OOFEM  2.4
OOFEM.org - Object Oriented Finite Element Solver
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 - 2013 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 
45 namespace oofem {
47  //
48  // constructor
49  //
50 {
52  linearElasticMaterial = NULL;
53 }
54 
55 
57 //
58 // destructor
59 //
60 {}
61 
62 int
64 //
65 // returns whether receiver supports given mode
66 //
67 {
68  return mode == _3dMat || mode == _PlaneStress || mode == _PlaneStrain;
69 }
70 
71 void
73  const FloatArray &totalStrain,
74  TimeStep *tStep)
75 //
76 // returns real stress vector in 3d stress space of receiver according to
77 // previous level of stress and current
78 // strain increment, the only way, how to correctly update gp records
79 //
80 {
81  // g = global direction
82  // l = local (crack) direction
83 
84  // stress and strain transformation matrices
85  FloatMatrix epsG2L, sigmaL2G;
86  // stiffness matrices (total, elastic, cracking)
87  FloatMatrix D, De, Dcr;
88  FloatMatrix principalDirs;
89 
90  // STRESSES
91  FloatArray stressIncrement_l, stressVector_g, trialStress_g, trialStress_l, sigmaResid, sigmaElast_l, sigmaCrack_l;
92 
93  // STRAINS
94  FloatArray reducedStrain_g, reducedStrain_l, strainIncrement_g, strainIncrement_l, crackStrain, crackStrainIncrement, elasticStrain_l;
95 
96  FCMMaterialStatus *status = static_cast< FCMMaterialStatus * >( this->giveStatus(gp) );
97  MaterialMode mMode = gp->giveMaterialMode();
98 
99  int nCr = status->giveNumberOfCracks();
100  int nMaxCr = status->giveMaxNumberOfCracks(gp);
101 
102  double maxErr;
103  int iterLimit;
104  double maxTau;
105 
106  // for the bisection method:
107  int index, indexCount;
108  bool exitWhileCond = false;
109  bool cancelledStrainFlag = false;
110 
111  bool plus, minus;
112  FloatArray crackStrainPlus, crackStrainMinus, residPlus, residMinus; //local
113 
114  index = indexCount = 0;
115 
116  this->initTempStatus(gp);
117 
118  // get elastic stiffness matrix
119  this->giveStiffnessMatrix(De, ElasticStiffness, gp, tStep);
120 
121 
122  // NOW INDEPENDENTLY OF THE HISTORY GET PRINCIPAL STRESSES / STRESSES IN THE CRACK DIRS
123 
124  // ELASTIC CASE - NO CRACKING SO FAR (in prevs steps)
125  if ( nCr == 0 ) {
126  this->giveStressDependentPartOfStrainVector(reducedStrain_g, gp, totalStrain, tStep, VM_Total);
127  trialStress_g.beProductOf(De, reducedStrain_g);
128 
129  this->computePrincipalValDir(trialStress_l, principalDirs, trialStress_g, principal_stress);
130 
131  // REMAINS ELASTIC or element cracking is prevented from above
132  if ( !this->isStrengthExceeded( principalDirs, gp, tStep, 1, trialStress_l.at(1) ) ) {
133  answer = trialStress_g;
134  status->letTempStrainVectorBe(totalStrain);
135  status->letTempStressVectorBe(trialStress_g);
136  return;
137 
138  // STARTS CRACKING - 1st crack
139  } else {
140  this->initializeCrack(gp, principalDirs, 1);
141 
142  for ( int iCrack = 2; iCrack <= min(nMaxCr, this->nAllowedCracks); iCrack++ ) {
143  if ( this->isStrengthExceeded( principalDirs, gp, tStep, iCrack, trialStress_l.at(iCrack) ) ) { // for "nonlocal model"
144  this->initializeCrack(gp, principalDirs, iCrack);
145  }
146  }
147  }
148 
149  // CRACKING OCCURED IN PREVIOUS STEPS
150  } else {
151  // check other principal directions - if strength is not exceeded, take transformation matrices from previous analysi
152  if ( nCr < nMaxCr ) {
153  // equilibrated global stress
154  stressVector_g = status->giveStressVector();
155  // get strain increment
156  this->giveStressDependentPartOfStrainVector(reducedStrain_g, gp, totalStrain, tStep, VM_Incremental);
157  strainIncrement_g.beDifferenceOf( reducedStrain_g, status->giveStrainVector() );
158 
159  // stress increment in crack direction
160  trialStress_g.beProductOf(De, strainIncrement_g);
161  trialStress_g.add(stressVector_g);
162 
163  principalDirs = status->giveCrackDirs();
164 
165  for ( int iCrack = nCr + 1; iCrack <= min(nMaxCr, this->nAllowedCracks); iCrack++ ) {
166  // condition to guarantee that third crack won't be initiated if second crack does not exist
167  if ( ( status->giveNumberOfTempCracks() + 1 ) < iCrack ) {
168  break;
169  }
170 
171  if ( !this->checkStrengthCriterion(principalDirs, trialStress_g, gp, tStep, iCrack) ) { // = strength of composite
172  // second and third crack is now initialized
173  this->initializeCrack(gp, principalDirs, iCrack);
174  }
175  } // end for
176  } // end nMaxCr
177  } // cracking in previous steps
178 
179  sigmaL2G = status->giveL2GStressVectorTransformationMtrx();
180  epsG2L = status->giveG2LStrainVectorTransformationMtrx();
181 
182 
183  // PREPARE STRAIN AND STRESS INCREMENT IN CRACK DIRECTIONS
184 
185  if ( ( nCr == 0 ) || ( nCr == nMaxCr ) ) { // not to do the same job twice
186  // get strain increment
187  this->giveStressDependentPartOfStrainVector(reducedStrain_g, gp, totalStrain, tStep, VM_Incremental);
188  strainIncrement_g.beDifferenceOf( reducedStrain_g, status->giveStrainVector() );
189  }
190 
191 
192  // strain and strain increment in crack direction
193  strainIncrement_l = strainIncrement_g;
194  strainIncrement_l.rotatedWith(epsG2L, 'n'); // from global to local
195 
196  reducedStrain_l = reducedStrain_g;
197  reducedStrain_l.rotatedWith(epsG2L, 'n'); // from global to local
198 
199  // stress increment in crack direction
200  stressIncrement_l.beProductOf(De, strainIncrement_l);
201 
202 
203  // SIMILARLY TO RCM2: THE STRESS INCREMENT IN ELASTIC
204  // AND CRACKING UNIT MUST BE EQUAL - ITERATE
205 
206  crackStrain = status->giveCrackStrainVector();
207 
208  iterLimit = 20;
209 
210  for ( int iter = 1; iter <= iterLimit; iter++ ) {
211  this->giveLocalCrackedStiffnessMatrix(Dcr, TangentStiffness, gp, tStep);
212 
213  if ( iter == 1 ) {
214  // residuum
215  sigmaResid = stressIncrement_l;
216  }
217 
218  D = De;
219  D.add(Dcr);
220 
221  D.solveForRhs(sigmaResid, crackStrainIncrement);
222 
223  // update and store cracking strain
224  crackStrain.add(crackStrainIncrement);
225 
226  status->setTempCrackStrainVector(crackStrain);
227 
228  // update statuses
229  this->updateCrackStatus(gp);
230 
231  // compute and compare stress in elastic and cracking unit:
232  // EL
233  elasticStrain_l = reducedStrain_l;
234  elasticStrain_l.subtract(crackStrain);
235  sigmaElast_l.beProductOf(De, elasticStrain_l);
236  // CR
237  this->giveLocalCrackedStiffnessMatrix(Dcr, SecantStiffness, gp, tStep);
238  // sigmaCrack_l.beProductOf(Dr, crackStrain); // matrix is diagonal
239  sigmaCrack_l.resize( crackStrain.giveSize() );
240  for ( int i = 1; i <= crackStrain.giveSize(); i++ ) {
241  sigmaCrack_l.at(i) = Dcr.at(i, i) * crackStrain.at(i);
242 
243  switch ( mMode ) {
244  case _PlaneStress:
245 
246  if ( i == 3 ) {
247  maxTau = this->maxShearStress(gp, 6);
248 
249  if ( sigmaCrack_l.at(i) > maxTau ) {
250  sigmaCrack_l.at(i) = maxTau;
251  }
252  }
253  break;
254 
255  case _3dMat:
256 
257  if ( i >= 4 ) {
258  maxTau = this->maxShearStress(gp, i);
259 
260  if ( sigmaCrack_l.at(i) > maxTau ) {
261  sigmaCrack_l.at(i) = maxTau;
262  }
263  }
264  break;
265 
266  case _PlaneStrain: // check
267 
268  if ( i == 4 ) {
269  maxTau = this->maxShearStress(gp, 6);
270 
271  if ( sigmaCrack_l.at(i) > maxTau ) {
272  sigmaCrack_l.at(i) = maxTau;
273  }
274  }
275  break;
276 
277  default:
278  OOFEM_ERROR( "Material mode %s not supported", __MaterialModeToString(mMode) );
279  }
280  }
281 
282  // residuum
283  sigmaResid = sigmaElast_l;
284  sigmaResid.subtract(sigmaCrack_l);
285 
286 
287  maxErr = 0.;
288  for ( int i = 1; i <= sigmaResid.giveSize(); i++ ) {
289  if ( fabs( sigmaResid.at(i) ) > maxErr ) {
290  maxErr = fabs( sigmaResid.at(i) );
291  }
292  }
293 
294 
295  if ( maxErr < fcm_TOLERANCE * this->giveTensileStrength(gp) ) {
296  break;
297  }
298  }
299 
300 
301  // modified iteration method - first gradient and then bisection method
302 
303  while ( maxErr > fcm_TOLERANCE * this->giveTensileStrength(gp) ) {
304  for ( int i = 1; i <= sigmaResid.giveSize(); i++ ) {
305  if ( fabs( sigmaResid.at(i) ) == maxErr ) {
306  // found the same index as before -> do not do the same job over, exit while
307  if ( index == i ) {
308  exitWhileCond = true;
309  break;
310  } else {
311  // the crack strain index is not the same as before - clear the following flag
312  cancelledStrainFlag = false;
313  }
314 
315  index = i;
316 
317  break;
318  }
319  }
320 
321  if ( exitWhileCond ) {
322  break;
323  }
324 
325  iterLimit = 100;
326 
327  plus = false;
328  minus = false;
329 
330  crackStrainPlus.zero();
331  crackStrainMinus.zero();
332  residPlus.zero();
333  residMinus.zero();
334 
335  indexCount++;
336 
337  if ( indexCount > 10 ) {
338  OOFEM_WARNING("Fixed crack model: Local equilibrium not reached!, max. stress error %f", maxErr);
339  break;
340  }
341 
342 
343  for ( int iter = 1; iter <= iterLimit; iter++ ) {
344  if ( iter == iterLimit ) {
345  OOFEM_WARNING("Fixed crack model: Local equilibrium not reached!, max. stress error %f", maxErr);
346  }
347 
348 
349  if ( ( iter == 1 ) && ( indexCount == 1 ) ) {
350  crackStrain.zero();
351  crackStrain.add( status->giveCrackStrainVector() );
352  status->setTempCrackStrainVector(crackStrain);
353  //sigmaResid.zero();
354  sigmaResid = stressIncrement_l;
355  this->updateCrackStatus(gp);
356  }
357 
358  // DO THE FOLLOWING SECTION ONLY UNTIL "PLUS" AND "MINUS" CRACK STRAIN IS FOUND
359 
360  if ( ( !plus ) || ( !minus ) ) {
361  this->giveLocalCrackedStiffnessMatrix(Dcr, TangentStiffness, gp, tStep);
362 
363  D = De;
364  D.add(Dcr);
365 
366  D.solveForRhs(sigmaResid, crackStrainIncrement);
367 
368  // update and store cracking strain
369  crackStrain.add(crackStrainIncrement);
370 
371  status->setTempCrackStrainVector(crackStrain);
372 
373  // update statuses
374  this->updateCrackStatus(gp);
375 
376  // compute and compare stress in elastic and cracking unit:
377  // EL
378  elasticStrain_l = reducedStrain_l;
379  elasticStrain_l.subtract(crackStrain);
380  sigmaElast_l.beProductOf(De, elasticStrain_l);
381  // CR
382  this->giveLocalCrackedStiffnessMatrix(Dcr, SecantStiffness, gp, tStep);
383  // sigmaCrack_l.beProductOf(Dr, crackStrain); // matrix is diagonal
384  sigmaCrack_l.resize( crackStrain.giveSize() );
385  for ( int i = 1; i <= crackStrain.giveSize(); i++ ) {
386  sigmaCrack_l.at(i) = Dcr.at(i, i) * crackStrain.at(i);
387 
388  switch ( mMode ) {
389  case _PlaneStress:
390 
391  if ( i == 3 ) {
392  maxTau = this->maxShearStress(gp, 6);
393 
394  if ( sigmaCrack_l.at(i) > maxTau ) {
395  sigmaCrack_l.at(i) = maxTau;
396  }
397  }
398  break;
399 
400  case _3dMat:
401 
402  if ( i >= 4 ) {
403  maxTau = this->maxShearStress(gp, i);
404 
405  if ( sigmaCrack_l.at(i) > maxTau ) {
406  sigmaCrack_l.at(i) = maxTau;
407  }
408  }
409  break;
410 
411  case _PlaneStrain: // check
412 
413  if ( i == 4 ) {
414  maxTau = this->maxShearStress(gp, 6);
415 
416  if ( sigmaCrack_l.at(i) > maxTau ) {
417  sigmaCrack_l.at(i) = maxTau;
418  }
419  }
420  break;
421 
422  default:
423  OOFEM_ERROR( "Material mode %s not supported", __MaterialModeToString(mMode) );
424  }
425  }
426 
427  // residuum
428  sigmaResid = sigmaElast_l;
429  sigmaResid.subtract(sigmaCrack_l);
430 
431  if ( sigmaResid.at(index) > 0 ) {
432  if ( !plus ) { // first positive sigma resid
433  residPlus = sigmaResid;
434  } else {
435  if ( sigmaResid.at(index) < residPlus.at(index) ) { // smaller error
436  residPlus = sigmaResid;
437  }
438  }
439 
440  crackStrainPlus = crackStrain;
441  plus = true;
442  } else { // sigmaResid < 0
443  if ( !minus ) { // first negative sigma resid
444  residMinus = sigmaResid;
445  } else {
446  if ( sigmaResid.at(index) > residMinus.at(index) ) { // smaller error
447  residMinus = sigmaResid;
448  }
449  }
450 
451  crackStrainMinus = crackStrain;
452  minus = true;
453  }
454 
455 
456  maxErr = 0.;
457  for ( int i = 1; i <= sigmaResid.giveSize(); i++ ) {
458  if ( fabs( sigmaResid.at(i) ) > maxErr ) {
459  maxErr = fabs( sigmaResid.at(i) );
460  }
461  }
462 
463 
464  if ( fabs( sigmaResid.at(index) ) < fcm_TOLERANCE * this->giveTensileStrength(gp) ) {
465  break;
466  }
467  } else { // we have "plus" and "minus" cracking and the bisection method can start
468  crackStrain.zero();
469  crackStrain.add(crackStrainMinus);
470  crackStrain.add(crackStrainPlus);
471  crackStrain.times(0.5);
472 
473  status->setTempCrackStrainVector(crackStrain);
474  // update statuses - to have correct stiffnesses (unlo-relo vs. softening)
475  this->updateCrackStatus(gp);
476 
477 
478 
479  // compute and compare stress in elastic and cracking unit:
480  // EL
481  elasticStrain_l = reducedStrain_l;
482  elasticStrain_l.subtract(crackStrain);
483  sigmaElast_l.beProductOf(De, elasticStrain_l);
484  // CR
485  this->giveLocalCrackedStiffnessMatrix(Dcr, SecantStiffness, gp, tStep);
486  // sigmaCrack_l.beProductOf(Dr, crackStrain); // matrix is diagonal
487  sigmaCrack_l.resize( crackStrain.giveSize() );
488  for ( int i = 1; i <= crackStrain.giveSize(); i++ ) {
489  sigmaCrack_l.at(i) = Dcr.at(i, i) * crackStrain.at(i);
490 
491  switch ( mMode ) {
492  case _PlaneStress:
493 
494  if ( i == 3 ) {
495  maxTau = this->maxShearStress(gp, 6);
496 
497  if ( sigmaCrack_l.at(i) > maxTau ) {
498  sigmaCrack_l.at(i) = maxTau;
499  }
500  }
501  break;
502 
503  case _3dMat:
504 
505  if ( i >= 4 ) {
506  maxTau = this->maxShearStress(gp, i);
507 
508  if ( sigmaCrack_l.at(i) > maxTau ) {
509  sigmaCrack_l.at(i) = maxTau;
510  }
511  }
512  break;
513 
514  case _PlaneStrain: // check
515 
516  if ( i == 4 ) {
517  maxTau = this->maxShearStress(gp, 6);
518 
519  if ( sigmaCrack_l.at(i) > maxTau ) {
520  sigmaCrack_l.at(i) = maxTau;
521  }
522  }
523  break;
524 
525  default:
526  OOFEM_ERROR( "Material mode %s not supported", __MaterialModeToString(mMode) );
527  }
528  }
529 
530  // residuum
531  sigmaResid = sigmaElast_l;
532  sigmaResid.subtract(sigmaCrack_l);
533 
534 
535  maxErr = 0.;
536  for ( int i = 1; i <= sigmaResid.giveSize(); i++ ) {
537  if ( fabs( sigmaResid.at(i) ) > maxErr ) {
538  maxErr = fabs( sigmaResid.at(i) );
539  }
540  }
541 
542  // condition for the given stress component
543  if ( fabs( sigmaResid.at(index) ) < fcm_TOLERANCE * this->giveTensileStrength(gp) ) {
544  break;
545  }
546 
547 
548  if ( sigmaResid.at(index) > 0 ) {
549  if ( sigmaResid.at(index) < residPlus.at(index) ) { // smaller error
550  residPlus = sigmaResid;
551  }
552 
553  crackStrainPlus = crackStrain;
554  } else { // sigmaResid < 0
555  if ( sigmaResid.at(index) > residMinus.at(index) ) { // smaller error
556  residMinus = sigmaResid;
557  }
558 
559  crackStrainMinus = crackStrain;
560  }
561 
562 
563  if ( fabs( crackStrainMinus.at(index) - crackStrainPlus.at(index) ) < 1.e-18 ) {
564  crackStrain.zero();
565  crackStrain.add(crackStrainMinus);
566  crackStrain.add(crackStrainPlus);
567  crackStrain.times(0.5);
568 
569  // if ( fabs (crackStrain.at(index) < 1.e-18 ) ) {
570  crackStrain.at(index) = 0.;
571  cancelledStrainFlag = true;
572  // }
573 
574 #if DEBUG
575  if ( cancelledStrainFlag ) {
576  OOFEM_WARNING("Fixed crack model: cracking strain component %d set to zero", index);
577  }
578 #endif
579 
580  break;
581  }
582  } // plus-minus condition
583  } // iter looop
584  } // maxErr exceeded
585 
586 
587  if ( ( maxErr >= fcm_TOLERANCE * this->giveTensileStrength(gp) ) && ( !cancelledStrainFlag ) ) {
588  OOFEM_WARNING("Fixed crack model: Local equilibrium not reached!, max. stress error %f", maxErr);
589  }
590 
591  // set the final-equilibrated crack strain and store the correct corresponding statuses
592  status->setTempCrackStrainVector(crackStrain);
593  // update statuses - to have correct stiffnesses
594  this->updateCrackStatus(gp);
595 
596 
597  // correct non-zeros, statuses, etc...
598  for ( int i = 1; i <= nMaxCr; i++ ) {
599 #if DEBUG
600  // check for NaNs
601  if ( crackStrain.at(i) != crackStrain.at(i) ) {
602  OOFEM_ERROR( "crack strain is NaN: %f", crackStrain.at(i) );
603  }
604 #endif
605 
606  if ( ( status->giveTempCrackStatus(i) == pscm_NONE ) && ( crackStrain.at(i) <= 0. ) ) {
607  status->setTempCrackStrain(i, 0.);
608  } else if ( ( status->giveTempCrackStatus(i) == pscm_JUST_INIT ) && ( crackStrain.at(i) <= fcm_SMALL_STRAIN ) ) {
609  if ( i + 1 > status->giveMaxNumberOfCracks(gp) ) {
610  status->setTempCrackStatus(i, pscm_NONE);
611  status->setTempCrackStrain(i, 0.);
612  if ( i == 1 ) {
613  principalDirs.zero();
614  status->setCrackDirs(principalDirs);
615  }
616 
617  // can't change crack if the subsequent crack exists
618  } else if ( status->giveTempCrackStatus(i + 1) != pscm_NONE ) {} else {
619  status->setTempCrackStatus(i, pscm_NONE);
620  status->setTempCrackStrain(i, 0.);
621  if ( i == 1 ) {
622  principalDirs.zero();
623  status->setCrackDirs(principalDirs);
624  }
625  }
626  } else if ( crackStrain.at(i) < 0. ) {
627  status->setTempCrackStrain(i, 0.);
628  }
629  }
630 
631  // at first crack initiation there should not be any other cracking strains
632  if ( ( status->giveTempCrackStatus(1) == pscm_JUST_INIT ) && ( status->giveTempCrackStatus(2) == pscm_NONE ) ) {
633  status->setTempCrackStrain(2, 0.);
634  status->setTempCrackStrain(3, 0.);
635  }
636 
637 
638  // convert from local to global coordinates
639  sigmaElast_l.rotatedWith(sigmaL2G, 'n');
640 
641  answer = sigmaElast_l;
642  status->letTempStrainVectorBe(totalStrain);
643  status->letTempStressVectorBe(sigmaElast_l);
644 
645  return;
646 }
647 
648 
649 void
651 {
652  MaterialMode mMode = gp->giveMaterialMode();
653  FCMMaterialStatus *status = static_cast< FCMMaterialStatus * >( this->giveStatus(gp) );
654 
655  FloatArray crackVector(3);
656  crackVector.zero();
657 
658  FloatMatrix epsG2L, sigmaG2L, epsL2G, sigmaL2G;
659 
660  int nMaxCracks;
661  double Le;
662 
663  for ( int i = 1; i <= base.giveNumberOfRows(); i++ ) {
664  crackVector.at(i) = base.at(i, nCrack);
665  }
666 
667  Le = this->giveCharacteristicElementLength(gp, crackVector);
668  status->setCharLength(nCrack, Le);
669 
670  this->checkSnapBack(gp, nCrack);
671 
672  status->setTempCrackStatus(nCrack, pscm_JUST_INIT);
673  //status->setTempCrackStatus(nCrack, pscm_SOFTENING);
674 
675  // change the base crack vector and transformation matrices
676  // when the initiated crack is not the last possible one
677 
678  // return for the second crack in plane stress - base is set and transformation matrices are given by the first crack
679  if ( ( mMode == _PlaneStress ) && ( nCrack == 2 ) ) {
680  return;
681  }
682 
683  nMaxCracks = status->giveMaxNumberOfCracks(gp);
684 
685  if ( nCrack <= nMaxCracks ) {
686  // make sure that the matrix of base vectors is right-handed
687  if ( nMaxCracks == 3 ) {
688  base.at(1, 3) = base.at(2, 1) * base.at(3, 2) - base.at(3, 1) * base.at(2, 2);
689  base.at(2, 3) = base.at(3, 1) * base.at(1, 2) - base.at(1, 1) * base.at(3, 2);
690  base.at(3, 3) = base.at(1, 1) * base.at(2, 2) - base.at(2, 1) * base.at(1, 2);
691  }
692 
693  // status->setTempCrackDirs(base);
694  status->setCrackDirs(base);
695 
696  if ( mMode == _PlaneStress ) {
697  this->givePlaneStressVectorTranformationMtrx(sigmaG2L, base, 0);
698  this->give2DStrainVectorTranformationMtrx(epsG2L, base, 0);
699 
700  this->givePlaneStressVectorTranformationMtrx(sigmaL2G, base, 1);
701  this->give2DStrainVectorTranformationMtrx(epsL2G, base, 1);
702  } else {
703  this->giveStressVectorTranformationMtrx(sigmaG2L, base, 0);
704  this->giveStrainVectorTranformationMtrx(epsG2L, base, 0);
705 
706  this->giveStressVectorTranformationMtrx(sigmaL2G, base, 1);
707  this->giveStrainVectorTranformationMtrx(epsL2G, base, 1);
708  }
709 
710  status->setG2LStressVectorTransformationMtrx(sigmaG2L);
711  status->setG2LStrainVectorTransformationMtrx(epsG2L);
712 
713  status->setL2GStressVectorTransformationMtrx(sigmaL2G);
714  status->setL2GStrainVectorTransformationMtrx(epsL2G);
715  }
716 }
717 
718 
719 bool
721  FCMMaterialStatus *status = static_cast< FCMMaterialStatus * >( this->giveStatus(gp) );
722 
723  if ( icrack >= 4 ) {
724  OOFEM_ERROR("Unexpected crack number");
725  }
726 
727  if ( ( status->giveTempCrackStatus(icrack) != pscm_NONE ) && ( status->giveTempCrackStatus(icrack) != pscm_CLOSED ) ) {
728  return false;
729  } else {
730  return true;
731  }
732 }
733 
734 
735 
736 
737 bool
739  FCMMaterialStatus *status = static_cast< FCMMaterialStatus * >( this->giveStatus(gp) );
740 
741  int normal_1, normal_2;
742 
743  if ( i == 4 ) { // y-z
744  normal_1 = 2;
745  normal_2 = 3;
746  } else if ( i == 5 ) { // x-z
747  normal_1 = 1;
748  normal_2 = 3;
749  } else if ( i == 6 ) { // x-y
750  normal_1 = 1;
751  normal_2 = 2;
752  } else {
753  OOFEM_ERROR("Unexpected number for shear stress (must be either 4, 5 or 6).");
754  normal_1 = normal_2 = 0;
755  }
756 
757  if ( ( status->giveTempCrackStatus(normal_1) != pscm_NONE ) && ( status->giveTempCrackStatus(normal_1) != pscm_CLOSED ) ) {
758  return false;
759  } else if ( ( status->giveTempCrackStatus(normal_2) != pscm_NONE ) && ( status->giveTempCrackStatus(normal_2) != pscm_CLOSED ) ) {
760  return false;
761  } else {
762  return true;
763  }
764 }
765 
766 
767 
768 double
770  FCMMaterialStatus *status = static_cast< FCMMaterialStatus * >( this->giveStatus(gp) );
771 
772  double crackOpening, N;
773 
774  if ( i > status->giveNumberOfTempCracks() ) {
775  crackOpening = 0;
776  } else {
777  crackOpening = max(status->giveCharLength(i) * status->giveTempCrackStrain(i), 0.);
778  N = this->giveNumberOfCracksInDirection(gp, i);
779  crackOpening /= N;
780  }
781 
782  return crackOpening;
783 }
784 
785 
786 
787 double
789  FCMMaterialStatus *status = static_cast< FCMMaterialStatus * >( this->giveStatus(gp) );
790 
791  double crackOpening, N;
792 
793  if ( i > status->giveNumberOfTempCracks() ) {
794  crackOpening = 0;
795  } else {
796  crackOpening = max(status->giveCharLength(i) * status->giveMaxCrackStrain(i), 0.);
797  N = this->giveNumberOfCracksInDirection(gp, i);
798  crackOpening /= N;
799  }
800 
801  return crackOpening;
802 }
803 
804 
805 
806 double
808  MaterialMode mMode = gp->giveMaterialMode();
809 
810  FCMMaterialStatus *status = static_cast< FCMMaterialStatus * >( this->giveStatus(gp) );
811 
812  // total shear slip
813  double slip = 0;
814 
815  // number of cracks in the direction of iCrack
816  int nCracks = this->giveNumberOfCracksInDirection(gp, icrack);
817 
818  // shear directions on icrack plane
819  int dir_j, dir_k;
820 
821  // cracking shear strains
822  double gamma_cr_ij, gamma_cr_ik;
823 
824  // factor for redistribution of shear crack strain
825  double factor_ij, factor_ik;
826 
827  double u_ij, u_ik;
828 
829  if ( status->giveTempCrackStatus(icrack) == pscm_NONE ) { // crack not initiated
830  return slip;
831  }
832 
833  if ( icrack > 3 ) {
834  OOFEM_ERROR("Unexpected value of index i (4, 5, 6 permitted only)");
835  }
836 
837 
838  if ( ( mMode == _PlaneStress ) || ( mMode == _PlaneStrain ) ) {
839  // get shear strain
840  if ( mMode == _PlaneStress ) {
841  gamma_cr_ij = status->giveTempCrackStrain(3);
842  } else {
843  gamma_cr_ij = status->giveTempCrackStrain(6);
844  }
845 
846  factor_ij = 1.;
847 
848  // determine if it necessary to split the cracking strain
849  // both cracks exist in 2D
850 
851  if ( ( icrack == 2 ) || ( status->giveTempCrackStatus(2) != pscm_NONE ) ) {
852  if ( icrack == 1 ) {
853  dir_j = 2;
854  } else { // icrack == 2
855  dir_j = 1;
856  }
857 
858 
859  // 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.
860 
861  factor_ij = this->computeShearStiffnessRedistributionFactor(gp, icrack, dir_j);
862  }
863 
864  slip = factor_ij * fabs(gamma_cr_ij) * status->giveCharLength(icrack) / nCracks;
865  } else if ( mMode == _3dMat ) {
866  if ( icrack == 1 ) {
867  dir_j = 2;
868  dir_k = 3;
869 
870  gamma_cr_ij = status->giveTempCrackStrain(6);
871  gamma_cr_ik = status->giveTempCrackStrain(5);
872  } else if ( icrack == 2 ) {
873  dir_j = 1;
874  dir_k = 3;
875 
876  gamma_cr_ij = status->giveTempCrackStrain(6);
877  gamma_cr_ik = status->giveTempCrackStrain(4);
878  } else { // icrack == 3
879  dir_j = 1;
880  dir_k = 2;
881 
882  gamma_cr_ij = status->giveTempCrackStrain(5);
883  gamma_cr_ik = status->giveTempCrackStrain(4);
884  }
885 
886 
887  factor_ij = 1.;
888  factor_ik = 1.;
889 
890 
891  if ( ( status->giveTempCrackStatus(dir_j) != pscm_NONE ) || ( status->giveTempCrackStatus(dir_k) != pscm_NONE ) ) {
892  if ( status->giveTempCrackStatus(dir_j) != pscm_NONE ) {
893  factor_ij = this->computeShearStiffnessRedistributionFactor(gp, icrack, dir_j);
894  }
895 
896  if ( status->giveTempCrackStatus(dir_k) != pscm_NONE ) {
897  factor_ik = this->computeShearStiffnessRedistributionFactor(gp, icrack, dir_k);
898  }
899  }
900 
901  u_ij = factor_ij * fabs(gamma_cr_ij) * status->giveCharLength(icrack) / nCracks;
902  u_ik = factor_ik * fabs(gamma_cr_ik) * status->giveCharLength(icrack) / nCracks;
903 
904  slip = sqrt( pow(u_ij, 2) + pow(u_ik, 2) );
905  } else {
906  OOFEM_ERROR( "Material mode %s not supported", __MaterialModeToString(mMode) );
907  }
908 
909  return slip;
910 }
911 
912 
913 bool
914 FCMMaterial :: isStrengthExceeded(const FloatMatrix &base, GaussPoint *gp, TimeStep *tStep, int iCrack, double trialStress) {
915  if ( trialStress > this->giveTensileStrength(gp) ) {
916  return true;
917  } else {
918  return false;
919  }
920 }
921 
922 
923 
924 
925 double
926 FCMMaterial :: computeShearStiffnessRedistributionFactor(GaussPoint *gp, int ithCrackPlane, int jthCrackDirection) {
927  double factor_ij;
928  double D2_i, D2_j;
929 
930  D2_i = this->computeD2ModulusForCrack(gp, ithCrackPlane);
931  D2_j = this->computeD2ModulusForCrack(gp, jthCrackDirection);
932 
933  factor_ij = D2_j / ( D2_i + D2_j );
934 
935  return factor_ij;
936 }
937 
938 
939 bool
940 FCMMaterial :: checkStrengthCriterion(FloatMatrix &newBase, const FloatArray &globalStress, GaussPoint *gp, TimeStep *tStep, int nCrack) {
941  FCMMaterialStatus *status = static_cast< FCMMaterialStatus * >( this->giveStatus(gp) );
942 
943  double sigX, sigY, tau, sig2;
944  FloatArray trialStress, planeStress, princStress, crackingStrain, shearStrains, newShearStrains;
945  FloatMatrix sigmaG2L, princCrackDir, oldBase, princCrackDirExt;
946 
947  sigmaG2L = status->giveG2LStressVectorTransformationMtrx();
948 
949  trialStress = globalStress;
950 
951  // rotate stress to local coordinates
952  trialStress.rotatedWith(sigmaG2L, 'n');
953 
954  if ( status->giveMaxNumberOfCracks(gp) < 3 ) { // plane stress
955  // test material strength but keep crack coordinates
956  if ( this->isStrengthExceeded( newBase, gp, tStep, nCrack, trialStress.at(nCrack) ) ) {
957  return false;
958  } else {
959  return true;
960  }
961  } else if ( nCrack == 3 ) { // 3D but is the last crack
962  if ( this->isStrengthExceeded( newBase, gp, tStep, nCrack, trialStress.at(nCrack) ) ) {
963  return false;
964  } else {
965  return true;
966  }
967  } else { // second crack in 3D case
968  // if the stress-strength criterion is violated
969  // the crack coordinates have to be changed
970  // in-plane shear strain remains unchanged
971 
972  // compute second principal stress cheap - if passed, compute eigenvectors and eigenvalues
973  sigX = trialStress.at(2); //1st normal stress in crack plane
974  sigY = trialStress.at(3); //2nd normal stress in crack plane
975  tau = trialStress.at(4); // shear stress in crack plane
976 
977  sig2 = ( sigX + sigY ) / 2. + sqrt( ( sigX - sigY ) * ( sigX - sigY ) / 4. + tau * tau );
978 
979  if ( this->isStrengthExceeded(newBase, gp, tStep, 2, sig2) ) {
980  oldBase = newBase;
981  planeStress.resize(3);
982  planeStress.zero();
983 
984  planeStress.at(1) = trialStress.at(2); //1st normal stress in crack plane
985  planeStress.at(2) = trialStress.at(3); //2nd normal stress in crack plane
986  planeStress.at(3) = trialStress.at(4); //shear in crack plane
987 
988 
989  // compute principal stresses on the already existing crack plane
990  this->computePrincipalValDir(princStress, princCrackDir, planeStress, principal_stress);
991 
992  // right-hand orientation
993  princCrackDir.at(1, 2) = -1. * princCrackDir.at(2, 1);
994  princCrackDir.at(2, 2) = princCrackDir.at(1, 1);
995 
996  // establish vector of crack directions on the crack plane
997  // the first vector is the normal direction, [1, 0, 0]
998 
999  princCrackDirExt.resize(3, 3);
1000  princCrackDirExt.zero();
1001  princCrackDirExt.at(1, 1) = 1.;
1002 
1003  for ( int i = 1; i <= 2; i++ ) {
1004  for ( int j = 1; j <= 2; j++ ) {
1005  princCrackDirExt.at(i + 1, j + 1) = princCrackDir.at(i, j);
1006  }
1007  }
1008 
1009  // new crack base vector in global coordinates
1010  newBase.beProductOf(oldBase, princCrackDirExt);
1011 
1012  /*
1013  * // algorithm which gives the same results - rotation of axes in space
1014  * double theta;
1015  *
1016  * if ( princCrackDir.at(1,1) > 0. ) {
1017  * theta = asin( princCrackDir.at(2,1) );
1018  * } else {
1019  * theta = M_PI - asin ( princCrackDir.at(2,1) );
1020  * }
1021  *
1022  * if (theta < 0.) {
1023  * theta += 2 * M_PI;
1024  * }
1025  *
1026  *
1027  * FloatArray baseVector2;
1028  * FloatArray baseVector3;
1029  *
1030  * FloatArray newBaseVector2;
1031  * FloatArray newBaseVector3;
1032  *
1033  * baseVector2.resize(3);
1034  * baseVector3.resize(3);
1035  *
1036  * newBaseVector2.resize(3);
1037  * newBaseVector3.resize(3);
1038  *
1039  * for (int i = 1; i <= 3; i++) {
1040  * baseVector2.at(i) = oldBase.at(i,2);
1041  * baseVector3.at(i) = oldBase.at(i,3);
1042  * }
1043  *
1044  * FloatMatrix K;
1045  * FloatMatrix R;
1046  * FloatMatrix help;
1047  *
1048  * K.resize(3,3);
1049  * K.zero();
1050  *
1051  * K.at(1,2) = -1. * oldBase.at(3,1);
1052  * K.at(1,3) = oldBase.at(2,1);
1053  * K.at(2,1) = oldBase.at(3,1);
1054  * K.at(2,3) = -1. * oldBase.at(1,1);
1055  * K.at(3,1) = -1. * oldBase.at(2,1);
1056  * K.at(3,2) = oldBase.at(1,1);
1057  *
1058  * R.resize(3,3);
1059  * R.beUnitMatrix();
1060  *
1061  * help = K;
1062  * help.times( sin(theta) );
1063  *
1064  * R.add(help);
1065  *
1066  * help.beProductOf (K, K);
1067  * help.times( 1.-cos(theta) );
1068  *
1069  * R.add(help);
1070  *
1071  * newBaseVector2.beProductOf(R, baseVector2);
1072  * newBaseVector3.beProductOf(R, baseVector3);
1073  */
1074 
1075 
1076  /*
1077  * // the same way - R evaluated symbolically in Matlab
1078  * // c = cos, s = sin
1079  * // k = axis of rotation
1080  * R.at(1,1) = 1. + (c-1.) * k2*k2 + (c-1.) * k3*k3;
1081  * R.at(1,2) = -k3*s - k1*k2 * (c-1.);
1082  * R.at(1,3) = k2*s - k1*k3 * (c-1.);
1083  *
1084  * R.at(2,1) = k3*s - k1*k2 * (c-1.);
1085  * R.at(2,2) = 1. + (c-1.) * k1*k1 + (c-1.) * k3*k3;
1086  * R.at(2,3) = -k1*s - k2*k3 * (c-1.);
1087  *
1088  * R.at(3,1) = -k2*s - k1*k3 * (c-1.);
1089  * R.at(3,2) = k1*s - k2*k3 * (c-1.);
1090  * R.at(3,3) = 1. + (c-1.) * k1*k1 + (c-1.) * k2*k2;
1091  */
1092 
1093 
1094  // modify shear strains
1095 
1096  crackingStrain = status->giveCrackStrainVector();
1097 
1098  shearStrains.resize(2);
1099  shearStrains.at(1) = crackingStrain.at(6); //aka x-y
1100  shearStrains.at(2) = crackingStrain.at(5); //aka x-z
1101 
1102  // gamma_xy_l [ c s ] gamma_xy_g
1103  // = [ ] *
1104  // gamma_xz_l [ -s c ] gamma_xz_g
1105 
1106  newShearStrains.beTProductOf(princCrackDir, shearStrains);
1107 
1108  crackingStrain.at(5) = newShearStrains.at(2); // aka x-z
1109  crackingStrain.at(6) = newShearStrains.at(1); // aka x-y
1110 
1111  status->setCrackStrainVector(crackingStrain);
1112  status->setTempCrackStrainVector(crackingStrain);
1113 
1114  // modify maximum shear strains
1115  // get max shear strains in the original configuration
1116  shearStrains.at(1) = status->giveMaxCrackStrain(6); //aka x-y
1117  shearStrains.at(2) = status->giveMaxCrackStrain(5); //aka x-z
1118 
1119  // gamma_xy_l [ c s ] gamma_xy_g
1120  // = [ ] *
1121  // gamma_xz_l [ -s c ] gamma_xz_g
1122 
1123  newShearStrains.beTProductOf(princCrackDir, shearStrains);
1124 
1125  status->setMaxCrackStrain( 5, max( newShearStrains.at(2), crackingStrain.at(5) ) ); // aka x-z
1126  status->setMaxCrackStrain( 6, max( newShearStrains.at(1), crackingStrain.at(6) ) ); // aka x-y
1127 
1128  status->setTempMaxCrackStrain( 5, max( newShearStrains.at(2), crackingStrain.at(5) ) ); // aka x-z
1129  status->setTempMaxCrackStrain( 6, max( newShearStrains.at(1), crackingStrain.at(6) ) ); // aka x-y
1130 
1131  return false;
1132  } else { // strength not reached
1133  return true;
1134  }
1135  } // number of crack in given stress-state
1136 }
1137 
1138 
1139 double
1141 {
1142  return gp->giveElement()->giveCharacteristicLength(crackPlaneNormal);
1143 }
1144 
1145 
1146 void
1148 {
1149  FCMMaterialStatus *status = static_cast< FCMMaterialStatus * >( this->giveStatus(gp) );
1150 
1151  int nMaxCracks = status->giveMaxNumberOfCracks(gp);
1152 
1153  FloatArray crackStrain, maxCrackStrain;
1154  crackStrain = status->giveTempCrackStrainVector();
1155  maxCrackStrain = status->giveMaxCrackStrainVector();
1156 
1157 
1158  // loop over NORMAL components of the crack strain vector
1159  // statuses are only for normal components
1160  for ( int i = 1; i <= nMaxCracks; i++ ) {
1161  // changes in status only in previously cracked (or at least initiated in this step) material
1162  if ( status->giveTempCrackStatus(i) != pscm_NONE ) {
1163  // temp crack bigger or equal than max strain so far
1164  if ( crackStrain.at(i) >= maxCrackStrain.at(i) ) {
1165  status->setTempMaxCrackStrain( i, crackStrain.at(i) );
1166 
1167  if ( status->giveCrackStatus(i) != pscm_NONE ) { //already existing crack in prev step
1168  status->setTempCrackStatus(i, pscm_SOFTENING);
1169  } else if ( ( fabs( crackStrain.at(i) ) <= 1.e-18 ) && ( status->giveCrackStatus(i) == pscm_NONE ) ) { // just initiated crack with zero strain
1170  if ( i == nMaxCracks ) {
1171  status->setTempCrackStatus(i, pscm_NONE);
1172  } else if ( nMaxCracks > i ) { // curent crack is not the last possible one
1173  if ( status->giveTempCrackStatus(i + 1) == pscm_NONE ) { // the next crack is not initiated
1174  status->setTempCrackStatus(i, pscm_NONE); // reset the current one to NONE
1175  } else {
1176  status->setTempCrackStatus(i, pscm_JUST_INIT); // else keep as just initiated
1177  }
1178  } else {
1179  status->setTempCrackStatus(i, pscm_JUST_INIT);
1180  }
1181  } else {
1182  status->setTempCrackStatus(i, pscm_JUST_INIT);
1183  }
1184 
1185  // closed previously
1186  } else if ( crackStrain.at(i) <= 0. ) {
1187  status->setTempCrackStatus(i, pscm_CLOSED);
1188 
1189  if ( status->giveMaxCrackStrain(i) == 0. ) { //not existing crack in previous step that wants to close?
1190  status->setTempMaxCrackStrain(i, 0.);
1191  // status->setTempCrackStatus(i, pscm_NONE);
1192  // CLOSED status should have similar behavior and it is safer. The status can be changed within one iteration loop on the local scale
1193  // the crack can be changed from CLOSED to NONE during overall updating
1194  // status->setTempCrackStatus(i, pscm_CLOSED);
1195  }
1196 
1197  // unloading--reloading
1198  } else if ( crackStrain.at(i) < maxCrackStrain.at(i) ) {
1199  status->setTempCrackStatus(i, pscm_UNLO_RELO);
1200  } else {
1201  OOFEM_ERROR("Unexpected value of cracking strain");
1202  }
1203  }
1204  }
1205 
1206  // SHEAR components
1207  for ( int i = nMaxCracks + 1; i <= crackStrain.giveSize(); i++ ) {
1208  if ( fabs( crackStrain.at(i) ) > maxCrackStrain.at(i) ) {
1209  status->setTempMaxCrackStrain( i, fabs( crackStrain.at(i) ) );
1210  }
1211  }
1212 }
1213 
1214 
1215 
1216 
1217 
1218 void
1220  MatResponseMode rMode, GaussPoint *gp,
1221  TimeStep *tStep)
1222 {
1223  int dim, j;
1224  FloatMatrix Dcr;
1225  IntArray indx;
1227 
1228  dim = indx.giveSize();
1229 
1230  Dcr.resize(dim, dim);
1231 
1232  for ( int i = 1; i <= dim; i++ ) {
1233  j = indx.at(i);
1234 
1235  if ( j <= 3 ) {
1236  Dcr.at(i, i) = this->giveCrackingModulus(rMode, gp, j);
1237  } else {
1238  Dcr.at(i, i) = this->computeTotalD2Modulus(gp, j);
1239  }
1240  }
1241 
1242  answer = Dcr;
1243 }
1244 
1245 
1246 
1247 
1248 
1249 void
1251  MatResponseMode rMode,
1252  GaussPoint *gp,
1253  TimeStep *tStep)
1254 //
1255 // returns effective material stiffness matrix in full form
1256 // for gp stress strain mode
1257 //
1258 {
1259  FCMMaterialStatus *status = static_cast< FCMMaterialStatus * >( this->giveStatus(gp) );
1260  StructuralMaterial *lMat = static_cast< StructuralMaterial * >( this->giveLinearElasticMaterial() );
1261 
1262  MaterialMode mMode = gp->giveMaterialMode();
1263 
1264  int numberOfActiveCracks, nMaxCracks;
1265  double overallElasticStiffness;
1266 
1267  FloatMatrix D, De, DeHelp, Dcr, DcrHelp, DcrHelp2, inv, stiffnessL2G;
1268 
1269  numberOfActiveCracks = status->giveNumberOfTempCracks();
1270 
1271  // ELASTIC MATRIX
1272  if ( ( rMode == ElasticStiffness ) || ( numberOfActiveCracks == 0 ) ) {
1273  lMat->giveStiffnessMatrix(D, rMode, gp, tStep);
1274 
1275  overallElasticStiffness = this->computeOverallElasticStiffness();
1276  if ( overallElasticStiffness != ( this->give('E', gp) ) ) {
1277  D.times( overallElasticStiffness / ( this->give('E', gp) ) );
1278  }
1279 
1280  answer = D;
1281  return;
1282  }
1283 
1284  // SECANT OR TANGENT MATRIX - first in local direction
1285  nMaxCracks = status->giveMaxNumberOfCracks(gp);
1286  lMat->giveStiffnessMatrix(De, rMode, gp, tStep);
1287 
1288  overallElasticStiffness = this->computeOverallElasticStiffness();
1289  if ( overallElasticStiffness != ( this->give('E', gp) ) ) {
1290  De.times( overallElasticStiffness / ( this->give('E', gp) ) );
1291  }
1292 
1293  // extract 1x1 / 2x2 / 3x3 submatrix
1294  De.resizeWithData(nMaxCracks, nMaxCracks);
1295 
1296  // shrink to square matrix number of cracks x number of cracks
1297  DeHelp = De;
1298  DeHelp.resizeWithData(numberOfActiveCracks, numberOfActiveCracks);
1299 
1300  Dcr.resize(numberOfActiveCracks, numberOfActiveCracks);
1301  Dcr.zero();
1302 
1303  for ( int i = 1; i <= numberOfActiveCracks; i++ ) {
1304  Dcr.at(i, i) = this->giveCrackingModulus(rMode, gp, i);
1305  }
1306 
1307  Dcr.add(DeHelp);
1308  inv.beInverseOf(Dcr);
1309  inv.resizeWithData(nMaxCracks, nMaxCracks);
1310 
1311 
1312  DcrHelp.beProductOf(De, inv); // De (De + Dcr)^-1
1313  DcrHelp2.beProductOf(DcrHelp, De); // De (De + Dcr)^-1 De
1314 
1315  D.zero();
1316  D.resize(nMaxCracks, nMaxCracks);
1317  D.add(De);
1318  D.subtract(DcrHelp2); // De - De (De + Dcr)^-1 De
1319 
1320 
1321  // resize add shear moduli on diagonal
1322  switch ( mMode ) {
1323  case _PlaneStress:
1324  D.resizeWithData(3, 3);
1325  D.at(3, 3) = this->computeEffectiveShearModulus(gp, 6);
1326 
1327  break;
1328 
1329  case _3dMat:
1330  D.resizeWithData(6, 6);
1331  for ( int i = 4; i <= 6; i++ ) {
1332  D.at(i, i) = this->computeEffectiveShearModulus(gp, i);
1333  }
1334  break;
1335 
1336  case _PlaneStrain:
1337  D.resizeWithData(6, 6);
1338  D.at(6, 6) = this->computeEffectiveShearModulus(gp, 6);
1339  break;
1340 
1341  default:
1342  OOFEM_ERROR( "Material mode %s not supported", __MaterialModeToString(mMode) );
1343  }
1344 
1345  // transform stiffnes to global c.s
1346  // to transform stiffness from LOCAL to GLOBAL coordinate system
1347  // the transformation matrix is the same as for strain transformation
1348  // from GLOBAL to LOCAL coordinate system
1349  stiffnessL2G = status->giveG2LStrainVectorTransformationMtrx();
1350  D.rotatedWith(stiffnessL2G, 'n');
1351 
1352 
1353  switch ( mMode ) {
1354  case _PlaneStress:
1355  answer = D;
1356  return;
1357 
1358  break;
1359  case _3dMat:
1360  answer = D;
1361  return;
1362 
1363  break;
1364  case _PlaneStrain:
1366  return;
1367 
1368  break;
1369  default:
1370  OOFEM_ERROR( "Material mode %s not supported", __MaterialModeToString(mMode) );
1371  return;
1372  }
1373 }
1374 
1375 
1376 
1377 double
1379 {
1380  double D2 = 0.;
1381  int crackA, crackB;
1382  double D2_1, D2_2;
1383 
1384  if ( this->isIntactForShear(gp, shearDirection) ) {
1386  } else {
1387  if ( shearDirection == 4 ) {
1388  crackA = 2;
1389  crackB = 3;
1390  } else if ( shearDirection == 5 ) {
1391  crackA = 1;
1392  crackB = 3;
1393  } else if ( shearDirection == 6 ) {
1394  crackA = 1;
1395  crackB = 2;
1396  } else {
1397  OOFEM_ERROR("Unexpected value of index i (4, 5, 6 permitted only)");
1398  crackA = crackB = 0;
1399  }
1400 
1401  if ( ( this->isIntact(gp, crackA) ) || ( this->isIntact(gp, crackB) ) ) {
1402  if ( this->isIntact(gp, crackA) ) {
1403  D2 = this->computeD2ModulusForCrack(gp, crackB);
1404  } else {
1405  D2 = this->computeD2ModulusForCrack(gp, crackA);
1406  }
1407  } else {
1408  D2_1 = this->computeD2ModulusForCrack(gp, crackA);
1409  D2_2 = this->computeD2ModulusForCrack(gp, crackB);
1410 
1411  if ( multipleCrackShear ) {
1412  // serial coupling of stiffnesses 1/D = 1/D1 + 1/D2
1413  D2 = D2_1 * D2_2 / ( D2_1 + D2_2 );
1414  } else {
1415  D2 = min(D2_1, D2_2);
1416  }
1417  }
1418  }
1419 
1420  return D2;
1421 }
1422 
1425 {
1426  IRResultType result; // Required by IR_GIVE_FIELD macro
1427 
1429  if ( result != IRRT_OK ) {
1430  return result;
1431  }
1432 
1433  this->nAllowedCracks = 3;
1435 
1436 
1437  this->crackSpacing = -1.;
1438  if ( ir->hasField(_IFT_FCM_crackSpacing) ) {
1440  }
1441 
1442  this->multipleCrackShear = false;
1444  this->multipleCrackShear = true;
1445  }
1446 
1447  int ecsm = 0;
1449  switch ( ecsm ) {
1450  case 1: ecsMethod = ECSM_SquareRootOfArea;
1451  break;
1453  break;
1454  case 3: ecsMethod = ECSM_Oliver1;
1455  break;
1456  case 4: ecsMethod = ECSM_Oliver1modified;
1457  break;
1458  default: ecsMethod = ECSM_Projection;
1459  }
1460 
1461 
1462  return IRRT_OK;
1463 }
1464 
1465 
1466 double
1467 FCMMaterial :: give(int aProperty, GaussPoint *gp)
1468 {
1469  return linearElasticMaterial->give(aProperty, gp);
1470 }
1471 
1472 
1473 
1474 double
1476 {
1477  return crackSpacing;
1478 }
1479 
1480 
1481 double
1483 {
1484  FCMMaterialStatus *status = static_cast< FCMMaterialStatus * >( this->giveStatus(gp) );
1485  double spacing, L, N;
1486 
1487  L = status->giveCharLength(iCrack);
1488  spacing = this->giveCrackSpacing();
1489 
1490  if ( ( spacing > L ) || ( spacing < 0. ) ) {
1491  N = 1;
1492  } else {
1493  N = ( L / spacing );
1494  }
1495 
1496  return N;
1497 }
1498 
1499 
1500 double
1502 {
1503  double N;
1504  int dir_1, dir_2;
1505 
1506  if ( i == 4 ) {
1507  dir_1 = 2;
1508  dir_2 = 3;
1509  } else if ( i == 5 ) {
1510  dir_1 = 1;
1511  dir_2 = 3;
1512  } else if ( i == 6 ) {
1513  dir_1 = 1;
1514  dir_2 = 2;
1515  } else {
1516  OOFEM_ERROR("Unexpected value of index i (4, 5, 6 permitted only)");
1517  dir_1 = dir_2 = 0;
1518  }
1519 
1520  N = max( this->giveNumberOfCracksInDirection(gp, dir_1), this->giveNumberOfCracksInDirection(gp, dir_2) );
1521 
1522  return N;
1523 }
1524 
1525 
1526 int
1528 {
1529  FCMMaterialStatus *status = static_cast< FCMMaterialStatus * >( this->giveStatus(gp) );
1530  double width;
1531  int index;
1532 
1533 
1534  // first/dominant crack vector
1535  if ( type == IST_CrackVector ) {
1536  answer.resize(3);
1537  answer.zero();
1538 
1539  // MAX CRACK
1540  width = 0.;
1541  index = 1;
1542  for ( int i = 1; i <= status->giveNumberOfCracks(); i++ ) {
1543  if ( status->giveCharLength(i) * status->giveCrackStrain(i) / this->giveNumberOfCracksInDirection(gp, i) > width ) {
1544  width = status->giveCharLength(i) * status->giveCrackStrain(i) / this->giveNumberOfCracksInDirection(gp, i);
1545  index = i;
1546  }
1547  }
1548 
1549  for ( int i = 1; i <= status->giveCrackDirs().giveNumberOfRows(); i++ ) {
1550  answer.at(i) = status->giveCrackDirs().at(i, index);
1551  }
1552 
1553  return 1;
1554  } else if ( type == IST_2ndCrackVector ) {
1555  answer.resize(3);
1556  answer.zero();
1557  index = 2;
1558 
1559  if ( status->giveNumberOfCracks() <= 1 ) {
1560  return 1;
1561  } else if ( status->giveNumberOfCracks() == 2 ) {
1562  if ( status->giveCharLength(2) * status->giveCrackStrain(2) / this->giveNumberOfCracksInDirection(gp, 2) > status->giveCharLength(1) * status->giveCrackStrain(1) / this->giveNumberOfCracksInDirection(gp, 1) ) {
1563  index = 1;
1564  }
1565  } else { // 3 cracks
1566  if ( status->giveCharLength(2) * status->giveCrackStrain(2) / this->giveNumberOfCracksInDirection(gp, 2) > status->giveCharLength(1) * status->giveCrackStrain(1) / this->giveNumberOfCracksInDirection(gp, 1) ) { // #2 is biggest
1567  if ( status->giveCharLength(3) * status->giveCrackStrain(3) / this->giveNumberOfCracksInDirection(gp, 3) > status->giveCharLength(1) * status->giveCrackStrain(1) / this->giveNumberOfCracksInDirection(gp, 1) ) {
1568  index = 3;
1569  } else {
1570  index = 1;
1571  }
1572  } else { // #1 is biggest
1573  if ( status->giveCharLength(3) * status->giveCrackStrain(3) / this->giveNumberOfCracksInDirection(gp, 3) > status->giveCharLength(2) * status->giveCrackStrain(2) / this->giveNumberOfCracksInDirection(gp, 2) ) {
1574  index = 3;
1575  } else {
1576  index = 2;
1577  }
1578  }
1579  }
1580 
1581  for ( int i = 1; i <= status->giveCrackDirs().giveNumberOfRows(); i++ ) {
1582  answer.at(i) = status->giveCrackDirs().at(i, index);
1583  }
1584 
1585  return 1;
1586  } else if ( type == IST_3rdCrackVector ) {
1587  answer.resize(3);
1588  answer.zero();
1589  index = 3;
1590 
1591  if ( status->giveNumberOfCracks() <= 2 ) {
1592  return 1;
1593  } else { // 3 cracks
1594  if ( status->giveCharLength(2) * status->giveCrackStrain(2) / this->giveNumberOfCracksInDirection(gp, 2) > status->giveCharLength(1) * status->giveCrackStrain(1) / this->giveNumberOfCracksInDirection(gp, 1) ) { // #2 is biggest
1595  if ( status->giveCharLength(3) * status->giveCrackStrain(3) / this->giveNumberOfCracksInDirection(gp, 3) > status->giveCharLength(1) * status->giveCrackStrain(1) / this->giveNumberOfCracksInDirection(gp, 1) ) {
1596  index = 1;
1597  } else {
1598  index = 3;
1599  }
1600  } else { // #1 is biggest
1601  if ( status->giveCharLength(3) * status->giveCrackStrain(3) / this->giveNumberOfCracksInDirection(gp, 3) > status->giveCharLength(2) * status->giveCrackStrain(2) / this->giveNumberOfCracksInDirection(gp, 2) ) {
1602  index = 2;
1603  } else {
1604  index = 3;
1605  }
1606  }
1607  }
1608 
1609  for ( int i = 1; i <= status->giveCrackDirs().giveNumberOfRows(); i++ ) {
1610  answer.at(i) = status->giveCrackDirs().at(i, index);
1611  }
1612 
1613  return 1;
1614 
1615 
1616  // width of a first/dominant crack
1617  } else if ( type == IST_CrackWidth ) {
1618  answer.resize(1);
1619  answer.zero();
1620 
1621  // MAX WIDTH
1622  width = 0.;
1623  for ( int i = 1; i <= status->giveNumberOfCracks(); i++ ) {
1624  width = max( width, status->giveCharLength(i) * status->giveCrackStrain(i) / this->giveNumberOfCracksInDirection(gp, i) );
1625  }
1626  answer.at(1) = width;
1627  return 1;
1628  } else if ( type == IST_2ndCrackWidth ) {
1629  answer.resize(1);
1630  answer.zero();
1631  index = 2;
1632  if ( status->giveNumberOfCracks() <= 1 ) {
1633  return 1;
1634  } else if ( status->giveNumberOfCracks() == 2 ) {
1635  if ( status->giveCharLength(2) * status->giveCrackStrain(2) / this->giveNumberOfCracksInDirection(gp, 2) > status->giveCharLength(1) * status->giveCrackStrain(1) / this->giveNumberOfCracksInDirection(gp, 1) ) {
1636  index = 1;
1637  }
1638  } else { // 3 cracks
1639  if ( status->giveCharLength(2) * status->giveCrackStrain(2) / this->giveNumberOfCracksInDirection(gp, 2) > status->giveCharLength(1) * status->giveCrackStrain(1) / this->giveNumberOfCracksInDirection(gp, 1) ) { // #2 is biggest
1640  if ( status->giveCharLength(3) * status->giveCrackStrain(3) / this->giveNumberOfCracksInDirection(gp, 3) > status->giveCharLength(1) * status->giveCrackStrain(1) / this->giveNumberOfCracksInDirection(gp, 1) ) {
1641  index = 3;
1642  } else {
1643  index = 1;
1644  }
1645  } else { // #1 is biggest
1646  if ( status->giveCharLength(3) * status->giveCrackStrain(3) / this->giveNumberOfCracksInDirection(gp, 3) > status->giveCharLength(2) * status->giveCrackStrain(2) / this->giveNumberOfCracksInDirection(gp, 2) ) {
1647  index = 3;
1648  } else {
1649  index = 2;
1650  }
1651  }
1652  }
1653 
1654  width = status->giveCharLength(index) * status->giveCrackStrain(index) / this->giveNumberOfCracksInDirection(gp, index);
1655 
1656  answer.at(1) = width;
1657  return 1;
1658  } else if ( type == IST_CrackDirs ) {
1659  const FloatMatrix &help = status->giveCrackDirs();
1660  answer.resize(9);
1661  for ( int i = 1; i <= 3; i++ ) {
1662  answer.at(i) = help.at(1, i);
1663  answer.at(i + 3) = help.at(2, i);
1664  answer.at(i + 6) = help.at(3, i);
1665  }
1666 
1667  return 1;
1668  } else if ( type == IST_CrackStatuses ) {
1669  const IntArray &crackStatus = status->giveCrackStatus();
1670  answer.resize(3);
1671  for ( int i = 1; i <= 3; i++ ) {
1672  answer.at(i) = crackStatus.at(i);
1673  }
1674 
1675  return 1;
1676  } else if ( type == IST_CrackStrainTensor ) {
1677  FloatArray crackStrain = status->giveCrackStrainVector();
1679  // from local to global
1680  crackStrain.rotatedWith(epsL2G, 'n');
1681 
1682  StructuralMaterial :: giveFullSymVectorForm( answer, crackStrain, gp->giveMaterialMode() );
1683 
1684  return 1;
1685  } else {
1686  return StructuralMaterial :: giveIPValue(answer, gp, type, tStep);
1687  }
1688 }
1689 
1690 void
1692  MatResponseMode mode,
1693  GaussPoint *gp,
1694  TimeStep *tStep)
1695 {
1696  this->giveMaterialStiffnessMatrix(answer, mode, gp, tStep);
1697 }
1698 
1699 
1700 void
1702  MatResponseMode mode,
1703  GaussPoint *gp,
1704  TimeStep *tStep)
1705 
1706 {
1707  this->giveMaterialStiffnessMatrix(answer, mode, gp, tStep);
1708 }
1709 
1710 
1711 void
1713  MatResponseMode mode,
1714  GaussPoint *gp,
1715  TimeStep *tStep)
1716 
1717 {
1718  this->giveMaterialStiffnessMatrix(answer, mode, gp, tStep);
1719 }
1720 
1721 
1723  StructuralMaterialStatus(n, d, gp),
1724  crackStatuses(), tempCrackStatuses(),
1725  maxCrackStrains(), tempMaxCrackStrains(),
1726  crackStrainVector(), tempCrackStrainVector(),
1727  crackDirs(),
1728  charLengths(),
1729  transMatrix_G2Lstress(), transMatrix_G2Lstrain(),
1730  transMatrix_L2Gstress(), transMatrix_L2Gstrain()
1731 {
1732  // resize in constructor according to stress-state
1733  this->nMaxCracks = 0;
1734  this->nMaxCracks = giveMaxNumberOfCracks(gp);
1735 
1737  crackStatuses.zero();
1739 
1740  charLengths.resize(this->nMaxCracks);
1741  charLengths.zero();
1742 
1743  crackDirs.resize(this->nMaxCracks, this->nMaxCracks);
1744  crackDirs.zero();
1745  for ( int i = 1; i <= this->nMaxCracks; i++ ) {
1746  crackDirs.at(i, i) = 1.0;
1747  }
1748 
1749 
1750  if ( this->nMaxCracks == 2 ) { //plane stress
1754 
1757 
1758 
1762  } else {
1766 
1769 
1770 
1774  }
1775 }
1776 
1777 
1779 { }
1780 
1781 
1782 
1783 void
1785 {
1786  int i;
1787  char s [ 11 ];
1788 
1790  fprintf(file, "status { ");
1791  if ( this->giveNumberOfCracks() > 0 ) {
1792  for ( i = 1; i <= crackDirs.giveNumberOfColumns(); i++ ) {
1793  switch ( crackStatuses.at(i) ) {
1794  case pscm_NONE:
1795  strcpy(s, "NONE");
1796  break;
1797  case pscm_JUST_INIT:
1798  strcpy(s, "JUST_INIT");
1799  break;
1800  case pscm_SOFTENING:
1801  strcpy(s, "SOFTENING");
1802  break;
1803  case pscm_UNLO_RELO:
1804  strcpy(s, "UNLO_RELO");
1805  break;
1806  case pscm_CLOSED:
1807  strcpy(s, "CLOSED");
1808  break;
1809  default:
1810  strcpy(s, "UNKNOWN");
1811  break;
1812  }
1813 
1814  fprintf(file, "crack %d {status %s, crackplane is normal to { ", i, s);
1815 
1816  for ( int j = 1; j <= crackDirs.giveNumberOfRows(); j++ ) {
1817  fprintf( file, "%f ", crackDirs.at(j, i) );
1818  }
1819 
1820  fprintf(file, "}}");
1821  ;
1822  }
1823  }
1824 
1825  fprintf(file, " }\n");
1826 }
1827 
1828 int
1830 //
1831 // return number of existing cracks
1832 //
1833 {
1834  int answer = 0;
1835 
1836  for ( int i = 1; i <= crackStatuses.giveSize(); i++ ) {
1837  if ( crackStatuses.at(i) != pscm_NONE ) {
1838  answer++;
1839  }
1840  }
1841 
1842  return answer;
1843 }
1844 
1845 
1846 int
1848 //
1849 // return number of existing temp cracks
1850 //
1851 {
1852  int answer = 0;
1853 
1854  for ( int i = 1; i <= tempCrackStatuses.giveSize(); i++ ) {
1855  if ( tempCrackStatuses.at(i) != pscm_NONE ) {
1856  answer++;
1857  }
1858  }
1859 
1860  return answer;
1861 }
1862 
1863 int
1865 //
1866 // return number of maximum allowable cracks
1867 //
1868 {
1869  int nCr = 0;
1870  IntArray indx;
1871 
1872  if ( this->nMaxCracks == 0 ) {
1874 
1875  for ( int i = 1; i <= 3; i++ ) {
1876  if ( indx.contains(i) ) {
1877  nCr++;
1878  }
1879  }
1880 
1881  this->nMaxCracks = nCr;
1882  }
1883 
1884  return this->nMaxCracks;
1885 }
1886 
1887 
1888 
1889 
1890 void
1892 //
1893 // initializes temp variables according to variables form previous equlibrium state.
1894 //
1895 {
1897 
1901 }
1902 
1903 
1904 
1905 void
1907 //
1908 // updates variables (nonTemp variables describing situation at previous equilibrium state)
1909 // after a new equilibrium state has been reached
1910 // temporary variables correspond to newly reched equilibrium.
1911 //
1912 {
1914 
1915 
1918 
1919  for ( int i = 1; i <= crackStrainVector.giveSize(); i++ ) {
1920  if ( crackStrainVector.at(i) < 0. ) {
1921  crackStrainVector.at(i) = 0.;
1922  }
1923  }
1924 
1925  // updating of statuses has to be done more carefully.
1926  // consider a crack which does not exist in the previous step and ends as closed in the end of this step
1927  // this crack is naturally treated as "NONE" in the following steps
1928 
1929  for ( int i = 1; i <= crackStatuses.giveSize(); i++ ) {
1930  if ( ( tempCrackStatuses.at(i) == pscm_CLOSED ) && ( crackStatuses.at(i) == pscm_NONE ) ) {
1931  // no other crack so this one can be set as non-existing
1932  if ( i + 1 > nMaxCracks ) {
1933  crackStatuses.at(i) = pscm_NONE;
1934 
1935 
1936  // be sure that in the second and third crack does not exist, if it does we have to copy CLOSED status
1937  } else if ( tempCrackStatuses.at(i + 1) != pscm_NONE ) {
1939  } else {
1940  crackStatuses.at(i) = pscm_NONE;
1941  }
1942  } else {
1944  }
1945  }
1946 }
1947 
1948 
1949 
1952 //
1953 // saves full information stored in this Status
1954 // no temp variables stored
1955 //
1956 {
1957  contextIOResultType iores;
1958 
1959  // save parent class status
1960  if ( ( iores = StructuralMaterialStatus :: saveContext(stream, mode, obj) ) != CIO_OK ) {
1961  THROW_CIOERR(iores);
1962  }
1963 
1964  // write a raw data
1965 
1966  if ( ( iores = crackStatuses.storeYourself(stream) ) != CIO_OK ) {
1967  THROW_CIOERR(iores);
1968  }
1969 
1970  if ( ( iores = maxCrackStrains.storeYourself(stream) ) != CIO_OK ) {
1971  THROW_CIOERR(iores);
1972  }
1973 
1974  if ( ( iores = crackDirs.storeYourself(stream) ) != CIO_OK ) {
1975  THROW_CIOERR(iores);
1976  }
1977 
1978  if ( ( iores = charLengths.storeYourself(stream) ) != CIO_OK ) {
1979  THROW_CIOERR(iores);
1980  }
1981 
1982  if ( ( iores = crackStrainVector.storeYourself(stream) ) != CIO_OK ) {
1983  THROW_CIOERR(iores);
1984  }
1985 
1986  if ( ( iores = crackStrainVector.storeYourself(stream) ) != CIO_OK ) {
1987  THROW_CIOERR(iores);
1988  }
1989 
1990  if ( ( iores = transMatrix_G2Lstrain.storeYourself(stream) ) != CIO_OK ) {
1991  THROW_CIOERR(iores);
1992  }
1993 
1994  if ( ( iores = transMatrix_G2Lstress.storeYourself(stream) ) != CIO_OK ) {
1995  THROW_CIOERR(iores);
1996  }
1997 
1998  if ( ( iores = transMatrix_L2Gstrain.storeYourself(stream) ) != CIO_OK ) {
1999  THROW_CIOERR(iores);
2000  }
2001 
2002  if ( ( iores = transMatrix_L2Gstress.storeYourself(stream) ) != CIO_OK ) {
2003  THROW_CIOERR(iores);
2004  }
2005 
2006  return CIO_OK;
2007 }
2008 
2011 //
2012 // restores full information stored in stream to this Status
2013 //
2014 {
2015  contextIOResultType iores;
2016 
2017  // read parent class status
2018  if ( ( iores = StructuralMaterialStatus :: restoreContext(stream, mode, obj) ) != CIO_OK ) {
2019  THROW_CIOERR(iores);
2020  }
2021 
2022  // read raw data
2023 
2024  if ( ( iores = crackStatuses.restoreYourself(stream) ) != CIO_OK ) {
2025  THROW_CIOERR(iores);
2026  }
2027 
2028  if ( ( iores = maxCrackStrains.restoreYourself(stream) ) != CIO_OK ) {
2029  THROW_CIOERR(iores);
2030  }
2031 
2032  if ( ( iores = crackDirs.restoreYourself(stream) ) != CIO_OK ) {
2033  THROW_CIOERR(iores);
2034  }
2035 
2036  if ( ( iores = charLengths.restoreYourself(stream) ) != CIO_OK ) {
2037  THROW_CIOERR(iores);
2038  }
2039 
2040  if ( ( iores = crackStrainVector.restoreYourself(stream) ) != CIO_OK ) {
2041  THROW_CIOERR(iores);
2042  }
2043 
2044  if ( ( iores = transMatrix_G2Lstrain.restoreYourself(stream) ) != CIO_OK ) {
2045  THROW_CIOERR(iores);
2046  }
2047 
2048  if ( ( iores = transMatrix_G2Lstress.restoreYourself(stream) ) != CIO_OK ) {
2049  THROW_CIOERR(iores);
2050  }
2051 
2052  if ( ( iores = transMatrix_L2Gstrain.restoreYourself(stream) ) != CIO_OK ) {
2053  THROW_CIOERR(iores);
2054  }
2055 
2056  if ( ( iores = transMatrix_L2Gstress.restoreYourself(stream) ) != CIO_OK ) {
2057  THROW_CIOERR(iores);
2058  }
2059 
2060  return CIO_OK; // return succes
2061 }
2062 } // end namespace oofem
#define pscm_JUST_INIT
Definition: fcm.h:55
bool contains(int value) const
Definition: intarray.h:283
void setG2LStressVectorTransformationMtrx(FloatMatrix t)
sets transformation matrix for stress transformation from global to local coordinate system ...
Definition: fcm.h:140
virtual void updateYourself(TimeStep *tStep)
Update equilibrium history variables according to temp-variables.
Definition: structuralms.C:96
InternalStateType
Type representing the physical meaning of element or constitutive model internal variable.
#define fcm_BIGNUMBER
Definition: fcm.h:61
MaterialMode giveMaterialMode()
Returns corresponding material mode of receiver.
Definition: gausspoint.h:191
contextIOResultType storeYourself(DataStream &stream) const
Stores array to output stream.
Definition: intarray.C:289
int giveNumberOfColumns() const
Returns number of columns of receiver.
Definition: floatmatrix.h:158
bool multipleCrackShear
if true = takes shear compliance of all cracks, false = only dominant crack contribution, default value is false
Definition: fcm.h:274
virtual double computeShearSlipOnCrack(GaussPoint *gp, int i)
computes total shear slip on a given crack plane (i = 1, 2, 3); the slip is computed from the tempora...
Definition: fcm.C:807
void subtract(const FloatArray &src)
Subtracts array src to receiver.
Definition: floatarray.C:258
static int giveVoigtSymVectorMask(IntArray &answer, MaterialMode mmode)
The same as giveVoigtVectorMask but returns a mask corresponding to a symmetric second order tensor...
virtual IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description stored in input record.
Definition: fcm.C:1424
void letTempStrainVectorBe(const FloatArray &v)
Assigns tempStrainVector to given vector v.
Definition: structuralms.h:137
virtual double computeOverallElasticStiffness(void)
returns overall Young&#39;s modulus
Definition: fcm.h:344
FCMMaterial(int n, Domain *d)
Definition: fcm.C:46
static void givePlaneStressVectorTranformationMtrx(FloatMatrix &answer, const FloatMatrix &base, bool transpose=false)
Computes 2d stress vector transformation matrix from standard vector transformation matrix...
FloatMatrix crackDirs
Storing direction of cracks (crack normals) in columwise format.
Definition: fcm.h:78
int nMaxCracks
number of maximum possible cracks (optional parameter)
Definition: fcm.h:92
virtual bool isStrengthExceeded(const FloatMatrix &base, GaussPoint *gp, TimeStep *tStep, int iCrack, double trialStress)
compares trial stress with strength. Returns true if the strength is exceeded. Function oveloaded in ...
Definition: fcm.C:914
virtual MaterialStatus * giveStatus(GaussPoint *gp) const
Returns material status of receiver in given integration point.
Definition: material.C:244
void setL2GStrainVectorTransformationMtrx(FloatMatrix s)
sets transformation matrix for stress transformation from global to local coordinate system ...
Definition: fcm.h:146
virtual void givePlaneStressStiffMtrx(FloatMatrix &answer, MatResponseMode mmode, GaussPoint *gp, TimeStep *tStep)
Method for computing plane stress stiffness matrix of receiver.
Definition: fcm.C:1701
virtual int giveIPValue(FloatArray &answer, GaussPoint *gp, InternalStateType type, TimeStep *tStep)
Returns the integration point corresponding value in Reduced form.
Definition: fcm.C:1527
void setMaxCrackStrain(int icrack, double val)
sets value of the maximum crack strain for the i-th crack (equilibrated value)
Definition: fcm.h:109
GaussPoint * gp
Associated integration point.
Class and object Domain.
Definition: domain.h:115
virtual void givePlaneStrainStiffMtrx(FloatMatrix &answer, MatResponseMode mmode, GaussPoint *gp, TimeStep *tStep)
Method for computing plane strain stiffness matrix of receiver.
Definition: fcm.C:1712
For computing principal stresses.
double crackSpacing
value of crack spacing (allows to "have" more parallel cracks in one direction if the element size ex...
Definition: fcm.h:322
#define _IFT_FCM_ecsm
Definition: fcm.h:49
virtual double giveCharacteristicElementLength(GaussPoint *gp, const FloatArray &crackPlaneNormal)
returns characteristic element length in given direction
Definition: fcm.C:1140
virtual bool isIntact(GaussPoint *gp, int icrack)
returns true for closed or no crack (i = 1, 2, 3)
Definition: fcm.C:720
bool solveForRhs(const FloatArray &b, FloatArray &answer, bool transpose=false)
Solves the system of linear equations .
Definition: floatmatrix.C:1112
int nAllowedCracks
allowed number of cracks (user-defined)
Definition: fcm.h:268
This class implements associated Material Status to FCMMaterial (fixed crack material).
Definition: fcm.h:68
void setTempCrackStatus(int icrack, int val)
sets temporary value of status for of the i-th crack
Definition: fcm.h:121
const FloatMatrix & giveG2LStressVectorTransformationMtrx()
returns transformation matrix for stress transformation from global to local coordinate system ...
Definition: fcm.h:149
virtual void give3dMaterialStiffnessMatrix(FloatMatrix &answer, MatResponseMode mode, GaussPoint *gp, TimeStep *tStep)
Computes full 3d material stiffness matrix at given integration point, time, respecting load history ...
Definition: fcm.C:1691
#define fcm_TOLERANCE
Definition: fcm.h:62
The purpose of DataStream abstract class is to allow to store/restore context to different streams...
Definition: datastream.h:54
virtual contextIOResultType saveContext(DataStream &stream, ContextMode mode, void *obj=NULL)
saves current context(state) into stream
Definition: fcm.C:1951
contextIOResultType storeYourself(DataStream &stream) const
Definition: floatarray.C:872
virtual int giveIPValue(FloatArray &answer, GaussPoint *gp, InternalStateType type, TimeStep *tStep)
Returns the integration point corresponding value in Reduced form.
virtual double maxShearStress(GaussPoint *gp, int i)=0
computes the maximum value of the shear stress; if the shear stress exceeds this value, it is cropped
void zero()
Sets all component to zero.
Definition: intarray.C:52
double & at(int i)
Coefficient access function.
Definition: floatarray.h:131
int max(int i, int j)
Returns bigger value form two given decimals.
Definition: mathfem.h:71
#define pscm_CLOSED
Definition: fcm.h:58
virtual double giveCrackSpacing(void)
returns either user-provided value of crack spacing or a value computed from composition ...
Definition: fcm.C:1475
virtual bool checkStrengthCriterion(FloatMatrix &newBase, const FloatArray &globalStress, GaussPoint *gp, TimeStep *tStep, int nCrack)
checks if the globalStress does not exceed strength in the direction of newBase for n-th crack ...
Definition: fcm.C:940
const FloatMatrix & giveG2LStrainVectorTransformationMtrx()
sets transformation matrix for strain transformation from global to local coordinate system ...
Definition: fcm.h:151
FloatArray charLengths
Characteristic lengths computed from the crack orientation and element geometry.
Definition: fcm.h:80
This class implements a structural material status information.
Definition: structuralms.h:65
#define pscm_NONE
Definition: fcm.h:54
FloatArray maxCrackStrains
Max. crack strain reached in the entire previous history.
Definition: fcm.h:74
virtual bool hasField(InputFieldType id)=0
Returns true if record contains field identified by idString keyword.
virtual int giveNumberOfCracks() const
returns number of cracks from the previous time step (equilibrated value)
Definition: fcm.C:1829
Element * giveElement()
Returns corresponding element to receiver.
Definition: gausspoint.h:188
FloatArray tempCrackStrainVector
Definition: fcm.h:76
const FloatMatrix & giveCrackDirs()
returns crack directions
Definition: fcm.h:168
double giveCharLength(int icrack) const
returns characteristic length associated with i-th crack direction
Definition: fcm.h:158
static void giveStrainVectorTranformationMtrx(FloatMatrix &answer, const FloatMatrix &base, bool transpose=false)
Computes 3d strain vector transformation matrix from standard vector transformation matrix...
ElementCharSizeMethod ecsMethod
Method used for evaluation of characteristic element size.
Definition: fcm.h:271
#define _IFT_FCM_nAllowedCracks
Definition: fcm.h:46
virtual void initializeCrack(GaussPoint *gp, FloatMatrix &base, int nCrack)
Definition: fcm.C:650
virtual ~FCMMaterial()
Definition: fcm.C:56
void giveStressDependentPartOfStrainVector(FloatArray &answer, GaussPoint *gp, const FloatArray &reducedStrainVector, TimeStep *tStep, ValueModeType mode)
Method for subtracting from reduced space strain vector its stress-independent parts (caused by tempe...
MaterialMode
Type representing material mode of integration point.
Definition: materialmode.h:89
virtual void printOutputAt(FILE *file, TimeStep *tStep)
Print receiver&#39;s output to given stream.
Definition: structuralms.C:73
static void computePrincipalValDir(FloatArray &answer, FloatMatrix &dir, const FloatArray &s, stressStrainPrincMode mode)
Computes principal values and directions of stress or strain vector.
virtual double computeD2ModulusForCrack(GaussPoint *gp, int icrack)=0
shear modulus for a given crack plane (1, 2, 3)
Class implementing an array of integers.
Definition: intarray.h:61
int & at(int i)
Coefficient access function.
Definition: intarray.h:103
MatResponseMode
Describes the character of characteristic material matrix.
virtual double computeMaxNormalCrackOpening(GaussPoint *gp, int i)
uses maximum equilibrated cracking strain and characteristic length to obtain the maximum reached cra...
Definition: fcm.C:788
#define THROW_CIOERR(e)
Definition: contextioerr.h:61
virtual double giveCrackingModulus(MatResponseMode rMode, GaussPoint *gp, int i)=0
returns stiffness in the normal direction of the i-th crack
virtual void updateCrackStatus(GaussPoint *gp)
updates crack statuses
Definition: fcm.C:1147
virtual void giveRealStressVector(FloatArray &answer, GaussPoint *gp, const FloatArray &reducedStrain, TimeStep *tStep)
Computes the real stress vector for given total strain and integration point.
Definition: fcm.C:72
contextIOResultType storeYourself(DataStream &stream) const
Definition: floatmatrix.C:1852
void rotatedWith(FloatMatrix &r, char mode)
Returns the receiver a rotated according the change-of-base matrix r.
Definition: floatarray.C:799
virtual double computeTotalD2Modulus(GaussPoint *gp, int i)
shear modulus for a given shear direction (4, 5, 6)
Definition: fcm.C:1378
void beDifferenceOf(const FloatArray &a, const FloatArray &b)
Sets receiver to be a - b.
Definition: floatarray.C:341
static void giveFullSymVectorForm(FloatArray &answer, const FloatArray &vec, MaterialMode matMode)
Converts the reduced unsymmetric Voigt vector (2nd order tensor) to full form.
const FloatArray & giveMaxCrackStrainVector()
returns vector with maximum cracking strains (max 3 components)
Definition: fcm.h:105
void setCrackDirs(FloatMatrix a)
sets matrix with crack directions (normal vectors)
Definition: fcm.h:172
virtual double give(int aProperty, GaussPoint *gp)
Returns the value of material property &#39;aProperty&#39;.
Definition: fcm.C:1467
virtual void initTempStatus()
initializes temporary status
Definition: fcm.C:1891
const char * __MaterialModeToString(MaterialMode _value)
Definition: cltypes.C:314
virtual IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description stored in input record.
virtual void giveStiffnessMatrix(FloatMatrix &answer, MatResponseMode mode, GaussPoint *gp, TimeStep *tStep)
Computes the stiffness matrix for giveRealStressVector of receiver in given integration point...
double giveTempCrackStrain(int icrack) const
returns i-th component of the crack strain vector (temporary)
Definition: fcm.h:132
virtual void checkSnapBack(GaussPoint *gp, int crack)=0
checks possible snap-back
IntArray tempCrackStatuses
Definition: fcm.h:72
const FloatArray & giveTempCrackStrainVector()
return temporary crack strain vector (max 6 components)
Definition: fcm.h:128
void setTempCrackStrain(int icrack, double val)
sets temporary value of i-th cracking strain (max 6 components)
Definition: fcm.h:136
void setG2LStrainVectorTransformationMtrx(FloatMatrix s)
sets transformation matrix for strain transformation from global to local coordinate system ...
Definition: fcm.h:142
#define OOFEM_ERROR(...)
Definition: error.h:61
const FloatMatrix & giveL2GStressVectorTransformationMtrx()
sets transformation matrix for stress transformation from local to global coordinate system ...
Definition: fcm.h:153
void resizeWithData(int, int)
Checks size of receiver towards requested bounds.
Definition: floatmatrix.C:1369
virtual int giveNumberOfTempCracks() const
returns temporary number of cracks
Definition: fcm.C:1847
void times(double f)
Multiplies receiver by factor f.
Definition: floatmatrix.C:1594
virtual double computeEffectiveShearModulus(GaussPoint *gp, int i)=0
returns Geff which is necessary in the global stiffness matrix
void rotatedWith(const FloatMatrix &r, char mode= 'n')
Returns the receiver &#39;a&#39; transformed using give transformation matrix r.
Definition: floatmatrix.C:1557
contextIOResultType restoreYourself(DataStream &stream)
Definition: floatarray.C:895
virtual contextIOResultType restoreContext(DataStream &stream, ContextMode mode, void *obj=NULL)
restores context(state) from stream
Definition: fcm.C:2010
void beProductOf(const FloatMatrix &aMatrix, const FloatArray &anArray)
Receiver becomes the result of the product of aMatrix and anArray.
Definition: floatarray.C:676
FloatArray tempMaxCrackStrains
Definition: fcm.h:74
#define N(p, q)
Definition: mdm.C:367
virtual int hasMaterialModeCapability(MaterialMode mode)
Tests if material supports material mode.
Definition: fcm.C:63
FloatMatrix transMatrix_G2Lstress
transformation matrix converting stress from global to local coordinate system
Definition: fcm.h:83
void beTProductOf(const FloatMatrix &aMatrix, const FloatArray &anArray)
Receiver becomes the result of the product of aMatrix^T and anArray.
Definition: floatarray.C:708
double at(int i, int j) const
Coefficient access function.
Definition: floatmatrix.h:176
void resize(int n)
Checks size of receiver towards requested bounds.
Definition: intarray.C:124
void setTempMaxCrackStrain(int icrack, double val)
sets value of the maximum crack strain for the i-th crack (temporary value)
Definition: fcm.h:114
virtual contextIOResultType restoreContext(DataStream &stream, ContextMode mode, void *obj=NULL)
Restores the receiver state previously written in stream.
Definition: structuralms.C:157
contextIOResultType restoreYourself(DataStream &stream)
Restores array from image on stream.
Definition: intarray.C:305
void setTempCrackStrainVector(FloatArray a)
sets temporary vector of cracking strains (max 6 components)
Definition: fcm.h:134
virtual void giveMaterialStiffnessMatrix(FloatMatrix &answer, MatResponseMode, GaussPoint *gp, TimeStep *tStep)
Definition: fcm.C:1250
void subtract(const FloatMatrix &a)
Subtracts matrix from the receiver.
Definition: floatmatrix.C:1084
contextIOResultType restoreYourself(DataStream &stream)
Definition: floatmatrix.C:1877
Class representing vector of real numbers.
Definition: floatarray.h:82
virtual double giveNumberOfCracksForShearDirection(GaussPoint *gp, int i)
returns number of cracks for given shear direction (i = 4, 5, 6) which is treated as the maximum of t...
Definition: fcm.C:1501
const IntArray & giveTempCrackStatus()
returns vector of temporary crack statuses
Definition: fcm.h:117
Implementation of matrix containing floating point numbers.
Definition: floatmatrix.h:94
const FloatArray & giveCrackStrainVector() const
return equilibrated crack strain vector (max 6 components)
Definition: fcm.h:126
IRResultType
Type defining the return values of InputRecord reading operations.
Definition: irresulttype.h:47
static void giveStressVectorTranformationMtrx(FloatMatrix &answer, const FloatMatrix &base, bool transpose=false)
Computes 3d stress vector transformation matrix from standard vector transformation matrix...
void letTempStressVectorBe(const FloatArray &v)
Assigns tempStressVector to given vector v.
Definition: structuralms.h:135
const FloatMatrix & giveL2GStrainVectorTransformationMtrx()
sets transformation matrix for stress transformation from global to local coordinate system ...
Definition: fcm.h:155
const FloatArray & giveStressVector() const
Returns the const pointer to receiver&#39;s stress vector.
Definition: structuralms.h:107
#define pscm_UNLO_RELO
Definition: fcm.h:57
void resize(int rows, int cols)
Checks size of receiver towards requested bounds.
Definition: floatmatrix.C:1358
Class representing the general Input Record.
Definition: inputrecord.h:101
#define _IFT_FCM_multipleCrackShear
Definition: fcm.h:48
void add(const FloatMatrix &a)
Adds matrix to the receiver.
Definition: floatmatrix.C:1023
IntArray crackStatuses
crack statuses (none, just initialized, softenin, unlo-relo, closed)
Definition: fcm.h:72
void setL2GStressVectorTransformationMtrx(FloatMatrix t)
sets transformation matrix for stress transformation from local to global coordinate system ...
Definition: fcm.h:144
IsotropicLinearElasticMaterial * linearElasticMaterial
Definition: fcm.h:200
double giveMaxCrackStrain(int icrack)
returns maximum crack strain for the i-th crack (equilibrated value)
Definition: fcm.h:107
double giveCrackStrain(int icrack) const
returns i-th component of the crack strain vector (equilibrated)
Definition: fcm.h:130
virtual bool isIntactForShear(GaussPoint *gp, int i)
returns true for closed or no cracks associated to given shear direction (i = 4, 5, 6)
Definition: fcm.C:738
void zero()
Zeroes all coefficients of receiver.
Definition: floatarray.C:658
virtual double giveCharacteristicLength(const FloatArray &normalToCrackPlane)
Returns the size of element in the given direction, in some cases adjusted (e.g.
Definition: element.h:874
virtual int giveMaxNumberOfCracks(GaussPoint *gp)
returns maximum number of cracks associated with current mode
Definition: fcm.C:1864
virtual double giveTensileStrength(GaussPoint *gp)=0
comutes tensile strength
virtual double giveNumberOfCracksInDirection(GaussPoint *gp, int iCrack)
returns number of fictiotious parallel cracks in the direction of i-th crack
Definition: fcm.C:1482
virtual void printOutputAt(FILE *file, TimeStep *tStep)
prints the output into the output file
Definition: fcm.C:1784
virtual void updateYourself(TimeStep *tStep)
replaces equilibrated values with temporary values
Definition: fcm.C:1906
void times(double s)
Multiplies receiver with scalar.
Definition: floatarray.C:818
long ContextMode
Context mode (mask), defining the type of information written/read to/from context.
Definition: contextmode.h:43
Abstract base class for all "structural" constitutive models.
void zero()
Zeroes all coefficient of receiver.
Definition: floatmatrix.C:1326
virtual void initTempStatus()
Initializes the temporary internal variables, describing the current state according to previously re...
Definition: structuralms.C:108
FloatArray crackStrainVector
Components of crack strain vector (normal as well as shear).
Definition: fcm.h:76
static void give2DStrainVectorTranformationMtrx(FloatMatrix &answer, const FloatMatrix &base, bool transpose=false)
Computes 2d strain vector transformation matrix from standard vector transformation matrix...
int min(int i, int j)
Returns smaller value from two given decimals.
Definition: mathfem.h:59
virtual double give(int aProperty, GaussPoint *gp)
Returns the value of material property &#39;aProperty&#39;.
#define IR_GIVE_OPTIONAL_FIELD(__ir, __value, __id)
Macro facilitating the use of input record reading methods.
Definition: inputrecord.h:78
void beProductOf(const FloatMatrix &a, const FloatMatrix &b)
Assigns to the receiver product of .
Definition: floatmatrix.C:337
void setCharLength(int icrack, double val)
sets characteristic length for i-th crack
Definition: fcm.h:166
FloatMatrix transMatrix_G2Lstrain
transformation matrix converting strain from global to local coordinate system
Definition: fcm.h:85
int giveSize() const
Definition: intarray.h:203
FloatMatrix transMatrix_L2Gstrain
transformation matrix converting strain from local to global coordinate system
Definition: fcm.h:89
int giveSize() const
Returns the size of receiver.
Definition: floatarray.h:218
virtual double computeNormalCrackOpening(GaussPoint *gp, int i)
uses temporary cracking strain and characteristic length to obtain the crack opening ...
Definition: fcm.C:769
virtual contextIOResultType saveContext(DataStream &stream, ContextMode mode, void *obj=NULL)
Stores receiver state to output stream.
Definition: structuralms.C:133
the oofem namespace is to define a context or scope in which all oofem names are defined.
void setCrackStrainVector(FloatArray a)
sets equilibrated vector of cracking strains (max 6 components)
Definition: fcm.h:138
#define IR_GIVE_FIELD(__ir, __value, __id)
Macro facilitating the use of input record reading methods.
Definition: inputrecord.h:69
virtual void giveLocalCrackedStiffnessMatrix(FloatMatrix &answer, MatResponseMode rMode, GaussPoint *gp, TimeStep *tStep)
returns local stiffness matrix of the crack
Definition: fcm.C:1219
virtual void initTempStatus(GaussPoint *gp)
Initializes temporary variables stored in integration point status at the beginning of new time step...
Definition: material.C:267
void beInverseOf(const FloatMatrix &src)
Modifies receiver to become inverse of given parameter.
Definition: floatmatrix.C:835
const FloatArray & giveStrainVector() const
Returns the const pointer to receiver&#39;s strain vector.
Definition: structuralms.h:105
#define pscm_SOFTENING
Definition: fcm.h:56
int giveNumberOfRows() const
Returns number of rows of receiver.
Definition: floatmatrix.h:156
FCMMaterialStatus(int n, Domain *d, GaussPoint *g)
Definition: fcm.C:1722
int giveCrackStatus(int icrack) const
return equilibrated value of status associated with i-th crack direction
Definition: fcm.h:123
#define fcm_SMALL_STRAIN
Definition: fcm.h:60
Class representing integration point in finite element program.
Definition: gausspoint.h:93
#define OOFEM_WARNING(...)
Definition: error.h:62
Class representing solution step.
Definition: timestep.h:80
void add(const FloatArray &src)
Adds array src to receiver.
Definition: floatarray.C:156
IsotropicLinearElasticMaterial * giveLinearElasticMaterial()
Definition: fcm.h:216
static void giveReducedSymMatrixForm(FloatMatrix &answer, const FloatMatrix &full, MaterialMode matMode)
Converts the full unsymmetric Voigt matrix (4th order tensor) to reduced form.
FloatMatrix transMatrix_L2Gstress
transformation matrix converting stress from local to global coordinate system
Definition: fcm.h:87
#define _IFT_FCM_crackSpacing
Definition: fcm.h:47
void resize(int s)
Resizes receiver towards requested size.
Definition: floatarray.C:631
virtual ~FCMMaterialStatus()
Definition: fcm.C:1778
virtual double computeShearStiffnessRedistributionFactor(GaussPoint *gp, int ithCrackPlane, int jthCrackDirection)
function calculating ratio used to split shear slips on two crack planes
Definition: fcm.C:926

This page is part of the OOFEM documentation. Copyright (c) 2011 Borek Patzak
Project e-mail: info@oofem.org
Generated at Tue Jan 2 2018 20:07:28 for OOFEM by doxygen 1.8.11 written by Dimitri van Heesch, © 1997-2011