OOFEM  2.4
OOFEM.org - Object Oriented Finite Element Solver
intmatisodamage.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 "intmatisodamage.h"
36 #include "gausspoint.h"
37 #include "floatmatrix.h"
38 #include "floatarray.h"
39 #include "mathfem.h"
40 #include "datastream.h"
41 #include "contextioerr.h"
42 #include "classfactory.h"
43 #include "dynamicinputrecord.h"
44 
45 namespace oofem {
46 REGISTER_Material(IntMatIsoDamage);
47 
49 {
50  maxOmega = 0.999999;
51  semiExplicit = false;
52 }
53 
54 
56 { }
57 
58 
59 void
61  const FloatArray &jump,
62  TimeStep *tStep)
63 {
64  IntMatIsoDamageStatus *status = static_cast< IntMatIsoDamageStatus * >( this->giveStatus(gp) );
65 
66  double f, equivJump, tempKappa = 0.0, omega = 0.0;
67 
68  // compute equivalent strain
69  this->computeEquivalentJump(equivJump, jump);
70 
71  // compute value of loading function
72  f = equivJump - status->giveKappa();
73 
74  if ( f <= 0.0 ) {
75  // no damage evolution
76  tempKappa = status->giveKappa();
77  omega = status->giveDamage();
78  } else {
79  // damage evolution
80  tempKappa = equivJump;
81  // evaluate damage parameter
82  this->computeDamageParam(omega, tempKappa);
83  }
84 
85  double strength = 1.0 - min(omega, maxOmega);
86 
87  if ( semiExplicit ) {
88  double oldOmega = status->giveDamage();
89  strength = 1.0 - min(oldOmega, maxOmega);
90  }
91 
92  //answer = {ks * jump.at(1) * strength, ks * jump.at(2) * strength, kn * jump.at(3)};
93  answer = {kn * jump.at(1), ks * jump.at(2) * strength, ks * jump.at(3) * strength};
94  // damage in tension only
95  if ( jump.at(1) >= 0 ) {
96  answer.at(1) *= strength;
97  }
98 
99  // update gp
100  status->letTempJumpBe(jump);
101  status->letTempTractionBe(answer);
102 
103  status->setTempKappa(tempKappa);
104  status->setTempDamage(omega);
105 }
106 
107 
108 void
110 {
111  IntMatIsoDamageStatus *status = static_cast< IntMatIsoDamageStatus * > ( this->giveStatus ( gp ) );
112  this->giveEngTraction_3d(answer, gp, jump, tStep);
113  status->letTempFirstPKTractionBe(answer);
114 
115 }
116 
117 
118 void
120  GaussPoint *gp, TimeStep *tStep)
121 {
122  IntMatIsoDamageStatus *status = static_cast< IntMatIsoDamageStatus * >( this->giveStatus(gp) );
123 
124  // assemble eleastic stiffness
125  answer.resize(2, 2);
126  answer.zero();
127  answer.at(1, 1) = kn;
128  answer.at(2, 2) = ks;
129 
130  if ( rMode == ElasticStiffness ) {
131  return;
132  }
133 
134  const FloatArray &jump3d = status->giveTempJump();
135  FloatArray jump2d = {jump3d.at(1),jump3d.at(2)};
136  double om = min(status->giveTempDamage(), maxOmega);
137  double un = jump2d.at(1);
138 
139  if ( rMode == SecantStiffness ) {
140  answer.at(2, 2) *= 1.0 - om;
141  // damage in tension only
142  if ( un >= 0 ) {
143  answer.at(1, 1) *= 1.0 - om;
144  }
145 
146  return;
147  } else { // Tangent Stiffness
148  answer.at(2, 2) *= 1.0 - om;
149  // damage in tension only
150  if ( un >= 0 ) {
151  answer.at(1, 1) *= 1.0 - om;
152 
153 #if 0
154  se.beProductOf(answer, jump2d);
155  double omega, omega_plus;
156  computeDamageParam(omega, un);
157  computeDamageParam(omega_plus, un + 1.0e-8 );
158  double dom = (omega_plus - omega)/ 1.0e-8;
159 
160  // d( (1-omega)*D*j ) / dj = (1-omega)D - D*j openprod domega/dj
161  double fac = ft*(e0 - un)/gf;
162  dom = e0*exp(fac) /(un*un + 1.0e-9) + e0*ft*exp(fac) / (gf*un + 1.0e-9);
163 
164 
165  dom = -( -e0 / (un * un+1.0e-9) * exp( -( ft / gf ) * ( un - e0 ) ) + e0 / (un+1.0e-9) * exp( -( ft / gf ) * ( un - e0 ) ) * ( -( ft / gf ) ) );
166  if ( ( om > 0. ) && ( status->giveTempKappa() > status->giveKappa() ) ) {
167  answer.at(1, 2) -= se.at(1) * dom;
168  answer.at(2, 2) -= se.at(2) * dom;
169  }
170 #endif
171  }
172  }
173 }
174 
175 
176 void
178  GaussPoint *gp, TimeStep *tStep)
179 {
180  IntMatIsoDamageStatus *status = static_cast< IntMatIsoDamageStatus * >( this->giveStatus(gp) );
181 
182  // assemble eleastic stiffness
183  answer.resize(3, 3);
184  answer.zero();
185  answer.at(1, 1) = kn;
186  answer.at(2, 2) = ks;
187  answer.at(3, 3) = ks;
188 
189  if ( rMode == ElasticStiffness ) {
190  return;
191  }
192 
193  const FloatArray &jump3d = status->giveTempJump();
194  double om = min(status->giveTempDamage(), maxOmega);
195  double un = jump3d.at(1);
196 
197  if ( rMode == SecantStiffness ) {
198  // Secant stiffness
199  answer.at(2, 2) *= 1.0 - om;
200  answer.at(3, 3) *= 1.0 - om;
201  // damage in tension only
202  if ( un >= 0 ) {
203  answer.at(1, 1) *= 1.0 - om;
204  }
205 
206  return;
207  } else {
208  // Tangent Stiffness
209  answer.at(2, 2) *= 1.0 - om;
210  answer.at(3, 3) *= 1.0 - om;
211  // damage in tension only
212  if ( un >= 0 ) {
213  answer.at(1, 1) *= 1.0 - om;
214  }
215  return;
216  }
217 }
218 
219 
220 int
222 {
223  IntMatIsoDamageStatus *status = static_cast< IntMatIsoDamageStatus * >( this->giveStatus(gp) );
224  if ( type == IST_DamageTensor ) {
225  answer.resize(6);
226  answer.zero();
227  answer.at(1) = answer.at(2) = answer.at(3) = status->giveDamage();
228  return 1;
229  } else if ( type == IST_DamageTensorTemp ) {
230  answer.resize(6);
231  answer.zero();
232  answer.at(1) = answer.at(2) = answer.at(3) = status->giveTempDamage();
233  return 1;
234  } else if ( type == IST_PrincipalDamageTensor ) {
235  answer.resize(3);
236  answer.at(1) = answer.at(2) = answer.at(3) = status->giveDamage();
237  return 1;
238  } else if ( type == IST_PrincipalDamageTempTensor ) {
239  answer.resize(3);
240  answer.at(1) = answer.at(2) = answer.at(3) = status->giveTempDamage();
241  return 1;
242  } else if ( type == IST_MaxEquivalentStrainLevel ) {
243  answer.resize(1);
244  answer.at(1) = status->giveKappa();
245  return 1;
246  } else {
247  return StructuralInterfaceMaterial :: giveIPValue(answer, gp, type, tStep);
248  }
249 }
250 
251 
254 {
255  IRResultType result; // Required by IR_GIVE_FIELD macro
256 
259 
262  this->e0 = ft / kn;
263 
264  //Set limit on the maximum isotropic damage parameter if needed
265  maxOmega = 0.999999;
267  maxOmega = min(maxOmega, 0.999999);
268  maxOmega = max(maxOmega, 0.0);
269 
271 }
272 
273 
274 void
276 {
278 
279  input.setField(this->kn, _IFT_IntMatIsoDamage_kn);
280  input.setField(this->ks, _IFT_IntMatIsoDamage_ks);
281 
282  input.setField(this->ft, _IFT_IntMatIsoDamage_ft);
283  input.setField(this->gf, _IFT_IntMatIsoDamage_gf);
284 
286 }
287 
288 
289 void
291 {
292  kappa = jump.at(1);
293 }
294 
295 void
296 IntMatIsoDamage :: computeDamageParam(double &omega, double kappa)
297 {
298  if ( kappa > this->e0 ) {
299  omega = 1.0 - ( this->e0 / kappa ) * exp( -( ft / gf ) * ( kappa - e0 ) );
300  } else {
301  omega = 0.0;
302  }
303 }
304 
305 
307 {
308  kappa = tempKappa = 0.0;
309  damage = tempDamage = 0.0;
310 }
311 
312 
314 { }
315 
316 
317 void
319 {
321  fprintf(file, "status { ");
322  if ( this->damage > 0.0 ) {
323  fprintf(file, "kappa %f, damage %f ", this->kappa, this->damage);
324  }
325 
326  fprintf(file, "}\n");
327 }
328 
329 
330 void
332 {
334  this->tempKappa = this->kappa;
335  this->tempDamage = this->damage;
336 }
337 
338 void
340 {
342  this->kappa = this->tempKappa;
343  this->damage = this->tempDamage;
344 }
345 
346 
349 {
350  contextIOResultType iores;
351 
352  // save parent class status
353  if ( ( iores = StructuralInterfaceMaterialStatus :: saveContext(stream, mode, obj) ) != CIO_OK ) {
354  THROW_CIOERR(iores);
355  }
356 
357  // write a raw data
358  if ( !stream.write(kappa) ) {
360  }
361 
362  if ( !stream.write(damage) ) {
364  }
365 
366  return CIO_OK;
367 }
368 
371 {
372  contextIOResultType iores;
373 
374  // read parent class status
375  if ( ( iores = StructuralInterfaceMaterialStatus :: restoreContext(stream, mode, obj) ) != CIO_OK ) {
376  THROW_CIOERR(iores);
377  }
378 
379  // read raw data
380  if ( !stream.read(kappa) ) {
382  }
383 
384  if ( !stream.read(damage) ) {
386  }
387 
388  return CIO_OK;
389 }
390 
391 } // end namespace oofem
InternalStateType
Type representing the physical meaning of element or constitutive model internal variable.
void setField(int item, InputFieldType id)
virtual IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description stored in input record.
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
#define _IFT_IntMatIsoDamage_maxOmega
double giveKappa()
Returns the last equilibrated scalar measure of the largest jump level.
The purpose of DataStream abstract class is to allow to store/restore context to different streams...
Definition: datastream.h:54
double & at(int i)
Coefficient access function.
Definition: floatarray.h:131
int max(int i, int j)
Returns bigger value form two given decimals.
Definition: mathfem.h:71
double gf
Fracture energy.
virtual contextIOResultType saveContext(DataStream &stream, ContextMode mode, void *obj=NULL)
Stores receiver state to output stream.
virtual ~IntMatIsoDamageStatus()
Destructor.
virtual void giveInputRecord(DynamicInputRecord &input)
Setups the input record string of receiver.
double kn
Elastic properties (normal moduli).
General IO error.
double tempKappa
Non-equilibrated scalar measure of the largest equivalent displacement.
double giveDamage()
Returns the last equilibrated damage level.
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.
virtual contextIOResultType restoreContext(DataStream &stream, ContextMode mode, void *obj=NULL)
Restores the receiver state previously written in stream.
#define THROW_CIOERR(e)
Definition: contextioerr.h:61
double tempDamage
Non-equilibrated damage level of material.
virtual void updateYourself(TimeStep *tStep)
Update equilibrium history variables according to temp-variables.
IntMatIsoDamage(int n, Domain *d)
Constructor.
virtual int giveIPValue(FloatArray &answer, GaussPoint *gp, InternalStateType type, TimeStep *tStep)
Returns the integration point corresponding value in Reduced form.
virtual void computeDamageParam(double &omega, double kappa)
computes the value of damage parameter omega, based on given value of equivalent strain.
virtual int write(const int *data, int count)=0
Writes count integer values from array pointed by data.
double giveTempDamage()
Returns the temp. damage level.
double e0
Limit elastic deformation.
virtual void updateYourself(TimeStep *tStep)
Update equilibrium history variables according to temp-variables.
virtual void computeEquivalentJump(double &kappa, const FloatArray &jump)
Computes the equivalent jump measure from given jump vector (full form).
virtual contextIOResultType restoreContext(DataStream &stream, ContextMode mode, void *obj=NULL)
Restores the receiver state previously written in stream.
void letTempJumpBe(FloatArray v)
Assigns tempJump to given vector v.
double ks
Shear moduli.
IntMatIsoDamageStatus(int n, Domain *d, GaussPoint *g)
Constructor.
double kappa
Scalar measure of the largest equivalent displacement ever reached in material.
This class implements a structural interface material status information.
#define _IFT_IntMatIsoDamage_ks
virtual void give2dStiffnessMatrix_Eng(FloatMatrix &answer, MatResponseMode rMode, GaussPoint *gp, TimeStep *tStep)
double at(int i, int j) const
Coefficient access function.
Definition: floatmatrix.h:176
virtual void give3dStiffnessMatrix_Eng(FloatMatrix &answer, MatResponseMode rMode, GaussPoint *gp, TimeStep *tStep)
virtual void giveEngTraction_3d(FloatArray &answer, GaussPoint *gp, const FloatArray &jump, TimeStep *tStep)
void letTempFirstPKTractionBe(FloatArray v)
Assigns tempFirstPKTraction to given vector v.
virtual void printOutputAt(FILE *file, TimeStep *tStep)
Print receiver&#39;s output to given stream.
Class representing vector of real numbers.
Definition: floatarray.h:82
virtual void giveFirstPKTraction_3d(FloatArray &answer, GaussPoint *gp, const FloatArray &jump, const FloatMatrix &F, TimeStep *tStep)
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 ~IntMatIsoDamage()
Destructor.
void setTempKappa(double newKappa)
Sets the temp scalar measure of the largest strain level to given value.
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
virtual IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description stored in input record.
virtual int giveIPValue(FloatArray &answer, GaussPoint *gp, InternalStateType type, TimeStep *tStep)
Returns the integration point corresponding value in Reduced form.
void zero()
Zeroes all coefficients of receiver.
Definition: floatarray.C:658
virtual contextIOResultType saveContext(DataStream &stream, ContextMode mode, void *obj=NULL)
Stores receiver state to output stream.
Class representing the a dynamic Input Record.
long ContextMode
Context mode (mask), defining the type of information written/read to/from context.
Definition: contextmode.h:43
double ft
Tension strength.
virtual void initTempStatus()
Initializes the temporary internal variables, describing the current state according to previously re...
void zero()
Zeroes all coefficient of receiver.
Definition: floatmatrix.C:1326
int min(int i, int j)
Returns smaller value from two given decimals.
Definition: mathfem.h:59
double giveTempKappa()
Returns the temp. scalar measure of the largest jump level.
void setTempDamage(double newDamage)
Sets the temp damage level to given value.
REGISTER_Material(DummyMaterial)
void letTempTractionBe(FloatArray v)
Assigns tempTraction to given vector v.
#define IR_GIVE_OPTIONAL_FIELD(__ir, __value, __id)
Macro facilitating the use of input record reading methods.
Definition: inputrecord.h:78
virtual void printOutputAt(FILE *file, TimeStep *tStep)
Print receiver&#39;s output to given stream.
Abstract base class for all "structural" interface models.
the oofem namespace is to define a context or scope in which all oofem names are defined.
#define _IFT_IntMatIsoDamage_gf
double damage
Damage level of material.
#define IR_GIVE_FIELD(__ir, __value, __id)
Macro facilitating the use of input record reading methods.
Definition: inputrecord.h:69
const FloatArray & giveTempJump() const
Returns the const pointer to receiver&#39;s temporary jump.
virtual void initTempStatus()
Initializes the temporary internal variables, describing the current state according to previously re...
#define _IFT_IntMatIsoDamage_ft
double maxOmega
Maximum limit on omega. The purpose is elimination of a too compliant material which may cause conver...
#define _IFT_IntMatIsoDamage_kn
Class representing integration point in finite element program.
Definition: gausspoint.h:93
virtual void giveInputRecord(DynamicInputRecord &input)
Setups the input record string of receiver.
Class representing solution step.
Definition: timestep.h:80
This class implements the InterfaceMaterialStatus associated with IntMatIsoDamage.
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