OOFEM  2.4
OOFEM.org - Object Oriented Finite Element Solver
maxwellChM.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 "mathfem.h"
36 #include "maxwellChM.h"
37 #include "floatarray.h"
38 #include "floatmatrix.h"
39 #include "gausspoint.h"
40 #include "timestep.h"
41 #include "contextioerr.h"
42 
43 namespace oofem {
45 { }
46 
47 
48 void
50 {
51  int rSize;
52  FloatArray rhs(this->nUnits), discreteRelaxFunctionVal;
53  FloatMatrix A(this->nUnits, this->nUnits);
54 
55  const FloatArray &rTimes = this->giveDiscreteTimes();
56  rSize = rTimes.giveSize();
57 
58  // compute discrete values of the relaxation function at times rTimes
59  // from the creep function (by numerically solving integral equations)
60  //
61  // a direct call to the relaxation function could be used, if available
62  this->computeDiscreteRelaxationFunction(discreteRelaxFunctionVal,
63  rTimes,
64  tPrime,
65  tPrime,
66  gp,
67  tStep);
68 
69  // assemble the matrix of the set of linear equations
70  // for computing the optimal moduli
71  for ( int i = 1; i <= this->nUnits; i++ ) {
72  double taui = this->giveCharTime(i);
73  for ( int j = 1; j <= this->nUnits; j++ ) {
74  double tauj = this->giveCharTime(j);
75  double sum = 0;
76  for ( int r = 1; r <= rSize; r++ ) {
77  double tti = pow( ( tPrime + rTimes.at(r) ) / taui, giveCharTimeExponent(i) ) -
78  pow( tPrime / taui, giveCharTimeExponent(i) );
79  double ttj = pow( ( tPrime + rTimes.at(r) ) / tauj, giveCharTimeExponent(j) ) -
80  pow( tPrime / tauj, giveCharTimeExponent(j) );
81  sum += exp(-tti - ttj);
82  }
83 
84  A.at(i, j) = sum;
85  }
86 
87  // assemble rhs
88  double sumRhs = 0;
89  for ( int r = 1; r <= rSize; r++ ) {
90  double tti = pow( ( tPrime + rTimes.at(r) ) / taui, giveCharTimeExponent(i) ) -
91  pow( tPrime / taui, giveCharTimeExponent(i) );
92  sumRhs += exp(-tti) * discreteRelaxFunctionVal.at(r);
93  }
94 
95  rhs.at(i) = sumRhs;
96  }
97 
98  // solve the linear system
99  A.solveForRhs(rhs, answer);
100 }
101 
102 
103 
104 double
106 {
107  /*
108  * This function returns the incremental modulus for the given time increment.
109  * The modulus may also depend on the specimen geometry (gp - dependence).
110  *
111  * Note: time -1 refers to the previous time.
112  */
113  double E = 0.0;
114 
116 
117  if ( (tStep->giveIntrinsicTime() < this->castingTime) ) {
118  OOFEM_ERROR("Attempted to evaluate E modulus at time lower than casting time");
119  }
120 
121  double tPrime = relMatAge + ( tStep->giveTargetTime() - 0.5 * tStep->giveTimeIncrement() ) / timeFactor;
122  this->updateEparModuli(tPrime, gp, tStep);
123 
124  for ( int mu = 1; mu <= nUnits; mu++ ) {
125  double deltaYmu = tStep->giveTimeIncrement() / timeFactor / this->giveCharTime(mu);
126  if ( deltaYmu <= 0.0 ) {
127  deltaYmu = 1.e-3;
128  }
129 
130  deltaYmu = pow( deltaYmu, this->giveCharTimeExponent(mu) );
131 
132  double lambdaMu = ( 1.0 - exp(-deltaYmu) ) / deltaYmu;
133  double Emu = this->giveEparModulus(mu); // previously updated by updateEparModuli
134  E += lambdaMu * Emu;
135  }
136 
137  return E;
138 }
139 
140 
141 
142 
143 void
145  GaussPoint *gp, TimeStep *tStep, ValueModeType mode)
146 //
147 // computes the strain due to creep at constant stress during the increment
148 // (in fact, the INCREMENT of creep strain is computed for mode == VM_Incremental)
149 //
150 {
151  FloatArray help, reducedAnswer;
152  FloatArray sigmaMu;
153  FloatMatrix B;
154  MaxwellChainMaterialStatus *status = static_cast< MaxwellChainMaterialStatus * >( this->giveStatus(gp) );
155 
156  if ( (tStep->giveIntrinsicTime() < this->castingTime) ) {
157  OOFEM_ERROR("Attempted to evaluate creep strain for time lower than casting time");
158  }
159 
160  if ( mode == VM_Incremental ) {
161  this->giveUnitComplianceMatrix(B, gp, tStep);
162  reducedAnswer.resize( B.giveNumberOfRows() );
163  reducedAnswer.zero();
164 
165  for ( int mu = 1; mu <= nUnits; mu++ ) {
166  double deltaYmu = tStep->giveTimeIncrement() / timeFactor / this->giveCharTime(mu);
167  deltaYmu = pow( deltaYmu, this->giveCharTimeExponent(mu) );
168 
169  sigmaMu = status->giveHiddenVarsVector(mu); // JB
170 
171  if ( sigmaMu.giveSize() ) {
172  help.beProductOf(B, sigmaMu); // B can be moved before sum !!!
173  help.times( 1.0 - exp(-deltaYmu) );
174  reducedAnswer.add(help);
175  }
176  }
177 
178  double E = this->giveEModulus(gp, tStep);
179  // E = this->giveIncrementalModulus(gp, tStep);
180  reducedAnswer.times(1.0 / E);
181 
182  answer = reducedAnswer;
183  } else {
184  /* error - total mode not implemented yet */
185  OOFEM_ERROR("mode is not supported");
186  }
187 }
188 
189 
190 void
192 {
193  RheoChainMaterial :: giveRealStressVector(answer, gp, reducedStrain, tStep);
194 
195  // Computes hidden variables and stores them as temporary
196  this->computeHiddenVars(gp, tStep);
197 }
198 
199 
200 void
202 {
203  /*
204  * Updates hidden variables used to effectively trace the load history
205  */
206  FloatArray help, deltaEps0, help1;
207  FloatArray muthHiddenVarsVector;
208 
209  FloatMatrix Binv;
211  static_cast< MaxwellChainMaterialStatus * >( this->giveStatus(gp) );
212 
213  this->giveUnitStiffnessMatrix(Binv, gp, tStep);
214  help = status->giveTempStrainVector();
215  help.subtract( status->giveStrainVector() );
216 
217  // Subtract the stress-independent part of strain
218  this->computeTrueStressIndependentStrainVector(deltaEps0, gp, tStep, VM_Incremental);
219  if ( deltaEps0.giveSize() ) {
220  help.subtract(deltaEps0);
221  }
222 
223  help1.beProductOf(Binv, help);
224 
226  double tPrime = relMatAge + ( tStep->giveTargetTime() - 0.5 * tStep->giveTimeIncrement() ) / timeFactor;
227  this->updateEparModuli(tPrime, gp, tStep);
228 
229  for ( int mu = 1; mu <= nUnits; mu++ ) {
230  double deltaYmu = tStep->giveTimeIncrement() / timeFactor / this->giveCharTime(mu);
231  deltaYmu = pow( deltaYmu, this->giveCharTimeExponent(mu) );
232 
233  double lambdaMu = ( 1.0 - exp(-deltaYmu) ) / deltaYmu;
234  double Emu = this->giveEparModulus(mu);
235 
236  muthHiddenVarsVector = status->giveHiddenVarsVector(mu);
237  help = help1;
238  help.times(lambdaMu * Emu);
239  if ( muthHiddenVarsVector.giveSize() ) {
240  muthHiddenVarsVector.times( exp(-deltaYmu) );
241  muthHiddenVarsVector.add(help);
242  status->letTempHiddenVarsVectorBe(mu, muthHiddenVarsVector);
243  } else {
244  status->letTempHiddenVarsVectorBe(mu, help);
245  }
246  }
247 }
248 
249 
252 {
253  return new MaxwellChainMaterialStatus(1, this->giveDomain(), gp, nUnits);
254 }
255 
256 
259 {
261 }
262 
263 /****************************************************************************************/
264 
266  GaussPoint *g, int nunits) :
267  RheoChainMaterialStatus(n, d, g, nunits) { }
268 
269 
270 void
272 {
274 }
275 
276 
277 void
279 {
281 }
282 
283 void
285 {
287 }
288 
291 {
292  contextIOResultType iores;
293 
294  if ( ( iores = RheoChainMaterialStatus :: saveContext(stream, mode, obj) ) != CIO_OK ) {
295  THROW_CIOERR(iores);
296  }
297 
298  return CIO_OK;
299 }
300 
301 
304 {
306 
307  if ( iores != CIO_OK ) {
308  THROW_CIOERR(iores);
309  }
310 
311  return CIO_OK;
312 }
313 } // end namespace oofem
double relMatAge
Physical age of the material at simulation time = 0.
Definition: rheoChM.h:138
void subtract(const FloatArray &src)
Subtracts array src to receiver.
Definition: floatarray.C:258
virtual MaterialStatus * CreateStatus(GaussPoint *gp) const
Creates new copy of associated status and inserts it into given integration point.
Definition: maxwellChM.C:251
virtual void updateYourself(TimeStep *tStep)
Update equilibrium history variables according to temp-variables.
Definition: maxwellChM.C:278
virtual MaterialStatus * giveStatus(GaussPoint *gp) const
Returns material status of receiver in given integration point.
Definition: material.C:244
double giveCharTime(int) const
Access to the characteristic time of a given unit.
Definition: rheoChM.C:397
Class and object Domain.
Definition: domain.h:115
virtual contextIOResultType restoreContext(DataStream &stream, ContextMode mode, void *obj=NULL)
Restores the receiver state previously written in stream.
Definition: maxwellChM.C:303
virtual void printOutputAt(FILE *file, TimeStep *tStep)
Print receiver&#39;s output to given stream.
Definition: maxwellChM.C:271
bool solveForRhs(const FloatArray &b, FloatArray &answer, bool transpose=false)
Solves the system of linear equations .
Definition: floatmatrix.C:1112
virtual void updateYourself(TimeStep *tStep)
Update equilibrium history variables according to temp-variables.
Definition: rheoChM.C:781
virtual IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description stored in input record.
Definition: maxwellChM.C:258
virtual void printOutputAt(FILE *file, TimeStep *tStep)
Print receiver&#39;s output to given stream.
Definition: rheoChM.C:751
The purpose of DataStream abstract class is to allow to store/restore context to different streams...
Definition: datastream.h:54
void computeTrueStressIndependentStrainVector(FloatArray &answer, GaussPoint *gp, TimeStep *tStep, ValueModeType mode)
Computes, for the given integration point, the strain vector induced by stress-independent internal p...
Definition: rheoChM.C:344
double & at(int i)
Coefficient access function.
Definition: floatarray.h:131
ValueModeType
Type representing the mode of UnknownType or CharType, or similar types.
Definition: valuemodetype.h:78
virtual void giveRealStressVector(FloatArray &answer, GaussPoint *gp, const FloatArray &reducedStrain, TimeStep *tStep)
Computes the real stress vector for given total strain and integration point.
Definition: rheoChM.C:81
virtual contextIOResultType saveContext(DataStream &stream, ContextMode mode, void *obj=NULL)
Stores receiver state to output stream.
Definition: rheoChM.C:802
This class implements associated Material Status to RheoChainMaterial.
Definition: rheoChM.h:71
virtual contextIOResultType restoreContext(DataStream &stream, ContextMode mode, void *obj=NULL)
Restores the receiver state previously written in stream.
Definition: rheoChM.C:829
double giveTargetTime()
Returns target time.
Definition: timestep.h:146
This class implements a rheologic chain model describing a viscoelastic material. ...
Definition: rheoChM.h:130
virtual void giveEigenStrainVector(FloatArray &answer, GaussPoint *gp, TimeStep *tStep, ValueModeType mode)
Computes, for the given integration point, the strain vector induced by the stress history (typically...
Definition: maxwellChM.C:144
MaxwellChainMaterialStatus(int n, Domain *d, GaussPoint *g, int nunits)
Definition: maxwellChM.C:265
virtual void computeCharCoefficients(FloatArray &answer, double tPrime, GaussPoint *gp, TimeStep *tStep)
This function computes the moduli of individual Maxwell units such that the corresponding Dirichlet s...
Definition: maxwellChM.C:49
const FloatArray & giveDiscreteTimes()
Definition: rheoChM.C:260
#define THROW_CIOERR(e)
Definition: contextioerr.h:61
double giveEparModulus(int iChain)
Access to partial modulus of a given unit.
Definition: rheoChM.C:303
void letTempHiddenVarsVectorBe(int i, FloatArray &valueArray)
Definition: rheoChM.C:739
double giveTimeIncrement()
Returns solution step associated time increment.
Definition: timestep.h:150
void giveUnitStiffnessMatrix(FloatMatrix &answer, GaussPoint *gp, TimeStep *tStep)
Evaluation of elastic stiffness matrix for unit Young&#39;s modulus.
Definition: rheoChM.C:267
double timeFactor
Scaling factor transforming the simulation time units into days (gives the number of simulation time ...
Definition: rheoChM.h:167
#define E(p)
Definition: mdm.C:368
void giveUnitComplianceMatrix(FloatMatrix &answer, GaussPoint *gp, TimeStep *tStep)
Evaluation of elastic compliance matrix for unit Young&#39;s modulus.
Definition: rheoChM.C:286
virtual void giveRealStressVector(FloatArray &answer, GaussPoint *gp, const FloatArray &reducedStrain, TimeStep *tStep)
Computes the real stress vector for given total strain and integration point.
Definition: maxwellChM.C:191
virtual double giveEModulus(GaussPoint *gp, TimeStep *tStep)
Evaluation of the incremental modulus.
Definition: maxwellChM.C:105
#define OOFEM_ERROR(...)
Definition: error.h:61
void beProductOf(const FloatMatrix &aMatrix, const FloatArray &anArray)
Receiver becomes the result of the product of aMatrix and anArray.
Definition: floatarray.C:676
void computeDiscreteRelaxationFunction(FloatArray &answer, const FloatArray &tSteps, double t0, double tr, GaussPoint *gp, TimeStep *tSte)
Evaluation of the relaxation function at given times.
Definition: rheoChM.C:193
const FloatArray & giveTempStrainVector() const
Returns the const pointer to receiver&#39;s temporary strain vector.
Definition: structuralms.h:115
double at(int i, int j) const
Coefficient access function.
Definition: floatmatrix.h:176
double giveIntrinsicTime()
Returns intrinsic time, e.g. time in which constitutive model is evaluated.
Definition: timestep.h:148
FloatArray & giveHiddenVarsVector(int i)
Definition: rheoChM.h:99
Abstract base class representing a material status information.
Definition: matstatus.h:84
virtual void updateEparModuli(double tPrime, GaussPoint *gp, TimeStep *tStep)
Update of partial moduli of individual chain units.
Definition: rheoChM.C:313
virtual double giveCharTimeExponent(int i)
Exponent to be used with the char time of a given unit, usually = 1.0.
Definition: rheoChM.h:356
MaxwellChainMaterial(int n, Domain *d)
Definition: maxwellChM.C:44
Class representing vector of real numbers.
Definition: floatarray.h:82
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
Class representing the general Input Record.
Definition: inputrecord.h:101
void zero()
Zeroes all coefficients of receiver.
Definition: floatarray.C:658
virtual void initTempStatus()
Initializes the temporary internal variables, describing the current state according to previously re...
Definition: maxwellChM.C:284
virtual void initTempStatus()
Initializes the temporary internal variables, describing the current state according to previously re...
Definition: rheoChM.C:796
void times(double s)
Multiplies receiver with scalar.
Definition: floatarray.C:818
long ContextMode
Context mode (mask), defining the type of information written/read to/from context.
Definition: contextmode.h:43
virtual IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description stored in input record.
Definition: rheoChM.C:583
void computeHiddenVars(GaussPoint *gp, TimeStep *tStep)
Definition: maxwellChM.C:201
Domain * giveDomain() const
Definition: femcmpnn.h:100
double castingTime
Casting time.
Definition: material.h:113
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.
const FloatArray & giveStrainVector() const
Returns the const pointer to receiver&#39;s strain vector.
Definition: structuralms.h:105
int giveNumberOfRows() const
Returns number of rows of receiver.
Definition: floatmatrix.h:156
int nUnits
Number of (Maxwell or Kelvin) units in the rheologic chain.
Definition: rheoChM.h:136
This class implements associated Material Status to MaxwellChainMaterial.
Definition: maxwellChM.h:44
virtual contextIOResultType saveContext(DataStream &stream, ContextMode mode, void *obj=NULL)
Stores receiver state to output stream.
Definition: maxwellChM.C:290
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