OOFEM  2.4
OOFEM.org - Object Oriented Finite Element Solver
concretedpm2.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 "concretedpm2.h"
36 #include "floatarray.h"
37 #include "floatmatrix.h"
38 #include "../sm/Materials/structuralms.h"
39 #include "gausspoint.h"
40 #include "intarray.h"
41 #include "datastream.h"
42 #include "contextioerr.h"
43 #include "timestep.h"
44 #include "../sm/Materials/structuralmaterial.h"
46 #include "../sm/CrossSections/structuralcrosssection.h"
47 #include "mathfem.h"
48 #include "classfactory.h"
49 #include <limits>
50 
51 namespace oofem {
52 REGISTER_Material(ConcreteDPM2);
53 
55  StructuralMaterialStatus(n, d, gp),
56  plasticStrain(6),
57  tempPlasticStrain(6)
58 {
63 
64  kappaP = tempKappaP = 0.;
65 
69 
72 
77 
78  alpha = tempAlpha = 0.;
79 
82 
83  deltaLambda = 0.;
85  rateFactor = 1.;
87 
88 #ifdef keep_track_of_dissipated_energy
89  stressWork = tempStressWork = 0.0;
90  dissWork = tempDissWork = 0.0;
91 #endif
92 }
93 
95 { }
96 
97 void
99 {
100  // Call the function of the parent class to initialize the variables defined there.
102 
103  //Plasticity part
105  tempKappaP = kappaP;
106 
107  tempAlpha = alpha;
108 
109  //Damage part
113 
116 
121 
124 
128 #ifdef keep_track_of_dissipated_energy
131 #endif
132 }
133 
134 void
136 {
137  // Call corresponding function of the parent class to update
138  // variables defined there.
140 
141  // update variables defined in ConcreteDPM2Status
142 
143  alpha = tempAlpha;
144 
145  //Plasticity part
147  kappaP = tempKappaP;
148 
149  //Damage part
153 
156 
159 
162 
165 
167 
169 
171 #ifdef keep_track_of_dissipated_energy
174 #endif
175 }
176 
177 void
179 {
180  // Call corresponding function of the parent class to print
182 
183  fprintf(file, "\tstatus { ");
184 
185  // print status flag
186  switch ( state_flag ) {
188  fprintf(file, "Elastic, ");
189  break;
191  fprintf(file, "Unloading, ");
192  break;
194  fprintf(file, "Plastic, ");
195  break;
197  fprintf(file, "Damage, ");
198  break;
200  fprintf(file, "PlasticDamage, ");
201  break;
202  }
203 
204  // print plastic strain vector and inelastic strain vector
205  FloatArray plasticStrainVector = this->givePlasticStrain();
206  FloatArray inelasticStrainVector = strainVector;
207  inelasticStrainVector.subtract(plasticStrainVector);
208  inelasticStrainVector.times(damageTension);
209  inelasticStrainVector.add(plasticStrainVector);
210  inelasticStrainVector.times(le);
211 
212  fprintf(file, " plastic ");
213  for ( auto &val : plasticStrainVector ) {
214  fprintf( file, " %.10e", val );
215  }
216 
217  fprintf(file, " inelastic ");
218  for ( auto &val : inelasticStrainVector ) {
219  fprintf( file, " %.10e", val );
220  }
221 
222  fprintf(file, " equivStrain %.10e,", equivStrain);
223 
224  fprintf(file, " kappaDTension %.10e,", kappaDTension);
225 
226  fprintf(file, " kappaDCompression %.10e,", kappaDCompression);
227 
228  fprintf(file, " kappaP %.10e,", kappaP);
229 
230  fprintf(file, " kappaDTensionOne %.10e,", kappaDTensionOne);
231 
232  fprintf(file, " kappaDCompressionOne %.10e,", kappaDCompressionOne);
233 
234  fprintf(file, " kappaDTensionTwo %.10e,", kappaDTensionTwo);
235 
236  fprintf(file, " kappaDCompressionTwo %.10e,", kappaDCompressionTwo);
237 
238  fprintf(file, " damageTension %.10e,", damageTension);
239 
240  fprintf(file, " damageCompression %.10e,", damageCompression);
241 
242  fprintf(file, " alpha %.10e,", this->alpha);
243 
244 #ifdef keep_track_of_dissipated_energy
245  fprintf(file, " dissW %g, freeE %g, stressW %g ", this->dissWork, ( this->stressWork ) - ( this->dissWork ), this->stressWork);
246 #endif
247  fprintf(file, "}\n");
248 }
249 
252 {
253  contextIOResultType iores;
254 
255  // save parent class status
256  if ( ( iores = StructuralMaterialStatus :: saveContext(stream, mode, obj) ) != CIO_OK ) {
257  THROW_CIOERR(iores);
258  }
259 
260  if ( ( iores = plasticStrain.storeYourself(stream) ) != CIO_OK ) {
261  THROW_CIOERR(iores);
262  }
263 
264  if ( !stream.write(kappaP) ) {
266  }
267 
268  if ( !stream.write(deltaLambda) ) {
270  }
271 
272  if ( !stream.write(dFDKappa) ) {
274  }
275 
276  if ( !stream.write(le) ) {
278  }
279 
280  if ( !stream.write(alpha) ) {
282  }
283 
284  if ( !stream.write(equivStrainTension) ) {
286  }
287 
288  if ( !stream.write(equivStrainCompression) ) {
290  }
291 
292  if ( !stream.write(kappaDTension) ) {
294  }
295 
296  if ( !stream.write(kappaDCompression) ) {
298  }
299 
300  if ( !stream.write(kappaDCompressionOne) ) {
302  }
303 
304  if ( !stream.write(kappaDTensionTwo) ) {
306  }
307 
308  if ( !stream.write(kappaDCompressionTwo) ) {
310  }
311 
312 
313  if ( !stream.write(damageTension) ) {
315  }
316 
317  if ( !stream.write(damageCompression) ) {
319  }
320 
321  if ( !stream.write(deltaEquivStrain) ) {
323  }
324 
325  if ( !stream.write(rateFactor) ) {
327  }
328 
329  if ( !stream.write(rateStrain) ) {
331  }
332 
333 
334  if ( !stream.write(state_flag) ) {
336  }
337 
338 #ifdef keep_track_of_dissipated_energy
339  if ( !stream.write(stressWork) ) {
341  }
342 
343  if ( !stream.write(dissWork) ) {
345  }
346 
347 #endif
348  return CIO_OK;
349 }
350 
351 
354 {
355  contextIOResultType iores;
356 
357  // read parent class status
358  if ( ( iores = StructuralMaterialStatus :: restoreContext(stream, mode, obj) ) != CIO_OK ) {
359  THROW_CIOERR(iores);
360  }
361 
362  // read raw data
363  if ( ( iores = plasticStrain.restoreYourself(stream) ) != CIO_OK ) {
364  THROW_CIOERR(iores);
365  }
366 
367  if ( !stream.read(kappaP) ) {
369  }
370 
371  if ( !stream.read(deltaLambda) ) {
373  }
374 
375  if ( !stream.read(dFDKappa) ) {
377  }
378 
379  if ( !stream.read(le) ) {
381  }
382 
383  if ( !stream.read(alpha) ) {
385  }
386 
387  if ( !stream.read(equivStrainTension) ) {
389  }
390 
391  if ( !stream.read(equivStrainCompression) ) {
393  }
394 
395  if ( !stream.read(kappaDTension) ) {
397  }
398 
399  if ( !stream.read(kappaDCompression) ) {
401  }
402 
403  if ( !stream.read(kappaDCompressionOne) ) {
405  }
406 
407  if ( !stream.read(kappaDTensionTwo) ) {
409  }
410 
411  if ( !stream.read(kappaDCompressionTwo) ) {
413  }
414 
415 
416  if ( !stream.read(damageTension) ) {
418  }
419 
420  if ( !stream.read(damageCompression) ) {
422  }
423 
424  if ( !stream.read(deltaEquivStrain) ) {
426  }
427 
428  if ( !stream.read(rateFactor) ) {
430  }
431 
432  if ( !stream.read(rateStrain) ) {
434  }
435 
436 
437  if ( !stream.read(state_flag) ) {
439  }
440 
441 
442 #ifdef keep_track_of_dissipated_energy
443  if ( !stream.read(stressWork) ) {
445  }
446 
447  if ( !stream.read(dissWork) ) {
449  }
450 
451 #endif
452  return CIO_OK;
453 }
454 
455 #ifdef keep_track_of_dissipated_energy
456 void
458 {
459  FloatArray tempTotalstrain = tempStrainVector;
460  FloatArray totalstrain = strainVector;
461 
462  //Calculate increase or decrease of total strain tensor during iteration/step
463  FloatArray deltaTotalStrain = tempTotalstrain;
464  deltaTotalStrain.subtract(totalstrain);
465 
466  //Ask for stress tensor at step
467  FloatArray stress = tempStressVector;
468  int n = stress.giveSize();
469  //Calculate increase/decrease in total work
470  double dSW = ( tempStressVector.dotProduct(deltaTotalStrain, n) + stressVector.dotProduct(deltaTotalStrain, n) ) / 2.;
471 
472  double tempStressWork = this->giveTempStressWork() + dSW;
473 
474  //Calculate temporary elastic strain
475  FloatArray tempElasticStrain = tempTotalstrain;
476  tempElasticStrain.subtract(tempPlasticStrain);
477 
478  //Calculate elastically stored energy density
479  double We = tempStressVector.dotProduct(tempElasticStrain, n) / 2.;
480 
481  // dissipative work density
482  tempDissWork = tempStressWork - We;
483 
484  // to avoid extremely small negative dissipation due to round-off error
485  // (note: gf is the dissipation density at complete failure, per unit volume)
486  //Rough estimation of gf which is ok for this purpose
487 
488  if ( fabs(tempDissWork) < 1.e-12 * gf ) {
489  tempDissWork = 0.;
490  }
491 
492  this->setTempStressWork(tempStressWork);
494 }
495 #endif
496 
497 // ********************************
498 // *** CLASS DYNAMIC CONCRETE ***
499 // ********************************
500 
501 #define IDM_ITERATION_LIMIT 1.e-8
502 
504  StructuralMaterial(n, d)
505 {
506  yieldTol = 0.;
507  yieldTolDamage = 0.;
508  newtonIter = 0;
509 
511 }
512 
514 {
515  delete linearElasticMaterial;
516 }
517 
520 {
521  // Required by IR_GIVE_FIELD macro
522  IRResultType result;
523 
524  // call the corresponding service for the linear elastic material
526  if ( result != IRRT_OK ) {
527  return result;
528  }
529 
531  if ( result != IRRT_OK ) {
532  return result;
533  }
534 
535  //isotropic flag
536  isotropicFlag = 0;
538 
539  double value;
540  // elastic parameters
543  propertyDictionary.add('E', eM);
544  propertyDictionary.add('n', nu);
545 
547  propertyDictionary.add(tAlpha, value);
548 
549  gM = eM / ( 2. * ( 1. + nu ) );
550  kM = eM / ( 3. * ( 1. - 2. * nu ) );
551 
554 
555  this->e0 = this->ft / this->eM;
556 
557  // default parameters
558  this->ecc = 0.525;
560  yieldHardInitial = 0.3;
562  //Inclination at transition point
563  yieldHardPrimePeak = 0.5;
565 
566  if ( yieldHardPrimePeak < 0 ) {
567  yieldHardPrimePeak = 0.;
568  OOFEM_WARNING("kPrimePeak cannot be less than zero");
569  } else if ( yieldHardPrimePeak > ( 1. - yieldHardInitial ) ) {
571  OOFEM_WARNING("kPrimePeak cannot be greater than 1.-kinit");
572  }
573 
574  AHard = 8.e-2;
576  BHard = 3.e-3;
578  CHard = 2.;
580  DHard = 1.e-6;
582  dilationConst = 0.85;
584 
585  softeningType = 1; //0-Linear softening; 1-Bilinear softening; 2-Exponential
587 
588  if ( softeningType > 2 ) {
589  OOFEM_WARNING("softening type not implemented");
590  return IRRT_BAD_FORMAT;
591  }
592 
594 
595  if ( softeningType == 1 ) {
596  this->ftOne = 0.3 * this->ft;
598  this->wfOne = 0.15 * this->wf;
600  }
601 
602  this->efCompression = 100.e-6;
604 
605  this->ASoft = 15;
607 
608  helem = 0.;
610 
611 
612  //Compute m
613  m = 3. * ( pow(fc, 2.) - pow(ft, 2.) ) / ( fc * ft ) * ecc / ( ecc + 1. );
614 
615  //Compute default value of dilationConst
616  yieldTol = 1.e-6;
618 
619  yieldTolDamage = yieldTol * 10.;
620 
621  newtonIter = 100;
623 
624  //parameters for rate dependence
625  strainRateFlag = 0;
627 
628  timeFactor = 1.;
629  if ( strainRateFlag == 1 ) {
632  }
633 
634  return IRRT_OK;
635 }
636 
637 
638 void
640  GaussPoint *gp,
641  const FloatArray &strainVector,
642  TimeStep *tStep)
643 {
644  ConcreteDPM2Status *status = static_cast< ConcreteDPM2Status * >( this->giveStatus(gp) );
645 
646  // Initialize temp variables for this gauss point
647  this->initTempStatus(gp);
648 
649  //Calculate strain rate
650  //Time step
651  double deltaTime = 1.;
652  if ( strainRateFlag == 1 ) {
653  if ( tStep->giveTimeIncrement() == 0 ) { //Problem with the first step. For some reason the time increment is zero
654  deltaTime = this->timeFactor;
655  } else {
656  deltaTime = tStep->giveTimeIncrement() * this->timeFactor;
657  }
658  } else {
659  if ( tStep->giveTimeIncrement() == 0 ) {
660  deltaTime = 1.;
661  } else {
662  deltaTime = tStep->giveTimeIncrement();
663  }
664  }
665 
666  FloatMatrix D;
667  this->giveLinearElasticMaterial()->give1dStressStiffMtrx(D, ElasticStiffness, gp, tStep);
668 
669  // perform plasticity return
670  performPlasticityReturn(gp, D, strainVector);
671 
672  // compute elastic strains and trial stress
673  FloatArray effectiveStress;
674  FloatArray elasticStrain = strainVector;
675  elasticStrain.subtract( status->giveTempPlasticStrain() );
676  effectiveStress.beProductOf(D, elasticStrain);
677 
678 
679  FloatArray effectiveStressTension;
680  FloatArray effectiveStressCompression;
681  double alpha = 0.;
682  if ( effectiveStress.at(1) >= 0 ) { //1D tensile stress state
683  alpha = 0.;
684  effectiveStressTension = effectiveStress;
685  effectiveStressCompression.zero();
686  } else if ( effectiveStress.at(1) < 0 ) { //1D compressive stress state
687  alpha = 1.;
688  effectiveStressTension.zero();
689  effectiveStressCompression = effectiveStress;
690  }
691 
692  FloatArray damages;
693  computeDamage(damages, strainVector, D, deltaTime, gp, tStep, alpha);
694 
695  //Split damage in a tension and compression part
696 
697  if ( isotropicFlag == 0 ) { //Default
698  effectiveStressTension.times( 1. - damages.at(1) );
699  effectiveStressCompression.times( 1. - damages.at(2) );
700  answer = effectiveStressTension;
701  answer.add(effectiveStressCompression);
702  } else { //Consider only tensile damage. Reduction to a fully isotropic model
703  answer = effectiveStress;
704  answer.times( 1. - damages.at(1) );
705  }
706 
707  status->letTempStrainVectorBe(strainVector);
708  status->letTempAlphaBe(alpha);
709  status->letTempStressVectorBe(answer);
710 #ifdef keep_track_of_dissipated_energy
711  double gf = pow(ft, 2) / this->eM; //rough estimation only for this purpose
712  status->computeWork(gp, gf);
713 #endif
714  assignStateFlag(gp);
715 }
716 
717 
718 void
720  GaussPoint *gp,
721  const FloatArray &strainVector,
722  TimeStep *tStep)
723 {
724  ConcreteDPM2Status *status = static_cast< ConcreteDPM2Status * >( this->giveStatus(gp) );
725 
726  // Initialize temp variables for this gauss point
727  status->initTempStatus();
728 
729  status->letTempStrainVectorBe(strainVector);
730 
731  //Calculate strain rate
732  //Time step
733  double deltaTime = 1.;
734  if ( strainRateFlag == 1 ) {
735  if ( tStep->giveTimeIncrement() == 0 ) { //Problem with the first step. For some reason the time increment is zero
736  deltaTime = this->timeFactor;
737  } else {
738  deltaTime = tStep->giveTimeIncrement() * this->timeFactor;
739  }
740  } else {
741  if ( tStep->giveTimeIncrement() == 0 ) {
742  deltaTime = 1.;
743  } else {
744  deltaTime = tStep->giveTimeIncrement();
745  }
746  }
747 
748  FloatMatrix D;
749  this->giveLinearElasticMaterial()->give3dMaterialStiffnessMatrix(D, ElasticStiffness, gp, tStep);
750 
751  // perform plasticity return
752  performPlasticityReturn(gp, D, strainVector);
753 
754 
755  // compute elastic strains and trial stress
756  FloatArray effectiveStress;
757  FloatArray elasticStrain = strainVector;
758  elasticStrain.subtract( status->giveTempPlasticStrain() );
759  effectiveStress.beProductOf(D, elasticStrain);
760 
761  FloatArray effectiveStressTension;
762  FloatArray effectiveStressCompression;
763  double alpha;
764 
765  alpha = computeAlpha(effectiveStressTension, effectiveStressCompression, effectiveStress);
766 
767  FloatArray damages;
768  computeDamage(damages, strainVector, D, deltaTime, gp, tStep, alpha);
769 
770  //Split damage in a tension and compression part
771 
772  if ( isotropicFlag == 0 ) { //Default
773  effectiveStressTension.times( 1. - damages.at(1) );
774  effectiveStressCompression.times( 1. - damages.at(2) );
775  answer = effectiveStressTension;
776  answer.add(effectiveStressCompression);
777  } else { //Consider only tensile damage. Reduction to a fully isotropic model
778  answer = effectiveStress;
779  answer.times( 1. - damages.at(1) );
780  }
781 
782  status->letTempStrainVectorBe(strainVector);
783  status->letTempAlphaBe(alpha);
784  status->letTempStressVectorBe(answer);
785 #ifdef keep_track_of_dissipated_energy
786  double gf = pow(ft, 2) / this->eM; //rough estimation only for this purpose
787  status->computeWork(gp, gf);
788 #endif
789  assignStateFlag(gp);
790 }
791 
792 
793 void
795  const FloatArray &strain,
796  const FloatMatrix &D,
797  double deltaTime,
798  GaussPoint *gp,
799  TimeStep *tStep,
800  double tempAlpha)
801 {
802  ConcreteDPM2Status *status = static_cast< ConcreteDPM2Status * >( this->giveStatus(gp) );
803 
804  double tempEquivStrain;
805  double deltaPlasticStrainNorm;
806  double tempDamageTension = 0.0;
807  double tempDamageCompression = 0.0;
808 
809  double tempKappaDTension = 0.0, tempKappaDCompression = 0.0;
810  double tempKappaDTensionOne = 0.0, tempKappaDTensionTwo = 0.0;
811  double tempKappaDCompressionOne = 0.0, tempKappaDCompressionTwo = 0.0;
812 
813  double rateFactor;
814 
815  int unAndReloadingFlag = 0;
816  double minEquivStrain = 0.;
817 
818  unAndReloadingFlag = checkForUnAndReloading(tempEquivStrain, minEquivStrain, D, gp);
819 
820  if ( ( status->giveDamageTension() == 0. ) && ( status->giveDamageCompression() == 0. ) ) {
821  rateFactor = computeRateFactor(tempAlpha, deltaTime, gp, tStep);
822  } else {
823  rateFactor = status->giveRateFactor();
824  }
825 
826  //Compute equivalent strains for tension and compression
827  double tempEquivStrainTension = 0.;
828  double tempEquivStrainCompression = 0.;
829 
830  tempEquivStrainTension = status->giveEquivStrainTension() + ( tempEquivStrain - status->giveEquivStrain() ) / rateFactor;
831 
832  if ( unAndReloadingFlag == 0 ) { //Standard way
833  tempEquivStrainCompression = status->giveEquivStrainCompression() + ( tempAlpha * ( tempEquivStrain - status->giveEquivStrain() ) ) / rateFactor;
834  } else {
835  tempEquivStrainCompression = status->giveEquivStrainCompression() + status->giveAlpha() * ( minEquivStrain - status->giveEquivStrain() ) / rateFactor + ( tempAlpha * ( tempEquivStrain - minEquivStrain ) ) / rateFactor;
836  }
837 
838 
839  //If damage threshold is exceeded determine the rate factor from the previous step
840  if ( ( tempEquivStrainTension > e0 || tempEquivStrainCompression > e0 ) &&
841  ( ( status->giveDamageTension() == 0. ) && ( status->giveDamageCompression() == 0. ) ) && !tStep->isTheFirstStep() ) {
842  //Rate factor from last step
843  rateFactor = status->giveRateFactor();
844 
845  tempEquivStrainTension = status->giveEquivStrainTension() + ( tempEquivStrain - status->giveEquivStrain() ) / rateFactor;
846  if ( unAndReloadingFlag == 0 ) { //Standard way
847  tempEquivStrainCompression = status->giveEquivStrainCompression() + ( tempAlpha * ( tempEquivStrain - status->giveEquivStrain() ) ) / rateFactor;
848  } else {
849  tempEquivStrainCompression = status->giveEquivStrainCompression() + status->giveAlpha() * ( minEquivStrain - status->giveEquivStrain() ) / rateFactor + ( tempAlpha * ( tempEquivStrain - minEquivStrain ) ) / rateFactor;
850  }
851  }
852 
853  status->letTempRateFactorBe(rateFactor);
854 
855  double fTension = tempEquivStrainTension - status->giveKappaDTension();
856  double fCompression = tempEquivStrainCompression - status->giveKappaDCompression();
857 
858  //Normalize the fs
859  fTension = fTension / e0;
860  fCompression = fCompression / e0;
861 
862  double ductilityMeasure = computeDuctilityMeasureDamage(strain, gp);
863  double deltaPlasticStrainNormTension, deltaPlasticStrainNormCompression;
864 
865  if ( fTension < -yieldTolDamage && fCompression < -yieldTolDamage ) {
866  //Neither tension nor compression is active
867 
868  tempKappaDTension = status->giveKappaDTension();
869  tempKappaDTensionOne = status->giveKappaDTensionOne();
870  tempKappaDTensionTwo = status->giveKappaDTensionTwo();
871 
872  tempKappaDCompression = status->giveKappaDCompression();
873  tempKappaDCompressionOne = status->giveKappaDCompressionOne();
874  tempKappaDCompressionTwo = status->giveKappaDCompressionTwo();
875 
876  tempDamageTension = status->giveDamageTension();
877  tempDamageCompression = status->giveDamageCompression();
878  } else if ( fTension >= -yieldTolDamage && fCompression < -yieldTolDamage ) { //Only tension is active
879  //Update tension history variables
880  tempKappaDTension = tempEquivStrainTension;
881  deltaPlasticStrainNorm = computeDeltaPlasticStrainNormTension(tempKappaDTension, status->giveKappaDTension(), gp);
882  tempKappaDTensionOne = status->giveKappaDTensionOne() + deltaPlasticStrainNorm / ductilityMeasure / rateFactor;
883  tempKappaDTensionTwo = status->giveKappaDTensionTwo() + ( tempKappaDTension - status->giveKappaDTension() ) / ductilityMeasure;
884 
885  //Nothing changes for compression history variables
886  tempKappaDCompression = status->giveKappaDCompression();
887  tempKappaDCompressionOne = status->giveKappaDCompressionOne();
888  tempKappaDCompressionTwo = status->giveKappaDCompressionTwo();
889 
890  //Initialise damage with tensile history variable
891  this->initDamaged(tempKappaDTension, strain, gp);
892 
893  tempDamageTension = computeDamageParamTension( tempKappaDTension, tempKappaDTensionOne, tempKappaDTensionTwo, status->giveLe(), status->giveDamageTension() );
894 
895  tempDamageCompression = status->giveDamageCompression();
896  } else if ( fTension < -yieldTolDamage && fCompression >= -yieldTolDamage ) {
897  //Only compression is active
898 
899  //Nothing changes for the history variables in tension
900  tempKappaDTension = status->giveKappaDTension();
901  tempKappaDTensionOne = status->giveKappaDTensionOne();
902  tempKappaDTensionTwo = status->giveKappaDTensionTwo();
903 
904  //Update compression history variables
905  tempKappaDCompression = tempEquivStrainCompression;
906  deltaPlasticStrainNormCompression = computeDeltaPlasticStrainNormCompression(tempAlpha, tempKappaDCompression, status->giveKappaDCompression(), gp);
907  tempKappaDCompressionOne = status->giveKappaDCompressionOne() + deltaPlasticStrainNormCompression / ( ductilityMeasure * rateFactor );
908  tempKappaDCompressionTwo = status->giveKappaDCompressionTwo() + ( tempKappaDCompression - status->giveKappaDCompression() ) / ductilityMeasure;
909 
910  //Determine damage parameters
911  tempDamageTension = status->giveDamageTension();
912  tempDamageCompression = computeDamageParamCompression( tempKappaDCompression, tempKappaDCompressionOne, tempKappaDCompressionTwo, status->giveDamageCompression() );
913  } else if ( fTension >= -yieldTolDamage && fCompression >= -yieldTolDamage ) {
914  //Both tension and compression is active
915 
916  //Update tension history variables
917  tempKappaDTension = tempEquivStrainTension;
918  deltaPlasticStrainNormTension = computeDeltaPlasticStrainNormTension(tempKappaDTension, status->giveKappaDTension(), gp);
919  tempKappaDTensionOne = status->giveKappaDTensionOne() + deltaPlasticStrainNormTension / ( ductilityMeasure * rateFactor );
920  tempKappaDTensionTwo = status->giveKappaDTensionTwo() + ( tempKappaDTension - status->giveKappaDTension() ) / ductilityMeasure;
921 
922  //Update the compression history variables
923  tempKappaDCompression = tempEquivStrainCompression;
924  deltaPlasticStrainNormCompression =
925  computeDeltaPlasticStrainNormCompression(tempAlpha, tempKappaDCompression, status->giveKappaDCompression(), gp);
926  tempKappaDCompressionOne = status->giveKappaDCompressionOne() + deltaPlasticStrainNormCompression / ( ductilityMeasure * rateFactor );
927  tempKappaDCompressionTwo = status->giveKappaDCompressionTwo() +
928  ( tempKappaDCompression - status->giveKappaDCompression() ) / ductilityMeasure;
929 
930  //Determine the damage parameters
931  this->initDamaged(tempKappaDTension, strain, gp);
932 
933  tempDamageTension = computeDamageParamTension( tempKappaDTension, tempKappaDTensionOne, tempKappaDTensionTwo, status->giveLe(), status->giveDamageTension() );
934 
935  tempDamageCompression = computeDamageParamCompression( tempKappaDCompression, tempKappaDCompressionOne, tempKappaDCompressionTwo, status->giveDamageCompression() );
936  }
937 
938  //Write all temp history variables to the status
939  status->letTempEquivStrainBe(tempEquivStrain);
940 
941  //Tension
942  status->letTempEquivStrainTensionBe(tempEquivStrainTension);
943  status->letTempKappaDTensionBe(tempKappaDTension);
944  status->letTempKappaDTensionOneBe(tempKappaDTensionOne);
945  status->letTempKappaDTensionTwoBe(tempKappaDTensionTwo);
946  status->letTempDamageTensionBe(tempDamageTension);
947 
948  //Compression
949  status->letTempEquivStrainCompressionBe(tempEquivStrainCompression);
950  status->letTempKappaDCompressionBe(tempKappaDCompression);
951  status->letTempKappaDCompressionOneBe(tempKappaDCompressionOne);
952  status->letTempKappaDCompressionTwoBe(tempKappaDCompressionTwo);
953  status->letTempDamageCompressionBe(tempDamageCompression);
954 
955  answer.resize(2);
956  answer.at(1) = tempDamageTension;
957  answer.at(2) = tempDamageCompression;
958 }
959 
960 int
962  double &minEquivStrain,
963  const FloatMatrix &D,
964  GaussPoint *gp)
965 {
966  ConcreteDPM2Status *status = static_cast< ConcreteDPM2Status * >( this->giveStatus(gp) );
967 
968  double sigEffective, rhoEffective, thetaEffective;
969 
970  //Access old and new strains
971  FloatArray oldStrain = status->giveStrainVector();
972  FloatArray strain = status->giveTempStrainVector();
973 
974  //Compute the temp equivalent strain
975  FloatArray tempEffectiveStress;
976  FloatArray tempElasticStrain = strain;
977  tempElasticStrain.subtract( status->giveTempPlasticStrain() );
978  tempEffectiveStress.beProductOf(D, tempElasticStrain);
979  computeTrialCoordinates(tempEffectiveStress, sigEffective, rhoEffective, thetaEffective);
980  tempEquivStrain = computeEquivalentStrain(sigEffective, rhoEffective, thetaEffective);
981  //Get the equivalent strain from the status
982  double equivStrain = status->giveEquivStrain();
983 
984  //Compute the increment of effective stress
985  FloatArray effectiveStress;
986  FloatArray elasticStrain = oldStrain;
987  for ( int i = 1; i <= elasticStrain.giveSize(); ++i ) elasticStrain.at(i) -= status->givePlasticStrain().at(i);
988  //elasticStrain.subtract( status->givePlasticStrain() );
989  effectiveStress.beProductOf(D, elasticStrain);
990  FloatArray deltaEffectiveStress;
991  deltaEffectiveStress = tempEffectiveStress;
992  deltaEffectiveStress.subtract(effectiveStress);
993 
994  //Compute equivalent strains for stress state slightly greater than the effective stress and smaller than the temp effective stress
995  FloatArray intermediateEffectiveStress;
996  intermediateEffectiveStress = effectiveStress;
997  //For slightly more than effective stress
998  intermediateEffectiveStress.add(0.01 * deltaEffectiveStress);
999  computeTrialCoordinates(intermediateEffectiveStress, sigEffective, rhoEffective, thetaEffective);
1000  double equivStrainPlus = computeEquivalentStrain(sigEffective, rhoEffective, thetaEffective);
1001  //For slightly less than temp effective stress
1002  intermediateEffectiveStress = effectiveStress;
1003  intermediateEffectiveStress.add(0.99 * deltaEffectiveStress);
1004  computeTrialCoordinates(intermediateEffectiveStress, sigEffective, rhoEffective, thetaEffective);
1005  double tempEquivStrainMinus = computeEquivalentStrain(sigEffective, rhoEffective, thetaEffective);
1006 
1007  //Check for unloading and reloading in the same step
1008  int unloadingFlag = 0;
1009  minEquivStrain = equivStrain;
1010  double midEquivStrain = 0.;
1011 
1012  unloadingFlag = 0;
1013  if ( ( equivStrain > equivStrainPlus && tempEquivStrain > tempEquivStrainMinus ) && ( fabs(equivStrainPlus - equivStrain) > yieldTolDamage / 100. && fabs(tempEquivStrainMinus - tempEquivStrain) > yieldTolDamage / 100. ) ) {
1014  unloadingFlag = 1;
1015  //Unloading and reloading takes place. Find the minimum equivalent strain by subincrementing the effective stress increment
1016  for ( double k = 1.0; k <= 100.0; k = k + 1.0 ) {
1017  intermediateEffectiveStress = effectiveStress;
1018  intermediateEffectiveStress.add(k / 100. * deltaEffectiveStress);
1019  computeTrialCoordinates(intermediateEffectiveStress, sigEffective, rhoEffective, thetaEffective);
1020  midEquivStrain = computeEquivalentStrain(sigEffective, rhoEffective, thetaEffective);
1021 
1022  if ( midEquivStrain <= minEquivStrain ) {
1023  minEquivStrain = midEquivStrain;
1024  } else {
1025  return unloadingFlag;
1026  }
1027  }
1028  }
1029  return unloadingFlag;
1030 }
1031 
1032 
1033 double
1035  double deltaTime,
1036  GaussPoint *gp,
1037  TimeStep *tStep)
1038 {
1039  if ( strainRateFlag == 0 ) {
1040  return 1;
1041  }
1042 
1043  ConcreteDPM2Status *status = static_cast< ConcreteDPM2Status * >( this->giveStatus(gp) );
1044 
1045  //Access old and new strain
1046  FloatArray strain = status->giveTempStrainVector();
1047 
1048  //Determine the principal values of the strain
1049  FloatArray principalStrain;
1051  //Determine max and min value;
1052  double maxStrain = -1.e20, minStrain = 1.e20;
1053  for ( int k = 1; k <= principalStrain.giveSize(); k++ ) {
1054  //maximum
1055  if ( principalStrain.at(k) > maxStrain ) {
1056  maxStrain = principalStrain.at(k);
1057  }
1058 
1059  //minimum
1060  if ( principalStrain.at(k) < minStrain ) {
1061  minStrain = principalStrain.at(k);
1062  }
1063  }
1064 
1065  //Evaluate the equivalent strains
1066  double strainRate;
1067  double oldRateStrain = status->giveRateStrain();
1068  if ( 1. - alpha > DYNCON_TOL ) { //Tension
1069  strainRate = ( maxStrain - oldRateStrain ) / deltaTime;
1070  status->letTempRateStrainBe(maxStrain);
1071  } else { //Compression
1072  strainRate = ( minStrain - oldRateStrain ) / deltaTime;
1073  status->letTempRateStrainBe(minStrain);
1074  }
1075 
1076  //For tension
1077  double deltaS = 1. / ( 1. + 8. * this->fc / this->fcZero );
1078  double betaS = exp( ( 6 * deltaS - 2. ) * log(10.) );
1079 
1080 
1081  double strainRateZeroTension = 1.e-6;
1082  double strainRateZeroCompression = -30.e-6;
1083 
1084  //For compression
1085  // Fip-Model code 1990 expressions modified to take into account that we have an equivalent strain!
1086  double alphaS = 1. / ( 5. + 9 * this->fc / this->fcZero );
1087  double gammaS = exp( ( 6.156 * alphaS - 2. ) * log(10.0) );
1088 
1089  double strainRateRatioTension, strainRateRatioCompression;
1090 
1091  double rateFactorTension = 1.;
1092  double rateFactorCompression = 1.;
1093 
1094  strainRateRatioTension = strainRate / strainRateZeroTension;
1095 
1096  //Tension
1097  if ( strainRate < 30.e-6 ) {
1098  rateFactorTension = 1.;
1099  } else if ( 30.e-6 < strainRate && strainRate < 1 ) {
1100  rateFactorTension = pow(strainRateRatioTension, deltaS);
1101  } else {
1102  rateFactorTension = betaS * pow(strainRateRatioTension, 0.333);
1103  }
1104 
1105  //Compression
1106  strainRateRatioCompression = strainRate / strainRateZeroCompression;
1107  if ( strainRate > -30.e-6 ) {
1108  rateFactorCompression = 1.;
1109  } else if ( -30.e-6 > strainRate && strainRate > -30 ) {
1110  rateFactorCompression = pow(strainRateRatioCompression, 1.026 * alphaS);
1111  } else {
1112  rateFactorCompression = gammaS * pow(strainRateRatioCompression, 0.333);
1113  }
1114 
1115  double rateFactor = ( 1. - alpha ) * rateFactorTension + alpha * rateFactorCompression;
1116 
1117  return rateFactor;
1118 }
1119 
1120 
1121 
1122 double
1124 {
1125  ConcreteDPM2Status *status = static_cast< ConcreteDPM2Status * >( this->giveStatus(gp) );
1126 
1127  const FloatArray &tempPlasticStrain = status->giveTempPlasticStrain();
1128  const FloatArray &plasticStrain = status->givePlasticStrain();
1129 
1130  FloatArray deltaPlasticStrain = tempPlasticStrain;
1131  for ( int i = 1; i <= deltaPlasticStrain.giveSize(); ++i ) deltaPlasticStrain.at(i) -= plasticStrain.at(i);
1132  //deltaPlasticStrain.subtract(plasticStrain);
1133 
1134  double deltaPlasticStrainNorm = 0;
1135 
1136  //Distinguish pre-peak, peak and post-peak
1137 
1138  double factor = 0.;
1139  if ( tempKappaD < e0 * ( 1. - yieldTolDamage ) ) {
1140  deltaPlasticStrainNorm = 0.;
1141  } else if ( tempKappaD > e0 * ( 1. - yieldTolDamage ) && kappaD < e0 * ( 1. - yieldTolDamage ) ) {
1142  factor = ( 1. - ( e0 - kappaD ) / ( tempKappaD - kappaD ) );
1143  deltaPlasticStrain.times(factor);
1144  deltaPlasticStrainNorm = deltaPlasticStrain.computeNorm();
1145  } else {
1146  deltaPlasticStrainNorm = deltaPlasticStrain.computeNorm();
1147  }
1148 
1149  return deltaPlasticStrainNorm;
1150 }
1151 
1152 
1153 
1154 double
1155 ConcreteDPM2 :: computeDeltaPlasticStrainNormCompression(double tempAlpha, double tempKappaD, double kappaD, GaussPoint *gp)
1156 {
1157  ConcreteDPM2Status *status = static_cast< ConcreteDPM2Status * >( this->giveStatus(gp) );
1158 
1159  const FloatArray &tempPlasticStrain = status->giveTempPlasticStrain();
1160  const FloatArray &plasticStrain = status->givePlasticStrain();
1161 
1162  FloatArray deltaPlasticStrain;
1163  deltaPlasticStrain.add(tempAlpha, tempPlasticStrain);
1164  for ( int i = 1; i <= deltaPlasticStrain.giveSize(); ++i ) deltaPlasticStrain.at(i) -= tempAlpha * plasticStrain.at(i);
1165  //deltaPlasticStrain.add(-tempAlpha, plasticStrain);
1166 
1167  double deltaPlasticStrainNorm = 0;
1168 
1169  //Distinguish pre-peak, peak and post-peak
1170  if ( tempKappaD < e0 * ( 1. - yieldTolDamage ) ) {
1171  deltaPlasticStrainNorm = 0.;
1172  } else if ( tempKappaD > e0 * ( 1. - yieldTolDamage ) && kappaD < e0 * ( 1. - yieldTolDamage ) ) {
1173  double factor = ( 1. - ( e0 - kappaD ) / ( tempKappaD - kappaD ) );
1174  deltaPlasticStrain.times(factor);
1175  deltaPlasticStrainNorm = deltaPlasticStrain.computeNorm();
1176  } else {
1177  deltaPlasticStrainNorm = deltaPlasticStrain.computeNorm();
1178  }
1179 
1180  double tempKappaP = status->giveTempKappaP();
1181  double yieldHardTwo = computeHardeningTwo(tempKappaP);
1182  double extraFactor;
1183  if ( rho < 1.e-16 ) {
1184  extraFactor = this->ft * yieldHardTwo * sqrt(2. / 3.) / 1.e-16 / sqrt( 1. + 2. * pow(this->dilationConst, 2.) );
1185  } else {
1186  extraFactor = this->ft * yieldHardTwo * sqrt(2. / 3.) / this->rho / sqrt( 1. + 2. * pow(this->dilationConst, 2.) );
1187  }
1188 
1189  return deltaPlasticStrainNorm * extraFactor;
1190 }
1191 
1192 
1193 double
1195  double rho,
1196  double theta)
1197 {
1198  double rFunction = ( 4. * ( 1. - pow(this->ecc, 2.) ) * pow(cos(theta), 2.) + pow(2. * this->ecc - 1., 2.) ) / ( 2. * ( 1. - pow(this->ecc, 2.) ) * cos(theta) + ( 2. * this->ecc - 1. ) * sqrt(4. * ( 1. - pow(this->ecc, 2.) ) * pow(cos(theta), 2.) + 5. * pow(this->ecc, 2.) - 4. * this->ecc) );
1199 
1200  double pHelp = -this->m * ( rho * rFunction / ( sqrt(6.) * fc ) + sig / fc );
1201 
1202  double qHelp = -3. / 2. * pow(rho, 2.) / pow(this->fc, 2.);
1203 
1204  double help = -0.5 * pHelp + sqrt(pow(pHelp, 2.) / 4. - qHelp);
1205 
1206  double tempEquivStrain = 0.;
1207  if ( help > 0 ) {
1208  tempEquivStrain = help * e0;
1209  }
1210  return tempEquivStrain;
1211 }
1212 
1213 
1214 double
1215 ConcreteDPM2 :: computeDamageParamTension(double equivStrain, double kappaOne, double kappaTwo, double le, double omegaOld)
1216 {
1217  double omega = 0.;
1218 
1219  double help;
1220  if ( equivStrain > e0 * ( 1. - yieldTolDamage ) ) {
1221  if ( softeningType == 0 ) { //linear
1222  omega = ( this->eM * equivStrain * this->wf - this->ft * this->wf + this->ft * kappaOne * le ) /
1223  ( this->eM * equivStrain * this->wf - this->ft * le * kappaTwo );
1224  } else if ( softeningType == 1 ) { //bilinear: Calculate damage parameter for both parts of bilinear curve and check which fulfils limits.
1225  omega = ( this->eM * equivStrain * this->wfOne - this->ft * this->wfOne - ( this->ftOne - this->ft ) * kappaOne * le ) /
1226  ( this->eM * equivStrain * this->wfOne + ( this->ftOne - this->ft ) * le * kappaTwo );
1227  help = le * kappaOne + le * omega * kappaTwo;
1228 
1229  if ( help >= 0. && help < this->wfOne ) {
1230  return omega;
1231  }
1232 
1233  omega = ( this->eM * equivStrain * ( this->wf - this->wfOne ) - this->ftOne * ( this->wf - this->wfOne ) +
1234  this->ftOne * kappaOne * le - this->ftOne * this->wfOne ) /
1235  ( this->eM * equivStrain * ( this->wf - this->wfOne ) - this->ftOne * le * kappaTwo );
1236  help = le * kappaOne + le * omega * kappaTwo;
1237 
1238  if ( help > this->wfOne && help < this->wf ) {
1239  return omega;
1240  }
1241  } else if ( softeningType == 2 ) { //exponential: Iterative solution
1242  omega = 1.; //Initial guess
1243  double residual = 0.;
1244  double dResidualDOmega = 0.;
1245  int nite = 0;
1246 
1247  do {
1248  nite++;
1249 
1250 
1251  residual = ( 1 - omega ) * eM * equivStrain - ft *exp(-le * ( omega * kappaTwo + kappaOne ) / wf);
1252  dResidualDOmega = -eM * equivStrain + ft * le * kappaTwo / wf *exp(-le * ( omega * kappaTwo + kappaOne ) / wf);
1253 
1254  omega -= residual / dResidualDOmega;
1255  if ( nite > newtonIter ) {
1256  OOFEM_ERROR("algorithm not converging");
1257  }
1258  } while ( fabs(residual / ft) >= 1.e-8 );
1259  }
1260  } else {
1261  omega = 0.;
1262  }
1263 
1264 
1265  if ( omega > 1. ) {
1266  omega = 1.;
1267  }
1268 
1269  if ( omega < 0. || omega < omegaOld ) {
1270  omega = omegaOld;
1271  }
1272 
1273 
1274  return omega;
1275 }
1276 
1277 double
1278 ConcreteDPM2 :: computeDamageParamCompression(double equivStrain, double kappaOne, double kappaTwo, double omegaOld)
1279 {
1280  if ( isotropicFlag == 1 ) {
1281  return 0.;
1282  }
1283 
1284 
1285  double omega = 1.;
1286  int nite = 0;
1287  double residual = 0.;
1288  double dResidualDOmega = 0.;
1289  if ( equivStrain > e0 * ( 1. - yieldTolDamage ) ) {
1290  do {
1291  nite++;
1292 
1293  residual = ( 1. - omega ) * this->eM * equivStrain - this->ft *exp(-( kappaOne + omega * kappaTwo ) / this->efCompression);
1294  dResidualDOmega = -this->eM * equivStrain + this->ft * kappaTwo / this->efCompression *exp(-( kappaOne + omega * kappaTwo ) / this->efCompression);
1295 
1296  omega -= residual / dResidualDOmega;
1297  if ( nite > newtonIter ) {
1298  OOFEM_ERROR("algorithm not converging");
1299  }
1300  } while ( fabs(residual / ft) >= 1.e-8 );
1301  } else {
1302  omega = 0.;
1303  }
1304 
1305  if ( omega > 1. ) {
1306  omega = 1.;
1307  }
1308  if ( omega < omegaOld || omega < 0. ) {
1309  omega = omegaOld;
1310  }
1311 
1312  return omega;
1313 }
1314 
1315 
1316 
1317 void
1319  const FloatArray &strain,
1320  GaussPoint *gp)
1321 {
1322  if ( kappaD <= 0. ) {
1323  return;
1324  }
1325 
1326  int indx = 1;
1327  double le;
1328  FloatArray principalStrains, crackPlaneNormal(3);
1329  FloatMatrix principalDir;
1330  ConcreteDPM2Status *status = static_cast< ConcreteDPM2Status * >( this->giveStatus(gp) );
1331 
1332  if ( helem > 0. ) {
1333  status->setLe(helem);
1334  } else if ( strain.giveSize() == 1 ) {
1335  le = gp->giveElement()->computeLength();
1336  status->setLe(le);
1337  } else if ( status->giveDamageTension() == 0. && status->giveDamageCompression() == 0. ) {
1338  StructuralMaterial :: computePrincipalValDir(principalStrains, principalDir, strain, principal_strain);
1339  // find index of max positive principal strain
1340  for ( int i = 2; i <= 3; i++ ) {
1341  if ( principalStrains.at(i) > principalStrains.at(indx) ) {
1342  indx = i;
1343  }
1344  }
1345 
1346  for ( int i = 1; i <= 3; i++ ) {
1347  crackPlaneNormal.at(i) = principalDir.at(i, indx);
1348  }
1349 
1350  // evaluate the projected element size
1351  le = gp->giveElement()->giveCharacteristicLength(crackPlaneNormal);
1352  if ( le == 0. ) {
1353  le = gp->giveElement()->computeMeanSize();
1354  }
1355 
1356  // store le in the corresponding status
1357  status->setLe(le);
1358  } else if ( status->giveLe() == 0. ) {
1359  // this happens if the status is initialized from a file
1360  // with nonzero damage
1361  // le determined as square root of element area or cube root of el. volume
1362  le = gp->giveElement()->computeMeanSize();
1363  status->setLe(le);
1364  }
1365 }
1366 
1367 
1368 double
1370 {
1371  //Angle in uniaxial compression is atan(1./sqrt(6.))=0.387597
1372  double alphaZero = 0.40824829;
1373 
1374  double Rs = 0;
1375  if ( this->sig < 0. ) {
1376  if ( this->rho > 1.e-16 ) {
1377  Rs = -this->sig / ( alphaZero * this->rho );
1378  } else { //Need to set a dummy valye
1379  Rs = -this->sig * 1.e16 / alphaZero;
1380  }
1381  } else {
1382  Rs = 0;
1383  }
1384 
1385  return 1. + ( ASoft - 1. ) * Rs; // ductilityMeasure
1386 }
1387 
1388 
1389 void
1391  const FloatMatrix &D,
1392  const FloatArray &strain)
1393 {
1394  ConcreteDPM2Status *status = static_cast< ConcreteDPM2Status * >( this->giveStatus(gp) );
1395 
1396  FloatMatrix C;
1397  C.beInverseOf(D);
1398 
1399  //get temp plastic strain and tempKappa
1400  FloatArray tempPlasticStrain = status->givePlasticStrain();
1401  double tempKappaP = status->giveKappaP();
1402 
1403  // compute elastic strains and trial stress
1404  FloatArray effectiveStress;
1405  FloatArray elasticStrain = strain;
1406  for ( int i = 1; i <= elasticStrain.giveSize(); ++i ) elasticStrain.at(i) -= tempPlasticStrain.at(i);
1407  effectiveStress.beProductOf(D, elasticStrain);
1408 
1409  //Compute trial coordinates
1410  computeTrialCoordinates(effectiveStress, this->sig, this->rho, this->thetaTrial);
1411 
1412  double yieldValue;
1413 
1414  FloatArray convergedStrain;
1415 
1416  //What happens here to thermal strains that are subtracted first?
1417  FloatArray oldStrain = status->giveStrainVector();
1418  FloatArray tempStrain;
1419  FloatArray deltaStrain;
1420 
1421  // introduce a strange subincrementation flag
1422  int subIncrementFlag = 0;
1423 
1424  double apexStress = 0.;
1425  int subincrementcounter = 0;
1426  //Attempt to implement subincrementation
1427  // initialize variables
1428  subIncrementFlag = 0;
1429  convergedStrain = oldStrain;
1430  tempStrain = strain;
1431  deltaStrain.beDifferenceOf(strain, oldStrain);
1432  //To get into the loop
1434  while ( returnResult == RR_NotConverged || subIncrementFlag == 1 ) {
1435  elasticStrain = tempStrain;
1436  for ( int i = 1; i <= elasticStrain.giveSize(); ++i ) elasticStrain.at(i) -= tempPlasticStrain.at(i);
1437  //elasticStrain.subtract(tempPlasticStrain);
1438  effectiveStress.beProductOf(D, elasticStrain);
1439  double theta;
1440  computeTrialCoordinates(effectiveStress, this->sig, this->rho, theta);
1441  yieldValue = computeYieldValue(this->sig, this->rho, this->thetaTrial, tempKappaP);
1442 
1443  apexStress = 0.;
1444 
1445  if ( yieldValue > 0. ) {
1446  checkForVertexCase(apexStress, sig, tempKappaP, strain.giveSize() == 1);
1448  tempKappaP = performVertexReturn(effectiveStress, apexStress, tempKappaP, gp);
1449  status->letTempKappaPBe(tempKappaP);
1450  }
1451  if ( returnType == RT_Regular ) {
1452  tempKappaP = performRegularReturn(effectiveStress, tempKappaP, gp);
1453  status->letTempKappaPBe(tempKappaP);
1454  }
1455  } else {
1457  tempPlasticStrain = status->givePlasticStrain();
1458  status->letTempPlasticStrainBe(tempPlasticStrain);
1459  status->letTempKappaPBe(tempKappaP);
1460  break;
1461  }
1462 
1463  if ( returnResult == RR_NotConverged ) {
1464  subincrementcounter++;
1465  if ( subincrementcounter > 10 ) {
1466  OOFEM_LOG_INFO( "Unstable element %d \n", gp->giveElement()->giveGlobalNumber() );
1467  OOFEM_LOG_INFO( "Old strain vector %g %g %g %g %g %g \n", oldStrain.at(1), oldStrain.at(2), oldStrain.at(3), oldStrain.at(4), oldStrain.at(5), oldStrain.at(6) );
1468  FloatArray help = status->giveTempPlasticStrain();
1469  FloatArray help1;
1470  double sig1, rho1, theta1;
1471  oldStrain.subtract(help);
1472  help1.beProductOf(D, oldStrain);
1473  OOFEM_LOG_INFO( "Old plastic strain vector %g %g %g %g %g %g \n", help.at(1), help.at(2), help.at(3), help.at(4), help.at(5), help.at(6) );
1474  OOFEM_LOG_INFO( "New strain vector %g %g %g %g %g %g \n", strain.at(1), strain.at(2), strain.at(3), strain.at(4), strain.at(5), strain.at(6) );
1475  computeTrialCoordinates(effectiveStress, this->sig, this->rho, theta);
1476  computeTrialCoordinates(help1, sig1, rho1, theta1);
1477  yieldValue = computeYieldValue(this->sig, this->rho, this->thetaTrial, tempKappaP);
1478  OOFEM_LOG_INFO("OLD Sig %g rho %g theta %g \n", sig1, rho1, theta1);
1479  OOFEM_LOG_INFO("NEW Sig %g rho %g theta %g \n", sig, rho, theta);
1481  OOFEM_LOG_INFO("Vertex case apexstress %g\n", apexStress);
1482  } else {
1483  OOFEM_LOG_INFO("Regular case %g \n", 15.18);
1484  }
1485  OOFEM_LOG_INFO("KappaP old %g new %g yieldfun %g\n", status->giveTempKappaP(), tempKappaP, yieldValue);
1486  OOFEM_ERROR("Could not reach convergence with small deltaStrain, giving up.");
1487  } else if ( subincrementcounter > 9 && tempKappaP < 1. ) {
1488  tempKappaP = 1.;
1489  status->letTempKappaPBe(tempKappaP);
1490  }
1491 
1492  subIncrementFlag = 1;
1493  deltaStrain.times(0.5);
1494  tempStrain = convergedStrain;
1495  tempStrain.add(deltaStrain);
1496  } else if ( returnResult == RR_Converged && subIncrementFlag == 0 ) {
1497  elasticStrain.beProductOf(C, effectiveStress);
1498  tempPlasticStrain = strain;
1499  tempPlasticStrain.subtract(elasticStrain);
1500  status->letTempPlasticStrainBe(tempPlasticStrain);
1501  } else if ( returnResult == RR_Converged && subIncrementFlag == 1 ) {
1502  OOFEM_LOG_INFO("Subincrementation %d required\n", subincrementcounter);
1503  subincrementcounter = 0;
1504  elasticStrain.beProductOf(C, effectiveStress);
1505  tempPlasticStrain = tempStrain;
1506  tempPlasticStrain.subtract(elasticStrain);
1507  status->letTempPlasticStrainBe(tempPlasticStrain);
1508 
1509  subIncrementFlag = 0;
1511  convergedStrain = tempStrain;
1512  deltaStrain.beDifferenceOf(strain, convergedStrain);
1513  tempStrain = strain;
1514  }
1515  }
1516 }
1517 
1518 
1519 bool
1521  double sig,
1522  double tempKappa,
1523  bool mat1d)
1524 {
1525  if ( mat1d ) {
1527  return false;
1528  }
1529 
1530  answer = 0.;
1531  if ( sig > 0. ) {
1533  } else if ( sig < 0. && tempKappa < 1. ) {
1535  } else {
1537  }
1538  return false;
1539 }
1540 
1541 
1542 
1543 double
1545  double apexStress, double tempKappaP,
1546  GaussPoint *gp)
1547 {
1548  FloatArray deviatoricStressTrial;
1549  double sigTrial;
1550  double yieldValue = 0.;
1551  double yieldValueMid = 0.;
1552  double sig2 = 0.;
1553  double dSig;
1554  double sigMid;
1555  double sigAnswer;
1556  double ratioPotential;
1557 
1558  sigTrial = computeDeviatoricVolumetricSplit(deviatoricStressTrial, effectiveStress);
1559  double rhoTrial = computeSecondCoordinate(deviatoricStressTrial);
1560 
1561  double kappaInitial = tempKappaP;
1562 
1563  sig2 = apexStress;
1564 
1565  tempKappaP =
1566  computeTempKappa(kappaInitial, sigTrial, rhoTrial, sigTrial);
1567 
1568  yieldValue =
1569  computeYieldValue(sigTrial, 0., 0., tempKappaP);
1570 
1571  tempKappaP =
1572  computeTempKappa(kappaInitial, sigTrial, rhoTrial, sig2);
1573 
1574  yieldValueMid =
1575  computeYieldValue(sig2, 0., 0., tempKappaP);
1576 
1577  if ( yieldValue * yieldValueMid >= 0. ) {
1580  return kappaInitial;
1581  }
1582 
1583  if ( yieldValue < 0.0 ) {
1584  dSig = sig2 - sigTrial;
1585  sigAnswer = sig2;
1586  } else {
1587  dSig = sigTrial - sig2;
1588  sigAnswer = sig2;
1589  }
1590 
1591  for ( int j = 0; j < 250; j++ ) {
1592  dSig = 0.5 * dSig;
1593 
1594  sigMid = sigAnswer + dSig;
1595 
1596 
1597  tempKappaP =
1598  computeTempKappa(kappaInitial, sigTrial, rhoTrial, sigMid);
1599 
1600  yieldValueMid = computeYieldValue(sigMid, 0., 0., tempKappaP);
1601 
1602  if ( yieldValueMid <= 0. ) {
1603  sigAnswer = sigMid;
1604  }
1605 
1606  if ( fabs(yieldValueMid) < yieldTol && yieldValueMid <= 0. ) {
1607  // //Put tempKappaP in the status
1608  // status->letTempKappaPBe(tempKappaP);
1609 
1610  ratioPotential =
1611  computeRatioPotential(sigAnswer, tempKappaP);
1612 
1613 
1614  double ratioTrial = rhoTrial / ( sigTrial - sigAnswer );
1615 
1616  if ( ( ( ( ratioPotential >= ratioTrial ) && returnType == RT_Tension ) ) ||
1617  ( ( ratioPotential <= ratioTrial ) && returnType == RT_Compression ) ) {
1618  for ( int i = 0; i < 3; i++ ) {
1619  effectiveStress(i) = sigAnswer;
1620  }
1621 
1622  for ( int i = 3; i < effectiveStress.giveSize(); i++ ) {
1623  effectiveStress(i) = 0.;
1624  }
1626  return tempKappaP;
1627  } else {
1630  return kappaInitial;
1631  }
1632  }
1633  }
1634 
1635  for ( int i = 0; i < 3; i++ ) {
1636  effectiveStress(i) = sigAnswer;
1637  }
1638 
1639  for ( int i = 3; i < effectiveStress.giveSize(); i++ ) {
1640  effectiveStress(i) = 0.;
1641  }
1644  return tempKappaP;
1645 }
1646 
1647 
1648 double
1650  double sigTrial,
1651  double rhoTrial,
1652  double sig)
1653 {
1654  //This function is called, if stress state is in vertex case
1655  double equivalentDeltaPlasticStrain;
1656  rho = 0.;
1657  equivalentDeltaPlasticStrain = sqrt( 1. / 9. * pow( ( sigTrial - sig ) / ( kM ), 2. ) +
1658  pow(rhoTrial / ( 2. * gM ), 2.) );
1659 
1660  double thetaVertex = M_PI / 3.;
1661  double ductilityMeasure = computeDuctilityMeasure(sig, rho, thetaVertex);
1662 
1663  return kappaInitial + equivalentDeltaPlasticStrain / ductilityMeasure;
1664 }
1665 
1666 
1667 double
1669  double rho,
1670  double theta)
1671 {
1672  double thetaConst = pow(2. * cos(theta), 2.);
1673  double ductilityMeasure;
1674  double x = -( sig + fc / 3 ) / fc;
1675  if ( x < 0. ) {
1676  /*Introduce exponential help function which results in a smooth
1677  * transition. */
1678  double EHard = BHard - DHard;
1679  double FHard = ( BHard - DHard ) * CHard / ( AHard - BHard );
1680  ductilityMeasure = ( EHard * exp(x / FHard) + DHard ) / thetaConst;
1681  } else {
1682  ductilityMeasure = ( AHard + ( BHard - AHard ) * exp( -x / ( CHard ) ) ) / thetaConst;
1683  }
1684 
1685  return ductilityMeasure;
1686 }
1687 
1688 
1689 double
1691  double tempKappa)
1692 {
1693  //compute yieldHard and yieldSoft
1694  double yieldHardOne = computeHardeningOne(tempKappa);
1695  double yieldHardTwo = computeHardeningTwo(tempKappa);
1696 
1697  //Compute dilation parameter
1698  double AGParam = this->ft * yieldHardTwo * 3 / this->fc + m / 2;
1699  double BGParam =
1700  yieldHardTwo / 3. * ( 1. + this->ft / this->fc ) /
1701  ( log(AGParam) + log(this->dilationConst + 1.) - log(2 * this->dilationConst - 1.) - log(3. * yieldHardTwo + this->m / 2) );
1702 
1703  double R = ( sig - ft / 3. * yieldHardTwo ) / fc / BGParam;
1704  double mQ = AGParam * exp(R);
1705 
1706  double Bl = sig / fc + rho / ( fc * sqrt(6.) );
1707 
1708  double Al = ( 1. - yieldHardOne ) * pow(Bl, 2.) + sqrt(3. / 2.) * rho / fc;
1709 
1710  double dgdsig = 4. * ( 1. - yieldHardOne ) / fc * Al * Bl + pow(yieldHardOne, 2.) * mQ / fc;
1711 
1712  double dgdrho = Al / ( sqrt(6.) * fc ) * ( 4. * ( 1. - yieldHardOne ) * Bl + 6. ) +
1713 
1714  m *pow(yieldHardOne, 2.) / ( sqrt(6.) * fc );
1715 
1716  return dgdrho / dgdsig * 3. * ( 1. - 2. * nu ) / ( 1. + nu );
1717 }
1718 
1719 
1720 double
1722  double kappaP,
1723  GaussPoint *gp)
1724 {
1725  FloatArray trialStress, deviatoricTrialStress;
1726  FloatArray residuals, residualsNorm, deltaIncrement, dGDInv;
1727  FloatMatrix jacobian;
1728  FloatArray unknowns;
1729 
1730  double yieldValue;
1731  double dKappaDDeltaLambda;
1732  double deltaLambda = 0.;
1733  double normOfResiduals = 0.;
1734  int iterationCount = 0;
1735  bool mode3d = effectiveStress.giveSize() > 1;
1736 
1737  //Define stressVariables
1738  double trialSig, trialRho;
1739 
1740  trialStress = effectiveStress;
1741 
1742  //compute the principal directions of the stress
1743  FloatArray helpStress;
1744  FloatMatrix stressPrincipalDir;
1745 
1746  //compute invariants from stress state
1747  if ( mode3d ) {
1748  StructuralMaterial :: computePrincipalValDir(helpStress, stressPrincipalDir, trialStress, principal_stress);
1749  trialSig = computeDeviatoricVolumetricSplit(deviatoricTrialStress, trialStress);
1750  trialRho = computeSecondCoordinate(deviatoricTrialStress);
1751  } else { //1d case
1752  double angle; // this variable is used only as an input to the function computeTrialCoordinates and is not important and it is already calculated
1753  computeTrialCoordinates(trialStress, trialSig, trialRho, angle);
1754  }
1755 
1756  sig = trialSig;
1757  rho = trialRho;
1758 
1759 
1760  // Starting guess:
1761  double tempKappaP = kappaP;
1762 
1763  //initialise unknowns
1764  if ( mode3d ) {
1765  residuals.resize(4);
1766  unknowns.resize(4);
1767  unknowns.at(1) = trialSig;
1768  unknowns.at(2) = trialRho;
1769  unknowns.at(3) = tempKappaP;
1770  unknowns.at(4) = 0.;
1771  } else { //1D case
1772  residuals.resize(3);
1773  unknowns.resize(3);
1774  unknowns.at(1) = trialSig * 3.; // It is calculated as the volumetric stress in this case sigma/3
1775  unknowns.at(2) = tempKappaP;
1776  unknowns.at(3) = 0.;
1777  }
1778 
1779  // Look at the magnitudes of the residuals. You have to scale the yieldValue down.
1780  yieldValue = computeYieldValue(sig, rho, thetaTrial, tempKappaP);
1781 
1782  //initiate residuals
1783  residuals.zero();
1784  residuals.at( residuals.giveSize() ) = yieldValue; //store in the last element of the array
1785 
1786  normOfResiduals = 1.; //just to get into the loop
1787 
1788  /* N.R. iteration for finding the correct plastic return which is found when the norm of the residuals are equal to zero*/
1789 
1790  while ( normOfResiduals > yieldTol ) {
1791  iterationCount++;
1792  if ( iterationCount == newtonIter ) {
1794  return kappaP;
1795  }
1796 
1797  residualsNorm = residuals;
1798  if ( effectiveStress.giveSize() > 1 ) {
1799  //Normalize residuals. Think about it more.
1800  residualsNorm.at(1) /= this->kM;
1801  residualsNorm.at(2) /= 2. * this->gM;
1802  } else { //1D case
1803  residualsNorm.at(1) /= this->eM;
1804  }
1805 
1806  normOfResiduals = residualsNorm.computeNorm();
1807 
1808  if ( std :: isnan(normOfResiduals) ) {
1810  return kappaP;
1811  }
1812 
1813  if ( normOfResiduals > yieldTol ) {
1814  // Test to run newton iteration using inverse of Jacobian
1815  if ( mode3d ) {
1816  computeJacobian(jacobian, sig, rho, tempKappaP, deltaLambda, gp);
1817 
1818  if ( !jacobian.solveForRhs(residuals, deltaIncrement) ) {
1820  return kappaP;
1821  }
1822 
1823  //compute Unknowns
1824  unknowns.subtract(deltaIncrement);
1825 
1826  if ( unknowns.at(4) <= 0. ) { //Keep deltaLambda greater than zero!
1827  unknowns.at(4) = 0.;
1828  }
1829 
1830  if ( unknowns.at(2) <= 0. ) { //Keep rho greater than zero!
1831  unknowns.at(2) = 0.;
1832  }
1833 
1834  if ( unknowns.at(3) - kappaP <= 0. ) { //Keep deltaKappa greater than zero!
1835  unknowns.at(3) = kappaP;
1836  }
1837 
1838  //compute residuals
1839  sig = unknowns.at(1);
1840  rho = unknowns.at(2);
1841 
1842  tempKappaP = unknowns.at(3);
1843  deltaLambda = unknowns.at(4);
1844 
1845  /* Compute the mVector holding the derivatives of the g function and the hardening function*/
1846  computeDGDInv(dGDInv, sig, rho, tempKappaP);
1847  dKappaDDeltaLambda = computeDKappaDDeltaLambda(sig, rho, tempKappaP);
1848 
1849  residuals.at(1) = sig - trialSig + this->kM *deltaLambda *dGDInv.at(1);
1850  residuals.at(2) = rho - trialRho + ( 2. * this->gM ) * deltaLambda * dGDInv.at(2);
1851  residuals.at(3) = -tempKappaP + kappaP + deltaLambda * dKappaDDeltaLambda;
1852  residuals.at(4) = computeYieldValue(sig, rho, thetaTrial, tempKappaP);
1853  } else {
1854  compute1dJacobian(jacobian, 3. * sig, tempKappaP, deltaLambda, gp);
1855 
1856  if ( !jacobian.solveForRhs(residuals, deltaIncrement) ) {
1858  return kappaP;
1859  }
1860 
1861  //compute Unknowns
1862  unknowns.subtract(deltaIncrement);
1863 
1864  if ( unknowns.at(3) <= 0. ) { //Keep deltaLambda greater equal than zero!
1865  unknowns.at(3) = 0.;
1866  }
1867 
1868  if ( unknowns.at(2) - kappaP <= 0. ) { //Keep deltaKappa greater equal than zero!
1869  unknowns.at(2) = kappaP;
1870  }
1871 
1872  //compute residuals
1873  sig = unknowns.at(1) / 3.;
1874  rho = unknowns.at(1) * sqrt(2. / 3.); //for the 1d case
1875 
1876  tempKappaP = unknowns.at(2);
1877  deltaLambda = unknowns.at(3);
1878 
1879  /* Compute the mVector holding the derivatives of the g function and the hardening function*/
1880  double dginv = computeDGDInv1d(unknowns.at(1), tempKappaP);
1881  dKappaDDeltaLambda = computeDKappaDDeltaLambda1d(unknowns.at(1), tempKappaP);
1882 
1883  residuals.at(1) = 3. * ( sig - trialSig ) + this->eM *deltaLambda * dginv;
1884  residuals.at(2) = -tempKappaP + kappaP + deltaLambda * dKappaDDeltaLambda;
1885  residuals.at(3) = computeYieldValue(sig, rho, thetaTrial, tempKappaP);
1886  }
1887  }
1888  }
1889 
1891  if ( mode3d ) {
1892  FloatArray stressPrincipal(6);
1893  stressPrincipal.zero();
1894 
1895  stressPrincipal(0) = sig + sqrt(2. / 3.) * rho * cos(thetaTrial);
1896  stressPrincipal(1) = sig + sqrt(2. / 3.) * rho * cos(thetaTrial - 2. * M_PI / 3.);
1897  stressPrincipal(2) = sig + sqrt(2. / 3.) * rho * cos(thetaTrial + 2. * M_PI / 3.);
1898  transformStressVectorTo(effectiveStress, stressPrincipalDir, stressPrincipal, 1);
1899  } else {
1900  effectiveStress.at(1) = sig * 3;
1901  }
1902 
1903  return tempKappaP;
1904 }
1905 
1906 void
1908  double totalsigma,
1909  double kappa,
1910  double deltaLambda,
1911  GaussPoint *gp)
1912 {
1913  double dFDInv = computeDFDInv1d(totalsigma, kappa);
1914  double dGDInv = computeDGDInv1d(totalsigma, kappa);
1915  double dDGDDInv = computeDDGDDInv1d(totalsigma, kappa);
1916  double dKappaDDeltaLambda = computeDKappaDDeltaLambda(totalsigma, 1, kappa);
1917  double dFDKappa = computeDFDKappa(totalsigma, 1, kappa, true);
1918  double dDGDInvDKappa = computeDDGDInvDKappa1d(totalsigma, kappa);
1919  double dDKappaDDeltaLambdaDKappa = computeDDKappaDDeltaLambdaDKappa1d(totalsigma, kappa);
1920  double dDKappaDDeltaLambdaDInv = computeDDKappaDDeltaLambdaDInv1d(totalsigma, kappa);
1921 
1922  answer.resize(3, 3);
1923  /* Compute matrix*/
1924  /* 1st row */
1925  answer.at(1, 1) = 1. + this->eM * deltaLambda * dDGDDInv;
1926  answer.at(1, 2) = this->eM * deltaLambda * dDGDInvDKappa;
1927  answer.at(1, 3) = this->eM * dGDInv;
1928  /* 2nd row */
1929  answer.at(2, 1) = deltaLambda * dDKappaDDeltaLambdaDInv;
1930  answer.at(2, 2) = deltaLambda * dDKappaDDeltaLambdaDKappa - 1.;
1931  answer.at(2, 3) = dKappaDDeltaLambda;
1932  /* 3rd row */
1933  answer.at(3, 1) = dFDInv;
1934  answer.at(3, 2) = dFDKappa;
1935  answer.at(3, 3) = 0.;
1936 }
1937 
1938 
1939 void
1941  double sig,
1942  double rho,
1943  double kappa,
1944  double deltaLambda,
1945  GaussPoint *gp)
1946 {
1947  FloatArray dFDInv;
1948  computeDFDInv(dFDInv, sig, rho, kappa);
1949 
1950  FloatArray dGDInv;
1951  computeDGDInv(dGDInv, sig, rho, kappa);
1952 
1953  FloatMatrix dDGDDInv;
1954  computeDDGDDInv(dDGDDInv, sig, rho, kappa);
1955 
1956  double dKappaDDeltaLambda = computeDKappaDDeltaLambda(sig, rho, kappa);
1957 
1958  double dFDKappa = computeDFDKappa(sig, rho, kappa, false);
1959 
1960  FloatArray dDGDInvDKappa;
1961  computeDDGDInvDKappa(dDGDInvDKappa, sig, rho, kappa);
1962 
1963 
1964  double dDKappaDDeltaLambdaDKappa;
1965  dDKappaDDeltaLambdaDKappa = computeDDKappaDDeltaLambdaDKappa(sig, rho, kappa);
1966 
1967 
1968  FloatArray dDKappaDDeltaLambdaDInv;
1969  computeDDKappaDDeltaLambdaDInv(dDKappaDDeltaLambdaDInv, sig, rho, kappa);
1970 
1971 
1972  answer.resize(4, 4);
1973  /* Compute matrix*/
1974  answer.at(1, 1) = 1. + this->kM *deltaLambda *dDGDDInv.at(1, 1);
1975  answer.at(1, 2) = this->kM * deltaLambda * dDGDDInv.at(1, 2);
1976  answer.at(1, 3) = this->kM * deltaLambda * dDGDInvDKappa.at(1);
1977  answer.at(1, 4) = this->kM * dGDInv.at(1);
1978 
1979  answer.at(2, 1) = 2. *this->gM *deltaLambda *dDGDDInv.at(2, 1);
1980  answer.at(2, 2) = 1. + 2. *this->gM *deltaLambda *dDGDDInv.at(2, 2);
1981  answer.at(2, 3) = 2. *this->gM *deltaLambda *dDGDInvDKappa.at(2);
1982  answer.at(2, 4) = 2. *this->gM *dGDInv.at(2);
1983 
1984  answer.at(3, 1) = deltaLambda * dDKappaDDeltaLambdaDInv.at(1);
1985  answer.at(3, 2) = deltaLambda * dDKappaDDeltaLambdaDInv.at(2);
1986  answer.at(3, 3) = deltaLambda * dDKappaDDeltaLambdaDKappa - 1.;
1987  answer.at(3, 4) = dKappaDDeltaLambda;
1988 
1989  answer.at(4, 1) = dFDInv.at(1);
1990  answer.at(4, 2) = dFDInv.at(2);
1991  answer.at(4, 3) = dFDKappa;
1992  answer.at(4, 4) = 0.;
1993 }
1994 
1995 
1996 
1997 
1998 double
2000  double rho,
2001  double theta,
2002  double tempKappa) const
2003 {
2004  //compute yieldHard
2005  double yieldHardOne = computeHardeningOne(tempKappa);
2006  double yieldHardTwo = computeHardeningTwo(tempKappa);
2007 
2008 
2009  // compute elliptic function r
2010  double rFunction = ( 4. * ( 1. - pow(ecc, 2.) ) * pow(cos(theta), 2.) +
2011  pow( ( 2. * ecc - 1. ), 2. ) ) /
2012  ( 2. * ( 1. - pow(ecc, 2.) ) * cos(theta) +
2013  ( 2. * ecc - 1. ) * sqrt(4. * ( 1. - pow(ecc, 2.) ) * pow(cos(theta), 2.)
2014  + 5. * pow(ecc, 2.) - 4. * ecc) );
2015 
2016  //compute help function Al
2017  double Al = ( 1. - yieldHardOne ) * pow( ( sig / fc + rho / ( sqrt(6.) * fc ) ), 2. ) +
2018  sqrt(3. / 2.) * rho / fc;
2019 
2020  //Compute yield equation
2021  return pow(Al, 2.) +
2022  pow(yieldHardOne, 2.) * yieldHardTwo * m * ( sig / fc + rho * rFunction / ( sqrt(6.) * fc ) ) -
2023  pow(yieldHardOne, 2.) * pow(yieldHardTwo, 2.);
2024 }
2025 
2026 
2027 double
2029  double rho,
2030  double tempKappa,
2031  bool mode1d)
2032 {
2033  double dFDKappa;
2034  double theta = thetaTrial;
2035  //compute yieldHard and yieldSoft
2036  double yieldHardOne = computeHardeningOne(tempKappa);
2037  double yieldHardTwo = computeHardeningTwo(tempKappa);
2038  // compute the derivative of the hardening and softening laws
2039  double dYieldHardOneDKappa = computeHardeningOnePrime(tempKappa);
2040  double dYieldHardTwoDKappa = computeHardeningTwoPrime(tempKappa);
2041  //compute elliptic function r
2042  double rFunction =
2043  ( 4. * ( 1. - pow(ecc, 2) ) * pow(cos(theta), 2) + pow( ( 2. * ecc - 1. ), 2 ) ) /
2044  ( 2 * ( 1. - pow(ecc, 2) ) * cos(theta) + ( 2. * ecc - 1. ) * sqrt(4. * ( 1. - pow(ecc, 2) ) * pow(cos(theta), 2) + 5. * pow(ecc, 2) - 4. * ecc) );
2045 
2046 
2047  if ( !mode1d ) {
2048  //compute help functions Al, Bl
2049  double Al = ( 1. - yieldHardOne ) * pow( ( sig / fc + rho / ( sqrt(6.) * fc ) ), 2. ) + sqrt(3. / 2.) * rho / fc;
2050 
2051  double Bl = sig / fc + rho / ( fc * sqrt(6.) );
2052 
2053  double dFDYieldHardOne = -2. *Al *pow(Bl, 2.)
2054  + 2. * yieldHardOne * yieldHardTwo * m * ( sig / fc + rho * rFunction / ( sqrt(6.) * fc ) ) - 2. *yieldHardOne *pow(yieldHardTwo, 2.);
2055 
2056  double dFDYieldHardTwo = pow(yieldHardOne, 2.) * m * ( sig / fc + rho * rFunction / ( sqrt(6.) * fc ) ) - 2. *yieldHardTwo *pow(yieldHardOne, 2.);
2057 
2058  // compute dFDKappa
2059  dFDKappa = dFDYieldHardOne * dYieldHardOneDKappa +
2060  dFDYieldHardTwo * dYieldHardTwoDKappa;
2061  } else { //1d case
2062  dFDKappa = -2 * pow(2 * sig / 3 / fc, 2) * ( sig / fc + pow(2 / 3 * sig / fc, 2) * ( 1 - yieldHardOne ) ) * dYieldHardOneDKappa +
2063  ( 1 + rFunction ) * m * sig / 3 / fc * ( dYieldHardOneDKappa * 2. * yieldHardOne * yieldHardTwo + dYieldHardTwoDKappa * yieldHardOne ) -
2064  2 * ( yieldHardOne * pow(yieldHardTwo, 2) * dYieldHardOneDKappa + yieldHardTwo * pow(yieldHardOne, 2) * dYieldHardTwoDKappa );
2065  }
2066 
2067  /*
2068  * set dFDKappa to zero, if it becomes greater than zero.
2069  * dFDKappa can only be negative or zero in the converged state for
2070  * the case of hardenig and perfect plasticity. For trial stresses, however,
2071  * it might be psoitive, which may lead to convergency problems. Therefore,
2072  * it is set to zero in this cases.
2073  */
2074  if ( dFDKappa > 0. ) {
2075  dFDKappa = 0.;
2076  }
2077 
2078  return dFDKappa;
2079 }
2080 
2081 double
2083  double tempKappa) const
2084 {
2085  double theta = thetaTrial;
2086  //compute yieldHard
2087  double yieldHardOne = computeHardeningOne(tempKappa);
2088  double yieldHardTwo = computeHardeningTwo(tempKappa);
2089  //compute elliptic function r
2090  double rFunction = ( 4. * ( 1. - pow(ecc, 2) ) * pow(cos(theta), 2) + pow( ( 2. * ecc - 1. ), 2 ) ) /
2091  ( 2. * ( 1. - pow(ecc, 2) ) * cos(theta) + ( 2. * ecc - 1. ) * sqrt(4. * ( 1. - pow(ecc, 2) ) * pow(cos(theta), 2) + 5. * pow(ecc, 2) - 4. * ecc) );
2092  return 2 * ( 1 / fc + 8 * sigma / pow(3 * fc, 2) * ( 1 - yieldHardOne ) ) * ( sigma / fc + pow(2 * sigma / 3 / fc, 2) * ( 1 - yieldHardOne ) ) + ( 1 + rFunction ) * m / ( 3 * fc ) * pow(yieldHardOne, 2) * yieldHardTwo;
2093 }
2094 
2095 
2096 void
2098  double sig,
2099  double rho,
2100  double tempKappa) const
2101 {
2102  double theta = thetaTrial;
2103 
2104  //compute yieldHard
2105  double yieldHardOne = computeHardeningOne(tempKappa);
2106  double yieldHardTwo = computeHardeningTwo(tempKappa);
2107 
2108  //compute elliptic function r
2109  double rFunction = ( 4. * ( 1. - ecc * ecc ) * cos(theta) * cos(theta) + ( 2. * ecc - 1. ) * ( 2. * ecc - 1. ) ) /
2110  ( 2. * ( 1. - ecc * ecc ) * cos(theta) + ( 2. * ecc - 1. ) * sqrt(4. * ( 1. - ecc * ecc ) * cos(theta) * cos(theta)
2111  + 5. * ecc * ecc - 4. * ecc) );
2112 
2113  //compute help functions AL, BL
2114  double AL = ( 1. - yieldHardOne ) * pow( ( sig / fc + rho / ( sqrt(6.) * fc ) ), 2. ) + sqrt(3. / 2.) * rho / fc;
2115  double BL = sig / fc + rho / ( fc * sqrt(6.) );
2116 
2117  //compute dfdsig
2118  double dfdsig = 4. * ( 1. - yieldHardOne ) / fc * AL * BL + yieldHardTwo *pow(yieldHardOne, 2.) * m / fc;
2119 
2120  //compute dfdrho
2121  double dfdrho = AL / ( sqrt(6.) * fc ) * ( 4. * ( 1. - yieldHardOne ) * BL + 6. ) + rFunction *m *yieldHardTwo *pow(yieldHardOne, 2.) / ( sqrt(6.) * fc );
2122 
2123  answer.resize(2);
2124  answer(0) = dfdsig;
2125  answer(1) = dfdrho;
2126 }
2127 
2128 
2129 double
2131 {
2132  double equivalentDGDStress = fabs(computeDGDInv1d(sig, tempKappa)); // In 1D is the absolute value
2133  double ductilityMeasure = computeDuctilityMeasure(sig / 3., sqrt(2. / 3.) * sig, this->thetaTrial);
2134  return equivalentDGDStress / ductilityMeasure; // dKappaDDeltaLambda
2135 }
2136 
2137 
2138 double
2140  double rho,
2141  double tempKappa)
2142 {
2143  double equivalentDGDStress;
2144  double ductilityMeasure;
2145 
2146  FloatArray dGDInv;
2147  computeDGDInv(dGDInv, sig, rho, tempKappa);
2148 
2149  equivalentDGDStress = sqrt( 1. / 3. * pow(dGDInv(0), 2.) + pow(dGDInv(1), 2.) );
2150 
2151  ductilityMeasure = computeDuctilityMeasure(sig, rho, this->thetaTrial);
2152 
2153  return equivalentDGDStress / ductilityMeasure; // dKappaDDeltaLambda
2154 }
2155 
2156 double
2158 {
2159  double dGDInv;
2160  double dDGDDInv;
2161 
2162  //Compute first and second derivative of plastic potential
2163  dGDInv = computeDGDInv1d(sigma, tempKappa);
2164  dDGDDInv = computeDDGDDInv1d(sigma, tempKappa);
2165 
2166 
2167  //computeDuctilityMeasure
2168  double ductilityMeasure = computeDuctilityMeasure(sigma / 3., sigma * sqrt(2. / 3.), this->thetaTrial);
2169 
2170  // compute the derivative of
2171  double dDuctilityMeasureDInv = computeDDuctilityMeasureDInv1d(sigma, tempKappa);
2172  if ( dGDInv < 0 ) {
2173  dDGDDInv = -dDGDDInv;
2174  dGDInv = -dGDInv;
2175  }
2176 
2177  return dDGDDInv / ductilityMeasure - dGDInv * dDuctilityMeasureDInv / pow(ductilityMeasure, 2);
2178 }
2179 
2180 
2181 
2182 void
2184  double sig,
2185  double rho,
2186  double tempKappa)
2187 {
2188  double equivalentDGDStress;
2189  FloatArray dGDInv;
2190  FloatMatrix dDGDDInv;
2191  FloatArray dEquivalentDGDStressDInv(2);
2192 
2193  //Compute first and second derivative of plastic potential
2194  computeDGDInv(dGDInv, sig, rho, tempKappa);
2195  computeDDGDDInv(dDGDDInv, sig, rho, tempKappa);
2196 
2197  //Compute equivalentDGDStress
2198  equivalentDGDStress = sqrt( 1. / 3. * pow(dGDInv(0), 2.) + pow(dGDInv(1), 2.) );
2199 
2200  //computeDuctilityMeasure
2201  double ductilityMeasure = computeDuctilityMeasure(sig, rho, this->thetaTrial);
2202 
2203  //Compute dEquivalentDGDStressDInv
2204  dEquivalentDGDStressDInv(0) =
2205  ( 2. / 3. * dGDInv(0) * dDGDDInv(0, 0) + 2. * dGDInv(1) * dDGDDInv(1, 0) ) / ( 2. * equivalentDGDStress );
2206  dEquivalentDGDStressDInv(1) =
2207  ( 2. / 3. * dGDInv(0) * dDGDDInv(0, 1) + 2. * dGDInv(1) * dDGDDInv(1, 1) ) / ( 2. * equivalentDGDStress );
2208 
2209 
2210  // compute the derivative of
2211  FloatArray dDuctilityMeasureDInv(2);
2212  computeDDuctilityMeasureDInv(dDuctilityMeasureDInv, sig, rho, tempKappa);
2213 
2214  answer.resize(2);
2215  answer(0) = ( dEquivalentDGDStressDInv(0) * ductilityMeasure - equivalentDGDStress * dDuctilityMeasureDInv(0) ) / pow(ductilityMeasure, 2.);
2216 
2217  answer(1) = ( dEquivalentDGDStressDInv(1) * ductilityMeasure - equivalentDGDStress * dDuctilityMeasureDInv(1) ) / pow(ductilityMeasure, 2.);
2218 }
2219 
2220 
2221 double
2223 {
2224  double equivalentDGDStress, dEquivalentDGDStressDKappa, ductilityMeasure;
2225 
2226  FloatArray dGDInv, dDGDInvDKappa;
2227 
2228  //Compute first and second derivative of plastic potential
2229  computeDGDInv(dGDInv, sig, rho, tempKappa);
2230  computeDDGDInvDKappa(dDGDInvDKappa, sig, rho, tempKappa);
2231 
2232  equivalentDGDStress = sqrt( 1. / 3. * pow(dGDInv(0), 2.) +
2233  pow(dGDInv(1), 2.) );
2234 
2235  ductilityMeasure = computeDuctilityMeasure(sig, rho, this->thetaTrial); //computeDuctilityMeasure
2236  //Compute dEquivalentDGDStressDKappa
2237  dEquivalentDGDStressDKappa =
2238  ( 2. / 3. * dGDInv(0) * dDGDInvDKappa(0) + 2. * dGDInv(1) * dDGDInvDKappa(1) ) / ( 2. * equivalentDGDStress );
2239 
2240 #if 0
2241  // compute the derivative of
2242  double dDuctilityMeasureDKappa = 0.;
2243 
2245  double dDKappaDDeltaLambdaDKappa =
2246  ( dEquivalentDGDStressDKappa * ductilityMeasure -
2247  equivalentDGDStress * dDuctilityMeasureDKappa ) / pow(ductilityMeasure, 2.);
2248 #endif
2249  return dEquivalentDGDStressDKappa / ductilityMeasure; // dDKappaDDeltaLambdaDKappa
2250 }
2251 
2252 
2253 double
2255  double tempKappa)
2256 {
2257  double equivalentDGDStress, dEquivalentDGDStressDKappa, ductilityMeasure;
2258 
2259  equivalentDGDStress = computeDGDInv1d(sig, tempKappa);
2260  dEquivalentDGDStressDKappa = computeDDGDInvDKappa1d(sig, tempKappa);
2261  if ( equivalentDGDStress < 0 ) {
2262  dEquivalentDGDStressDKappa = ( -1 ) * dEquivalentDGDStressDKappa; //We are differentiating the absolute value of the first derivative of G with respect to stress
2263  }
2264 
2265  ductilityMeasure = computeDuctilityMeasure(sig / 3., sig * sqrt(2. / 3.), this->thetaTrial);
2266 
2267  return dEquivalentDGDStressDKappa / ductilityMeasure; // dDKappaDDeltaLambdaDKappa
2268 }
2269 
2270 
2271 double
2273  double tempKappa)
2274 {
2275  double thetaConst = pow(2. * cos(thetaTrial), 2.);
2276  double x = -( sigma + fc ) / ( 3 * fc ); //R hardening variable
2277  double dXDSigma = -1. / ( 3. * fc );
2278  if ( x < 0. ) {
2279  /* Introduce exponential help function which results in a
2280  * smooth transition. */
2281  double EHard = BHard - DHard;
2282  double FHard = ( BHard - DHard ) * CHard / ( AHard - BHard );
2283 
2284  double dDuctilityMeasureDX = EHard / FHard *exp(x / FHard) / thetaConst;
2285  return dDuctilityMeasureDX * dXDSigma;
2286  } else {
2287  double dDuctilityMeasureDX = ( AHard - BHard ) / CHard / thetaConst *exp(-x / CHard);
2288  return dDuctilityMeasureDX * dXDSigma;
2289  }
2290 }
2291 
2292 void
2294  double sig,
2295  double rho,
2296  double tempKappa)
2297 {
2298  double thetaConst = pow(2. * cos(thetaTrial), 2.);
2299  double x = ( -( sig + fc / 3. ) ) / fc;
2300 
2301  if ( x < 0. ) {
2302  double dXDSig = -1. / fc;
2303  /* Introduce exponential help function which results in a
2304  * smooth transition. */
2305  double EHard = BHard - DHard;
2306  double FHard = ( BHard - DHard ) * CHard / ( AHard - BHard );
2307 
2308  double dDuctilityMeasureDX = EHard / FHard *exp(x / FHard) / thetaConst;
2309  answer = {dDuctilityMeasureDX * dXDSig, 0.};
2310  } else {
2311  double dXDSig = -1. / fc;
2312  double dDuctilityMeasureDX = -( BHard - AHard ) / ( CHard ) / thetaConst *exp( -x / ( CHard ) );
2313  answer = {dDuctilityMeasureDX * dXDSig, 0.};
2314  }
2315 }
2316 
2317 
2318 double
2320  double tempKappa)
2321 {
2322  //compute yieldHard and yieldSoft
2323  double yieldHardOne = computeHardeningOne(tempKappa);
2324  double yieldHardTwo = computeHardeningTwo(tempKappa);
2325  double AGParam = this->ft * yieldHardTwo * 3 / this->fc + m / 2;
2326  double BGParam = yieldHardTwo / 3. * ( 1. + this->ft / this->fc ) / ( log(AGParam) + log(this->dilationConst + 1.) - log(2 * this->dilationConst - 1.) - log(3. * yieldHardTwo + this->m / 2) );
2327  double R = ( sigma - yieldHardTwo * ft ) / ( 3 * fc * BGParam );
2328  double mQ = AGParam * exp(R) / 3;
2329  return 2 * ( 1 / fc + 8 * sigma / pow(3 * fc, 2) * ( 1 - yieldHardOne ) ) * ( sigma / fc + pow(2 * sigma / ( 3 * fc ), 2) * ( 1 - yieldHardOne ) ) + pow(yieldHardOne, 2) / fc * ( m / 3 + mQ );
2330 }
2331 
2332 
2333 void
2335  double sig,
2336  double rho,
2337  double tempKappa)
2338 {
2339  //compute yieldHard and yieldSoft
2340  double yieldHardOne = computeHardeningOne(tempKappa);
2341  double yieldHardTwo = computeHardeningTwo(tempKappa);
2342 
2343  //Compute dilation parameter
2344  double AGParam = this->ft * yieldHardTwo * 3 / this->fc + m / 2;
2345  double BGParam =
2346  yieldHardTwo / 3. * ( 1. + this->ft / this->fc ) /
2347  ( log(AGParam) + log(this->dilationConst + 1.) - log(2 * this->dilationConst - 1.) - log(3. * yieldHardTwo + this->m / 2) );
2348 
2349  double R = ( sig - ft / 3. * yieldHardTwo ) / fc / BGParam;
2350  double mQ = AGParam * exp(R);
2351 
2352  double Bl = sig / fc + rho / ( fc * sqrt(6.) );
2353 
2354  double Al = ( 1. - yieldHardOne ) * pow(Bl, 2.) + sqrt(3. / 2.) * rho / fc;
2355 
2356  double dgdsig = 4. * ( 1. - yieldHardOne ) / fc * Al * Bl + pow(yieldHardOne, 2.) * mQ / fc;
2357 
2358  double dgdrho = Al / ( sqrt(6.) * fc ) * ( 4. * ( 1. - yieldHardOne ) * Bl + 6. ) +
2359  m *pow(yieldHardOne, 2.) / ( sqrt(6.) * fc );
2360 
2361  answer = {dgdsig, dgdrho};
2362 }
2363 
2364 double
2366  double tempKappa)
2367 {
2368  //compute yieldHard and yieldSoft
2369  double yieldHardOne = computeHardeningOne(tempKappa);
2370  double yieldHardTwo = computeHardeningTwo(tempKappa);
2371 
2372  double dYieldHardOneDKappa = computeHardeningOnePrime(tempKappa);
2373  double dYieldHardTwoDKappa = computeHardeningTwoPrime(tempKappa);
2374 
2375  //Compute dilation parameter
2376  double AGParam = this->ft * yieldHardTwo * 3 / this->fc + m / 2;
2377  double BGParam =
2378  yieldHardTwo / 3. * ( 1. + this->ft / this->fc ) /
2379  ( log(AGParam) + log(this->dilationConst + 1.) - log(2 * this->dilationConst - 1.) - log(3. * yieldHardTwo + this->m / 2) );
2380 
2381  double R = ( sigma - ft * yieldHardTwo ) / ( 3 * fc * BGParam );
2382  double mQ = AGParam * exp(R) / 3;
2383  //Compute the derivative of mQ with respect to kappa
2384 
2385  //Derivative of AGParam
2386  double dAGParamDKappa = dYieldHardTwoDKappa * 3. * this->ft / this->fc;
2387 
2388  //Derivative of BGParam
2389  double BGParamTop = yieldHardTwo / 3. * ( 1. + this->ft / this->fc );
2390  double BGParamBottom = ( log(AGParam) + log(this->dilationConst + 1.) - log(2 * this->dilationConst - 1.) - log(3. * yieldHardTwo + this->m / 2) );
2391  double dBGParamTopDKappa1 = dYieldHardTwoDKappa * ( 1 + ft / fc ) / 3;
2392  double dBGParamBottomDKappa1 = BGParamBottom;
2393  double dBGParamTopDKappa2 = BGParamTop * ( dAGParamDKappa / AGParam - 3 * dYieldHardTwoDKappa / ( m / 2 + 3 * yieldHardTwo ) );
2394  double dBGParamBottomDKappa2 = pow(BGParamBottom, 2);
2395 
2396  double dBGParamDKappa = dBGParamTopDKappa1 / dBGParamBottomDKappa1 - dBGParamTopDKappa2 / dBGParamBottomDKappa2;
2397  double dMQDKappa = 1 / 3 * exp(R) * ( dAGParamDKappa - AGParam * ( ( sigma - ft * yieldHardTwo ) * dBGParamDKappa / ( 3 * fc * pow(BGParam, 2) ) + ft * dYieldHardTwoDKappa / 3 / fc / BGParam ) );
2398 
2399  return -8 / 9 * pow(sigma / fc, 2) * ( 1 / fc + 8 * sigma / pow(3 * fc, 2) * ( 1 - yieldHardOne ) ) * dYieldHardOneDKappa - sigma *pow(4 / 3 / fc, 2) * ( sigma / fc + pow(2 / 3 * sigma / fc, 2) * ( 1 - yieldHardOne ) ) * dYieldHardOneDKappa + 2 * dYieldHardOneDKappa * yieldHardOne / fc * ( this->m / 3 + mQ ) + yieldHardOne / fc * dMQDKappa;
2400 }
2401 
2402 
2403 void
2405  double sig,
2406  double rho,
2407  double tempKappa)
2408 {
2409  //Compute dilation parameter
2410 
2411  //compute yieldHard and yieldSoft
2412  double yieldHardOne = computeHardeningOne(tempKappa);
2413  double yieldHardTwo = computeHardeningTwo(tempKappa);
2414 
2415  double dYieldHardOneDKappa = computeHardeningOnePrime(tempKappa);
2416  double dYieldHardTwoDKappa = computeHardeningTwoPrime(tempKappa);
2417 
2418  //Compute dilation parameter
2419  double AGParam = this->ft * yieldHardTwo * 3 / this->fc + m / 2;
2420  double BGParam =
2421  yieldHardTwo / 3. * ( 1. + this->ft / this->fc ) /
2422  ( log(AGParam) + log(this->dilationConst + 1.) - log(2 * this->dilationConst - 1.) - log(3. * yieldHardTwo + this->m / 2) );
2423 
2424  double R = ( sig - ft / 3. * yieldHardTwo ) / ( fc * BGParam );
2425  double mQ = AGParam * exp(R);
2426 
2427  //Compute the derivative of mQ with respect to kappa
2428 
2429  //Derivative of AGParam
2430  double dAGParamDKappa = dYieldHardTwoDKappa * 3. * this->ft / this->fc;
2431 
2432  //Derivative of BGParam
2433  double BGParamTop = yieldHardTwo / 3. * ( 1. + this->ft / this->fc );
2434  double BGParamBottom = ( log(AGParam) + log(this->dilationConst + 1.) - log(2 * this->dilationConst - 1.) - log(3. * yieldHardTwo + this->m / 2) );
2435 
2436  double dBGParamTopDKappa = dYieldHardTwoDKappa / 3. * ( 1. + this->ft / this->fc );
2437  double dBGParamBottomDKappa = 1./AGParam*dAGParamDKappa - 3. * dYieldHardTwoDKappa / ( 3 * yieldHardTwo + m / 2. );
2438  double dBGParamDKappa = ( dBGParamTopDKappa * BGParamBottom - BGParamTop * dBGParamBottomDKappa ) / pow(BGParamBottom, 2.);
2439 
2440  //Derivative of R
2441  double RTop = ( sig - ft / 3. * yieldHardTwo );
2442  double RBottom = fc * BGParam;
2443  double dRTopDKappa = -this->ft / 3. * dYieldHardTwoDKappa;
2444  double dRBottomDKappa = this->fc * dBGParamDKappa;
2445  double dRDKappa = ( dRTopDKappa * RBottom - RTop * dRBottomDKappa ) / pow(RBottom, 2.);
2446 
2447  double dMQDKappa = dAGParamDKappa * exp(R) + AGParam *dRDKappa *exp(R);
2448 
2449  double Bl = sig / fc + rho / ( fc * sqrt(6.) );
2450 
2451  double Al = ( 1. - yieldHardOne ) * pow(Bl, 2.) + sqrt(3. / 2.) * rho / fc;
2452 
2453  double dAlDYieldHard = -pow(Bl, 2.);
2454 
2455  const double dDGDSigDKappa =
2456  ( -4. * Al * Bl / fc + 4. * ( 1 - yieldHardOne ) / fc * dAlDYieldHard * Bl ) * dYieldHardOneDKappa +
2457  dYieldHardOneDKappa * 2 * yieldHardOne * mQ / fc + pow(yieldHardOne,2.) * dMQDKappa / fc;
2458 
2459  const double dDGDRhoDKappa =
2460  ( dAlDYieldHard / ( sqrt(6.) * fc ) * ( 4. * ( 1. - yieldHardOne ) * Bl + 6. ) -
2461  4. * Al / ( sqrt(6.) * fc ) * Bl + m / ( sqrt(6.) * fc ) * 2 * yieldHardOne) * dYieldHardOneDKappa;
2462 
2463  answer = {dDGDSigDKappa, dDGDRhoDKappa};
2464 }
2465 
2466 double
2468  double tempKappa)
2469 {
2470  //compute yieldHardOne and yieldSoft
2471  double yieldHardOne = computeHardeningOne(tempKappa);
2472  double yieldHardTwo = computeHardeningTwo(tempKappa);
2473  double AGParam = this->ft * yieldHardTwo * 3 / this->fc + m / 2;
2474  double BGParam =
2475  yieldHardTwo / 3. * ( 1. + this->ft / this->fc ) /
2476  ( log(AGParam) + log(this->dilationConst + 1.) - log(2 * this->dilationConst - 1.) - log(3. * yieldHardTwo + this->m / 2) );
2477  double R = ( sigma - ft * yieldHardTwo ) / ( 3 * fc * BGParam );
2478  double dMQDSigma = AGParam / ( 9 * BGParam * fc ) * exp(R);
2479  return 2 * pow(1 / fc + 8 * sigma / pow(3 * fc, 2) * ( 1 - yieldHardOne ), 2) + pow(4 / 3 / fc, 2) * ( sigma / fc + pow(2 / 3 * sigma / fc, 2) * ( 1 - yieldHardOne ) ) * ( 1 - yieldHardOne ) + pow(yieldHardOne, 2) / fc * dMQDSigma;
2480 }
2481 
2482 void
2484  double sig,
2485  double rho,
2486  double tempKappa)
2487 {
2488  //compute yieldHardOne and yieldSoft
2489  double yieldHardOne = computeHardeningOne(tempKappa);
2490  double yieldHardTwo = computeHardeningTwo(tempKappa);
2491 
2492  //CoQpute dilation parameter
2493  double AGParam = this->ft * yieldHardTwo * 3 / this->fc + m / 2;
2494  double BGParam =
2495  yieldHardTwo / 3. * ( 1. + this->ft / this->fc ) /
2496  ( log(AGParam) + log(this->dilationConst + 1.) - log(2 * this->dilationConst - 1.) - log(3. * yieldHardTwo + this->m / 2) );
2497 
2498  double R = ( sig - ft / 3. * yieldHardTwo ) / fc / BGParam;
2499 
2500  double dMQDSig = AGParam / ( BGParam * fc ) * exp(R);
2501 
2502  //compute help parameter Al and Bl and the corresponding derivatives
2503  double Bl = sig / fc + rho / ( fc * sqrt(6.) );
2504 
2505  double Al = ( 1. - yieldHardOne ) * pow(Bl, 2.) +
2506  sqrt(3. / 2.) * rho / fc;
2507 
2508  double dAlDSig = 2. * ( 1. - yieldHardOne ) * Bl / fc;
2509  double dBlDSig = 1. / fc;
2510 
2511  double dAlDRho = 2. * ( 1. - yieldHardOne ) * Bl / ( fc * sqrt(6.) ) + sqrt(3. / 2.) / fc;
2512  double dBlDRho = 1. / ( fc * sqrt(6.) );
2513 
2514  //compute second derivatives of g
2515  double ddgddSig = 4. * ( 1. - yieldHardOne ) / fc * ( dAlDSig * Bl + Al * dBlDSig ) +
2516  pow(yieldHardOne, 2.) * dMQDSig / fc;
2517 
2518  double ddgddRho = dAlDRho / ( sqrt(6.) * fc ) * ( 4. * ( 1. - yieldHardOne ) * Bl + 6. ) +
2519  Al * dBlDRho * 4. * ( 1. - yieldHardOne ) / ( sqrt(6.) * fc );
2520 
2521  double ddgdSigdRho = 4. * ( 1. - yieldHardOne ) / fc * ( dAlDRho * Bl + Al * dBlDRho );
2522 
2523  double ddgdRhodSig = dAlDSig / ( sqrt(6.) * fc ) * ( 4. * ( 1. - yieldHardOne ) * Bl + 6. )
2524  + Al / ( sqrt(6.) * fc ) * ( 4. * ( 1. - yieldHardOne ) * dBlDSig );
2525 
2526  answer.resize(2, 2);
2527  answer(0, 0) = ddgddSig;
2528  answer(0, 1) = ddgdSigdRho;
2529  answer(1, 0) = ddgdRhodSig;
2530  answer(1, 1) = ddgddRho;
2531 }
2532 
2533 
2534 
2535 double
2537  FloatArray &effectiveStressCompression,
2538  FloatArray &effectiveStress)
2539 {
2540  FloatMatrix stressPrincipalDir;
2541  FloatArray principalStress;
2542  StructuralMaterial :: computePrincipalValDir(principalStress, stressPrincipalDir, effectiveStress, principal_stress);
2543 
2544  //Split the principal values in a tension and a compression part
2545  FloatArray principalStressTension(6);
2546  FloatArray principalStressCompression(6);
2547 
2548  for ( int i = 1; i <= principalStress.giveSize(); i++ ) {
2549  if ( principalStress.at(i) >= 0 ) {
2550  principalStressTension.at(i) = principalStress.at(i);
2551  principalStressCompression.at(i) = 0.;
2552  } else {
2553  principalStressTension.at(i) = 0.;
2554  principalStressCompression.at(i) = principalStress.at(i);
2555  }
2556  }
2557 
2558  //Transform the tension and compression principal stresses back to the original coordinate system
2559 
2560  //Take care of type of stress state for tension
2561  transformStressVectorTo(effectiveStressTension, stressPrincipalDir, principalStressTension, 1);
2562 
2563  //Take care of type of stress state for compression
2564  transformStressVectorTo(effectiveStressCompression, stressPrincipalDir, principalStressCompression, 1);
2565 
2566  //Determine the two factors from the stress
2567 
2568  double squareNormOfPrincipalStress = 0.;
2569  for ( int i = 1; i <= principalStress.giveSize(); i++ ) {
2570  squareNormOfPrincipalStress += pow(principalStress.at(i), 2.);
2571  }
2572 
2573  double alphaTension = 0.;
2574 
2575  if ( squareNormOfPrincipalStress > 0 ) {
2576  for ( int i = 1; i <= principalStress.giveSize(); i++ ) {
2577  alphaTension += principalStressTension.at(i) *
2578  ( principalStressTension.at(i) + principalStressCompression.at(i) ) / squareNormOfPrincipalStress;
2579  }
2580  }
2581 
2582  return 1. - alphaTension;
2583 }
2584 
2585 
2586 double
2588 {
2589  if ( kappa <= 0. ) {
2590  return yieldHardInitial;
2591  } else if ( kappa > 0. && kappa < 1. ) {
2592  return
2593  ( 1. - yieldHardInitial - yieldHardPrimePeak ) * pow(kappa, 3.)
2594  - ( 3. * ( 1. - yieldHardInitial ) - 3. * yieldHardPrimePeak ) * pow(kappa, 2.)
2595  + ( 3. * ( 1. - yieldHardInitial ) - 2. * yieldHardPrimePeak ) * kappa
2596  + yieldHardInitial;
2597  } else {
2598  return 1.;
2599  }
2600 }
2601 
2602 
2603 double
2605 {
2606  if ( kappa <= 0. ) {
2607  return 3. * ( 1 - yieldHardInitial ) - 2. * yieldHardPrimePeak;
2608  } else if ( kappa >= 0. && kappa < 1. ) {
2609  return
2610  3. * ( 1. - yieldHardInitial - yieldHardPrimePeak ) * pow(kappa, 2.)
2611  - 2. * ( 3. * ( 1. - yieldHardInitial ) - 3. * yieldHardPrimePeak ) * kappa
2612  + ( 3. * ( 1. - yieldHardInitial ) - 2. * yieldHardPrimePeak );
2613  } else {
2614  return 0.;
2615  }
2616 }
2617 
2618 
2619 double
2621 {
2622  if ( kappa <= 0. ) {
2623  return 1.;
2624  } else if ( kappa > 0. && kappa < 1. ) {
2625  return 1.;
2626  } else {
2627  return 1. + ( kappa - 1. ) * yieldHardPrimePeak;
2628  }
2629 }
2630 
2631 
2632 double
2634 {
2635  if ( kappa <= 0. ) {
2636  return 0.;
2637  } else if ( kappa >= 0. && kappa < 1. ) {
2638  return 0.;
2639  } else {
2640  return yieldHardPrimePeak;
2641  }
2642 }
2643 
2644 void
2646 {
2647  ConcreteDPM2Status *status = static_cast< ConcreteDPM2Status * >( this->giveStatus(gp) );
2648  if ( mode == ElasticStiffness ) {
2649  answer.resize(1, 1);
2650  answer.at(1, 1) = eM;
2651  } else if ( mode == SecantStiffness ) {
2652  double om;
2653  const FloatArray &strain = status->giveTempStrainVector();
2654  if ( strain.at(1) >= 0 || isotropicFlag == 1 ) {
2655  om = status->giveTempDamageTension();
2656  } else {
2657  om = status->giveTempDamageCompression();
2658  }
2659 
2660  if ( om > 0.999999 ) {
2661  om = 0.999999;
2662  }
2663 
2664  answer.resize(1, 1);
2665  answer.at(1, 1) = eM;
2666  answer.times(1.0 - om);
2667  } else {
2668  OOFEM_WARNING("unknown type of stiffness (tangent stiffness not implemented for 1d). Elastic stiffness used!\n");
2669  answer.resize(1, 1);
2670  answer.at(1, 1) = eM;
2671  }
2672 }
2673 
2674 void
2676  MatResponseMode mode,
2677  GaussPoint *gp,
2678  TimeStep *tStep)
2679 {
2680  if ( mode == ElasticStiffness ) {
2681  this->giveLinearElasticMaterial()->give3dMaterialStiffnessMatrix(answer, mode, gp, tStep);
2682  } else if ( mode == SecantStiffness ) {
2683  this->compute3dSecantStiffness(answer, gp, tStep);
2684  } else if ( mode == TangentStiffness ) {
2685  OOFEM_WARNING("unknown type of stiffness (tangent stiffness not implemented). Elastic stiffness used!");
2686  this->giveLinearElasticMaterial()->give3dMaterialStiffnessMatrix(answer, mode, gp, tStep);
2687  }
2688 }
2689 
2690 void
2692  GaussPoint *gp,
2693  TimeStep *tStep)
2694 {
2695  ConcreteDPM2Status *status = static_cast< ConcreteDPM2Status * >( giveStatus(gp) );
2696 
2697  // Damage parameters
2698  double omegaTension = min(status->giveTempDamageTension(), 0.999999);
2699 
2700  this->giveLinearElasticMaterial()->give3dMaterialStiffnessMatrix(answer, ElasticStiffness, gp, tStep);
2701 
2702  if ( isotropicFlag == 1 ) {
2703  answer.times(1. - omegaTension);
2704  return;
2705  }
2706 
2707 
2708  FloatArray effectiveStress;
2709  FloatArray elasticStrain = status->giveTempStrainVector();
2710  elasticStrain.subtract( status->giveTempPlasticStrain() );
2711  effectiveStress.beProductOf(answer, elasticStrain);
2712 
2713  //Calculate the principal values of the effective stress
2714  FloatArray principalStress;
2715  StructuralMaterial :: computePrincipalValues(principalStress, effectiveStress, principal_stress);
2716 
2717  //exclude two special cases.
2718  if ( principalStress.containsOnlyZeroes() ) {
2719  return;
2720  }
2721 
2722  for ( int i = 1; i <= 3; i++ ) {
2723  if ( principalStress.at(i) < -DYNCON_TOL ) {
2724  return;
2725  }
2726  }
2727 
2728  answer.times(1. - omegaTension);
2729 }
2730 
2731 
2732 
2733 void
2734 ConcreteDPM2 :: computeTrialCoordinates(const FloatArray &stress, double &sigNew, double &rhoNew, double &thetaNew)
2735 {
2736  if ( stress.giveSize() == 1 ) { //1d case
2737  sigNew = stress[0] / 3.;
2738  rhoNew = stress[0] * sqrt(2. / 3.);
2739  if ( sigNew >= 0 ) {
2740  thetaNew = 0.;
2741  } else {
2742  thetaNew = M_PI / 6;
2743  }
2744  } else {
2745  FloatArray deviatoricStress;
2746  sigNew = computeDeviatoricVolumetricSplit(deviatoricStress, stress);
2747  rhoNew = computeSecondCoordinate(deviatoricStress);
2748  thetaNew = computeThirdCoordinate(deviatoricStress);
2749  }
2750 }
2751 
2752 
2753 void
2755 {
2756  ConcreteDPM2Status *status = static_cast< ConcreteDPM2Status * >( this->giveStatus(gp) );
2757 
2758  //Get kappaD from status to define state later on
2759  double damageTension = status->giveDamageTension();
2760  double damageCompression = status->giveDamageCompression();
2761  double tempDamageTension = status->giveTempDamageTension();
2762  double tempDamageCompression = status->giveTempDamageCompression();
2763  double kappaP = status->giveKappaP();
2764  double tempKappaP = status->giveTempKappaP();
2765 
2766  if ( tempKappaP > kappaP ) {
2767  if ( tempDamageTension > damageTension || tempDamageTension == 1. ||
2768  tempDamageCompression > damageCompression || tempDamageCompression == 1. ) {
2770  } else {
2772  }
2773  } else {
2774  const int state_flag = status->giveStateFlag();
2775  if ( tempDamageTension > damageTension || tempDamageTension == 1. ||
2776  tempDamageCompression > damageCompression || tempDamageCompression == 1. ) {
2778  } else {
2779  if ( state_flag == ConcreteDPM2Status :: ConcreteDPM2_Elastic ) {
2781  } else {
2783  }
2784  }
2785  }
2786 }
2787 
2788 void
2790  const FloatArray &stress) const
2791 {
2792  int size = 6;
2793  //compute volumetric deviatoric split
2794  FloatArray deviatoricStress;
2795  computeDeviatoricVolumetricSplit(deviatoricStress, stress);
2796  double rho = computeSecondCoordinate(deviatoricStress);
2797 
2798  //compute the derivative of J2 with respect to the stress
2799  FloatArray dJ2DStress;
2800  dJ2DStress = deviatoricStress;
2801  for ( int i = 3; i < size; i++ ) {
2802  dJ2DStress(i) = deviatoricStress(i) * 2.0;
2803  }
2804 
2805  //compute the derivative of rho with respect to stress
2806  answer = dJ2DStress;
2807  answer.times(1. / rho);
2808 }
2809 
2810 void
2812 {
2813  int size = 6;
2814  for ( int i = 0; i < 3; i++ ) {
2815  answer(i) = 1. / 3.;
2816  }
2817 
2818  for ( int i = 3; i < size; i++ ) {
2819  answer(i) = 0.;
2820  }
2821 }
2822 
2823 
2824 void
2826  const FloatArray &stress) const
2827 
2828 {
2829  int size = 6;
2830 
2831  //compute volumetric deviatoric split
2832  FloatArray deviatoricStress;
2833  computeDeviatoricVolumetricSplit(deviatoricStress, stress);
2834  double rho = computeSecondCoordinate(deviatoricStress);
2835 
2836 
2837  //compute first dericative of J2
2838  FloatArray dJ2dstress;
2839  dJ2dstress = deviatoricStress;
2840  for ( int i = 3; i < deviatoricStress.giveSize(); i++ ) {
2841  dJ2dstress(i) = deviatoricStress(i) * 2.;
2842  }
2843 
2844  //compute second derivative of J2
2845  FloatMatrix ddJ2ddstress(size, size);
2846  ddJ2ddstress.zero();
2847  for ( int i = 0; i < size; i++ ) {
2848  if ( i < 3 ) {
2849  ddJ2ddstress(i, i) = 2. / 3.;
2850  }
2851 
2852  if ( i > 2 ) {
2853  ddJ2ddstress(i, i) = 2.;
2854  }
2855  }
2856 
2857  ddJ2ddstress(0, 1) = -1. / 3.;
2858  ddJ2ddstress(0, 2) = -1. / 3.;
2859  ddJ2ddstress(1, 0) = -1. / 3.;
2860  ddJ2ddstress(1, 2) = -1. / 3.;
2861  ddJ2ddstress(2, 0) = -1. / 3.;
2862  ddJ2ddstress(2, 1) = -1. / 3.;
2863 
2864  //compute the second derivative of rho
2865  answer = ddJ2ddstress;
2866  answer.times(1. / rho);
2867  //compute square of the first derivative of J2
2868  answer.plusDyadUnsym(dJ2dstress, dJ2dstress, -1. / ( rho * rho * rho ));
2869 }
2870 
2871 int
2873  GaussPoint *gp,
2874  InternalStateType type,
2875  TimeStep *tStep)
2876 {
2877  ConcreteDPM2Status *status = static_cast< ConcreteDPM2Status * >( this->giveStatus(gp) );
2878 
2879  switch ( type ) {
2880  case IST_PlasticStrainTensor:
2881  answer = status->givePlasticStrain();
2882  return 1;
2883 
2884  case IST_DamageTensor:
2885  answer.resize(6);
2886  answer.zero();
2887  answer.at(1) = status->giveDamageTension();
2888  answer.at(2) = status->giveDamageCompression();
2889  answer.at(3) = status->giveDissWork();
2890  return 1;
2891 
2892 #ifdef keep_track_of_dissipated_energy
2893  case IST_StressWorkDensity:
2894  answer.resize(1);
2895  answer.at(1) = status->giveStressWork();
2896  return 1;
2897 
2898  case IST_DissWorkDensity:
2899  answer.resize(1);
2900  answer.at(1) = status->giveDissWork();
2901  return 1;
2902 
2903  case IST_FreeEnergyDensity:
2904  answer.resize(1);
2905  answer.at(1) = status->giveStressWork() - status->giveDissWork();
2906  return 1;
2907 
2908 #endif
2909  default:
2910  return StructuralMaterial :: giveIPValue(answer, gp, type, tStep);
2911  }
2912 }
2913 
2916 {
2918 }
2919 } //end of namespace
void computeDDRhoDDStress(FloatMatrix &answer, const FloatArray &stress) const
Computes the seconfd derivative of rho with the respect to the stress.
virtual void updateYourself(TimeStep *tStep)
Update equilibrium history variables according to temp-variables.
Definition: structuralms.C:96
#define _IFT_ConcreteDPM2_softeningType
Definition: concretedpm2.h:68
InternalStateType
Type representing the physical meaning of element or constitutive model internal variable.
#define _IFT_ConcreteDPM2_isoflag
Definition: concretedpm2.h:74
static void computePrincipalValues(FloatArray &answer, const FloatArray &s, stressStrainPrincMode mode)
Common functions for convenience.
void subtract(const FloatArray &src)
Subtracts array src to receiver.
Definition: floatarray.C:258
double mQ
Dilation parameter of the plastic potential.
Definition: concretedpm2.h:646
double kM
Elastic bulk modulus.
Definition: concretedpm2.h:659
#define _IFT_ConcreteDPM2_kinit
Definition: concretedpm2.h:56
double timeFactor
Input parameter which simulates a loading rate. Only for debugging purposes.
Definition: concretedpm2.h:688
double giveKappaDCompressionTwo() const
Get the compression hardening variable two of the damage model from the material status.
Definition: concretedpm2.h:258
void letTempPlasticStrainBe(const FloatArray &v)
Assign the temp value of deviatoric plastic strain.
Definition: concretedpm2.h:427
void letTempStrainVectorBe(const FloatArray &v)
Assigns tempStrainVector to given vector v.
Definition: structuralms.h:137
virtual void printOutputAt(FILE *file, TimeStep *tStep)
Print receiver&#39;s output to given stream.
Definition: concretedpm2.C:178
void computeDSigDStress(FloatArray &answer) const
Computes the derivative of sig with respect to the stress.
virtual void give1dStressStiffMtrx(FloatMatrix &answer, MatResponseMode mmode, GaussPoint *gp, TimeStep *tStep)
Method for computing 1d stiffness matrix of receiver.
void letTempKappaDCompressionOneBe(double v)
Assign the temp value of the hardening variable of the damage model.
Definition: concretedpm2.h:469
double computeDGDInv1d(double sig, double tempKappa)
Compute derivative the palstic potential function with respect to the stress state.
void letTempKappaDTensionBe(double v)
Assign the temp value of the rate factor of the damage model.
Definition: concretedpm2.h:448
double giveKappaDCompression() const
Get the temp value of the hardening variable of the damage model from the material status...
Definition: concretedpm2.h:382
double giveKappaDTension() const
Get the temp value of the hardening variable of the damage model from the material status...
Definition: concretedpm2.h:370
void computeDGDInv(FloatArray &answer, double sig, double rho, double tempKappa)
Compute derivative the palstic potential function with respect to the stress state.
virtual MaterialStatus * giveStatus(GaussPoint *gp) const
Returns material status of receiver in given integration point.
Definition: material.C:244
virtual double computeDamageParamTension(double equivStrain, double kappaOne, double kappaTwo, double le, double omegaOld)
Compute damage parameter.
MaterialStatus * CreateStatus(GaussPoint *gp) const
Creates new copy of associated status and inserts it into given integration point.
#define _IFT_ConcreteDPM2_dhard
Definition: concretedpm2.h:60
GaussPoint * gp
Associated integration point.
Class and object Domain.
Definition: domain.h:115
#define _IFT_ConcreteDPM2_helem
Definition: concretedpm2.h:73
double fcZero
This parameter is needed for the rate dependence. It should be read in if rate dependence is consider...
Definition: concretedpm2.h:691
virtual double computeLength()
Computes the length (zero for all but 1D geometries)
Definition: element.C:1129
double computeRatioPotential(double sig, double tempKappa)
This function computes the ratio of the volumetric and deviatoric component of the flow direction...
double computeRateFactor(double alpha, double timeFactor, GaussPoint *gp, TimeStep *deltaTime)
This function computes the rate factor which is used to take into account the strain rate dependence ...
For computing principal stresses.
For computing principal strains from engineering strains.
double giveKappaDCompressionOne() const
Get the compression hardening variable one of the damage model from the material status.
Definition: concretedpm2.h:240
virtual ~ConcreteDPM2Status()
Destructor.
Definition: concretedpm2.C:94
double computeMeanSize()
Computes the size of the element defined as its length.
Definition: element.C:1078
double performRegularReturn(FloatArray &stress, double kappaP, GaussPoint *gp)
Perform regular stress return for the plasticity model, i.e.
bool solveForRhs(const FloatArray &b, FloatArray &answer, bool transpose=false)
Solves the system of linear equations .
Definition: floatmatrix.C:1112
#define _IFT_ConcreteDPM2_chard
Definition: concretedpm2.h:59
void computeJacobian(FloatMatrix &answer, double sig, double rho, double tempKappa, double deltaLambda, GaussPoint *gp)
Compute jacobian for 2D(plane strain) and 3d cases.
const FloatArray & givePlasticStrain() const
Get the plastic strain deviator from the material status.
Definition: concretedpm2.h:196
#define tAlpha
Definition: matconst.h:67
double m
Friction parameter of the yield surface.
Definition: concretedpm2.h:643
static double computeSecondCoordinate(const FloatArray &s)
FloatArray tempStrainVector
Temporary strain vector in reduced form (to find balanced state)
Definition: structuralms.h:75
#define _IFT_ConcreteDPM2_asoft
Definition: concretedpm2.h:62
#define _IFT_ConcreteDPM2_dilation
Definition: concretedpm2.h:61
double eM
Elastic Young&#39;s modulus.
Definition: concretedpm2.h:655
void letTempKappaDCompressionBe(double v)
Assign the temp value of the rate factor of the damage model.
Definition: concretedpm2.h:455
int newtonIter
Maximum number of iterations for stress return.
Definition: concretedpm2.h:682
The purpose of DataStream abstract class is to allow to store/restore context to different streams...
Definition: datastream.h:54
double giveEquivStrainCompression() const
Get the compression equivalent strain from the material status.
Definition: concretedpm2.h:284
double ASoft
Parameter of the ductilityMeasure of the damage model.
Definition: concretedpm2.h:619
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.
double & at(int i)
Coefficient access function.
Definition: floatarray.h:131
int giveGlobalNumber() const
Definition: element.h:1059
int checkForUnAndReloading(double &tempEquivStrain, double &minEquivStrain, const FloatMatrix &D, GaussPoint *gp)
Check for un- and reloading in the damage part.
Definition: concretedpm2.C:961
void letTempAlphaBe(double v)
Definition: concretedpm2.h:329
This class implements a structural material status information.
Definition: structuralms.h:65
void computeWork(GaussPoint *gp, double ft)
Computes the increment of total stress work and of dissipated work (gf is the dissipation density per...
Definition: concretedpm2.C:457
LinearElasticMaterial * giveLinearElasticMaterial()
Definition: concretedpm2.h:707
Updated Lagrange.
Definition: fmode.h:45
double performVertexReturn(FloatArray &stress, double apexStress, double tempKappaP, GaussPoint *gp)
Perform stress return for vertex case of the plasticity model, i.e.
FloatArray stressVector
Equilibrated stress vector in reduced form.
Definition: structuralms.h:71
bool isnan(double x)
Definition: mathfem.h:103
Dictionary propertyDictionary
Property dictionary.
Definition: material.h:105
void setTempStressWork(double w)
Sets the density of total work of stress on strain increments to given value.
Definition: concretedpm2.h:557
void compute3dSecantStiffness(FloatMatrix &answer, GaussPoint *gp, TimeStep *tStep)
LinearElasticMaterial * linearElasticMaterial
Pointer for linear elastic material.
Definition: concretedpm2.h:652
virtual IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description stored in input record.
double yieldHardInitial
Parameter of the hardening law of the plasticity model.
Definition: concretedpm2.h:625
#define _IFT_ConcreteDPM2_ecc
Definition: concretedpm2.h:55
Element * giveElement()
Returns corresponding element to receiver.
Definition: gausspoint.h:188
General IO error.
void computeDDKappaDDeltaLambdaDInv(FloatArray &answer, double sig, double rho, double tempKappa)
Computes the mixed derivative of the hardening parameter kappa with respect to the plastic multiplier...
void letTempDamageCompressionBe(double v)
Assign the temp value of the compressive damage variable of the damage model.
Definition: concretedpm2.h:497
void letTempKappaDTensionTwoBe(double v)
Assign the temp value of the second tension hardening variable of the damage model.
Definition: concretedpm2.h:476
void plusDyadUnsym(const FloatArray &a, const FloatArray &b, double dV)
Adds to the receiver the product .
Definition: floatmatrix.C:811
virtual void updateYourself(TimeStep *tStep)
Update equilibrium history variables according to temp-variables.
Definition: concretedpm2.C:135
void computeDDGDDInv(FloatMatrix &answer, double sig, double rho, double tempKappa)
Here, the second derivative of the plastic potential with respect to the invariants sig and rho are c...
double giveTempDamageCompression() const
Get the temp value of the hardening variable of the damage model from the material status...
Definition: concretedpm2.h:400
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.
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.
double helem
Element size (to be used in fracture energy approach (crack band).
Definition: concretedpm2.h:649
double giveTempKappaP() const
Get the temp value of the hardening variable of the plasticity model from the material status...
Definition: concretedpm2.h:361
#define THROW_CIOERR(e)
Definition: contextioerr.h:61
virtual contextIOResultType saveContext(DataStream &, ContextMode mode, void *obj=NULL)
Stores receiver state to output stream.
Definition: concretedpm2.C:251
double computeDDuctilityMeasureDInv1d(double sigma, double tempKappa)
Compute derivative the ductility measure with respect to the stress state.
virtual double computeDuctilityMeasure(double sig, double rho, double theta)
Compute the ductility measure based on the stress state.
void computeTrialCoordinates(const FloatArray &stress, double &sig, double &rho, double &theta)
Compute the trial coordinates.
void letTempStateFlagBe(const int v)
Assign the temp value of the state flag.
Definition: concretedpm2.h:545
double AHard
Parameter of the ductilityMeasure of the plasticity model.
Definition: concretedpm2.h:607
double giveDissWork()
Returns the density of dissipated work.
Definition: concretedpm2.h:559
double giveAlpha() const
Definition: concretedpm2.h:373
void letTempRateStrainBe(double v)
Definition: concretedpm2.h:326
void beDifferenceOf(const FloatArray &a, const FloatArray &b)
Sets receiver to be a - b.
Definition: floatarray.C:341
void computeDFDInv(FloatArray &answer, double sig, double rho, double tempKappa) const
Computes the derivative of the yield surface with respect to the invariants sig and rho...
double giveTimeIncrement()
Returns solution step associated time increment.
Definition: timestep.h:150
double giveRateFactor() const
Get the rate factor of the damage model from the material status.
Definition: concretedpm2.h:311
bool isTheFirstStep()
Check if receiver is first step.
Definition: timestep.C:134
double tempStressWork
Non-equilibrated density of total work done by stresses on strain increments.
Definition: concretedpm2.h:169
#define _IFT_ConcreteDPM2_newtoniter
Definition: concretedpm2.h:65
double computeDuctilityMeasureDamage(const FloatArray &strain, GaussPoint *gp)
Compute the ductility measure for the damage model.
virtual void giveRealStressVector_3d(FloatArray &answer, GaussPoint *gp, const FloatArray &strainVector, TimeStep *tStep)
Default implementation relies on giveRealStressVector for second Piola-Kirchoff stress.
Definition: concretedpm2.C:719
const FloatArray & giveTempPlasticStrain() const
Get the temp value of the full plastic strain vector from the material status.
Definition: concretedpm2.h:347
virtual int write(const int *data, int count)=0
Writes count integer values from array pointed by data.
double computeDDKappaDDeltaLambdaDInv1d(double sigma, double tempKappa)
void letTempDamageTensionBe(double v)
Assign the temp value of the tensile damage variable of the damage model.
Definition: concretedpm2.h:490
double giveEquivStrain() const
Get the equivalent strain from the material status.
Definition: concretedpm2.h:267
virtual IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description stored in input record.
double computeDKappaDDeltaLambda(double sig, double rho, double tempKappa)
Compute the derivative of kappa with respect of delta lambda based on the stress state and the harden...
#define _IFT_IsotropicLinearElasticMaterial_talpha
#define _IFT_IsotropicLinearElasticMaterial_n
double yieldHardPrimePeak
Parameter of the hardening law of the plasticity model.
Definition: concretedpm2.h:622
#define OOFEM_LOG_INFO(...)
Definition: logger.h:127
#define M_PI
Definition: mathfem.h:52
virtual void initTempStatus()
Initializes the temporary internal variables, describing the current state according to previously re...
Definition: concretedpm2.C:98
double dotProduct(const FloatArray &x) const
Computes the dot product (or inner product) of receiver and argument.
Definition: floatarray.C:463
void assignStateFlag(GaussPoint *gp)
Assign state flag.
double wfOne
Control parameter for the bilinear softening law in tension.
Definition: concretedpm2.h:670
This class implements an isotropic linear elastic material in a finite element problem.
#define _IFT_ConcreteDPM2_timeFactor
Definition: concretedpm2.h:72
#define OOFEM_ERROR(...)
Definition: error.h:61
double thetaTrial
Lode angle of the trial stress..
Definition: concretedpm2.h:639
void initDamaged(double kappa, const FloatArray &strain, GaussPoint *gp)
Initialize the characteristic length, if damage is not yet activated Set the increase factor for the ...
double tempDissWork
Non-equilibrated density of dissipated work.
Definition: concretedpm2.h:173
double nu
Elastic poisson&#39;s ration.
Definition: concretedpm2.h:661
double BHard
Parameter of the ductilityMeasure of the plasticity model.
Definition: concretedpm2.h:609
#define _IFT_ConcreteDPM2_rateFlag
Definition: concretedpm2.h:71
virtual int giveIPValue(FloatArray &answer, GaussPoint *gp, InternalStateType type, TimeStep *tStep)
Returns the integration point corresponding value in Reduced form.
void computeDDuctilityMeasureDInv(FloatArray &answer, double sig, double rho, double tempKappa)
Compute derivative the ductility measure with respect to the stress state.
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 ...
virtual void giveRealStressVector_1d(FloatArray &answer, GaussPoint *gp, const FloatArray &totalStrain, TimeStep *tStep)
Default implementation relies on giveRealStressVector_StressControl.
Definition: concretedpm2.C:639
double dilationConst
Control parameter for te volumetric plastic flow of the plastic potential.
Definition: concretedpm2.h:628
double computeDFDInv1d(double sigma, double tempKappa) const
void times(double f)
Multiplies receiver by factor f.
Definition: floatmatrix.C:1594
void setLe(double ls)
Sets the characteristic length.
Definition: concretedpm2.h:538
void computeDDGDInvDKappa(FloatArray &answer, double sig, double rho, double tempKappa)
Here, the mixed derivative of the plastic potential with respect to the invariants and the hardening ...
double computeDeltaPlasticStrainNormCompression(double tempAlpha, double tempKappaD, double kappaD, GaussPoint *gp)
double efCompression
Control parameter for the exponential softening law.
Definition: concretedpm2.h:664
contextIOResultType restoreYourself(DataStream &stream)
Definition: floatarray.C:895
#define _IFT_ConcreteDPM2_wf
Definition: concretedpm2.h:66
double giveRateStrain() const
Definition: concretedpm2.h:323
#define _IFT_ConcreteDPM2_hp
Definition: concretedpm2.h:63
double computeDFDKappa(double sig, double rho, double tempKappa, bool mode1d)
Compute the derivative of the yield surface with respect to the hardening variable based on the stres...
void beProductOf(const FloatMatrix &aMatrix, const FloatArray &anArray)
Receiver becomes the result of the product of aMatrix and anArray.
Definition: floatarray.C:676
#define _IFT_IsotropicLinearElasticMaterial_e
void letTempRateFactorBe(double v)
Assign the temp value of the rate factor of the damage model.
Definition: concretedpm2.h:504
void computeDamage(FloatArray &answer, const FloatArray &strain, const FloatMatrix &D, double timeFactor, GaussPoint *gp, TimeStep *tStep, double alpha)
Compute damage parameters.
Definition: concretedpm2.C:794
bool checkForVertexCase(double &answer, double sig, double tempKappa, bool mode1d)
Check if the trial stress state falls within the vertex region of the plasticity model at the apex of...
#define DYNCON_TOL
Definition: concretedpm2.h:47
const FloatArray & giveTempStrainVector() const
Returns the const pointer to receiver&#39;s temporary strain vector.
Definition: structuralms.h:115
double giveTempStressWork()
Returns the temp density of total work of stress on strain increments.
Definition: concretedpm2.h:555
double computeTempKappa(double kappaInitial, double sigTrial, double rhoTrial, double sig)
Compute tempKappa.
ConcreteDPM2_ReturnType returnType
Definition: concretedpm2.h:594
virtual void give1dStressStiffMtrx(FloatMatrix &answer, MatResponseMode mode, GaussPoint *gp, TimeStep *tStep)
Method for computing 1d stiffness matrix of receiver.
double at(int i, int j) const
Coefficient access function.
Definition: floatmatrix.h:176
double stressWork
Density of total work done by stresses on strain increments.
Definition: concretedpm2.h:167
ConcreteDPM2Status(int n, Domain *d, GaussPoint *gp)
Constructor.
Definition: concretedpm2.C:54
virtual double computeDamageParamCompression(double equivStrain, double kappaOne, double kappaTwo, double omegaOld)
double wf
Control parameter for the linear/bilinear softening law in tension.
Definition: concretedpm2.h:667
virtual contextIOResultType restoreContext(DataStream &stream, ContextMode mode, void *obj=NULL)
Restores the receiver state previously written in stream.
Definition: structuralms.C:157
double computeDDGDDInv1d(double sigma, double tempKappa)
Here, the second derivative of the plastic potential with respect to the invariants sig and rho are c...
double yieldTol
yield tolerance for the plasticity model.
Definition: concretedpm2.h:676
double giveKappaP() const
Get the hardening variable of the plasticity model.
Definition: concretedpm2.h:224
double giveDamageTension() const
Get the tension damage variable of the damage model from the material status.
Definition: concretedpm2.h:292
void setTempDissWork(double w)
Sets the density of dissipated work to given value.
Definition: concretedpm2.h:563
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 ...
Abstract base class representing a material status information.
Definition: matstatus.h:84
Pair * add(int aKey, double value)
Adds a new Pair with given keyword and value into receiver.
Definition: dictionary.C:81
static double computeThirdCoordinate(const FloatArray &s)
Class representing vector of real numbers.
Definition: floatarray.h:82
This class implements the material status associated to ConcreteDPM2.
Definition: concretedpm2.h:82
double computeHardeningTwo(double tempKappa) const
Compute the value of the hardening function based on the hardening variable.
virtual ~ConcreteDPM2()
Destructor.
Definition: concretedpm2.C:513
double computeDDKappaDDeltaLambdaDKappa1d(double sig, double tempKappa)
FloatArray strainVector
Equilibrated strain vector in reduced form.
Definition: structuralms.h:69
Implementation of matrix containing floating point numbers.
Definition: floatmatrix.h:94
IRResultType
Type defining the return values of InputRecord reading operations.
Definition: irresulttype.h:47
double ftOne
Control parameter for the bilinear softening law.
Definition: concretedpm2.h:673
virtual contextIOResultType restoreContext(DataStream &, ContextMode mode, void *obj=NULL)
Restores the receiver state previously written in stream.
Definition: concretedpm2.C:353
void letTempStressVectorBe(const FloatArray &v)
Assigns tempStressVector to given vector v.
Definition: structuralms.h:135
#define _IFT_ConcreteDPM2_yieldtol
Definition: concretedpm2.h:64
double computeYieldValue(double sig, double rho, double theta, double tempKappa) const
Compute the yield value based on stress and hardening variable.
double DHard
Parameter of the ductilityMeasure of the plasticity model.
Definition: concretedpm2.h:613
double giveDamageCompression() const
Get the compressive damage variable of the damage model from the material status. ...
Definition: concretedpm2.h:301
double computeNorm() const
Computes the norm (or length) of the vector.
Definition: floatarray.C:840
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
double computeDDKappaDDeltaLambdaDKappa(double sig, double rho, double tempKappa)
Computes the derivative of the evolution law of the hardening parameter kappa with respect to the har...
void performPlasticityReturn(GaussPoint *gp, const FloatMatrix &D, const FloatArray &strain)
Perform stress return of the plasticity model and compute history variables.
double computeAlpha(FloatArray &effectiveStressTension, FloatArray &effectiveStressCompression, FloatArray &effectiveStress)
double rho
Length of the deviatoric stress.
Definition: concretedpm2.h:636
int softeningType
Type of softening function used.
Definition: concretedpm2.h:685
double computeHardeningOnePrime(double tempKappa) const
Compute the derivative of the hardening function based on the hardening parameter.
double giveTempDamageTension() const
Get the temp value of the hardening variable of the damage model from the material status...
Definition: concretedpm2.h:391
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
double yieldTolDamage
yield tolerance for the damage model.
Definition: concretedpm2.h:679
double giveLe()
Gives the characteristic length.
Definition: concretedpm2.h:532
int state_flag
Indicates the state (i.e. elastic, unloading, plastic, damage, vertex) of the Gauss point...
Definition: concretedpm2.h:163
double fc
Parameters of the yield surface of the plasticity model. fc is the uniaxial compressive strength...
Definition: concretedpm2.h:600
void times(double s)
Multiplies receiver with scalar.
Definition: floatarray.C:818
int strainRateFlag
Flag which signals if strainRate effects should be considered.
Definition: concretedpm2.h:694
#define _IFT_ConcreteDPM2_ft
Definition: concretedpm2.h:54
static double computeDeviatoricVolumetricSplit(FloatArray &dev, const FloatArray &s)
Computes split of receiver into deviatoric and volumetric part.
long ContextMode
Context mode (mask), defining the type of information written/read to/from context.
Definition: contextmode.h:43
#define _IFT_ConcreteDPM2_wfOne
Definition: concretedpm2.h:70
double rateStrain
Strains that are used for calculation of strain rates.
Definition: concretedpm2.h:159
Abstract base class for all "structural" constitutive models.
#define _IFT_ConcreteDPM2_fcZero
Definition: concretedpm2.h:53
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
Domain * giveDomain() const
Definition: femcmpnn.h:100
void letTempKappaPBe(double v)
Assign the temp value of the hardening variable of the plasticity model.
Definition: concretedpm2.h:441
void compute1dJacobian(FloatMatrix &answer, double totalsigma, double tempKappa, double deltaLambda, GaussPoint *gp)
Compute jacobian for 1D case.
double giveStressWork()
Returns the density of total work of stress on strain increments.
Definition: concretedpm2.h:553
double computeDKappaDDeltaLambda1d(double sig, double tempKappa)
int min(int i, int j)
Returns smaller value from two given decimals.
Definition: mathfem.h:59
double gM
Elastic shear modulus.
Definition: concretedpm2.h:657
double computeHardeningOne(double tempKappa) const
Compute the value of the hardening function based on the hardening variable.
#define _IFT_ConcreteDPM2_ftOne
Definition: concretedpm2.h:69
int giveStateFlag() const
Get the state flag from the material status.
Definition: concretedpm2.h:336
void letTempEquivStrainCompressionBe(double v)
Assign the temp value of the rate factor of the damage model.
Definition: concretedpm2.h:525
REGISTER_Material(DummyMaterial)
double computeDeltaPlasticStrainNormTension(double tempKappaD, double kappaD, GaussPoint *gp)
Compute equivalent strain value.
#define IR_GIVE_OPTIONAL_FIELD(__ir, __value, __id)
Macro facilitating the use of input record reading methods.
Definition: inputrecord.h:78
double sig
Volumetric stress.
Definition: concretedpm2.h:633
#define _IFT_ConcreteDPM2_bhard
Definition: concretedpm2.h:58
ConcreteDPM2_ReturnResult returnResult
Definition: concretedpm2.h:597
void letTempKappaDCompressionTwoBe(double v)
Assign the temp value of the second compression hardening variable of the damage model.
Definition: concretedpm2.h:483
int giveSize() const
Returns the size of receiver.
Definition: floatarray.h:218
double computeHardeningTwoPrime(double tempKappa) const
Compute the derivative of the hardening function based on the hardening parameter.
bool containsOnlyZeroes() const
Returns nonzero if all coefficients of the receiver are 0, else returns zero.
Definition: floatarray.C:646
virtual contextIOResultType saveContext(DataStream &stream, ContextMode mode, void *obj=NULL)
Stores receiver state to output stream.
Definition: structuralms.C:133
#define _IFT_ConcreteDPM2_ahard
Definition: concretedpm2.h:57
#define _IFT_ConcreteDPM2_efc
Definition: concretedpm2.h:67
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(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
void computeDRhoDStress(FloatArray &answer, const FloatArray &stress) const
Computes the derivative of rho with respect to the stress.
virtual double computeEquivalentStrain(double sig, double rho, double theta)
const FloatArray & giveStrainVector() const
Returns the const pointer to receiver&#39;s strain vector.
Definition: structuralms.h:105
#define _IFT_ConcreteDPM2_fc
Definition: concretedpm2.h:52
FloatArray tempStressVector
Temporary stress vector in reduced form (increments are used mainly in nonlinear analysis) ...
Definition: structuralms.h:73
static void transformStressVectorTo(FloatArray &answer, const FloatMatrix &base, const FloatArray &stressVector, bool transpose=false)
Transforms 3d stress vector into another coordinate system.
void letTempEquivStrainBe(double v)
Assign the temp value of the rate factor of the damage model.
Definition: concretedpm2.h:511
Class representing integration point in finite element program.
Definition: gausspoint.h:93
#define OOFEM_WARNING(...)
Definition: error.h:62
double computeDDGDInvDKappa1d(double sigma, double tempKappa)
virtual IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description stored in input record.
Definition: concretedpm2.C:519
ConcreteDPM2(int n, Domain *d)
Constructor.
Definition: concretedpm2.C:503
Class representing solution step.
Definition: timestep.h:80
double CHard
Parameter of the ductilityMeasure of the plasticity model.
Definition: concretedpm2.h:611
void add(const FloatArray &src)
Adds array src to receiver.
Definition: floatarray.C:156
double giveKappaDTensionTwo() const
Get the tension hardening variable two of the damage model from the material status.
Definition: concretedpm2.h:249
double dissWork
Density of dissipated work.
Definition: concretedpm2.h:171
double giveKappaDTensionOne() const
Get the hardening variable of the damage model from the material status.
Definition: concretedpm2.h:232
void resize(int s)
Resizes receiver towards requested size.
Definition: floatarray.C:631
void letTempKappaDTensionOneBe(double v)
Assign the temp value of the hardening variable of the damage model.
Definition: concretedpm2.h:462
double giveEquivStrainTension() const
Get the tension equivalent strain from the material status.
Definition: concretedpm2.h:275
void letTempEquivStrainTensionBe(double v)
Assign the temp value of the rate factor of the damage model.
Definition: concretedpm2.h:518

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