OOFEM  2.4
OOFEM.org - Object Oriented Finite Element Solver
tutorialmaterial.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 "tutorialmaterial.h"
36 #include "gausspoint.h"
37 #include "floatmatrix.h"
38 #include "floatarray.h"
39 #include "dynamicinputrecord.h"
40 #include "classfactory.h"
42 #include "mathfem.h"
43 
44 namespace oofem {
45 REGISTER_Material(TutorialMaterial);
46 
48 {}
49 
51 {}
52 
55 {
56  IRResultType result; // Required by IR_GIVE_FIELD macro
57 
58  result = D.initializeFrom(ir);
59  if ( result != IRRT_OK ) return result;
60 
62 
64 
66 }
67 
68 
69 void
71 {
73  D.giveInputRecord(ir);
74 
77 }
78 
79 
82 {
83  return new TutorialMaterialStatus(1, this->giveDomain(), gp);
84 }
85 
86 
87 void
89  const FloatArray &totalStrain, TimeStep *tStep)
90 {
91  FloatArray strain;
92  TutorialMaterialStatus *status = static_cast< TutorialMaterialStatus * >( this->giveStatus(gp) );
93 
94  // subtract stress thermal expansion
95  this->giveStressDependentPartOfStrainVector_3d(strain, gp, totalStrain, tStep, VM_Total);
96 
97  FloatArray trialElastStrain;
98  trialElastStrain.beDifferenceOf(strain, status->givePlasticStrain());
99 
100  // Compute trial stress = sig_tr = sig_old + E*delta_eps
101  FloatMatrix elasticStiffness;
102  D.give3dMaterialStiffnessMatrix(elasticStiffness, TangentStiffness, gp, tStep);
103  FloatArray trialStress;
104  trialStress.beProductOf(elasticStiffness, trialElastStrain);
105 
106  FloatArray sphTrialStress, devTrialStress;
107  computeSphDevPartOf(trialStress, sphTrialStress, devTrialStress);
108 
109  double J2 = this->computeSecondStressInvariant(devTrialStress);
110  double effectiveTrialStress = sqrt(3 * J2);
111 
112 #if 1
113  double temperatureScaling = 1.0;
114 #else
115  double temperature;
116  FloatArray et;
117  static_cast< StructuralElement *>(gp->giveIntegrationRule()->giveElement())->computeResultingIPTemperatureAt(et, tStep, gp, VM_Total);
118  temperature = et.at(1) + 800;
119 
120  double temperatureScaling = temperature <= 400 ? 1.0 : 1.0 - (temperature - 400) * 0.5 / 400;
121 #endif
122 
123  // evaluate the yield surface
124  double k = status->giveK();
125  double phiTrial = effectiveTrialStress - ( temperatureScaling * this->sig0 + temperatureScaling * H * k );
126 
127  if ( phiTrial < 0.0 ) { // elastic
128  answer = trialStress;
129 
130  status->letTempPlasticStrainBe(status->givePlasticStrain());
131  } else { // plastic loading
132  double G = D.giveShearModulus();
133  double mu = phiTrial / ( 3.0 * G + temperatureScaling * H ); // plastic multiplier
134  answer = devTrialStress * ( 1.0 - 3.0*G*mu/effectiveTrialStress); // radial return
135  answer.add(sphTrialStress);
136  k += mu;
137 
138  FloatArray plasticStrain = status->givePlasticStrain();
139  FloatArray dPlStrain;
140  applyDeviatoricElasticCompliance(dPlStrain, devTrialStress, 0.5);
141  plasticStrain.add(mu * 3. / (2. * effectiveTrialStress), dPlStrain);
142  status->letTempPlasticStrainBe(plasticStrain);
143  }
144 
145  // Store the temporary values for the given iteration
146  status->letTempStrainVectorBe(totalStrain);
147  status->letTempStressVectorBe(answer);
148  status->letTempKBe(k);
149  status->letTempDevTrialStressBe(devTrialStress);
150 }
151 
152 
153 void
155 {
156  double volumetricPart = ( sigV.at(1) + sigV.at(2) + sigV.at(3) ) / 3.0;
157  sigSph = {volumetricPart, volumetricPart, volumetricPart, 0.0, 0.0, 0.0};
158  sigDev.beDifferenceOf(sigV, sigSph);
159 }
160 
161 
162 void
164 {
165  answer.resize(6,6);
166  answer.zero();
167  answer.at(1,1) = answer.at(2,2) = answer.at(3,3) = 2.0/3.0;
168  answer.at(1,2) = answer.at(1,3) = answer.at(2,3) = -1.0/3.0;
169  answer.at(2,1) = answer.at(3,1) = answer.at(3,2) = -1.0/3.0;
170  answer.at(4,4) = answer.at(5,5) = answer.at(6,6) = 1.0/2.0;
171 }
172 
173 
174 void
176 {
177  TutorialMaterialStatus *status = static_cast< TutorialMaterialStatus * >( this->giveStatus(gp) );
178  FloatArray devTrialStress = status->giveTempDevTrialStress();
179 
180  double J2 = this->computeSecondStressInvariant(devTrialStress);
181  double effectiveTrialStress = sqrt(3 * J2);
182 
183 #if 1
184  double temperatureScaling = 1.0;
185 #else
186  double temperature;
187  FloatArray et;
188  static_cast< StructuralElement *>(gp->giveIntegrationRule()->giveElement())->computeResultingIPTemperatureAt(et, tStep, gp, VM_Total);
189  temperature = et.at(1) + 800;
190 
191  double temperatureScaling = temperature <= 400 ? 1.0 : 1.0 - (temperature - 400) * 0.5 / 400;
192 #endif
193 
194  // evaluate the yield surface
195  double k = status->giveK();
196  double phiTrial = effectiveTrialStress - ( temperatureScaling * this->sig0 + temperatureScaling * H * k );
197 
198  FloatMatrix elasticStiffness;
199  D.give3dMaterialStiffnessMatrix(elasticStiffness, ElasticStiffness, gp, tStep);
200 
201  if ( phiTrial < 0.0 ) { // elastic
202  answer = elasticStiffness;
203  } else { // plastic loading
204  double G = D.giveShearModulus();
205  // E_t = elasticStiffness - correction
206  // correction = 2.0 * G * ( 2.0 * G / h *( sig0 + kappa ) / sigtre *openProd(nu,nu) + mu*3*G/sigtre *Idev);
207  double h = 3.0 * G + temperatureScaling * H;
208  double mu = phiTrial / h; // plasic multiplier
209 
210  FloatArray nu = (3.0/2.0 / effectiveTrialStress ) * devTrialStress;
211  FloatMatrix Idev, correction;
213 
214  correction.plusDyadUnsym(nu, nu, 2.0 * G / h * ( temperatureScaling * this->sig0 + temperatureScaling * H * k ));
215  correction.add(mu * 3.0 * G, Idev);
216  correction.times(2.0 * G / effectiveTrialStress);
217  answer = elasticStiffness;
218  answer.subtract(correction);
219  }
220 }
221 
222 
223 void
225 {
226  double alpha = D.give(tAlpha, gp);
227  answer = {alpha, alpha, alpha, 0., 0., 0.};
228 }
229 
230 
231 int
233 {
234  TutorialMaterialStatus *status = static_cast< TutorialMaterialStatus * >( this->giveStatus(gp) );
235  if ( type == IST_PlasticStrainTensor ) {
236  answer = status->givePlasticStrain();
237  return 1;
238  } else {
239  return StructuralMaterial :: giveIPValue(answer, gp, type, tStep);
240  }
241 }
242 
243 
244 //=============================================================================
245 
246 
248  StructuralMaterialStatus(n, d, g),
249  tempPlasticStrain(6), plasticStrain(6),
250  tempDevTrialStress(6),
251  tempK(0.), k(0.)
252 {
253  strainVector.resize(6);
254  stressVector.resize(6);
257 }
258 
259 void
261 {
262  //StructuralMaterialStatus :: initTempStatus();
263 
264  // reset temp vars.
269 
270 
272  tempK = k;
275 }
276 
277 
278 void
280 {
281  // Copy the temporary values to the convered ones.
282  // This method is called after equilibrium has been reached and should
283  // set all...
285 
287  k = tempK;
288  // deviatoric trial stress is not really a state variable and was used not to repeat some code...
289 }
290 
291 } // end namespace oofem
FloatArray PVector
Equilibrated first Piola-Kirchhoff stress vector.
Definition: structuralms.h:78
const FloatArray & giveTempDevTrialStress()
virtual void updateYourself(TimeStep *tStep)
Update equilibrium history variables according to temp-variables.
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.
void setField(int item, InputFieldType id)
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
Class and object Domain.
Definition: domain.h:115
static void applyDeviatoricElasticCompliance(FloatArray &strain, const FloatArray &stress, double EModulus, double nu)
#define _IFT_TutorialMaterial_yieldstress
#define tAlpha
Definition: matconst.h:67
FloatArray tempStrainVector
Temporary strain vector in reduced form (to find balanced state)
Definition: structuralms.h:75
static void computeSphDevPartOf(const FloatArray &sigV, FloatArray &sigSph, FloatArray &sigDev)
virtual int giveIPValue(FloatArray &answer, GaussPoint *gp, InternalStateType type, TimeStep *tStep)
Returns the integration point corresponding value in Reduced form.
double & at(int i)
Coefficient access function.
Definition: floatarray.h:131
FloatArray tempFVector
Temporary deformation gradient in reduced form (to find balanced state)
Definition: structuralms.h:88
This class implements a structural material status information.
Definition: structuralms.h:65
FloatArray stressVector
Equilibrated stress vector in reduced form.
Definition: structuralms.h:71
Element * giveElement()
Returns reference to element containing receiver.
const FloatArray & givePlasticStrain()
double sig0
Initial (uniaxial) yield stress.
void plusDyadUnsym(const FloatArray &a, const FloatArray &b, double dV)
Adds to the receiver the product .
Definition: floatmatrix.C:811
#define _IFT_TutorialMaterial_hardeningmoduli
TutorialMaterial(int n, Domain *d)
double giveShearModulus()
Returns the shear elastic modulus .
MatResponseMode
Describes the character of characteristic material matrix.
void beDifferenceOf(const FloatArray &a, const FloatArray &b)
Sets receiver to be a - b.
Definition: floatarray.C:341
virtual IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description stored in input record.
Abstract base class for all "structural" finite elements.
FloatArray tempPlasticStrain
Temporary plastic strain (the given iteration)
virtual IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description stored in input record.
virtual void giveRealStressVector_3d(FloatArray &answer, GaussPoint *gp, const FloatArray &reducedE, TimeStep *tStep)
Default implementation relies on giveRealStressVector for second Piola-Kirchoff stress.
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 ...
FloatArray tempPVector
Temporary first Piola-Kirchhoff stress vector (to find balanced state)
Definition: structuralms.h:80
virtual void initTempStatus()
Initializes the temporary internal variables, describing the current state according to previously re...
FloatArray plasticStrain
Last equilibriated plastic strain (end of last time step)
void times(double f)
Multiplies receiver by factor f.
Definition: floatmatrix.C:1594
void letTempDevTrialStressBe(const FloatArray &values)
void giveStressDependentPartOfStrainVector_3d(FloatArray &answer, GaussPoint *gp, const FloatArray &reducedStrainVector, TimeStep *tStep, ValueModeType mode)
virtual void giveThermalDilatationVector(FloatArray &answer, GaussPoint *gp, TimeStep *tStep)
Returns a vector of coefficients of thermal dilatation in direction of each material principal (local...
void beProductOf(const FloatMatrix &aMatrix, const FloatArray &anArray)
Receiver becomes the result of the product of aMatrix and anArray.
Definition: floatarray.C:676
double at(int i, int j) const
Coefficient access function.
Definition: floatmatrix.h:176
void letTempPlasticStrainBe(const FloatArray &values)
static void giveDeviatoricProjectionMatrix(FloatMatrix &answer)
IntegrationRule * giveIntegrationRule()
Returns corresponding integration rule to receiver.
Definition: gausspoint.h:186
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
virtual void give3dMaterialStiffnessMatrix(FloatMatrix &answer, MatResponseMode mode, GaussPoint *gp, TimeStep *tStep)
Computes full 3d material stiffness matrix at given integration point, time, respecting load history ...
virtual IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description stored in input record.
FloatArray strainVector
Equilibrated strain vector in reduced form.
Definition: structuralms.h:69
Implementation of matrix containing floating point numbers.
Definition: floatmatrix.h:94
IRResultType
Type defining the return values of InputRecord reading operations.
Definition: irresulttype.h:47
void letTempStressVectorBe(const FloatArray &v)
Assigns tempStressVector to given vector v.
Definition: structuralms.h:135
virtual void giveInputRecord(DynamicInputRecord &input)
Setups the input record string of receiver.
void resize(int rows, int cols)
Checks size of receiver towards requested bounds.
Definition: floatmatrix.C:1358
static double computeSecondStressInvariant(const FloatArray &s)
Class representing the general Input Record.
Definition: inputrecord.h:101
FloatArray FVector
Equilibrated deformation gradient in reduced form.
Definition: structuralms.h:86
void add(const FloatMatrix &a)
Adds matrix to the receiver.
Definition: floatmatrix.C:1023
void zero()
Zeroes all coefficients of receiver.
Definition: floatarray.C:658
virtual void giveInputRecord(DynamicInputRecord &ir)
Setups the input record string of receiver.
Class representing the a dynamic Input Record.
IsotropicLinearElasticMaterial D
Abstract base class for all "structural" constitutive models.
void zero()
Zeroes all coefficient of receiver.
Definition: floatmatrix.C:1326
Domain * giveDomain() const
Definition: femcmpnn.h:100
virtual void giveInputRecord(DynamicInputRecord &input)
Setups the input record string of receiver.
virtual double give(int aProperty, GaussPoint *gp)
Returns the value of material property &#39;aProperty&#39;.
REGISTER_Material(DummyMaterial)
the oofem namespace is to define a context or scope in which all oofem names are defined.
virtual int giveIPValue(FloatArray &answer, GaussPoint *gp, InternalStateType type, TimeStep *tStep)
Returns the integration point corresponding value in Reduced form.
#define IR_GIVE_FIELD(__ir, __value, __id)
Macro facilitating the use of input record reading methods.
Definition: inputrecord.h:69
FloatArray tempStressVector
Temporary stress vector in reduced form (increments are used mainly in nonlinear analysis) ...
Definition: structuralms.h:73
Class representing integration point in finite element program.
Definition: gausspoint.h:93
double H
Hardening modulus.
Class representing solution step.
Definition: timestep.h:80
void add(const FloatArray &src)
Adds array src to receiver.
Definition: floatarray.C:156
TutorialMaterialStatus(int n, Domain *d, GaussPoint *g)
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:32 for OOFEM by doxygen 1.8.11 written by Dimitri van Heesch, © 1997-2011