OOFEM  2.4
OOFEM.org - Object Oriented Finite Element Solver
mplasticmaterial.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 "mplasticmaterial.h"
36 #include "gausspoint.h"
37 #include "floatmatrix.h"
38 #include "floatarray.h"
39 #include "intarray.h"
40 #include "datastream.h"
41 #include "contextioerr.h"
42 
43 namespace oofem {
44 #define YIELD_TOL 1.e-4
45 #define RES_TOL 1.e-4
46 #define PLASTIC_MATERIAL_MAX_ITERATIONS 90
47 
49  //
50  // constructor
51  //
52 {
53  linearElasticMaterial = NULL;
56 }
57 
58 
60 //
61 // destructor
62 //
63 {
64  delete linearElasticMaterial;
65 }
66 
67 
68 int
70 //
71 // returns whether receiver supports given mode
72 //
73 {
74  return mode == _3dMat ||
75  mode == _1dMat ||
76  //<RESTRICTED_SECTION>
77  mode == _PlaneStress ||
78  //</RESTRICTED_SECTION>
79  mode == _PlaneStrain;
80 }
81 
82 
85 {
86  return new MPlasticMaterialStatus(1, this->giveDomain(), gp, this->giveSizeOfReducedHardeningVarsVector(gp));
87 }
88 
89 
90 void
92  GaussPoint *gp,
93  const FloatArray &totalStrain,
94  TimeStep *tStep)
95 //
96 // returns real stress vector in 3d stress space of receiver according to
97 // previous level of stress and current
98 // strain increment, the only way, how to correctly update gp records
99 //
100 // completely formulated in Reduced stress-strain space
101 {
102  FloatArray strainSpaceHardeningVariables;
103  FloatArray fullStressVector;
104  FloatArray strainVectorR, plasticStrainVectorR;
105  FloatArray helpVec;
106  IntArray activeConditionMap(this->nsurf);
107  FloatArray gamma;
108 
109  MPlasticMaterialStatus *status = static_cast< MPlasticMaterialStatus * >( this->giveStatus(gp) );
110 
111  this->initTempStatus(gp);
112 
113  // subtract stress independent part
114  // note: eigenStrains (temperature) is not contained in mechanical strain stored in gp
115  // therefore it is necessary to subtract always the total eigen strain value
116  this->giveStressDependentPartOfStrainVector(strainVectorR, gp, totalStrain,
117  tStep, VM_Total);
118 
119  /*
120  * stress return algorithm
121  */
122 
123  //if (0) this->cuttingPlaneReturn (fullStressVector, activeConditionMap, gamma, form, gp, totalStrain, tStep);
124  //else this->closestPointReturn (fullStressVector, activeConditionMap, gamma, form, gp, totalStrain, tStep);
125  if ( rmType == mpm_ClosestPoint ) {
126  this->closestPointReturn(fullStressVector, activeConditionMap, gamma, gp, strainVectorR,
127  plasticStrainVectorR, strainSpaceHardeningVariables, tStep);
128  } else {
129  this->cuttingPlaneReturn(fullStressVector, activeConditionMap, gamma, gp, strainVectorR,
130  plasticStrainVectorR, strainSpaceHardeningVariables, tStep);
131  }
132 
133 
134  status->letTempStrainVectorBe(totalStrain);
135  StructuralMaterial :: giveReducedSymVectorForm( helpVec, fullStressVector, gp->giveMaterialMode() );
136  status->letTempStressVectorBe(helpVec);
137 
138  status->letTempPlasticStrainVectorBe(plasticStrainVectorR);
139  status->letTempStrainSpaceHardeningVarsVectorBe(strainSpaceHardeningVariables);
140 
141  status->setTempGamma(gamma);
142  status->setTempActiveConditionMap(activeConditionMap);
143 
144 
145  // update state flag
146  int newState, state = status->giveStateFlag();
147  bool yieldFlag = false;
148  for ( int i = 1; i <= nsurf; i++ ) {
149  if ( gamma.at(i) > 0. ) {
150  yieldFlag = true;
151  }
152  }
153 
154  if ( yieldFlag ) {
155  newState = MPlasticMaterialStatus :: PM_Yielding; // test if plastic loading occur
156  }
157  // no plastic loading - check for unloading
158  else if ( ( state == MPlasticMaterialStatus :: PM_Yielding ) || ( state == MPlasticMaterialStatus :: PM_Unloading ) ) {
160  } else {
162  }
163 
164  status->letTempStateFlagBe(newState);
165 
166  answer = status->giveTempStressVector();
167 }
168 
169 
170 void
172  IntArray &activeConditionMap,
173  FloatArray &gamma,
174  GaussPoint *gp,
175  const FloatArray &totalStrain,
176  FloatArray &plasticStrainVectorR,
177  FloatArray &strainSpaceHardeningVariables,
178  TimeStep *tStep)
179 {
180  FloatArray fullStressVector;
181  FloatArray elasticStrainVectorR;
182  FloatArray fullStressSpaceHardeningVars, residualVectorR;
183 
184  FloatArray helpVector, helpVector2;
185  FloatArray dgamma, tempGamma;
186  FloatMatrix elasticModuli, hardeningModuli, consistentModuli;
187  FloatMatrix elasticModuliInverse, hardeningModuliInverse;
188  FloatMatrix helpMtrx;
189  FloatMatrix gmat;
190  std :: vector< FloatArray > yieldGradVec(this->nsurf), loadGradVec(this->nsurf), * yieldGradVecPtr, * loadGradVecPtr;
191  FloatArray rhs;
192  double yieldValue;
193  int nIterations = 0;
194  int i, j, strSize, totSize;
195  int elastic, restart, actSurf, indx;
196  bool yieldConsistency;
197 
198 
199  if ( this->plType == associatedPT ) {
200  yieldGradVecPtr = loadGradVecPtr = & yieldGradVec;
201  } else {
202  yieldGradVecPtr = & yieldGradVec;
203  loadGradVecPtr = & loadGradVec;
204  }
205 
206  MPlasticMaterialStatus *status = static_cast< MPlasticMaterialStatus * >( this->giveStatus(gp) );
207 
208  plasticStrainVectorR = status->givePlasticStrainVector();
209  strainSpaceHardeningVariables = status->giveStrainSpaceHardeningVars();
210 
211  dgamma.resize(nsurf);
212  dgamma.zero();
213  gamma.resize(nsurf);
214  gamma.zero();
215 
216  strSize = totalStrain.giveSize(); // size of reducedStrain Vector
217  totSize = strSize + strainSpaceHardeningVariables.giveSize();
218 
219  // compute elastic moduli and its inverse
220  this->computeReducedElasticModuli(elasticModuli, gp, tStep);
221  elasticModuliInverse.beInverseOf(elasticModuli);
222  // delete elasticModuli;
223 
224  //
225  // Elastic Predictor
226  //
227 
228  elasticStrainVectorR = totalStrain;
229  elasticStrainVectorR.subtract(plasticStrainVectorR);
230  // stress vector in full form due to computational convinience
231  //if (fullStressVector) delete fullStressVector;
232  this->computeTrialStressIncrement(fullStressVector, gp, elasticStrainVectorR, tStep);
233  this->computeStressSpaceHardeningVars(fullStressSpaceHardeningVars, gp, strainSpaceHardeningVariables);
234 
235  //
236  // check for plastic process
237  //
238  elastic = 1;
239  actSurf = 0;
240  activeConditionMap.zero();
241  for ( i = 1; i <= this->nsurf; i++ ) {
242  yieldValue = this->computeYieldValueAt(gp, i, fullStressVector, fullStressSpaceHardeningVars);
243  if ( yieldValue > YIELD_TOL ) {
244  elastic = 0;
245  actSurf += 1;
246  activeConditionMap.at(i) = actSurf;
247  }
248  }
249 
250  if ( !elastic ) {
251  /*
252  * plastic multisurface closest-point corrector
253  */
254 
255 
256  do {
257  do { // restart loop
258  elasticStrainVectorR = totalStrain;
259  elasticStrainVectorR.subtract(plasticStrainVectorR);
260  // stress vector in full form due to computational convinience
261  //if (fullStressVector) delete fullStressVector;
262  this->computeTrialStressIncrement(fullStressVector, gp, elasticStrainVectorR, tStep);
263  this->computeStressSpaceHardeningVars(fullStressSpaceHardeningVars, gp, strainSpaceHardeningVariables);
264 
265  // compute gradients
266  for ( i = 1; i <= nsurf; i++ ) {
267  this->computeGradientVector(yieldGradVec [ i - 1 ], yieldFunction, i, gp, fullStressVector, fullStressSpaceHardeningVars);
268  }
269 
270  if ( this->plType == nonassociatedPT ) {
271  for ( i = 1; i <= nsurf; i++ ) {
272  this->computeGradientVector(loadGradVec [ i - 1 ], loadFunction, i, gp, fullStressVector, fullStressSpaceHardeningVars);
273  }
274  }
275 
276 
277  //this->computeGradientVector (gradientVectorR, gp, &fullStressVector, fullStressSpaceHardeningVars);
278  this->computeResidualVector(residualVectorR, gp, gamma, activeConditionMap, plasticStrainVectorR,
279  strainSpaceHardeningVariables, * loadGradVecPtr);
280 
281  // check convergence
282  yieldConsistency = true;
283  for ( i = 1; i <= nsurf; i++ ) {
284  if ( activeConditionMap.at(i) ) {
285  yieldValue = this->computeYieldValueAt(gp, i, fullStressVector, fullStressSpaceHardeningVars);
286  if ( ( yieldValue >= YIELD_TOL ) ) {
287  yieldConsistency = false;
288  }
289  }
290  }
291 
292  if ( yieldConsistency && ( residualVectorR.computeNorm() < RES_TOL ) ) {
293  answer = fullStressVector;
294  printf(" (%d iterations)", nIterations);
295  return;
296  }
297 
298  // compute consistent tangent moduli
299  this->computeHardeningReducedModuli(hardeningModuli, gp, strainSpaceHardeningVariables, tStep);
300  if ( hardeningModuli.giveNumberOfRows() ) {
301  hardeningModuliInverse.beInverseOf(hardeningModuli);
302  // delete hardeningModuli;
303  } else {
304  hardeningModuliInverse.clear();
305  }
306 
307  this->computeAlgorithmicModuli(consistentModuli, gp, elasticModuliInverse,
308  hardeningModuliInverse, gamma, activeConditionMap,
309  fullStressVector, fullStressSpaceHardeningVars);
310 
311  rhs.resize(actSurf);
312  gmat.resize(actSurf, actSurf);
313  for ( i = 1; i <= nsurf; i++ ) {
314  if ( activeConditionMap.at(i) ) {
315  yieldValue = this->computeYieldValueAt(gp, i, fullStressVector, fullStressSpaceHardeningVars);
316  helpVector.beTProductOf(consistentModuli, ( * yieldGradVecPtr ) [ i - 1 ]);
317 
318  for ( j = 1; j <= this->nsurf; j++ ) {
319  if ( activeConditionMap.at(j) ) {
320  gmat.at(i, j) = ( * loadGradVecPtr ) [ j - 1 ].dotProduct(helpVector);
321  }
322  }
323 
324  rhs.at(i) = yieldValue - residualVectorR.dotProduct(helpVector);
325  }
326  }
327 
328  // obtain increment to consistency parameter
329  gmat.solveForRhs(rhs, helpVector);
330 
331  // assign gammas
332  for ( i = 1; i <= this->nsurf; i++ ) {
333  if ( ( indx = activeConditionMap.at(i) ) ) {
334  dgamma.at(i) = helpVector.at(indx);
335  } else {
336  dgamma.at(i) = 0.0;
337  }
338  }
339 
340  tempGamma = gamma;
341  tempGamma.add(dgamma);
342 
343  //
344  // check active constraint
345  //
346 
347  restart = 0;
348  for ( i = 1; i <= this->nsurf; i++ ) {
349  if ( tempGamma.at(i) < 0.0 ) {
350  restart = 1;
351  }
352  }
353 
354  if ( restart ) {
355  // build up new activeConditionMap and restart the iterartion from begining
356  actSurf = 0;
357  activeConditionMap.zero();
358  for ( i = 1; i <= this->nsurf; i++ ) {
359  if ( tempGamma.at(i) > 0.0 ) {
360  actSurf += 1;
361  activeConditionMap.at(i) = actSurf;
362  }
363  }
364 
365  continue;
366  } else {
367  break;
368  }
369  } while ( 1 ); // end restart loop
370 
371  // obtain incremental plastic strains and internal variables
372  // residualVectorR vector and gradVec array are changed
373  // but they are computed once again when new iteration starts
374  for ( i = 1; i <= nsurf; i++ ) {
375  if ( activeConditionMap.at(i) ) {
376  ( * loadGradVecPtr ) [ i - 1 ].times( dgamma.at(i) );
377  residualVectorR.add( ( * loadGradVecPtr ) [ i - 1 ] );
378  }
379  }
380 
381  helpVector.beProductOf(consistentModuli, residualVectorR);
382  this->computeDiagModuli(helpMtrx, gp, elasticModuliInverse, hardeningModuliInverse);
383  helpVector2.beProductOf(helpMtrx, helpVector);
384 
385  // Update state variables and consistency parameter
386  for ( i = 1; i <= strSize; i++ ) {
387  plasticStrainVectorR.at(i) += helpVector2.at(i);
388  }
389 
390  for ( i = strSize + 1; i <= totSize; i++ ) {
391  strainSpaceHardeningVariables.at(i - strSize) += helpVector2.at(i);
392  }
393 
394  gamma.add(dgamma);
395 
396  // increment iteration count
397  nIterations++;
398 
399  if ( nIterations > PLASTIC_MATERIAL_MAX_ITERATIONS ) {
400  OOFEM_WARNING("local equlibrium not reached in %d iterations\nElement %d, gp %d, continuing",
402  answer = fullStressVector;
403  break;
404  }
405  } while ( 1 );
406  } else {
407  answer = fullStressVector;
408  }
409 }
410 
411 
412 void
414  IntArray &activeConditionMap,
415  FloatArray &gamma,
416  GaussPoint *gp,
417  const FloatArray &totalStrain,
418  FloatArray &plasticStrainVectorR,
419  FloatArray &strainSpaceHardeningVariables,
420  TimeStep *tStep)
421 {
422  FloatArray elasticStrainVectorR;
423  FloatArray fullStressVector, fullStressSpaceHardeningVars;
424  FloatArray helpVector, trialStressIncrement, rhs;
425  FloatMatrix elasticModuli, hardeningModuli, dmat, gmatInv;
426  std :: vector< FloatArray > yieldGradVec(this->nsurf), loadGradVec(this->nsurf), * yieldGradVecPtr, * loadGradVecPtr;
427  FloatArray dgamma(this->nsurf);
428  double yieldValue;
429  int nIterations = 0;
430  int size, sizeR, i, j, elastic, restart, actSurf, indx, iindx, jindx;
431 
432  MPlasticMaterialStatus *status = static_cast< MPlasticMaterialStatus * >( this->giveStatus(gp) );
433 
434  if ( this->plType == associatedPT ) {
435  yieldGradVecPtr = loadGradVecPtr = & yieldGradVec;
436  } else {
437  yieldGradVecPtr = & yieldGradVec;
438  loadGradVecPtr = & loadGradVec;
439  }
440 
441  // huhu:
442  plasticStrainVectorR = status->givePlasticStrainVector();
443  strainSpaceHardeningVariables = status->giveStrainSpaceHardeningVars();
444 
445  dgamma.resize(nsurf);
446  dgamma.zero();
447  gamma.resize(nsurf);
448  gamma.zero();
449  // compute elastic moduli
450  this->computeReducedElasticModuli(elasticModuli, gp, tStep);
451  // compute hardening moduli and its inverse
452  this->computeHardeningReducedModuli(hardeningModuli, gp, strainSpaceHardeningVariables, tStep);
453 
454  // initialize activeConditionMap
455  elasticStrainVectorR = totalStrain;
456  elasticStrainVectorR.subtract(plasticStrainVectorR);
457  // stress vector in full form due to computational convinience
458  this->computeTrialStressIncrement(fullStressVector, gp, elasticStrainVectorR, tStep);
459  StructuralMaterial :: giveReducedSymVectorForm( trialStressIncrement, fullStressVector, gp->giveMaterialMode() );
460  trialStressIncrement.subtract( status->giveStressVector() );
461  this->computeStressSpaceHardeningVars(fullStressSpaceHardeningVars, gp, strainSpaceHardeningVariables);
462 
463  elastic = 1;
464  actSurf = 0;
465  activeConditionMap.zero();
466  for ( i = 1; i <= this->nsurf; i++ ) {
467  yieldValue = this->computeYieldValueAt(gp, i, fullStressVector, fullStressSpaceHardeningVars);
468  if ( yieldValue > YIELD_TOL ) {
469  elastic = 0;
470  actSurf += 1;
471  activeConditionMap.at(i) = actSurf;
472  }
473  }
474 
475  if ( !elastic ) {
476  //
477  // compute consistent moduli
478  //
479  sizeR = elasticModuli.giveNumberOfRows();
480  size = sizeR + hardeningModuli.giveNumberOfRows();
481 
482  dmat.resize(size, size);
483  dmat.zero();
484 
485  for ( i = 1; i <= sizeR; i++ ) {
486  for ( j = 1; j <= sizeR; j++ ) {
487  dmat.at(i, j) = elasticModuli.at(i, j);
488  }
489  }
490 
491  for ( i = sizeR + 1; i <= size; i++ ) {
492  for ( j = sizeR + 1; j <= size; j++ ) {
493  dmat.at(i, j) = hardeningModuli.at(i - sizeR, j - sizeR);
494  }
495  }
496 
497  do { // local equilibrium loop
498  // update hardening moduli and its inverse
499  this->computeHardeningReducedModuli(hardeningModuli, gp, strainSpaceHardeningVariables, tStep);
500  for ( i = sizeR + 1; i <= size; i++ ) {
501  for ( j = sizeR + 1; j <= size; j++ ) {
502  dmat.at(i, j) = hardeningModuli.at(i - sizeR, j - sizeR);
503  }
504  }
505 
506  do { // loop for restart
507  gmatInv.resize(actSurf, actSurf);
508  gmatInv.zero();
509  rhs.resize(actSurf);
510  rhs.zero();
511 
512  // compute gradients
513  for ( i = 1; i <= nsurf; i++ ) {
514  this->computeGradientVector(yieldGradVec [ i - 1 ], yieldFunction, i, gp, fullStressVector, fullStressSpaceHardeningVars);
515  }
516 
517  if ( this->plType == nonassociatedPT ) {
518  for ( i = 1; i <= nsurf; i++ ) {
519  this->computeGradientVector(loadGradVec [ i - 1 ], loadFunction, i, gp, fullStressVector, fullStressSpaceHardeningVars);
520  }
521  }
522 
523 
524 
525  for ( i = 1; i <= nsurf; i++ ) {
526  if ( ( iindx = activeConditionMap.at(i) ) ) {
527  rhs.at(iindx) = this->computeYieldValueAt(gp, i, fullStressVector, fullStressSpaceHardeningVars);
528  helpVector.beTProductOf(dmat, ( * yieldGradVecPtr ) [ i - 1 ]);
529 
530  for ( j = 1; j <= this->nsurf; j++ ) {
531  if ( ( jindx = activeConditionMap.at(j) ) ) {
532  gmatInv.at(iindx, jindx) = ( * loadGradVecPtr ) [ j - 1 ].dotProduct(helpVector);
533  }
534  }
535  }
536  }
537 
538  /*
539  * // compute plastic multipliers for active conds.
540  * //computee gmatInv
541  * gradMat.resize(actSurf, size);
542  * for (i=1; i<=nsurf; i++) {
543  * if ((iindx = activeConditionMap.at(i))) {
544  * rhs.at (iindx) = this-> computeYieldValueAt (gp, i, fullStressVector, fullStressSpaceHardeningVars);
545  * for (j=1; j<=size; j++) {
546  * gradMat.at(iindx, j) = (*yieldGradVecPtr)[i-1].at(j);
547  * }
548  * }
549  * }
550  * helpMtrx.resize(actSurf, size);
551  * helpMtrx.beProductOf (gradMat, dmat);
552  * gmatInv.beProductTOf (helpMtrx, gradMat);
553  */
554  // solve for plastic multiplier increments
555  gmatInv.solveForRhs(rhs, helpVector);
556  // assign gammas
557  for ( i = 1; i <= this->nsurf; i++ ) {
558  if ( ( indx = activeConditionMap.at(i) ) ) {
559  dgamma.at(i) = helpVector.at(indx);
560  } else {
561  dgamma.at(i) = 0.0;
562  }
563  }
564 
565 
566 
567  restart = 0;
568  // check for active conditions
569  for ( i = 1; i <= this->nsurf; i++ ) {
570  if ( ( gamma.at(i) + dgamma.at(i) ) < 0.0 ) {
571  gamma.at(i) = 0.0;
572  restart = 1;
573  }
574  }
575 
576  if ( restart ) {
577  // build up new activeConditionMap and restart the iterartion from begining
578  actSurf = 0;
579  activeConditionMap.zero();
580  for ( i = 1; i <= this->nsurf; i++ ) {
581  if ( ( gamma.at(i) + dgamma.at(i) ) > 0.0 ) {
582  actSurf += 1;
583  activeConditionMap.at(i) = actSurf;
584  }
585  }
586 
587  continue;
588  } else {
589  break;
590  }
591  } while ( 1 ); // end restart loop
592 
593  // compute plastic strain vector
594  helpVector.resize(sizeR);
595  helpVector.zero();
596  for ( i = 1; i <= this->nsurf; i++ ) {
597  if ( activeConditionMap.at(i) ) {
598  ( * loadGradVecPtr ) [ i - 1 ].times( dgamma.at(i) );
599  for ( j = 1; j <= sizeR; j++ ) {
600  helpVector.at(j) += ( * loadGradVecPtr ) [ i - 1 ].at(j);
601  }
602  }
603  }
604 
605  plasticStrainVectorR.add(helpVector);
606  // update strain space hardening vector
607  helpVector.resize(size - sizeR);
608  helpVector.zero();
609  for ( i = 1; i <= this->nsurf; i++ ) {
610  if ( activeConditionMap.at(i) ) {
611  //gradVec[i-1].times(gamma.at(i));
612  for ( j = sizeR + 1; j <= size; j++ ) {
613  helpVector.at(j - sizeR) += ( * loadGradVecPtr ) [ i - 1 ].at(j);
614  }
615  }
616  }
617 
618  strainSpaceHardeningVariables.add(helpVector);
619  // update plastic multipliers
620  gamma.add(dgamma);
621 
622  elasticStrainVectorR = totalStrain;
623  elasticStrainVectorR.subtract(plasticStrainVectorR);
624  // stress vector in full form due to computational convinience
625  this->computeTrialStressIncrement(fullStressVector, gp, elasticStrainVectorR, tStep);
626  this->computeStressSpaceHardeningVars(fullStressSpaceHardeningVars, gp, strainSpaceHardeningVariables);
627 
628  elastic = 1;
629  for ( i = 1; i <= this->nsurf; i++ ) {
630  yieldValue = this->computeYieldValueAt(gp, i, fullStressVector, fullStressSpaceHardeningVars);
631  // check for end of iteration
632  if ( yieldValue > YIELD_TOL ) {
633  elastic = 0;
634  }
635  }
636 
637  if ( elastic ) {
638  break;
639  }
640 
641  // increment iteration count
642  nIterations++;
643 
644  if ( nIterations > PLASTIC_MATERIAL_MAX_ITERATIONS ) {
645  char buff [ 200 ], buff1 [ 150 ];
646  sprintf(buff, "GiveRealStressVector: local equlibrium not reached in %d iterations", PLASTIC_MATERIAL_MAX_ITERATIONS);
647  sprintf( buff1, "Element %d, gp %d, continuing", gp->giveElement()->giveNumber(), gp->giveNumber() );
648  OOFEM_WARNING(buff, buff1);
649  //nIterations = 0; goto huhu;
650  break;
651  }
652  } while ( 1 );
653  }
654 
655  printf(" (%d iterations)", nIterations);
656  answer = fullStressVector;
657 }
658 
659 
660 void
662  GaussPoint *gp,
663  const FloatArray &fullStressVector,
664  const FloatArray &fullStressSpaceHardeningVars)
665 {
666  /*
667  * Computes gradient vector R in reduced form.
668  * R = {\der{\sigma}{f_{n+1}}, \der{q}{f_{n+1}}}^T
669  * Gradient vector contains partial derivatives of yield function
670  * with respect to stresses and with respect to
671  * strain space hardening variables
672  *
673  * Note: variables with R postfix are in reduced stress-space.
674  */
675  FloatArray stressGradient, stressGradientR;
676  FloatArray stressSpaceHardVarGradient;
677  int i, isize, size;
678 
679  this->computeStressGradientVector(stressGradient, ftype, isurf, gp, fullStressVector,
680  fullStressSpaceHardeningVars);
681 
682  StructuralMaterial :: giveReducedSymVectorForm( stressGradientR, stressGradient, gp->giveMaterialMode() );
683 
684  this->computeStressSpaceHardeningVarsReducedGradient(stressSpaceHardVarGradient, ftype, isurf, gp,
685  fullStressVector, fullStressSpaceHardeningVars);
686 
687  isize = stressGradientR.giveSize();
688  if ( stressSpaceHardVarGradient.isNotEmpty() ) {
689  size = isize + stressSpaceHardVarGradient.giveSize();
690  } else {
691  size = isize;
692  }
693 
694  answer.resize(size);
695  for ( i = 1; i <= isize; i++ ) {
696  answer.at(i) = stressGradientR.at(i);
697  }
698 
699  for ( i = isize + 1; i <= size; i++ ) {
700  answer.at(i) = stressSpaceHardVarGradient.at(i - isize);
701  }
702 }
703 
704 
705 void
707  const IntArray &activeConditionMap, const FloatArray &plasticStrainVectorR,
708  const FloatArray &strainSpaceHardeningVariables,
709  std :: vector< FloatArray > &gradientVectorR)
710 {
711  /* Computes Residual vector for closest point projection algorithm */
712 
713  FloatArray oldPlasticStrainVectorR, oldStrainSpaceHardeningVariables;
714  int i, j, isize, size;
715  MPlasticMaterialStatus *status = static_cast< MPlasticMaterialStatus * >( this->giveStatus(gp) );
716 
717  isize = plasticStrainVectorR.giveSize();
718  size = gradientVectorR [ 0 ].giveSize();
719 
720  answer.resize(size);
721  oldPlasticStrainVectorR = status->givePlasticStrainVector();
722  oldStrainSpaceHardeningVariables = status->giveStrainSpaceHardeningVars();
723 
724  for ( i = 1; i <= isize; i++ ) {
725  answer.at(i) = oldPlasticStrainVectorR.at(i) - plasticStrainVectorR.at(i);
726  }
727 
728  for ( i = isize + 1; i <= size; i++ ) {
729  answer.at(i) = oldStrainSpaceHardeningVariables.at(i - isize) - strainSpaceHardeningVariables.at(i - isize);
730  }
731 
732  for ( i = 1; i <= this->nsurf; i++ ) {
733  if ( activeConditionMap.at(i) ) {
734  for ( j = 1; j <= size; j++ ) {
735  answer.at(j) += gamma.at(i) * gradientVectorR [ i - 1 ].at(j);
736  }
737  }
738  }
739 }
740 
741 
742 void
744  const FloatArray &elasticStrainVectorR,
745  TimeStep *tStep)
746 {
747  /* Computes the full trial elastic stress vector */
748 
749  FloatMatrix de;
750  FloatArray reducedAnswer;
751 
752  this->computeReducedElasticModuli(de, gp, tStep);
753  //this->giveLinearElasticMaterial()->giveCharacteristicMatrix(de, TangentStiffness, gp, tStep);
754  reducedAnswer.beProductOf(de, elasticStrainVectorR);
755  StructuralMaterial :: giveFullSymVectorForm( answer, reducedAnswer, gp->giveMaterialMode() );
756 }
757 
758 
759 void
761  GaussPoint *gp,
762  const FloatMatrix &elasticModuliInverse,
763  const FloatMatrix &hardeningModuliInverse,
764  const FloatArray &gamma,
765  const IntArray &activeConditionMap,
766  const FloatArray &fullStressVector,
767  const FloatArray &fullStressSpaceHardeningVars)
768 {
769  /* returns consistent moduli in reduced form.
770  * Note: elasticModuli and hardeningModuli will be inverted
771  *
772  * stressVector in full form
773  */
774  FloatMatrix gradientMatrix;
775  FloatMatrix helpInverse;
776  int i, j, isize, size;
777 
778  isize = elasticModuliInverse.giveNumberOfRows();
779  if ( this->hasHardening() ) {
780  size = isize + hardeningModuliInverse.giveNumberOfRows();
781  } else {
782  size = isize;
783  }
784 
785  // assemble consistent moduli
786  helpInverse.resize(size, size);
787  helpInverse.zero();
788  for ( i = 1; i <= nsurf; i++ ) {
789  if ( activeConditionMap.at(i) ) {
790  // ask gradientMatrix
791  this->computeReducedGradientMatrix(gradientMatrix, i, gp, fullStressVector,
792  fullStressSpaceHardeningVars);
793 
794  gradientMatrix.times( gamma.at(i) );
795  helpInverse.add(gradientMatrix);
796  }
797  }
798 
799  for ( i = 1; i <= isize; i++ ) {
800  for ( j = 1; j <= isize; j++ ) {
801  helpInverse.at(i, j) += elasticModuliInverse.at(i, j);
802  }
803  }
804 
805  if ( hardeningModuliInverse.giveNumberOfRows() > 0 ) {
806  for ( i = isize + 1; i <= size; i++ ) {
807  for ( j = isize + 1; j <= size; j++ ) {
808  helpInverse.at(i, j) += hardeningModuliInverse.at(i - isize, j - isize);
809  }
810  }
811  }
812 
813  if ( hasHardening() && ( hardeningModuliInverse.giveNumberOfRows() == 0 ) ) {
814  FloatMatrix help;
815  help.beInverseOf(helpInverse);
816  answer.resize( isize + fullStressSpaceHardeningVars.giveSize(), isize + fullStressSpaceHardeningVars.giveSize() );
817  answer.zero();
818  for ( i = 1; i <= size; i++ ) {
819  for ( j = 1; j <= size; j++ ) {
820  answer.at(i, j) = help.at(i, j);
821  }
822  }
823  } else {
824  answer.beInverseOf(helpInverse);
825  }
826 }
827 
828 // ----------------------------------------------------------------------------//
829 
830 
831 void
833  MatResponseMode mode,
834  GaussPoint *gp,
835  TimeStep *tStep)
836 {
837  //
838  // returns receiver material matrix for given reached state
839  // described by temp variables.
840  //
841 
842  FloatMatrix consistentModuli, elasticModuli, hardeningModuli;
843  FloatMatrix elasticModuliInverse, hardeningModuliInverse;
844  FloatMatrix gmatInv, gmat, gradMat, helpMtrx, helpMtrx2, nmat, sbm1;
845  FloatArray stressVector, fullStressVector;
846  FloatArray stressSpaceHardeningVars;
847  FloatArray strainSpaceHardeningVariables;
848  std :: vector< FloatArray > yieldGradVec(this->nsurf), loadGradVec(this->nsurf), * yieldGradVecPtr, * loadGradVecPtr;
849 
850  IntArray activeConditionMap;
851  FloatArray gamma;
852  int size, sizeR, i, j, iindx, actSurf = 0;
853 
854  MPlasticMaterialStatus *status = static_cast< MPlasticMaterialStatus * >( this->giveStatus(gp) );
855 
856  if ( this->plType == associatedPT ) {
857  yieldGradVecPtr = loadGradVecPtr = & yieldGradVec;
858  } else {
859  yieldGradVecPtr = & yieldGradVec;
860  loadGradVecPtr = & loadGradVec;
861  }
862 
863  // ask for plastic consistency parameter
864  gamma = status->giveTempGamma();
865  activeConditionMap = status->giveTempActiveConditionMap();
866  //
867  // check for elastic cases
868  //
870  this->giveLinearElasticMaterial()->giveStiffnessMatrix(answer, TangentStiffness, gp, tStep);
871  return;
872  }
873 
874  //
875  // plastic case
876  //
877  // determine number of active surfaces
878  for ( i = 1; i <= nsurf; i++ ) {
879  if ( activeConditionMap.at(i) ) {
880  actSurf++;
881  }
882  }
883 
884  // compute elastic moduli and its inverse
885  this->computeReducedElasticModuli(elasticModuli, gp, tStep);
886  elasticModuliInverse.beInverseOf(elasticModuli);
887  sizeR = elasticModuliInverse.giveNumberOfRows();
888 
889  // compute hardening moduli and its inverse
890  this->computeHardeningReducedModuli(hardeningModuli, gp, strainSpaceHardeningVariables, tStep);
891  if ( hardeningModuli.giveNumberOfRows() ) {
892  hardeningModuliInverse.beInverseOf(hardeningModuli);
893  } else {
894  hardeningModuliInverse.clear();
895  }
896 
897  stressVector = status->giveStressVector();
898  StructuralMaterial :: giveFullSymVectorForm( fullStressVector, stressVector, gp->giveMaterialMode() );
899  strainSpaceHardeningVariables = status->giveStrainSpaceHardeningVars();
900  this->computeStressSpaceHardeningVars(stressSpaceHardeningVars, gp, strainSpaceHardeningVariables);
901 
902 
903  //
904  // compute consistent moduli
905  //
906  this->computeAlgorithmicModuli(consistentModuli, gp, elasticModuliInverse, hardeningModuliInverse, gamma,
907  activeConditionMap, fullStressVector, stressSpaceHardeningVars);
908 
909  //computee gmatInv
910  for ( i = 1; i <= nsurf; i++ ) {
911  this->computeGradientVector(yieldGradVec [ i - 1 ], yieldFunction, i, gp, fullStressVector, stressSpaceHardeningVars);
912  }
913 
914  if ( this->plType == nonassociatedPT ) {
915  for ( i = 1; i <= nsurf; i++ ) {
916  this->computeGradientVector(loadGradVec [ i - 1 ], loadFunction, i, gp, fullStressVector, stressSpaceHardeningVars);
917  }
918  }
919 
920 
921  size = consistentModuli.giveNumberOfColumns();
922  gradMat.resize(actSurf, size);
923  for ( i = 1; i <= nsurf; i++ ) {
924  if ( ( iindx = activeConditionMap.at(i) ) ) {
925  for ( j = 1; j <= size; j++ ) {
926  gradMat.at(iindx, j) = ( * yieldGradVecPtr ) [ i - 1 ].at(j);
927  }
928  }
929  }
930 
931  if ( this->plType == associatedPT ) {
932  helpMtrx.resize(actSurf, size);
933  helpMtrx.beProductTOf(consistentModuli, gradMat);
934  gmatInv.beProductOf(gradMat, helpMtrx);
935  gmat.beInverseOf(gmatInv);
936 
937  nmat.resize(sizeR, size);
938  sbm1.beSubMatrixOf(consistentModuli, 1, sizeR, 1, size);
939  nmat.beProductTOf(sbm1, gradMat);
940  helpMtrx.beProductOf(nmat, gmat);
941  sbm1.beSubMatrixOf(consistentModuli, 1, size, 1, sizeR);
942  nmat.beProductOf(gradMat, sbm1);
943 
944  helpMtrx2.beProductOf(helpMtrx, nmat);
945 
946  answer.beSubMatrixOf(consistentModuli, 1, sizeR, 1, sizeR);
947  answer.subtract(helpMtrx2);
948  } else {
949  FloatMatrix lgradMat(actSurf, size);
950  for ( i = 1; i <= nsurf; i++ ) {
951  if ( ( iindx = activeConditionMap.at(i) ) ) {
952  for ( j = 1; j <= size; j++ ) {
953  lgradMat.at(iindx, j) = ( * loadGradVecPtr ) [ i - 1 ].at(j);
954  }
955  }
956  }
957 
958  helpMtrx.resize(actSurf, size);
959  helpMtrx.beProductTOf(consistentModuli, lgradMat);
960  gmatInv.beProductOf(gradMat, helpMtrx);
961  gmat.beInverseOf(gmatInv);
962 
963  nmat.resize(sizeR, size);
964  sbm1.beSubMatrixOf(consistentModuli, 1, sizeR, 1, size);
965  nmat.beProductTOf(sbm1, lgradMat);
966  helpMtrx.beProductOf(nmat, gmat);
967  sbm1.beSubMatrixOf(consistentModuli, 1, size, 1, sizeR);
968  nmat.beProductOf(gradMat, sbm1);
969  helpMtrx2.beProductOf(helpMtrx, nmat);
970 
971  answer.beSubMatrixOf(consistentModuli, 1, sizeR, 1, sizeR);
972  answer.subtract(helpMtrx2);
973  }
974 }
975 
976 
977 void
979  MatResponseMode mode,
980  GaussPoint *gp,
981  TimeStep *tStep)
982 {
983  //
984  // returns receiver material matrix for given reached state
985  // described by temp variables.
986  //
987 
988  FloatMatrix dmat, elasticModuli, hardeningModuli;
989  FloatMatrix gmatInv, gmat, gradMat, helpMtrx, helpMtrx2;
990  FloatArray stressVector, fullStressVector;
991  FloatArray stressSpaceHardeningVars;
992  FloatArray strainSpaceHardeningVariables;
993  std :: vector< FloatArray > yieldGradVec(this->nsurf), loadGradVec(this->nsurf);
994 
995  IntArray activeConditionMap;
996  int size, sizeR, i, j, iindx, actSurf = 0;
997 
998  MPlasticMaterialStatus *status = static_cast< MPlasticMaterialStatus * >( this->giveStatus(gp) );
999 
1000  // ask for plastic consistency parameter
1001  activeConditionMap = status->giveTempActiveConditionMap();
1002  //
1003  // check for elastic cases
1004  //
1006  this->giveLinearElasticMaterial()->giveStiffnessMatrix(answer, TangentStiffness, gp, tStep);
1007  return;
1008  }
1009 
1010  //
1011  // plastic case
1012  //
1013  // determine number of active surfaces
1014  for ( i = 1; i <= nsurf; i++ ) {
1015  if ( activeConditionMap.at(i) ) {
1016  actSurf++;
1017  }
1018  }
1019 
1020  // compute elastic moduli and its inverse
1021  this->computeReducedElasticModuli(elasticModuli, gp, tStep);
1022  sizeR = elasticModuli.giveNumberOfRows();
1023 
1024  // compute hardening moduli and its inverse
1025  this->computeHardeningReducedModuli(hardeningModuli, gp, strainSpaceHardeningVariables, tStep);
1026 
1027  stressVector = status->giveStressVector();
1028  StructuralMaterial :: giveFullSymVectorForm( fullStressVector, stressVector, gp->giveMaterialMode() );
1029  strainSpaceHardeningVariables = status->giveStrainSpaceHardeningVars();
1030  this->computeStressSpaceHardeningVars(stressSpaceHardeningVars, gp, strainSpaceHardeningVariables);
1031 
1032 
1033  //
1034  // compute consistent moduli
1035  //
1036  size = elasticModuli.giveNumberOfRows() + hardeningModuli.giveNumberOfRows();
1037  dmat.resize(size, size);
1038  dmat.zero();
1039 
1040  for ( i = 1; i <= sizeR; i++ ) {
1041  for ( j = 1; j <= sizeR; j++ ) {
1042  dmat.at(i, j) = elasticModuli.at(i, j);
1043  }
1044  }
1045 
1046  for ( i = sizeR + 1; i <= size; i++ ) {
1047  for ( j = sizeR + 1; j <= size; j++ ) {
1048  dmat.at(i, j) = hardeningModuli.at(i - sizeR, j - sizeR);
1049  }
1050  }
1051 
1052 
1053  //computee gmatInv
1054  for ( i = 1; i <= nsurf; i++ ) {
1055  this->computeGradientVector(yieldGradVec [ i - 1 ], yieldFunction, i, gp, fullStressVector, stressSpaceHardeningVars);
1056  }
1057 
1058  if ( this->plType == nonassociatedPT ) {
1059  for ( i = 1; i <= nsurf; i++ ) {
1060  this->computeGradientVector(loadGradVec [ i - 1 ], loadFunction, i, gp, fullStressVector, stressSpaceHardeningVars);
1061  }
1062  }
1063 
1064  gradMat.resize(actSurf, size);
1065  for ( i = 1; i <= nsurf; i++ ) {
1066  if ( ( iindx = activeConditionMap.at(i) ) ) {
1067  for ( j = 1; j <= size; j++ ) {
1068  gradMat.at(iindx, j) = yieldGradVec [ i - 1 ].at(j);
1069  }
1070  }
1071  }
1072 
1073  if ( this->plType == associatedPT ) {
1074  helpMtrx.resize(actSurf, size);
1075  helpMtrx.beProductOf(gradMat, dmat);
1076  gmatInv.beProductTOf(helpMtrx, gradMat);
1077  gmat.beInverseOf(gmatInv);
1078 
1079  helpMtrx.beProductOf(gradMat, elasticModuli);
1080  helpMtrx2.beProductOf(gmat, helpMtrx);
1081  answer.beTProductOf(helpMtrx, helpMtrx2);
1082  answer.negated();
1083  answer.add(elasticModuli);
1084  } else {
1085  FloatMatrix lgradMat(actSurf, size);
1086  for ( i = 1; i <= nsurf; i++ ) {
1087  if ( ( iindx = activeConditionMap.at(i) ) ) {
1088  for ( j = 1; j <= size; j++ ) {
1089  lgradMat.at(iindx, j) = loadGradVec [ i - 1 ].at(j);
1090  }
1091  }
1092  }
1093 
1094  helpMtrx.resize(actSurf, size);
1095  helpMtrx.beProductOf(gradMat, dmat);
1096  gmatInv.beProductTOf(helpMtrx, lgradMat);
1097  gmat.beInverseOf(gmatInv);
1098 
1099  helpMtrx.beProductOf(gradMat, elasticModuli);
1100  helpMtrx2.beProductOf(gmat, helpMtrx);
1101  helpMtrx.beProductTOf(elasticModuli, lgradMat);
1102  answer.beTProductOf(helpMtrx, helpMtrx2);
1103  answer.negated();
1104  answer.add(elasticModuli);
1105  }
1106 }
1107 
1108 
1109 void
1111  GaussPoint *gp, FloatMatrix &elasticModuliInverse,
1112  FloatMatrix &hardeningModuliInverse)
1113 {
1114  //
1115  // assembles diagonal moduli from elasticModuliInverse and hardeningModuliInverse
1116  //
1117  int size1, size2, i, j;
1118 
1119  size1 = elasticModuliInverse.giveNumberOfRows();
1120  if ( hardeningModuliInverse.giveNumberOfRows() ) {
1121  size2 = size1 + hardeningModuliInverse.giveNumberOfRows();
1122  } else {
1123  size2 = size1;
1124  }
1125 
1126  answer.resize(size2, size2);
1127  answer.zero();
1128 
1129  for ( i = 1; i <= size1; i++ ) {
1130  for ( j = 1; j <= size1; j++ ) {
1131  answer.at(i, j) = elasticModuliInverse.at(i, j);
1132  }
1133  }
1134 
1135  for ( i = size1 + 1; i <= size2; i++ ) {
1136  for ( j = size1 + 1; j <= size2; j++ ) {
1137  answer.at(i, j) = hardeningModuliInverse.at(i - size1, j - size1);
1138  }
1139  }
1140 }
1141 
1142 
1143 void
1145  GaussPoint *gp,
1146  TimeStep *tStep)
1147 { /* Returns elastic moduli in reduced stress-strain space*/
1149  ElasticStiffness,
1150  gp, tStep);
1151 }
1152 
1153 
1154 // overloaded from structural material
1155 
1156 void
1158  MatResponseMode mode,
1159  GaussPoint *gp,
1160  TimeStep *tStep)
1161 //
1162 //
1163 //
1164 // computes full 3d constitutive matrix for case of 3d stress-strain state.
1165 // it returns elasto-plastic stiffness material matrix.
1166 // if strainIncrement == NULL a loading is assumed
1167 // for detailed description see (W.F.Chen: Plasticity in Reinforced Concrete, McGraw-Hill, 1982,
1168 // chapter 6.)
1169 //
1170 // if derived material would like to implement failure behaviour
1171 // it must redefine basic Give3dMaterialStiffnessMatrix function
1172 // in order to take possible failure (tension cracking) into account
1173 //
1174 //
1175 //
1176 {
1177  MaterialMode originalMode = gp->giveMaterialMode();
1178  if ( originalMode != _3dMat ) {
1179  OOFEM_ERROR("Different stressStrain mode encountered");
1180  }
1181 
1182  // we can force 3d response, and we obtain correct 3d tangent matrix,
1183  // but in fact, stress integration algorithm will not work
1184  // because in stress integration algorithm we are unable to recognize
1185  // which reduction from 3d case should be performed to obtain correct result.
1186  // so for new stressStrain state, instead of programming 3d reduction,
1187  // you should enhance imposeConstraints functions for ne state, and
1188  // then programming simple inteface function for you stressstrain state
1189  // calling GiveMaterailStiffenssMatrix, which imposes constrains correctly.
1190  if ( mode == ElasticStiffness ) {
1191  this->giveLinearElasticMaterial()->give3dMaterialStiffnessMatrix(answer, mode, gp, tStep);
1192  } else if ( rmType == mpm_ClosestPoint ) {
1193  this->giveConsistentStiffnessMatrix(answer, mode, gp, tStep);
1194  } else {
1195  this->giveElastoPlasticStiffnessMatrix(answer, mode, gp, tStep);
1196  }
1197 }
1198 
1199 
1200 void
1202  MatResponseMode mode,
1203  GaussPoint *gp,
1204  TimeStep *tStep)
1205 
1206 //
1207 // returns receiver's 2dPlaneStressMtrx
1208 // (2dPlaneStres ==> sigma_z = tau_xz = tau_yz = 0.)
1209 //
1210 // standard method from Material Class overloaded, because no inversion is needed.
1211 // the reduction from 3d case will not work
1212 // this implementation should be faster.
1213 {
1214  if ( mode == ElasticStiffness ) {
1215  this->giveLinearElasticMaterial()->giveStiffnessMatrix(answer, mode, gp, tStep);
1216  } else if ( rmType == mpm_ClosestPoint ) {
1217  this->giveConsistentStiffnessMatrix(answer, mode, gp, tStep);
1218  } else {
1219  this->giveElastoPlasticStiffnessMatrix(answer, mode, gp, tStep);
1220  }
1221 }
1222 
1223 
1224 void
1226  MatResponseMode mode,
1227  GaussPoint *gp,
1228  TimeStep *tStep)
1229 
1230 //
1231 // return receiver's 2dPlaneStrainMtrx constructed from
1232 // general 3dMatrialStiffnessMatrix
1233 // (2dPlaneStrain ==> eps_z = gamma_xz = gamma_yz = 0.)
1234 //
1235 {
1236  if ( mode == ElasticStiffness ) {
1237  this->giveLinearElasticMaterial()->giveStiffnessMatrix(answer, mode, gp, tStep);
1238  } else if ( rmType == mpm_ClosestPoint ) {
1239  this->giveConsistentStiffnessMatrix(answer, mode, gp, tStep);
1240  } else {
1241  this->giveElastoPlasticStiffnessMatrix(answer, mode, gp, tStep);
1242  }
1243 }
1244 
1245 
1246 void
1248  MatResponseMode mode,
1249  GaussPoint *gp,
1250  TimeStep *tStep)
1251 
1252 //
1253 // returns receiver's 1dMaterialStiffnessMAtrix
1254 // (1d case ==> sigma_y = sigma_z = tau_yz = tau_zx = tau_xy = 0.)
1255 {
1256  if ( mode == ElasticStiffness ) {
1257  this->giveLinearElasticMaterial()->giveStiffnessMatrix(answer, mode, gp, tStep);
1258  } else if ( rmType == mpm_ClosestPoint ) {
1259  this->giveConsistentStiffnessMatrix(answer, mode, gp, tStep);
1260  } else {
1261  this->giveElastoPlasticStiffnessMatrix(answer, mode, gp, tStep);
1262  }
1263 }
1264 
1265 
1266 void
1268  MatResponseMode mode,
1269  GaussPoint *gp,
1270  TimeStep *tStep)
1271 //
1272 // returns receiver's 2dBeamLayerStiffMtrx.
1273 // (2dPlaneStres ==> sigma_z = tau_xz = tau_yz = 0.)
1274 //
1275 // standard method from Material Class overloaded, because no inversion is needed.
1276 // the reduction from 3d case will not work
1277 // this implementation should be faster.
1278 {
1279  if ( mode == ElasticStiffness ) {
1280  this->giveLinearElasticMaterial()->giveStiffnessMatrix(answer, mode, gp, tStep);
1281  } else if ( rmType == mpm_ClosestPoint ) {
1282  this->giveConsistentStiffnessMatrix(answer, mode, gp, tStep);
1283  } else {
1284  this->giveElastoPlasticStiffnessMatrix(answer, mode, gp, tStep);
1285  }
1286 }
1287 
1288 
1289 void
1291  MatResponseMode mode,
1292  GaussPoint *gp,
1293  TimeStep *tStep)
1294 //
1295 // returns receiver's 2dPlateLayerMtrx
1296 // (2dPlaneStres ==> sigma_z = tau_xz = tau_yz = 0.)
1297 //
1298 // standard method from Material Class overloaded, because no inversion is needed.
1299 // the reduction from 3d case will not work
1300 // this implementation should be faster.
1301 {
1302  if ( mode == ElasticStiffness ) {
1303  this->giveLinearElasticMaterial()->giveStiffnessMatrix(answer, mode, gp, tStep);
1304  } else if ( rmType == mpm_ClosestPoint ) {
1305  this->giveConsistentStiffnessMatrix(answer, mode, gp, tStep);
1306  } else {
1307  this->giveElastoPlasticStiffnessMatrix(answer, mode, gp, tStep);
1308  }
1309 }
1310 
1311 
1312 void
1314  MatResponseMode mode,
1315  GaussPoint *gp,
1316  TimeStep *tStep)
1317 //
1318 // returns receiver's Fiber
1319 // (1dFiber ==> sigma_y = sigma_z = tau_yz = 0.)
1320 //
1321 // standard method from Material Class overloaded, because no inversion is needed.
1322 // the reduction from 3d case will not work
1323 // this implementation should be faster.
1324 {
1325  if ( mode == ElasticStiffness ) {
1326  this->giveLinearElasticMaterial()->giveStiffnessMatrix(answer, mode, gp, tStep);
1327  } else if ( rmType == mpm_ClosestPoint ) {
1328  this->giveConsistentStiffnessMatrix(answer, mode, gp, tStep);
1329  } else {
1330  this->giveElastoPlasticStiffnessMatrix(answer, mode, gp, tStep);
1331  }
1332 }
1333 
1334 
1335 int
1337 {
1338  MPlasticMaterialStatus *status = static_cast< MPlasticMaterialStatus * >( this->giveStatus(gp) );
1339  if ( type == IST_PlasticStrainTensor ) {
1340  const FloatArray &ep = status->givePlasticStrainVector();
1343  return 1;
1344  } else if ( type == IST_PrincipalPlasticStrainTensor ) {
1345  FloatArray st;
1346 
1347  const FloatArray &s = status->givePlasticStrainVector();
1350 
1351  this->computePrincipalValues(answer, st, principal_strain);
1352  return 1;
1353  } else {
1354  return StructuralMaterial :: giveIPValue(answer, gp, type, tStep);
1355  }
1356 }
1357 
1358 
1360  StructuralMaterialStatus(n, d, g), plasticStrainVector(), tempPlasticStrainVector(),
1361  strainSpaceHardeningVarsVector(statusSize), tempStrainSpaceHardeningVarsVector(statusSize),
1362  gamma(),
1363  tempGamma()
1364 {
1366 }
1367 
1368 
1370 { }
1371 
1372 
1373 void
1375 {
1377  fprintf(file, "status { ");
1380  fprintf(file, " Yielding, ");
1381  } else {
1382  fprintf(file, " Unloading, ");
1383  }
1384 
1385  fprintf(file, " plastic strains ");
1386  for ( auto &val : plasticStrainVector ) {
1387  fprintf( file, " %.4e", val );
1388  }
1389 
1391  fprintf(file, ", strain space hardening vars ");
1392  for ( auto &val : strainSpaceHardeningVarsVector ) {
1393  fprintf( file, " %.4e", val );
1394  }
1395  }
1396  }
1397 
1398  fprintf(file, "}\n");
1399 }
1400 
1401 
1403 //
1404 // initializes temp variables according to variables form previous equlibrium state.
1405 //
1406 {
1408 
1409  if ( plasticStrainVector.giveSize() == 0 ) {
1412  }
1413 
1415 
1417 
1418 
1419  tempGamma = gamma;
1421 }
1422 
1423 
1424 void
1426 //
1427 // updates variables (nonTemp variables describing situation at previous equilibrium state)
1428 // after a new equilibrium state has been reached
1429 // temporary variables are having values corresponding to newly reched equilibrium.
1430 //
1431 {
1433 
1436 
1438  gamma = tempGamma;
1440 }
1441 
1442 
1445 //
1446 // saves full information stored in this Status
1447 // no temp variables stored
1448 //
1449 {
1450  contextIOResultType iores;
1451 
1452  // save parent class status
1453  if ( ( iores = StructuralMaterialStatus :: saveContext(stream, mode, obj) ) != CIO_OK ) {
1454  THROW_CIOERR(iores);
1455  }
1456 
1457  // write a raw data
1458  if ( ( iores = plasticStrainVector.storeYourself(stream) ) != CIO_OK ) {
1459  THROW_CIOERR(iores);
1460  }
1461 
1462  if ( ( iores = strainSpaceHardeningVarsVector.storeYourself(stream) ) != CIO_OK ) {
1463  THROW_CIOERR(iores);
1464  }
1465 
1466  if ( !stream.write(state_flag) ) {
1468  }
1469 
1470  if ( ( iores = gamma.storeYourself(stream) ) != CIO_OK ) {
1471  THROW_CIOERR(iores);
1472  }
1473 
1474  if ( ( iores = activeConditionMap.storeYourself(stream) ) != CIO_OK ) {
1475  THROW_CIOERR(iores);
1476  }
1477 
1478 
1479  return CIO_OK;
1480 }
1481 
1482 
1485 //
1486 // restores full information stored in stream to this Status
1487 //
1488 {
1489  contextIOResultType iores;
1490 
1491  // read parent class status
1492  if ( ( iores = StructuralMaterialStatus :: restoreContext(stream, mode, obj) ) != CIO_OK ) {
1493  THROW_CIOERR(iores);
1494  }
1495 
1496  if ( ( iores = plasticStrainVector.restoreYourself(stream) ) != CIO_OK ) {
1497  THROW_CIOERR(iores);
1498  }
1499 
1500  if ( ( iores = strainSpaceHardeningVarsVector.restoreYourself(stream) ) != CIO_OK ) {
1501  THROW_CIOERR(iores);
1502  }
1503 
1504  if ( !stream.read(state_flag) ) {
1506  }
1507 
1508  if ( ( iores = gamma.restoreYourself(stream) ) != CIO_OK ) {
1509  THROW_CIOERR(iores);
1510  }
1511 
1512  if ( ( iores = activeConditionMap.restoreYourself(stream) ) != CIO_OK ) {
1513  THROW_CIOERR(iores);
1514  }
1515 
1516 
1517  return CIO_OK; // return succes
1518 }
1519 } // end namespace oofem
static int giveSizeOfVoigtSymVector(MaterialMode mmode)
Returns the size of symmetric part of a reduced stress/strain vector according to given mode...
virtual void computeReducedGradientMatrix(FloatMatrix &answer, int isurf, GaussPoint *gp, const FloatArray &stressVector, const FloatArray &stressSpaceHardeningVars)=0
virtual void updateYourself(TimeStep *tStep)
Update equilibrium history variables according to temp-variables.
Definition: structuralms.C:96
InternalStateType
Type representing the physical meaning of element or constitutive model internal variable.
static void computePrincipalValues(FloatArray &answer, const FloatArray &s, stressStrainPrincMode mode)
Common functions for convenience.
MaterialMode giveMaterialMode()
Returns corresponding material mode of receiver.
Definition: gausspoint.h:191
contextIOResultType storeYourself(DataStream &stream) const
Stores array to output stream.
Definition: intarray.C:289
int giveNumberOfColumns() const
Returns number of columns of receiver.
Definition: floatmatrix.h:158
void subtract(const FloatArray &src)
Subtracts array src to receiver.
Definition: floatarray.C:258
virtual void computeStressSpaceHardeningVarsReducedGradient(FloatArray &answer, functType ftype, int isurf, GaussPoint *gp, const FloatArray &stressVector, const FloatArray &stressSpaceHardeningVars)=0
virtual void giveConsistentStiffnessMatrix(FloatMatrix &answer, MatResponseMode, GaussPoint *gp, TimeStep *tStep)
void letTempStrainVectorBe(const FloatArray &v)
Assigns tempStrainVector to given vector v.
Definition: structuralms.h:137
virtual MaterialStatus * giveStatus(GaussPoint *gp) const
Returns material status of receiver in given integration point.
Definition: material.C:244
virtual void computeHardeningReducedModuli(FloatMatrix &answer, GaussPoint *gp, const FloatArray &strainSpaceHardeningVariables, TimeStep *tStep)=0
GaussPoint * gp
Associated integration point.
Class and object Domain.
Definition: domain.h:115
FloatArray plasticStrainVector
Plastic strain vector.
const FloatArray & giveTempGamma()
For computing principal strains from engineering strains.
functType
Type that allows to distinguish between yield function and loading function.
bool solveForRhs(const FloatArray &b, FloatArray &answer, bool transpose=false)
Solves the system of linear equations .
Definition: floatmatrix.C:1112
void letTempPlasticStrainVectorBe(FloatArray v)
virtual void updateYourself(TimeStep *tStep)
Update equilibrium history variables according to temp-variables.
The purpose of DataStream abstract class is to allow to store/restore context to different streams...
Definition: datastream.h:54
contextIOResultType storeYourself(DataStream &stream) const
Definition: floatarray.C:872
virtual int giveIPValue(FloatArray &answer, GaussPoint *gp, InternalStateType type, TimeStep *tStep)
Returns the integration point corresponding value in Reduced form.
void zero()
Sets all component to zero.
Definition: intarray.C:52
double & at(int i)
Coefficient access function.
Definition: floatarray.h:131
void beSubMatrixOf(const FloatMatrix &src, int topRow, int bottomRow, int topCol, int bottomCol)
Assigns to the receiver the sub-matrix of another matrix.
Definition: floatmatrix.C:962
enum oofem::MPlasticMaterial::plastType plType
virtual contextIOResultType saveContext(DataStream &stream, ContextMode mode, void *obj=NULL)
Stores receiver state to output stream.
virtual void givePlaneStrainStiffMtrx(FloatMatrix &answer, MatResponseMode, GaussPoint *gp, TimeStep *tStep)
Method for computing plane strain stiffness matrix of receiver.
This class implements a structural material status information.
Definition: structuralms.h:65
void computeGradientVector(FloatArray &answer, functType ftype, int isurf, GaussPoint *gp, const FloatArray &fullStressVector, const FloatArray &fullStressSpaceHardeningVars)
void computeDiagModuli(FloatMatrix &answer, GaussPoint *gp, FloatMatrix &elasticModuliInverse, FloatMatrix &hardeningModuliInverse)
virtual void computeReducedElasticModuli(FloatMatrix &answer, GaussPoint *gp, TimeStep *tStep)
void computeResidualVector(FloatArray &answer, GaussPoint *gp, const FloatArray &gamma, const IntArray &activeConditionMap, const FloatArray &plasticStrainVectorR, const FloatArray &strainSpaceHardeningVariables, std::vector< FloatArray > &gradVec)
LinearElasticMaterial * linearElasticMaterial
Reference to bulk (undamaged) material.
void cuttingPlaneReturn(FloatArray &answer, IntArray &activeConditionMap, FloatArray &gamma, GaussPoint *gp, const FloatArray &totalStrain, FloatArray &plasticStrainR, FloatArray &strainSpaceHardeningVariables, TimeStep *tStep)
Element * giveElement()
Returns corresponding element to receiver.
Definition: gausspoint.h:188
General IO error.
void negated()
Changes sign of receiver values.
Definition: floatmatrix.C:1605
void giveStressDependentPartOfStrainVector(FloatArray &answer, GaussPoint *gp, const FloatArray &reducedStrainVector, TimeStep *tStep, ValueModeType mode)
Method for subtracting from reduced space strain vector its stress-independent parts (caused by tempe...
MaterialMode
Type representing material mode of integration point.
Definition: materialmode.h:89
LinearElasticMaterial * giveLinearElasticMaterial()
Returns reference to undamaged (bulk) material.
const FloatArray & givePlasticStrainVector() const
Returns the equilibrated strain vector.
virtual void printOutputAt(FILE *file, TimeStep *tStep)
Print receiver&#39;s output to given stream.
Definition: structuralms.C:73
Class implementing an array of integers.
Definition: intarray.h:61
int & at(int i)
Coefficient access function.
Definition: intarray.h:103
virtual int read(int *data, int count)=0
Reads count integer values into array pointed by data.
MatResponseMode
Describes the character of characteristic material matrix.
#define THROW_CIOERR(e)
Definition: contextioerr.h:61
const IntArray & giveTempActiveConditionMap()
virtual void give1dStressStiffMtrx(FloatMatrix &answer, MatResponseMode, GaussPoint *gp, TimeStep *tStep)
Method for computing 1d stiffness matrix of receiver.
virtual void computeStressGradientVector(FloatArray &answer, functType ftype, int isurf, GaussPoint *gp, const FloatArray &stressVector, const FloatArray &stressSpaceHardeningVars)=0
virtual int hasMaterialModeCapability(MaterialMode mode)
Tests if material supports material mode.
static void giveFullSymVectorForm(FloatArray &answer, const FloatArray &vec, MaterialMode matMode)
Converts the reduced unsymmetric Voigt vector (2nd order tensor) to full form.
virtual void givePlateLayerStiffMtrx(FloatMatrix &answer, MatResponseMode, GaussPoint *gp, TimeStep *tStep)
Method for computing 2d plate layer stiffness matrix of receiver.
FloatArray strainSpaceHardeningVarsVector
Strain space hardening variables.
void letTempStrainSpaceHardeningVarsVectorBe(FloatArray v)
virtual int write(const int *data, int count)=0
Writes count integer values from array pointed by data.
virtual int giveIPValue(FloatArray &answer, GaussPoint *gp, InternalStateType type, TimeStep *tStep)
Returns the integration point corresponding value in Reduced form.
const FloatArray & giveStrainSpaceHardeningVars() const
Returns the equilibrated hardening variable vector.
int nsurf
Number of yield surfaces.
virtual void giveStiffnessMatrix(FloatMatrix &answer, MatResponseMode mode, GaussPoint *gp, TimeStep *tStep)
Computes the stiffness matrix for giveRealStressVector of receiver in given integration point...
virtual void givePlaneStressStiffMtrx(FloatMatrix &answer, MatResponseMode, GaussPoint *gp, TimeStep *tStep)
Method for computing plane stress stiffness matrix of receiver.
static void giveReducedSymVectorForm(FloatArray &answer, const FloatArray &vec, MaterialMode matMode)
Converts the full unsymmetric Voigt vector (2nd order tensor) to reduced form.
double dotProduct(const FloatArray &x) const
Computes the dot product (or inner product) of receiver and argument.
Definition: floatarray.C:463
#define OOFEM_ERROR(...)
Definition: error.h:61
int giveNumber()
Returns number of receiver.
Definition: gausspoint.h:184
virtual int giveSizeOfReducedHardeningVarsVector(GaussPoint *) const
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 ...
void times(double f)
Multiplies receiver by factor f.
Definition: floatmatrix.C:1594
contextIOResultType restoreYourself(DataStream &stream)
Definition: floatarray.C:895
virtual void computeStressSpaceHardeningVars(FloatArray &answer, GaussPoint *gp, const FloatArray &strainSpaceHardeningVariables)=0
virtual void give2dBeamLayerStiffMtrx(FloatMatrix &answer, MatResponseMode, GaussPoint *gp, TimeStep *tStep)
Method for computing 2d beam layer stiffness matrix of receiver.
void beProductOf(const FloatMatrix &aMatrix, const FloatArray &anArray)
Receiver becomes the result of the product of aMatrix and anArray.
Definition: floatarray.C:676
virtual MaterialStatus * CreateStatus(GaussPoint *gp) const
Creates new copy of associated status and inserts it into given integration point.
virtual void give3dMaterialStiffnessMatrix(FloatMatrix &answer, MatResponseMode, GaussPoint *gp, TimeStep *tStep)
Computes full 3d material stiffness matrix at given integration point, time, respecting load history ...
void beTProductOf(const FloatMatrix &aMatrix, const FloatArray &anArray)
Receiver becomes the result of the product of aMatrix^T and anArray.
Definition: floatarray.C:708
double at(int i, int j) const
Coefficient access function.
Definition: floatmatrix.h:176
enum oofem::MPlasticMaterial::ReturnMappingAlgoType rmType
virtual contextIOResultType restoreContext(DataStream &stream, ContextMode mode, void *obj=NULL)
Restores the receiver state previously written in stream.
Definition: structuralms.C:157
This class implements associated Material Status to MPlasticMaterial.
contextIOResultType restoreYourself(DataStream &stream)
Restores array from image on stream.
Definition: intarray.C:305
void subtract(const FloatMatrix &a)
Subtracts matrix from the receiver.
Definition: floatmatrix.C:1084
Abstract base class representing a material status information.
Definition: matstatus.h:84
Class representing vector of real numbers.
Definition: floatarray.h:82
Implementation of matrix containing floating point numbers.
Definition: floatmatrix.h:94
void setTempActiveConditionMap(IntArray v)
void letTempStressVectorBe(const FloatArray &v)
Assigns tempStressVector to given vector v.
Definition: structuralms.h:135
IntArray activeConditionMap
Active set of yield functions (needed for algorithmic stiffness).
const FloatArray & giveStressVector() const
Returns the const pointer to receiver&#39;s stress vector.
Definition: structuralms.h:107
void closestPointReturn(FloatArray &answer, IntArray &activeConditionMap, FloatArray &gamma, GaussPoint *gp, const FloatArray &totalStrain, FloatArray &plasticStrainR, FloatArray &strainSpaceHardeningVariables, TimeStep *tStep)
MPlasticMaterial(int n, Domain *d)
double computeNorm() const
Computes the norm (or length) of the vector.
Definition: floatarray.C:840
void resize(int rows, int cols)
Checks size of receiver towards requested bounds.
Definition: floatmatrix.C:1358
void add(const FloatMatrix &a)
Adds matrix to the receiver.
Definition: floatmatrix.C:1023
void computeAlgorithmicModuli(FloatMatrix &answer, GaussPoint *gp, const FloatMatrix &elasticModuliInverse, const FloatMatrix &hardeningModuliInverse, const FloatArray &gamma, const IntArray &activeConditionMap, const FloatArray &fullStressVector, const FloatArray &fullStressSpaceHardeningVars)
int state_flag
Yield function status indicator.
void zero()
Zeroes all coefficients of receiver.
Definition: floatarray.C:658
virtual void giveElastoPlasticStiffnessMatrix(FloatMatrix &answer, MatResponseMode mode, GaussPoint *gp, TimeStep *tStep)
void beTProductOf(const FloatMatrix &a, const FloatMatrix &b)
Assigns to the receiver product of .
Definition: floatmatrix.C:367
virtual void printOutputAt(FILE *file, TimeStep *tStep)
Print receiver&#39;s output to given stream.
long ContextMode
Context mode (mask), defining the type of information written/read to/from context.
Definition: contextmode.h:43
Abstract base class for all "structural" constitutive models.
#define RES_TOL
virtual double computeYieldValueAt(GaussPoint *gp, int isurf, const FloatArray &stressVector, const FloatArray &stressSpaceHardeningVars)=0
void beProductTOf(const FloatMatrix &a, const FloatMatrix &b)
Assigns to the receiver product of .
Definition: floatmatrix.C:397
void zero()
Zeroes all coefficient of receiver.
Definition: floatmatrix.C:1326
virtual void initTempStatus()
Initializes the temporary internal variables, describing the current state according to previously re...
Definition: structuralms.C:108
Domain * giveDomain() const
Definition: femcmpnn.h:100
FloatArray gamma
Consistency parameter values (needed for algorithmic stiffness).
void beProductOf(const FloatMatrix &a, const FloatMatrix &b)
Assigns to the receiver product of .
Definition: floatmatrix.C:337
const FloatArray & giveTempStressVector() const
Returns the const pointer to receiver&#39;s temporary stress vector.
Definition: structuralms.h:117
int giveSize() const
Returns the size of receiver.
Definition: floatarray.h:218
virtual contextIOResultType saveContext(DataStream &stream, ContextMode mode, void *obj=NULL)
Stores receiver state to output stream.
Definition: structuralms.C:133
#define PLASTIC_MATERIAL_MAX_ITERATIONS
the oofem namespace is to define a context or scope in which all oofem names are defined.
bool isNotEmpty() const
Returns true if receiver is not empty.
Definition: floatarray.h:220
void clear()
Sets size of receiver to be an empty matrix. It will have zero rows and zero columns size...
Definition: floatmatrix.h:516
virtual void giveFiberStiffMtrx(FloatMatrix &answer, MatResponseMode, GaussPoint *gp, TimeStep *tStep)
Method for computing 1d fiber stiffness matrix of receiver.
virtual void initTempStatus(GaussPoint *gp)
Initializes temporary variables stored in integration point status at the beginning of new time step...
Definition: material.C:267
int giveNumber() const
Definition: femcmpnn.h:107
void beInverseOf(const FloatMatrix &src)
Modifies receiver to become inverse of given parameter.
Definition: floatmatrix.C:835
#define YIELD_TOL
int giveNumberOfRows() const
Returns number of rows of receiver.
Definition: floatmatrix.h:156
virtual contextIOResultType restoreContext(DataStream &stream, ContextMode mode, void *obj=NULL)
Restores the receiver state previously written in stream.
Class representing integration point in finite element program.
Definition: gausspoint.h:93
virtual void giveRealStressVector(FloatArray &answer, GaussPoint *, const FloatArray &, TimeStep *)
Computes the real stress vector for given total strain and integration point.
#define OOFEM_WARNING(...)
Definition: error.h:62
Class representing solution step.
Definition: timestep.h:80
virtual void computeTrialStressIncrement(FloatArray &answer, GaussPoint *gp, const FloatArray &strainIncrement, TimeStep *tStep)
void add(const FloatArray &src)
Adds array src to receiver.
Definition: floatarray.C:156
virtual void initTempStatus()
Initializes the temporary internal variables, describing the current state according to previously re...
MPlasticMaterialStatus(int n, Domain *d, GaussPoint *g, int statusSize)
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