OOFEM 3.0
Loading...
Searching...
No Matches
fcm.C
Go to the documentation of this file.
1/*
2 *
3 * ##### ##### ###### ###### ### ###
4 * ## ## ## ## ## ## ## ### ##
5 * ## ## ## ## #### #### ## # ##
6 * ## ## ## ## ## ## ## ##
7 * ## ## ## ## ## ## ## ##
8 * ##### ##### ## ###### ## ##
9 *
10 *
11 * OOFEM : Object Oriented Finite Element Code
12 *
13 * Copyright (C) 1993 - 2025 Borek Patzak
14 *
15 *
16 *
17 * Czech Technical University, Faculty of Civil Engineering,
18 * Department of Structural Mechanics, 166 29 Prague, Czech Republic
19 *
20 * This library is free software; you can redistribute it and/or
21 * modify it under the terms of the GNU Lesser General Public
22 * License as published by the Free Software Foundation; either
23 * version 2.1 of the License, or (at your option) any later version.
24 *
25 * This program is distributed in the hope that it will be useful,
26 * but WITHOUT ANY WARRANTY; without even the implied warranty of
27 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
28 * Lesser General Public License for more details.
29 *
30 * You should have received a copy of the GNU Lesser General Public
31 * License along with this library; if not, write to the Free Software
32 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
33 */
34
35#include "fcm.h"
36#include "gausspoint.h"
37#include "floatmatrix.h"
38#include "floatarray.h"
39#include "contextioerr.h"
40#include "mathfem.h"
41#include "stressvector.h"
42
43#include <cstring>
44
45namespace oofem {
46FCMMaterial :: FCMMaterial(int n, Domain *d) : StructuralMaterial(n, d),
47 linearElasticMaterial(n, d),
48 ecsMethod(ECSM_Unknown)
49{}
50
51
52bool
53FCMMaterial :: hasMaterialModeCapability(MaterialMode mode) const
54{
55 return mode == _3dMat || mode == _PlaneStress || mode == _PlaneStrain;
56}
57
58
59void
60FCMMaterial :: giveRealStressVector(FloatArray &answer, GaussPoint *gp,
61 const FloatArray &totalStrain,
62 TimeStep *tStep) const
63//
64// returns real stress vector in 3d stress space of receiver according to
65// previous level of stress and current
66// strain increment, the only way, how to correctly update gp records
67//
68{
69 // g = global direction
70 // l = local (crack) direction
71 // normal = normal components only (in direction of crack planes, "l")
72
73 // stress and strain transformation matrices
74 FloatMatrix epsG2L, sigmaL2G, sigmaG2L;
75 // stiffness matrices (total, elastic, cracking)
76 FloatMatrix D, De, Dcr, DeRed;
77 FloatMatrix principalDirs;
78
79 // STRESSES
80 FloatArray stressIncrement_l, stressVector_g, trialStress_g, trialStress_l, sigmaResid, sigmaElast_l, sigmaCrack_l, sigmaResid_old;
81 FloatArray stressVector_l;
82
83 // STRAINS
84 FloatArray reducedStrain_g, strainIncrement_g, strainIncrement_l, crackStrain, crackStrainIncrement, elasticStrain_l, elasticStrainIncrement_l; // FloatArray reducedStrain_l
85
86 FCMMaterialStatus *status = static_cast< FCMMaterialStatus * >( this->giveStatus(gp) );
87 MaterialMode mMode = gp->giveMaterialMode();
88
89 int nCr = status->giveNumberOfCracks();
90 int nMaxCr = status->giveMaxNumberOfCracks(gp);
91
92 double maxErr = 0.;
93 double maxTau = 0.;
94
95 double G;
96
97 double DII;
98
99 double gamma_cr, d_gamma_cr;
100 double tau_el, tau_cr;
101 double d_tau;
102 double d_tau_old = 0.;
103 bool illinoisFlag = false;
104
105 int iterLimitGradient = 20;
106 int iterLimitNormal = 100;
107 int iterLimitShear = 100;
108 int indexLimit = 10;
109
110 FloatArray tempNormalCrackStrain, normalCrackStrain;
111 // FloatArray normalStrainVector_l;
112 FloatArray normalStrainIncrement_l;
113 FloatArray normalStressVector_l;
114 FloatArray normalCrackStrainIncrement;
115 FloatArray normalElasticStrainIncrement_l;
116 FloatArray normalElasticStrain_l;
117
118
119 // for the bisection method:
120 int index, indexOld, indexCount;
121
122 bool changeIndexFlag = false;
123
124 bool plus, minus;
125 FloatArray crackStrainPlus, crackStrainMinus, residPlus, residMinus; //local
126
127 FloatArray helpRF;
128
129 int iterN;
130
131 plus = false;
132 minus = false;
133
134 sigmaResid_old.resize(nMaxCr);
135
136 index = indexOld = indexCount = 0;
137
138
139 this->initTempStatus(gp);
140
141 if ( !this->isActivated(tStep) ) {
142 FloatArray zeros;
143 zeros.resize( StructuralMaterial :: giveSizeOfVoigtSymVector( gp->giveMaterialMode() ) );
144 zeros.zero();
145 status->letTempStrainVectorBe(zeros);
146 status->letTempStressVectorBe(zeros);
147 answer = zeros;
148 return;
149 }
150
151 // get elastic stiffness matrix
152 this->giveStiffnessMatrix(De, ElasticStiffness, gp, tStep);
153
154 // equilibrated global stress
155 stressVector_g = status->giveStressVector();
156
157
158 this->giveStressDependentPartOfStrainVector(reducedStrain_g, gp, totalStrain, tStep, VM_Incremental);
159 strainIncrement_g.beDifferenceOf( reducedStrain_g, status->giveStrainVector() );
160
161 trialStress_g.beProductOf(De, strainIncrement_g);
162 trialStress_g.add(stressVector_g);
163
164 // ELASTIC CASE - NO CRACKING SO FAR (in prevs steps)
165 if ( nCr == 0 ) {
166
167 this->computePrincipalValDir(trialStress_l, principalDirs, trialStress_g, principal_stress);
168
169 // REMAINS ELASTIC or element cracking is prevented from above
170 if ( !this->isStrengthExceeded( principalDirs, gp, tStep, 1, trialStress_l.at(1) ) ) {
171 answer = trialStress_g;
172 status->letTempStrainVectorBe(totalStrain);
173 status->letTempStressVectorBe(trialStress_g);
174 return;
175
176 // STARTS CRACKING - 1st crack
177 } else {
178 this->initializeCrack(gp, tStep, principalDirs, 1);
179
180 for ( int iCrack = 2; iCrack <= min(nMaxCr, this->nAllowedCracks); iCrack++ ) {
181 if ( this->isStrengthExceeded( principalDirs, gp, tStep, iCrack, trialStress_l.at(iCrack) ) ) { // for "nonlocal model"
182 this->initializeCrack(gp, tStep, principalDirs, iCrack);
183 }
184 }
185 }
186
187 // CRACKING OCCURED IN PREVIOUS STEPS
188 } else {
189 // check other principal directions - if strength is not exceeded, take transformation matrices from previous analysi
190 if ( nCr < nMaxCr ) {
191
192 principalDirs = status->giveCrackDirs();
193
194 for ( int iCrack = nCr + 1; iCrack <= min(nMaxCr, this->nAllowedCracks); iCrack++ ) {
195 // condition to guarantee that third crack won't be initiated if second crack does not exist
196 if ( ( status->giveNumberOfTempCracks() + 1 ) < iCrack ) {
197 break;
198 }
199
200 if ( !this->checkStrengthCriterion(principalDirs, trialStress_g, gp, tStep, iCrack) ) { // = strength of composite
201 // second and third crack is now initialized
202 this->initializeCrack(gp, tStep, principalDirs, iCrack);
203 }
204 } // end for
205 } // end nMaxCr
206 } // cracking in previous steps
207
208
209 epsG2L = status->giveG2LStrainVectorTransformationMtrx();
210 sigmaG2L = status->giveG2LStressVectorTransformationMtrx();
211 sigmaL2G = status->giveL2GStressVectorTransformationMtrx();
212
213
214 // PREPARE STRAIN AND STRESS INCREMENT IN CRACK DIRECTIONS
215
216 // strain and strain increment in crack direction
217 strainIncrement_l.beProductOf(epsG2L, strainIncrement_g); // from global to local
218
219 // stress increment in crack direction
220 stressIncrement_l.beProductOf(De, strainIncrement_l);
221
222 // (equilibrated) stress in crack direction
223
224 stressVector_l.beProductOf(sigmaG2L, status->giveStressVector() );
225
226
227 // SIMILARLY TO RCM2: THE STRESS INCREMENT IN ELASTIC
228 // AND CRACKING UNIT MUST BE EQUAL - ITERATE
229
230 // average shear modulus
231 G = computeOverallElasticShearModulus(gp, tStep);
232
233 DeRed = De;
234 DeRed.resizeWithData(nMaxCr,nMaxCr);
235
236 crackStrain = status->giveCrackStrainVector();
237
238 tempNormalCrackStrain = status->giveCrackStrainVector();
239 tempNormalCrackStrain.resize(nMaxCr);
240
241 normalCrackStrain = tempNormalCrackStrain;
242
243
244 normalStressVector_l = stressVector_l;
245 normalStressVector_l.resize(nMaxCr);
246
247 normalStrainIncrement_l = strainIncrement_l;
248 normalStrainIncrement_l.resizeWithValues(nMaxCr);
249
250 for ( int iterG = 1; iterG <= iterLimitGlobal; iterG++ ) {
251
252 if ( iterG == 1 ) {
253 // residuum
254 sigmaResid = stressIncrement_l;
255 }
256
257 for ( int indexCount = 1; indexCount <= indexLimit; indexCount++ ) {
258
259 if (indexCount > 1 ) {
260 if ( !changeIndexFlag ) {// same index as before - do not even try to fix it
261 OOFEM_WARNING("FCM: maximum error on the same component, max. stress error %f, index %d, indexCount %d, element %d\n", maxErr, index, indexCount, gp->giveElement()->giveNumber() );
262
263 break;
264 }
265 }
266
267 minus = plus = false;
268 illinoisFlag = false;
269 changeIndexFlag = false;
270 iterN = 0;
271 sigmaResid_old.zero();
272
273
274 do { // find convergence in the crack directions
275
276 iterN++;
277
278 if ( (iterN == 1) && (indexCount == 1) ) {
279 // shrink the total residuum
280 sigmaResid.resizeWithValues(nMaxCr);
281 }
282
283 this->giveNormalLocalCrackedStiffnessMatrix(Dcr, TangentStiffness, gp, tStep);
284
285 D = DeRed;
286 D.add(Dcr);
287
288 D.solveForRhs(sigmaResid, normalCrackStrainIncrement);
289
290 // update and store cracking strain
291 tempNormalCrackStrain.add(normalCrackStrainIncrement);
292 status->setTempNormalCrackStrainVector(tempNormalCrackStrain);
293
294 // update statuses
295 this->updateCrackStatus(gp);
296
297 // compute and compare stress in elastic and cracking unit:
298 // CR
299 this->giveNormalLocalCrackedStiffnessMatrix(Dcr, SecantStiffness, gp, tStep);
300
301 sigmaCrack_l.resize( nMaxCr );
302 for ( int i = 1; i <= nMaxCr; i++ ) {
303 sigmaCrack_l.at(i) = Dcr.at(i, i) * tempNormalCrackStrain.at(i);
304 }
305
306 // compute and compare stress in elastic and cracking unit:
307 // increment of elastic strain
308 // d_eps_el = d_eps - d_eps_cr
309 // d_eps_cr = eps_cr_n - eps_cr_n-1
310 normalElasticStrainIncrement_l = normalStrainIncrement_l;
311 normalElasticStrainIncrement_l.subtract(tempNormalCrackStrain);
312 normalElasticStrainIncrement_l.add(normalCrackStrain);
313
314 // d_sigma = Del * d_eps_el
315 sigmaElast_l.beProductOf(DeRed, normalElasticStrainIncrement_l);
316
317 // sigma_n = sigma_n-1 + d_sigma
318 sigmaElast_l.add(normalStressVector_l);
319
320 // residuum
321 sigmaResid = sigmaElast_l;
322 sigmaResid.subtract(sigmaCrack_l);
323
324 maxErr = 0.;
325 indexOld = index;
326
327 for ( int i = 1; i <= sigmaResid.giveSize(); i++ ) {
328 if ( fabs( sigmaResid.at(i) ) > maxErr ) {
329 maxErr = fabs( sigmaResid.at(i) );
330
331 // if (indexCount == 1) {
332 index = i;
333 // }
334 }
335 }
336
337
338
339 if ( maxErr <= fcm_TOLERANCE * this->giveTensileStrength(gp, tStep) ) {
340 break; // convergence criteria satisfied
341 }
342
343 // try only gradient method and do not slow down with the preparation for the bisection method
344 if ( indexCount == 1 ) {
345
346 if (iterN < iterLimitGradient/2.) {
347 continue;
348
349 }
350 }
351 // identify the biggest error only in the first run of the gradient method
352 //} else {
353
354 if (index != indexOld) {
355 minus = plus = false;
356 // crackStrainPlus.resize(nMaxCr);
357 /*crackStrainPlus.zero();
358 //crackStrainMinus.resize(nMaxCr);
359 crackStrainMinus.zero();
360 // residPlus.resize(nMaxCr);
361 residPlus.zero();
362 //residMinus.resize(nMaxCr);
363 residMinus.zero();*/
364 }
365
366
367
368 if ( sigmaResid.at(index) > 0 ) {
369 if ( !plus) { // first positive error
370 plus = true;
371 residPlus = sigmaResid;
372 crackStrainPlus = tempNormalCrackStrain;
373 } else if ( fabs(sigmaResid.at(index)) < fabs(residPlus.at(index)) ) { // smaller error
374 residPlus = sigmaResid;
375 crackStrainPlus = tempNormalCrackStrain;
376 }
377
378 } else {
379 if ( !minus) { // first negative error
380 minus = true;
381 residMinus = sigmaResid;
382 crackStrainMinus = tempNormalCrackStrain;
383 }
384
385 else if ( fabs(sigmaResid.at(index)) < fabs(residMinus.at(index)) ) { // smaller error
386 residMinus = sigmaResid;
387 crackStrainMinus = tempNormalCrackStrain;
388 }
389 }
390
391
392 // if ( (! plus) || (! minus) ) {
393
394 if ( ( plus) && ( minus) ) {
395 if ( ( indexCount == 1 ) && ( iterN < iterLimitGradient ) ) {
396 continue;
397 }
398
399
400
401 //( iterN < iterLimitGradient ) || (! plus) || (! minus) ) {
402
403 } else {
404 continue;
405 }
406
407 do { // perform bisection method
408
409 // check if the cracking strain at index is negative and try to fix it
410 tempNormalCrackStrain.zero();
411
412 if( (crackStrainMinus.at(index) < 0.) || (crackStrainPlus.at(index) < 0.) ) {
413 tempNormalCrackStrain.add(crackStrainMinus);
414 tempNormalCrackStrain.add(crackStrainPlus);
415 tempNormalCrackStrain.times(0.5);
416
417 } else {
418 // OOFEM_WARNING("Trying regula-falsi iterations\n");
419
420 // try regula-falsi
421 // eps_cr_new = (eps_cr_min * resid_plus - eps_cr_plus * resid_minus) / (resid_plus - resid_minus)
422 // https://en.wikipedia.org/wiki/False_position_method#The_Regula_Falsi_.28False_Position.29_Method
423
424 if (illinoisFlag) { // same sign of the error (at least) twice in a row
425
426 if ( sigmaResid.at(index) == residPlus.at(index) ) { // error is positive
427
428 tempNormalCrackStrain.add(crackStrainMinus);
429 tempNormalCrackStrain.times(residPlus.at(index) );
430
431 helpRF = crackStrainPlus;
432 helpRF.times( 0.5 * residMinus.at(index) );
433
434 tempNormalCrackStrain.subtract(helpRF);
435 tempNormalCrackStrain.times( 1. / ( residPlus.at(index) - 0.5 * residMinus.at(index) ) );
436
437 } else { // error is negative
438
439 tempNormalCrackStrain.add(crackStrainMinus);
440 tempNormalCrackStrain.times( 0.5 * residPlus.at(index) );
441
442 helpRF = crackStrainPlus;
443 helpRF.times( residMinus.at(index) );
444
445 tempNormalCrackStrain.subtract(helpRF);
446 tempNormalCrackStrain.times( 1. / ( 0.5 * residPlus.at(index) - residMinus.at(index) ) );
447
448 }
449
450 } else { // different sign of the error (consecutive errors)
451 // perform standard regula-falsi
452
453 tempNormalCrackStrain.add(crackStrainMinus);
454 tempNormalCrackStrain.times(residPlus.at(index) );
455
456 helpRF = crackStrainPlus;
457 helpRF.times(residMinus.at(index) );
458
459 tempNormalCrackStrain.subtract(helpRF);
460 tempNormalCrackStrain.times( 1. / ( residPlus.at(index) - residMinus.at(index) ) );
461 }
462 }
463
464
465 // too small strain and the crack has not existed before
466 if ( ( fabs (crackStrainMinus.at(index) ) <= fcm_THRESHOLD_CRACK_STRAIN ) &&
467 ( fabs (crackStrainPlus.at(index) ) <= fcm_THRESHOLD_CRACK_STRAIN ) &&
468 ( status->giveCrackStatus(index) == 0 ) ) {
469 tempNormalCrackStrain.at(index) = 0.;
470 }
471
472 // replace negative cracking strain by a zero strain
473 if (tempNormalCrackStrain.at(index) < 0.) {
474 tempNormalCrackStrain.at(index) = 0.;
475 }
476
477 status->setTempNormalCrackStrainVector(tempNormalCrackStrain);
478 // update statuses - to have correct stiffnesses (unlo-relo vs. softening)
479 this->updateCrackStatus(gp);
480
481 // compute and compare stress in elastic and cracking unit:
482 // CR
483 this->giveNormalLocalCrackedStiffnessMatrix(Dcr, SecantStiffness, gp, tStep);
484 sigmaCrack_l.resize( nMaxCr );
485 for ( int i = 1; i <= nMaxCr; i++ ) {
486 sigmaCrack_l.at(i) = Dcr.at(i, i) * tempNormalCrackStrain.at(i);
487 }
488
489 // increment of elastic strain
490 // d_eps_el = d_eps - d_eps_cr
491 // d_eps_cr = eps_cr_n - eps_cr_n-1
492 normalElasticStrainIncrement_l = normalStrainIncrement_l;
493 normalElasticStrainIncrement_l.subtract(tempNormalCrackStrain);
494 normalElasticStrainIncrement_l.add(normalCrackStrain);
495
496 // d_sigma = Del * d_eps_el
497 sigmaElast_l.beProductOf(DeRed, normalElasticStrainIncrement_l);
498
499 // sigma_n = sigma_n-1 + d_sigma
500 sigmaElast_l.add(normalStressVector_l);
501
502 // residuum
503 sigmaResid = sigmaElast_l;
504 sigmaResid.subtract(sigmaCrack_l);
505
506
507 /* if ( (status->giveTempCrackStatus(index) == 4) && (tempNormalCrackStrain.at(index) < fcm_SMALL_STRAIN) ) {
508 sigmaResid.at(index) = 0.;
509 }*/
510
511
512
513 // adjust residuum if the crack does not exist
514 /* if ( status->giveTempCrackStatus(index) == 0 ) {
515 sigmaResid.at(index) = 0.;
516 }*/
517
518 // condition for the given stress component
519 if ( fabs( sigmaResid.at(index) ) < fcm_TOLERANCE * this->giveTensileStrength(gp, tStep) ) {
520 changeIndexFlag = true;
521 break;
522 }
523
524
525 if ( (sigmaResid.at(index) >= 0.) && ( sigmaResid.at(index) < residPlus.at(index) ) ) {
526 residPlus = sigmaResid;
527 crackStrainPlus = tempNormalCrackStrain;
528
529 } else if ( (sigmaResid.at(index) <= 0.) && ( fabs(sigmaResid.at(index)) < fabs(residMinus.at(index) ) ) ) {
530
531 residMinus = sigmaResid;
532 crackStrainMinus = tempNormalCrackStrain;
533
534 } else { // error is increasing -> change index
535 changeIndexFlag = true;
536 }
537
538
539 // cracking strain at index component is close to zero
540 if ( ( ( fabs( crackStrainMinus.at(index) - crackStrainPlus.at(index) ) < 1.e-16 ) &&
541 ( max( fabs( crackStrainMinus.at(index) ), fabs( crackStrainPlus.at(index) ) ) < 1.e-14 ) ) ||
542 ( status->giveTempCrackStatus(index) == pscm_NONE ) ) {
543 tempNormalCrackStrain.zero();
544 tempNormalCrackStrain.add(crackStrainMinus);
545 tempNormalCrackStrain.add(crackStrainPlus);
546 tempNormalCrackStrain.times(0.5);
547
548 tempNormalCrackStrain.at(index) = 0.;
549
550#if DEBUG
551 // OOFEM_WARNING("Fixed crack model: cracking strain component %d set to zero", index);
552#endif
553
554
555 status->setTempNormalCrackStrainVector(tempNormalCrackStrain);
556
557 // update statuses
558 this->updateCrackStatus(gp);
559
560 // compute and compare stress in elastic and cracking unit:
561 for ( int i_count = 1; i_count <= 2; i_count++ ) {
562
563 // EL
564 // increment of elastic strain
565 // d_eps_el = d_eps - d_eps_cr
566 // d_eps_cr = eps_cr_n - eps_cr_n-1
567 normalElasticStrainIncrement_l = normalStrainIncrement_l;
568 normalElasticStrainIncrement_l.subtract(tempNormalCrackStrain);
569 normalElasticStrainIncrement_l.add(normalCrackStrain);
570 // d_sigma = Del * d_eps_el
571 sigmaElast_l.beProductOf(DeRed, normalElasticStrainIncrement_l);
572
573 // eliminating stress error at index component
574 tempNormalCrackStrain.at(index) = sigmaElast_l.at(index) / this->giveCrackingModulus(SecantStiffness, gp, tStep, index);
575 status->setTempNormalCrackStrainVector(tempNormalCrackStrain);
576
577 // update statuses
578 this->updateCrackStatus(gp);
579 }
580
581 // stress in cracking unit
582 this->giveNormalLocalCrackedStiffnessMatrix(Dcr, SecantStiffness, gp, tStep);
583 sigmaCrack_l.resize( nMaxCr );
584 for ( int i = 1; i <= nMaxCr; i++ ) {
585 sigmaCrack_l.at(i) = Dcr.at(i, i) * tempNormalCrackStrain.at(i);
586 }
587
588 // reevaluate the residuum
589 sigmaResid = sigmaElast_l;
590 sigmaResid.subtract(sigmaCrack_l);
591
592 } // end - crack strain is close to zero
593
594 // evaluate maximum error
595 indexOld = index;
596
597 maxErr = 0.;
598 for ( int i = 1; i <= sigmaResid.giveSize(); i++ ) {
599 if ( fabs( sigmaResid.at(i) ) > maxErr ) {
600 maxErr = fabs( sigmaResid.at(i) );
601 index = i;
602 }
603 }
604
605 // biggest error on different index
606 if (indexOld != index) {
607 changeIndexFlag = true;
608 illinoisFlag = false;
609
610 } else { // evaluate
611
612 if ( sgn(sigmaResid.at(index) ) == sgn(sigmaResid_old.at(index) ) ) {
613 illinoisFlag = true;
614 } else {
615 illinoisFlag = false;
616 }
617
618 sigmaResid_old = sigmaResid;
619
620 }
621
622
623 // end loop bisection method
624 } while ( (!changeIndexFlag) && (iterN <= iterLimitNormal) && ( fabs( sigmaResid.at(index) ) < fcm_TOLERANCE * this->giveTensileStrength(gp, tStep) ) );
625
626
627 // end normal loop
628 } while ( ( !changeIndexFlag ) && ( maxErr > fcm_TOLERANCE * this->giveTensileStrength(gp, tStep) ) && (iterN <= iterLimitNormal) );
629
630 if ( maxErr <= fcm_TOLERANCE * this->giveTensileStrength(gp, tStep) ) {
631 break; // break index loop
632 }
633
634 if ( indexCount >= indexLimit/2.) {
635 OOFEM_WARNING("FCM: slowly converging equilibrium in crack direction!, max. stress error %f, index %d, indexCount %d, element %d\n", maxErr, index, indexCount, gp->giveElement()->giveNumber() );
636 }
637
638 if ( indexCount >= indexLimit ) {
639 OOFEM_WARNING("FCM: slowly converging equilibrium in crack direction!, max. stress error %f, index %d, indexCount %d, element %d\n", maxErr, index, indexCount, gp->giveElement()->giveNumber() );
640 break;
641 }
642
643
644 } // end index loop in normal direction
645
646 // FIND EQUILIBRIUM FOR SHEAR
647 tau_cr = 0.;
648 DII = 0.;
649
650
651 double d_gamma_el;
652
653 for ( int i = nMaxCr+1; i <= crackStrain.giveSize(); i++ ) {
654
655 plus = false;
656 minus = false;
657
658 crackStrainPlus.resize(1);
659 crackStrainPlus.zero();
660 crackStrainMinus.resize(1);
661 crackStrainMinus.zero();
662 residPlus.resize(1);
663 residPlus.zero();
664 residMinus.resize(1);
665 residMinus.zero();
666
667
668 d_tau = G * strainIncrement_l.at(i);
669 gamma_cr = status->giveCrackStrain(i);
670
671 for ( int iterS = 1; iterS <= iterLimitShear; iterS++ ) {
672
673 switch ( mMode ) {
674 case _PlaneStress:
675
676 DII = computeTotalD2Modulus(gp, tStep, 6);
677 break;
678
679 case _PlaneStrain: // check
680 DII = computeTotalD2Modulus(gp, tStep, 6);
681 break;
682
683 case _3dMat:
684 DII = computeTotalD2Modulus(gp, tStep, i);
685 break;
686
687 default:
688 OOFEM_ERROR( "Material mode %s not supported", __MaterialModeToString(mMode) );
689 }
690
691
692 if ( ( !plus ) || ( !minus ) ) {
693
694 if (iterS == 1) {
695 d_gamma_cr = ( stressVector_l.at(i) - DII * status->giveCrackStrain(i) + d_tau ) / ( G + DII );
696 } else {
697 d_gamma_cr = ( d_tau ) / ( G + DII );
698 }
699
700 gamma_cr += d_gamma_cr;
701 } else {
702 // same sign of the cracking strain
703 if (sgn( crackStrainMinus.at(1) ) == sgn( crackStrainPlus.at(1) ) ) {
704
705 if (illinoisFlag) { // same sign of the error (at least) twice in a row
706
707 if (d_tau == residPlus.at(1) ) { // error is positive
708 gamma_cr = ( crackStrainMinus.at(1) * residPlus.at(1) - 0.5 * crackStrainPlus.at(1) * residMinus.at(1) ) / ( residPlus.at(1) - 0.5 * residMinus.at(1) );
709
710 } else {
711 gamma_cr = ( 0.5 * crackStrainMinus.at(1) * residPlus.at(1) - crackStrainPlus.at(1) * residMinus.at(1) ) / ( 0.5 * residPlus.at(1) - residMinus.at(1) );
712 }
713
714 // different sign of the error (consecutive errors)
715 } else {
716 gamma_cr = ( crackStrainMinus.at(1) * residPlus.at(1) - crackStrainPlus.at(1) * residMinus.at(1) ) / ( residPlus.at(1) - residMinus.at(1) );
717 }
718 } else { // different sign of the cracking strain
719 gamma_cr = ( crackStrainMinus.at(1) + crackStrainPlus.at(1) ) / 2.;
720 }
721 }
722
723
724 status->setTempCrackStrain(i, gamma_cr);
725 this->updateCrackStatus(gp);
726
727 switch ( mMode ) {
728 case _PlaneStress:
729
730 maxTau = this->maxShearStress(gp, tStep, 6);
731 tau_cr = computeTotalD2Modulus(gp, tStep, 6) * gamma_cr;
732 break;
733
734 case _PlaneStrain: // check
735
736 maxTau = this->maxShearStress(gp, tStep, 6);
737 tau_cr = computeTotalD2Modulus(gp, tStep, 6) * gamma_cr;
738 break;
739
740 case _3dMat:
741
742 maxTau = this->maxShearStress(gp, tStep, i);
743 tau_cr = computeTotalD2Modulus(gp, tStep, i) * gamma_cr;
744 break;
745
746 default:
747 OOFEM_ERROR( "Material mode %s not supported", __MaterialModeToString(mMode) );
748 }
749
750 if( fabs(tau_cr) > maxTau) {
751
752 if ( ( !plus ) || ( !minus ) ) {
753
754 // d_eps_el = (tau_max - tau_n-1)/G
755 d_gamma_el = ( sgn(tau_cr)* maxTau - stressVector_l.at(i) ) / G;
756
757 // eps_cr_n = eps_cr_n-1 + d_eps - d_eps_el
758 gamma_cr = status->giveCrackStrain(i) + strainIncrement_l.at(i) - d_gamma_el;
759
760 // maxTau can be dependent on the actual cracking strain
761 // (such as in the case of FRC)
762
763 } else {
764 d_gamma_el = strainIncrement_l.at(i) - (gamma_cr - status->giveCrackStrain(i) );
765 }
766
767 status->setTempCrackStrain(i, gamma_cr);
768 this->updateCrackStatus(gp);
769
770 // reevaluate TauMax
771 switch ( mMode ) {
772 case _PlaneStress:
773
774 maxTau = this->maxShearStress(gp, tStep, 6);
775 break;
776
777 case _PlaneStrain: // check
778
779 maxTau = this->maxShearStress(gp, tStep, 6);
780 break;
781
782 case _3dMat:
783
784 maxTau = this->maxShearStress(gp, tStep, i);
785 break;
786
787 default:
788 OOFEM_ERROR( "Material mode %s not supported", __MaterialModeToString(mMode) );
789 }
790
791 tau_el = stressVector_l.at(i) + d_gamma_el * G;
792
793 tau_cr = min( maxTau, DII * fabs(gamma_cr) );
794 tau_cr *= sgn(gamma_cr);
795
796 } else { // tau is smaller than tau_max
797
798 // old shear stress
799 tau_el = stressVector_l.at(i);
800 // shear stress increment
801 tau_el += G * ( strainIncrement_l.at(i) - gamma_cr + status->giveCrackStrain(i) );
802 }
803
804 // compute residuum = sig_el - sig_cr
805 d_tau = tau_el - tau_cr;
806
807 if ( fabs( d_tau ) < fcm_TOLERANCE * this->giveTensileStrength(gp, tStep) ) {
808 break;
809
810 } else {
811
812 if (d_tau < 0.) {
813 minus = true;
814
815 residMinus.at(1) = d_tau;
816 crackStrainMinus.at(1) = gamma_cr;
817
818 } else {
819 plus = true;
820
821 residPlus.at(1) = d_tau;
822 crackStrainPlus.at(1) = gamma_cr;
823 }
824
825 // compare signs of the error
826 if (sgn(d_tau_old) == sgn(d_tau) ) { // same sign of the error as before
827 illinoisFlag = true;
828
829 } else {
830 illinoisFlag = false;
831 }
832
833 d_tau_old = d_tau;
834
835 }
836
837#if DEBUG
838 if ( iterS > iterLimitShear*0.3) {
839 OOFEM_WARNING("Fixed crack model: Local equilibrium in shear not converging!, max. stress error %f, index %d, iter %d, element %d\n", fabs(d_tau), i, iterS, gp->giveElement()->giveNumber() );
840 }
841#endif
842
843
844 if ( iterS == iterLimitShear ) {
845 OOFEM_WARNING("Fixed crack model: Local equilibrium in shear direction(s) not reached!, max. stress error %f, element %d\n", fabs(d_tau), gp->giveElement()->giveNumber() );
846 }
847
848 } // i-th shear component equlibrium loop
849
850 } // shear components loop
851
852 // compute final stress residuum
853 crackStrain = status->giveTempCrackStrainVector();
854
855 // el
856 elasticStrainIncrement_l = strainIncrement_l;
857 elasticStrainIncrement_l.subtract(crackStrain);
858 elasticStrainIncrement_l.add(status->giveCrackStrainVector() );
859
860 sigmaElast_l.beProductOf(De, elasticStrainIncrement_l);
861 sigmaElast_l.add(stressVector_l);
862
863 this->giveTotalLocalCrackedStiffnessMatrix(Dcr, SecantStiffness, gp, tStep);
864
865 sigmaCrack_l.resize( crackStrain.giveSize() );
866
867 for ( int i = 1; i <= crackStrain.giveSize(); i++ ) {
868 sigmaCrack_l.at(i) = Dcr.at(i, i) * crackStrain.at(i);
869
870 if (this->isThisShearComponent(gp, i) ) {
871
872 switch ( mMode ) {
873 case _PlaneStress:
874 maxTau = this->maxShearStress(gp, tStep, 6);
875 break;
876
877 case _PlaneStrain:
878 maxTau = this->maxShearStress(gp, tStep, 6);
879 break;
880
881 case _3dMat:
882 maxTau = this->maxShearStress(gp, tStep, i);
883 break;
884
885 default:
886 OOFEM_ERROR( "Material mode %s not supported", __MaterialModeToString(mMode) );
887 }
888
889
890 if ( fabs( sigmaCrack_l.at(i) ) > maxTau ) {
891 sigmaCrack_l.at(i) = sgn(sigmaCrack_l.at(i)) * maxTau;
892 }
893 }
894 }
895
896 // residuum
897 sigmaResid = sigmaElast_l;
898 sigmaResid.subtract(sigmaCrack_l);
899
900 maxErr = 0.;
901 for ( int i = 1; i <= sigmaResid.giveSize(); i++ ) {
902 if ( fabs( sigmaResid.at(i) ) > maxErr ) {
903 maxErr = fabs( sigmaResid.at(i) );
904 }
905 }
906
907 if ( maxErr < fcm_TOLERANCE * this->giveTensileStrength(gp, tStep) ) {
908 break;
909 }
910
911 if ( iterG == iterLimitGlobal ) {
912 OOFEM_WARNING("Fixed crack model: Local equilibrium not reached!, max. stress error %f, element %d\n", maxErr, gp->giveElement()->giveNumber() );
913 }
914
915 if ( iterG > iterLimitGlobal/3 ) {
916 OOFEM_WARNING("Slowly converning (normal vs. shear interaction)!, iter %d/ iter limit %d, max. stress error %f, element %d\n", iterG, iterLimitGlobal, maxErr, gp->giveElement()->giveNumber() );
917 }
918
919 } // global loop
920
921
922 // correct non-zeros, statuses, etc...
923 for ( int i = 1; i <= nMaxCr; i++ ) {
924#if DEBUG
925 // check for NaNs
926 if ( crackStrain.at(i) != crackStrain.at(i) ) {
927 OOFEM_ERROR( "crack strain is NaN: %f", crackStrain.at(i) );
928 }
929#endif
930
931 if ( ( status->giveTempCrackStatus(i) == pscm_NONE ) && ( crackStrain.at(i) <= 0. ) ) {
932 status->setTempCrackStrain(i, 0.);
933 } else if ( ( status->giveTempCrackStatus(i) == pscm_JUST_INIT ) && ( crackStrain.at(i) <= fcm_SMALL_STRAIN ) ) {
934 if ( i + 1 > status->giveMaxNumberOfCracks(gp) ) {
935 status->setTempCrackStatus(i, pscm_NONE);
936 status->setTempCrackStrain(i, 0.);
937 if ( i == 1 ) {
938 principalDirs.zero();
939 status->setCrackDirs(principalDirs);
940 }
941
942 // can't change crack if the subsequent crack exists
943 } else if ( status->giveTempCrackStatus(i + 1) != pscm_NONE ) {} else {
944 status->setTempCrackStatus(i, pscm_NONE);
945 status->setTempCrackStrain(i, 0.);
946 if ( i == 1 ) {
947 principalDirs.zero();
948 status->setCrackDirs(principalDirs);
949 }
950 }
951 } else if ( crackStrain.at(i) < 0. ) {
952 status->setTempCrackStrain(i, 0.);
953 }
954 }
955
956 // at first crack initiation there should not be any other cracking strains
957 if ( ( status->giveTempCrackStatus(1) == pscm_JUST_INIT ) && ( status->giveTempCrackStatus(2) == pscm_NONE ) ) {
958 status->setTempCrackStrain(2, 0.);
959 status->setTempCrackStrain(3, 0.);
960 }
961
962 // convert from local to global coordinates
963 answer.beProductOf(sigmaL2G, sigmaElast_l);
964 status->letTempStrainVectorBe(totalStrain);
965 status->letTempStressVectorBe(answer);
966
967 return;
968}
969
970
971void
972FCMMaterial :: initializeCrack(GaussPoint *gp, TimeStep *tStep, FloatMatrix &base, int nCrack) const
973{
974 MaterialMode mMode = gp->giveMaterialMode();
975 FCMMaterialStatus *status = static_cast< FCMMaterialStatus * >( this->giveStatus(gp) );
976
977 FloatArray crackVector(3);
978 crackVector.zero();
979
980 FloatMatrix epsG2L, sigmaG2L, epsL2G, sigmaL2G;
981
982 int nMaxCracks;
983 double Le;
984
985 for ( int i = 1; i <= base.giveNumberOfRows(); i++ ) {
986 crackVector.at(i) = base.at(i, nCrack);
987 }
988
989 Le = this->giveCharacteristicElementLength(gp, crackVector);
990 status->setCharLength(nCrack, Le);
991
992 this->checkSnapBack(gp, tStep, nCrack);
993
994 status->setTempCrackStatus(nCrack, pscm_JUST_INIT);
995 //status->setTempCrackStatus(nCrack, pscm_SOFTENING);
996
997 // change the base crack vector and transformation matrices
998 // when the initiated crack is not the last possible one
999
1000 // return for the second crack in plane stress - base is set and transformation matrices are given by the first crack
1001 if ( ( mMode == _PlaneStress ) && ( nCrack == 2 ) ) {
1002 return;
1003 }
1004
1005 nMaxCracks = status->giveMaxNumberOfCracks(gp);
1006
1007 if ( nCrack <= nMaxCracks ) {
1008 // make sure that the matrix of base vectors is right-handed
1009 if ( nMaxCracks == 3 ) {
1010 base.at(1, 3) = base.at(2, 1) * base.at(3, 2) - base.at(3, 1) * base.at(2, 2);
1011 base.at(2, 3) = base.at(3, 1) * base.at(1, 2) - base.at(1, 1) * base.at(3, 2);
1012 base.at(3, 3) = base.at(1, 1) * base.at(2, 2) - base.at(2, 1) * base.at(1, 2);
1013 }
1014
1015 // status->setTempCrackDirs(base);
1016 status->setCrackDirs(base);
1017
1018 if ( mMode == _PlaneStress ) {
1019 sigmaG2L = this->givePlaneStressVectorTranformationMtrx(base, 0);
1020 epsG2L = this->give2DStrainVectorTranformationMtrx(base, 0);
1021
1022 sigmaL2G = this->givePlaneStressVectorTranformationMtrx(base, 1);
1023 epsL2G = this->give2DStrainVectorTranformationMtrx(base, 1);
1024 } else {
1025 sigmaG2L = this->giveStressVectorTranformationMtrx(base, 0);
1026 epsG2L = this->giveStrainVectorTranformationMtrx(base, 0);
1027
1028 sigmaL2G = this->giveStressVectorTranformationMtrx(base, 1);
1029 epsL2G = this->giveStrainVectorTranformationMtrx(base, 1);
1030 }
1031
1032 status->setG2LStressVectorTransformationMtrx(sigmaG2L);
1033 status->setG2LStrainVectorTransformationMtrx(epsG2L);
1034
1035 status->setL2GStressVectorTransformationMtrx(sigmaL2G);
1036 status->setL2GStrainVectorTransformationMtrx(epsL2G);
1037 }
1038}
1039
1040
1041bool
1042FCMMaterial :: isIntact(GaussPoint *gp, int icrack) const {
1043 FCMMaterialStatus *status = static_cast< FCMMaterialStatus * >( this->giveStatus(gp) );
1044
1045 if ( icrack >= 4 ) {
1046 OOFEM_ERROR("Unexpected crack number");
1047 }
1048
1049 if ( status->giveCrackStatus(icrack) != pscm_NONE ) {
1050 return false;
1051 } else if ( ( status->giveTempCrackStatus(icrack) != pscm_NONE ) && ( status->giveTempCrackStatus(icrack) != pscm_CLOSED ) ) {
1052 return false;
1053 } else {
1054 return true;
1055 }
1056}
1057
1058
1059bool
1060FCMMaterial :: isIntactForShear(GaussPoint *gp, int i) const {
1061 FCMMaterialStatus *status = static_cast< FCMMaterialStatus * >( this->giveStatus(gp) );
1062
1063 int normal_1, normal_2;
1064
1065 if ( i == 4 ) { // y-z
1066 normal_1 = 2;
1067 normal_2 = 3;
1068 } else if ( i == 5 ) { // x-z
1069 normal_1 = 1;
1070 normal_2 = 3;
1071 } else if ( i == 6 ) { // x-y
1072 normal_1 = 1;
1073 normal_2 = 2;
1074 } else {
1075 OOFEM_ERROR("Unexpected number for shear stress (must be either 4, 5 or 6).");
1076 }
1077
1078 // any crack had existed in previous steps - newly inserted condition
1079 if ( (status->giveCrackStatus(normal_1) != pscm_NONE ) || (status->giveCrackStatus(normal_2) != pscm_NONE ) ) {
1080 return false;
1081
1082 // } else if ( ( status->giveTempCrackStatus(normal_1) != pscm_NONE ) && ( status->giveTempCrackStatus(normal_1) != pscm_CLOSED ) ) {
1083 // return false;
1084 // } else if ( ( status->giveTempCrackStatus(normal_2) != pscm_NONE ) && ( status->giveTempCrackStatus(normal_2) != pscm_CLOSED ) ) {
1085 // return false;
1086 } else {
1087 return true;
1088 }
1089}
1090
1091
1092
1093bool
1094FCMMaterial :: isThisShearComponent(GaussPoint *gp, int component) const {
1095
1096 MaterialMode mMode = gp->giveMaterialMode();
1097
1098 switch ( mMode ) {
1099
1100 case _PlaneStress:
1101 if (component == 3) {
1102 return true;
1103 }
1104 break;
1105
1106 case _PlaneStrain:
1107 if (component == 4) {
1108 return true;
1109 }
1110 break;
1111
1112 case _3dMat:
1113
1114 if (component >=4 ) {
1115 return true;
1116 }
1117 break;
1118
1119 default:
1120 OOFEM_ERROR("Unsupported material mode");
1121 }
1122
1123 return false;
1124}
1125
1126
1127
1128
1129double
1130FCMMaterial :: computeNormalCrackOpening(GaussPoint *gp, int i) const {
1131 FCMMaterialStatus *status = static_cast< FCMMaterialStatus * >( this->giveStatus(gp) );
1132
1133 double crackOpening, N;
1134
1135 if ( i > status->giveNumberOfTempCracks() ) {
1136 crackOpening = 0;
1137 } else {
1138 crackOpening = max(status->giveCharLength(i) * status->giveTempCrackStrain(i), 0.);
1139 N = this->giveNumberOfCracksInDirection(gp, i);
1140 crackOpening /= N;
1141 }
1142
1143 return crackOpening;
1144}
1145
1146
1147
1148double
1149FCMMaterial :: computeMaxNormalCrackOpening(GaussPoint *gp, TimeStep *tStep, int i) const {
1150 FCMMaterialStatus *status = static_cast< FCMMaterialStatus * >( this->giveStatus(gp) );
1151
1152 double crackOpening, N;
1153
1154 if ( i > status->giveNumberOfTempCracks() ) {
1155 crackOpening = 0;
1156 } else {
1157 crackOpening = max(status->giveCharLength(i) * status->giveMaxCrackStrain(i), 0.);
1158 N = this->giveNumberOfCracksInDirection(gp, i);
1159 crackOpening /= N;
1160 }
1161
1162 return crackOpening;
1163}
1164
1165
1166
1167double
1168FCMMaterial :: computeShearSlipOnCrack(GaussPoint *gp, TimeStep *tStep, int icrack) const{
1169 MaterialMode mMode = gp->giveMaterialMode();
1170
1171 FCMMaterialStatus *status = static_cast< FCMMaterialStatus * >( this->giveStatus(gp) );
1172
1173 // total shear slip
1174 double slip = 0;
1175
1176 // number of cracks in the direction of iCrack
1177 int nCracks = this->giveNumberOfCracksInDirection(gp, icrack);
1178
1179 // shear directions on icrack plane
1180 int dir_j, dir_k;
1181
1182 // cracking shear strains
1183 double gamma_cr_ij, gamma_cr_ik;
1184
1185 // factor for redistribution of shear crack strain
1186 double factor_ij, factor_ik;
1187
1188 double u_ij, u_ik;
1189
1190 if ( status->giveTempCrackStatus(icrack) == pscm_NONE ) { // crack not initiated
1191 return slip;
1192 }
1193
1194 if ( icrack > 3 ) {
1195 OOFEM_ERROR("Unexpected value of index i (4, 5, 6 permitted only)");
1196 }
1197
1198
1199 if ( ( mMode == _PlaneStress ) || ( mMode == _PlaneStrain ) ) {
1200 // get shear strain
1201 if ( mMode == _PlaneStress ) {
1202 gamma_cr_ij = status->giveTempCrackStrain(3);
1203 } else {
1204 gamma_cr_ij = status->giveTempCrackStrain(6);
1205 }
1206
1207 factor_ij = 1.;
1208
1209 // determine if it necessary to split the cracking strain
1210 // both cracks exist in 2D
1211
1212 if ( ( icrack == 2 ) || ( status->giveTempCrackStatus(2) != pscm_NONE ) ) {
1213 if ( icrack == 1 ) {
1214 dir_j = 2;
1215 } else { // icrack == 2
1216 dir_j = 1;
1217 }
1218
1219
1220 // well this is quite unfortunate. The problem is that the shear slip should be redistributed according to the D2 moduli for the individual crack. In reality this is a big problem for the FRC-FCM problem beacuse it would have led to recursive calling (D2 modulus evaluates damage and damage is computed from shear slip etc.
1221
1222 factor_ij = this->computeShearStiffnessRedistributionFactor(gp, tStep, icrack, dir_j);
1223 }
1224
1225 slip = factor_ij * fabs(gamma_cr_ij) * status->giveCharLength(icrack) / nCracks;
1226 } else if ( mMode == _3dMat ) {
1227 if ( icrack == 1 ) {
1228 dir_j = 2;
1229 dir_k = 3;
1230
1231 gamma_cr_ij = status->giveTempCrackStrain(6);
1232 gamma_cr_ik = status->giveTempCrackStrain(5);
1233 } else if ( icrack == 2 ) {
1234 dir_j = 1;
1235 dir_k = 3;
1236
1237 gamma_cr_ij = status->giveTempCrackStrain(6);
1238 gamma_cr_ik = status->giveTempCrackStrain(4);
1239 } else { // icrack == 3
1240 dir_j = 1;
1241 dir_k = 2;
1242
1243 gamma_cr_ij = status->giveTempCrackStrain(5);
1244 gamma_cr_ik = status->giveTempCrackStrain(4);
1245 }
1246
1247
1248 factor_ij = 1.;
1249 factor_ik = 1.;
1250
1251
1252 if ( ( status->giveTempCrackStatus(dir_j) != pscm_NONE ) || ( status->giveTempCrackStatus(dir_k) != pscm_NONE ) ) {
1253 if ( status->giveTempCrackStatus(dir_j) != pscm_NONE ) {
1254 factor_ij = this->computeShearStiffnessRedistributionFactor(gp, tStep, icrack, dir_j);
1255 }
1256
1257 if ( status->giveTempCrackStatus(dir_k) != pscm_NONE ) {
1258 factor_ik = this->computeShearStiffnessRedistributionFactor(gp, tStep, icrack, dir_k);
1259 }
1260 }
1261
1262 u_ij = factor_ij * fabs(gamma_cr_ij) * status->giveCharLength(icrack) / nCracks;
1263 u_ik = factor_ik * fabs(gamma_cr_ik) * status->giveCharLength(icrack) / nCracks;
1264
1265 slip = sqrt( pow(u_ij, 2) + pow(u_ik, 2) );
1266 } else {
1267 OOFEM_ERROR( "Material mode %s not supported", __MaterialModeToString(mMode) );
1268 }
1269
1270 return slip;
1271}
1272
1273
1274bool
1275FCMMaterial :: isStrengthExceeded(const FloatMatrix &base, GaussPoint *gp, TimeStep *tStep, int iCrack, double trialStress) const {
1276 if ( trialStress > this->giveTensileStrength(gp, tStep) ) {
1277 return true;
1278 } else {
1279 return false;
1280 }
1281}
1282
1283
1284
1285
1286double
1287FCMMaterial :: computeShearStiffnessRedistributionFactor(GaussPoint *gp, TimeStep *tStep, int ithCrackPlane, int jthCrackDirection) const {
1288 double factor_ij;
1289 double D2_i, D2_j;
1290
1291 D2_i = this->computeD2ModulusForCrack(gp, tStep, ithCrackPlane);
1292 D2_j = this->computeD2ModulusForCrack(gp, tStep, jthCrackDirection);
1293
1294 factor_ij = D2_j / ( D2_i + D2_j );
1295
1296 return factor_ij;
1297}
1298
1299
1300bool
1301FCMMaterial :: checkStrengthCriterion(FloatMatrix &newBase, const FloatArray &globalStress, GaussPoint *gp, TimeStep *tStep, int nCrack) const {
1302 FCMMaterialStatus *status = static_cast< FCMMaterialStatus * >( this->giveStatus(gp) );
1303
1304 double sigX, sigY, tau, sig2;
1305 FloatArray trialStress, planeStress, princStress, crackingStrain, shearStrains, newShearStrains;
1306 FloatMatrix sigmaG2L, princCrackDir, oldBase, princCrackDirExt;
1307
1308 sigmaG2L = status->giveG2LStressVectorTransformationMtrx();
1309
1310 // rotate stress to local coordinates
1311 // trialStress = globalStress;
1312 // trialStress.rotatedWith(sigmaG2L, 'n');
1313 trialStress.beProductOf(sigmaG2L, globalStress);
1314
1315 if ( status->giveMaxNumberOfCracks(gp) < 3 ) { // plane stress
1316 // test material strength but keep crack coordinates
1317 if ( this->isStrengthExceeded( newBase, gp, tStep, nCrack, trialStress.at(nCrack) ) ) {
1318 return false;
1319 } else {
1320 return true;
1321 }
1322 } else if ( nCrack == 3 ) { // 3D but is the last crack
1323 if ( this->isStrengthExceeded( newBase, gp, tStep, nCrack, trialStress.at(nCrack) ) ) {
1324 return false;
1325 } else {
1326 return true;
1327 }
1328 } else { // second crack in 3D case
1329 // if the stress-strength criterion is violated
1330 // the crack coordinates have to be changed
1331 // in-plane shear strain remains unchanged
1332
1333 // compute second principal stress cheap - if passed, compute eigenvectors and eigenvalues
1334 sigX = trialStress.at(2); //1st normal stress in crack plane
1335 sigY = trialStress.at(3); //2nd normal stress in crack plane
1336 tau = trialStress.at(4); // shear stress in crack plane
1337
1338 sig2 = ( sigX + sigY ) / 2. + sqrt( ( sigX - sigY ) * ( sigX - sigY ) / 4. + tau * tau );
1339
1340 if ( this->isStrengthExceeded(newBase, gp, tStep, 2, sig2) ) {
1341 oldBase = newBase;
1342 planeStress.resize(3);
1343 planeStress.zero();
1344
1345 planeStress.at(1) = trialStress.at(2); //1st normal stress in crack plane
1346 planeStress.at(2) = trialStress.at(3); //2nd normal stress in crack plane
1347 planeStress.at(3) = trialStress.at(4); //shear in crack plane
1348
1349
1350 // compute principal stresses on the already existing crack plane
1351 this->computePrincipalValDir(princStress, princCrackDir, planeStress, principal_stress);
1352
1353 // right-hand orientation
1354 princCrackDir.at(1, 2) = -1. * princCrackDir.at(2, 1);
1355 princCrackDir.at(2, 2) = princCrackDir.at(1, 1);
1356
1357 // establish vector of crack directions on the crack plane
1358 // the first vector is the normal direction, [1, 0, 0]
1359
1360 princCrackDirExt.resize(3, 3);
1361 princCrackDirExt.zero();
1362 princCrackDirExt.at(1, 1) = 1.;
1363
1364 for ( int i = 1; i <= 2; i++ ) {
1365 for ( int j = 1; j <= 2; j++ ) {
1366 princCrackDirExt.at(i + 1, j + 1) = princCrackDir.at(i, j);
1367 }
1368 }
1369
1370 // new crack base vector in global coordinates
1371 newBase.beProductOf(oldBase, princCrackDirExt);
1372
1373 /*
1374 * // algorithm which gives the same results - rotation of axes in space
1375 * double theta;
1376 *
1377 * if ( princCrackDir.at(1,1) > 0. ) {
1378 * theta = asin( princCrackDir.at(2,1) );
1379 * } else {
1380 * theta = M_PI - asin ( princCrackDir.at(2,1) );
1381 * }
1382 *
1383 * if (theta < 0.) {
1384 * theta += 2 * M_PI;
1385 * }
1386 *
1387 *
1388 * FloatArray baseVector2;
1389 * FloatArray baseVector3;
1390 *
1391 * FloatArray newBaseVector2;
1392 * FloatArray newBaseVector3;
1393 *
1394 * baseVector2.resize(3);
1395 * baseVector3.resize(3);
1396 *
1397 * newBaseVector2.resize(3);
1398 * newBaseVector3.resize(3);
1399 *
1400 * for (int i = 1; i <= 3; i++) {
1401 * baseVector2.at(i) = oldBase.at(i,2);
1402 * baseVector3.at(i) = oldBase.at(i,3);
1403 * }
1404 *
1405 * FloatMatrix K;
1406 * FloatMatrix R;
1407 * FloatMatrix help;
1408 *
1409 * K.resize(3,3);
1410 * K.zero();
1411 *
1412 * K.at(1,2) = -1. * oldBase.at(3,1);
1413 * K.at(1,3) = oldBase.at(2,1);
1414 * K.at(2,1) = oldBase.at(3,1);
1415 * K.at(2,3) = -1. * oldBase.at(1,1);
1416 * K.at(3,1) = -1. * oldBase.at(2,1);
1417 * K.at(3,2) = oldBase.at(1,1);
1418 *
1419 * R.resize(3,3);
1420 * R.beUnitMatrix();
1421 *
1422 * help = K;
1423 * help.times( sin(theta) );
1424 *
1425 * R.add(help);
1426 *
1427 * help.beProductOf (K, K);
1428 * help.times( 1.-cos(theta) );
1429 *
1430 * R.add(help);
1431 *
1432 * newBaseVector2.beProductOf(R, baseVector2);
1433 * newBaseVector3.beProductOf(R, baseVector3);
1434 */
1435
1436
1437 /*
1438 * // the same way - R evaluated symbolically in Matlab
1439 * // c = cos, s = sin
1440 * // k = axis of rotation
1441 * R.at(1,1) = 1. + (c-1.) * k2*k2 + (c-1.) * k3*k3;
1442 * R.at(1,2) = -k3*s - k1*k2 * (c-1.);
1443 * R.at(1,3) = k2*s - k1*k3 * (c-1.);
1444 *
1445 * R.at(2,1) = k3*s - k1*k2 * (c-1.);
1446 * R.at(2,2) = 1. + (c-1.) * k1*k1 + (c-1.) * k3*k3;
1447 * R.at(2,3) = -k1*s - k2*k3 * (c-1.);
1448 *
1449 * R.at(3,1) = -k2*s - k1*k3 * (c-1.);
1450 * R.at(3,2) = k1*s - k2*k3 * (c-1.);
1451 * R.at(3,3) = 1. + (c-1.) * k1*k1 + (c-1.) * k2*k2;
1452 */
1453
1454
1455 // modify shear strains
1456
1457 crackingStrain = status->giveCrackStrainVector();
1458
1459 shearStrains.resize(2);
1460 shearStrains.at(1) = crackingStrain.at(6); //aka x-y
1461 shearStrains.at(2) = crackingStrain.at(5); //aka x-z
1462
1463 // gamma_xy_l [ c s ] gamma_xy_g
1464 // = [ ] *
1465 // gamma_xz_l [ -s c ] gamma_xz_g
1466
1467 newShearStrains.beTProductOf(princCrackDir, shearStrains);
1468
1469 crackingStrain.at(5) = newShearStrains.at(2); // aka x-z
1470 crackingStrain.at(6) = newShearStrains.at(1); // aka x-y
1471
1472 status->setCrackStrainVector(crackingStrain);
1473 status->setTempCrackStrainVector(crackingStrain);
1474
1475 // modify maximum shear strains
1476 // get max shear strains in the original configuration
1477 shearStrains.at(1) = status->giveMaxCrackStrain(6); //aka x-y
1478 shearStrains.at(2) = status->giveMaxCrackStrain(5); //aka x-z
1479
1480 // gamma_xy_l [ c s ] gamma_xy_g
1481 // = [ ] *
1482 // gamma_xz_l [ -s c ] gamma_xz_g
1483
1484 newShearStrains.beTProductOf(princCrackDir, shearStrains);
1485
1486 status->setMaxCrackStrain( 5, max( newShearStrains.at(2), crackingStrain.at(5) ) ); // aka x-z
1487 status->setMaxCrackStrain( 6, max( newShearStrains.at(1), crackingStrain.at(6) ) ); // aka x-y
1488
1489 status->setTempMaxCrackStrain( 5, max( newShearStrains.at(2), crackingStrain.at(5) ) ); // aka x-z
1490 status->setTempMaxCrackStrain( 6, max( newShearStrains.at(1), crackingStrain.at(6) ) ); // aka x-y
1491
1492 return false;
1493 } else { // strength not reached
1494 return true;
1495 }
1496 } // number of crack in given stress-state
1497}
1498
1499
1500double
1501FCMMaterial :: giveCharacteristicElementLength(GaussPoint *gp, const FloatArray &crackPlaneNormal) const
1502{
1503 return gp->giveElement()->giveCharacteristicLength(crackPlaneNormal);
1504}
1505
1506
1507void
1508FCMMaterial :: updateCrackStatus(GaussPoint *gp) const
1509{
1510 FCMMaterialStatus *status = static_cast< FCMMaterialStatus * >( this->giveStatus(gp) );
1511
1512 int nMaxCracks = status->giveMaxNumberOfCracks(gp);
1513
1514 FloatArray crackStrain, maxCrackStrain;
1515 crackStrain = status->giveTempCrackStrainVector();
1516 maxCrackStrain = status->giveMaxCrackStrainVector();
1517
1518
1519 // loop over NORMAL components of the crack strain vector
1520 // statuses are only for normal components
1521 for ( int i = 1; i <= nMaxCracks; i++ ) {
1522 // changes in status only in previously cracked (or at least initiated in this step) material
1523 if ( status->giveTempCrackStatus(i) != pscm_NONE ) {
1524 // temp crack bigger or equal than max strain so far
1525 if ( crackStrain.at(i) >= maxCrackStrain.at(i) ) {
1526 status->setTempMaxCrackStrain( i, crackStrain.at(i) );
1527
1528 if ( status->giveCrackStatus(i) != pscm_NONE ) { //already existing crack in prev step
1529 status->setTempCrackStatus(i, pscm_SOFTENING);
1530 } else if ( ( fabs( crackStrain.at(i) ) <= 1.e-18 ) && ( status->giveCrackStatus(i) == pscm_NONE ) ) { // just initiated crack with zero strain
1531 if ( i == nMaxCracks ) {
1532 status->setTempCrackStatus(i, pscm_NONE);
1533 } else if ( nMaxCracks > i ) { // curent crack is not the last possible one
1534 if ( status->giveTempCrackStatus(i + 1) == pscm_NONE ) { // the next crack is not initiated
1535 status->setTempCrackStatus(i, pscm_NONE); // reset the current one to NONE
1536 } else {
1537 status->setTempCrackStatus(i, pscm_JUST_INIT); // else keep as just initiated
1538 }
1539 } else {
1540 status->setTempCrackStatus(i, pscm_JUST_INIT);
1541 }
1542 } else {
1543 status->setTempCrackStatus(i, pscm_JUST_INIT);
1544 }
1545
1546 // closed previously
1547 } else if ( crackStrain.at(i) <= 0. ) {
1548
1549 if (status->giveCrackStatus(i) != pscm_NONE ) {
1550 status->setTempCrackStatus(i, pscm_CLOSED);
1551 } else {
1552 status->setTempCrackStatus(i, pscm_NONE);
1553 status->setTempMaxCrackStrain(i, 0.);
1554 }
1555
1556
1557 // if ( status->giveMaxCrackStrain(i) == 0. ) { //not existing crack in previous step that wants to close?
1558
1559
1560
1561 // status->setTempMaxCrackStrain(i, 0.);
1562
1563 // status->setTempCrackStatus(i, pscm_NONE);
1564
1565 // status->setTempCrackStatus(i, pscm_NONE);
1566 // CLOSED status should have similar behavior and it is safer. The status can be changed within one iteration loop on the local scale
1567 // the crack can be changed from CLOSED to NONE during overall updating
1568 // status->setTempCrackStatus(i, pscm_CLOSED);
1569 //}
1570
1571 // unloading--reloading
1572 } else if ( crackStrain.at(i) < maxCrackStrain.at(i) ) {
1573 status->setTempCrackStatus(i, pscm_UNLO_RELO);
1574 } else {
1575 OOFEM_ERROR("Unexpected value of %d-th cracking strain %f ", i, crackStrain.at(i) );
1576 }
1577 }
1578 }
1579
1580 // SHEAR components
1581 for ( int i = nMaxCracks + 1; i <= crackStrain.giveSize(); i++ ) {
1582 if ( fabs( crackStrain.at(i) ) > maxCrackStrain.at(i) ) {
1583 status->setTempMaxCrackStrain( i, fabs( crackStrain.at(i) ) );
1584 }
1585 }
1586}
1587
1588
1589
1590
1591
1592void
1593FCMMaterial :: giveTotalLocalCrackedStiffnessMatrix(FloatMatrix &answer,
1594 MatResponseMode rMode, GaussPoint *gp,
1595 TimeStep *tStep) const
1596{
1597 int dim, j;
1598 FloatMatrix Dcr;
1599 IntArray indx;
1600 StructuralMaterial :: giveVoigtSymVectorMask( indx, gp->giveMaterialMode() );
1601
1602 dim = indx.giveSize();
1603
1604 Dcr.resize(dim, dim);
1605
1606 for ( int i = 1; i <= dim; i++ ) {
1607 j = indx.at(i);
1608
1609 if ( j <= 3 ) {
1610 Dcr.at(i, i) = this->giveCrackingModulus(rMode, gp, tStep, j);
1611 } else {
1612 Dcr.at(i, i) = this->computeTotalD2Modulus(gp, tStep, j);
1613 }
1614 }
1615
1616 answer = Dcr;
1617}
1618
1619
1620
1621void
1622FCMMaterial :: giveNormalLocalCrackedStiffnessMatrix(FloatMatrix &answer,
1623 MatResponseMode rMode, GaussPoint *gp,
1624 TimeStep *tStep) const
1625{
1626
1627 FloatMatrix Dcr;
1628 FCMMaterialStatus *status = static_cast< FCMMaterialStatus * >( this->giveStatus(gp) );
1629 int nMaxCr = status->giveMaxNumberOfCracks(gp);
1630
1631 Dcr.resize(nMaxCr, nMaxCr);
1632
1633 for ( int i = 1; i <= nMaxCr; i++ ) {
1634 Dcr.at(i, i) = this->giveCrackingModulus(rMode, gp, tStep, i);
1635 }
1636
1637 answer = Dcr;
1638}
1639
1640
1641
1642void
1643FCMMaterial :: giveMaterialStiffnessMatrix(FloatMatrix &answer,
1644 MatResponseMode rMode,
1645 GaussPoint *gp,
1646 TimeStep *tStep) const
1647//
1648// returns effective material stiffness matrix in full form
1649// for gp stress strain mode
1650//
1651{
1652 FCMMaterialStatus *status = static_cast< FCMMaterialStatus * >( this->giveStatus(gp) );
1653 // StructuralMaterial *lMat = static_cast< StructuralMaterial * >( this->giveLinearElasticMaterial() );
1654
1655 MaterialMode mMode = gp->giveMaterialMode();
1656
1657 int numberOfActiveCracks, nMaxCracks;
1658 double overallElasticStiffness;
1659 double G = this->computeOverallElasticShearModulus(gp, tStep);
1660
1661 FloatMatrix D, De, C, stiffnessL2G;
1662
1663 FloatMatrix DeHelp, Dcr, DcrHelp, DcrHelp2;
1664
1665 numberOfActiveCracks = status->giveNumberOfTempCracks();
1666
1667 // ELASTIC MATRIX
1668 linearElasticMaterial.giveStiffnessMatrix(De, rMode, gp, tStep);
1669 // lMat->giveStiffnessMatrix(De, rMode, gp, tStep);
1670
1671 overallElasticStiffness = this->computeOverallElasticStiffness(gp, tStep);
1672 if ( overallElasticStiffness != ( this->give('E', gp) ) ) {
1673 De.times( overallElasticStiffness / ( this->give('E', gp) ) );
1674 }
1675
1676 if ( ( rMode == ElasticStiffness ) || ( numberOfActiveCracks == 0 ) ) {
1677 answer = De;
1678 return;
1679 }
1680
1681 // SECANT OR TANGENT MATRIX - first in local direction
1682 nMaxCracks = status->giveMaxNumberOfCracks(gp);
1683 De.resizeWithData(nMaxCracks,nMaxCracks);
1684
1685
1686 double E_crack;
1687
1688 if ( rMode == SecantStiffness ) {
1689
1690 C.beInverseOf(De);
1691
1692 for ( int i = 1; i <= numberOfActiveCracks; i++ ) {
1693
1694 E_crack = this->giveCrackingModulus(rMode, gp, tStep, i);
1695
1696 C.at(i, i) += 1. / E_crack;
1697
1698 }
1699
1700
1701 D.beInverseOf(C);
1702
1703 } else if ( rMode == TangentStiffness ) {
1704
1705 // if ( ( rMode == SecantStiffness ) || ( rMode == TangentStiffness ) ) {
1706
1707 // shrink to square matrix number of cracks x number of cracks
1708 DeHelp = De;
1709 DeHelp.resizeWithData(numberOfActiveCracks, numberOfActiveCracks);
1710
1711 Dcr.resize(numberOfActiveCracks, numberOfActiveCracks);
1712 Dcr.zero();
1713
1714 for ( int i = 1; i <= numberOfActiveCracks; i++ ) {
1715 if ( status->giveTempCrackStrain(i) > fcm_THRESHOLD_CRACK_STRAIN ) {
1716 Dcr.at(i, i) = this->giveCrackingModulus(rMode, gp, tStep, i);
1717 } else {
1718 Dcr.at(i, i) = fcm_BIGNUMBER * overallElasticStiffness;
1719 }
1720 }
1721
1722 Dcr.add(DeHelp);
1723 C.beInverseOf(Dcr);
1724 C.resizeWithData(nMaxCracks, nMaxCracks);
1725
1726 DcrHelp.beProductOf(De, C); // De (De + Dcr)^-1
1727 DcrHelp2.beProductOf(DcrHelp, De); // De (De + Dcr)^-1 De
1728
1729 D.zero();
1730 D.resize(nMaxCracks, nMaxCracks);
1731 D.add(De);
1732 D.subtract(DcrHelp2); // De - De (De + Dcr)^-1 De
1733 }
1734
1735
1736 if ( this->normalCoeffNumer > 0. ) {
1737 // apply minimum stiffness
1738 double nu = this->givePoissonsRatio();
1739
1740 for ( int i = 1; i <= nMaxCracks; i++ ) {
1741 for ( int j = 1; j <= nMaxCracks; j++ ) {
1742 if (i == j) { // diagonal terms
1743 D.at(i, j) = max( D.at(i,j), this->normalCoeffNumer * overallElasticStiffness );
1744 } else {
1745 D.at(i, j) = max( D.at(i,j), this->normalCoeffNumer * overallElasticStiffness * nu/(1-nu) );
1746 }
1747 }
1748 }
1749 }
1750
1751
1752
1753
1754
1755 // resize add shear moduli on diagonal
1756 switch ( mMode ) {
1757 case _PlaneStress:
1758 D.resizeWithData(3, 3);
1759
1760 if ( max( status->giveCrackStrain(1), status->giveCrackStrain(2) ) > fcm_THRESHOLD_CRACK_STRAIN ) {
1761 D.at(3, 3) = this->computeEffectiveShearModulus(gp, tStep, 6);
1762 D.at(3, 3) = max( D.at(3,3), this->shearCoeffNumer * G );
1763 } else {
1764 D.at(3, 3) = G;
1765 }
1766
1767 // debug:
1768 /*double help;
1769 help = this->shearCoeffNumer * G;
1770 if (help > D.at(3, 3) ) {
1771 D.at(3, 3) = help;
1772 }*/
1773
1774 break;
1775
1776 case _3dMat:
1777 D.resizeWithData(6, 6);
1778 for ( int i = 4; i <= 6; i++ ) {
1779 if (i == 4) {
1780 if ( max( status->giveCrackStrain(2), status->giveCrackStrain(3) ) > fcm_THRESHOLD_CRACK_STRAIN ) {
1781 D.at(i, i) = this->computeEffectiveShearModulus(gp, tStep, i);
1782 D.at(i, i) = max( D.at(i,i), this->shearCoeffNumer * G );
1783 } else {
1784 D.at(i, i) = G;
1785 }
1786
1787 } else if ( i == 5 ) {
1788 if ( max( status->giveCrackStrain(1), status->giveCrackStrain(3) ) > fcm_THRESHOLD_CRACK_STRAIN ) {
1789 D.at(i, i) = this->computeEffectiveShearModulus(gp, tStep, i);
1790 D.at(i, i) = max( D.at(i,i), this->shearCoeffNumer * G );
1791 } else {
1792 D.at(i, i) = G;
1793 }
1794
1795 } else { // i == 6
1796 if ( max( status->giveCrackStrain(1), status->giveCrackStrain(2) ) > fcm_THRESHOLD_CRACK_STRAIN ) {
1797 D.at(i, i) = this->computeEffectiveShearModulus(gp, tStep, i);
1798 D.at(i, i) = max( D.at(i,i), this->shearCoeffNumer * G );
1799 } else {
1800 D.at(i, i) = G;
1801 }
1802 }
1803 } // end for
1804
1805 break;
1806
1807 case _PlaneStrain:
1808 D.resizeWithData(6, 6);
1809
1810 if ( max( status->giveCrackStrain(1), status->giveCrackStrain(2) ) > fcm_THRESHOLD_CRACK_STRAIN ) {
1811 D.at(6, 6) = this->computeEffectiveShearModulus(gp, tStep, 6);
1812 D.at(6, 6) = max( D.at(6,6), this->shearCoeffNumer * G );
1813 } else {
1814 D.at(6, 6) = G;
1815 }
1816
1817 break;
1818
1819 default:
1820 OOFEM_ERROR( "Material mode %s not supported", __MaterialModeToString(mMode) );
1821 }
1822
1823 // transform stiffnes to global c.s
1824 // to transform stiffness from LOCAL to GLOBAL coordinate system
1825 // the transformation matrix is the same as for strain transformation
1826 // from GLOBAL to LOCAL coordinate system
1827 stiffnessL2G = status->giveG2LStrainVectorTransformationMtrx();
1828 D.rotatedWith(stiffnessL2G, 'n');
1829
1830
1831 switch ( mMode ) {
1832 case _PlaneStress:
1833 answer = D;
1834 return;
1835
1836 break;
1837 case _3dMat:
1838 answer = D;
1839 return;
1840
1841 break;
1842 case _PlaneStrain:
1843 StructuralMaterial :: giveReducedSymMatrixForm( answer, D, gp->giveMaterialMode() );
1844 return;
1845
1846 break;
1847 default:
1848 OOFEM_ERROR( "Material mode %s not supported", __MaterialModeToString(mMode) );
1849 return;
1850 }
1851}
1852
1853
1854
1855double
1856FCMMaterial :: computeTotalD2Modulus(GaussPoint *gp, TimeStep *tStep, int shearDirection) const
1857{
1858 double D2 = 0.;
1859 int crackA, crackB;
1860 double D2_1, D2_2;
1861
1862 if ( this->isIntactForShear(gp, shearDirection) ) {
1863 D2 = this->computeOverallElasticStiffness(gp, tStep) * fcm_BIGNUMBER;
1864 } else {
1865 if ( shearDirection == 4 ) {
1866 crackA = 2;
1867 crackB = 3;
1868 } else if ( shearDirection == 5 ) {
1869 crackA = 1;
1870 crackB = 3;
1871 } else if ( shearDirection == 6 ) {
1872 crackA = 1;
1873 crackB = 2;
1874 } else {
1875 OOFEM_ERROR("Unexpected value of index i (4, 5, 6 permitted only)");
1876 }
1877
1878 if ( ( this->isIntact(gp, crackA) ) || ( this->isIntact(gp, crackB) ) ) {
1879 if ( this->isIntact(gp, crackA) ) {
1880 D2 = this->computeD2ModulusForCrack(gp, tStep, crackB);
1881 } else {
1882 D2 = this->computeD2ModulusForCrack(gp, tStep, crackA);
1883 }
1884 } else {
1885 D2_1 = this->computeD2ModulusForCrack(gp, tStep, crackA);
1886 D2_2 = this->computeD2ModulusForCrack(gp, tStep, crackB);
1887
1888 if ( multipleCrackShear ) {
1889 // serial coupling of stiffnesses 1/D = 1/D1 + 1/D2
1890 D2 = D2_1 * D2_2 / ( D2_1 + D2_2 );
1891 } else {
1892 D2 = min(D2_1, D2_2);
1893 }
1894 }
1895
1896 // adjust stiffness according to the maximum allowable stress
1897 /*FCMMaterialStatus *status = static_cast< FCMMaterialStatus * >( this->giveStatus(gp) );
1898 MaterialMode mMode = gp->giveMaterialMode();
1899 double D2_help, gamma_max;
1900 switch ( mMode ) {
1901 case _PlaneStress:
1902
1903 gamma_max = status->giveMaxCrackStrain(3);
1904 break;
1905
1906 case _PlaneStrain: // check
1907 gamma_max = status->giveMaxCrackStrain(3);
1908 break;
1909
1910 case _3dMat:
1911 gamma_max = status->giveMaxCrackStrain(shearDirection);
1912 break;
1913
1914 default:
1915 OOFEM_ERROR( "Material mode %s not supported", __MaterialModeToString(mMode) );
1916 }
1917 D2_help = this->maxShearStress(gp, tStep, shearDirection) / gamma_max;
1918 D2 = min(D2, D2_help);
1919 */
1920 }
1921
1922 return D2;
1923}
1924
1925
1926double
1927FCMMaterial :: computeNumerD2Modulus(GaussPoint *gp, TimeStep *tStep, int shearDirection) const
1928{
1929 double D2 = 0.;
1930 int crackA, crackB;
1931 double D2_1, D2_2;
1932
1933 if ( this->isIntactForShear(gp, shearDirection) ) {
1934 D2 = this->computeOverallElasticStiffness(gp, tStep) * fcm_BIGNUMBER;
1935 } else {
1936 if ( shearDirection == 4 ) {
1937 crackA = 2;
1938 crackB = 3;
1939 } else if ( shearDirection == 5 ) {
1940 crackA = 1;
1941 crackB = 3;
1942 } else if ( shearDirection == 6 ) {
1943 crackA = 1;
1944 crackB = 2;
1945 } else {
1946 OOFEM_ERROR("Unexpected value of index i (4, 5, 6 permitted only)");
1947 }
1948
1949 if ( ( this->isIntact(gp, crackA) ) || ( this->isIntact(gp, crackB) ) ) {
1950 if ( this->isIntact(gp, crackA) ) {
1951 D2 = this->computeNumerD2ModulusForCrack(gp, tStep, crackB);
1952 } else {
1953 D2 = this->computeNumerD2ModulusForCrack(gp, tStep, crackA);
1954 }
1955 } else {
1956 D2_1 = this->computeNumerD2ModulusForCrack(gp, tStep, crackA);
1957 D2_2 = this->computeNumerD2ModulusForCrack(gp, tStep, crackB);
1958
1959 if ( multipleCrackShear ) {
1960 // serial coupling of stiffnesses 1/D = 1/D1 + 1/D2
1961 D2 = D2_1 * D2_2 / ( D2_1 + D2_2 );
1962 } else {
1963 D2 = min(D2_1, D2_2);
1964 }
1965 }
1966 }
1967
1968 return D2;
1969}
1970
1971
1972
1973void
1974FCMMaterial :: initializeFrom(InputRecord &ir)
1975{
1976 StructuralMaterial :: initializeFrom(ir);
1977 linearElasticMaterial.initializeFrom(ir);
1978
1979 IR_GIVE_OPTIONAL_FIELD(ir,iterLimitGlobal,"iterlimitglobal");
1980
1981 this->nAllowedCracks = 3;
1982 IR_GIVE_OPTIONAL_FIELD(ir, nAllowedCracks, _IFT_FCM_nAllowedCracks);
1983
1984
1985 this->crackSpacing = -1.;
1986 if ( ir.hasField(_IFT_FCM_crackSpacing) ) {
1987 IR_GIVE_FIELD(ir, crackSpacing, _IFT_FCM_crackSpacing);
1988 }
1989
1990 this->multipleCrackShear = false;
1991 if ( ir.hasField(_IFT_FCM_multipleCrackShear) ) {
1992 this->multipleCrackShear = true;
1993 }
1994
1995 int ecsm = 0;
1996 IR_GIVE_OPTIONAL_FIELD(ir, ecsm, _IFT_FCM_ecsm);
1997 switch ( ecsm ) {
1998 case 1: ecsMethod = ECSM_SquareRootOfArea;
1999 break;
2000 case 2: ecsMethod = ECSM_ProjectionCentered;
2001 break;
2002 case 3: ecsMethod = ECSM_Oliver1;
2003 break;
2004 case 4: ecsMethod = ECSM_Oliver1modified;
2005 break;
2006 default: ecsMethod = ECSM_Projection;
2007 }
2008
2009 this->shearCoeffNumer = -1.;
2010 IR_GIVE_OPTIONAL_FIELD(ir, shearCoeffNumer, _IFT_FCM_shearCoeffNumer);
2011 this->normalCoeffNumer = -1. * fcm_BIGNUMBER;
2012 IR_GIVE_OPTIONAL_FIELD(ir, normalCoeffNumer, _IFT_FCM_normalCoeffNumer);
2013}
2014
2015
2016double
2017FCMMaterial :: give(int aProperty, GaussPoint *gp) const
2018{
2019 return linearElasticMaterial.give(aProperty, gp);
2020}
2021
2022
2023
2024double
2025FCMMaterial :: giveCrackSpacing(void) const
2026{
2027 return crackSpacing;
2028}
2029
2030
2031double
2032FCMMaterial :: giveNumberOfCracksInDirection(GaussPoint *gp, int iCrack) const
2033{
2034 FCMMaterialStatus *status = static_cast< FCMMaterialStatus * >( this->giveStatus(gp) );
2035 double spacing, L, N;
2036
2037 L = status->giveCharLength(iCrack);
2038 spacing = this->giveCrackSpacing();
2039
2040 if ( ( spacing > L ) || ( spacing < 0. ) ) {
2041 N = 1;
2042 } else {
2043 N = ( L / spacing );
2044 }
2045
2046 return N;
2047}
2048
2049
2050double
2051FCMMaterial :: giveNumberOfCracksForShearDirection(GaussPoint *gp, int i) const
2052{
2053 double N;
2054 int dir_1, dir_2;
2055
2056 if ( i == 4 ) {
2057 dir_1 = 2;
2058 dir_2 = 3;
2059 } else if ( i == 5 ) {
2060 dir_1 = 1;
2061 dir_2 = 3;
2062 } else if ( i == 6 ) {
2063 dir_1 = 1;
2064 dir_2 = 2;
2065 } else {
2066 OOFEM_ERROR("Unexpected value of index i (4, 5, 6 permitted only)");
2067 }
2068
2069 N = max( this->giveNumberOfCracksInDirection(gp, dir_1), this->giveNumberOfCracksInDirection(gp, dir_2) );
2070
2071 return N;
2072}
2073
2074
2075int
2076FCMMaterial :: giveIPValue(FloatArray &answer, GaussPoint *gp, InternalStateType type, TimeStep *tStep)
2077{
2078 FCMMaterialStatus *status = static_cast< FCMMaterialStatus * >( this->giveStatus(gp) );
2079 double width;
2080 int index;
2081
2082
2083 // first/dominant crack vector
2084 if ( type == IST_CrackVector ) {
2085 answer.resize(3);
2086 answer.zero();
2087
2088 // MAX CRACK
2089 width = 0.;
2090 index = 1;
2091 for ( int i = 1; i <= status->giveNumberOfCracks(); i++ ) {
2092 if ( status->giveCharLength(i) * status->giveCrackStrain(i) / this->giveNumberOfCracksInDirection(gp, i) > width ) {
2093 width = status->giveCharLength(i) * status->giveCrackStrain(i) / this->giveNumberOfCracksInDirection(gp, i);
2094 index = i;
2095 }
2096 }
2097
2098 for ( int i = 1; i <= status->giveCrackDirs().giveNumberOfRows(); i++ ) {
2099 answer.at(i) = status->giveCrackDirs().at(i, index);
2100 }
2101
2102 return 1;
2103 } else if ( type == IST_2ndCrackVector ) {
2104 answer.resize(3);
2105 answer.zero();
2106 index = 2;
2107
2108 if ( status->giveNumberOfCracks() <= 1 ) {
2109 return 1;
2110 } else if ( status->giveNumberOfCracks() == 2 ) {
2111 if ( status->giveCharLength(2) * status->giveCrackStrain(2) / this->giveNumberOfCracksInDirection(gp, 2) > status->giveCharLength(1) * status->giveCrackStrain(1) / this->giveNumberOfCracksInDirection(gp, 1) ) {
2112 index = 1;
2113 }
2114 } else { // 3 cracks
2115 if ( status->giveCharLength(2) * status->giveCrackStrain(2) / this->giveNumberOfCracksInDirection(gp, 2) > status->giveCharLength(1) * status->giveCrackStrain(1) / this->giveNumberOfCracksInDirection(gp, 1) ) { // #2 is biggest
2116 if ( status->giveCharLength(3) * status->giveCrackStrain(3) / this->giveNumberOfCracksInDirection(gp, 3) > status->giveCharLength(1) * status->giveCrackStrain(1) / this->giveNumberOfCracksInDirection(gp, 1) ) {
2117 index = 3;
2118 } else {
2119 index = 1;
2120 }
2121 } else { // #1 is biggest
2122 if ( status->giveCharLength(3) * status->giveCrackStrain(3) / this->giveNumberOfCracksInDirection(gp, 3) > status->giveCharLength(2) * status->giveCrackStrain(2) / this->giveNumberOfCracksInDirection(gp, 2) ) {
2123 index = 3;
2124 } else {
2125 index = 2;
2126 }
2127 }
2128 }
2129
2130 for ( int i = 1; i <= status->giveCrackDirs().giveNumberOfRows(); i++ ) {
2131 answer.at(i) = status->giveCrackDirs().at(i, index);
2132 }
2133
2134 return 1;
2135 } else if ( type == IST_3rdCrackVector ) {
2136 answer.resize(3);
2137 answer.zero();
2138 index = 3;
2139
2140 if ( status->giveNumberOfCracks() <= 2 ) {
2141 return 1;
2142 } else { // 3 cracks
2143 if ( status->giveCharLength(2) * status->giveCrackStrain(2) / this->giveNumberOfCracksInDirection(gp, 2) > status->giveCharLength(1) * status->giveCrackStrain(1) / this->giveNumberOfCracksInDirection(gp, 1) ) { // #2 is biggest
2144 if ( status->giveCharLength(3) * status->giveCrackStrain(3) / this->giveNumberOfCracksInDirection(gp, 3) > status->giveCharLength(1) * status->giveCrackStrain(1) / this->giveNumberOfCracksInDirection(gp, 1) ) {
2145 index = 1;
2146 } else {
2147 index = 3;
2148 }
2149 } else { // #1 is biggest
2150 if ( status->giveCharLength(3) * status->giveCrackStrain(3) / this->giveNumberOfCracksInDirection(gp, 3) > status->giveCharLength(2) * status->giveCrackStrain(2) / this->giveNumberOfCracksInDirection(gp, 2) ) {
2151 index = 2;
2152 } else {
2153 index = 3;
2154 }
2155 }
2156 }
2157
2158 for ( int i = 1; i <= status->giveCrackDirs().giveNumberOfRows(); i++ ) {
2159 answer.at(i) = status->giveCrackDirs().at(i, index);
2160 }
2161
2162 return 1;
2163
2164
2165 // width of a first/dominant crack
2166 } else if ( type == IST_CrackWidth ) {
2167 answer.resize(1);
2168 answer.zero();
2169
2170 // MAX WIDTH
2171 width = 0.;
2172 for ( int i = 1; i <= status->giveNumberOfCracks(); i++ ) {
2173 width = max( width, status->giveCharLength(i) * status->giveCrackStrain(i) / this->giveNumberOfCracksInDirection(gp, i) );
2174 }
2175 answer.at(1) = width;
2176 return 1;
2177 } else if ( type == IST_2ndCrackWidth ) {
2178 answer.resize(1);
2179 answer.zero();
2180 index = 2;
2181 if ( status->giveNumberOfCracks() <= 1 ) {
2182 return 1;
2183 } else if ( status->giveNumberOfCracks() == 2 ) {
2184 if ( status->giveCharLength(2) * status->giveCrackStrain(2) / this->giveNumberOfCracksInDirection(gp, 2) > status->giveCharLength(1) * status->giveCrackStrain(1) / this->giveNumberOfCracksInDirection(gp, 1) ) {
2185 index = 1;
2186 }
2187 } else { // 3 cracks
2188 if ( status->giveCharLength(2) * status->giveCrackStrain(2) / this->giveNumberOfCracksInDirection(gp, 2) > status->giveCharLength(1) * status->giveCrackStrain(1) / this->giveNumberOfCracksInDirection(gp, 1) ) { // #2 is biggest
2189 if ( status->giveCharLength(3) * status->giveCrackStrain(3) / this->giveNumberOfCracksInDirection(gp, 3) > status->giveCharLength(1) * status->giveCrackStrain(1) / this->giveNumberOfCracksInDirection(gp, 1) ) {
2190 index = 3;
2191 } else {
2192 index = 1;
2193 }
2194 } else { // #1 is biggest
2195 if ( status->giveCharLength(3) * status->giveCrackStrain(3) / this->giveNumberOfCracksInDirection(gp, 3) > status->giveCharLength(2) * status->giveCrackStrain(2) / this->giveNumberOfCracksInDirection(gp, 2) ) {
2196 index = 3;
2197 } else {
2198 index = 2;
2199 }
2200 }
2201 }
2202
2203 width = status->giveCharLength(index) * status->giveCrackStrain(index) / this->giveNumberOfCracksInDirection(gp, index);
2204
2205 answer.at(1) = width;
2206 return 1;
2207 } else if ( type == IST_CrackDirs ) {
2208 const FloatMatrix &help = status->giveCrackDirs();
2209 answer.resize(9);
2210 for ( int i = 1; i <= 3; i++ ) {
2211 answer.at(i) = help.at(1, i);
2212 answer.at(i + 3) = help.at(2, i);
2213 answer.at(i + 6) = help.at(3, i);
2214 }
2215
2216 return 1;
2217
2218 } else if ( type == IST_CrackStatuses ) {
2219 answer.resize(3);
2220 answer.zero();
2221 for ( int i = 1; i <= status->giveNumberOfCracks(); i++ ) {
2222 answer.at(i) = status->giveCrackStatus(i);
2223 }
2224 return 1;
2225
2226 } else if ( type == IST_CrackStatusesTemp ) {
2227 answer.resize(3);
2228 answer.zero();
2229 for ( int i = 1; i <= status->giveNumberOfTempCracks(); i++ ) {
2230 answer.at(i) = status->giveTempCrackStatus(i);
2231 }
2232 return 1;
2233
2234 } else if ( type == IST_CrackStrainTensor ) {
2235 // FloatArray crackStrain = status->giveCrackStrainVector();
2236 FloatArray crackStrain;
2237 FloatMatrix epsL2G = status->giveL2GStrainVectorTransformationMtrx();
2238 // from local to global
2239 // crackStrain.rotatedWith(epsL2G, 'n');
2240 crackStrain.beProductOf(epsL2G, status->giveCrackStrainVector() );
2241
2242 StructuralMaterial :: giveFullSymVectorForm( answer, crackStrain, gp->giveMaterialMode() );
2243
2244 return 1;
2245
2246 } else if ( type == IST_CrackSlip ) {
2247 answer.resize(1);
2248 answer.zero();
2249
2250 /* if ( ( mMode == _PlaneStress ) || ( mMode == _PlaneStrain ) ) {
2251 answer = this->computeShearSlipOnCrack(gp, tStep, 1);
2252 } else {*/
2253
2254 int nMaxCracks = status->giveNumberOfTempCracks();
2255 for ( int i = 1; i <= nMaxCracks; i++ ) {
2256 answer.at(1) = max(answer.at(1), this->computeShearSlipOnCrack(gp, tStep, i) );
2257 }
2258 // }
2259 return 1;
2260
2261 } else {
2262 return StructuralMaterial :: giveIPValue(answer, gp, type, tStep);
2263 }
2264}
2265
2266FloatMatrixF<6,6>
2267FCMMaterial :: give3dMaterialStiffnessMatrix(MatResponseMode mode,
2268 GaussPoint *gp,
2269 TimeStep *tStep) const
2270{
2271 FloatMatrix answer;
2272 const_cast<FCMMaterial*>(this)->giveMaterialStiffnessMatrix(answer, mode, gp, tStep);
2273 return answer;
2274}
2275
2276
2277FloatMatrixF<3,3>
2278FCMMaterial :: givePlaneStressStiffMtrx(MatResponseMode mode,
2279 GaussPoint *gp,
2280 TimeStep *tStep) const
2281{
2282 FloatMatrix answer;
2283 const_cast<FCMMaterial*>(this)->giveMaterialStiffnessMatrix(answer, mode, gp, tStep);
2284 return answer;
2285}
2286
2287
2288FloatMatrixF<4,4>
2289FCMMaterial :: givePlaneStrainStiffMtrx(MatResponseMode mode,
2290 GaussPoint *gp,
2291 TimeStep *tStep) const
2292
2293{
2294 FloatMatrix answer;
2295 const_cast<FCMMaterial*>(this)->giveMaterialStiffnessMatrix(answer, mode, gp, tStep);
2296 return answer;
2297}
2298
2299
2300FCMMaterialStatus :: FCMMaterialStatus(GaussPoint *gp) :
2301 StructuralMaterialStatus(gp),
2302 crackStatuses(), tempCrackStatuses(),
2303 maxCrackStrains(), tempMaxCrackStrains(),
2304 crackStrainVector(), tempCrackStrainVector(),
2305 crackDirs(),
2306 charLengths(),
2307 transMatrix_G2Lstress(), transMatrix_G2Lstrain(),
2308 transMatrix_L2Gstress(), transMatrix_L2Gstrain()
2309{
2310 // resize in constructor according to stress-state
2311 this->nMaxCracks = giveMaxNumberOfCracks(gp);
2312
2313 crackStatuses.resize(this->nMaxCracks);
2314 crackStatuses.zero();
2315 tempCrackStatuses = crackStatuses;
2316
2317 charLengths.resize(this->nMaxCracks);
2318 charLengths.zero();
2319
2320 crackDirs.resize(this->nMaxCracks, this->nMaxCracks);
2321 crackDirs.zero();
2322 for ( int i = 1; i <= this->nMaxCracks; i++ ) {
2323 crackDirs.at(i, i) = 1.0;
2324 }
2325
2326
2327 if ( this->nMaxCracks == 2 ) { //plane stress
2328 maxCrackStrains.resize(3);
2329 maxCrackStrains.zero();
2330 tempMaxCrackStrains = maxCrackStrains;
2331
2332 crackStrainVector = maxCrackStrains;
2333 tempCrackStrainVector = maxCrackStrains;
2334
2335
2336 transMatrix_G2Lstress.resize(3, 3);
2337 transMatrix_G2Lstress.zero();
2338 transMatrix_L2Gstrain = transMatrix_L2Gstress = transMatrix_G2Lstrain = transMatrix_G2Lstress;
2339 } else {
2340 maxCrackStrains.resize(6);
2341 maxCrackStrains.zero();
2342 tempMaxCrackStrains = maxCrackStrains;
2343
2344 crackStrainVector = maxCrackStrains;
2345 tempCrackStrainVector = maxCrackStrains;
2346
2347
2348 transMatrix_G2Lstress.resize(6, 6);
2349 transMatrix_G2Lstress.zero();
2350 transMatrix_L2Gstrain = transMatrix_L2Gstress = transMatrix_G2Lstrain = transMatrix_G2Lstress;
2351 }
2352}
2353
2354
2355void
2356FCMMaterialStatus :: printOutputAt(FILE *file, TimeStep *tStep) const
2357{
2358 char s [ 11 ];
2359
2360 StructuralMaterialStatus :: printOutputAt(file, tStep);
2361 fprintf(file, "status { ");
2362 if ( this->giveNumberOfCracks() > 0 ) {
2363 for ( int i = 1; i <= crackDirs.giveNumberOfColumns(); i++ ) {
2364 switch ( crackStatuses.at(i) ) {
2365 case pscm_NONE:
2366 strcpy(s, "NONE");
2367 break;
2368 case pscm_JUST_INIT:
2369 strcpy(s, "JUST_INIT");
2370 break;
2371 case pscm_SOFTENING:
2372 strcpy(s, "SOFTENING");
2373 break;
2374 case pscm_UNLO_RELO:
2375 strcpy(s, "UNLO_RELO");
2376 break;
2377 case pscm_CLOSED:
2378 strcpy(s, "CLOSED");
2379 break;
2380 default:
2381 strcpy(s, "UNKNOWN");
2382 break;
2383 }
2384
2385 fprintf(file, "crack %d {status %s, crackplane is normal to { ", i, s);
2386
2387 for ( int j = 1; j <= crackDirs.giveNumberOfRows(); j++ ) {
2388 fprintf( file, "%f ", crackDirs.at(j, i) );
2389 }
2390
2391 fprintf(file, "} ");
2392
2393 fprintf(file, "crack_width %.4e", this->giveCharLength(i) * this->giveCrackStrain(i) );
2394
2395 fprintf(file, " } ");
2396 }
2397 }
2398
2399 fprintf(file, " }\n");
2400}
2401
2402int
2403FCMMaterialStatus :: giveNumberOfCracks() const
2404//
2405// return number of existing cracks
2406//
2407{
2408 int answer = 0;
2409
2410 for ( int i = 1; i <= crackStatuses.giveSize(); i++ ) {
2411 if ( crackStatuses.at(i) != pscm_NONE ) {
2412 answer++;
2413 }
2414 }
2415
2416 return answer;
2417}
2418
2419
2420int
2421FCMMaterialStatus :: giveNumberOfTempCracks() const
2422//
2423// return number of existing temp cracks
2424//
2425{
2426 int answer = 0;
2427
2428 for ( int i = 1; i <= tempCrackStatuses.giveSize(); i++ ) {
2429 if ( tempCrackStatuses.at(i) != pscm_NONE ) {
2430 answer++;
2431 }
2432 }
2433
2434 return answer;
2435}
2436
2437int
2438FCMMaterialStatus :: giveMaxNumberOfCracks(GaussPoint *gp)
2439//
2440// return number of maximum allowable cracks
2441//
2442{
2443 int nCr = 0;
2444 IntArray indx;
2445
2446 if ( this->nMaxCracks == 0 ) {
2447 StructuralMaterial :: giveVoigtSymVectorMask( indx, gp->giveMaterialMode() );
2448
2449 for ( int i = 1; i <= 3; i++ ) {
2450 if ( indx.contains(i) ) {
2451 nCr++;
2452 }
2453 }
2454
2455 this->nMaxCracks = nCr;
2456 }
2457
2458 return this->nMaxCracks;
2459}
2460
2461
2462void
2463FCMMaterialStatus :: setTempNormalCrackStrainVector(FloatArray tempNormalCrackStrain)
2464//
2465// return number of maximum allowable cracks
2466//
2467{
2468
2469 for ( int i = 1; i <= tempNormalCrackStrain.giveSize(); i++ ) {
2470 tempCrackStrainVector.at(i) = tempNormalCrackStrain.at(i);
2471 }
2472
2473}
2474
2475
2476void
2477FCMMaterialStatus :: initTempStatus()
2478//
2479// initializes temp variables according to variables form previous equlibrium state.
2480//
2481{
2482 StructuralMaterialStatus :: initTempStatus();
2483
2484 tempCrackStatuses = crackStatuses;
2485 tempMaxCrackStrains = maxCrackStrains;
2486 tempCrackStrainVector = crackStrainVector;
2487}
2488
2489
2490
2491void
2492FCMMaterialStatus :: updateYourself(TimeStep *tStep)
2493//
2494// updates variables (nonTemp variables describing situation at previous equilibrium state)
2495// after a new equilibrium state has been reached
2496// temporary variables correspond to newly reched equilibrium.
2497//
2498{
2499 StructuralMaterialStatus :: updateYourself(tStep);
2500
2501
2502 maxCrackStrains = tempMaxCrackStrains;
2503 crackStrainVector = tempCrackStrainVector;
2504
2505 // for ( int i = 1; i <= crackStrainVector.giveSize(); i++ ) {
2506 for ( int i = 1; i <= this->nMaxCracks; i++ ) { // loop only in normal directions!
2507 if ( crackStrainVector.at(i) < 0. ) {
2508 crackStrainVector.at(i) = 0.;
2509 }
2510 }
2511
2512 // updating of statuses has to be done more carefully.
2513 // consider a crack which does not exist in the previous step and ends as closed in the end of this step
2514 // this crack is naturally treated as "NONE" in the following steps
2515
2516 for ( int i = 1; i <= crackStatuses.giveSize(); i++ ) {
2517 if ( ( tempCrackStatuses.at(i) == pscm_CLOSED ) && ( crackStatuses.at(i) == pscm_NONE ) ) {
2518 // no other crack so this one can be set as non-existing
2519 if ( i + 1 > nMaxCracks ) {
2520 crackStatuses.at(i) = pscm_NONE;
2521
2522
2523 // be sure that in the second and third crack does not exist, if it does we have to copy CLOSED status
2524 } else if ( tempCrackStatuses.at(i + 1) != pscm_NONE ) {
2525 crackStatuses.at(i) = tempCrackStatuses.at(i);
2526 } else {
2527 crackStatuses.at(i) = pscm_NONE;
2528 }
2529 } else {
2530 crackStatuses.at(i) = tempCrackStatuses.at(i);
2531 }
2532 }
2533}
2534
2535
2536
2537void
2538FCMMaterialStatus :: saveContext(DataStream &stream, ContextMode mode)
2539{
2540 StructuralMaterialStatus :: saveContext(stream, mode);
2541
2542 contextIOResultType iores;
2543 if ( ( iores = crackStatuses.storeYourself(stream) ) != CIO_OK ) {
2544 THROW_CIOERR(iores);
2545 }
2546
2547 if ( ( iores = maxCrackStrains.storeYourself(stream) ) != CIO_OK ) {
2548 THROW_CIOERR(iores);
2549 }
2550
2551 if ( ( iores = crackDirs.storeYourself(stream) ) != CIO_OK ) {
2552 THROW_CIOERR(iores);
2553 }
2554
2555 if ( ( iores = charLengths.storeYourself(stream) ) != CIO_OK ) {
2556 THROW_CIOERR(iores);
2557 }
2558
2559 if ( ( iores = crackStrainVector.storeYourself(stream) ) != CIO_OK ) {
2560 THROW_CIOERR(iores);
2561 }
2562
2563 if ( ( iores = crackStrainVector.storeYourself(stream) ) != CIO_OK ) {
2564 THROW_CIOERR(iores);
2565 }
2566
2567 if ( ( iores = transMatrix_G2Lstrain.storeYourself(stream) ) != CIO_OK ) {
2568 THROW_CIOERR(iores);
2569 }
2570
2571 if ( ( iores = transMatrix_G2Lstress.storeYourself(stream) ) != CIO_OK ) {
2572 THROW_CIOERR(iores);
2573 }
2574
2575 if ( ( iores = transMatrix_L2Gstrain.storeYourself(stream) ) != CIO_OK ) {
2576 THROW_CIOERR(iores);
2577 }
2578
2579 if ( ( iores = transMatrix_L2Gstress.storeYourself(stream) ) != CIO_OK ) {
2580 THROW_CIOERR(iores);
2581 }
2582}
2583
2584void
2585FCMMaterialStatus :: restoreContext(DataStream &stream, ContextMode mode)
2586{
2587 StructuralMaterialStatus :: restoreContext(stream, mode);
2588
2589 contextIOResultType iores;
2590 if ( ( iores = crackStatuses.restoreYourself(stream) ) != CIO_OK ) {
2591 THROW_CIOERR(iores);
2592 }
2593
2594 if ( ( iores = maxCrackStrains.restoreYourself(stream) ) != CIO_OK ) {
2595 THROW_CIOERR(iores);
2596 }
2597
2598 if ( ( iores = crackDirs.restoreYourself(stream) ) != CIO_OK ) {
2599 THROW_CIOERR(iores);
2600 }
2601
2602 if ( ( iores = charLengths.restoreYourself(stream) ) != CIO_OK ) {
2603 THROW_CIOERR(iores);
2604 }
2605
2606 if ( ( iores = crackStrainVector.restoreYourself(stream) ) != CIO_OK ) {
2607 THROW_CIOERR(iores);
2608 }
2609
2610 if ( ( iores = transMatrix_G2Lstrain.restoreYourself(stream) ) != CIO_OK ) {
2611 THROW_CIOERR(iores);
2612 }
2613
2614 if ( ( iores = transMatrix_G2Lstress.restoreYourself(stream) ) != CIO_OK ) {
2615 THROW_CIOERR(iores);
2616 }
2617
2618 if ( ( iores = transMatrix_L2Gstrain.restoreYourself(stream) ) != CIO_OK ) {
2619 THROW_CIOERR(iores);
2620 }
2621
2622 if ( ( iores = transMatrix_L2Gstress.restoreYourself(stream) ) != CIO_OK ) {
2623 THROW_CIOERR(iores);
2624 }
2625}
2626} // end namespace oofem
N
#define N(a, b)
oofem::Element::giveCharacteristicLength
virtual double giveCharacteristicLength(const FloatArray &normalToCrackPlane)
Definition element.h:938
oofem::FCMMaterialStatus::setCharLength
void setCharLength(int icrack, double val)
sets characteristic length for i-th crack
Definition fcm.h:170
oofem::FCMMaterialStatus::setTempMaxCrackStrain
void setTempMaxCrackStrain(int icrack, double val)
sets value of the maximum crack strain for the i-th crack (temporary value)
Definition fcm.h:116
oofem::FCMMaterialStatus::transMatrix_G2Lstress
FloatMatrix transMatrix_G2Lstress
transformation matrix converting stress from global to local coordinate system
Definition fcm.h:86
oofem::FCMMaterialStatus::setL2GStressVectorTransformationMtrx
void setL2GStressVectorTransformationMtrx(FloatMatrix t)
sets transformation matrix for stress transformation from local to global coordinate system
Definition fcm.h:148
oofem::FCMMaterialStatus::giveG2LStressVectorTransformationMtrx
const FloatMatrix & giveG2LStressVectorTransformationMtrx()
returns transformation matrix for stress transformation from global to local coordinate system
Definition fcm.h:153
oofem::FCMMaterialStatus::giveL2GStrainVectorTransformationMtrx
const FloatMatrix & giveL2GStrainVectorTransformationMtrx()
sets transformation matrix for stress transformation from global to local coordinate system
Definition fcm.h:159
oofem::FCMMaterialStatus::setTempNormalCrackStrainVector
void setTempNormalCrackStrainVector(FloatArray a)
sets temporary vector of cracking strains (normal components)
Definition fcm.C:2463
oofem::FCMMaterialStatus::tempMaxCrackStrains
FloatArray tempMaxCrackStrains
Definition fcm.h:77
oofem::FCMMaterialStatus::transMatrix_L2Gstress
FloatMatrix transMatrix_L2Gstress
transformation matrix converting stress from local to global coordinate system
Definition fcm.h:90
oofem::FCMMaterialStatus::crackStatuses
IntArray crackStatuses
crack statuses (none, just initialized, softening, unload-reload, closed)
Definition fcm.h:75
oofem::FCMMaterialStatus::giveCharLength
double giveCharLength(int icrack) const
returns characteristic length associated with i-th crack direction
Definition fcm.h:162
oofem::FCMMaterialStatus::giveCrackStrain
double giveCrackStrain(int icrack) const
returns i-th component of the crack strain vector (equilibrated)
Definition fcm.h:132
oofem::FCMMaterialStatus::setMaxCrackStrain
void setMaxCrackStrain(int icrack, double val)
sets value of the maximum crack strain for the i-th crack (equilibrated value)
Definition fcm.h:111
oofem::FCMMaterialStatus::setCrackStrainVector
void setCrackStrainVector(FloatArray a)
sets equilibrated vector of cracking strains (max 6 components)
Definition fcm.h:142
oofem::FCMMaterialStatus::charLengths
FloatArray charLengths
Characteristic lengths computed from the crack orientation and element geometry.
Definition fcm.h:83
oofem::FCMMaterialStatus::crackDirs
FloatMatrix crackDirs
Storing direction of cracks (crack normals) in columwise format.
Definition fcm.h:81
oofem::FCMMaterialStatus::giveTempCrackStatus
const IntArray & giveTempCrackStatus()
returns vector of temporary crack statuses
Definition fcm.h:119
oofem::FCMMaterialStatus::tempCrackStrainVector
FloatArray tempCrackStrainVector
Definition fcm.h:79
oofem::FCMMaterialStatus::setTempCrackStrain
void setTempCrackStrain(int icrack, double val)
sets temporary value of i-th cracking strain (max 6 components)
Definition fcm.h:140
oofem::FCMMaterialStatus::giveG2LStrainVectorTransformationMtrx
const FloatMatrix & giveG2LStrainVectorTransformationMtrx()
sets transformation matrix for strain transformation from global to local coordinate system
Definition fcm.h:155
oofem::FCMMaterialStatus::giveMaxCrackStrainVector
const FloatArray & giveMaxCrackStrainVector()
returns vector with maximum cracking strains (max 3 components)
Definition fcm.h:107
oofem::FCMMaterialStatus::giveL2GStressVectorTransformationMtrx
const FloatMatrix & giveL2GStressVectorTransformationMtrx()
sets transformation matrix for stress transformation from local to global coordinate system
Definition fcm.h:157
oofem::FCMMaterialStatus::giveCrackStatus
int giveCrackStatus(int icrack) const
return equilibrated value of status associated with i-th crack direction
Definition fcm.h:125
oofem::FCMMaterialStatus::maxCrackStrains
FloatArray maxCrackStrains
Max. crack strain reached in the entire previous history.
Definition fcm.h:77
oofem::FCMMaterialStatus::giveNumberOfCracks
virtual int giveNumberOfCracks() const
returns number of cracks from the previous time step (equilibrated value)
Definition fcm.C:2403
oofem::FCMMaterialStatus::setG2LStrainVectorTransformationMtrx
void setG2LStrainVectorTransformationMtrx(FloatMatrix s)
sets transformation matrix for strain transformation from global to local coordinate system
Definition fcm.h:146
oofem::FCMMaterialStatus::setTempCrackStrainVector
void setTempCrackStrainVector(FloatArray a)
sets temporary vector of cracking strains (max 6 components)
Definition fcm.h:136
oofem::FCMMaterialStatus::setTempCrackStatus
void setTempCrackStatus(int icrack, int val)
sets temporary value of status for of the i-th crack
Definition fcm.h:123
oofem::FCMMaterialStatus::setCrackDirs
void setCrackDirs(FloatMatrix a)
sets matrix with crack directions (normal vectors)
Definition fcm.h:176
oofem::FCMMaterialStatus::tempCrackStatuses
IntArray tempCrackStatuses
Definition fcm.h:75
oofem::FCMMaterialStatus::giveNumberOfTempCracks
virtual int giveNumberOfTempCracks() const
returns temporary number of cracks
Definition fcm.C:2421
oofem::FCMMaterialStatus::transMatrix_L2Gstrain
FloatMatrix transMatrix_L2Gstrain
transformation matrix converting strain from local to global coordinate system
Definition fcm.h:92
oofem::FCMMaterialStatus::setG2LStressVectorTransformationMtrx
void setG2LStressVectorTransformationMtrx(FloatMatrix t)
sets transformation matrix for stress transformation from global to local coordinate system
Definition fcm.h:144
oofem::FCMMaterialStatus::crackStrainVector
FloatArray crackStrainVector
Components of crack strain vector (normal as well as shear).
Definition fcm.h:79
oofem::FCMMaterialStatus::giveMaxCrackStrain
double giveMaxCrackStrain(int icrack)
returns maximum crack strain for the i-th crack (equilibrated value)
Definition fcm.h:109
oofem::FCMMaterialStatus::giveTempCrackStrain
double giveTempCrackStrain(int icrack) const
returns i-th component of the crack strain vector (temporary)
Definition fcm.h:134
oofem::FCMMaterialStatus::setL2GStrainVectorTransformationMtrx
void setL2GStrainVectorTransformationMtrx(FloatMatrix s)
sets transformation matrix for stress transformation from global to local coordinate system
Definition fcm.h:150
oofem::FCMMaterialStatus::giveCrackDirs
const FloatMatrix & giveCrackDirs()
returns crack directions
Definition fcm.h:172
oofem::FCMMaterialStatus::nMaxCracks
int nMaxCracks
number of maximum possible cracks (optional parameter)
Definition fcm.h:95
oofem::FCMMaterialStatus::giveTempCrackStrainVector
const FloatArray & giveTempCrackStrainVector()
return temporary crack strain vector (max 6 components)
Definition fcm.h:130
oofem::FCMMaterialStatus::giveMaxNumberOfCracks
virtual int giveMaxNumberOfCracks(GaussPoint *gp)
returns maximum number of cracks associated with current mode
Definition fcm.C:2438
oofem::FCMMaterialStatus::transMatrix_G2Lstrain
FloatMatrix transMatrix_G2Lstrain
transformation matrix converting strain from global to local coordinate system
Definition fcm.h:88
oofem::FCMMaterialStatus::giveCrackStrainVector
const FloatArray & giveCrackStrainVector() const
return equilibrated crack strain vector (max 6 components)
Definition fcm.h:128
oofem::FCMMaterial::givePoissonsRatio
virtual double givePoissonsRatio() const
returns Poisson's ratio
Definition fcm.h:383
oofem::FCMMaterial::giveCharacteristicElementLength
virtual double giveCharacteristicElementLength(GaussPoint *gp, const FloatArray &crackPlaneNormal) const
returns characteristic element length in given direction
Definition fcm.C:1501
oofem::FCMMaterial::updateCrackStatus
virtual void updateCrackStatus(GaussPoint *gp) const
updates crack statuses
Definition fcm.C:1508
oofem::FCMMaterial::multipleCrackShear
bool multipleCrackShear
if true = takes shear compliance of all cracks, false = only dominant crack contribution,...
Definition fcm.h:281
oofem::FCMMaterial::giveTensileStrength
virtual double giveTensileStrength(GaussPoint *gp, TimeStep *tStep) const =0
comutes tensile strength
oofem::FCMMaterial::checkSnapBack
virtual void checkSnapBack(GaussPoint *gp, TimeStep *tStep, int crack) const =0
checks possible snap-back
oofem::FCMMaterial::shearCoeffNumer
double shearCoeffNumer
minimum ratio of effective shear modulus vs. shear modulus - just for numerical purpose
Definition fcm.h:346
oofem::FCMMaterial::isIntactForShear
virtual bool isIntactForShear(GaussPoint *gp, int i) const
returns true for closed or no cracks associated to given shear direction (i = 4, 5,...
Definition fcm.C:1060
oofem::FCMMaterial::iterLimitGlobal
int iterLimitGlobal
Definition fcm.h:283
oofem::FCMMaterial::giveMaterialStiffnessMatrix
virtual void giveMaterialStiffnessMatrix(FloatMatrix &answer, MatResponseMode, GaussPoint *gp, TimeStep *tStep) const
Definition fcm.C:1643
oofem::FCMMaterial::computeTotalD2Modulus
virtual double computeTotalD2Modulus(GaussPoint *gp, TimeStep *tStep, int i) const
shear modulus for a given shear direction (4, 5, 6)
Definition fcm.C:1856
oofem::FCMMaterial::giveTotalLocalCrackedStiffnessMatrix
virtual void giveTotalLocalCrackedStiffnessMatrix(FloatMatrix &answer, MatResponseMode rMode, GaussPoint *gp, TimeStep *tStep) const
returns local stiffness matrix in the cracks' direction (total according to the material mode)
Definition fcm.C:1593
oofem::FCMMaterial::isIntact
virtual bool isIntact(GaussPoint *gp, int icrack) const
returns true for closed or no crack (i = 1, 2, 3)
Definition fcm.C:1042
oofem::FCMMaterial::giveCrackSpacing
virtual double giveCrackSpacing(void) const
returns either user-provided value of crack spacing or a value computed from composition
Definition fcm.C:2025
oofem::FCMMaterial::giveNumberOfCracksInDirection
virtual double giveNumberOfCracksInDirection(GaussPoint *gp, int iCrack) const
returns number of fictiotious parallel cracks in the direction of i-th crack
Definition fcm.C:2032
oofem::FCMMaterial::computeEffectiveShearModulus
virtual double computeEffectiveShearModulus(GaussPoint *gp, TimeStep *tStep, int i) const =0
returns Geff which is necessary in the global stiffness matrix
oofem::FCMMaterial::give
double give(int aProperty, GaussPoint *gp) const override
Definition fcm.C:2017
oofem::FCMMaterial::computeOverallElasticShearModulus
virtual double computeOverallElasticShearModulus(GaussPoint *gp, TimeStep *tStep) const
returns overall shear modulus
Definition fcm.h:380
oofem::FCMMaterial::computeOverallElasticStiffness
virtual double computeOverallElasticStiffness(GaussPoint *gp, TimeStep *tStep) const
returns overall Young's modulus
Definition fcm.h:377
oofem::FCMMaterial::normalCoeffNumer
double normalCoeffNumer
minimum ratio of effective normal stiffness vs. overall modulus - just for numerical purpose
Definition fcm.h:349
oofem::FCMMaterial::giveNormalLocalCrackedStiffnessMatrix
virtual void giveNormalLocalCrackedStiffnessMatrix(FloatMatrix &answer, MatResponseMode rMode, GaussPoint *gp, TimeStep *tStep) const
returns local stiffness matrix in the cracks' direction (only normal components)
Definition fcm.C:1622
oofem::FCMMaterial::checkStrengthCriterion
virtual bool checkStrengthCriterion(FloatMatrix &newBase, const FloatArray &globalStress, GaussPoint *gp, TimeStep *tStep, int nCrack) const
checks if the globalStress does not exceed strength in the direction of newBase for n-th crack
Definition fcm.C:1301
oofem::FCMMaterial::isStrengthExceeded
virtual bool isStrengthExceeded(const FloatMatrix &base, GaussPoint *gp, TimeStep *tStep, int iCrack, double trialStress) const
compares trial stress with strength. Returns true if the strength is exceeded. Function oveloaded in ...
Definition fcm.C:1275
oofem::FCMMaterial::maxShearStress
virtual double maxShearStress(GaussPoint *gp, TimeStep *tStep, int i) const =0
computes the maximum value of the shear stress; if the shear stress exceeds this value,...
oofem::FCMMaterial::isThisShearComponent
virtual bool isThisShearComponent(GaussPoint *gp, int component) const
returns true if current component is associated with shear
Definition fcm.C:1094
oofem::FCMMaterial::computeD2ModulusForCrack
virtual double computeD2ModulusForCrack(GaussPoint *gp, TimeStep *tStep, int icrack) const =0
shear modulus for a given crack plane (1, 2, 3)
oofem::FCMMaterial::linearElasticMaterial
IsotropicLinearElasticMaterial linearElasticMaterial
Definition fcm.h:199
oofem::FCMMaterial::initializeCrack
virtual void initializeCrack(GaussPoint *gp, TimeStep *tStep, FloatMatrix &base, int nCrack) const
Definition fcm.C:972
oofem::FCMMaterial::computeShearStiffnessRedistributionFactor
virtual double computeShearStiffnessRedistributionFactor(GaussPoint *gp, TimeStep *tStep, int ithCrackPlane, int jthCrackDirection) const
function calculating ratio used to split shear slips on two crack planes
Definition fcm.C:1287
oofem::FCMMaterial::ecsMethod
ElementCharSizeMethod ecsMethod
Method used for evaluation of characteristic element size.
Definition fcm.h:278
oofem::FCMMaterial::giveCrackingModulus
virtual double giveCrackingModulus(MatResponseMode rMode, GaussPoint *gp, TimeStep *tStep, int i) const =0
returns stiffness in the normal direction of the i-th crack
oofem::FCMMaterial::FCMMaterial
FCMMaterial(int n, Domain *d)
Definition fcm.C:46
oofem::FCMMaterial::crackSpacing
double crackSpacing
value of crack spacing (allows to "have" more parallel cracks in one direction if the element size ex...
Definition fcm.h:343
oofem::FCMMaterial::computeNumerD2ModulusForCrack
virtual double computeNumerD2ModulusForCrack(GaussPoint *gp, TimeStep *tStep, int icrack) const =0
shear modulus for numerical purpose (stiffness matrix) for a given crack plane (1,...
oofem::FCMMaterial::nAllowedCracks
int nAllowedCracks
allowed number of cracks (user-defined)
Definition fcm.h:275
oofem::FloatArray::resize
void resize(Index s)
Definition floatarray.C:94
oofem::FloatArray::at
double & at(Index i)
Definition floatarray.h:202
oofem::FloatArray::giveSize
Index giveSize() const
Returns the size of receiver.
Definition floatarray.h:261
oofem::FloatArray::resizeWithValues
void resizeWithValues(Index s, std::size_t allocChunk=0)
Definition floatarray.C:103
oofem::FloatArray::beDifferenceOf
void beDifferenceOf(const FloatArray &a, const FloatArray &b)
Definition floatarray.C:403
oofem::FloatArray::zero
void zero()
Zeroes all coefficients of receiver.
Definition floatarray.C:683
oofem::FloatArray::beProductOf
void beProductOf(const FloatMatrix &aMatrix, const FloatArray &anArray)
Definition floatarray.C:689
oofem::FloatArray::beTProductOf
void beTProductOf(const FloatMatrix &aMatrix, const FloatArray &anArray)
Definition floatarray.C:721
oofem::FloatArray::add
void add(const FloatArray &src)
Definition floatarray.C:218
oofem::FloatArray::subtract
void subtract(const FloatArray &src)
Definition floatarray.C:320
oofem::FloatArray::times
void times(double s)
Definition floatarray.C:834
oofem::FloatMatrix::times
void times(double f)
Definition floatmatrix.C:1155
oofem::FloatMatrix::rotatedWith
void rotatedWith(const FloatMatrix &r, char mode='n')
Definition floatmatrix.C:1612
oofem::FloatMatrix::resizeWithData
void resizeWithData(Index, Index)
Definition floatmatrix.C:91
oofem::FloatMatrix::add
void add(const FloatMatrix &a)
Definition floatmatrix.C:962
oofem::FloatMatrix::resize
void resize(Index rows, Index cols)
Definition floatmatrix.C:79
oofem::FloatMatrix::beProductOf
void beProductOf(const FloatMatrix &a, const FloatMatrix &b)
Definition floatmatrix.C:365
oofem::FloatMatrix::zero
void zero()
Zeroes all coefficient of receiver.
Definition floatmatrix.C:1064
oofem::FloatMatrix::beInverseOf
bool beInverseOf(const FloatMatrix &src)
Definition floatmatrix.C:768
oofem::FloatMatrix::giveNumberOfRows
int giveNumberOfRows() const
Returns number of rows of receiver.
Definition floatmatrix.h:200
oofem::FloatMatrix::at
double at(std::size_t i, std::size_t j) const
Definition floatmatrix.h:221
oofem::FloatMatrix::solveForRhs
bool solveForRhs(const FloatArray &b, FloatArray &answer, bool transpose=false)
Definition floatmatrix.C:1303
oofem::FloatMatrix::subtract
void subtract(const FloatMatrix &a)
Definition floatmatrix.C:1020
oofem::GaussPoint::giveMaterialMode
MaterialMode giveMaterialMode()
Returns corresponding material mode of receiver.
Definition gausspoint.h:190
oofem::GaussPoint::giveElement
Element * giveElement()
Returns corresponding element to receiver.
Definition gausspoint.h:187
oofem::InputRecord::hasField
virtual bool hasField(InputFieldType id)=0
Returns true if record contains field identified by idString keyword.
oofem::IntArray::contains
bool contains(int value) const
Definition intarray.h:292
oofem::IntArray::at
int & at(std::size_t i)
Definition intarray.h:104
oofem::IntArray::giveSize
int giveSize() const
Definition intarray.h:211
oofem::IntegrationPointStatus::gp
GaussPoint * gp
Associated integration point.
Definition integrationpointstatus.h:71
oofem::Material::giveStatus
virtual MaterialStatus * giveStatus(GaussPoint *gp) const
Definition material.C:206
oofem::Material::isActivated
virtual bool isActivated(TimeStep *tStep) const
Definition material.h:182
oofem::Material::initTempStatus
virtual void initTempStatus(GaussPoint *gp) const
Definition material.C:221
oofem::StructuralMaterialStatus::giveStrainVector
const FloatArray & giveStrainVector() const
Returns the const pointer to receiver's strain vector.
Definition structuralms.h:103
oofem::StructuralMaterialStatus::StructuralMaterialStatus
StructuralMaterialStatus(GaussPoint *g)
Constructor. Creates new StructuralMaterialStatus with IntegrationPoint g.
Definition structuralms.C:42
oofem::StructuralMaterialStatus::giveStressVector
const FloatArray & giveStressVector() const
Returns the const pointer to receiver's stress vector.
Definition structuralms.h:105
oofem::StructuralMaterialStatus::letTempStressVectorBe
void letTempStressVectorBe(const FloatArray &v)
Assigns tempStressVector to given vector v.
Definition structuralms.h:133
oofem::StructuralMaterialStatus::letTempStrainVectorBe
void letTempStrainVectorBe(const FloatArray &v)
Assigns tempStrainVector to given vector v.
Definition structuralms.h:135
oofem::StructuralMaterial::givePlaneStressVectorTranformationMtrx
static FloatMatrixF< 3, 3 > givePlaneStressVectorTranformationMtrx(const FloatMatrixF< 2, 2 > &base, bool transpose=false)
Definition structuralmaterial.C:1859
oofem::StructuralMaterial::giveStressVectorTranformationMtrx
static FloatMatrixF< 6, 6 > giveStressVectorTranformationMtrx(const FloatMatrixF< 3, 3 > &base, bool transpose=false)
Definition structuralmaterial.C:1798
oofem::StructuralMaterial::StructuralMaterial
StructuralMaterial(int n, Domain *d)
Definition structuralmaterial.C:66
oofem::StructuralMaterial::giveStiffnessMatrix
virtual void giveStiffnessMatrix(FloatMatrix &answer, MatResponseMode mode, GaussPoint *gp, TimeStep *tStep) const
Definition structuralmaterial.C:623
oofem::StructuralMaterial::computePrincipalValDir
static void computePrincipalValDir(FloatArray &answer, FloatMatrix &dir, const FloatArray &s, stressStrainPrincMode mode)
Definition structuralmaterial.C:1337
oofem::StructuralMaterial::giveStrainVectorTranformationMtrx
static FloatMatrixF< 6, 6 > giveStrainVectorTranformationMtrx(const FloatMatrixF< 3, 3 > &base, bool transpose=false)
Definition structuralmaterial.C:1708
oofem::StructuralMaterial::giveStressDependentPartOfStrainVector
void giveStressDependentPartOfStrainVector(FloatArray &answer, GaussPoint *gp, const FloatArray &reducedStrainVector, TimeStep *tStep, ValueModeType mode) const
Definition structuralmaterial.C:772
oofem::StructuralMaterial::give2DStrainVectorTranformationMtrx
static FloatMatrixF< 3, 3 > give2DStrainVectorTranformationMtrx(const FloatMatrixF< 2, 2 > &base, bool transpose=false)
Definition structuralmaterial.C:1769
THROW_CIOERR
#define THROW_CIOERR(e)
Definition contextioerr.h:63
OOFEM_WARNING
#define OOFEM_WARNING(...)
Definition error.h:80
OOFEM_ERROR
#define OOFEM_ERROR(...)
Definition error.h:79
_IFT_FCM_multipleCrackShear
#define _IFT_FCM_multipleCrackShear
Definition fcm.h:48
_IFT_FCM_nAllowedCracks
#define _IFT_FCM_nAllowedCracks
Definition fcm.h:46
_IFT_FCM_shearCoeffNumer
#define _IFT_FCM_shearCoeffNumer
Definition fcm.h:50
pscm_NONE
#define pscm_NONE
Definition fcm.h:56
fcm_BIGNUMBER
#define fcm_BIGNUMBER
Definition fcm.h:63
_IFT_FCM_ecsm
#define _IFT_FCM_ecsm
Definition fcm.h:49
pscm_JUST_INIT
#define pscm_JUST_INIT
Definition fcm.h:57
_IFT_FCM_crackSpacing
#define _IFT_FCM_crackSpacing
Definition fcm.h:47
pscm_CLOSED
#define pscm_CLOSED
Definition fcm.h:60
pscm_UNLO_RELO
#define pscm_UNLO_RELO
Definition fcm.h:59
fcm_THRESHOLD_CRACK_STRAIN
#define fcm_THRESHOLD_CRACK_STRAIN
Definition fcm.h:65
fcm_TOLERANCE
#define fcm_TOLERANCE
Definition fcm.h:64
fcm_SMALL_STRAIN
#define fcm_SMALL_STRAIN
Definition fcm.h:62
pscm_SOFTENING
#define pscm_SOFTENING
Definition fcm.h:58
_IFT_FCM_normalCoeffNumer
#define _IFT_FCM_normalCoeffNumer
Definition fcm.h:51
IR_GIVE_OPTIONAL_FIELD
#define IR_GIVE_OPTIONAL_FIELD(__ir, __value, __id)
Definition inputrecord.h:75
IR_GIVE_FIELD
#define IR_GIVE_FIELD(__ir, __value, __id)
Definition inputrecord.h:67
oofem::ContextMode
long ContextMode
Definition contextmode.h:43
oofem::min
FloatArrayF< N > min(const FloatArrayF< N > &a, const FloatArrayF< N > &b)
Definition floatarrayf.h:522
oofem::max
FloatArrayF< N > max(const FloatArrayF< N > &a, const FloatArrayF< N > &b)
Definition floatarrayf.h:512
oofem::principal_stress
@ principal_stress
For computing principal stresses.
Definition stressstrainprincmode.h:46
oofem::sgn
double sgn(double i)
Returns the signum of given value (if value is < 0 returns -1, otherwise returns 1).
Definition mathfem.h:104
oofem::ECSM_SquareRootOfArea
@ ECSM_SquareRootOfArea
Definition internalstatetype.h:212
oofem::ECSM_ProjectionCentered
@ ECSM_ProjectionCentered
Definition internalstatetype.h:214
oofem::zeros
FloatArrayF< N > zeros()
For more readable code.
Definition floatarrayf.h:533
This page is part of the OOFEM-3.0 documentation. Copyright Copyright (C) 1994-2025 Borek Patzak Bořek Patzák
Project e-mail: oofem@fsv.cvut.cz
Generated at for OOFEM by doxygen 1.15.0 written by Dimitri van Heesch, © 1997-2011