OOFEM  2.4
OOFEM.org - Object Oriented Finite Element Solver
mpsdammat.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 "mpsdammat.h"
36 #include "mathfem.h"
37 #include "gausspoint.h"
38 #include "timestep.h"
39 #include "stressvector.h"
40 #include "engngm.h"
41 #include "contextioerr.h"
42 #include "datastream.h"
43 #include "classfactory.h"
44 
45 namespace oofem {
46 REGISTER_Material(MPSDamMaterial);
47 
48 /***************************************************************/
49 /************** MPSDamMaterialStatus *******************/
50 /***************************************************************/
51 
52 
53 
55  MPSMaterialStatus(n, d, g, nunits), effectiveStressVector(), tempEffectiveStressVector()
56 {
57  kappa = tempKappa = 0.0;
58  damage = tempDamage = 0.0;
59  charLength = 0.0;
61  crackVector.zero();
62 
63  var_e0 = var_gf = 0.;
64 
69 
70 #ifdef supplementary_info
71  crackWidth = 0.;
73 #endif
74 }
75 
76 
77 void
79 {
81 
82  this->tempKappa = this->kappa;
83  this->tempDamage = this->damage;
84 
86 
87  if ( !damage ) {
88  var_e0 = var_gf = 0.;
89  }
90 }
91 
92 void
94 {
96 
97  this->kappa = this->tempKappa;
98  this->damage = this->tempDamage;
99 
101 
102  if ( !damage ) {
103  var_e0 = var_gf = 0.;
104  }
105 }
106 
107 void
109 {
110  answer = crackVector;
111  answer.times(damage);
112 }
113 
114 
115 void
117 {
119 
120  fprintf(file, "damage status { ");
121  if ( this->kappa > 0 && this->damage <= 0 ) {
122  fprintf(file, "kappa %f", this->kappa);
123  } else if ( this->damage > 0.0 ) {
124  fprintf( file, "kappa %f, damage %f crackVector %f %f %f", this->kappa, this->damage, this->crackVector.at(1), this->crackVector.at(2), this->crackVector.at(3) );
125  }
126  fprintf(file, "}\n");
127 }
128 
129 
130 
133 //
134 // saves full information stored in this Status
135 //
136 {
137  contextIOResultType iores;
138 
139  if ( ( iores = MPSMaterialStatus :: saveContext(stream, mode, obj) ) != CIO_OK ) {
140  THROW_CIOERR(iores);
141  }
142 
143  if ( !stream.write(kappa) ) {
145  }
146 
147  if ( !stream.write(damage) ) {
149  }
150 
151  if ( !stream.write(charLength) ) {
153  }
154 
155  if ( ( iores = crackVector.storeYourself(stream) ) != CIO_OK ) {
156  THROW_CIOERR(iores);
157  }
158 
159  if ( !stream.write(var_e0) ) {
161  }
162 
163  if ( !stream.write(var_gf) ) {
165  }
166 
167  if ( ( iores = effectiveStressVector.storeYourself(stream) ) != CIO_OK ) {
168  THROW_CIOERR(iores);
169  }
170 
171  return CIO_OK;
172 }
173 
176 //
177 // restore the state variables from a stream
178 //
179 {
180  contextIOResultType iores;
181  if ( ( iores = MPSMaterialStatus :: restoreContext(stream, mode, obj) ) != CIO_OK ) {
182  THROW_CIOERR(iores);
183  }
184 
185  if ( !stream.read(kappa) ) {
186  return CIO_IOERR;
187  }
188 
189  if ( !stream.read(damage) ) {
190  return CIO_IOERR;
191  }
192 
193  if ( !stream.read(charLength) ) {
194  return CIO_IOERR;
195  }
196 
197  if ( ( iores = crackVector.restoreYourself(stream) ) != CIO_OK ) {
198  THROW_CIOERR(iores);
199  }
200 
201  if ( !stream.read(var_e0) ) {
202  return CIO_IOERR;
203  }
204 
205  if ( !stream.read(var_gf) ) {
206  return CIO_IOERR;
207  }
208 
209  if ( ( iores = effectiveStressVector.restoreYourself(stream) ) != CIO_OK ) {
210  THROW_CIOERR(iores);
211  }
212 
213 
214  return CIO_OK;
215 }
216 
217 // ***********************************************************************************
218 // *** CLASS Damage extension of MPS model for creep and shrinakge of concrete ***
219 // ***********************************************************************************
220 
221 
223 
224 {
225  maxOmega = 0.999999;
226 
229  // const_e0 = 0.;
230  const_gf = 0.;
231  checkSnapBack = 1; //snapback check by default
232 
233  E = -1.;
234 
235 }
236 
237 
238 int
240 //
241 // returns whether receiver supports given mode
242 //
243 {
244  return mode == _3dMat || mode == _PlaneStress || mode == _PlaneStrain || mode == _1dMat;
245 }
246 
247 
248 
251 {
252  IRResultType result; // Required by IR_GIVE_FIELD macro
253 
255 
256  this->isotropic = false;
258  this->isotropic = true;
259  }
260 
261  int damageLaw = 0;
263 
264  this->timeDepFracturing = false;
265 
267  this->timeDepFracturing = true;
268  //
270  // the same compressive strength as for the prediction using the B3 formulas
271 
272  this->gf28 = 0.;
273  this->ft28 = 0.;
274 
276 
278  if (gf28 < 0.) {
279  OOFEM_ERROR("Fracture energy at 28 days must be positive");
280  }
282 
283  if (ft28 < 0.) {
284  OOFEM_ERROR("Tensile strength at 28 days must be positive");
285  }
286 
287  } else {
291  }
292 
293  } else {
294 
295  //applies only in this class
296  switch ( damageLaw ) {
297 
298  case 0: // exponential softening - default
302  break;
303 
304  case 1: // linear softening law
308  break;
309 
310  case 6:
311  this->softType = ST_Disable_Damage;
312  break;
313 
314  default:
315  OOFEM_ERROR("Softening type number %d is unknown", damageLaw);
316  }
317  }
318 
320 
321  return IRRT_OK;
322 }
323 
324 
325 void
327 {
328 
329  if (this->E < 0.) { // initialize dummy elastic modulus E
330  this->E = 1. / MPSMaterial :: computeCreepFunction(28.01, 28., gp, tStep);
331  }
332 
333  MPSDamMaterialStatus *status = static_cast< MPSDamMaterialStatus * >( this->giveStatus(gp) );
334 
335  MaterialMode mode = gp->giveMaterialMode();
336 
337  FloatArray principalStress;
338  FloatMatrix principalDir;
339 
340  StressVector tempEffectiveStress(mode);
341  StressVector tempNominalStress(mode);
342 
343  double f, sigma1, kappa, tempKappa = 0.0, omega = 0.0;
344 
345  // effective stress computed by the viscoelastic MPS material
346  MPSMaterial :: giveRealStressVector(tempEffectiveStress, gp, totalStrain, tStep);
347 
348  if ( !this->isActivated(tStep) ) {
349  FloatArray help;
351  help.zero();
352 
353  status->letTempStrainVectorBe(help);
354  status->letTempStressVectorBe(help);
355  status->letTempViscoelasticStressVectorBe(help);
356 
357 #ifdef supplementary_info
358  status->setResidualTensileStrength(0.);
359  status->setCrackWidth(0.);
360  status->setTempKappa(0.);
361  status->setTempDamage(0.);
362 #endif
363 
364  return;
365  }
366 
367  tempEffectiveStress.computePrincipalValDir(principalStress, principalDir);
368 
369  sigma1 = 0.;
370  if ( principalStress.at(1) > 0.0 ) {
371  sigma1 = principalStress.at(1);
372  }
373 
374  kappa = sigma1 / this->E;
375 
376  f = kappa - status->giveKappa();
377 
378  if ( f <= 0.0 ) { // damage does not grow
379  tempKappa = status->giveKappa();
380  omega = status->giveDamage();
381 
382 #ifdef supplementary_info
383  double residualStrength = 0.;
384  double e0;
385 
386  if ( ( this->timeDepFracturing ) && ( this->givee0(gp) == 0. ) ) {
387  this->initDamagedFib(gp, tStep);
388  }
389 
390  e0 = this->givee0(gp);
391 
392  if ( omega == 0. ) {
393  residualStrength = E * e0; // undamaged material
394  } else {
395  double gf, ef, wf = 0., Le;
396  gf = this->givegf(gp);
397  Le = status->giveCharLength();
398 
399  if ( softType == ST_Exponential_Cohesive_Crack ) { // exponential softening
400  wf = gf / this->E / e0; // wf is the crack opening
401  } else if ( softType == ST_Linear_Cohesive_Crack ) { // linear softening law
402  wf = 2. * gf / this->E / e0; // wf is the crack opening
403  } else {
404  OOFEM_ERROR("Gf unsupported for softening type softType = %d", softType);
405  }
406 
407  ef = wf / Le; //ef is the fracturing strain
408 
409  if ( this->softType == ST_Linear_Cohesive_Crack ) {
410  residualStrength = E * e0 * ( ef - tempKappa ) / ( ef - e0 );
411  } else if ( this->softType == ST_Exponential_Cohesive_Crack ) {
412  residualStrength = E * e0 * exp(-1. * ( tempKappa - e0 ) / ef);
413  } else {
414  OOFEM_ERROR("Unknown softening type for cohesive crack model.");
415  }
416  }
417 
418  if ( status ) {
419  status->setResidualTensileStrength(residualStrength);
420  }
421 
422 #endif
423  } else {
424  // damage grows
425  tempKappa = kappa;
426 
427  FloatArray crackPlaneNormal(3);
428  for ( int i = 1; i <= principalStress.giveSize(); i++ ) {
429  crackPlaneNormal.at(i) = principalDir.at(i, 1);
430  }
431  this->initDamaged(tempKappa, crackPlaneNormal, gp, tStep);
432  this->computeDamage(omega, tempKappa, gp);
433  }
434 
435  answer.zero();
436 
437  if ( omega > 0. ) {
438  tempNominalStress = tempEffectiveStress;
439 
440  if ( this->isotropic ) {
441  // convert effective stress to nominal stress
442  tempNominalStress.times(1. - omega);
443  answer.add(tempNominalStress);
444  } else {
445  // stress transformation matrix
446  FloatMatrix Tstress;
447 
448  // compute principal nominal stresses by multiplying effective stresses by damage
449  for ( int i = 1; i <= principalStress.giveSize(); i++ ) {
450  if ( principalStress.at(i) > 0. ) {
451  // convert principal effective stress to nominal stress
452  principalStress.at(i) *= ( 1. - omega );
453  }
454  }
455 
456  if ( mode == _PlaneStress ) {
457  principalStress.resizeWithValues(3);
458  givePlaneStressVectorTranformationMtrx(Tstress, principalDir, true);
459  } else {
460  principalStress.resizeWithValues(6);
461  giveStressVectorTranformationMtrx(Tstress, principalDir, true);
462  }
463 
464  principalStress.rotatedWith(Tstress, 'n');
465 
466 
467  if ( mode == _PlaneStress ) { // plane stress
468  answer.add(principalStress);
469  } else if ( this->giveSizeOfVoigtSymVector(mode) != this->giveSizeOfVoigtSymVector(_3dMat) ) { // mode = _PlaneStrain or axial symmetry
470  StressVector redFormStress(mode);
471  redFormStress.convertFromFullForm(principalStress, mode);
472  answer.add(redFormStress);
473  } else { // 3D
474  answer.add(principalStress);
475  }
476  }
477  } else {
478  answer.add(tempEffectiveStress);
479  }
480 
481 #ifdef supplementary_info
482 
483  if ( ( omega == 0. ) || ( sigma1 <= 0 ) ) {
484  status->setCrackWidth(0.);
485  } else {
486  FloatArray principalStrains;
487  this->computePrincipalValues(principalStrains, totalStrain, principal_strain);
488 
489  double crackWidth;
490  //case when the strain localizes into narrow band and after a while the stresses relax
491 
492  double strainWithoutTemperShrink = principalStrains.at(1);
493  strainWithoutTemperShrink -= status->giveTempThermalStrain();
494  strainWithoutTemperShrink -= status->giveTempDryingShrinkageStrain();
495  strainWithoutTemperShrink -= status->giveTempAutogenousShrinkageStrain();
496 
497 
498  crackWidth = status->giveCharLength() * omega * strainWithoutTemperShrink;
499 
500  status->setCrackWidth(crackWidth);
501  }
502 
503 #endif
504 
505  // update gp
506  status->letTempStrainVectorBe(totalStrain);
507  status->letTempStressVectorBe(answer);
508  status->letTempViscoelasticStressVectorBe(tempEffectiveStress);
509  status->setTempKappa(tempKappa);
510  status->setTempDamage(omega);
511 }
512 
513 
514 void
516 {
517  MPSDamMaterialStatus *status = ( MPSDamMaterialStatus * ) this->giveStatus(gp);
518 
519  if ( status->giveDamage() == 0. ) {
520  double tequiv;
521 
522  tequiv = this->computeEquivalentTime(gp, tStep, 0);
523 
524  double e0 = this->computeTensileStrength(tequiv) / this->E;
525  double gf = this->computeFractureEnergy(tequiv);
526 
527  status->sete0(e0);
528  status->setgf(gf);
529  }
530 }
531 
532 double
534 {
535  if ( this->timeDepFracturing ) {
536  MPSDamMaterialStatus *status = ( MPSDamMaterialStatus * ) this->giveStatus(gp);
537  return status->givee0();
538  } else {
539  return this->ft / this->E;
540  }
541 }
542 
543 double
545 {
546  if ( this->timeDepFracturing ) {
547  MPSDamMaterialStatus *status = ( MPSDamMaterialStatus * ) this->giveStatus(gp);
548  return status->givegf();
549  } else {
550  return this->const_gf;
551  }
552 }
553 
554 
555 double
557 {
558 
559  double fractureEnergy, fractureEnergy28;
560  double ftm, ftm28;
561 
562  // the fracture energy has the same time evolution as the tensile strength,
563  // the direct relation to the mean value of compressive strength according to Model Code
564  // highly overestimates the initial (early age) value of the fracture energy
565  //( fractureEnergy = 73. * fcm(t)^0.18 )
566 
567  // evolution of the mean compressive strength with respect to equivalent time/age/maturity
568  // returns fcm in MPa
569 
570  /* if (this->gf28 > 0.) {
571  double fcm28mod;
572  fcm28mod = pow ( this->gf28 * MPSMaterial :: stiffnessFactor / 73. , 1. / 0.18 );
573  fcm = exp( fib_s * ( 1. - sqrt(28. * MPSMaterial :: lambda0 / equivalentTime) ) ) * fcm28mod;
574 
575  } else {
576  fcm = exp( fib_s * ( 1. - sqrt(28. * MPSMaterial :: lambda0 / equivalentTime) ) ) * fib_fcm28;
577  }
578 
579  fractureEnergy = 73. * pow(fcm, 0.18) / MPSMaterial :: stiffnessFactor;
580  */
581 
582  // 1) read or estimate the 28-day value of fracture energy
583 
584  if (this->gf28 > 0.) {
585  fractureEnergy28 = this->gf28;
586  } else {
587  fractureEnergy28 = 73. * pow(fib_fcm28, 0.18) / MPSMaterial :: stiffnessFactor;
588  }
589 
590  // 2) compute the tensile strengh according to provided equivalent time
591 
592  ftm = this->computeTensileStrength(equivalentTime);
593  ftm28 = this->computeTensileStrength(28. * MPSMaterial :: lambda0);
594 
595  // 3) calculate the resulting fracture energy as gf28 * ft/ft28
596 
597  fractureEnergy = fractureEnergy28 * ftm / ftm28;
598 
599  return fractureEnergy;
600 }
601 
602 double
604 {
605  double fcm, ftm;
606 
607 
608  if (this->ft28 > 0.) {
609  double fcm28mod;
610  fcm28mod = pow ( this->ft28 * MPSMaterial :: stiffnessFactor / 0.3e6, 3./2. ) + 8.;
611  fcm = exp( fib_s * ( 1. - sqrt(28. * MPSMaterial :: lambda0 / equivalentTime) ) ) * fcm28mod;
612 
613  } else {
614  // returns fcm in MPa - formula 5.1-51, Table 5.1-9
615  fcm = exp( fib_s * ( 1. - sqrt(28. * MPSMaterial :: lambda0 / equivalentTime) ) ) * fib_fcm28;
616  }
617 
618 
619  // ftm adjusted according to the stiffnessFactor (MPa by default)
620  if ( fcm >= 58. ) {
621  ftm = 2.12 * log ( 1. + 0.1 * fcm ) * 1.e6 / MPSMaterial :: stiffnessFactor;
622  } else if ( fcm <= 20. ) {
623  ftm = 0.07862 * fcm * 1.e6 / MPSMaterial :: stiffnessFactor; // 12^(2/3) * 0.3 / 20 = 0.07862
624  } else {
625  ftm = 0.3 * pow(fcm - 8., 2. / 3.) * 1.e6 / MPSMaterial :: stiffnessFactor; //5.1-3a
626  }
627 
628  /*
629  if ( fcm >= 20. ) {
630  ftm = 0.3 * pow(fcm - 8., 2. / 3.) * 1.e6 / MPSMaterial :: stiffnessFactor; //5.1-3a
631  } else if ( fcm < 8. ) {
632  // upper formula does not hold for concretes with fcm < 8 MPa
633  ftm = 0.3 * pow(fcm, 2. / 3.) * 1.e6 / MPSMaterial :: stiffnessFactor;
634  } else {
635  // smooth transition
636  ftm = 0.3 * pow(fcm - ( 8. * ( fcm - 8. ) / ( 20. - 8. ) ), 2. / 3.) * 1.e6 / MPSMaterial :: stiffnessFactor;
637  }*/
638 
639  return ftm;
640 }
641 
642 
643 void
644 MPSDamMaterial :: computeDamage(double &omega, double kappa, GaussPoint *gp)
645 {
646  if ( this->softType == ST_Disable_Damage ) { //dummy material with no damage
647  omega = 0.;
648  } else {
649  computeDamageForCohesiveCrack(omega, kappa, gp);
650  }
651 }
652 
653 void
655 {
656  MPSDamMaterialStatus *status = NULL;
657 
658  omega = 0.0;
659  double e0 = this->givee0(gp);
660  double ef = 0.;
661 
662  if ( kappa > e0 ) {
663  double gf = this->givegf(gp);
664 
665  double Le;
666  double wf = 0.;
667  if ( softType == ST_Exponential_Cohesive_Crack ) { // exponential softening
668  wf = gf / this->E / e0; // wf is the crack opening
669  } else if ( softType == ST_Linear_Cohesive_Crack ) { // linear softening law
670  wf = 2. * gf / this->E / e0; // wf is the crack opening
671  } else {
672  OOFEM_ERROR("Gf unsupported for softening type softType = %d", softType);
673  }
674 
675  // MPSDamMaterialStatus *status = static_cast< MPSDamMaterialStatus * >( this->giveStatus(gp) );
676  status = static_cast< MPSDamMaterialStatus * >( this->giveStatus(gp) );
677 
678  Le = status->giveCharLength();
679  ef = wf / Le; //ef is the fracturing strain
680  if ( ef < e0 ) { //check that no snapback occurs
681  double minGf = 0.;
682  OOFEM_WARNING("ef %e < e0 %e, this leads to material snapback in element %d, characteristic length %f", ef, e0, gp->giveElement()->giveNumber(), Le);
683 
684  if ( softType == ST_Exponential_Cohesive_Crack ) { //exponential softening
685  minGf = this->E * e0 * e0 * Le;
686  } else if ( softType == ST_Linear_Cohesive_Crack ) { //linear softening law
687  minGf = this->E * e0 * e0 * Le / 2.;
688  } else {
689  OOFEM_WARNING("Gf unsupported for softening type softType = %d", softType);
690  }
691 
692  OOFEM_WARNING("Material number %d, decrease e0, or increase Gf from %f to Gf=%f", this->giveNumber(), gf, minGf);
693  if ( checkSnapBack ) {
694  OOFEM_ERROR("");
695  }
696  }
697 
698  if ( this->softType == ST_Linear_Cohesive_Crack ) {
699  if ( kappa < ef ) {
700  omega = ( ef / kappa ) * ( kappa - e0 ) / ( ef - e0 );
701  } else {
702  omega = 1.0; //maximum omega (maxOmega) is adjusted just for stiffness matrix in isodamagemodel.C
703  }
704  } else if ( this->softType == ST_Exponential_Cohesive_Crack ) {
705  // exponential cohesive crack - iteration needed
706  double R, Lhs, help;
707  int nite = 0;
708  // iteration to achieve objectivity
709  // we are looking for a state in which the elastic stress is equal to
710  // the stress from crack-opening relation
711  // ef has now the meaning of strain
712  do {
713  nite++;
714  help = omega * kappa / ef;
715  R = ( 1. - omega ) * kappa - e0 *exp(-help); //residuum
716  Lhs = kappa - e0 *exp(-help) * kappa / ef; //- dR / (d omega)
717  omega += R / Lhs;
718  if ( nite > 40 ) {
719  OOFEM_ERROR("algorithm not converging");
720  }
721  } while ( fabs(R) >= e0 * MPSDAMMAT_ITERATION_LIMIT );
722  } else {
723  OOFEM_ERROR("Unknown softening type for cohesive crack model.");
724  }
725 
726  if ( omega > 1.0 ) {
727  OOFEM_ERROR("damage parameter is %f, which is greater than 1, snap-back problems", omega);
728  }
729 
730  if ( omega < 0.0 ) {
731  OOFEM_WARNING("damage parameter is %f, which is smaller than 0, snap-back problems", omega);
732  omega = 1.;
733  if ( checkSnapBack ) {
734  OOFEM_ERROR("");
735  }
736 
737  }
738  }
739 
740 
741 #ifdef supplementary_info
742  double residualStrength = 0.;
743 
744  if ( omega == 0. ) {
745  residualStrength = E * e0; // undamaged material
746  status = static_cast< MPSDamMaterialStatus * >( this->giveStatus(gp) );
747  } else {
748  if ( this->softType == ST_Linear_Cohesive_Crack ) {
749  residualStrength = E * e0 * ( ef - kappa ) / ( ef - e0 );
750  } else if ( this->softType == ST_Exponential_Cohesive_Crack ) {
751  residualStrength = E * e0 * exp(-1. * ( kappa - e0 ) / ef);
752  } else {
753  OOFEM_ERROR("Unknown softening type for cohesive crack model.");
754  }
755  }
756 
757  if ( status ) {
758  status->setResidualTensileStrength(residualStrength);
759  }
760 
761 #endif
762 }
763 
764 void
765 MPSDamMaterial :: initDamaged(double kappa, FloatArray &principalDirection, GaussPoint *gp, TimeStep *tStep)
766 {
767  double le = 0.;
768 
769  MPSDamMaterialStatus *status = static_cast< MPSDamMaterialStatus * >( this->giveStatus(gp) );
770 
771  if ( this->timeDepFracturing ) {
772  this->initDamagedFib(gp, tStep);
773  }
774 
775  double e0 = this->givee0(gp);
776  double gf = this->givegf(gp);
777  double wf = 0.;
778 
779  if ( softType == ST_Disable_Damage ) {
780  return;
781  } else if ( softType == ST_Exponential_Cohesive_Crack ) { // exponential softening
782  wf = gf / E / e0; // wf is the crack opening
783  } else if ( softType == ST_Linear_Cohesive_Crack ) { // linear softening law
784  wf = 2. * gf / E / e0; // wf is the crack opening
785  } else {
786  OOFEM_ERROR("Gf unsupported for softening type softType = %d", softType);
787  }
788 
789  if ( ( kappa > e0 ) && ( status->giveDamage() == 0. ) ) {
790  status->setCrackVector(principalDirection);
791 
792  le = gp->giveElement()->giveCharacteristicSize(gp, principalDirection, ecsMethod);
793  status->setCharLength(le);
794 
795  if ( gf != 0. && e0 >= ( wf / le ) ) { // case for a given fracture energy
796  OOFEM_WARNING("Fracturing strain %e is lower than the elastic strain e0=%e, possible snap-back. Element number %d, wf %e, le %e. Increase fracturing strain or decrease element size by at least %f", wf / le, e0, gp->giveElement()->giveLabel(), wf, le, e0/(wf/le) );
797  if ( checkSnapBack ) {
798  OOFEM_ERROR("");
799  }
800  }
801  }
802 }
803 
804 
805 
808 /*
809  * creates a new material status corresponding to this class
810  */
811 {
812  return new MPSDamMaterialStatus(1, this->giveDomain(), gp, nUnits);
813 }
814 
815 
816 void
818  MatResponseMode mode,
819  GaussPoint *gp,
820  TimeStep *tStep)
821 {
822  RheoChainMaterial :: give3dMaterialStiffnessMatrix(answer, ElasticStiffness, gp, tStep);
823 
824  if ( mode == ElasticStiffness || ( mode == SecantStiffness && !this->isotropic ) ) {
825  return;
826  }
827 
828  double tempDamage;
829  MPSDamMaterialStatus *status = ( MPSDamMaterialStatus * ) this->giveStatus(gp);
830  tempDamage = min(status->giveTempDamage(), this->maxOmega);
831  answer.times(1.0 - tempDamage);
832 }
833 
834 
835 void
837  MatResponseMode mode,
838  GaussPoint *gp,
839  TimeStep *tStep)
840 {
841  RheoChainMaterial :: givePlaneStressStiffMtrx(answer, ElasticStiffness, gp, tStep);
842 
843  if ( mode == ElasticStiffness || ( mode == SecantStiffness && !this->isotropic ) ) {
844  return;
845  }
846 
847  double tempDamage;
848  MPSDamMaterialStatus *status = ( MPSDamMaterialStatus * ) this->giveStatus(gp);
849  tempDamage = min(status->giveTempDamage(), this->maxOmega);
850  answer.times(1.0 - tempDamage);
851 }
852 
853 void
855  MatResponseMode mode,
856  GaussPoint *gp,
857  TimeStep *tStep)
858 {
859  RheoChainMaterial :: givePlaneStrainStiffMtrx(answer, ElasticStiffness, gp, tStep);
860 
861  if ( mode == ElasticStiffness || ( mode == SecantStiffness && !this->isotropic ) ) {
862  return;
863  }
864 
865  double tempDamage;
866  MPSDamMaterialStatus *status = ( MPSDamMaterialStatus * ) this->giveStatus(gp);
867  tempDamage = min(status->giveTempDamage(), this->maxOmega);
868  answer.times(1.0 - tempDamage);
869 }
870 
871 
872 void
874  MatResponseMode mode,
875  GaussPoint *gp,
876  TimeStep *tStep)
877 {
878  RheoChainMaterial :: give1dStressStiffMtrx(answer, ElasticStiffness, gp, tStep);
879 
880  if ( mode == ElasticStiffness || ( mode == SecantStiffness && !this->isotropic ) ) {
881  return;
882  }
883 
884  double tempDamage;
885  MPSDamMaterialStatus *status = ( MPSDamMaterialStatus * ) this->giveStatus(gp);
886  tempDamage = min(status->giveTempDamage(), this->maxOmega);
887  answer.times(1.0 - tempDamage);
888 }
889 
890 
891 int
893 {
894  MPSDamMaterialStatus *status = static_cast< MPSDamMaterialStatus * >( this->giveStatus(gp) );
895  if ( type == IST_DamageScalar ) {
896  answer.resize(1);
897  answer.zero();
898  answer.at(1) = status->giveDamage();
899  return 1;
900  } else if ( type == IST_DamageTensor ) {
901  answer.resize(6);
902  answer.zero();
903  answer.at(1) = answer.at(2) = answer.at(3) = status->giveDamage();
904  return 1;
905  } else if ( type == IST_PrincipalDamageTensor ) {
906  answer.resize(3);
907  answer.zero();
908  answer.at(1) = status->giveDamage();
909  return 1;
910  } else if ( type == IST_DamageTensorTemp ) {
911  answer.resize(6);
912  answer.zero();
913  answer.at(1) = answer.at(2) = answer.at(3) = status->giveTempDamage();
914  return 1;
915  } else if ( type == IST_PrincipalDamageTempTensor ) {
916  answer.resize(3);
917  answer.zero();
918  answer.at(1) = status->giveTempDamage();
919  return 1;
920  } else if ( type == IST_CharacteristicLength ) {
921  answer.resize(1);
922  answer.zero();
923  answer.at(1) = status->giveCharLength();
924  return 1;
925  } else if ( type == IST_CrackVector ) {
926  answer.resize(3);
927  answer.zero();
928  status->giveCrackVector(answer);
929  return 1;
930  } else if ( type == IST_CrackWidth ) {
931  answer.resize(1);
932  answer.zero();
933  answer.at(1) = status->giveCrackWidth();
934  return 1;
935  } else if ( type == IST_ResidualTensileStrength ) {
936  answer.resize(1);
937  answer.zero();
938  answer.at(1) = status->giveResidualTensileStrength();
939  return 1;
940  } else if ( type == IST_TensileStrength ) {
941  MPSDamMaterialStatus *status = static_cast< MPSDamMaterialStatus * >( this->giveStatus(gp) );
942  double tequiv = status->giveEquivalentTime();
943  answer.resize(1);
944  answer.zero();
945  if (tequiv >= this->castingTime) {
946  answer.at(1) = this->computeTensileStrength(tequiv);
947  }
948  return 1;
949  } else if ( type == IST_CrackIndex ) {
950  //ratio of real principal stress / strength. 1 if damage already occured.
951  FloatArray principalStress;
952  MPSDamMaterialStatus *status = static_cast< MPSDamMaterialStatus * >( this->giveStatus(gp) );
953  answer.resize(1);
954  answer.zero();
955  if ( status->giveDamage()>0. ){
956  answer.at(1)=1.;
957  return 1;
958  }
959  //FloatArray effectiveStress = status->giveTempViscoelasticStressVector();
960  //StructuralMaterial :: computePrincipalValues(principalStress, effectiveStress, principal_stress);
961  StructuralMaterial :: giveIPValue(principalStress, gp, IST_PrincipalStressTensor, tStep);
962  double tequiv = status->giveEquivalentTime();
963  if (tequiv >= this->castingTime) {
964  double ft = this->computeTensileStrength(tequiv);
965  if (ft > 1.e-20 && principalStress.at(1)>1.e-20){
966  answer.at(1) = principalStress.at(1)/ft;
967  }
968  }
969  return 1;
970  } else {
971  return MPSMaterial :: giveIPValue(answer, gp, type, tStep);
972  }
973 
974  return 1; // to make the compiler happy
975 }
976 } // end namespace oofem
double computeEquivalentTime(GaussPoint *gp, TimeStep *tStep, int option)
Computes equivalent time at given time step and GP.
Definition: mps.C:1542
double giveDamage()
Returns the last equilibrated damage level.
Definition: mpsdammat.h:119
static int giveSizeOfVoigtSymVector(MaterialMode mmode)
Returns the size of symmetric part of a reduced stress/strain vector according to given mode...
#define _IFT_MPSDamMaterial_isotropic
Definition: mpsdammat.h:53
virtual contextIOResultType restoreContext(DataStream &stream, ContextMode mode, void *obj=NULL)
Restores the receiver state previously written in stream.
Definition: mps.C:164
#define _IFT_MPSDamMaterial_checkSnapBack
Definition: mpsdammat.h:58
InternalStateType
Type representing the physical meaning of element or constitutive model internal variable.
static void computePrincipalValues(FloatArray &answer, const FloatArray &s, stressStrainPrincMode mode)
Common functions for convenience.
MaterialMode giveMaterialMode()
Returns corresponding material mode of receiver.
Definition: gausspoint.h:191
void letTempViscoelasticStressVectorBe(FloatArray v)
Assigns tempStressVector to given vector v.
Definition: mpsdammat.h:109
virtual void updateYourself(TimeStep *tStep)
Update equilibrium history variables according to temp-variables.
Definition: mpsdammat.C:93
void letTempStrainVectorBe(const FloatArray &v)
Assigns tempStrainVector to given vector v.
Definition: structuralms.h:137
static void givePlaneStressVectorTranformationMtrx(FloatMatrix &answer, const FloatMatrix &base, bool transpose=false)
Computes 2d stress vector transformation matrix from standard vector transformation matrix...
void setCharLength(double length)
Sets characteristic length to given value.
Definition: mpsdammat.h:128
virtual MaterialStatus * giveStatus(GaussPoint *gp) const
Returns material status of receiver in given integration point.
Definition: material.C:244
virtual bool isActivated(TimeStep *tStep)
Definition: rheoChM.h:292
double ft
Equivalent strain at stress peak (or a similar parameter).
Definition: mpsdammat.h:201
GaussPoint * gp
Associated integration point.
Class and object Domain.
Definition: domain.h:115
MPSDamMaterialStatus(int n, Domain *d, GaussPoint *g, int nunits)
Definition: mpsdammat.C:54
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: rheoChM.C:459
For computing principal strains from engineering strains.
ElementCharSizeMethod ecsMethod
Method used for evaluation of characteristic element size.
Definition: mpsdammat.h:223
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: mpsdammat.C:817
void initDamaged(double kappa, FloatArray &totalStrainVector, GaussPoint *gp, TimeStep *tStep)
Abstract service allowing to perform some initialization, when damage first appear.
Definition: mpsdammat.C:765
double giveKappa()
Returns the last equilibrated scalar measure of the largest strain level.
Definition: mpsdammat.h:113
The purpose of DataStream abstract class is to allow to store/restore context to different streams...
Definition: datastream.h:54
void setTempKappa(double newKappa)
Sets the temp scalar measure of the largest strain level to given value.
Definition: mpsdammat.h:117
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.
Specialization of a floating point array for representing a stress state.
Definition: stressvector.h:48
virtual void give1dStressStiffMtrx(FloatMatrix &answer, MatResponseMode mmode, GaussPoint *gp, TimeStep *tStep)
Method for computing 1d stiffness matrix of receiver.
Definition: rheoChM.C:521
double & at(int i)
Coefficient access function.
Definition: floatarray.h:131
#define _IFT_MPSDamMaterial_gf
Definition: mpsdammat.h:60
void computeDamageForCohesiveCrack(double &omega, double kappa, GaussPoint *gp)
computes the value of damage parameter omega, based on a given value of equivalent strain...
Definition: mpsdammat.C:654
void computePrincipalValDir(FloatArray &answer, FloatMatrix &dir) const
Computes principal values and directions of receiver vector.
Definition: stressvector.C:182
virtual IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description stored in input record.
Definition: mps.C:192
virtual void printOutputAt(FILE *file, TimeStep *tStep)
Print receiver&#39;s output to given stream.
Definition: kelvinChSolM.C:275
virtual bool hasField(InputFieldType id)=0
Returns true if record contains field identified by idString keyword.
double giveEquivalentTime()
Returns equivalent time.
Definition: mps.h:160
double giveTempThermalStrain(void)
Definition: rheoChM.h:116
virtual MaterialStatus * CreateStatus(GaussPoint *gp) const
Creates new copy of associated status and inserts it into given integration point.
Definition: mpsdammat.C:807
double giveTempDamage()
Returns the temp. damage level.
Definition: mpsdammat.h:121
Element * giveElement()
Returns corresponding element to receiver.
Definition: gausspoint.h:188
General IO error.
void setResidualTensileStrength(double src)
Definition: mpsdammat.h:146
double kappa
Scalar measure of the largest strain level ever reached in material.
Definition: mpsdammat.h:80
virtual void printOutputAt(FILE *file, TimeStep *tStep)
Print receiver&#39;s output to given stream.
Definition: mpsdammat.C:116
MaterialMode
Type representing material mode of integration point.
Definition: materialmode.h:89
#define _IFT_MPSDamMaterial_fib_s
Definition: mpsdammat.h:51
MatResponseMode
Describes the character of characteristic material matrix.
virtual int read(int *data, int count)=0
Reads count integer values into array pointed by data.
#define _IFT_MPSDamMaterial_timedepfracturing
Definition: mpsdammat.h:50
virtual int giveIPValue(FloatArray &answer, GaussPoint *gp, InternalStateType type, TimeStep *tStep)
Returns the integration point corresponding value in Reduced form.
Definition: mpsdammat.C:892
MPSDamMaterial(int n, Domain *d)
Definition: mpsdammat.C:222
#define THROW_CIOERR(e)
Definition: contextioerr.h:61
void rotatedWith(FloatMatrix &r, char mode)
Returns the receiver a rotated according the change-of-base matrix r.
Definition: floatarray.C:799
double giveCrackWidth(void)
Definition: mpsdammat.h:145
double E
dummy Young&#39;s modulus
Definition: mpsdammat.h:192
double lambda0
constant equal to one day in time units of analysis (eg. 86400 if the analysis runs in seconds) ...
Definition: mps.h:214
double ft28
28-day value of tensile strength. Used only with "timedepfracturing"
Definition: mpsdammat.h:226
void setCrackWidth(double src)
Definition: mpsdammat.h:144
FloatArray effectiveStressVector
Equilibrated stress vector in reduced form.
Definition: mpsdammat.h:76
virtual int write(const int *data, int count)=0
Writes count integer values from array pointed by data.
void giveCrackVector(FloatArray &answer)
Returns crack vector stored in receiver. This is useful for plotting cracks as a vector field (paravi...
Definition: mpsdammat.C:108
void setCrackVector(FloatArray cv)
Sets crack vector to given value. This is useful for plotting cracks as a vector field (paraview etc...
Definition: mpsdammat.h:132
virtual int hasMaterialModeCapability(MaterialMode mode)
Tests if material supports material mode.
Definition: mpsdammat.C:239
double const_gf
Determines the softening -> corresponds to the initial fracture energy.
Definition: mpsdammat.h:209
int checkSnapBack
Check possible snap back flag.
Definition: mpsdammat.h:217
This class implements the extended B3 model for concrete creep and shrinkage based on the microprestr...
Definition: mps.h:204
void setTempDamage(double newDamage)
Sets the temp damage level to given value.
Definition: mpsdammat.h:123
virtual contextIOResultType saveContext(DataStream &stream, ContextMode mode, void *obj=NULL)
Stores receiver state to output stream.
Definition: mpsdammat.C:132
double stiffnessFactor
scaling factor 1. for Pa, 1.e6 for MPa - only for empirical formulas - q1-q4 and ft and gf ...
Definition: mps.h:251
virtual contextIOResultType restoreContext(DataStream &stream, ContextMode mode, void *obj=NULL)
Restores the receiver state previously written in stream.
Definition: mpsdammat.C:175
double var_gf
hydration-degree dependent fracture energy
Definition: mpsdammat.h:96
double tempDamage
Non-equilibrated damage level of material.
Definition: mpsdammat.h:86
#define MPSDAMMAT_ITERATION_LIMIT
Definition: mpsdammat.h:68
#define OOFEM_ERROR(...)
Definition: error.h:61
#define _IFT_MPSDamMaterial_ft
Definition: mpsdammat.h:59
double gf28
28-day value of fracture energy. Used only with "timedepfracturing"
Definition: mpsdammat.h:228
#define _IFT_MPSDamMaterial_ft28
Definition: mpsdammat.h:62
void times(double f)
Multiplies receiver by factor f.
Definition: floatmatrix.C:1594
double giveTempAutogenousShrinkageStrain(void)
Definition: mps.h:182
FloatArray crackVector
Crack orientation normalized to damage magnitude. This is useful for plotting cracks as a vector fiel...
Definition: mpsdammat.h:91
contextIOResultType restoreYourself(DataStream &stream)
Definition: floatarray.C:895
double giveTempDryingShrinkageStrain(void)
Definition: mps.h:178
virtual IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description stored in input record.
Definition: mpsdammat.C:250
#define _IFT_MPSDamMaterial_damageLaw
Definition: mpsdammat.h:57
double maxOmega
Maximum limit on omega. The purpose is elimination of a too compliant material which may cause conver...
Definition: mpsdammat.h:195
void resizeWithValues(int s, int allocChunk=0)
Checks size of receiver towards requested bounds.
Definition: floatarray.C:615
double at(int i, int j) const
Coefficient access function.
Definition: floatmatrix.h:176
#define _IFT_MPSDamMaterial_gf28
Definition: mpsdammat.h:63
double givee0(GaussPoint *gp)
Definition: mpsdammat.C:533
virtual void givePlaneStressStiffMtrx(FloatMatrix &answer, MatResponseMode mmode, GaussPoint *gp, TimeStep *tStep)
Method for computing plane stress stiffness matrix of receiver.
Definition: rheoChM.C:486
double giveResidualTensileStrength(void)
Definition: mpsdammat.h:147
Abstract base class representing a material status information.
Definition: matstatus.h:84
double giveCharLength()
Returns characteristic length stored in receiver.
Definition: mpsdammat.h:126
Class representing vector of real numbers.
Definition: floatarray.h:82
virtual double giveCharacteristicSize(GaussPoint *gp, FloatArray &normalToCrackPlane, ElementCharSizeMethod method)
Returns characteristic element size for a given integration point and given direction.
Definition: element.h:901
virtual void givePlaneStrainStiffMtrx(FloatMatrix &answer, MatResponseMode mmode, GaussPoint *gp, TimeStep *tStep)
Method for computing plane strain stiffness matrix of receiver.
Definition: rheoChM.C:503
Implementation of matrix containing floating point numbers.
Definition: floatmatrix.h:94
double damage
Damage level of material.
Definition: mpsdammat.h:84
IRResultType
Type defining the return values of InputRecord reading operations.
Definition: irresulttype.h:47
This class implements associated Material Status to MPSMaterial, which corresponds to a model for hum...
Definition: mps.h:97
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
virtual double computeFractureEnergy(double equivalentTime)
Definition: mpsdammat.C:556
virtual double computeTensileStrength(double equivalentTime)
Definition: mpsdammat.C:603
Class representing the general Input Record.
Definition: inputrecord.h:101
void zero()
Zeroes all coefficients of receiver.
Definition: floatarray.C:658
virtual double computeCreepFunction(double t, double t_prime, GaussPoint *gp, TimeStep *tStep)
Evaluation of the basic creep compliance function - can be used to compute elastic modulus in derived...
Definition: mps.C:582
FloatArray tempEffectiveStressVector
Temporary stress vector in reduced form (increments are used mainly in nonlinear analysis) ...
Definition: mpsdammat.h:78
virtual void updateYourself(TimeStep *tStep)
Update equilibrium history variables according to temp-variables.
Definition: mps.C:81
void convertFromFullForm(const FloatArray &reducedform, MaterialMode mode)
Assign to receiver the reduced form of given vector.
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
void initDamagedFib(GaussPoint *gp, TimeStep *tStep)
Definition: mpsdammat.C:515
double tempKappa
Non-equilibrated scalar measure of the largest strain level.
Definition: mpsdammat.h:82
Domain * giveDomain() const
Definition: femcmpnn.h:100
virtual void initTempStatus()
Initializes the temporary internal variables, describing the current state according to previously re...
Definition: mpsdammat.C:78
int min(int i, int j)
Returns smaller value from two given decimals.
Definition: mathfem.h:59
#define _IFT_MPSMaterial_fc
Definition: mps.h:49
REGISTER_Material(DummyMaterial)
double castingTime
Casting time.
Definition: material.h:113
#define IR_GIVE_OPTIONAL_FIELD(__ir, __value, __id)
Macro facilitating the use of input record reading methods.
Definition: inputrecord.h:78
double givegf(GaussPoint *gp)
Definition: mpsdammat.C:544
void computeDamage(double &omega, double kappa, GaussPoint *gp)
Computes the value of damage omega, based on given value of equivalent strain.
Definition: mpsdammat.C:644
int giveSize() const
Returns the size of receiver.
Definition: floatarray.h:218
int giveLabel() const
Definition: element.h:1055
the oofem namespace is to define a context or scope in which all oofem names are defined.
#define IR_GIVE_FIELD(__ir, __value, __id)
Macro facilitating the use of input record reading methods.
Definition: inputrecord.h:69
virtual void initTempStatus()
Initializes the temporary internal variables, describing the current state according to previously re...
Definition: mps.C:114
int giveNumber() const
Definition: femcmpnn.h:107
int nUnits
Number of (Maxwell or Kelvin) units in the rheologic chain.
Definition: rheoChM.h:136
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: mps.C:430
double charLength
Characteristic length.
Definition: mpsdammat.h:89
Class representing integration point in finite element program.
Definition: gausspoint.h:93
#define OOFEM_WARNING(...)
Definition: error.h:62
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: mpsdammat.C:326
Class representing solution step.
Definition: timestep.h:80
virtual void givePlaneStressStiffMtrx(FloatMatrix &answer, MatResponseMode mode, GaussPoint *gp, TimeStep *tStep)
Method for computing plane stress stiffness matrix of receiver.
Definition: mpsdammat.C:836
void add(const FloatArray &src)
Adds array src to receiver.
Definition: floatarray.C:156
void sete0(double e0)
Definition: mpsdammat.h:134
virtual void givePlaneStrainStiffMtrx(FloatMatrix &answer, MatResponseMode mode, GaussPoint *gp, TimeStep *tStep)
Method for computing plane strain stiffness matrix of receiver.
Definition: mpsdammat.C:854
double var_e0
hydration-degree dependent equivalent strain at stress peak
Definition: mpsdammat.h:94
void setgf(double gf)
Definition: mpsdammat.h:135
virtual int giveIPValue(FloatArray &answer, GaussPoint *gp, InternalStateType type, TimeStep *tStep)
Returns the integration point corresponding value in Reduced form.
Definition: mps.C:1577
virtual void give1dStressStiffMtrx(FloatMatrix &answer, MatResponseMode mode, GaussPoint *gp, TimeStep *tStep)
Method for computing 1d stiffness matrix of receiver.
Definition: mpsdammat.C:873
virtual contextIOResultType saveContext(DataStream &stream, ContextMode mode, void *obj=NULL)
Stores receiver state to output stream.
Definition: mps.C:140
void resize(int s)
Resizes receiver towards requested size.
Definition: floatarray.C:631
SofteningType softType
Parameter specifying the type of softening (damage law).
Definition: mpsdammat.h:220

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:30 for OOFEM by doxygen 1.8.11 written by Dimitri van Heesch, © 1997-2011