OOFEM  2.4
OOFEM.org - Object Oriented Finite Element Solver
latticedamage2d.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 "latticedamage2d.h"
36 #include "gausspoint.h"
37 #include "floatmatrix.h"
38 #include "floatarray.h"
39 #include "datastream.h"
40 #include "contextioerr.h"
41 #include "mathfem.h"
44 #include "staggeredproblem.h"
45 #include "classfactory.h"
46 #ifdef __TM_MODULE
47  #include "../tm/latticetransportelement.h"
48 #endif
49 
50 #include <cstdlib>
51 
52 namespace oofem {
53 REGISTER_Material(LatticeDamage2d);
54 
56 { }
57 
58 
60 { }
61 
62 int
64 {
65  return mode == _2dLattice;
66 }
67 
68 
71 {
72  IRResultType result; // Required by IR_GIVE_FIELD macro
73 
75  if ( result != IRRT_OK ) return result;
76 
78  if ( result != IRRT_OK ) return result;
79 
80  double value = 0.;
83 
85 
86  cAlpha = 0;
88 
89 
90  //factor which relates the shear stiffness to the normal stiffness. Default is 1
91  alphaOne = 1.;
94 
95  //Parameter which is used for the definition of the moment.
96  alphaTwo = 1.;
98 
100 
101  localRandomType = 0; //Default: No local random field
103  if ( localRandomType == 1 ) { //Gaussian random generator
106  }
107 
108  e0Mean = 0;
110 
113 
114  softeningType = 0;
116 
117  if ( softeningType == 1 || softeningType == 3 ) { //linear or exponential softening
119  } else if ( softeningType == 2 ) { //bilinear softening
121  wfOne = 0.15 * wf;
123  e0OneMean = 0.3 * e0Mean;
125  } else {
126  OOFEM_WARNING("Unknown softening type");
127  return IRRT_BAD_FORMAT;
128  }
129 
130  this->biotCoefficient = 0.;
132 
133  this->biotType = 0;
135 
136  return IRRT_OK;
137 }
138 
139 void
140 LatticeDamage2d :: computeEquivalentStrain(double &tempEquivStrain, const FloatArray &strain, GaussPoint *gp, TimeStep *tStep)
141 {
142  if ( strain.isEmpty() ) {
143  tempEquivStrain = 0.;
144  return;
145  }
146 
147  //LatticeDamage2dStatus *status = static_cast< LatticeDamage2dStatus * >( this->giveStatus(gp) );
148  const double e0 = this->give(e0_ID, gp) * this->e0Mean;
149 
150  double paramA, paramB, paramC;
151 
152  paramA = 0.5 * ( e0 + ec * e0 );
153  paramB = ( coh * e0 ) / sqrt( 1. - pow( ( ec * e0 - e0 ) / ( e0 + ec * e0 ), 2. ) );
154  paramC = 0.5 * ( this->ec * e0 - e0 );
155 
156  //double equivStrain = status->giveEquivalentStrain();
157  tempEquivStrain = sqrt( pow(this->alphaOne * strain.at(2) / paramB, 2.) + pow( ( strain.at(1) + paramC ) / paramA, 2. ) ) * paramA - paramC;
158 }
159 
160 void
161 LatticeDamage2d :: computeDamageParam(double &omega, double tempKappa, const FloatArray &strain, GaussPoint *gp)
162 {
163  double le = ( static_cast< LatticeStructuralElement * >( gp->giveElement() ) )->giveLength();
164  const double e0 = this->give(e0_ID, gp) * this->e0Mean;
165  omega = 0.0;
166 
167  int nite = 0;
168  double R, Lhs, Ft, help;
169 
170  if ( softeningType == 1 ) { //linear
171  if ( tempKappa >= e0 && tempKappa < this->wf / le ) {
172  //linear stress-crack opening relation
173  //check if input parameter make sense
174  if ( this->wf / le <= e0 ) {
175  OOFEM_ERROR("e0>wf/Le \n Possible solutions: Increase fracture energy or reduce element size");
176  }
177 
178  omega = ( 1. - e0 / tempKappa ) / ( 1. - e0 / ( this->wf / le ) );
179  } else if ( tempKappa >= this->wf / le ) {
180  omega = 1.;
181  } else {
182  omega = 0.;
183  }
184  } else if ( softeningType == 2 ) { //bilinear softening
185  double helpStrain = 0.3 * e0;
186 
187  //Check if input parameter make sense
188  if ( e0 > wfOne / le ) {
189  OOFEM_ERROR("parameter wf1 is too small");
190  } else if ( wfOne / le > this->wf / le ) {
191  OOFEM_ERROR("parameter wf is too small");
192  }
193 
194  //
195  if ( tempKappa > e0 ) {
196  omega = ( 1 - e0 / tempKappa ) / ( ( helpStrain - e0 ) / this->wfOne * le + 1. );
197 
198  if ( omega * tempKappa * le > 0 && omega * tempKappa * le < this->wfOne ) {
199  return;
200  } else {
201  omega = ( 1. - helpStrain / tempKappa - helpStrain * this->wfOne / ( tempKappa * ( this->wf - this->wfOne ) ) ) / ( 1. - helpStrain * le / ( this->wf - this->wfOne ) );
202 
203  if ( omega * tempKappa * le >= this->wfOne && omega * tempKappa * le < this->wf ) {
204  return;
205  }
206  }
207 
208  if ( omega > 1. ) {
209  omega = 1.;
210  } else if ( omega < 0 ) {
211  omega = 0;
212  }
213  } else {
214  omega = 0.;
215  }
216  } else if ( softeningType == 3 ) { //exponential softening
217  // iteration to achieve objectivity
218  // we are finding state, where elastic stress is equal to
219  // stress from crack-opening relation (wf = wf characterizes the carc opening diagram)
220 
221  if ( tempKappa <= e0 ) {
222  omega = 0.0;
223  } else {
224  omega = 0.0;
225  Ft = this->eNormal * e0;
226  do {
227  nite++;
228  help = le * omega * tempKappa / this->wf;
229  R = ( 1. - omega ) * eNormal * tempKappa - Ft *exp(-help);
230  Lhs = eNormal * tempKappa - Ft *exp(-help) * le * tempKappa / this->wf;
231  omega += R / Lhs;
232  if ( nite > 40 ) {
233  OOFEM_ERROR("algorithm not converging");
234  }
235  } while ( fabs(R) >= 1.e-4 );
236 
237  if ( ( omega > 1.0 ) || ( omega < 0.0 ) ) {
238  OOFEM_ERROR("internal error");
239  }
240  }
241  } else {
242  OOFEM_ERROR("Unknown softening type");
243  }
244 }
245 
246 
247 void
249  GaussPoint *gp,
250  TimeStep *tStep,
251  ValueModeType mode)
252 {
253  FloatArray et;
254  LatticeStructuralElement *lselem = static_cast< LatticeStructuralElement * >( gp->giveElement() );
255 
256  answer.resize(3);
257  answer.zero();
258 
259  if ( tStep->giveIntrinsicTime() < this->castingTime ) {
260  return;
261  }
262 
263  if ( lselem ) {
264  lselem->computeResultingIPTemperatureAt(et, tStep, gp, mode);
265  }
266 
267  if ( et.giveSize() == 0 ) {
268  answer.clear();
269  return;
270  }
271 
272  if ( et.giveSize() < 1 ) {
273  OOFEM_ERROR("Bad format of TemperatureLoad");
274  exit(1);
275  }
276 
277  if ( mode == VM_Total ) {
278  // compute temperature difference
279  double deltaTemperature = et.at(1) - this->referenceTemperature;
280  answer.at(1) = this->give(tAlpha, gp) * deltaTemperature;
281  } else {
282  answer.at(1) = this->give(tAlpha, gp) * et.at(1);
283  }
284 
285  double length = ( static_cast< LatticeStructuralElement * >( gp->giveElement() ) )->giveLength();
286 
287  answer.at(1) += this->cAlpha * et.at(1) / length;
288  return;
289 }
290 
293 {
295 }
296 
299 {
300  MaterialStatus *status = static_cast< MaterialStatus * >( gp->giveMaterialStatus() );
301  if ( status == NULL ) {
302  // create a new one
303  status = this->CreateStatus(gp);
304 
305  if ( status != NULL ) {
306  gp->setMaterialStatus( status, this->giveNumber() );
307  this->_generateStatusVariables(gp);
308  }
309  }
310 
311  return status;
312 }
313 
314 
315 void
317  GaussPoint *gp,
318  const FloatArray &totalStrain,
319  TimeStep *tStep)
320 //
321 // returns real stress vector in 3d stress space of receiver according to
322 // previous level of stress and current
323 // strain increment, the only way, how to correctly update gp records
324 //
325 {
326  LatticeDamage2dStatus *status = static_cast< LatticeDamage2dStatus * >( this->giveStatus(gp) );
327 
328  const double e0 = this->give(e0_ID, gp) * this->e0Mean;
329  status->setE0(e0);
330 
331  FloatArray reducedStrain;
332 
333  double f, equivStrain, tempKappa, omega = 0.;
334 
335  this->initTempStatus(gp);
336  reducedStrain = totalStrain;
337 
338  // subtract stress independent part
339  this->giveStressDependentPartOfStrainVector(reducedStrain, gp, totalStrain, tStep, VM_Total);
340 
341  // compute equivalent strain
342  this->computeEquivalentStrain(equivStrain, reducedStrain, gp, tStep);
343 
344 
345  // compute value of loading function if strainLevel crit apply
346  f = equivStrain - status->giveKappa();
347 
348 
349  if ( f <= 0.0 ) {
350  // damage does not grow
351  tempKappa = status->giveKappa();
352  omega = status->giveDamage();
353  if ( status->giveCrackFlag() != 0 ) {
354  status->setTempCrackFlag(2);
355  } else {
356  status->setTempCrackFlag(0);
357  }
358  } else {
359  // damage grows
360  tempKappa = equivStrain;
361 
362  // evaluate damage parameter
363  this->computeDamageParam(omega, tempKappa, reducedStrain, gp);
364  if ( omega > 0 ) {
365  status->setTempCrackFlag(1);
366  }
367  }
368 
369  //Compute the stress
370  answer.resize(3);
371  answer.zero();
372  answer.at(1) = ( 1. - omega ) * eNormal * reducedStrain.at(1);
373  answer.at(2) = ( 1. - omega ) * eShear * reducedStrain.at(2);
374  answer.at(3) = ( 1. - omega ) * eTorsion * reducedStrain.at(3);
375 
376  //Add the water pressure to the normal component of the stress
377  IntArray coupledModels;
378  double waterPressure = 0.;
379 
380 #ifdef __TM_MODULE
382  ( static_cast< StaggeredProblem * >( domain->giveEngngModel()->giveMasterEngngModel() ) )->giveCoupledModels(coupledModels);
383  int couplingFlag = ( static_cast< LatticeStructuralElement * >( gp->giveElement() ) )->giveCouplingFlag();
384 
385  if ( couplingFlag == 1 && coupledModels.at(2) != 0 && !tStep->isTheFirstStep() ) {
386  IntArray couplingNumbers;
387 
388  static_cast< LatticeStructuralElement * >( gp->giveElement() ) ->giveCouplingNumbers(couplingNumbers);
389  LatticeTransportElement *coupledElement;
390  coupledElement = static_cast< LatticeTransportElement * >(domain->giveEngngModel()->giveMasterEngngModel()->giveSlaveProblem( coupledModels.at(2) )->giveDomain(1)->giveElement(couplingNumbers.at(1)) );
391  waterPressure = coupledElement->givePressure();
392  }
393  }
394 #endif
395 
396  //Compute dissipation
397  double tempDissipation = status->giveDissipation();
398  double tempDeltaDissipation = computeDeltaDissipation(omega, reducedStrain, gp, tStep);
399 
400  tempDissipation += tempDeltaDissipation;
401 
402  //Compute crack width
403  double le = ( static_cast< LatticeStructuralElement * >( gp->giveElement() ) )->giveLength();
404  double crackWidth = omega * sqrt( pow(reducedStrain.at(1), 2.) + pow(reducedStrain.at(2), 2.) ) * le;
405 
406  //Calculate the the bio coefficient;
407  double biot = 0.;
408  if ( this->biotType == 0 ) {
409  biot = this->biotCoefficient;
410  } else if ( this->biotType == 1 ) {
411  biot = computeBiot(omega, tempKappa, le);
412  } else {
413  OOFEM_ERROR("Unknown biot type\n");
414  }
415 
416  answer.at(1) = answer.at(1) + biot * waterPressure;
417 
418 
419  //Set all temp values
420  status->setTempDissipation(tempDissipation);
421  status->setTempDeltaDissipation(tempDeltaDissipation);
422 
423  status->setTempEquivalentStrain(equivStrain);
424  status->letTempStrainVectorBe(totalStrain);
425  status->letTempReducedStrainBe(reducedStrain);
426  status->letTempStressVectorBe(answer);
427  status->setTempKappa(tempKappa);
428  status->setTempDamage(omega);
429 
430  status->setTempNormalStress( answer.at(1) );
431  status->setTempCrackWidth(crackWidth);
432 
433  status->setBiotCoefficientInStatus(biot);
434 }
435 
436 double
438  double kappa,
439  double le)
440 {
441 
442  double biot = 0.;
443 
444  if ( this->softeningType == 1 || this->softeningType == 3 ) {
445  if ( omega == 0 ) {
446  biot = this->biotCoefficient;
447  } else if ( omega*kappa*le > 0 && omega*kappa*le < this->wf ) {
448  biot = this->biotCoefficient + (1.-biotCoefficient)*omega*kappa*le/this->wf;
449  } else {
450  biot = 1.;
451  }
452  } else {
453  OOFEM_ERROR("Wrong stype for btype=1. Only linear and exponential softening considered so far\n");
454  }
455 
456  return biot;
457 }
458 
459 
460 double
462  FloatArray &reducedStrain,
463  GaussPoint *gp,
464  TimeStep *tStep)
465 {
466  LatticeDamage2dStatus *status = static_cast< LatticeDamage2dStatus * >( this->giveStatus(gp) );
467  double length = ( static_cast< LatticeStructuralElement * >( gp->giveElement() ) )->giveLength();
468  const double e0 = this->give(e0_ID, gp) * this->e0Mean;
469 
470  FloatArray reducedStrainOld;
471 
472  reducedStrainOld = status->giveReducedStrain();
473  double omegaOld = status->giveDamage();
474  double deltaOmega;
475 
476  FloatArray crackOpeningOld(3);
477  crackOpeningOld.times(omegaOld);
478  crackOpeningOld.times(length);
479  FloatArray intermediateStrain(3);
480 
481  double tempDeltaDissipation = 0.;
482  double deltaTempDeltaDissipation = 0.;
483 
484  double intermediateOmega = 0;
485  FloatArray oldIntermediateStrain(3);
486  oldIntermediateStrain = reducedStrainOld;
487  double oldIntermediateOmega = omegaOld;
488  deltaOmega = ( omega - omegaOld );
489  double testDissipation =
490  0.5 * length * ( pow( ( reducedStrain(0) + reducedStrainOld(0) ) / 2., 2. ) * eNormal +
491  pow( ( reducedStrain(1) + reducedStrainOld(1) ) / 2., 2. ) * eShear +
492  pow( ( reducedStrain(2) + reducedStrainOld(2) ) / 2., 2. ) * eTorsion ) * deltaOmega;
493 
494 
495  double intervals = 0.;
496 
497  double referenceGf = 0;
498 
499  if ( softeningType == 1 ) {
500  referenceGf = e0 * eNormal * this->wf / 2.;
501  } else { //This is for the exponential law. Should also implement it for the bilinear one.
502  referenceGf = e0 * eNormal * this->wf;
503  }
504 
505  if ( testDissipation / ( referenceGf ) > 0.01 ) {
506  intervals = 1000. * testDissipation / referenceGf;
507  } else {
508  intervals = 1.;
509  }
510 
511  if ( intervals > 1000. ) {
512  intervals = 1000.;
513  }
514 
515  double oldKappa = status->giveKappa();
516  double f, equivStrain;
517  if ( deltaOmega > 0 ) {
518  for ( int k = 0; k < intervals; k++ ) {
519  intermediateStrain(0) = reducedStrainOld(0) + ( k + 1 ) / intervals * ( reducedStrain(0) - reducedStrainOld(0) );
520  intermediateStrain(1) = reducedStrainOld(1) + ( k + 1 ) / intervals * ( reducedStrain(1) - reducedStrainOld(1) );
521  intermediateStrain(2) = reducedStrainOld(2) + ( k + 1 ) / intervals * ( reducedStrain(2) - reducedStrainOld(2) );
522  this->computeEquivalentStrain(equivStrain, intermediateStrain, gp, tStep);
523  f = equivStrain - oldKappa;
524  if ( f > 0 ) {
525  this->computeDamageParam(intermediateOmega, equivStrain, intermediateStrain, gp);
526  deltaOmega = ( intermediateOmega - oldIntermediateOmega );
527  deltaTempDeltaDissipation =
528  0.5 * length * ( pow( ( intermediateStrain(0) + oldIntermediateStrain(0) ) / 2., 2. ) * eNormal +
529  pow( ( intermediateStrain(1) + oldIntermediateStrain(1) ) / 2., 2. ) * eShear +
530  pow( ( intermediateStrain(2) + oldIntermediateStrain(2) ) / 2., 2. ) * eTorsion ) * deltaOmega;
531 
532  oldKappa = equivStrain;
533  oldIntermediateOmega = intermediateOmega;
534  } else {
535  deltaTempDeltaDissipation = 0.;
536  }
537 
538  tempDeltaDissipation += deltaTempDeltaDissipation;
539  oldIntermediateStrain = intermediateStrain;
540  }
541  } else {
542  tempDeltaDissipation = 0.;
543  }
544 
545  if ( tempDeltaDissipation >= 2. * referenceGf ) {
546  tempDeltaDissipation = 2. * referenceGf;
547  }
548 
549  return tempDeltaDissipation;
550 }
551 
552 
554 {
555  param.resize(3);
556  param.at(1) = localRandomType;
557 
558  if ( localRandomType == 1 ) { //Gaussian
559  param.at(2) = coefficientOfVariation;
560  } else {
561  OOFEM_ERROR("Unknown local random type:\n randomtype 1 = Gaussian");
562  }
563 }
564 
565 
566 Interface *
568 {
569  return NULL;
570 }
571 
572 
573 void
575  MatResponseMode rMode,
576  GaussPoint *gp, TimeStep *tStep)
577 {
578  MaterialMode mMode = gp->giveMaterialMode();
579  switch ( mMode ) {
580  case _2dLattice:
581 
582  if ( rMode == ElasticStiffness ) {
583  this->giveElasticStiffnessMatrix(answer, gp, tStep);
584  } else if ( rMode == SecantStiffness ) {
585  this->giveSecantStiffnessMatrix(answer, gp, tStep);
586  } else if ( rMode == TangentStiffness ) {
587  this->giveSecantStiffnessMatrix(answer, gp, tStep);
588  } else {
589  OOFEM_ERROR("Unsupported stiffness mode");
590  }
591 
592  break;
593  default:
594  OOFEM_ERROR("unknown material mode");
595  }
596 }
597 
598 
599 void
601  GaussPoint *gp,
602  TimeStep *tStep)
603 {
604  LatticeDamage2dStatus *status = static_cast< LatticeDamage2dStatus * >( this->giveStatus(gp) );
605 
606  double omega = status->giveTempDamage();
607 
608  if ( omega > 1. - 1.e-9 ) { //To avoid convergence problems
609  omega = 1. - 1.e-9;
610  }
611 
612  /* Returns elastic moduli in reduced stress-strain space*/
613  answer.resize(3, 3);
614  answer.zero();
615  answer.at(1, 1) = ( 1 - omega ) * eNormal;
616  answer.at(2, 2) = ( 1 - omega ) * eShear;
617  answer.at(3, 3) = ( 1 - omega ) * eTorsion;
618 }
619 
620 void
622  GaussPoint *gp,
623  TimeStep *tStep)
624 {
625  /* Returns elastic moduli in reduced stress-strain space*/
626  answer.resize(3, 3);
627  answer.zero();
628  answer.at(1, 1) = eNormal;
629  answer.at(2, 2) = eShear;
630  answer.at(3, 3) = eTorsion;
631 }
632 
633 
634 void
636  GaussPoint *gp, TimeStep *tStep)
637 {
638  answer.resize(3);
639  answer.zero();
640  answer.at(1) = this->give(tAlpha, gp);
641 }
642 
643 
644 double
646 {
647  double answer;
648  if ( static_cast< LatticeDamage2dStatus * >( this->giveStatus(gp) )->_giveProperty(aProperty, answer) ) {
649  return answer;
650  } else if ( aProperty == e0_ID ) {
651  return 1.;
652  } else if ( aProperty == ef_ID ) {
653  return 1.;
654  } else {
655  return StructuralMaterial :: give(aProperty, gp);
656  }
657 }
658 
659 
660 int
662  GaussPoint *gp,
663  InternalStateType type,
664  TimeStep *tStep)
665 {
666  LatticeDamage2dStatus *status = static_cast< LatticeDamage2dStatus * >( this->giveStatus(gp) );
667  if ( type == IST_CrackStatuses ) {
668  answer.resize(1);
669  answer.at(1) = status->giveCrackFlag();
670  return 1;
671  } else if ( type == IST_DamageScalar ) {
672  answer.resize(1);
673  answer.at(1) = status->giveDamage();
674  return 1;
675  } else if ( type == IST_DamageTensor ) {
676  answer.resize(6);
677  answer.at(1) = answer.at(2) = answer.at(3) = status->giveDamage();
678  return 1;
679  } else if ( type == IST_DissWork ) {
680  answer.resize(1);
681  answer.at(1) = status->giveDissipation();
682  return 1;
683  } else if ( type == IST_DeltaDissWork ) {
684  answer.resize(1);
685  answer.at(1) = status->giveDeltaDissipation();
686  return 1;
687  } else if ( type == IST_CrackWidth ) {
688  answer.resize(1);
689  answer.at(1) = status->giveCrackWidth();
690  return 1;
691  } else {
692  return StructuralMaterial :: giveIPValue(answer, gp, type, tStep);
693  }
694 }
695 
697  LatticeMaterialStatus(n, d, g), RandomMaterialStatusExtensionInterface(), reducedStrain(3), tempReducedStrain(3)
698 {
699  le = 0.0;
703  oldNormalStress = 0;
704  e0 = 0.;
705  damage = tempDamage = 0.;
707  kappa = tempKappa = 0.;
710  biot = 0.;
711 }
712 
713 void
715 //
716 // initializes temp variables according to variables form previous equlibrium state.
717 // builds new crackMap
718 //
719 {
721  this->tempReducedStrain = this->reducedStrain;
722  this->tempKappa = this->kappa;
723  this->tempEquivStrain = this->equivStrain;
724  this->tempDamage = this->damage;
725  this->tempDissipation = this->dissipation;
727  this->temp_crack_flag = this->crack_flag;
728  this->tempCrackWidth = this->crackWidth;
729  this->tempNormalStress = this->normalStress;
730  this->updateFlag = 0;
731 }
732 
733 void
735 {
737  fprintf(file, "status { ");
738  fprintf(file, "reduced strains ");
739  for ( auto &val : reducedStrain ) {
740  fprintf( file, "%.4e ", val );
741  }
742 
743  fprintf(file, "kappa %f, equivStrain %f, damage %f, dissipation %f, deltaDissipation %f, e0 %f, crack_flag %d, crackWidth %.8e, biotCoeff %f", this->kappa, this->equivStrain, this->damage, this->dissipation, this->deltaDissipation, this->e0, this->crack_flag, this->crackWidth, this->biot);
744  fprintf(file, "}\n");
745 }
746 
747 
748 void
750 {
752 
753  this->reducedStrain = this->tempReducedStrain;
754  this->kappa = this->tempKappa;
755  this->equivStrain = this->tempEquivStrain;
756  this->damage = this->tempDamage;
757  this->dissipation = this->tempDissipation;
759  this->crack_flag = this->temp_crack_flag;
760  this->crackWidth = this->tempCrackWidth;
761  //Need to preserve oldValues for coupling;
762  this->oldNormalStress = this->normalStress;
763  this->normalStress = this->tempNormalStress;
764  //set updateflag
765  this->updateFlag = 1;
766 }
767 
768 Interface *
770 {
772  return static_cast< RandomMaterialStatusExtensionInterface * >(this);
773  } else {
774  return NULL;
775  }
776 }
777 
778 void
780 {
781  tempNormalStress = val;
782 }
783 
784 void
786 {
787  oldNormalStress = val;
788 }
789 
790 
793 //
794 // saves full information stored in this Status
795 // no temp variables stored
796 //
797 {
798  contextIOResultType iores;
799  // save parent class status
800  if ( ( iores = StructuralMaterialStatus :: saveContext(stream, mode, obj) ) != CIO_OK ) {
801  THROW_CIOERR(iores);
802  }
803 
804 
805 
806  // write a raw data
807  if ( ( iores = reducedStrain.storeYourself(stream) ) != CIO_OK ) {
808  THROW_CIOERR(iores);
809  }
810 
811  // write a raw data
812  if ( !stream.write(kappa) ) {
814  }
815 
816  if ( !stream.write(equivStrain) ) {
818  }
819 
820  if ( !stream.write(damage) ) {
822  }
823 
824  if ( !stream.write(le) ) {
826  }
827 
828  if ( !stream.write(dissipation) ) {
830  }
831 
832  if ( !stream.write(deltaDissipation) ) {
834  }
835 
836  if ( !stream.write(crack_flag) ) {
838  }
839 
840  if ( !stream.write(crackWidth) ) {
842  }
843 
844 
845  if ( !stream.write(e0) ) {
847  }
848 
849  if ( !stream.write(biot) ) {
851  }
852 
853  return CIO_OK;
854 }
855 
856 void
858 {
859  temp_crack_flag = val;
860 }
861 
862 
863 
864 void
866 {
867  tempCrackWidth = val;
868 }
869 
870 
871 void
873 {
874  e0 = variable;
875 }
876 
877 
878 void
880 {
881  biot = variable;
882 }
883 
884 int
886 {
887  if ( crack_flag != 0 ) { }
888 
889  return crack_flag;
890 }
891 
892 
895 //
896 // restores full information stored in stream to this Status
897 //
898 {
899  contextIOResultType iores;
900  // read parent class status
901  if ( ( iores = StructuralMaterialStatus :: restoreContext(stream, mode, obj) ) != CIO_OK ) {
902  THROW_CIOERR(iores);
903  }
904 
905  if ( ( iores = reducedStrain.restoreYourself(stream) ) != CIO_OK ) {
906  THROW_CIOERR(iores);
907  }
908 
909  // read raw data
910  if ( !stream.read(kappa) ) {
912  }
913 
914  if ( !stream.read(equivStrain) ) {
916  }
917 
918  if ( !stream.read(damage) ) {
920  }
921 
922  if ( !stream.read(le) ) {
924  }
925 
926  if ( !stream.read(dissipation) ) {
928  }
929 
930  if ( !stream.read(deltaDissipation) ) {
932  }
933 
934  if ( !stream.read(crack_flag) ) {
936  }
937 
938  if ( !stream.read(crackWidth) ) {
940  }
941 
942  if ( !stream.read(e0) ) {
944  }
945 
946  if ( !stream.read(biot) ) {
948  }
949 
950 
951  return CIO_OK;
952 }
953 } // end namespace oofem
virtual EngngModel * giveSlaveProblem(int i)
Returns i-th slave problem.
Definition: engngm.h:1082
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.
MaterialMode giveMaterialMode()
Returns corresponding material mode of receiver.
Definition: gausspoint.h:191
virtual void giveRealStressVector(FloatArray &answer, GaussPoint *, const FloatArray &, TimeStep *)
Computes the real stress vector for given total strain and integration point.
virtual void computeResultingIPTemperatureAt(FloatArray &answer, TimeStep *tStep, GaussPoint *gp, ValueModeType mode)
Computes at given time (tStep) the the resulting temperature component array.
void letTempStrainVectorBe(const FloatArray &v)
Assigns tempStrainVector to given vector v.
Definition: structuralms.h:137
void setTempDamage(double newDamage)
Sets the temp damage level to given value.
double tempNormalStress
nonequilibrated normal stress
double coefficientOfVariation
Coefficient variation of the Gaussian distribution.
Class and object Domain.
Definition: domain.h:115
IntegrationPointStatus * setMaterialStatus(IntegrationPointStatus *ptr, int n)
Sets Material status managed by receiver.
Definition: gausspoint.h:213
double tempKappa
Non-equilibrated scalar measure of the largest strain level.
virtual IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description stored in input record.
Domain * domain
Link to domain object, useful for communicating with other FEM components.
Definition: femcmpnn.h:82
virtual contextIOResultType saveContext(DataStream &stream, ContextMode mode, void *obj=NULL)
Stores receiver state to output stream.
double giveKappa()
Returns the last equilibrated scalar measure of the largest strain level.
#define tAlpha
Definition: matconst.h:67
int temp_crack_flag
Non-equilibrated temp flag.
double tempDeltaDissipation
Non-equilibrated increment of dissipation.
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.
LatticeDamage2dStatus(int n, Domain *d, GaussPoint *g)
Constructor.
double e0
Random material parameter stored in status, since each gp has a different value.
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
double computeDeltaDissipation(double omega, FloatArray &reducedStrain, GaussPoint *gp, TimeStep *tStep)
Compute increment of dissipation for post-processing reasons.
This class implements associated Material Status to LatticeDamage2d.
void clear()
Clears receiver (zero size).
Definition: floatarray.h:206
#define _IFT_LatticeDamage2d_calpha
LatticeDamage2d(int n, Domain *d)
Constructor.
double normalStress
equilibrated normal stress
virtual int giveIPValue(FloatArray &answer, GaussPoint *gp, InternalStateType type, TimeStep *tStep)
Returns the integration point corresponding value in Reduced form.
virtual ~LatticeDamage2d()
Destructor.
double deltaDissipation
Increment of dissipation.
EngngModel * giveEngngModel()
Returns engineering model to which receiver is associated.
Definition: domain.C:433
#define _IFT_LatticeDamage2d_localrandomtype
Dictionary propertyDictionary
Property dictionary.
Definition: material.h:105
double coh
Parameter setting ratio of shear and tensile strength.
void setBiotCoefficientInStatus(double variable)
Element * giveElement()
Returns corresponding element to receiver.
Definition: gausspoint.h:188
General IO error.
This class implements a base lattice material status.
void setTempEquivalentStrain(double newEquivStrain)
Sets the temp scalar measure of the largest strain level to given value.
virtual Interface * giveInterface(InterfaceType)
Interface requesting service.
double giveDeltaDissipation()
Returns the last equilibrated increment of dissipation.
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
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.
virtual void computeEquivalentStrain(double &kappa, const FloatArray &strain, GaussPoint *gp, TimeStep *tStep)
virtual void setVariableInStatus(double variable)
double giveTempDamage()
Returns the temp. damage level.
#define THROW_CIOERR(e)
Definition: contextioerr.h:61
#define _IFT_LatticeDamage2d_alphaTwo
double tempDissipation
Non-equilibrated dissipation..
bool isTheFirstStep()
Check if receiver is first step.
Definition: timestep.C:134
FloatArray reducedStrain
Reduced strain.
void _generateStatusVariables(GaussPoint *) const
Sets up (generates) the variables identified in randVariables array using generators given in randomV...
void letTempReducedStrainBe(FloatArray v)
Assign the temp value of plastic strain.
virtual int write(const int *data, int count)=0
Writes count integer values from array pointed by data.
void setTempNormalStress(double val)
Sets the temp normalStress.
virtual double give(int aProperty, GaussPoint *gp)
Returns the value of material property &#39;aProperty&#39;.
Definition: material.C:52
void setE0(double val)
Set random e0.
virtual IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description stored in input record.
double tempDamage
Non-equilibrated damage level of material.
#define _IFT_IsotropicLinearElasticMaterial_talpha
double tempEquivStrain
Non-equilibrated scalar measure of the strain.
double crackWidth
Crack width.
void setTempDissipation(double newDiss)
Sets the temp dissipation.
Element * giveElement(int n)
Service for accessing particular domain fe element.
Definition: domain.C:160
virtual MaterialStatus * CreateStatus(GaussPoint *gp) const
Creates new copy of associated status and inserts it into given integration point.
void printOutputAt(FILE *file, TimeStep *tStep)
Prints the receiver state to given stream.
double ec
Parameter for the elliptic equivalent strain function.
virtual void computeStressIndependentStrainVector(FloatArray &answer, GaussPoint *gp, TimeStep *tStep, ValueModeType mode)
Computes reduced strain vector in given integration point, generated by internal processes in materia...
#define OOFEM_ERROR(...)
Definition: error.h:61
double le
Characteristic length.
void setTempCrackWidth(double val)
Sets the temp crack width.
#define _IFT_LatticeDamage2d_softeningType
#define _IFT_LatticeDamage2d_alphaOne
void setTempKappa(double newKappa)
Sets the temp scalar measure of the largest strain level to given value.
#define _IFT_LatticeDamage2d_ec
double e0OneMean
Mean effective strain at sigma1.
Abstract base class for all random materials.
contextIOResultType restoreYourself(DataStream &stream)
Definition: floatarray.C:895
virtual double givePressure()
Returns the pressure.
FloatArray tempReducedStrain
Non-equilibrated reduced strain.
Implementation of general sequence (staggered) problem.
#define _IFT_LatticeDamage2d_bio
virtual Interface * giveInterface(InterfaceType)
Interface requesting service.
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 initTempStatus()
Initializes the temporary internal variables, describing the current state according to previously re...
Abstract base class for all random constitutive model statuses.
void setOldNormalStress(double val)
Sets the old normalStress.
void setTempCrackFlag(int val)
Sets the temp_crack_flag.
double giveIntrinsicTime()
Returns intrinsic time, e.g. time in which constitutive model is evaluated.
Definition: timestep.h:148
double equivStrain
Scalar measure of the strain.
EngngModel * giveMasterEngngModel()
Returns the master engnmodel.
Definition: engngm.h:516
double giveDissipation()
Returns the last equilibrated dissipation.
virtual contextIOResultType restoreContext(DataStream &stream, ContextMode mode, void *obj=NULL)
Restores the receiver state previously written in stream.
Definition: structuralms.C:157
virtual void updateYourself(TimeStep *)
Update equilibrium history variables according to temp-variables.
#define _IFT_LatticeDamage2d_wfOne
#define _IFT_LatticeDamage2d_e0OneMean
int crack_flag
The crack_flag indicates if the gp is cracked: crack_flag = 0 gp is uncracked crack_flag = 1 gp is cr...
Abstract base class representing a material status information.
Definition: matstatus.h:84
Pair * add(int aKey, double value)
Adds a new Pair with given keyword and value into receiver.
Definition: dictionary.C:81
Class representing vector of real numbers.
Definition: floatarray.h:82
double damage
Damage level of material.
int biotType
Parameter specifying how the biot coefficient changes with the crack opening.
Implementation of matrix containing floating point numbers.
Definition: floatmatrix.h:94
#define _IFT_LatticeDamage2d_coefficientOfVariation
double oldNormalStress
old normal stress
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
IntegrationPointStatus * giveMaterialStatus()
Returns reference to associated material status (NULL if not defined).
Definition: gausspoint.h:205
double referenceTemperature
Reference temperature (temperature, when material has been built into structure). ...
double kappa
Scalar measure of the largest strain level ever reached in material.
double alphaOne
Ratio of shear and normal modulus.
virtual void computeDamageParam(double &omega, double kappa, const FloatArray &strain, GaussPoint *gp)
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 void giveElasticStiffnessMatrix(FloatMatrix &answer, GaussPoint *gp, TimeStep *tStep)
double eNormal
Normal modulus.
virtual contextIOResultType restoreContext(DataStream &stream, ContextMode mode, void *obj=NULL)
Restores the receiver state previously written in stream.
#define _IFT_LatticeDamage2d_coh
Class Interface.
Definition: interface.h:82
void zero()
Zeroes all coefficients of receiver.
Definition: floatarray.C:658
double eShear
Shear modulus.
void times(double s)
Multiplies receiver with scalar.
Definition: floatarray.C:818
#define _IFT_LatticeDamage2d_e0Mean
virtual int hasMaterialModeCapability(MaterialMode mode)
Tests if material supports material mode.
long ContextMode
Context mode (mask), defining the type of information written/read to/from context.
Definition: contextmode.h:43
virtual void giveThermalDilatationVector(FloatArray &answer, GaussPoint *gp, TimeStep *tStep)
Returns a vector of coefficients of thermal dilatation in direction of each material principal (local...
Abstract base class for all "structural" constitutive models.
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
double localRandomType
Flag which chooses between no distribution (0) and Gaussian distribution (1)
InterfaceType
Enumerative type, used to identify interface type.
Definition: interfacetype.h:43
int softeningType
Parameter which determines the typ of the softeningFunction 1 = linear softening 2 = bilinear softeni...
#define ef_ID
Definition: matconst.h:84
double tempCrackWidth
Non-equilibrated crack width.
virtual double computeBiot(double omega, double kappa, double le)
virtual void giveRandomParameters(FloatArray &param)
double biot
Set biot coefficient.
REGISTER_Material(DummyMaterial)
double castingTime
Casting time.
Definition: material.h:113
#define _IFT_LatticeDamage2d_btype
This class implements the base of a special lattice element following the concepts orginally develope...
virtual MaterialStatus * giveStatus(GaussPoint *gp) const
Returns material status of receiver in given integration point.
#define IR_GIVE_OPTIONAL_FIELD(__ir, __value, __id)
Macro facilitating the use of input record reading methods.
Definition: inputrecord.h:78
double biotCoefficient
Parameter controlling the amount of fluid pressure added to the mechanical stress (Biot&#39;s coefficient...
int giveSize() const
Returns the size of receiver.
Definition: floatarray.h:218
virtual void giveSecantStiffnessMatrix(FloatMatrix &answer, GaussPoint *gp, TimeStep *tStep)
virtual contextIOResultType saveContext(DataStream &stream, ContextMode mode, void *obj=NULL)
Stores receiver state to output stream.
Definition: structuralms.C:133
double cAlpha
coefficient used for modelling eigendisplacements
the oofem namespace is to define a context or scope in which all oofem names are defined.
double e0Mean
Mean effective strain at peak.
double dissipation
Dissipation.
double wf
Determines the softening -> corresponds to threshold of crack opening (not strain) ...
#define IR_GIVE_FIELD(__ir, __value, __id)
Macro facilitating the use of input record reading methods.
Definition: inputrecord.h:69
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
double eTorsion
Torsion modulus.
#define _IFT_LatticeDamage2d_eNormal
virtual void giveStiffnessMatrix(FloatMatrix &answer, MatResponseMode mode, GaussPoint *gp, TimeStep *tStep)
Computes the stiffness matrix for giveRealStressVector of receiver in given integration point...
double giveDamage()
Returns the last equilibrated damage level.
Class representing integration point in finite element program.
Definition: gausspoint.h:93
#define OOFEM_WARNING(...)
Definition: error.h:62
IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description stored in input record.
double alphaTwo
Ratio of torsion and normal modulus.
void setTempDeltaDissipation(double newDiss)
Sets the temp. increment dissipation.
Class representing solution step.
Definition: timestep.h:80
const FloatArray & giveReducedStrain() const
Gives the old equilibrated value of plastic strain.
int giveCrackFlag()
Returns the crack_flag.
#define _IFT_LatticeDamage2d_wf
#define e0_ID
Definition: matconst.h:83
double giveCrackWidth()
Gives the last equilibrated crack width.
virtual double give(int aProperty, GaussPoint *gp)
Returns the value of material property &#39;aProperty&#39;.
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