OOFEM  2.4
OOFEM.org - Object Oriented Finite Element Solver
mplasticmaterial2.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 "mplasticmaterial2.h"
36 #include "gausspoint.h"
37 #include "floatmatrix.h"
38 #include "floatarray.h"
39 #include "intarray.h"
40 #include "mathfem.h"
41 #include "datastream.h"
42 #include "contextioerr.h"
43 
44 namespace oofem {
45 #define YIELD_TOL 1.e-6
46 #define YIELD_TOL_SECONDARY 1.e-4
47 #define RES_TOL 1.e-8
48 #define PLASTIC_MATERIAL_MAX_ITERATIONS 120
49 
50 MPlasticMaterial2 :: MPlasticMaterial2(int n, Domain *d) : StructuralMaterial(n, d)
51  //
52  // constructor
53  //
54 {
55  linearElasticMaterial = NULL;
56  rmType = mpm_ClosestPoint;
57  plType = associatedPT;
58  iterativeUpdateOfActiveConds = false; // safe
59 }
60 
61 
62 
63 MPlasticMaterial2 :: ~MPlasticMaterial2()
64 //
65 // destructor
66 //
67 {
68  delete linearElasticMaterial;
69 }
70 
71 
72 int
73 MPlasticMaterial2 :: hasMaterialModeCapability(MaterialMode mode)
74 //
75 // returns whether receiver supports given mode
76 //
77 {
78  return mode == _3dMat ||
79  mode == _1dMat ||
80  //<RESTRICTED_SECTION>
81  mode == _PlaneStress ||
82  //</RESTRICTED_SECTION>
83  mode == _PlaneStrain;
84 }
85 
86 
87 MaterialStatus *
88 MPlasticMaterial2 :: CreateStatus(GaussPoint *gp) const
89 {
90  return new MPlasticMaterial2Status(1, this->giveDomain(), gp, this->giveSizeOfReducedHardeningVarsVector(gp));
91 }
92 
93 
94 void
95 MPlasticMaterial2 :: giveRealStressVector(FloatArray &answer,
96  GaussPoint *gp,
97  const FloatArray &totalStrain,
98  TimeStep *tStep)
99 //
100 // returns real stress vector in 3d stress space of receiver according to
101 // previous level of stress and current
102 // strain increment, the only way, how to correctly update gp records
103 //
104 // completely formulated in Reduced stress-strain space
105 {
106  FloatArray strainSpaceHardeningVariables;
107  FloatArray fullStressVector;
108  FloatArray strainVectorR, plasticStrainVectorR;
109  FloatArray helpVec;
110  IntArray activeConditionMap(this->nsurf);
111  FloatArray gamma;
112 
113  MPlasticMaterial2Status *status = static_cast< MPlasticMaterial2Status * >( this->giveStatus(gp) );
114 
115  this->initTempStatus(gp);
116 
117  // subtract stress independent part
118  // note: eigenStrains (temperature) is not contained in mechanical strain stored in gp
119  // therefore it is necessary to subtract always the total eigen strain value
120  this->giveStressDependentPartOfStrainVector(strainVectorR, gp, totalStrain,
121  tStep, VM_Total);
122 
123  /*
124  * stress return algorithm
125  */
126 
127  if ( rmType == mpm_ClosestPoint ) {
128  this->closestPointReturn(fullStressVector, activeConditionMap, gamma, gp, strainVectorR,
129  plasticStrainVectorR, strainSpaceHardeningVariables, tStep);
130  } else {
131  this->cuttingPlaneReturn(fullStressVector, activeConditionMap, gamma, gp, strainVectorR,
132  plasticStrainVectorR, strainSpaceHardeningVariables, tStep);
133  }
134 
135  status->letTempStrainVectorBe(totalStrain);
136  StructuralMaterial :: giveReducedSymVectorForm( helpVec, fullStressVector, gp->giveMaterialMode() );
137 
138  //perform isotropic damage on effective stress
139  double omega = computeDamage(gp, strainSpaceHardeningVariables, tStep);
140  status->letTempDamageBe(omega);
141  helpVec.times(1. - omega);
142  //end of damage part
143 
144  status->letTempStressVectorBe(helpVec);
145 
146  status->letTempPlasticStrainVectorBe(plasticStrainVectorR);
147  status->letTempStrainSpaceHardeningVarsVectorBe(strainSpaceHardeningVariables);
148 
149  status->setTempGamma(gamma);
150  status->setTempActiveConditionMap(activeConditionMap);
151 
152 
153  // update state flag
154  int newState, state = status->giveStateFlag();
155  bool yieldFlag = false;
156  for ( int i = 1; i <= nsurf; i++ ) {
157  if ( gamma.at(i) > 0. ) {
158  yieldFlag = true;
159  }
160  }
161 
162  if ( yieldFlag ) {
163  newState = MPlasticMaterial2Status :: PM_Yielding; // test if plastic loading occur
164  }
165  // no plastic loading - check for unloading
166  else if ( ( state == MPlasticMaterial2Status :: PM_Yielding ) || ( state == MPlasticMaterial2Status :: PM_Unloading ) ) {
167  newState = MPlasticMaterial2Status :: PM_Unloading;
168  } else {
169  newState = MPlasticMaterial2Status :: PM_Elastic;
170  }
171 
172  status->letTempStateFlagBe(newState);
173 
174  StructuralMaterial :: giveReducedSymVectorForm( answer, fullStressVector, gp->giveMaterialMode() );
175 }
176 
177 double
178 MPlasticMaterial2 :: computeDamage(GaussPoint *gp, const FloatArray &strainSpaceHardeningVariables, TimeStep *tStep) {
179  return 0.;
180 }
181 
182 
183 void
184 MPlasticMaterial2 :: closestPointReturn(FloatArray &answer,
185  IntArray &activeConditionMap,
186  FloatArray &gamma,
187  GaussPoint *gp,
188  const FloatArray &totalStrain,
189  FloatArray &plasticStrainVectorR,
190  FloatArray &strainSpaceHardeningVariables,
191  TimeStep *tStep)
192 {
193  FloatArray fullStressVector;
194  FloatArray elasticStrainVectorR;
195  FloatArray residualVectorR;
196  FloatArray helpVector, helpVector2;
197  FloatArray dgamma, dgammared, tempGamma, dkappa;
198  FloatMatrix elasticModuli, consistentModuli;
199  FloatMatrix elasticModuliInverse;
200  FloatMatrix helpMtrx, helpMtrx2;
201  FloatMatrix ks, kl, lmat, rmat, gradientMatrix;
202  FloatMatrix gmat;
203  IntArray initialConditionMap;
204  std :: vector< FloatArray > yieldGradSigVec(this->nsurf), loadGradSigVec(this->nsurf), * loadGradSigVecPtr;
205  std :: vector< FloatArray > yieldGradKVec(this->nsurf), loadGradKVec(this->nsurf);
206  FloatArray rhs;
207  double yieldValue;
208  int nIterations = 0;
209  int strSize;
210  int elastic, restart, actSurf, indx;
211  bool yieldConsistency, init = true;
212  bool hasHardening = this->hasHardening() > 0;
213 #ifdef DEBUG
214  bool debug = false;
215 #endif
216 
217  this->clearPopulationSet();
218 
219  if ( this->plType == associatedPT ) {
220  loadGradSigVecPtr = & yieldGradSigVec;
221  } else {
222  loadGradSigVecPtr = & loadGradSigVec;
223  }
224 
225  MPlasticMaterial2Status *status = static_cast< MPlasticMaterial2Status * >( this->giveStatus(gp) );
226 
227 huhu: //label for goto
228 
229  plasticStrainVectorR = status->givePlasticStrainVector();
230  strainSpaceHardeningVariables = status->giveStrainSpaceHardeningVars();
231 
232  dgamma.resize(nsurf);
233  dgamma.zero();
234  gamma.resize(nsurf);
235  gamma.zero();
236 
237  strSize = totalStrain.giveSize(); // size of reducedStrain Vector
238 
239  // compute elastic moduli and its inverse
240  this->computeReducedElasticModuli(elasticModuli, gp, tStep);
241  elasticModuliInverse.beInverseOf(elasticModuli);
242 
243  //
244  // Elastic Predictor
245  //
246 
247  elasticStrainVectorR = totalStrain;
248  elasticStrainVectorR.subtract(plasticStrainVectorR);
249  // stress vector in full form due to computational convinience
250  //if (fullStressVector) delete fullStressVector;
251  this->computeTrialStressIncrement(fullStressVector, gp, elasticStrainVectorR, tStep);
252  //
253  // check for plastic process
254  //
255  elastic = 1;
256  actSurf = 0;
257  activeConditionMap.zero();
258  if ( init ) {
259  initialConditionMap.resize(nsurf);
260  initialConditionMap.zero();
261  }
262 
263  for ( int i = 1; i <= this->nsurf; i++ ) {
264  yieldValue = this->computeYieldValueAt(gp, i, fullStressVector, strainSpaceHardeningVariables);
265  if ( yieldValue > YIELD_TOL ) {
266  elastic = 0;
267  if ( init ) {
268  initialConditionMap.at(i) = 1;
269  }
270 
271  if ( actSurf < this->giveMaxNumberOfActiveYieldConds(gp) ) {
272  actSurf += 1;
273  activeConditionMap.at(i) = actSurf;
274  }
275  }
276  }
277 
278  init = false;
279 
280  if ( !elastic ) {
281  /*
282  * plastic multisurface closest-point corrector
283  */
284  this->addNewPopulation(activeConditionMap);
285 
286 #ifdef DEBUG
287  if ( debug ) {
288  printf("New activeConditionMap:");
289  activeConditionMap.printYourself();
290  printf("\n");
291  }
292 
293 #endif
294 
295  do {
296  do { // restart loop
297  // compute gradients
298  for ( int i = 1; i <= nsurf; i++ ) {
299  this->computeReducedStressGradientVector(yieldGradSigVec [ i - 1 ], yieldFunction, i, gp, fullStressVector,
300  strainSpaceHardeningVariables);
301  if ( hasHardening ) {
302  this->computeKGradientVector(yieldGradKVec [ i - 1 ], yieldFunction, i, gp, fullStressVector,
303  strainSpaceHardeningVariables);
304  }
305  }
306 
307  if ( this->plType == nonassociatedPT ) {
308  for ( int i = 1; i <= nsurf; i++ ) {
309  this->computeReducedStressGradientVector(loadGradSigVec [ i - 1 ], loadFunction, i, gp, fullStressVector,
310  strainSpaceHardeningVariables);
311  if ( hasHardening ) {
312  this->computeKGradientVector(loadGradKVec [ i - 1 ], loadFunction, i, gp, fullStressVector,
313  strainSpaceHardeningVariables);
314  }
315  }
316  }
317 
318 
319  //this->computeGradientVector (gradientVectorR, gp, &fullStressVector, fullStressSpaceHardeningVars);
320  this->computeResidualVector(residualVectorR, gp, gamma, activeConditionMap, plasticStrainVectorR,
321  strainSpaceHardeningVariables, * loadGradSigVecPtr);
322 #ifdef DEBUG
323  if ( debug ) {
324  printf("Convergence[actSet: ");
325  for ( int i = 1; i <= nsurf; i++ ) {
326  printf( "%d ", activeConditionMap.at(i) );
327  }
328 
329  printf("] yield_val (gamma) ");
330  }
331 
332 #endif
333 
334  // check convergence
335  yieldConsistency = true;
336  for ( int i = 1; i <= nsurf; i++ ) {
337  if ( activeConditionMap.at(i) ) {
338  yieldValue = this->computeYieldValueAt(gp, i, fullStressVector, strainSpaceHardeningVariables);
339  if ( fabs(yieldValue) >= YIELD_TOL ) {
340  yieldConsistency = false;
341  }
342 
343 #ifdef DEBUG
344  if ( debug ) {
345  printf( "%e(%e) ", yieldValue, gamma.at(i) );
346  }
347 
348 #endif
349  }
350  }
351 
352 #ifdef DEBUG
353  if ( debug ) {
354  printf("\n");
355  }
356 
357 #endif
358 
359  //if (yieldConsistency && (sqrt(dotProduct(residualVectorR, residualVectorR, residualVectorR.giveSize())) < RES_TOL)) {
360  if ( yieldConsistency ) {
361  answer = fullStressVector;
362  //printf (" (%d iterations)", nIterations);
363 
364 
365  restart = 0;
366  // check for active conditions
367 #ifdef DEBUG
368  if ( debug ) {
369  printf("Consistency: ");
370  }
371 
372 #endif
373  for ( int i = 1; i <= this->nsurf; i++ ) {
374  yieldValue = this->computeYieldValueAt(gp, i, fullStressVector, strainSpaceHardeningVariables);
375 
376 #ifdef DEBUG
377  if ( debug ) {
378  printf("%e ", yieldValue);
379  }
380 
381 #endif
382 
383  if ( ( gamma.at(i) < 0.0 ) || ( ( activeConditionMap.at(i) == 0 ) && ( yieldValue > YIELD_TOL_SECONDARY ) ) ) {
384  restart = 1;
385  //printf ("*");
386  }
387  }
388 
389 #ifdef DEBUG
390  if ( debug ) {
391  printf("\nStress: ");
392  for ( int i = 1; i <= fullStressVector.giveSize(); i++ ) {
393  printf( " %e", fullStressVector.at(i) );
394  }
395 
396  printf("\nKappa : ");
397  for ( int i = 1; i <= strainSpaceHardeningVariables.giveSize(); i++ ) {
398  printf( " %e", strainSpaceHardeningVariables.at(i) );
399  }
400 
401  printf("\n Restart %d\n", restart);
402  }
403 
404 #endif
405 
406  if ( restart ) {
407  IntArray newmap(nsurf);
408  actSurf = 0;
409  newmap.zero();
410  // buid new actsurf mask
411  for ( int i = 1; i <= this->nsurf; i++ ) {
412  yieldValue = this->computeYieldValueAt(gp, i, fullStressVector, strainSpaceHardeningVariables);
413  if ( ( ( activeConditionMap.at(i) ) && ( gamma.at(i) >= 0.0 ) ) || ( yieldValue > YIELD_TOL ) ) {
414  elastic = 0;
415  actSurf += 1;
416  newmap.at(i) = actSurf;
417  // if (actSurf == this->giveMaxNumberOfActiveYieldConds(gp)) break;
418  }
419  }
420 
421  if ( this->getNewPopulation(activeConditionMap, newmap, this->giveMaxNumberOfActiveYieldConds(gp), nsurf) == 0 ) {
422 #ifdef DEBUG
423  if ( !debug ) {
424  this->clearPopulationSet();
425  debug = true;
426  goto huhu;
427  } else {
428  printf("ep:");
429  plasticStrainVectorR = status->givePlasticStrainVector();
430  plasticStrainVectorR.printYourself();
431  printf("kp:");
432  strainSpaceHardeningVariables = status->giveStrainSpaceHardeningVars();
433  strainSpaceHardeningVariables.printYourself();
434  elasticStrainVectorR = totalStrain;
435  elasticStrainVectorR.subtract(plasticStrainVectorR);
436  printf("sg:");
437  this->computeTrialStressIncrement(fullStressVector, gp, elasticStrainVectorR, tStep);
438  fullStressVector.printYourself();
439 
440  OOFEM_ERROR("Internal Consistency error: all combinations of yield functions tried, no consistent return");
441  }
442 
443 #else
444  OOFEM_ERROR("Internal Consistency error: all combinations of yield functions tried, no consistent return");
445 #endif
446  }
447 
448  actSurf = 0;
449  for ( int i = 1; i <= nsurf; i++ ) {
450  if ( activeConditionMap.at(i) ) {
451  actSurf++;
452  }
453  }
454 
455  /*
456  * if (actSurf) {
457  * activeConditionMap=newmap;
458  * } else {
459  * restart = 0;
460  * // select any possible
461  * for (i=1; i<=this->nsurf; i++) {
462  * if (initialConditionMap.at(i)) {
463  * actSurf = 1;
464  * activeConditionMap.zero();
465  * activeConditionMap.at(i) = 1;
466  * initialConditionMap.at(i) = 0;
467  * restart = 1;
468  * break;
469  * }
470  * }
471  *
472  * if (!restart) {
473  * OOFEM_ERROR("Internal Consistency error");
474  * }
475  * }
476  */
477  plasticStrainVectorR = status->givePlasticStrainVector();
478  strainSpaceHardeningVariables = status->giveStrainSpaceHardeningVars();
479  elasticStrainVectorR = totalStrain;
480  elasticStrainVectorR.subtract(plasticStrainVectorR);
481  this->computeTrialStressIncrement(fullStressVector, gp, elasticStrainVectorR, tStep);
482  gamma.zero();
483  nIterations = 0;
484  int restart = 0;
485  if ( restart ) {
486  goto huhu;
487  }
488 
489  continue;
490  }
491 
492  if ( yieldConsistency && ( residualVectorR.computeNorm() < RES_TOL ) ) {
493  return;
494  }
495  }
496 
497  // compute consistent tangent moduli
498  this->computeAlgorithmicModuli(consistentModuli, gp, elasticModuliInverse,
499  gamma, activeConditionMap,
500  fullStressVector, strainSpaceHardeningVariables);
501 
502 
504 
505 
506  rhs.resize(actSurf);
507  rhs.zero();
508  gmat.resize(actSurf, actSurf);
509  gmat.zero();
510  lmat.resize(actSurf, strSize);
511  lmat.zero();
512  rmat.resize(strSize, actSurf);
513  rmat.zero();
514  if ( hasHardening ) {
515  this->computeReducedHardeningVarsSigmaGradient(ks, gp, activeConditionMap, fullStressVector,
516  strainSpaceHardeningVariables, gamma);
517  this->computeReducedHardeningVarsLamGradient(kl, gp, actSurf, activeConditionMap,
518  fullStressVector, strainSpaceHardeningVariables, gamma);
519  }
520 
521  for ( int i = 1; i <= nsurf; i++ ) {
522  if ( ( indx = activeConditionMap.at(i) ) ) {
523  if ( hasHardening ) {
524  helpVector.beTProductOf(ks, yieldGradKVec [ i - 1 ]);
525  helpVector.add(yieldGradSigVec [ i - 1 ]);
526  for ( int j = 1; j <= strSize; j++ ) {
527  lmat.at(indx, j) = helpVector.at(j);
528  }
529  } else {
530  for ( int j = 1; j <= strSize; j++ ) {
531  lmat.at(indx, j) = yieldGradSigVec [ i - 1 ].at(j);
532  }
533  }
534 
535  if ( hasHardening ) {
536  this->computeReducedSKGradientMatrix(gradientMatrix, i, gp, fullStressVector,
537  strainSpaceHardeningVariables);
538  helpMtrx.beProductOf(gradientMatrix, kl);
539  helpMtrx.times( gamma.at(i) );
540  rmat.add(helpMtrx);
541  }
542 
543  for ( int j = 1; j <= strSize; j++ ) {
544  rmat.at(j, indx) += ( * loadGradSigVecPtr ) [ i - 1 ].at(j);
545  }
546 
547  if ( hasHardening ) {
548  helpVector2.beTProductOf(kl, yieldGradKVec [ i - 1 ]);
549  for ( int j = 1; j <= actSurf; j++ ) {
550  gmat.at(indx, j) = ( -1.0 ) * helpVector2.at(j);
551  }
552  }
553 
554  yieldValue = this->computeYieldValueAt(gp, i, fullStressVector, strainSpaceHardeningVariables);
555  rhs.at(indx) = yieldValue;
556  //helpVector2.beTProductOf (consistentModuli, helpVector);
557  //rhs.at(indx)=yieldValue-dotProduct(helpVector2, residualVectorR, residualVectorR.giveSize());
558  }
559  }
560 
561  helpMtrx.beProductOf(lmat, consistentModuli);
562  helpMtrx2.beProductOf(helpMtrx, rmat);
563  gmat.add(helpMtrx2);
564  helpVector2.beProductOf(consistentModuli, residualVectorR);
565  helpVector.beProductOf(lmat, helpVector2);
566  rhs.subtract(helpVector);
567 
568 
569  // obtain increment to consistency parameter
570  gmat.solveForRhs(rhs, dgammared);
571 
572  // assign gammas
573  for ( int i = 1; i <= this->nsurf; i++ ) {
574  if ( ( indx = activeConditionMap.at(i) ) ) {
575  dgamma.at(i) = dgammared.at(indx);
576  } else {
577  dgamma.at(i) = 0.0;
578  }
579  }
580 
581  tempGamma = gamma;
582  tempGamma.add(dgamma);
583 
584 
585  if ( iterativeUpdateOfActiveConds ) {
586  //
587  // check active constraint
588  //
589 
590  restart = 0;
591  for ( int i = 1; i <= this->nsurf; i++ ) {
592  if ( tempGamma.at(i) < 0.0 ) {
593  restart = 1;
594  }
595  }
596 
597  if ( restart ) {
598  // build up new activeConditionMap and restart the iterartion from begining
599  actSurf = 0;
600  activeConditionMap.zero();
601  for ( int i = 1; i <= this->nsurf; i++ ) {
602  if ( tempGamma.at(i) > 0.0 ) {
603  actSurf += 1;
604  activeConditionMap.at(i) = actSurf;
605  }
606  }
607 
608  plasticStrainVectorR = status->givePlasticStrainVector();
609  strainSpaceHardeningVariables = status->giveStrainSpaceHardeningVars();
610  elasticStrainVectorR = totalStrain;
611  elasticStrainVectorR.subtract(plasticStrainVectorR);
612  this->computeTrialStressIncrement(fullStressVector, gp, elasticStrainVectorR, tStep);
613  gamma.zero();
614  nIterations = 0;
615  continue;
616  } else {
617  break;
618  }
619  } else {
620  break;
621  }
622  } while ( 1 ); // end restart loop
623 
624  // obtain incremental plastic strains and internal variables
625  // residualVectorR vector and gradVec array are changed
626  // but they are computed once again when new iteration starts
627 
628  helpVector2.resize(strSize);
629  // Update state variables and consistency parameter
630  for ( int i = 1; i <= nsurf; i++ ) {
631  if ( activeConditionMap.at(i) ) {
632  ( * loadGradSigVecPtr ) [ i - 1 ].times( dgamma.at(i) );
633  residualVectorR.add( ( * loadGradSigVecPtr ) [ i - 1 ] );
634 
635  if ( hasHardening ) {
636  this->computeReducedSKGradientMatrix(gradientMatrix, i, gp, fullStressVector,
637  strainSpaceHardeningVariables);
638  this->computeReducedHardeningVarsLamGradient(kl, gp, actSurf, activeConditionMap,
639  fullStressVector, strainSpaceHardeningVariables, gamma);
640  helpMtrx.beProductOf(gradientMatrix, kl);
641  helpMtrx.times( gamma.at(i) );
642  helpVector2.beProductOf(helpMtrx, dgammared);
643 
644  residualVectorR.add(helpVector2);
645  }
646  }
647  }
648 
649  helpVector.beProductOf(consistentModuli, residualVectorR);
650  helpVector2.beProductOf(elasticModuliInverse, helpVector);
651 
652  gamma.add(dgamma);
653  for ( int i = 1; i <= strSize; i++ ) {
654  plasticStrainVectorR.at(i) += helpVector2.at(i);
655  }
656 
657  elasticStrainVectorR = totalStrain;
658  elasticStrainVectorR.subtract(plasticStrainVectorR);
659  // stress vector in full form due to computational convinience
660  this->computeTrialStressIncrement(fullStressVector, gp, elasticStrainVectorR, tStep);
661 
662  if ( hasHardening ) {
663  this->computeStrainHardeningVarsIncrement(dkappa, gp, fullStressVector, gamma, helpVector2, activeConditionMap);
664  strainSpaceHardeningVariables = status->giveStrainSpaceHardeningVars();
665  strainSpaceHardeningVariables.add(dkappa);
666  }
667 
668  //this->computeStrainHardeningVarsIncrement(dkappa, fullStressVector, dgamma);
669  //strainSpaceHardeningVariables.add(dkappa);
670  //this-> computeTrialStressIncrement (fullStressVector, gp, elasticStrainVectorR, tStep);
671 
672  // increment iteration count
673  nIterations++;
674 
675  if ( nIterations > PLASTIC_MATERIAL_MAX_ITERATIONS ) {
676 #ifdef DEBUG
677  // issue warning
678  fprintf( stderr, "\nclosestPointReturn (mat no. %d): reached max number of iterations\n", this->giveNumber() );
679  fprintf(stderr, "activeSet: ");
680  for ( int i = 1; i <= nsurf; i++ ) {
681  fprintf( stderr, "%d ", activeConditionMap.at(i) );
682  }
683 
684  fprintf(stderr, "\n");
685  if ( !debug ) {
686  debug = true;
687  this->clearPopulationSet();
688  goto huhu;
689  }
690 
691 #endif
692 
693 
694  // try last resort select single function from initial active set
695  restart = 0;
696  bool smartCandidate = false;
697  int candidate = 0;
698  IntArray newmap(nsurf);
699  actSurf = 0;
700  newmap.zero();
701  // check for active conditions
702  for ( int i = 1; i <= this->nsurf; i++ ) {
703  if ( activeConditionMap.at(i) ) {
704  if ( gamma.at(i) < 0.0 ) {
705  smartCandidate = true;
706  } else {
707  candidate = i;
708  }
709  }
710  }
711 
712  if ( smartCandidate && candidate && initialConditionMap.at(candidate) ) {
713  actSurf = 1;
714  IntArray newMap(nsurf);
715  newMap.at(candidate) = 1;
716 
717  //initialConditionMap.at(candidate) = 0;
718  if ( this->getNewPopulation(activeConditionMap, newmap, this->giveMaxNumberOfActiveYieldConds(gp), nsurf) == 0 ) {
719 #ifdef DEBUG
720  if ( !debug ) {
721  this->clearPopulationSet();
722  debug = true;
723  goto huhu;
724  } else {
725  printf("ep:");
726  plasticStrainVectorR = status->givePlasticStrainVector();
727  plasticStrainVectorR.printYourself();
728  printf("kp:");
729  strainSpaceHardeningVariables = status->giveStrainSpaceHardeningVars();
730  strainSpaceHardeningVariables.printYourself();
731  elasticStrainVectorR = totalStrain;
732  elasticStrainVectorR.subtract(plasticStrainVectorR);
733  printf("sg:");
734  this->computeTrialStressIncrement(fullStressVector, gp, elasticStrainVectorR, tStep);
735  fullStressVector.printYourself();
736 
737  OOFEM_ERROR("Internal Consistency error: all combinations of yield functions tried, no consistent return");
738  }
739 
740 #else
741  OOFEM_ERROR("Internal Consistency error: all combinations of yield functions tried, no consistent return");
742 #endif
743  }
744 
745  actSurf = 0;
746  for ( int i = 1; i <= nsurf; i++ ) {
747  if ( activeConditionMap.at(i) ) {
748  actSurf++;
749  }
750  }
751 
752  restart = 1;
753  } else {
754  IntArray newMap(nsurf);
755  // select any possible
756  for ( int i = 1; i <= this->nsurf; i++ ) {
757  if ( initialConditionMap.at(i) ) {
758  actSurf = 1;
759  newMap.at(i) = 1;
760  //initialConditionMap.at(i) = 0;
761  if ( this->getNewPopulation(activeConditionMap, newmap, this->giveMaxNumberOfActiveYieldConds(gp), nsurf) == 0 ) {
762 #ifdef DEBUG
763  if ( !debug ) {
764  this->clearPopulationSet();
765  debug = true;
766  goto huhu;
767  } else {
768  printf("ep:");
769  plasticStrainVectorR = status->givePlasticStrainVector();
770  plasticStrainVectorR.printYourself();
771  printf("kp:");
772  strainSpaceHardeningVariables = status->giveStrainSpaceHardeningVars();
773  strainSpaceHardeningVariables.printYourself();
774  elasticStrainVectorR = totalStrain;
775  elasticStrainVectorR.subtract(plasticStrainVectorR);
776  printf("sg:");
777  this->computeTrialStressIncrement(fullStressVector, gp, elasticStrainVectorR, tStep);
778  fullStressVector.printYourself();
779 
780  OOFEM_ERROR("Internal Consistency error: all combinations of yield functions tried, no consistent return");
781  }
782 
783 #else
784  OOFEM_ERROR("Internal Consistency error: all combinations of yield functions tried, no consistent return");
785 #endif
786  }
787 
788  for ( actSurf = 0, i = 1; i <= nsurf; i++ ) {
789  if ( activeConditionMap.at(i) ) {
790  actSurf++;
791  }
792  }
793 
794  restart = 1;
795  break;
796  }
797  }
798  }
799 
800  if ( restart ) {
801  //fprintf(stderr,"===>LAST RESORT<====");
802 
803  plasticStrainVectorR = status->givePlasticStrainVector();
804  strainSpaceHardeningVariables = status->giveStrainSpaceHardeningVars();
805  elasticStrainVectorR = totalStrain;
806  elasticStrainVectorR.subtract(plasticStrainVectorR);
807  this->computeTrialStressIncrement(fullStressVector, gp, elasticStrainVectorR, tStep);
808  gamma.zero();
809  nIterations = 0;
810  int restart = 0;
811  if ( restart ) {
812  goto huhu;
813  }
814 
815  continue;
816  } else {
817  OOFEM_WARNING("local equlibrium not reached in %d iterations\nElement %d, gp %d, continuing",
818  PLASTIC_MATERIAL_MAX_ITERATIONS, gp->giveElement()->giveNumber(), gp->giveNumber() );
819  answer = fullStressVector;
820  // debug line
821  nIterations = 0;
822  goto huhu;
823  }
824  }
825  } while ( 1 );
826  } else {
827  answer = fullStressVector;
828  }
829 }
830 
831 void
832 MPlasticMaterial2 :: cuttingPlaneReturn(FloatArray &answer,
833  IntArray &activeConditionMap,
834  FloatArray &gamma,
835  GaussPoint *gp,
836  const FloatArray &totalStrain,
837  FloatArray &plasticStrainVectorR,
838  FloatArray &strainSpaceHardeningVariables,
839  TimeStep *tStep)
840 {
841  FloatArray elasticStrainVectorR;
842  FloatArray fullStressVector;
843  FloatArray helpVector, helpVector2, rhs, dkappa;
844  FloatMatrix elasticModuli, helpMtrx, helpMtrx2, gmat;
845  FloatMatrix kl, ks, lmat, rmat;
846  IntArray initialConditionMap;
847  std :: vector< FloatArray > yieldGradSigVec(this->nsurf), loadGradSigVec(this->nsurf), * loadGradSigVecPtr;
848  std :: vector< FloatArray > yieldGradKVec(this->nsurf), loadGradKVec(this->nsurf);
849  FloatArray dgamma(this->nsurf);
850  double yieldValue;
851  int nIterations = 0;
852  int strSize, i, j, elastic, restart, actSurf, indx;
853  bool yieldConsistency, init = true;
854  bool hasHardening = this->hasHardening() > 0;
855 
856  MPlasticMaterial2Status *status = static_cast< MPlasticMaterial2Status * >( this->giveStatus(gp) );
857 
858  if ( this->plType == associatedPT ) {
859  loadGradSigVecPtr = & yieldGradSigVec;
860  } else {
861  loadGradSigVecPtr = & loadGradSigVec;
862  }
863 
864 huhu: //label for goto
865  plasticStrainVectorR = status->givePlasticStrainVector();
866  strainSpaceHardeningVariables = status->giveStrainSpaceHardeningVars();
867 
868 
869  dgamma.resize(nsurf);
870  dgamma.zero();
871  gamma.resize(nsurf);
872  gamma.zero();
873  // compute elastic moduli
874  this->computeReducedElasticModuli(elasticModuli, gp, tStep);
875  strSize = elasticModuli.giveNumberOfRows();
876 
877  // initialize activeConditionMap
878  elasticStrainVectorR = totalStrain;
879  elasticStrainVectorR.subtract(plasticStrainVectorR);
880  // stress vector in full form due to computational convinience
881  this->computeTrialStressIncrement(fullStressVector, gp, elasticStrainVectorR, tStep);
882 
883  elastic = 1;
884  actSurf = 0;
885  activeConditionMap.zero();
886  if ( init ) {
887  initialConditionMap.resize(nsurf);
888  initialConditionMap.zero();
889  }
890 
891  for ( i = 1; i <= this->nsurf; i++ ) {
892  yieldValue = this->computeYieldValueAt(gp, i, fullStressVector, strainSpaceHardeningVariables);
893  if ( yieldValue > YIELD_TOL ) {
894  elastic = 0;
895  if ( init ) {
896  initialConditionMap.at(i) = 1;
897  }
898 
899  if ( actSurf < this->giveMaxNumberOfActiveYieldConds(gp) ) {
900  actSurf += 1;
901  activeConditionMap.at(i) = actSurf;
902  }
903  }
904  }
905 
906  init = false;
907 
908  if ( !elastic ) {
909  do { // local equilibrium loop
910  do { // loop for restart
911  // compute gradients
912  for ( i = 1; i <= nsurf; i++ ) {
913  this->computeReducedStressGradientVector(yieldGradSigVec [ i - 1 ], yieldFunction, i, gp, fullStressVector,
914  strainSpaceHardeningVariables);
915 
916  if ( hasHardening ) {
917  this->computeKGradientVector(yieldGradKVec [ i - 1 ], yieldFunction, i, gp, fullStressVector,
918  strainSpaceHardeningVariables);
919  }
920  }
921 
922  if ( this->plType == nonassociatedPT ) {
923  for ( i = 1; i <= nsurf; i++ ) {
924  this->computeReducedStressGradientVector(loadGradSigVec [ i - 1 ], loadFunction, i, gp, fullStressVector,
925  strainSpaceHardeningVariables);
926  if ( hasHardening ) {
927  this->computeKGradientVector(loadGradKVec [ i - 1 ], loadFunction, i, gp, fullStressVector,
928  strainSpaceHardeningVariables);
929  }
930  }
931  }
932 
933  rhs.resize(actSurf);
934  rhs.zero();
935  gmat.resize(actSurf, actSurf);
936  gmat.zero();
937  lmat.resize(actSurf, strSize);
938  lmat.zero();
939  rmat.resize(strSize, actSurf);
940  rmat.zero();
941 
942  if ( hasHardening ) {
943  this->computeReducedHardeningVarsSigmaGradient(ks, gp, activeConditionMap, fullStressVector,
944  strainSpaceHardeningVariables, gamma);
945  this->computeReducedHardeningVarsLamGradient(kl, gp, actSurf, activeConditionMap,
946  fullStressVector, strainSpaceHardeningVariables, gamma);
947  }
948 
949  for ( i = 1; i <= nsurf; i++ ) {
950  if ( ( indx = activeConditionMap.at(i) ) ) {
951  rhs.at(indx) = this->computeYieldValueAt(gp, i, fullStressVector, strainSpaceHardeningVariables);
952 
953  if ( hasHardening ) {
954  helpVector.beTProductOf(ks, yieldGradKVec [ i - 1 ]);
955  helpVector.add(yieldGradSigVec [ i - 1 ]);
956  for ( j = 1; j <= strSize; j++ ) {
957  lmat.at(indx, j) = helpVector.at(j);
958  }
959  } else {
960  for ( j = 1; j <= strSize; j++ ) {
961  lmat.at(indx, j) = yieldGradSigVec [ i - 1 ].at(j);
962  }
963  }
964 
965  for ( j = 1; j <= strSize; j++ ) {
966  rmat.at(j, indx) += ( * loadGradSigVecPtr ) [ i - 1 ].at(j);
967  }
968 
969  if ( hasHardening ) {
970  helpVector2.beTProductOf(kl, yieldGradKVec [ i - 1 ]);
971  for ( j = 1; j <= actSurf; j++ ) {
972  gmat.at(indx, j) = ( -1.0 ) * helpVector2.at(j);
973  }
974  }
975  }
976  }
977 
978  helpMtrx.beProductOf(lmat, elasticModuli);
979  helpMtrx2.beProductOf(helpMtrx, rmat);
980  gmat.add(helpMtrx2);
981 
982  // solve for plastic multiplier increments
983  gmat.solveForRhs(rhs, helpVector);
984  // assign gammas
985  for ( i = 1; i <= this->nsurf; i++ ) {
986  if ( ( indx = activeConditionMap.at(i) ) ) {
987  dgamma.at(i) = helpVector.at(indx);
988  } else {
989  dgamma.at(i) = 0.0;
990  }
991  }
992 
993  if ( iterativeUpdateOfActiveConds ) {
994  restart = 0;
995  // check for active conditions
996  for ( i = 1; i <= this->nsurf; i++ ) {
997  if ( ( gamma.at(i) + dgamma.at(i) ) < 0.0 ) {
998  gamma.at(i) = 0.0;
999  restart = 1;
1000  }
1001  }
1002 
1003  if ( restart ) {
1004  // build up new activeConditionMap and restart the iterartion from begining
1005  actSurf = 0;
1006  activeConditionMap.zero();
1007  for ( i = 1; i <= this->nsurf; i++ ) {
1008  if ( ( gamma.at(i) + dgamma.at(i) ) > 0.0 ) {
1009  actSurf += 1;
1010  activeConditionMap.at(i) = actSurf;
1011  }
1012  }
1013 
1014  plasticStrainVectorR = status->givePlasticStrainVector();
1015  strainSpaceHardeningVariables = status->giveStrainSpaceHardeningVars();
1016  elasticStrainVectorR = totalStrain;
1017  elasticStrainVectorR.subtract(plasticStrainVectorR);
1018  this->computeTrialStressIncrement(fullStressVector, gp, elasticStrainVectorR, tStep);
1019  gamma.zero();
1020  nIterations = 0;
1021  continue;
1022  } else {
1023  break;
1024  }
1025  } else {
1026  break;
1027  }
1028  } while ( 1 ); // end restart loop
1029 
1030  // compute plastic strain vector
1031  helpVector.resize(strSize);
1032  helpVector.zero();
1033  for ( i = 1; i <= this->nsurf; i++ ) {
1034  if ( activeConditionMap.at(i) ) {
1035  ( * loadGradSigVecPtr ) [ i - 1 ].times( dgamma.at(i) );
1036  for ( j = 1; j <= strSize; j++ ) {
1037  helpVector.at(j) += ( * loadGradSigVecPtr ) [ i - 1 ].at(j);
1038  }
1039  }
1040  }
1041 
1042  plasticStrainVectorR.add(helpVector);
1043 
1044  // update plastic multipliers
1045  gamma.add(dgamma);
1046 
1047  elasticStrainVectorR = totalStrain;
1048  elasticStrainVectorR.subtract(plasticStrainVectorR);
1049  // stress vector in full form due to computational convinience
1050  this->computeTrialStressIncrement(fullStressVector, gp, elasticStrainVectorR, tStep);
1051 
1052  if ( hasHardening ) {
1053  this->computeStrainHardeningVarsIncrement(dkappa, gp, fullStressVector, gamma, helpVector, activeConditionMap);
1054  strainSpaceHardeningVariables = status->giveStrainSpaceHardeningVars();
1055  strainSpaceHardeningVariables.add(dkappa);
1056  }
1057 
1058  // check convergence
1059  yieldConsistency = true;
1060  for ( i = 1; i <= nsurf; i++ ) {
1061  if ( activeConditionMap.at(i) ) {
1062  yieldValue = this->computeYieldValueAt(gp, i, fullStressVector, strainSpaceHardeningVariables);
1063  if ( fabs(yieldValue) >= YIELD_TOL ) {
1064  yieldConsistency = false;
1065  }
1066  }
1067  }
1068 
1069  if ( yieldConsistency ) {
1070  restart = 0;
1071  // check for active conditions
1072  for ( i = 1; i <= this->nsurf; i++ ) {
1073  yieldValue = this->computeYieldValueAt(gp, i, fullStressVector, strainSpaceHardeningVariables);
1074  if ( ( gamma.at(i) < 0.0 ) || ( ( activeConditionMap.at(i) == 0 ) && ( yieldValue > YIELD_TOL ) ) ) {
1075  restart = 1;
1076  //printf ("*");
1077  }
1078  }
1079 
1080  if ( !restart ) {
1081  break;
1082  } else {
1083  IntArray newmap(nsurf);
1084  actSurf = 0;
1085  newmap.zero();
1086  // buid new actsurf mask
1087  for ( i = 1; i <= this->nsurf; i++ ) {
1088  yieldValue = this->computeYieldValueAt(gp, i, fullStressVector, strainSpaceHardeningVariables);
1089  if ( ( ( activeConditionMap.at(i) ) && ( gamma.at(i) >= 0.0 ) ) || ( yieldValue > YIELD_TOL ) ) {
1090  elastic = 0;
1091  actSurf += 1;
1092  newmap.at(i) = actSurf;
1093  if ( actSurf == this->giveMaxNumberOfActiveYieldConds(gp) ) {
1094  break;
1095  }
1096  }
1097  }
1098 
1099  activeConditionMap = newmap;
1100 
1101  if ( actSurf ) {
1102  activeConditionMap = newmap;
1103  } else {
1104  restart = 0;
1105  // select any possible
1106  for ( i = 1; i <= this->nsurf; i++ ) {
1107  if ( initialConditionMap.at(i) ) {
1108  actSurf = 1;
1109  activeConditionMap.zero();
1110  activeConditionMap.at(i) = 1;
1111  initialConditionMap.at(i) = 0;
1112  restart = 1;
1113  break;
1114  }
1115  }
1116  }
1117 
1118 
1119 
1120 
1121  plasticStrainVectorR = status->givePlasticStrainVector();
1122  strainSpaceHardeningVariables = status->giveStrainSpaceHardeningVars();
1123  elasticStrainVectorR = totalStrain;
1124  elasticStrainVectorR.subtract(plasticStrainVectorR);
1125  this->computeTrialStressIncrement(fullStressVector, gp, elasticStrainVectorR, tStep);
1126  gamma.zero();
1127  nIterations = 0;
1128  int restart = 0;
1129  if ( restart ) {
1130  goto huhu;
1131  }
1132 
1133  continue;
1134  //int restart = 0;
1135  //if (restart) goto huhu;
1136  }
1137  }
1138 
1139  // increment iteration count
1140  nIterations++;
1141 
1142  if ( nIterations > PLASTIC_MATERIAL_MAX_ITERATIONS ) {
1143  // try last resort select single function from initial active set
1144  restart = 0;
1145  bool smartCandidate = false;
1146  int candidate = 0;
1147  IntArray newmap(nsurf);
1148  actSurf = 0;
1149  newmap.zero();
1150  // check for active conditions
1151  for ( i = 1; i <= this->nsurf; i++ ) {
1152  if ( activeConditionMap.at(i) ) {
1153  if ( gamma.at(i) < 0.0 ) {
1154  smartCandidate = true;
1155  } else {
1156  candidate = i;
1157  }
1158  }
1159  }
1160 
1161  if ( smartCandidate && candidate && initialConditionMap.at(candidate) ) {
1162  actSurf = 1;
1163  activeConditionMap.zero();
1164  activeConditionMap.at(candidate) = 1;
1165  initialConditionMap.at(candidate) = 0;
1166  restart = 1;
1167  } else {
1168  // select any possible
1169  for ( i = 1; i <= this->nsurf; i++ ) {
1170  if ( initialConditionMap.at(i) ) {
1171  actSurf = 1;
1172  activeConditionMap.zero();
1173  activeConditionMap.at(i) = 1;
1174  initialConditionMap.at(i) = 0;
1175  restart = 1;
1176  break;
1177  }
1178  }
1179  }
1180 
1181  if ( restart ) {
1182  //fprintf(stderr,"===>LAST RESORT<====");
1183 
1184  plasticStrainVectorR = status->givePlasticStrainVector();
1185  strainSpaceHardeningVariables = status->giveStrainSpaceHardeningVars();
1186  elasticStrainVectorR = totalStrain;
1187  elasticStrainVectorR.subtract(plasticStrainVectorR);
1188  this->computeTrialStressIncrement(fullStressVector, gp, elasticStrainVectorR, tStep);
1189  gamma.zero();
1190  nIterations = 0;
1191  int restart = 0;
1192  if ( restart ) {
1193  goto huhu;
1194  }
1195 
1196  continue;
1197  } else {
1198  OOFEM_WARNING("local equlibrium not reached in %d iterations\nElement %d, gp %d, continuing",
1199  PLASTIC_MATERIAL_MAX_ITERATIONS, gp->giveElement()->giveNumber(), gp->giveNumber() );
1200  answer = fullStressVector;
1201  // debug line
1202  nIterations = 0;
1203  goto huhu;
1204  }
1205  }
1206  } while ( 1 );
1207  }
1208 
1209  //printf (" (%d iterations)", nIterations);
1210  answer = fullStressVector;
1211 }
1212 
1213 
1214 
1215 void
1216 MPlasticMaterial2 :: computeReducedStressGradientVector(FloatArray &answer, functType ftype, int isurf,
1217  GaussPoint *gp, const FloatArray &stressVector,
1218  const FloatArray &strainSpaceHardeningVariables)
1219 {
1220  FloatArray stressGradient;
1221 
1222  this->computeStressGradientVector(stressGradient, ftype, isurf, gp, stressVector,
1223  strainSpaceHardeningVariables);
1224 
1225  StructuralMaterial :: giveReducedSymVectorForm( answer, stressGradient, gp->giveMaterialMode() );
1226 }
1227 
1228 void
1229 MPlasticMaterial2 :: computeResidualVector(FloatArray &answer, GaussPoint *gp, const FloatArray &gamma,
1230  const IntArray &activeConditionMap, const FloatArray &plasticStrainVectorR,
1231  const FloatArray &strainSpaceHardeningVariables,
1232  std :: vector< FloatArray > &gradientVectorR)
1233 {
1234  /* Computes Residual vector for closes point projection algorithm */
1235 
1236  FloatArray oldPlasticStrainVectorR;
1237  int size;
1238  MPlasticMaterial2Status *status = static_cast< MPlasticMaterial2Status * >( this->giveStatus(gp) );
1239 
1240  size = plasticStrainVectorR.giveSize();
1241 
1242  answer.resize(size);
1243  oldPlasticStrainVectorR = status->givePlasticStrainVector();
1244 
1245  for ( int i = 1; i <= size; i++ ) {
1246  answer.at(i) = oldPlasticStrainVectorR.at(i) - plasticStrainVectorR.at(i);
1247  }
1248 
1249  for ( int i = 1; i <= this->nsurf; i++ ) {
1250  if ( activeConditionMap.at(i) ) {
1251  for ( int j = 1; j <= size; j++ ) {
1252  answer.at(j) += gamma.at(i) * gradientVectorR [ i - 1 ].at(j);
1253  }
1254  }
1255  }
1256 }
1257 
1258 
1259 void
1260 MPlasticMaterial2 :: computeTrialStressIncrement(FloatArray &answer, GaussPoint *gp,
1261  const FloatArray &elasticStrainVectorR,
1262  TimeStep *tStep)
1263 {
1264  /* Computes the full trial elastic stress vector */
1265 
1266  FloatMatrix de;
1267  FloatArray reducedAnswer;
1268 
1269  this->computeReducedElasticModuli(de, gp, tStep);
1270  //this->giveLinearElasticMaterial()->giveCharacteristicMatrix(de, TangentStiffness, gp, tStep);
1271  reducedAnswer.beProductOf(de, elasticStrainVectorR);
1272  StructuralMaterial :: giveFullSymVectorForm( answer, reducedAnswer, gp->giveMaterialMode() );
1273 }
1274 
1275 
1276 void
1277 MPlasticMaterial2 :: computeAlgorithmicModuli(FloatMatrix &answer,
1278  GaussPoint *gp,
1279  const FloatMatrix &elasticModuliInverse,
1280  const FloatArray &gamma,
1281  const IntArray &activeConditionMap,
1282  const FloatArray &fullStressVector,
1283  const FloatArray &strainSpaceHardeningVariables)
1284 {
1285  /* returns consistent moduli in reduced form.
1286  * stressVector in full form
1287  */
1288  FloatMatrix gradientMatrix, ks;
1289  FloatMatrix help, helpInverse;
1290  int i, j, size;
1291  bool hasHardening = this->hasHardening() > 0;
1292 
1293  size = elasticModuliInverse.giveNumberOfRows();
1294 
1295  // assemble consistent moduli
1296  helpInverse.resize(size, size);
1297  helpInverse.zero();
1298  for ( i = 1; i <= nsurf; i++ ) {
1299  if ( activeConditionMap.at(i) ) {
1300  // ask gradientMatrix
1301  this->computeReducedSSGradientMatrix(gradientMatrix, i, gp, fullStressVector,
1302  strainSpaceHardeningVariables);
1303 
1304  gradientMatrix.times( gamma.at(i) );
1305  helpInverse.add(gradientMatrix);
1306 
1307  if ( hasHardening ) {
1308  this->computeReducedSKGradientMatrix(gradientMatrix, i, gp, fullStressVector,
1309  strainSpaceHardeningVariables);
1310 
1311  gradientMatrix.times( gamma.at(i) );
1312  this->computeReducedHardeningVarsSigmaGradient(ks, gp, activeConditionMap, fullStressVector, strainSpaceHardeningVariables, gamma);
1313  help.beProductOf(gradientMatrix, ks);
1314  helpInverse.add(help);
1315  }
1316  }
1317  }
1318 
1319  for ( i = 1; i <= size; i++ ) {
1320  for ( j = 1; j <= size; j++ ) {
1321  helpInverse.at(i, j) += elasticModuliInverse.at(i, j);
1322  }
1323  }
1324 
1325  // prevent too big numbers
1326  for ( i = 1; i <= size; i++ ) {
1327  if ( fabs( helpInverse.at(i, i) ) > 1.e16 ) {
1328  helpInverse.at(i, i) = sgn( helpInverse.at(i, i) ) * 1.e16;
1329  }
1330  }
1331 
1332  answer.beInverseOf(helpInverse);
1333 }
1334 
1335 // ----------------------------------------------------------------------------//
1336 
1337 
1338 void
1339 MPlasticMaterial2 :: giveConsistentStiffnessMatrix(FloatMatrix &answer,
1340  MatResponseMode mode,
1341  GaussPoint *gp,
1342  TimeStep *tStep)
1343 {
1344  //
1345  // returns receiver material matrix for given reached state
1346  // described by temp variables.
1347  //
1348 
1349  FloatMatrix consistentModuli, elasticModuli, umat, vmat;
1350  FloatMatrix elasticModuliInverse, kl, ks;
1351  FloatMatrix gradientMatrix, gmat, gmatInv, gradMat, helpMtrx, helpMtrx2, answerR;
1352  FloatArray gradientVector, stressVector, fullStressVector;
1353  FloatArray strainSpaceHardeningVariables, helpVector;
1354  std :: vector< FloatArray > yieldGradSigVec(this->nsurf), loadGradSigVec(this->nsurf), * loadGradSigVecPtr;
1355  std :: vector< FloatArray > yieldGradKVec(this->nsurf), loadGradKVec(this->nsurf);
1356  FloatArray helpVector2;
1357 
1358  IntArray activeConditionMap, mask;
1359  FloatArray gamma;
1360  int strSize, indx, actSurf = 0;
1361  bool hasHardening = this->hasHardening() > 0;
1362 
1363  MPlasticMaterial2Status *status = static_cast< MPlasticMaterial2Status * >( this->giveStatus(gp) );
1364 
1365  if ( this->plType == associatedPT ) {
1366  loadGradSigVecPtr = & yieldGradSigVec;
1367  } else {
1368  loadGradSigVecPtr = & loadGradSigVec;
1369  }
1370 
1371  // ask for plastic consistency parameter
1372  gamma = status->giveTempGamma();
1373  activeConditionMap = status->giveTempActiveConditionMap();
1374  //
1375  // check for elastic cases
1376  //
1377  if ( ( status->giveTempStateFlag() == MPlasticMaterial2Status :: PM_Elastic ) ||
1378  ( status->giveTempStateFlag() == MPlasticMaterial2Status :: PM_Unloading ) ) {
1379  this->giveStiffnessMatrix(answer, ElasticStiffness, gp, tStep);
1380  return;
1381  }
1382 
1383  //
1384  // plastic case
1385  //
1386  // determine number of active surfaces
1387  for ( int i = 1; i <= nsurf; i++ ) {
1388  if ( activeConditionMap.at(i) ) {
1389  actSurf++;
1390  }
1391  }
1392 
1393  // compute elastic moduli and its inverse
1394  this->computeReducedElasticModuli(elasticModuli, gp, tStep);
1395  elasticModuliInverse.beInverseOf(elasticModuli);
1396  strSize = elasticModuliInverse.giveNumberOfRows();
1397 
1398  stressVector = status->giveTempStressVector();
1399  StructuralMaterial :: giveFullSymVectorForm( fullStressVector, stressVector, gp->giveMaterialMode() );
1400  strainSpaceHardeningVariables = status->giveTempStrainSpaceHardeningVarsVector();
1401 
1402  //
1403  // compute consistent moduli
1404  //
1405  this->computeAlgorithmicModuli(consistentModuli, gp, elasticModuliInverse, gamma,
1406  activeConditionMap, fullStressVector, strainSpaceHardeningVariables);
1407 
1408  //computee gmatInv
1409  for ( int i = 1; i <= nsurf; i++ ) {
1410  this->computeReducedStressGradientVector(yieldGradSigVec [ i - 1 ], yieldFunction, i, gp, fullStressVector,
1411  strainSpaceHardeningVariables);
1412  if ( hasHardening ) {
1413  this->computeKGradientVector(yieldGradKVec [ i - 1 ], yieldFunction, i, gp, fullStressVector,
1414  strainSpaceHardeningVariables);
1415  }
1416  }
1417 
1418  if ( this->plType == nonassociatedPT ) {
1419  for ( int i = 1; i <= nsurf; i++ ) {
1420  this->computeReducedStressGradientVector(loadGradSigVec [ i - 1 ], loadFunction, i, gp, fullStressVector,
1421  strainSpaceHardeningVariables);
1422  if ( hasHardening ) {
1423  this->computeKGradientVector(loadGradKVec [ i - 1 ], loadFunction, i, gp, fullStressVector,
1424  strainSpaceHardeningVariables);
1425  }
1426  }
1427  }
1428 
1429 
1430  umat.resize(strSize, actSurf);
1431  umat.zero();
1432  vmat.resize(actSurf, strSize);
1433  vmat.zero();
1434  gmat.resize(actSurf, actSurf);
1435  gmat.zero();
1436 
1437  if ( hasHardening ) {
1438  this->computeReducedHardeningVarsSigmaGradient(ks, gp, activeConditionMap, fullStressVector,
1439  strainSpaceHardeningVariables, gamma);
1440  this->computeReducedHardeningVarsLamGradient(kl, gp, actSurf, activeConditionMap,
1441  fullStressVector, strainSpaceHardeningVariables, gamma);
1442  }
1443 
1444  for ( int i = 1; i <= nsurf; i++ ) {
1445  if ( ( indx = activeConditionMap.at(i) ) ) {
1446  if ( hasHardening ) {
1447  helpVector.beTProductOf(ks, yieldGradKVec [ i - 1 ]);
1448  helpVector.add(yieldGradSigVec [ i - 1 ]);
1449  for ( int j = 1; j <= strSize; j++ ) {
1450  vmat.at(indx, j) = helpVector.at(j);
1451  }
1452  } else {
1453  for ( int j = 1; j <= strSize; j++ ) {
1454  vmat.at(indx, j) = yieldGradSigVec [ i - 1 ].at(j);
1455  }
1456  }
1457 
1458  if ( hasHardening ) {
1459  this->computeReducedSKGradientMatrix(gradientMatrix, i, gp, fullStressVector,
1460  strainSpaceHardeningVariables);
1461  helpMtrx.beProductOf(gradientMatrix, kl);
1462  helpMtrx.times( gamma.at(i) );
1463  umat.add(helpMtrx);
1464  }
1465 
1466  for ( int j = 1; j <= strSize; j++ ) {
1467  umat.at(j, indx) += ( ( * loadGradSigVecPtr ) [ i - 1 ] ).at(j);
1468  }
1469 
1470 
1471  if ( hasHardening ) {
1472  helpVector2.beTProductOf(kl, yieldGradKVec [ i - 1 ]);
1473  for ( int j = 1; j <= actSurf; j++ ) {
1474  gmat.at(indx, j) = ( -1.0 ) * helpVector2.at(j);
1475  }
1476  }
1477  }
1478  }
1479 
1480  helpMtrx.beProductOf(vmat, consistentModuli); // V S
1481  helpMtrx2.beProductOf(helpMtrx, umat);
1482  gmat.add(helpMtrx2);
1484  gmatInv.beInverseOf(gmat);
1485 
1486 
1487  helpMtrx2.beProductOf(gmatInv, helpMtrx);
1488  helpMtrx.beProductOf(consistentModuli, umat); // S U
1489  answer.beProductOf(helpMtrx, helpMtrx2); // SUGVS
1490  answer.negated();
1491  answer.add(consistentModuli);
1492 }
1493 
1494 void
1495 MPlasticMaterial2 :: giveElastoPlasticStiffnessMatrix(FloatMatrix &answer,
1496  MatResponseMode mode,
1497  GaussPoint *gp,
1498  TimeStep *tStep)
1499 {
1500  //
1501  // returns receiver material matrix for given reached state
1502  // described by temp variables.
1503  //
1504 
1505  FloatMatrix elasticModuli, umat, vmat;
1506  FloatMatrix gmat, gmatInv, helpMtrx, helpMtrx2, kl, ks;
1507  FloatArray stressVector, fullStressVector;
1508  FloatArray strainSpaceHardeningVariables, helpVector, helpVector2;
1509  std :: vector< FloatArray > yieldGradSigVec(this->nsurf), loadGradSigVec(this->nsurf), * loadGradSigVecPtr;
1510  std :: vector< FloatArray > yieldGradKVec(this->nsurf), loadGradKVec(this->nsurf);
1511  FloatArray helpVec;
1512 
1513  IntArray activeConditionMap, mask;
1514  FloatArray gamma;
1515  int strSize, indx, actSurf = 0;
1516  bool hasHardening = this->hasHardening() > 0;
1517 
1518  MPlasticMaterial2Status *status = static_cast< MPlasticMaterial2Status * >( this->giveStatus(gp) );
1519 
1520  if ( this->plType == associatedPT ) {
1521  loadGradSigVecPtr = & yieldGradSigVec;
1522  } else {
1523  loadGradSigVecPtr = & loadGradSigVec;
1524  }
1525 
1526  // ask for plastic consistency parameter
1527  gamma = status->giveTempGamma();
1528  activeConditionMap = status->giveTempActiveConditionMap();
1529  //
1530  // check for elastic cases
1531  //
1532  if ( ( status->giveTempStateFlag() == MPlasticMaterial2Status :: PM_Elastic ) ||
1533  ( status->giveTempStateFlag() == MPlasticMaterial2Status :: PM_Unloading ) ) {
1534  this->giveStiffnessMatrix(answer, ElasticStiffness, gp, tStep);
1535  return;
1536  }
1537 
1538  //
1539  // plastic case
1540  //
1541  // determine number of active surfaces
1542  for ( int i = 1; i <= nsurf; i++ ) {
1543  if ( activeConditionMap.at(i) ) {
1544  actSurf++;
1545  }
1546  }
1547 
1548  // compute elastic moduli and its inverse
1549  this->computeReducedElasticModuli(elasticModuli, gp, tStep);
1550 
1551 
1552  stressVector = status->giveStressVector();
1553  StructuralMaterial :: giveFullSymVectorForm( fullStressVector, stressVector, gp->giveMaterialMode() );
1554  strainSpaceHardeningVariables = status->giveStrainSpaceHardeningVars();
1555  strSize = elasticModuli.giveNumberOfRows();
1556 
1557  //compute gmatInv
1558  for ( int i = 1; i <= nsurf; i++ ) {
1559  this->computeReducedStressGradientVector(yieldGradSigVec [ i - 1 ], yieldFunction, i, gp, fullStressVector,
1560  strainSpaceHardeningVariables);
1561  if ( hasHardening ) {
1562  this->computeKGradientVector(yieldGradKVec [ i - 1 ], yieldFunction, i, gp, fullStressVector,
1563  strainSpaceHardeningVariables);
1564  }
1565  }
1566 
1567  if ( this->plType == nonassociatedPT ) {
1568  for ( int i = 1; i <= nsurf; i++ ) {
1569  this->computeReducedStressGradientVector(loadGradSigVec [ i - 1 ], loadFunction, i, gp, fullStressVector,
1570  strainSpaceHardeningVariables);
1571  if ( hasHardening ) {
1572  this->computeKGradientVector(loadGradKVec [ i - 1 ], loadFunction, i, gp, fullStressVector,
1573  strainSpaceHardeningVariables);
1574  }
1575  }
1576  }
1577 
1578  umat.resize(strSize, actSurf);
1579  vmat.resize(actSurf, strSize);
1580  gmat.resize(actSurf, actSurf);
1581 
1582  if ( hasHardening ) {
1583  this->computeReducedHardeningVarsSigmaGradient(ks, gp, activeConditionMap, fullStressVector,
1584  strainSpaceHardeningVariables, gamma);
1585  this->computeReducedHardeningVarsLamGradient(kl, gp, actSurf, activeConditionMap,
1586  fullStressVector, strainSpaceHardeningVariables, gamma);
1587  }
1588 
1589  for ( int i = 1; i <= nsurf; i++ ) {
1590  if ( ( indx = activeConditionMap.at(i) ) ) {
1591  if ( hasHardening ) {
1592  helpVector.beTProductOf(ks, yieldGradKVec [ i - 1 ]);
1593  helpVector.add(yieldGradSigVec [ i - 1 ]);
1594  for ( int j = 1; j <= strSize; j++ ) {
1595  vmat.at(indx, j) = helpVector.at(j);
1596  }
1597  } else {
1598  for ( int j = 1; j <= strSize; j++ ) {
1599  vmat.at(indx, j) = yieldGradSigVec [ i - 1 ].at(j);
1600  }
1601  }
1602 
1603  for ( int j = 1; j <= strSize; j++ ) {
1604  umat.at(j, indx) += ( * loadGradSigVecPtr ) [ i - 1 ].at(j);
1605  }
1606 
1607 
1608  if ( hasHardening ) {
1609  helpVector2.beTProductOf(kl, yieldGradKVec [ i - 1 ]);
1610  for ( int j = 1; j <= actSurf; j++ ) {
1611  gmat.at(indx, j) = ( -1.0 ) * helpVector2.at(j);
1612  }
1613  }
1614  }
1615  }
1616 
1617  helpMtrx.beProductOf(vmat, elasticModuli); // V S
1618  helpMtrx2.beProductOf(helpMtrx, umat);
1619  gmat.add(helpMtrx2);
1621  gmatInv.beInverseOf(gmat);
1622 
1623  helpMtrx.beProductOf(gmatInv, vmat);
1624  helpMtrx2.beProductOf(umat, helpMtrx);
1625  answer.beProductOf(elasticModuli, helpMtrx2);
1626  answer.negated();
1627  answer.add(elasticModuli);
1628 }
1629 
1630 /*
1631  * void
1632  * MPlasticMaterial2 :: computeDiagModuli(FloatMatrix& answer,
1633  * GaussPoint *gp, FloatMatrix &elasticModuliInverse,
1634  * FloatMatrix &hardeningModuliInverse)
1635  * {
1636  * //
1637  * // assembles diagonal moduli from elasticModuliInverse and hardeningModuliInverse
1638  * //
1639  * int size1, size2, i,j;
1640  *
1641  * size1 = elasticModuliInverse.giveNumberOfRows();
1642  * if (hardeningModuliInverse.giveNumberOfRows()) size2=size1+hardeningModuliInverse.giveNumberOfRows();
1643  * else size2 = size1;
1644  *
1645  * answer.resize (size2, size2);
1646  * answer.zero();
1647  *
1648  * for (i=1; i<=size1; i++)
1649  * for (j=1; j<= size1; j++) answer.at(i,j) = elasticModuliInverse.at(i,j);
1650  *
1651  * for (i=size1+1; i<= size2; i++)
1652  * for (j=size1+1; j<= size2; j++) answer.at(i,j) = hardeningModuliInverse.at(i-size1, j-size1);
1653  * }
1654  */
1655 
1656 void
1657 MPlasticMaterial2 :: computeReducedElasticModuli(FloatMatrix &answer,
1658  GaussPoint *gp,
1659  TimeStep *tStep)
1660 { /* Returns elastic moduli in reduced stress-strain space*/
1661  this->giveLinearElasticMaterial()->giveStiffnessMatrix(answer, ElasticStiffness, gp, tStep);
1662 }
1663 
1664 
1665 // overloaded from structural material
1666 
1667 void
1668 MPlasticMaterial2 :: give3dMaterialStiffnessMatrix(FloatMatrix &answer,
1669  MatResponseMode mode,
1670  GaussPoint *gp,
1671  TimeStep *tStep)
1672 //
1673 //
1674 //
1675 // computes full 3d constitutive matrix for case of 3d stress-strain state.
1676 // it returns elasto-plastic stiffness material matrix.
1677 // if strainIncrement == NULL a loading is assumed
1678 // for detailed description see (W.F.Chen: Plasticity in Reinforced Concrete, McGraw-Hill, 1982,
1679 // chapter 6.)
1680 //
1681 // if derived material would like to implement failure behaviour
1682 // it must redefine basic Give3dMaterialStiffnessMatrix function
1683 // in order to take possible failure (tension cracking) into account
1684 //
1685 //
1686 //
1687 {
1688  MaterialMode originalMode = gp->giveMaterialMode();
1689  if ( originalMode != _3dMat ) {
1690  OOFEM_ERROR("Different stressStrain mode encountered");
1691  }
1692 
1693  // we can force 3d response, and we obtain correct 3d tangent matrix,
1694  // but in fact, stress integration algorithm will not work
1695  // because in stress integration algorithm we are unable to recognize
1696  // which reduction from 3d case should be performed to obtain correct result.
1697  // so for new stressStrain state, instead of programming 3d reduction,
1698  // you should enhance imposeConstraints functions for ne state, and
1699  // then programming simple inteface function for you stressstrain state
1700  // calling GiveMaterailStiffenssMatrix, which imposes constrains correctly.
1701  if ( mode == TangentStiffness ) {
1702  if ( rmType == mpm_ClosestPoint ) {
1703  this->giveConsistentStiffnessMatrix(answer, mode, gp, tStep);
1704  } else {
1705  this->giveElastoPlasticStiffnessMatrix(answer, mode, gp, tStep);
1706  }
1707  } else {
1708  this->giveLinearElasticMaterial()->give3dMaterialStiffnessMatrix(answer, mode, gp, tStep);
1709  }
1710 }
1711 
1712 
1713 void
1714 MPlasticMaterial2 :: givePlaneStressStiffMtrx(FloatMatrix &answer,
1715  MatResponseMode mode,
1716  GaussPoint *gp,
1717  TimeStep *tStep)
1718 
1719 //
1720 // returns receiver's 2dPlaneStressMtrx
1721 // (2dPlaneStres ==> sigma_z = tau_xz = tau_yz = 0.)
1722 //
1723 // standard method from Material Class overloaded, because no inversion is needed.
1724 // the reduction from 3d case will not work
1725 // this implementation should be faster.
1726 {
1727  if ( mode == TangentStiffness ) {
1728  if ( rmType == mpm_ClosestPoint ) {
1729  this->giveConsistentStiffnessMatrix(answer, mode, gp, tStep);
1730  } else {
1731  this->giveElastoPlasticStiffnessMatrix(answer, mode, gp, tStep);
1732  }
1733  } else {
1734  this->giveLinearElasticMaterial()->givePlaneStressStiffMtrx(answer, mode, gp, tStep);
1735  }
1736 }
1737 
1738 
1739 void
1740 MPlasticMaterial2 :: givePlaneStrainStiffMtrx(FloatMatrix &answer,
1741  MatResponseMode mode,
1742  GaussPoint *gp,
1743  TimeStep *tStep)
1744 
1745 //
1746 // return receiver's 2dPlaneStrainMtrx constructed from
1747 // general 3dMatrialStiffnessMatrix
1748 // (2dPlaneStrain ==> eps_z = gamma_xz = gamma_yz = 0.)
1749 //
1750 {
1751  if ( mode == TangentStiffness ) {
1752  if ( rmType == mpm_ClosestPoint ) {
1753  this->giveConsistentStiffnessMatrix(answer, mode, gp, tStep);
1754  } else {
1755  this->giveElastoPlasticStiffnessMatrix(answer, mode, gp, tStep);
1756  }
1757  } else {
1758  this->giveLinearElasticMaterial()->givePlaneStrainStiffMtrx(answer, mode, gp, tStep);
1759  }
1760 }
1761 
1762 
1763 void
1764 MPlasticMaterial2 :: give1dStressStiffMtrx(FloatMatrix &answer,
1765  MatResponseMode mode,
1766  GaussPoint *gp,
1767  TimeStep *tStep)
1768 
1769 //
1770 // returns receiver's 1dMaterialStiffnessMAtrix
1771 // (1d case ==> sigma_y = sigma_z = tau_yz = tau_zx = tau_xy = 0.)
1772 {
1773  if ( mode == TangentStiffness ) {
1774  if ( rmType == mpm_ClosestPoint ) {
1775  this->giveConsistentStiffnessMatrix(answer, mode, gp, tStep);
1776  } else {
1777  this->giveElastoPlasticStiffnessMatrix(answer, mode, gp, tStep);
1778  }
1779  } else {
1780  this->giveLinearElasticMaterial()->give1dStressStiffMtrx(answer, mode, gp, tStep);
1781  }
1782 }
1783 
1784 
1785 
1786 void
1787 MPlasticMaterial2 :: give2dBeamLayerStiffMtrx(FloatMatrix &answer,
1788  MatResponseMode mode,
1789  GaussPoint *gp,
1790  TimeStep *tStep)
1791 //
1792 // returns receiver's 2dBeamLayerStiffMtrx.
1793 // (2dPlaneStres ==> sigma_z = tau_xz = tau_yz = 0.)
1794 //
1795 // standard method from Material Class overloaded, because no inversion is needed.
1796 // the reduction from 3d case will not work
1797 // this implementation should be faster.
1798 {
1799  if ( mode == TangentStiffness ) {
1800  if ( rmType == mpm_ClosestPoint ) {
1801  this->giveConsistentStiffnessMatrix(answer, mode, gp, tStep);
1802  } else {
1803  this->giveElastoPlasticStiffnessMatrix(answer, mode, gp, tStep);
1804  }
1805  } else {
1806  this->giveLinearElasticMaterial()->give2dBeamLayerStiffMtrx(answer, mode, gp, tStep);
1807  }
1808 }
1809 
1810 
1811 
1812 void
1813 MPlasticMaterial2 :: givePlateLayerStiffMtrx(FloatMatrix &answer,
1814  MatResponseMode mode,
1815  GaussPoint *gp,
1816  TimeStep *tStep)
1817 //
1818 // returns receiver's 2dPlateLayerMtrx
1819 // (2dPlaneStres ==> sigma_z = tau_xz = tau_yz = 0.)
1820 //
1821 // standard method from Material Class overloaded, because no inversion is needed.
1822 // the reduction from 3d case will not work
1823 // this implementation should be faster.
1824 {
1825  if ( mode == TangentStiffness ) {
1826  if ( rmType == mpm_ClosestPoint ) {
1827  this->giveConsistentStiffnessMatrix(answer, mode, gp, tStep);
1828  } else {
1829  this->giveElastoPlasticStiffnessMatrix(answer, mode, gp, tStep);
1830  }
1831  } else {
1832  this->giveLinearElasticMaterial()->givePlateLayerStiffMtrx(answer, mode, gp, tStep);
1833  }
1834 }
1835 
1836 void
1837 MPlasticMaterial2 :: giveFiberStiffMtrx(FloatMatrix &answer,
1838  MatResponseMode mode,
1839  GaussPoint *gp,
1840  TimeStep *tStep)
1841 //
1842 // returns receiver's 1dFiber
1843 // (1dFiber ==> sigma_y = sigma_z = tau_yz = 0.)
1844 //
1845 // standard method from Material Class overloaded, because no inversion is needed.
1846 // the reduction from 3d case will not work
1847 // this implementation should be faster.
1848 {
1849  if ( mode == TangentStiffness ) {
1850  if ( rmType == mpm_ClosestPoint ) {
1851  this->giveConsistentStiffnessMatrix(answer, mode, gp, tStep);
1852  } else {
1853  this->giveElastoPlasticStiffnessMatrix(answer, mode, gp, tStep);
1854  }
1855  } else {
1856  this->giveLinearElasticMaterial()->giveFiberStiffMtrx(answer, mode, gp, tStep);
1857  }
1858 }
1859 
1860 
1861 int
1862 MPlasticMaterial2 :: giveIPValue(FloatArray &answer, GaussPoint *gp, InternalStateType type, TimeStep *tStep)
1863 {
1864  MPlasticMaterial2Status *status = static_cast< MPlasticMaterial2Status * >( this->giveStatus(gp) );
1865  if ( type == IST_PlasticStrainTensor ) {
1866  const FloatArray &ep = status->givePlasticStrainVector();
1868  StructuralMaterial :: giveFullSymVectorForm( answer, ep, gp->giveMaterialMode() );
1869  return 1;
1870  } else if ( type == IST_PrincipalPlasticStrainTensor ) {
1871  FloatArray st;
1872 
1873  const FloatArray &s = status->givePlasticStrainVector();
1874 
1876  StructuralMaterial :: giveFullSymVectorForm( st, s, gp->giveMaterialMode() );
1877 
1878  this->computePrincipalValues(answer, st, principal_strain);
1879  return 1;
1880  } else {
1881  return StructuralMaterial :: giveIPValue(answer, gp, type, tStep);
1882  }
1883 }
1884 
1885 long
1886 MPlasticMaterial2 :: getPopulationSignature(IntArray &mask)
1887 {
1888  long val = 0;
1889  int size = mask.giveSize();
1890  if ( size > ( int ) ( sizeof( long ) * 8 ) ) {
1891  printf("MPlasticMaterial2: too many yield conditions");
1892  }
1893 
1894  for ( int i = 1; i <= mask.giveSize(); i++ ) {
1895  if ( mask.at(i) ) {
1896  val |= ( 1L << ( i - 1 ) );
1897  }
1898  }
1899 
1900  return val;
1901 }
1902 
1903 int
1904 MPlasticMaterial2 :: testPopulation(long pop)
1905 {
1906  if ( populationSet.find(pop) != populationSet.end() ) {
1907  return 0;
1908  } else {
1909  return 1;
1910  }
1911 }
1912 
1913 void
1914 MPlasticMaterial2 :: clearPopulationSet()
1915 {
1916  populationSet.clear();
1917 }
1918 
1919 void
1920 MPlasticMaterial2 :: addNewPopulation(IntArray &mask)
1921 {
1922  long val = getPopulationSignature(mask);
1923  populationSet.insert(val);
1924 }
1925 
1926 
1927 int
1928 MPlasticMaterial2 :: getNewPopulation(IntArray &result, IntArray &candidateMask, int degree, int size)
1929 {
1930  long val = 0;
1931  int candDegree = 0, candSize = candidateMask.giveSize();
1932  for ( int i = 1; i <= candSize; i++ ) {
1933  if ( candidateMask.at(i) ) {
1934  candDegree++;
1935  }
1936  }
1937 
1938  result.resize(size);
1939 
1940  // first try candidate if compatible
1941  if ( candDegree <= degree ) {
1942  // generate candite popul value
1943  val = getPopulationSignature(candidateMask);
1944 
1945  if ( testPopulation(val) ) {
1946  // candite is o.k.
1947  populationSet.insert(val);
1948  result = candidateMask;
1949  return 1;
1950  }
1951  } else {
1952  // try to generate suitable candidate from candidateMask if possible
1953  // first construct array of candidates
1954  IntArray candidateArray(candDegree);
1955  IntArray ind(degree);
1956  for ( int i = 1, j = 1; i <= candSize; i++ ) {
1957  if ( candidateMask.at(i) ) {
1958  candidateArray.at(j++) = i;
1959  }
1960  }
1961 
1962  // first try first degree candidates
1963  int activeIndx = degree;
1964  for ( int i = 1; i <= degree; i++ ) {
1965  ind.at(i) = i;
1966  }
1967 
1968 
1969  do {
1970  result.zero();
1971  int ii, count = 1;
1972  for ( int i = 1; i <= degree; i++ ) {
1973  ii = candidateArray.at( ind.at(i) );
1974  if ( !result.at(ii) ) {
1975  result.at(ii) = count++;
1976  }
1977  }
1978 
1979  // order result array
1980  count = 1;
1981  for ( int i = 1; i <= size; i++ ) {
1982  if ( result.at(i) ) {
1983  result.at(i) = count++;
1984  }
1985  }
1986 
1987  val = getPopulationSignature(result);
1988  if ( testPopulation(val) ) {
1989  // candite is o.k.
1990  populationSet.insert(val);
1991  return 1;
1992  }
1993 
1994  bool cont = false;
1995  activeIndx = degree;
1996  do {
1997  if ( ind.at(activeIndx) < candDegree ) {
1998  ind.at(activeIndx)++;
1999  cont = false;
2000  } else {
2001  ind.at(activeIndx) = 1;
2002  activeIndx--;
2003  if ( activeIndx > 0 ) {
2004  cont = true;
2005  } else {
2006  cont = false;
2007  }
2008  }
2009  } while ( cont );
2010  } while ( activeIndx >= 1 );
2011  }
2012 
2013  // generate any suitable population
2014  IntArray ind(degree);
2015  // first try first degree candidates
2016  int activeIndx = degree;
2017  for ( int i = 1; i <= degree; i++ ) {
2018  ind.at(i) = 1;
2019  }
2020 
2021 
2022  do {
2023  result.zero();
2024  int ii, count = 1;
2025  for ( int i = 1; i <= degree; i++ ) {
2026  ii = ind.at(i);
2027  if ( !result.at(ii) ) {
2028  result.at(ii) = count++;
2029  }
2030  }
2031 
2032  // order result array
2033  count = 1;
2034  for ( int i = 1; i <= size; i++ ) {
2035  if ( result.at(i) ) {
2036  result.at(i) = count++;
2037  }
2038  }
2039 
2040  val = getPopulationSignature(result);
2041  if ( testPopulation(val) ) {
2042  // candite is o.k.
2043  populationSet.insert(val);
2044  return 1;
2045  }
2046 
2047  bool cont = false;
2048  activeIndx = degree;
2049  do {
2050  if ( ind.at(activeIndx) < size ) {
2051  ind.at(activeIndx)++;
2052  cont = false;
2053  } else {
2054  ind.at(activeIndx) = 1;
2055  activeIndx--;
2056  if ( activeIndx > 0 ) {
2057  cont = true;
2058  } else {
2059  cont = false;
2060  }
2061  }
2062  } while ( cont );
2063  } while ( activeIndx >= 1 );
2064 
2065  // there is no population left -> return 0;
2066  return 0;
2067 }
2068 
2069 
2070 //
2071 // SmearedCrackingMaterialStatus Class
2072 //
2073 
2074 MPlasticMaterial2Status :: MPlasticMaterial2Status(int n, Domain *d, GaussPoint *g, int statusSize) :
2075  StructuralMaterialStatus(n, d, g), plasticStrainVector(), tempPlasticStrainVector(),
2076  strainSpaceHardeningVarsVector(statusSize), tempStrainSpaceHardeningVarsVector(statusSize),
2077  state_flag(MPlasticMaterial2Status :: PM_Elastic),
2078  temp_state_flag(MPlasticMaterial2Status :: PM_Elastic),
2079  damage(0.),
2080  tempDamage(0.),
2081  gamma(),
2082  tempGamma()
2083 { }
2084 
2085 MPlasticMaterial2Status :: ~MPlasticMaterial2Status()
2086 { }
2087 
2088 void
2089 MPlasticMaterial2Status :: printOutputAt(FILE *file, TimeStep *tStep)
2090 {
2091  StructuralMaterialStatus :: printOutputAt(file, tStep);
2092  fprintf(file, "status { ");
2093  if ( ( state_flag == MPlasticMaterial2Status :: PM_Yielding ) || ( state_flag == MPlasticMaterial2Status :: PM_Unloading ) ) {
2094  if ( state_flag == MPlasticMaterial2Status :: PM_Yielding ) {
2095  fprintf(file, " Yielding,");
2096  } else {
2097  fprintf(file, " Unloading,");
2098  }
2099 
2100  fprintf(file, " Plastic strains");
2101  for ( auto &val : plasticStrainVector ) {
2102  fprintf( file, " %.4e", val );
2103  }
2104 
2105  if ( strainSpaceHardeningVarsVector.giveSize() ) {
2106  fprintf(file, " Strain space hardening vars");
2107  for ( auto &val : strainSpaceHardeningVarsVector ) {
2108  fprintf( file, " %.4e", val );
2109  }
2110  }
2111 
2112  fprintf(file, " ActiveConditionMap");
2113  for ( auto &val : activeConditionMap ) {
2114  fprintf( file, " %d", val );
2115  }
2116  }
2117 
2118  fprintf(file, "}\n");
2119 }
2120 
2121 
2122 void MPlasticMaterial2Status :: initTempStatus()
2123 //
2124 // initializes temp variables according to variables form previous equlibrium state.
2125 //
2126 {
2127  StructuralMaterialStatus :: initTempStatus();
2128 
2129  if ( plasticStrainVector.giveSize() == 0 ) {
2130  plasticStrainVector.resize( StructuralMaterial :: giveSizeOfVoigtSymVector( gp->giveMaterialMode() ) );
2131  plasticStrainVector.zero();
2132  }
2133 
2134  tempPlasticStrainVector = plasticStrainVector;
2135 
2136  tempStrainSpaceHardeningVarsVector = strainSpaceHardeningVarsVector;
2137 
2138 
2139  tempGamma = gamma;
2140  tempActiveConditionMap = activeConditionMap;
2141 }
2142 
2143 
2144 void
2145 MPlasticMaterial2Status :: updateYourself(TimeStep *tStep)
2146 //
2147 // updates variables (nonTemp variables describing situation at previous equilibrium state)
2148 // after a new equilibrium state has been reached
2149 // temporary variables are having values corresponding to newly reached equilibrium.
2150 //
2151 {
2152  StructuralMaterialStatus :: updateYourself(tStep);
2153 
2154  plasticStrainVector = tempPlasticStrainVector;
2155  strainSpaceHardeningVarsVector = tempStrainSpaceHardeningVarsVector;
2156 
2157  state_flag = temp_state_flag;
2158  gamma = tempGamma;
2159  activeConditionMap = tempActiveConditionMap;
2160 }
2161 
2162 
2163 
2164 
2165 contextIOResultType
2166 MPlasticMaterial2Status :: saveContext(DataStream &stream, ContextMode mode, void *obj)
2167 //
2168 // saves full information stored in this Status
2169 // no temp variables stored
2170 //
2171 {
2172  contextIOResultType iores;
2173 
2174  // save parent class status
2175  if ( ( iores = StructuralMaterialStatus :: saveContext(stream, mode, obj) ) != CIO_OK ) {
2176  THROW_CIOERR(iores);
2177  }
2178 
2179  // write a raw data
2180  if ( ( iores = plasticStrainVector.storeYourself(stream) ) != CIO_OK ) {
2181  THROW_CIOERR(iores);
2182  }
2183 
2184  if ( ( iores = strainSpaceHardeningVarsVector.storeYourself(stream) ) != CIO_OK ) {
2185  THROW_CIOERR(iores);
2186  }
2187 
2188  if ( !stream.write(state_flag) ) {
2189  THROW_CIOERR(CIO_IOERR);
2190  }
2191 
2192  if ( ( iores = gamma.storeYourself(stream) ) != CIO_OK ) {
2193  THROW_CIOERR(iores);
2194  }
2195 
2196  if ( ( iores = activeConditionMap.storeYourself(stream) ) != CIO_OK ) {
2197  THROW_CIOERR(iores);
2198  }
2199 
2200 
2201  return CIO_OK;
2202 }
2203 
2204 
2205 
2206 contextIOResultType
2207 MPlasticMaterial2Status :: restoreContext(DataStream &stream, ContextMode mode, void *obj)
2208 //
2209 // restores full information stored in stream to this Status
2210 //
2211 {
2212  contextIOResultType iores;
2213 
2214  // read parent class status
2215  if ( ( iores = StructuralMaterialStatus :: restoreContext(stream, mode, obj) ) != CIO_OK ) {
2216  THROW_CIOERR(iores);
2217  }
2218 
2219  if ( ( iores = plasticStrainVector.restoreYourself(stream) ) != CIO_OK ) {
2220  THROW_CIOERR(iores);
2221  }
2222 
2223  if ( ( iores = strainSpaceHardeningVarsVector.restoreYourself(stream) ) != CIO_OK ) {
2224  THROW_CIOERR(iores);
2225  }
2226 
2227  if ( !stream.read(state_flag) ) {
2228  THROW_CIOERR(CIO_IOERR);
2229  }
2230 
2231  if ( ( iores = gamma.restoreYourself(stream) ) != CIO_OK ) {
2232  THROW_CIOERR(iores);
2233  }
2234 
2235  if ( ( iores = activeConditionMap.restoreYourself(stream) ) != CIO_OK ) {
2236  THROW_CIOERR(iores);
2237  }
2238 
2239 
2240  return CIO_OK; // return succes
2241 }
2242 } // end namespace oofem
oofem::StructuralMaterial::giveSizeOfVoigtSymVector
static int giveSizeOfVoigtSymVector(MaterialMode mmode)
Returns the size of symmetric part of a reduced stress/strain vector according to given mode...
Definition: structuralmaterial.C:801
oofem::MPlasticMaterial2::computeReducedSSGradientMatrix
virtual void computeReducedSSGradientMatrix(FloatMatrix &gradientMatrix, int i, GaussPoint *gp, const FloatArray &fullStressVector, const FloatArray &strainSpaceHardeningVariables)=0
Computes second derivative of loading function with respect to stress.
oofem::MPlasticMaterial2::givePlateLayerStiffMtrx
virtual void givePlateLayerStiffMtrx(FloatMatrix &answer, MatResponseMode, GaussPoint *gp, TimeStep *tStep)
Method for computing 2d plate layer stiffness matrix of receiver.
Definition: mplasticmaterial2.C:1813
YIELD_TOL
#define YIELD_TOL
Definition: mplasticmaterial2.C:45
oofem::StructuralMaterialStatus::updateYourself
virtual void updateYourself(TimeStep *tStep)
Update equilibrium history variables according to temp-variables.
Definition: structuralms.C:96
oofem::InternalStateType
InternalStateType
Type representing the physical meaning of element or constitutive model internal variable.
Definition: internalstatetype.h:187
oofem::StructuralMaterial::computePrincipalValues
static void computePrincipalValues(FloatArray &answer, const FloatArray &s, stressStrainPrincMode mode)
Common functions for convenience.
Definition: structuralmaterial.C:1264
oofem::GaussPoint::giveMaterialMode
MaterialMode giveMaterialMode()
Returns corresponding material mode of receiver.
Definition: gausspoint.h:191
oofem::IntArray::storeYourself
contextIOResultType storeYourself(DataStream &stream) const
Stores array to output stream.
Definition: intarray.C:289
oofem::FloatArray::subtract
void subtract(const FloatArray &src)
Subtracts array src to receiver.
Definition: floatarray.C:258
oofem::MPlasticMaterial2Status::MPlasticMaterial2Status
MPlasticMaterial2Status(int n, Domain *d, GaussPoint *g, int statusSize)
Definition: mplasticmaterial2.C:2074
oofem::MPlasticMaterial2Status::gamma
FloatArray gamma
Consistency parameter values (needed for algorithmic stiffness).
Definition: mplasticmaterial2.h:92
oofem::StructuralMaterialStatus::letTempStrainVectorBe
void letTempStrainVectorBe(const FloatArray &v)
Assigns tempStrainVector to given vector v.
Definition: structuralms.h:137
oofem::StructuralMaterial::give1dStressStiffMtrx
virtual void give1dStressStiffMtrx(FloatMatrix &answer, MatResponseMode mmode, GaussPoint *gp, TimeStep *tStep)
Method for computing 1d stiffness matrix of receiver.
Definition: structuralmaterial.C:1109
oofem::IntArray::printYourself
void printYourself() const
Prints receiver on stdout.
Definition: intarray.C:225
oofem::Material::giveStatus
virtual MaterialStatus * giveStatus(GaussPoint *gp) const
Returns material status of receiver in given integration point.
Definition: material.C:244
oofem::IntegrationPointStatus::gp
GaussPoint * gp
Associated integration point.
Definition: integrationpointstatus.h:57
oofem::Domain
Class and object Domain.
Definition: domain.h:115
oofem::MPlasticMaterial2::computeReducedSKGradientMatrix
virtual void computeReducedSKGradientMatrix(FloatMatrix &gradientMatrix, int i, GaussPoint *gp, const FloatArray &fullStressVector, const FloatArray &strainSpaceHardeningVariables)=0
Computes second derivative of loading function with respect to stress and hardening vars...
oofem::MPlasticMaterial2::linearElasticMaterial
LinearElasticMaterial * linearElasticMaterial
Reference to bulk (undamaged) material.
Definition: mplasticmaterial2.h:179
oofem::MPlasticMaterial2::iterativeUpdateOfActiveConds
bool iterativeUpdateOfActiveConds
Flag indicating whether iterative update of a set of active yield conditions takes place...
Definition: mplasticmaterial2.h:188
oofem::principal_strain
For computing principal strains from engineering strains.
Definition: stressstrainprincmode.h:45
oofem::FloatMatrix::solveForRhs
bool solveForRhs(const FloatArray &b, FloatArray &answer, bool transpose=false)
Solves the system of linear equations .
Definition: floatmatrix.C:1112
oofem::StructuralMaterial::giveFiberStiffMtrx
virtual void giveFiberStiffMtrx(FloatMatrix &answer, MatResponseMode mmode, GaussPoint *gp, TimeStep *tStep)
Method for computing 1d fiber stiffness matrix of receiver.
Definition: structuralmaterial.C:1192
oofem::MPlasticMaterial2::computeReducedElasticModuli
virtual void computeReducedElasticModuli(FloatMatrix &answer, GaussPoint *gp, TimeStep *tStep)
Definition: mplasticmaterial2.C:1657
oofem::MPlasticMaterial2Status::initTempStatus
virtual void initTempStatus()
Initializes the temporary internal variables, describing the current state according to previously re...
Definition: mplasticmaterial2.C:2122
oofem::DataStream
The purpose of DataStream abstract class is to allow to store/restore context to different streams...
Definition: datastream.h:54
oofem::FloatArray::storeYourself
contextIOResultType storeYourself(DataStream &stream) const
Definition: floatarray.C:872
oofem::StructuralMaterial::giveIPValue
virtual int giveIPValue(FloatArray &answer, GaussPoint *gp, InternalStateType type, TimeStep *tStep)
Returns the integration point corresponding value in Reduced form.
Definition: structuralmaterial.C:2115
oofem::MPlasticMaterial2::nsurf
int nsurf
Number of yield surfaces.
Definition: mplasticmaterial2.h:181
oofem::IntArray::zero
void zero()
Sets all component to zero.
Definition: intarray.C:52
oofem::FloatArray::at
double & at(int i)
Coefficient access function.
Definition: floatarray.h:131
oofem::MPlasticMaterial2::computeReducedStressGradientVector
void computeReducedStressGradientVector(FloatArray &answer, functType ftype, int isurf, GaussPoint *gp, const FloatArray &stressVector, const FloatArray &strainSpaceHardeningVariables)
Definition: mplasticmaterial2.C:1216
oofem::StructuralMaterial::givePlaneStrainStiffMtrx
virtual void givePlaneStrainStiffMtrx(FloatMatrix &answer, MatResponseMode mmode, GaussPoint *gp, TimeStep *tStep)
Method for computing plane strain stiffness matrix of receiver.
Definition: structuralmaterial.C:1075
oofem::sgn
double sgn(double i)
Returns the signum of given value (if value is < 0 returns -1, otherwise returns 1) ...
Definition: mathfem.h:91
oofem::MPlasticMaterial2::computeTrialStressIncrement
virtual void computeTrialStressIncrement(FloatArray &answer, GaussPoint *gp, const FloatArray &strainIncrement, TimeStep *tStep)
Computes full-space trial stress increment (elastic).
Definition: mplasticmaterial2.C:1260
oofem::StructuralMaterialStatus
This class implements a structural material status information.
Definition: structuralms.h:65
oofem::MPlasticMaterial2::computeKGradientVector
virtual void computeKGradientVector(FloatArray &answer, functType ftype, int isurf, GaussPoint *gp, FloatArray &fullStressVector, const FloatArray &strainSpaceHardeningVariables)=0
Computes the derivative of yield/loading function with respect to vector.
oofem::MPlasticMaterial2Status::giveTempGamma
const FloatArray & giveTempGamma()
Definition: mplasticmaterial2.h:132
oofem::StructuralMaterial::givePlateLayerStiffMtrx
virtual void givePlateLayerStiffMtrx(FloatMatrix &answer, MatResponseMode mmode, GaussPoint *gp, TimeStep *tStep)
Method for computing 2d plate layer stiffness matrix of receiver.
Definition: structuralmaterial.C:1154
PLASTIC_MATERIAL_MAX_ITERATIONS
#define PLASTIC_MATERIAL_MAX_ITERATIONS
Definition: mplasticmaterial2.C:48
oofem::GaussPoint::giveElement
Element * giveElement()
Returns corresponding element to receiver.
Definition: gausspoint.h:188
oofem::CIO_IOERR
General IO error.
Definition: contextioresulttype.h:43
oofem::FloatMatrix::negated
void negated()
Changes sign of receiver values.
Definition: floatmatrix.C:1605
oofem::MPlasticMaterial2::plType
enum oofem::MPlasticMaterial2::plastType plType
oofem::MPlasticMaterial2Status::saveContext
virtual contextIOResultType saveContext(DataStream &stream, ContextMode mode, void *obj=NULL)
Stores receiver state to output stream.
Definition: mplasticmaterial2.C:2166
oofem::MPlasticMaterial2::cuttingPlaneReturn
void cuttingPlaneReturn(FloatArray &answer, IntArray &activeConditionMap, FloatArray &gamma, GaussPoint *gp, const FloatArray &totalStrain, FloatArray &plasticStrainR, FloatArray &strainSpaceHardeningVariables, TimeStep *tStep)
Definition: mplasticmaterial2.C:832
oofem::StructuralMaterial::giveStressDependentPartOfStrainVector
void giveStressDependentPartOfStrainVector(FloatArray &answer, GaussPoint *gp, const FloatArray &reducedStrainVector, TimeStep *tStep, ValueModeType mode)
Method for subtracting from reduced space strain vector its stress-independent parts (caused by tempe...
Definition: structuralmaterial.C:768
oofem::MaterialMode
MaterialMode
Type representing material mode of integration point.
Definition: materialmode.h:89
oofem::MPlasticMaterial2::givePlaneStressStiffMtrx
virtual void givePlaneStressStiffMtrx(FloatMatrix &answer, MatResponseMode, GaussPoint *gp, TimeStep *tStep)
Method for computing plane stress stiffness matrix of receiver.
Definition: mplasticmaterial2.C:1714
oofem::StructuralMaterialStatus::printOutputAt
virtual void printOutputAt(FILE *file, TimeStep *tStep)
Print receiver&#39;s output to given stream.
Definition: structuralms.C:73
oofem::MPlasticMaterial2Status
This class implements associated Material Status to MPlasticMaterial.
Definition: mplasticmaterial2.h:67
oofem::IntArray
Class implementing an array of integers.
Definition: intarray.h:61
oofem::IntArray::at
int & at(int i)
Coefficient access function.
Definition: intarray.h:103
oofem::DataStream::read
virtual int read(int *data, int count)=0
Reads count integer values into array pointed by data.
oofem::MatResponseMode
MatResponseMode
Describes the character of characteristic material matrix.
Definition: matresponsemode.h:64
oofem::MPlasticMaterial2::functType
functType
Type that allows to distinguish between yield function and loading function.
Definition: mplasticmaterial2.h:185
THROW_CIOERR
#define THROW_CIOERR(e)
Definition: contextioerr.h:61
oofem::MPlasticMaterial2::giveIPValue
virtual int giveIPValue(FloatArray &answer, GaussPoint *gp, InternalStateType type, TimeStep *tStep)
Returns the integration point corresponding value in Reduced form.
Definition: mplasticmaterial2.C:1862
oofem::MPlasticMaterial2Status::restoreContext
virtual contextIOResultType restoreContext(DataStream &stream, ContextMode mode, void *obj=NULL)
Restores the receiver state previously written in stream.
Definition: mplasticmaterial2.C:2207
RES_TOL
#define RES_TOL
Definition: mplasticmaterial2.C:47
oofem::StructuralMaterial::giveFullSymVectorForm
static void giveFullSymVectorForm(FloatArray &answer, const FloatArray &vec, MaterialMode matMode)
Converts the reduced unsymmetric Voigt vector (2nd order tensor) to full form.
Definition: structuralmaterial.C:2449
oofem::MPlasticMaterial2::populationSet
std::set< long > populationSet
Set for keeping record of generated populations of active yield conditions during return...
Definition: mplasticmaterial2.h:190
oofem::DataStream::write
virtual int write(const int *data, int count)=0
Writes count integer values from array pointed by data.
oofem::MPlasticMaterial2Status::letTempPlasticStrainVectorBe
void letTempPlasticStrainVectorBe(FloatArray v)
Definition: mplasticmaterial2.h:118
oofem::StructuralMaterial::giveStiffnessMatrix
virtual void giveStiffnessMatrix(FloatMatrix &answer, MatResponseMode mode, GaussPoint *gp, TimeStep *tStep)
Computes the stiffness matrix for giveRealStressVector of receiver in given integration point...
Definition: structuralmaterial.C:638
oofem::MPlasticMaterial2::giveRealStressVector
virtual void giveRealStressVector(FloatArray &answer, GaussPoint *, const FloatArray &, TimeStep *)
Computes the real stress vector for given total strain and integration point.
Definition: mplasticmaterial2.C:95
oofem::StructuralMaterial::giveReducedSymVectorForm
static void giveReducedSymVectorForm(FloatArray &answer, const FloatArray &vec, MaterialMode matMode)
Converts the full unsymmetric Voigt vector (2nd order tensor) to reduced form.
Definition: structuralmaterial.C:2497
oofem::MPlasticMaterial2Status::strainSpaceHardeningVarsVector
FloatArray strainSpaceHardeningVarsVector
Strain space hardening variables, e.g.
Definition: mplasticmaterial2.h:81
OOFEM_ERROR
#define OOFEM_ERROR(...)
Definition: error.h:61
oofem::MPlasticMaterial2::addNewPopulation
void addNewPopulation(IntArray &mask)
Definition: mplasticmaterial2.C:1920
oofem::GaussPoint::giveNumber
int giveNumber()
Returns number of receiver.
Definition: gausspoint.h:184
oofem::StructuralMaterial::give3dMaterialStiffnessMatrix
virtual void give3dMaterialStiffnessMatrix(FloatMatrix &answer, MatResponseMode mode, GaussPoint *gp, TimeStep *tStep)
Computes full 3d material stiffness matrix at given integration point, time, respecting load history ...
Definition: structuralmaterial.h:344
oofem::FloatMatrix::times
void times(double f)
Multiplies receiver by factor f.
Definition: floatmatrix.C:1594
oofem::FloatArray::restoreYourself
contextIOResultType restoreYourself(DataStream &stream)
Definition: floatarray.C:895
oofem::MPlasticMaterial2Status::state_flag
int state_flag
Yield function status indicator.
Definition: mplasticmaterial2.h:85
oofem::FloatArray::beProductOf
void beProductOf(const FloatMatrix &aMatrix, const FloatArray &anArray)
Receiver becomes the result of the product of aMatrix and anArray.
Definition: floatarray.C:676
oofem::MPlasticMaterial2::rmType
enum oofem::MPlasticMaterial2::ReturnMappingAlgoType rmType
oofem::MPlasticMaterial2::computeStressGradientVector
virtual void computeStressGradientVector(FloatArray &answer, functType ftype, int isurf, GaussPoint *gp, const FloatArray &stressVector, const FloatArray &strainSpaceHardeningVariables)=0
Computes the stress gradient of yield/loading function (df/d_sigma).
oofem::MPlasticMaterial2::giveLinearElasticMaterial
LinearElasticMaterial * giveLinearElasticMaterial()
Returns reference to undamaged (bulk) material.
Definition: mplasticmaterial2.h:202
oofem::FloatArray::beTProductOf
void beTProductOf(const FloatMatrix &aMatrix, const FloatArray &anArray)
Receiver becomes the result of the product of aMatrix^T and anArray.
Definition: floatarray.C:708
oofem::FloatMatrix::at
double at(int i, int j) const
Coefficient access function.
Definition: floatmatrix.h:176
oofem::IntArray::resize
void resize(int n)
Checks size of receiver towards requested bounds.
Definition: intarray.C:124
oofem::MPlasticMaterial2::getPopulationSignature
long getPopulationSignature(IntArray &mask)
Definition: mplasticmaterial2.C:1886
oofem::StructuralMaterialStatus::restoreContext
virtual contextIOResultType restoreContext(DataStream &stream, ContextMode mode, void *obj=NULL)
Restores the receiver state previously written in stream.
Definition: structuralms.C:157
oofem::IntArray::restoreYourself
contextIOResultType restoreYourself(DataStream &stream)
Restores array from image on stream.
Definition: intarray.C:305
YIELD_TOL_SECONDARY
#define YIELD_TOL_SECONDARY
Definition: mplasticmaterial2.C:46
oofem::MPlasticMaterial2::giveFiberStiffMtrx
virtual void giveFiberStiffMtrx(FloatMatrix &answer, MatResponseMode, GaussPoint *gp, TimeStep *tStep)
Method for computing 1d fiber stiffness matrix of receiver.
Definition: mplasticmaterial2.C:1837
oofem::MaterialStatus
Abstract base class representing a material status information.
Definition: matstatus.h:84
oofem::FloatArray
Class representing vector of real numbers.
Definition: floatarray.h:82
oofem::MPlasticMaterial2::computeReducedHardeningVarsSigmaGradient
virtual void computeReducedHardeningVarsSigmaGradient(FloatMatrix &answer, GaussPoint *gp, const IntArray &activeConditionMap, const FloatArray &fullStressVector, const FloatArray &strainSpaceHardeningVars, const FloatArray &gamma)=0
Computes derivative of vector with respect to stress.
oofem::MPlasticMaterial2Status::giveTempStrainSpaceHardeningVarsVector
const FloatArray & giveTempStrainSpaceHardeningVarsVector() const
Returns the actual (temp) hardening variable vector.
Definition: mplasticmaterial2.h:115
oofem::MPlasticMaterial2::computeAlgorithmicModuli
void computeAlgorithmicModuli(FloatMatrix &answer, GaussPoint *gp, const FloatMatrix &elasticModuliInverse, const FloatArray &gamma, const IntArray &activeConditionMap, const FloatArray &fullStressVector, const FloatArray &strainSpaceHardeningVariables)
Definition: mplasticmaterial2.C:1277
oofem::FloatMatrix
Implementation of matrix containing floating point numbers.
Definition: floatmatrix.h:94
oofem::MPlasticMaterial2Status::activeConditionMap
IntArray activeConditionMap
Active set of yield functions (needed for algorithmic stiffness).
Definition: mplasticmaterial2.h:94
oofem::MPlasticMaterial2::computeDamage
virtual double computeDamage(GaussPoint *gp, const FloatArray &strainSpaceHardeningVariables, TimeStep *tStep)
Definition: mplasticmaterial2.C:178
oofem::StructuralMaterialStatus::letTempStressVectorBe
void letTempStressVectorBe(const FloatArray &v)
Assigns tempStressVector to given vector v.
Definition: structuralms.h:135
oofem::MPlasticMaterial2Status::giveStrainSpaceHardeningVars
const FloatArray & giveStrainSpaceHardeningVars() const
Returns the equilibrated hardening variable vector.
Definition: mplasticmaterial2.h:113
oofem::StructuralMaterialStatus::giveStressVector
const FloatArray & giveStressVector() const
Returns the const pointer to receiver&#39;s stress vector.
Definition: structuralms.h:107
oofem::MPlasticMaterial2::give1dStressStiffMtrx
virtual void give1dStressStiffMtrx(FloatMatrix &answer, MatResponseMode, GaussPoint *gp, TimeStep *tStep)
Method for computing 1d stiffness matrix of receiver.
Definition: mplasticmaterial2.C:1764
oofem::FloatArray::computeNorm
double computeNorm() const
Computes the norm (or length) of the vector.
Definition: floatarray.C:840
oofem::FloatMatrix::resize
void resize(int rows, int cols)
Checks size of receiver towards requested bounds.
Definition: floatmatrix.C:1358
oofem::FloatArray::printYourself
virtual void printYourself() const
Print receiver on stdout.
Definition: floatarray.C:747
oofem::MPlasticMaterial2::CreateStatus
virtual MaterialStatus * CreateStatus(GaussPoint *gp) const
Creates new copy of associated status and inserts it into given integration point.
Definition: mplasticmaterial2.C:88
oofem::FloatMatrix::add
void add(const FloatMatrix &a)
Adds matrix to the receiver.
Definition: floatmatrix.C:1023
oofem::FloatArray::zero
void zero()
Zeroes all coefficients of receiver.
Definition: floatarray.C:658
oofem::MPlasticMaterial2Status::setTempActiveConditionMap
void setTempActiveConditionMap(IntArray v)
Definition: mplasticmaterial2.h:131
oofem::MPlasticMaterial2::closestPointReturn
void closestPointReturn(FloatArray &answer, IntArray &activeConditionMap, FloatArray &gamma, GaussPoint *gp, const FloatArray &totalStrain, FloatArray &plasticStrainR, FloatArray &strainSpaceHardeningVariables, TimeStep *tStep)
Definition: mplasticmaterial2.C:184
oofem::MPlasticMaterial2Status::letTempStrainSpaceHardeningVarsVectorBe
void letTempStrainSpaceHardeningVarsVectorBe(FloatArray v)
Definition: mplasticmaterial2.h:119
oofem::MPlasticMaterial2::computeReducedHardeningVarsLamGradient
virtual void computeReducedHardeningVarsLamGradient(FloatMatrix &answer, GaussPoint *gp, int actSurf, const IntArray &activeConditionMap, const FloatArray &fullStressVector, const FloatArray &strainSpaceHardeningVars, const FloatArray &gamma)=0
computes derivative of vector with respect to lambda vector
oofem::MPlasticMaterial2::MPlasticMaterial2
MPlasticMaterial2(int n, Domain *d)
Definition: mplasticmaterial2.C:50
oofem::FloatArray::times
void times(double s)
Multiplies receiver with scalar.
Definition: floatarray.C:818
oofem::ContextMode
long ContextMode
Context mode (mask), defining the type of information written/read to/from context.
Definition: contextmode.h:43
oofem::MPlasticMaterial2::givePlaneStrainStiffMtrx
virtual void givePlaneStrainStiffMtrx(FloatMatrix &answer, MatResponseMode, GaussPoint *gp, TimeStep *tStep)
Method for computing plane strain stiffness matrix of receiver.
Definition: mplasticmaterial2.C:1740
oofem::MPlasticMaterial2::computeStrainHardeningVarsIncrement
virtual void computeStrainHardeningVarsIncrement(FloatArray &answer, GaussPoint *gp, const FloatArray &stress, const FloatArray &dlambda, const FloatArray &dplasticStrain, const IntArray &activeConditionMap)=0
Computes the increment of strain-space hardening variables.
oofem::StructuralMaterial
Abstract base class for all "structural" constitutive models.
Definition: structuralmaterial.h:113
oofem::MPlasticMaterial2::computeResidualVector
void computeResidualVector(FloatArray &answer, GaussPoint *gp, const FloatArray &gamma, const IntArray &activeConditionMap, const FloatArray &plasticStrainVectorR, const FloatArray &strainSpaceHardeningVariables, std::vector< FloatArray > &gradVec)
Definition: mplasticmaterial2.C:1229
oofem::FloatMatrix::zero
void zero()
Zeroes all coefficient of receiver.
Definition: floatmatrix.C:1326
oofem::StructuralMaterialStatus::initTempStatus
virtual void initTempStatus()
Initializes the temporary internal variables, describing the current state according to previously re...
Definition: structuralms.C:108
oofem::MPlasticMaterial2Status::updateYourself
virtual void updateYourself(TimeStep *tStep)
Update equilibrium history variables according to temp-variables.
Definition: mplasticmaterial2.C:2145
oofem::FEMComponent::giveDomain
Domain * giveDomain() const
Definition: femcmpnn.h:100
oofem::StructuralMaterial::give2dBeamLayerStiffMtrx
virtual void give2dBeamLayerStiffMtrx(FloatMatrix &answer, MatResponseMode mmode, GaussPoint *gp, TimeStep *tStep)
Method for computing 2d beam layer stiffness matrix of receiver.
Definition: structuralmaterial.C:1130
oofem::MPlasticMaterial2::computeYieldValueAt
virtual double computeYieldValueAt(GaussPoint *gp, int isurf, const FloatArray &stressVector, const FloatArray &strainSpaceHardeningVariables)=0
Computes the value of yield function.
oofem::MPlasticMaterial2::give2dBeamLayerStiffMtrx
virtual void give2dBeamLayerStiffMtrx(FloatMatrix &answer, MatResponseMode, GaussPoint *gp, TimeStep *tStep)
Method for computing 2d beam layer stiffness matrix of receiver.
Definition: mplasticmaterial2.C:1787
oofem::MPlasticMaterial2Status::giveTempActiveConditionMap
const IntArray & giveTempActiveConditionMap()
Definition: mplasticmaterial2.h:130
oofem::MPlasticMaterial2::giveElastoPlasticStiffnessMatrix
virtual void giveElastoPlasticStiffnessMatrix(FloatMatrix &answer, MatResponseMode mode, GaussPoint *gp, TimeStep *tStep)
Definition: mplasticmaterial2.C:1495
oofem::MPlasticMaterial2::giveSizeOfReducedHardeningVarsVector
virtual int giveSizeOfReducedHardeningVarsVector(GaussPoint *) const
Definition: mplasticmaterial2.h:236
oofem::FloatMatrix::beProductOf
void beProductOf(const FloatMatrix &a, const FloatMatrix &b)
Assigns to the receiver product of .
Definition: floatmatrix.C:337
oofem::IntArray::giveSize
int giveSize() const
Definition: intarray.h:203
oofem::StructuralMaterialStatus::giveTempStressVector
const FloatArray & giveTempStressVector() const
Returns the const pointer to receiver&#39;s temporary stress vector.
Definition: structuralms.h:117
oofem::MPlasticMaterial2::giveMaxNumberOfActiveYieldConds
virtual int giveMaxNumberOfActiveYieldConds(GaussPoint *gp)=0
oofem::FloatArray::giveSize
int giveSize() const
Returns the size of receiver.
Definition: floatarray.h:218
oofem::StructuralMaterialStatus::saveContext
virtual contextIOResultType saveContext(DataStream &stream, ContextMode mode, void *obj=NULL)
Stores receiver state to output stream.
Definition: structuralms.C:133
oofem::MPlasticMaterial2::giveConsistentStiffnessMatrix
virtual void giveConsistentStiffnessMatrix(FloatMatrix &answer, MatResponseMode, GaussPoint *gp, TimeStep *tStep)
Definition: mplasticmaterial2.C:1339
oofem
the oofem namespace is to define a context or scope in which all oofem names are defined.
oofem::MPlasticMaterial2Status::plasticStrainVector
FloatArray plasticStrainVector
Plastic strain vector.
Definition: mplasticmaterial2.h:74
oofem::Material::initTempStatus
virtual void initTempStatus(GaussPoint *gp)
Initializes temporary variables stored in integration point status at the beginning of new time step...
Definition: material.C:267
oofem::FEMComponent::giveNumber
int giveNumber() const
Definition: femcmpnn.h:107
oofem::FloatMatrix::beInverseOf
void beInverseOf(const FloatMatrix &src)
Modifies receiver to become inverse of given parameter.
Definition: floatmatrix.C:835
oofem::StructuralMaterial::givePlaneStressStiffMtrx
virtual void givePlaneStressStiffMtrx(FloatMatrix &answer, MatResponseMode mmode, GaussPoint *gp, TimeStep *tStep)
Method for computing plane stress stiffness matrix of receiver.
Definition: structuralmaterial.C:1042
oofem::FloatMatrix::giveNumberOfRows
int giveNumberOfRows() const
Returns number of rows of receiver.
Definition: floatmatrix.h:156
oofem::MPlasticMaterial2::hasMaterialModeCapability
virtual int hasMaterialModeCapability(MaterialMode mode)
Tests if material supports material mode.
Definition: mplasticmaterial2.C:73
oofem::MPlasticMaterial2::getNewPopulation
int getNewPopulation(IntArray &result, IntArray &candidateMask, int degree, int size)
Definition: mplasticmaterial2.C:1928
oofem::GaussPoint
Class representing integration point in finite element program.
Definition: gausspoint.h:93
OOFEM_WARNING
#define OOFEM_WARNING(...)
Definition: error.h:62
oofem::MPlasticMaterial2Status::givePlasticStrainVector
const FloatArray & givePlasticStrainVector() const
Returns the equilibrated strain vector.
Definition: mplasticmaterial2.h:109
oofem::TimeStep
Class representing solution step.
Definition: timestep.h:80
oofem::FloatArray::add
void add(const FloatArray &src)
Adds array src to receiver.
Definition: floatarray.C:156
oofem::MPlasticMaterial2::hasHardening
virtual int hasHardening()=0
Indicates, whether receiver model has hardening/softening behavior or behaves according to perfect pl...
oofem::MPlasticMaterial2::give3dMaterialStiffnessMatrix
virtual void give3dMaterialStiffnessMatrix(FloatMatrix &answer, MatResponseMode, GaussPoint *gp, TimeStep *tStep)
Computes full 3d material stiffness matrix at given integration point, time, respecting load history ...
Definition: mplasticmaterial2.C:1668
oofem::MPlasticMaterial2Status::printOutputAt
virtual void printOutputAt(FILE *file, TimeStep *tStep)
Print receiver&#39;s output to given stream.
Definition: mplasticmaterial2.C:2089
oofem::FloatArray::resize
void resize(int s)
Resizes receiver towards requested size.
Definition: floatarray.C:631
This page is part of the OOFEM documentation. Copyright (c) 2011 Borek Patzak
Project e-mail: info@oofem.org
Generated at Tue Jan 2 2018 20:07:30 for OOFEM by doxygen 1.8.11 written by Dimitri van Heesch, © 1997-2011