OOFEM  2.4
OOFEM.org - Object Oriented Finite Element Solver
deidynamic.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 "../sm/EngineeringModels/deidynamic.h"
36 #include "timestep.h"
37 #include "dofmanager.h"
38 #include "element.h"
39 #include "dof.h"
40 #include "verbose.h"
41 #include "mathfem.h"
42 #include "classfactory.h"
43 #include "unknownnumberingscheme.h"
44 
45 namespace oofem {
46 #define ZERO_MASS 1.E-10 // unit dependent !!!!
47 
48 REGISTER_EngngModel(DEIDynamic);
49 
51 
53 // only one has reason for DEIDynamic
54 // - SolutionOfLinearEquations
55 
56 {
57  return NULL; // not necessary here - diagonal matrix is used-simple inversion
58 }
59 
60 
63 {
64  IRResultType result; // Required by IR_GIVE_FIELD macro
65 
66  IR_GIVE_FIELD(ir, dumpingCoef, _IFT_DEIDynamic_dumpcoef); // C = dumpingCoef * M
68 
70 }
71 
72 
74 // returns unknown quantity like displacement, velocity of equation eq in time t
75 // This function translates this request to numerical method language
76 {
77  int eq = dof->__giveEquationNumber();
78 #ifdef DEBUG
79  if ( eq == 0 ) {
80  OOFEM_ERROR("invalid equation number");
81  }
82 #endif
83 
84  if ( tStep != this->giveCurrentStep() ) {
85  OOFEM_ERROR("unknown time step encountered");
86  return 0.;
87  }
88 
89  switch ( mode ) {
90  case VM_Total:
91  return displacementVector.at(eq);
92 
93  case VM_Velocity:
94  return velocityVector.at(eq);
95 
96  case VM_Acceleration:
97  return accelerationVector.at(eq);
98 
99  default:
100  OOFEM_ERROR("Unknown is of undefined ValueModeType for this problem");
101  }
102 
103  return 0.0;
104 }
105 
106 
108 {
109  int istep = giveNumberOfFirstStep();
110  double totalTime = 0.;
111  StateCounterType counter = 1;
112 
113  if ( currentStep ) {
114  totalTime = currentStep->giveTargetTime() + deltaT;
115  istep = currentStep->giveNumber() + 1;
116  counter = currentStep->giveSolutionStateCounter() + 1;
117  }
118 
119  previousStep = std :: move(currentStep);
120  currentStep.reset( new TimeStep(istep, this, 1, totalTime, deltaT, counter) );
121 
122  return currentStep.get();
123 }
124 
125 
127 {
128  //
129  // creates system of governing eq's and solves them at given time step
130  //
131  // this is an explicit problem: we assemble governing equating at time t
132  // and solution is obtained for time t+dt
133  //
134  // first assemble problem at current time step to obtain results in following
135  // time step.
136  // and then print results for this step also.
137  // for first time step we need special start code
138  Domain *domain = this->giveDomain(1);
139  int nelem = domain->giveNumberOfElements();
140  IntArray loc;
141  Element *element;
142  int neq;
143  int n, init = 0;
144  double coeff, maxDt, maxOmi, maxOm = 0., maxOmEl, c1, c2, c3;
145  FloatMatrix charMtrx, charMtrx2, R;
146  FloatArray previousDisplacementVector;
147 
148 
150  if ( tStep->giveNumber() == giveNumberOfFirstStep() ) {
151  init = 1;
152 #ifdef VERBOSE
153  OOFEM_LOG_INFO("Assembling mass matrix\n");
154 #endif
155 
156  //
157  // first step assemble mass Matrix
158  //
159 
160  massMatrix.resize(neq);
161  massMatrix.zero();
163  for ( int i = 1; i <= nelem; i++ ) {
164  element = domain->giveElement(i);
165  element->giveLocationArray(loc, dn);
166  element->giveCharacteristicMatrix(charMtrx, LumpedMassMatrix, tStep);
167  // charMtrx.beLumpedOf(fullCharMtrx);
168  element->giveCharacteristicMatrix(charMtrx2, TangentStiffnessMatrix, tStep);
169  if ( charMtrx2.isNotEmpty() ) {
171  if ( element->giveRotationMatrix(R) ) {
172  charMtrx2.rotatedWith(R);
173  charMtrx.rotatedWith(R);
174  }
175  }
176 
177  //
178  // assemble it manually
179  //
180 #ifdef DEBUG
181  if ( loc.giveSize() != charMtrx.giveNumberOfRows() ) {
182  OOFEM_ERROR("dimension mismatch");
183  }
184 
185 #endif
186 
187  n = loc.giveSize();
188 
189  maxOmEl = 0.;
190  for ( int j = 1; j <= n; j++ ) {
191  if ( charMtrx.at(j, j) > ZERO_MASS ) {
192  maxOmi = charMtrx2.at(j, j) / charMtrx.at(j, j);
193  if ( init ) {
194  maxOmEl = ( maxOmEl > maxOmi ) ? ( maxOmEl ) : ( maxOmi );
195  }
196  }
197  }
198 
199  maxOm = ( maxOm > maxOmEl ) ? ( maxOm ) : ( maxOmEl );
200 
201  if (maxOmEl > ZERO_MASS) {
202  // update zero mass diagonal entries ; but skip massless elements
203  for ( int j = 1; j <= n; j++ ) {
204  int jj = loc.at(j);
205  if ( ( jj ) && ( charMtrx.at(j, j) <= ZERO_MASS ) ) {
206  charMtrx.at(j, j) = charMtrx2.at(j, j) / maxOmEl;
207  }
208  }
209  }
210 
211  for ( int j = 1; j <= n; j++ ) {
212  int jj = loc.at(j);
213  if ( jj ) {
214  massMatrix.at(jj) += charMtrx.at(j, j);
215  }
216  }
217  }
218 
219  // if init - try to determine the best deltaT
220  if ( init ) {
221  maxDt = 2 / sqrt(maxOm);
222  if ( deltaT > maxDt ) {
223  OOFEM_LOG_RELEVANT("DEIDynamic: deltaT reduced to %e\n", maxDt);
224  deltaT = maxDt;
225  tStep->setTimeIncrement(deltaT);
226  }
227  }
228 
229 
230  //
231  // special init step - compute displacements at tstep 0
232  //
237  velocityVector.resize(neq);
241 
242 
243  for ( auto &node : domain->giveDofManagers() ) {
244  for ( Dof *iDof: *node ) {
245  // ask for initial values obtained from
246  // bc (boundary conditions) and ic (initial conditions)
247  // now we are setting initial cond. for step -1.
248  if ( !iDof->isPrimaryDof() ) {
249  continue;
250  }
251 
252  int jj = iDof->__giveEquationNumber();
253  if ( jj ) {
254  nextDisplacementVector.at(jj) = iDof->giveUnknown(VM_Total, tStep);
255  // become displacementVector after init
256  velocityVector.at(jj) = iDof->giveUnknown(VM_Velocity, tStep);
257  // accelerationVector = iDof->giveUnknown(AccelerartionVector,tStep) ;
258  }
259  }
260  }
261 
263 
264  return;
265  } // end of init step
266 
267 #ifdef VERBOSE
268  OOFEM_LOG_INFO("Assembling right hand side\n");
269 #endif
270 
271 
272  c1 = ( 1. / ( deltaT * deltaT ) );
273  c2 = ( 1. / ( 2. * deltaT ) );
274  c3 = ( 2. / ( deltaT * deltaT ) );
275 
276  previousDisplacementVector = displacementVector;
278 
279  //
280  // assembling the element part of load vector
281  //
283  loadVector.zero();
285  VM_Total, EModelDefaultEquationNumbering(), domain);
287 
288 
289  //
290  // assembling additional parts of right hand side
291  //
293  for ( int i = 1; i <= nelem; i++ ) {
294  element = domain->giveElement(i);
295  element->giveLocationArray(loc, dn);
296  element->giveCharacteristicMatrix(charMtrx, TangentStiffnessMatrix, tStep);
297  if ( charMtrx.isNotEmpty() ) {
299  if ( element->giveRotationMatrix(R) ) {
300  charMtrx.rotatedWith(R);
301  }
302  }
303 
304  n = loc.giveSize();
305  for ( int j = 1; j <= n; j++ ) {
306  int jj = loc.at(j);
307  if ( jj ) {
308  for ( int k = 1; k <= n; k++ ) {
309  int kk = loc.at(k);
310  if ( kk ) {
311  loadVector.at(jj) -= charMtrx.at(j, k) * displacementVector.at(kk);
312  }
313  }
314 
315  //
316  // if init step - find minimum period of vibration in order to
317  // determine maximal admissible time step
318  //
319  //maxOmi = charMtrx.at(j,j)/massMatrix.at(jj) ;
320  //if (init) maxOm = (maxOm > maxOmi) ? (maxOm) : (maxOmi) ;
321  }
322  }
323  }
324 
325 
326 
327  for ( int j = 1; j <= neq; j++ ) {
328  coeff = massMatrix.at(j);
329  loadVector.at(j) += coeff * c3 * displacementVector.at(j) -
330  coeff * ( c1 - dumpingCoef * c2 ) *
331  previousDisplacementVector.at(j);
332  }
333 
334  //
335  // set-up numerical model
336  //
337  /* it is not necessary to call numerical method
338  * approach used here is not good, but effective enough
339  * inverse of diagonal mass matrix is done here
340  */
341  //
342  // call numerical model to solve arose problem - done locally here
343  //
344 #ifdef VERBOSE
345  OOFEM_LOG_RELEVANT( "Solving [step number %8d, time %15e]\n", tStep->giveNumber(), tStep->giveTargetTime() );
346 #endif
347  double prevD;
348 
349  for ( int i = 1; i <= neq; i++ ) {
350  prevD = previousDisplacementVector.at(i);
352  ( massMatrix.at(i) * ( c1 + dumpingCoef * c2 ) );
353  velocityVector.at(i) = nextDisplacementVector.at(i) - prevD;
356  2. * displacementVector.at(i) + prevD;
357  }
358 
360  velocityVector.times(c2);
361 }
362 
363 
364 void DEIDynamic :: printDofOutputAt(FILE *stream, Dof *iDof, TimeStep *tStep)
365 {
366  static char dofchar[] = "dva";
367  static ValueModeType dofmodes[] = {
368  VM_Total, VM_Velocity, VM_Acceleration
369  };
370 
371  iDof->printMultipleOutputAt(stream, tStep, dofchar, dofmodes, 3);
372 }
373 } // end namespace oofem
The representation of EngngModel default unknown numbering.
FloatArray velocityVector
Definition: deidynamic.h:79
std::unique_ptr< TimeStep > currentStep
Current time step.
Definition: engngm.h:231
Class and object Domain.
Definition: domain.h:115
virtual int giveNumberOfDomainEquations(int di, const UnknownNumberingScheme &num)
Returns number of equations for given domain in active (current time step) time step.
Definition: engngm.C:391
virtual IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description in input reader.
Definition: engngm.C:262
Class representing meta step.
Definition: metastep.h:62
std::unique_ptr< TimeStep > previousStep
Previous time step.
Definition: engngm.h:233
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
long StateCounterType
StateCounterType type used to indicate solution state.
virtual double giveUnknownComponent(ValueModeType type, TimeStep *tStep, Domain *d, Dof *dof)
Returns requested unknown.
Definition: deidynamic.C:73
virtual void printDofOutputAt(FILE *stream, Dof *iDof, TimeStep *tStep)
DOF printing routine.
Definition: deidynamic.C:364
This base class is an abstraction for numerical algorithm.
Definition: nummet.h:80
double giveTargetTime()
Returns target time.
Definition: timestep.h:146
void setTimeIncrement(double newDt)
Sets solution step time increment.
Definition: timestep.h:152
Abstract base class for all finite elements.
Definition: element.h:145
std::vector< std::unique_ptr< DofManager > > & giveDofManagers()
Definition: domain.h:400
virtual int __giveEquationNumber() const =0
Returns equation number of receiver, usually assigned by emodel.
int giveNumberOfElements() const
Returns number of elements in domain.
Definition: domain.h:434
virtual void giveCharacteristicMatrix(FloatMatrix &answer, CharType type, TimeStep *tStep)
Computes characteristic matrix of receiver of requested type in given time step.
Definition: element.C:569
virtual void printMultipleOutputAt(FILE *File, TimeStep *tStep, char *ch, ValueModeType *mode, int nite)
Prints Dof output (it prints value of unknown related to dof at given timeStep).
Definition: dof.C:92
REGISTER_EngngModel(ProblemSequence)
virtual int giveNumberOfFirstStep(bool force=false)
Returns number of first time step used by receiver.
Definition: deidynamic.h:102
Class implementing an array of integers.
Definition: intarray.h:61
int & at(int i)
Coefficient access function.
Definition: intarray.h:103
#define _IFT_DEIDynamic_dumpcoef
Definition: deidynamic.h:43
#define OOFEM_LOG_RELEVANT(...)
Definition: logger.h:126
FloatArray nextDisplacementVector
Definition: deidynamic.h:78
#define ZERO_MASS
Definition: deidynamic.C:46
bool isNotEmpty() const
Tests for empty matrix.
Definition: floatmatrix.h:162
#define OOFEM_LOG_INFO(...)
Definition: logger.h:127
int giveNumber()
Returns receiver&#39;s number.
Definition: timestep.h:129
Element * giveElement(int n)
Service for accessing particular domain fe element.
Definition: domain.C:160
double dumpingCoef
Definition: deidynamic.h:80
#define OOFEM_ERROR(...)
Definition: error.h:61
virtual ~DEIDynamic()
Definition: deidynamic.C:50
FloatArray loadVector
Definition: deidynamic.h:77
int updateSharedDofManagers(FloatArray &answer, const UnknownNumberingScheme &s, int ExchangeTag)
Exchanges necessary remote DofManagers data.
Definition: engngm.C:1957
FloatArray accelerationVector
Definition: deidynamic.h:79
Implementation for assembling external forces vectors in standard monolithic FE-problems.
virtual TimeStep * giveNextStep()
Returns next time step (next to current step) of receiver.
Definition: deidynamic.C:107
void giveLocationArray(IntArray &locationArray, const UnknownNumberingScheme &s, IntArray *dofIds=NULL) const
Returns the location array (array of code numbers) of receiver for given numbering scheme...
Definition: element.C:390
void rotatedWith(const FloatMatrix &r, char mode= 'n')
Returns the receiver &#39;a&#39; transformed using give transformation matrix r.
Definition: floatmatrix.C:1557
virtual bool giveRotationMatrix(FloatMatrix &answer)
Transformation matrices updates rotation matrix between element-local and primary DOFs...
Definition: element.C:259
double at(int i, int j) const
Coefficient access function.
Definition: floatmatrix.h:176
Initializes the variable VERBOSE, in order to get a few intermediate messages on screen: beginning an...
virtual void solveYourselfAt(TimeStep *tStep)
Solves problem for given time step.
Definition: deidynamic.C:126
Class representing vector of real numbers.
Definition: floatarray.h:82
virtual NumericalMethod * giveNumericalMethod(MetaStep *mStep)
Returns reference to receiver&#39;s numerical method.
Definition: deidynamic.C:52
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
FloatArray displacementVector
Definition: deidynamic.h:79
virtual void init()
Initializes the receiver state.
Definition: engngm.C:1883
Class representing the general Input Record.
Definition: inputrecord.h:101
void zero()
Zeroes all coefficients of receiver.
Definition: floatarray.C:658
void times(double s)
Multiplies receiver with scalar.
Definition: floatarray.C:818
#define _IFT_DEIDynamic_deltat
Definition: deidynamic.h:44
virtual IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description in input reader.
Definition: deidynamic.C:62
virtual TimeStep * giveCurrentStep(bool force=false)
Returns current time step.
Definition: engngm.h:683
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
Abstract class Dof represents Degree Of Freedom in finite element mesh.
Definition: dof.h:93
int giveNumberOfRows() const
Returns number of rows of receiver.
Definition: floatmatrix.h:156
void assembleVector(FloatArray &answer, TimeStep *tStep, const VectorAssembler &va, ValueModeType mode, const UnknownNumberingScheme &s, Domain *domain, FloatArray *eNorms=NULL)
Assembles characteristic vector of required type from dofManagers, element, and active boundary condi...
Definition: engngm.C:986
Class representing solution step.
Definition: timestep.h:80
FloatArray massMatrix
Definition: deidynamic.h:76
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:28 for OOFEM by doxygen 1.8.11 written by Dimitri van Heesch, © 1997-2011