OOFEM  2.4
OOFEM.org - Object Oriented Finite Element Solver
j2mplasticmaterial.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 "j2mplasticmaterial.h"
37 #include "gausspoint.h"
38 #include "floatmatrix.h"
39 #include "floatarray.h"
40 #include "intarray.h"
41 #include "mathfem.h"
42 #include "classfactory.h"
43 
44 namespace oofem {
45 REGISTER_Material(J2MPlasticMaterial);
46 
48 {
49  //
50  // constructor
51  //
54  this->nsurf = 1;
55 }
56 
58 {
59  //
60  // destructor
61  //
62 }
63 
66 {
67  IRResultType result; // Required by IR_GIVE_FIELD macro
68  double value;
69 
71  if ( result != IRRT_OK ) return result;
73  if ( result != IRRT_OK ) return result;
74 
76  k = value / sqrt(3.0);
77 
78  kinematicModuli = 0.0;
80 
81  isotropicModuli = 0.0;
83 
84  if ( fabs(kinematicModuli) > 1.e-12 ) {
86  }
87 
88  if ( fabs(isotropicModuli) > 1.e-12 ) {
90  }
91 
92  int rma = 0;
94  if ( rma == 0 ) {
95  this->rmType = mpm_ClosestPoint;
96  } else {
97  this->rmType = mpm_CuttingPlane;
98  }
99 
100  return IRRT_OK;
101 }
102 
103 
106 {
107  return new MPlasticMaterialStatus(1, this->giveDomain(), gp, this->giveSizeOfReducedHardeningVarsVector(gp));
108 }
109 
110 void
112  const FloatArray &strainSpaceHardeningVariables)
113 {
114  // in full stress strain space
115  int count = 0, size = this->giveSizeOfFullHardeningVarsVector(), isize, rSize;
116  IntArray mask;
117 
118  if ( !hasHardening() ) {
119  answer.clear();
120  return;
121  }
122 
123  answer.resize(size);
125  isize = mask.giveSize();
126  rSize = this->giveSizeOfReducedHardeningVarsVector(gp);
127 
128  /* kinematic hardening variables are first */
129  if ( this->kinematicHardeningFlag ) {
130  for ( int i = 1; i <= isize; i++ ) {
131  // to be consistent with equivalent plastic strain formulation
132  // we multiply by (sqrt(2.)*2./3.)
133  answer.at( mask.at(i) ) = ( sqrt(2.) * 2. / 3. ) * this->kinematicModuli * strainSpaceHardeningVariables.at(i);
134  }
135 
136  count = 6;
137  }
138 
139  if ( this->isotropicHardeningFlag ) {
140  answer.at(count + 1) = this->isotropicModuli *
141  strainSpaceHardeningVariables.at(rSize);
142  }
143 
144  answer.negated();
145 }
146 
147 
148 
149 double
151  const FloatArray &stressSpaceHardeningVars)
152 {
153  double f;
154  FloatArray helpVector, backStress;
155 
156  if ( this->kinematicHardeningFlag ) {
157  if ( stressVector.isNotEmpty() ) {
158  this->giveStressBackVector(backStress, stressSpaceHardeningVars);
159  helpVector = stressVector;
160  helpVector.add(backStress);
161  } else {
162  return -k;
163  }
164  } else {
165  helpVector = stressVector;
166  }
167 
168  f = sqrt( this->computeJ2InvariantAt(helpVector) );
169 
170  //if (this->kinematicHardeningFlag) delete helpVector;
171 
172  return f + sqrt(1. / 3.) * this->giveIsotropicHardeningVar(stressSpaceHardeningVars) - this->k;
173 }
174 
175 
176 void
178  const FloatArray &strainSpaceHardeningVariables,
179  TimeStep *tStep)
180 {
181  /* computes hardening moduli in reduced stress strain space (for kinematic back-stress)*/
182 
183  int size = this->giveSizeOfReducedHardeningVarsVector(gp);
184 
185  if ( !hasHardening() ) {
186  answer.clear();
187  return;
188  }
189 
190  answer.resize(size, size);
191  answer.zero();
192 
193  /* kinematic hardening variables are first */
194  if ( this->kinematicHardeningFlag ) {
196  for ( int i = 1; i <= ksize; i++ ) {
197  answer.at(i, i) = this->kinematicModuli;
198  }
199  }
200 
201  if ( this->isotropicHardeningFlag ) {
202  answer.at(size, size) = this->isotropicModuli;
203  }
204 }
205 
206 
207 void
209  const FloatArray &stressSpaceHardeningVars)
210 {
211  /* stress gradient of yield function in full stress - strain space */
212 
213  double f, ax, ay, az, sx, sy, sz;
214  FloatArray helpVector, backStress;
215 
216  answer.resize(6);
217  answer.zero();
218  if ( this->kinematicHardeningFlag ) {
219  if ( stressVector.isNotEmpty() ) {
220  this->giveStressBackVector(backStress, stressSpaceHardeningVars);
221  helpVector = stressVector;
222  helpVector.add(backStress);
223  } else {
224  return;
225  }
226  } else {
227  helpVector = stressVector;
228  }
229 
230  f = sqrt( this->computeJ2InvariantAt(helpVector) );
231  // check for yield value zero value
232  if ( fabs(f) < 1.e-6 ) {
233  return;
234  }
235 
236  ax = helpVector.at(1);
237  ay = helpVector.at(2);
238  az = helpVector.at(3);
239 
240  sx = ( 2. / 3. ) * ax - ( 1. / 3. ) * ay - ( 1. / 3. ) * az;
241  sy = ( 2. / 3. ) * ay - ( 1. / 3. ) * ax - ( 1. / 3. ) * az;
242  sz = ( 2. / 3. ) * az - ( 1. / 3. ) * ay - ( 1. / 3. ) * ax;
243 
244  answer.at(1) = 0.5 * sx / f;
245  answer.at(2) = 0.5 * sy / f;
246  answer.at(3) = 0.5 * sz / f;
247  answer.at(4) = helpVector.at(4) / f;
248  answer.at(5) = helpVector.at(5) / f;
249  answer.at(6) = helpVector.at(6) / f;
250 }
251 
252 void
254  const FloatArray &stressVector,
255  const FloatArray &stressSpaceHardeningVars)
256 {
257  /* computes stress space hardening gradient in reduced stress-strain space */
258 
259  int kcount = 0, size = this->giveSizeOfReducedHardeningVarsVector(gp);
260  //double f,ax,ay,az,sx,sy,sz;
261  FloatArray fullKinematicGradient, reducedKinematicGrad;
262 
263  if ( !hasHardening() ) {
264  answer.clear();
265  return;
266  }
267 
268  answer.resize(size);
269 
270  /* kinematic hardening variables first */
271  if ( this->kinematicHardeningFlag ) {
272  this->computeStressGradientVector(fullKinematicGradient, ftype, isurf, gp, stressVector, stressSpaceHardeningVars);
273  StructuralMaterial :: giveReducedSymVectorForm( reducedKinematicGrad, fullKinematicGradient, gp->giveMaterialMode() );
274 
275  kcount = reducedKinematicGrad.giveSize();
276  }
277 
278  if ( this->kinematicHardeningFlag ) {
279  for ( int i = 1; i <= kcount; i++ ) {
280  answer.at(i) = reducedKinematicGrad.at(i);
281  }
282  }
283 
284  if ( this->isotropicHardeningFlag ) {
285  answer.at(size) = sqrt(1. / 3.);
286  }
287 }
288 
289 
290 int
292 {
293  return ( this->kinematicHardeningFlag || this->isotropicHardeningFlag );
294 }
295 
296 
297 void
299  GaussPoint *gp,
300  const FloatArray &stressVector,
301  const FloatArray &stressSpaceHardeningVars)
302 {
303  int size;
304  int imask, jmask;
305  FloatArray helpVector, backStress, df(6);
306  IntArray mask;
307  double f, f32, f12, ax, ay, az;
308 
312 
313  answer.resize(size, size);
314  answer.zero();
315 
316 
317  if ( stressVector.giveSize() != 0 ) {
318  /* kinematic hardening variables first */
319  if ( this->kinematicHardeningFlag ) {
320  this->giveStressBackVector(backStress, stressSpaceHardeningVars);
321  helpVector = stressVector;
322  helpVector.add(backStress);
323  } else {
324  helpVector = stressVector;
325  }
326 
327  f = this->computeJ2InvariantAt(helpVector);
328  f12 = sqrt(f);
329  f32 = pow(f, 3. / 2.);
330 
331  ax = helpVector.at(1);
332  ay = helpVector.at(2);
333  az = helpVector.at(3);
334 
335  df.at(1) = ( 2. / 3. ) * ax - ( 1. / 3. ) * ay - ( 1. / 3. ) * az;
336  df.at(2) = ( 2. / 3. ) * ay - ( 1. / 3. ) * ax - ( 1. / 3. ) * az;
337  df.at(3) = ( 2. / 3. ) * az - ( 1. / 3. ) * ay - ( 1. / 3. ) * ax;
338  df.at(4) = 2. * helpVector.at(4);
339  df.at(5) = 2. * helpVector.at(5);
340  df.at(6) = 2. * helpVector.at(6);
341 
342  for ( int i = 1; i <= 3; i++ ) {
343  if ( ( imask = mask.at(i) ) == 0 ) {
344  continue;
345  }
346 
347  for ( int j = i; j <= 3; j++ ) {
348  if ( ( jmask = mask.at(j) ) == 0 ) {
349  continue;
350  }
351 
352  if ( i == j ) {
353  answer.at(imask, jmask) = -( 1. / 4. ) * ( 1. / f32 ) * df.at(i) * df.at(j) + 0.5 * ( 1. / f12 ) * ( 4. / 6 );
354  } else {
355  answer.at(imask, jmask) = -( 1. / 4. ) * ( 1. / f32 ) * df.at(i) * df.at(j) + 0.5 * ( 1. / f12 ) * ( -2. / 6 );
356  answer.at(jmask, imask) = -( 1. / 4. ) * ( 1. / f32 ) * df.at(i) * df.at(j) + 0.5 * ( 1. / f12 ) * ( -2. / 6 );
357  }
358  }
359  }
360 
361  for ( int i = 1; i <= 3; i++ ) {
362  if ( ( imask = mask.at(i) ) == 0 ) {
363  continue;
364  }
365 
366  for ( int j = 4; j <= 6; j++ ) {
367  if ( ( jmask = mask.at(j) ) == 0 ) {
368  continue;
369  }
370 
371  answer.at(imask, jmask) = -( 1. / 4. ) * ( 1. / f32 ) * df.at(i) * df.at(j);
372  answer.at(jmask, imask) = -( 1. / 4. ) * ( 1. / f32 ) * df.at(i) * df.at(j);
373  }
374  }
375 
376  for ( int i = 4; i <= 6; i++ ) {
377  if ( ( imask = mask.at(i) ) == 0 ) {
378  continue;
379  }
380 
381  for ( int j = i; j <= 6; j++ ) {
382  if ( ( jmask = mask.at(j) ) == 0 ) {
383  continue;
384  }
385 
386  if ( i == j ) {
387  answer.at(imask, jmask) = -( 1. / 4. ) * ( 1. / f32 ) * df.at(i) * df.at(j) + 0.5 * ( 1. / f12 ) * 2.;
388  } else {
389  answer.at(imask, jmask) = -( 1. / 4. ) * ( 1. / f32 ) * df.at(i) * df.at(j);
390  answer.at(jmask, imask) = -( 1. / 4. ) * ( 1. / f32 ) * df.at(i) * df.at(j);
391  }
392  }
393  }
394  }
395 
396  /* for isotropic hardening: the corresponding part of gradient matrix is zero valued */
397 }
398 
399 
400 void
402  GaussPoint *gp,
403  TimeStep *tStep)
404 {
405  /* Returns 3d elastic moduli */
406  this->giveLinearElasticMaterial()->give3dMaterialStiffnessMatrix(answer, ElasticStiffness, gp, tStep);
407 }
408 
409 
410 double
412 {
413  double answer;
414  double v1, v2, v3;
415 
416  if ( stressVector.isEmpty() ) {
417  return 0.0;
418  }
419 
420  v1 = ( ( stressVector.at(1) - stressVector.at(2) ) * ( stressVector.at(1) - stressVector.at(2) ) );
421  v2 = ( ( stressVector.at(2) - stressVector.at(3) ) * ( stressVector.at(2) - stressVector.at(3) ) );
422  v3 = ( ( stressVector.at(3) - stressVector.at(1) ) * ( stressVector.at(3) - stressVector.at(1) ) );
423 
424  answer = ( 1. / 6. ) * ( v1 + v2 + v3 ) + stressVector.at(4) * stressVector.at(4) +
425  stressVector.at(5) * stressVector.at(5) + stressVector.at(6) * stressVector.at(6);
426 
427  return answer;
428 }
429 
430 
431 int
433 {
434  /* Returns the size of hardening variables vector */
435  int size = 0;
436 
437  if ( kinematicHardeningFlag ) {
438  size += 6; /* size of full stress vector */
439  }
440 
441  if ( isotropicHardeningFlag ) {
442  size += 1; /* scalar value */
443  }
444 
445  return size;
446 }
447 
448 int
450 {
451  /* Returns the size of hardening variables vector */
452  int size = 0;
453 
454  if ( kinematicHardeningFlag ) {
456  }
457 
458  if ( isotropicHardeningFlag ) {
459  size += 1; /* scalar value */
460  }
461 
462  return size;
463 }
464 
465 
466 void
468  const FloatArray &stressSpaceHardeningVars)
469 {
470  /* returns part of hardening vector corresponding to kinematic hardening */
471  if ( this->kinematicHardeningFlag ) {
472  answer.resize(6);
473  for ( int i = 1; i <= 6; i++ ) {
474  answer.at(i) = stressSpaceHardeningVars.at(i);
475  }
476 
477  return;
478  }
479 
480  answer.clear();
481 }
482 
483 
484 double
486 {
487  /* returns value in hardening vector corresponding to isotropic hardening */
488  if ( !isotropicHardeningFlag ) {
489  return 0.;
490  } else if ( this->kinematicHardeningFlag ) {
491  return stressSpaceHardeningVars.at(7);
492  } else {
493  return stressSpaceHardeningVars.at(1);
494  }
495 }
496 } // 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 MaterialStatus * CreateStatus(GaussPoint *gp) const
Creates new copy of associated status and inserts it into given integration point.
MaterialMode giveMaterialMode()
Returns corresponding material mode of receiver.
Definition: gausspoint.h:191
static int giveVoigtSymVectorMask(IntArray &answer, MaterialMode mmode)
The same as giveVoigtVectorMask but returns a mask corresponding to a symmetric second order tensor...
J2MPlasticMaterial(int n, Domain *d)
Class and object Domain.
Definition: domain.h:115
virtual double computeYieldValueAt(GaussPoint *gp, int isurf, const FloatArray &stressVector, const FloatArray &stressSpaceHardeningVars)
functType
Type that allows to distinguish between yield function and loading function.
#define _IFT_J2MPlasticMaterial_khm
#define _IFT_J2MPlasticMaterial_ry
double & at(int i)
Coefficient access function.
Definition: floatarray.h:131
void clear()
Clears receiver (zero size).
Definition: floatarray.h:206
virtual IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description stored in input record.
double giveIsotropicHardeningVar(const FloatArray &stressSpaceHardeningVars)
LinearElasticMaterial * linearElasticMaterial
Reference to bulk (undamaged) material.
virtual IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description stored in input record.
LinearElasticMaterial * giveLinearElasticMaterial()
Returns reference to undamaged (bulk) material.
This class implements a general plastic material.
Class implementing an array of integers.
Definition: intarray.h:61
int & at(int i)
Coefficient access function.
Definition: intarray.h:103
#define _IFT_J2MPlasticMaterial_rma
int nsurf
Number of yield surfaces.
virtual IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description stored in input record.
static void giveReducedSymVectorForm(FloatArray &answer, const FloatArray &vec, MaterialMode matMode)
Converts the full unsymmetric Voigt vector (2nd order tensor) to reduced form.
virtual void computeReducedGradientMatrix(FloatMatrix &answer, int isurf, GaussPoint *gp, const FloatArray &stressVector, const FloatArray &stressSpaceHardeningVars)
This class implements an isotropic linear elastic material in a finite element problem.
virtual void computeStressGradientVector(FloatArray &answer, functType ftype, int isurf, GaussPoint *gp, const FloatArray &stressVector, const FloatArray &stressSpaceHardeningVars)
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 ...
bool isEmpty() const
Returns true if receiver is empty.
Definition: floatarray.h:222
double at(int i, int j) const
Coefficient access function.
Definition: floatmatrix.h:176
virtual void computeStressSpaceHardeningVarsReducedGradient(FloatArray &answer, functType ftype, int isurf, GaussPoint *gp, const FloatArray &stressVector, const FloatArray &stressSpaceHardeningVars)
enum oofem::MPlasticMaterial::ReturnMappingAlgoType rmType
This class implements associated Material Status to MPlasticMaterial.
Abstract base class representing a material status information.
Definition: matstatus.h:84
#define _IFT_J2MPlasticMaterial_ihm
Class representing vector of real numbers.
Definition: floatarray.h:82
virtual void computeStressSpaceHardeningVars(FloatArray &answer, GaussPoint *gp, const FloatArray &strainSpaceHardeningVariables)
Implementation of matrix containing floating point numbers.
Definition: floatmatrix.h:94
IRResultType
Type defining the return values of InputRecord reading operations.
Definition: irresulttype.h:47
virtual void compute3dElasticModuli(FloatMatrix &answer, GaussPoint *gp, TimeStep *tStep)
void resize(int rows, int cols)
Checks size of receiver towards requested bounds.
Definition: floatmatrix.C:1358
Class representing the general Input Record.
Definition: inputrecord.h:101
void zero()
Zeroes all coefficients of receiver.
Definition: floatarray.C:658
virtual void computeHardeningReducedModuli(FloatMatrix &answer, GaussPoint *gp, const FloatArray &strainSpaceHardeningVariables, TimeStep *tStep)
int giveSizeOfReducedHardeningVarsVector(GaussPoint *gp) const
static void giveInvertedVoigtVectorMask(IntArray &answer, MaterialMode mmode)
Gives the inverted version of giveVoigtVectorMask.
void zero()
Zeroes all coefficient of receiver.
Definition: floatmatrix.C:1326
void giveStressBackVector(FloatArray &answer, const FloatArray &stressSpaceHardeningVars)
Domain * giveDomain() const
Definition: femcmpnn.h:100
double computeJ2InvariantAt(const FloatArray &stressVector)
REGISTER_Material(DummyMaterial)
#define IR_GIVE_OPTIONAL_FIELD(__ir, __value, __id)
Macro facilitating the use of input record reading methods.
Definition: inputrecord.h:78
int giveSize() const
Definition: intarray.h:203
int giveSize() const
Returns the size of receiver.
Definition: floatarray.h:218
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
#define IR_GIVE_FIELD(__ir, __value, __id)
Macro facilitating the use of input record reading methods.
Definition: inputrecord.h:69
void negated()
Switches the sign of every coefficient of receiver.
Definition: floatarray.C:739
Class representing integration point in finite element program.
Definition: gausspoint.h:93
Class representing solution step.
Definition: timestep.h:80
void add(const FloatArray &src)
Adds array src to receiver.
Definition: floatarray.C:156
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:29 for OOFEM by doxygen 1.8.11 written by Dimitri van Heesch, © 1997-2011