OOFEM  2.4
OOFEM.org - Object Oriented Finite Element Solver
structuralmaterialevaluator.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 
36 #include "inputrecord.h"
37 #include "timestep.h"
38 #include "domain.h"
39 #include "gausspoint.h"
40 #include "../sm/Materials/structuralmaterial.h"
41 #include "../sm/Materials/structuralms.h"
42 #include "function.h"
43 #include "classfactory.h"
44 
45 #include <fstream>
46 
47 namespace oofem {
48 REGISTER_EngngModel(StructuralMaterialEvaluator);
49 
51 {
52  this->ndomains = 1;
53 }
54 
56 { }
57 
59 {
60  IRResultType result;
61 
62  this->deltaT = 1.0;
65 
66  //IR_GIVE_FIELD(ir, this->ndim, _IFT_StructuralMaterialEvaluator_nDimensions);
67 
71 
72  tolerance = 1.0;
73  if ( this->sControl.giveSize() > 0 ) {
75  }
76 
78 
79  this->suppressOutput = true;
80 
81  // Compute the strain control (everything not controlled by stress)
82  for ( int i = 1; i <= 6; ++i ) {
83  if ( !sControl.contains(i) ) {
85  }
86  }
87 
88  return IRRT_OK;
89 }
90 
91 
93 {
94  Domain *d = this->giveDomain(1);
95 
96  MaterialMode mode = _3dMat;
97  FloatArray initialStrain(6);
98  gps.clear();
99  gps.reserve(d->giveNumberOfMaterialModels());
100  for ( int i = 1; i <= d->giveNumberOfMaterialModels(); i++ ) {
101  std :: unique_ptr< GaussPoint > gp(new GaussPoint(nullptr, i, FloatArray(0), 1, mode));
102  gps.emplace_back( std :: move(gp) );
103  // Initialize the strain vector;
104  StructuralMaterialStatus *status = static_cast< StructuralMaterialStatus * >( d->giveMaterial(i)->giveStatus( gps[i-1].get() ) );
105  status->letStrainVectorBe(initialStrain);
106  }
107 
108  std :: string outname = this->giveOutputBaseFileName() + ".matdata";
109  this->outfile.open( outname.c_str() );
110 
112 
113  TimeStep *tStep = giveNextStep();
114 
115  // Note, strain == strain-rate (kept as strain for brevity)
116  int maxiter = 100; // User input?
117  FloatArray stressC, deltaStrain, strain, stress, res;
118  stressC.resize( sControl.giveSize() );
119  res.resize( sControl.giveSize() );
120 
121  FloatMatrix tangent, reducedTangent;
122  for ( int istep = 1; istep <= this->numberOfSteps; ++istep ) {
124  for ( int imat = 1; imat <= d->giveNumberOfMaterialModels(); ++imat ) {
125  GaussPoint *gp = gps[imat-1].get();
126  StructuralMaterial *mat = static_cast< StructuralMaterial * >( d->giveMaterial(imat) );
127  StructuralMaterialStatus *status = static_cast< StructuralMaterialStatus * >( mat->giveStatus(gp) );
128 
129  strain = status->giveStrainVector();
130  // Update the controlled parts
131  for ( int j = 1; j <= eControl.giveSize(); ++j ) {
132  int p = eControl.at(j);
133  strain.at(p) = d->giveFunction( cmpntFunctions.at(p) )->evaluateAtTime( tStep->giveIntrinsicTime() );
134  }
135 
136  for ( int j = 1; j <= sControl.giveSize(); ++j ) {
137  int p = sControl.at(j);
138  stressC.at(j) = d->giveFunction( cmpntFunctions.at(p) )->evaluateAtTime( tStep->giveIntrinsicTime() );
139  }
140 
141  //strain.add(-100, {6.27e-06, 6.27e-06, 6.27e-06, 0, 0, 0});
142  for ( int iter = 1; iter < maxiter; iter++ ) {
143 #if 0
144  // Debugging:
145  mat->give3dMaterialStiffnessMatrix(tangent, TangentStiffness, gp, tStep);
146  tangent.printYourself("# tangent");
147 
148  strain.zero();
149  mat->giveRealStressVector_3d(stress, gp, strain, tStep);
150  FloatArray strain2;
151  tangent.solveForRhs(stress, strain2);
152  strain2.printYourself("# thermal expansion");
153  break;
154 #endif
155 
156  strain.printYourself("Macro strain guess");
157  mat->giveRealStressVector_3d(stress, gp, strain, tStep);
158  for ( int j = 1; j <= sControl.giveSize(); ++j ) {
159  res.at(j) = stressC.at(j) - stress.at( sControl.at(j) );
160  }
161 
162  OOFEM_LOG_INFO("*** Time step: %d (t = %.2e), Material %d, Iteration: %d, Residual = %e (tolerance %.2e)\n",
163  istep, tStep->giveIntrinsicTime(), imat, iter, res.computeNorm(), tolerance);
164 
165  if ( res.computeNorm() <= tolerance ) {
166  break;
167  } else {
168  if ( tangent.giveNumberOfRows() == 0 || !keepTangent ) {
169  mat->give3dMaterialStiffnessMatrix(tangent, TangentStiffness, gp, tStep);
170  }
171 
172  // Pick out the stress-controlled part;
173  reducedTangent.beSubMatrixOf(tangent, sControl, sControl);
174 
175  // Update stress-controlled part of the strain
176  reducedTangent.solveForRhs(res, deltaStrain);
177  //deltaStrain.printYourself("deltaStrain");
178  for ( int j = 1; j <= sControl.giveSize(); ++j ) {
179  strain.at( sControl.at(j) ) += deltaStrain.at(j);
180  }
181  }
182  }
183 
184  if ( res.computeNorm() > tolerance ) {
185  OOFEM_WARNING("Residual did not converge!");
186  }
187 
188  // This material model has converged, so we update it and go on to the next.
189  gp->updateYourself(tStep);
190  }
191 
193  this->doStepOutput(tStep);
194  tStep = giveNextStep();
195  }
196 
198  this->outfile.close();
199 }
200 
202 {
203  Domain *d = this->giveDomain(1);
204  for ( auto &mat : d->giveMaterials() ) {
205  if ( !dynamic_cast< StructuralMaterial * >( mat.get() ) ) {
206  OOFEM_LOG_ERROR("Material %d is not a StructuralMaterial", mat->giveNumber());
207  return 0;
208  }
209  }
210 
212 }
213 
215 {
216  FloatArray outputValue;
217  Domain *d = this->giveDomain(1);
218  if ( tStep->isTheFirstStep() ) {
219  this->outfile << "# Time";
220  for ( int var : this->vars ) {
221  this->outfile << ", " << __InternalStateTypeToString( ( InternalStateType ) var );
222  }
223 
224  this->outfile << '\n';
225  }
226 
227  outfile << tStep->giveIntrinsicTime();
228  for ( int i = 1; i <= d->giveNumberOfMaterialModels(); i++ ) {
229  Material *mat = d->giveMaterial(i);
230  for ( int var : this->vars ) {
231  mat->giveIPValue(outputValue, gps[i-1].get(), ( InternalStateType ) var, tStep);
232  outfile << " " << outputValue;
233  }
234  }
235 
236  outfile << std :: endl;
237 }
238 
240 {
241  if ( !currentStep ) {
242  // first step -> generate initial step
243  //currentStep.reset( new TimeStep(*giveSolutionStepWhenIcApply()) );
244  currentStep.reset( new TimeStep(giveNumberOfTimeStepWhenIcApply(), this, 1, 0., this->deltaT, 0) );
245  }
246  previousStep = std :: move(currentStep);
247  currentStep.reset( new TimeStep(*previousStep, this->deltaT) );
248 
249  return currentStep.get();
250 }
251 } // end namespace oofem
bool contains(int value) const
Definition: intarray.h:283
EngngModelTimer timer
E-model timer.
Definition: engngm.h:267
InternalStateType
Type representing the physical meaning of element or constitutive model internal variable.
#define _IFT_StructuralMaterialEvaluator_outputVariables
Variables (from integration point) to be written.
std::unique_ptr< TimeStep > currentStep
Current time step.
Definition: engngm.h:231
virtual MaterialStatus * giveStatus(GaussPoint *gp) const
Returns material status of receiver in given integration point.
Definition: material.C:244
std::string giveOutputBaseFileName()
Returns base output file name to which extensions, like .out .vtu .osf should be added.
Definition: engngm.h:363
virtual void doStepOutput(TimeStep *tStep)
Prints the ouput of the solution step (using virtual this->printOutputAtservice) to the stream detemi...
Class and object Domain.
Definition: domain.h:115
#define OOFEM_LOG_ERROR(...)
Definition: logger.h:122
void letStrainVectorBe(const FloatArray &v)
Assigns strain vector to given vector v.
Definition: structuralms.h:125
bool solveForRhs(const FloatArray &b, FloatArray &answer, bool transpose=false)
Solves the system of linear equations .
Definition: floatmatrix.C:1112
virtual TimeStep * giveNextStep()
Returns next time step (next to current step) of receiver.
void startTimer(EngngModelTimerType t)
Definition: timer.h:128
#define _IFT_StructuralMaterialEvaluator_stressControl
Integer list of the stress components which are controlled.
std::unique_ptr< TimeStep > previousStep
Previous time step.
Definition: engngm.h:233
double & at(int i)
Coefficient access function.
Definition: floatarray.h:131
void beSubMatrixOf(const FloatMatrix &src, int topRow, int bottomRow, int topCol, int bottomCol)
Assigns to the receiver the sub-matrix of another matrix.
Definition: floatmatrix.C:962
This class implements a structural material status information.
Definition: structuralms.h:65
virtual int checkConsistency()
Allows programmer to test some receiver&#39;s internal data, before computation begins.
virtual bool hasField(InputFieldType id)=0
Returns true if record contains field identified by idString keyword.
void stopTimer(EngngModelTimerType t)
Definition: timer.h:129
MaterialMode
Type representing material mode of integration point.
Definition: materialmode.h:89
REGISTER_EngngModel(ProblemSequence)
int & at(int i)
Coefficient access function.
Definition: intarray.h:103
bool isTheFirstStep()
Check if receiver is first step.
Definition: timestep.C:134
virtual void updateYourself(TimeStep *tStep)
Updates internal state of receiver after finishing time step.
Definition: gausspoint.C:141
#define OOFEM_LOG_INFO(...)
Definition: logger.h:127
Material * giveMaterial(int n)
Service for accessing particular domain material model.
Definition: domain.C:281
virtual IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description in input reader.
std::vector< std::unique_ptr< Material > > & giveMaterials()
Definition: domain.h:344
int numberOfSteps
Total number of time steps.
Definition: engngm.h:209
int ndomains
Number of receiver domains.
Definition: engngm.h:205
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 ...
#define _IFT_StructuralMaterialEvaluator_componentFunctions
Integer list of time functions for each component.
int giveNumberOfMaterialModels() const
Returns number of material models in domain.
Definition: domain.h:436
Abstract base class for all material models.
Definition: material.h:95
double giveIntrinsicTime()
Returns intrinsic time, e.g. time in which constitutive model is evaluated.
Definition: timestep.h:148
#define _IFT_StructuralMaterialEvaluator_numberOfTimeSteps
virtual void giveRealStressVector_3d(FloatArray &answer, GaussPoint *gp, const FloatArray &reducedE, TimeStep *tStep)
Default implementation relies on giveRealStressVector for second Piola-Kirchoff stress.
Function * giveFunction(int n)
Service for accessing particular domain load time function.
Definition: domain.C:268
#define _IFT_StructuralMaterialEvaluator_deltat
Class representing vector of real numbers.
Definition: floatarray.h:82
virtual int checkConsistency()
Allows programmer to test some receiver&#39;s internal data, before computation begins.
Definition: engngm.h:995
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
std::vector< std::unique_ptr< GaussPoint > > gps
IntArray sControl
Time functions controlling each component of the deviatoric part of the stress.
int giveNumberOfTimeStepWhenIcApply()
Returns the time step number, when initial conditions should apply.
Definition: engngm.h:754
double computeNorm() const
Computes the norm (or length) of the vector.
Definition: floatarray.C:840
Class representing the general Input Record.
Definition: inputrecord.h:101
virtual void printYourself() const
Print receiver on stdout.
Definition: floatarray.C:747
void followedBy(const IntArray &b, int allocChunk=0)
Appends array b at the end of receiver.
Definition: intarray.C:145
#define _IFT_StructuralMaterialEvaluator_keepTangent
void zero()
Zeroes all coefficients of receiver.
Definition: floatarray.C:658
#define _IFT_StructuralMaterialEvaluator_tolerance
Tolerance for stress control.
void printYourself() const
Prints matrix to stdout. Useful for debugging.
Definition: floatmatrix.C:1458
Abstract base class for all "structural" constitutive models.
bool suppressOutput
Flag for suppressing output to file.
Definition: engngm.h:314
virtual void solveYourself()
Starts solution process.
const char * __InternalStateTypeToString(InternalStateType _value)
Definition: cltypes.C:298
Abstract base class representing the "problem" under consideration.
Definition: engngm.h:181
#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
the oofem namespace is to define a context or scope in which all oofem names are defined.
Domain * giveDomain(int n)
Service for accessing particular problem domain.
Definition: engngm.C:1720
#define IR_GIVE_FIELD(__ir, __value, __id)
Macro facilitating the use of input record reading methods.
Definition: inputrecord.h:69
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
StructuralMaterialEvaluator(int i, EngngModel *_master=NULL)
Class representing integration point in finite element program.
Definition: gausspoint.h:93
#define OOFEM_WARNING(...)
Definition: error.h:62
Class representing solution step.
Definition: timestep.h:80
virtual int giveIPValue(FloatArray &answer, GaussPoint *gp, InternalStateType type, TimeStep *tStep)
Returns the integration point corresponding value in Reduced form.
Definition: material.C:142
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:31 for OOFEM by doxygen 1.8.11 written by Dimitri van Heesch, © 1997-2011