OOFEM  2.4 OOFEM.org - Object Oriented Finite Element Solver
integrationrule.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
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 "integrationrule.h"
36 #include "material.h"
37 #include "crosssection.h"
38 #include "gausspoint.h"
39 #include "datastream.h"
40 #include "contextioerr.h"
41
42 namespace oofem {
43
44 IntegrationRule :: IntegrationRule(int n, Element *e, int startIndx, int endIndx, bool dynamic)
45 {
46  number = n;
47  elem = e;
48  firstLocalStrainIndx = startIndx;
49  lastLocalStrainIndx = endIndx;
50  isDynamic = dynamic;
52 }
53
55 {
56  number = n;
57  elem = e;
59  isDynamic = false;
61 }
62
63
65 {
66  this->clear();
67 }
68
69
70 void
72 {
73  for ( GaussPoint *gp: *this ) {
74  delete gp;
75  }
76
77  gaussPoints.clear();
78 }
79
80
81 GaussPoint *
83 {
84 # ifdef DEBUG
85  if ( ( i < 0 ) || ( i >= this->giveNumberOfIntegrationPoints() ) ) {
86  OOFEM_ERROR("request out of bounds (%d)", i);
87  }
88
89 # endif
90  return gaussPoints [ i ];
91 }
92
93
94 GaussPoint *
96 {
97  double mindist = -1.;
98  GaussPoint *minGp = NULL;
99  for ( GaussPoint *gp: *this ) {
100  double dist = lcoord.distance_square(gp->giveNaturalCoordinates());
101  if ( dist <= mindist || mindist < 0. ) {
102  mindist = dist;
103  minGp = gp;
104  }
105  }
106  return minGp;
107 }
108
109
110 void
112 // Performs end-of-step operations.
113 {
114  for ( GaussPoint *gp: *this ) {
115  gp->printOutputAt(file, tStep);
116  }
117 }
118
119 void
121 {
123  for ( GaussPoint *gp: *this ) {
124  gp->updateYourself(tStep);
125  }
126 }
127
128
131 {
132  //
133  // saves full context (saves state variables, that completely describe
134  // current state)
135  //
136
137  contextIOResultType iores;
138
139  int isdyn = isDynamic;
140  if ( !stream.write(isdyn) ) {
142  }
143
144  if ( isDynamic ) {
145  mode |= CM_Definition; // store definition if dynamic
146  }
147
148  if ( mode & CM_Definition ) {
149  int numberOfIntegrationPoints = (int)this->gaussPoints.size();
150  if ( !stream.write(numberOfIntegrationPoints) ) {
152  }
153
154  // write first and last integration indices
155  if ( !stream.write(firstLocalStrainIndx) ) {
157  }
158
159  if ( !stream.write(lastLocalStrainIndx) ) {
161  }
162  }
163
164  for ( GaussPoint *gp: *this ) {
165  if ( mode & CM_Definition ) {
166  // write gp weight, coordinates, element number, and material mode
167  double dval = gp->giveWeight();
168  if ( !stream.write(dval) ) {
170  }
171
172  if ( ( iores = gp->giveNaturalCoordinates().storeYourself(stream) ) != CIO_OK ) {
173  THROW_CIOERR(iores);
174  }
175
176  //int ival = gp->giveElement()->giveNumber();
177  //if (!stream.write(&ival,1)) THROW_CIOERR(CIO_IOERR);
178  int mmode = gp->giveMaterialMode();
179  if ( !stream.write(mmode) ) {
181  }
182  }
183
184  // write gp data
185  if ( ( iores = gp->giveCrossSection()->saveIPContext(stream, mode, gp) ) != CIO_OK ) {
186  THROW_CIOERR(iores);
187  }
188  }
189
190  return CIO_OK;
191 }
192
195 {
196  //
197  // restores full element context (saves state variables, that completely describe
198  // current state)
199  //
200
201  contextIOResultType iores;
202  int size;
203
204  int isdyn;
205  if ( !stream.read(isdyn) ) {
207  }
208
209  isDynamic = ( bool ) isdyn;
210
211  if ( isDynamic ) {
212  mode |= CM_Definition; // store definition if dynamic
213  }
214
215  if ( mode & CM_Definition ) {
216  if ( !stream.read(size) ) {
218  }
219
220  // read first and last integration indices
221  if ( !stream.read(firstLocalStrainIndx) ) {
223  }
224
225  if ( !stream.read(lastLocalStrainIndx) ) {
227  }
228
229  this->clear();
230
231  this->gaussPoints.resize(size);
232  }
233
234  int i = 1;
235  for ( GaussPoint *&gp: *this ) {
236  if ( mode & CM_Definition ) {
238  double w;
239  if ( !stream.read(w) ) {
241  }
242
244  FloatArray c;
245  if ( ( iores = c.restoreYourself(stream) ) != CIO_OK ) {
246  THROW_CIOERR(iores);
247  }
248
249  // restore element and material mode
250  //int n;
252  MaterialMode m;
253  int _m;
254  if ( !stream.read(_m) ) {
256  }
257
258  m = ( MaterialMode ) _m;
260
261  gp = new GaussPoint(this, i, std :: move(c), w, m);
262  i++;
263  }
264
266  if ( ( iores = gp->giveCrossSection()->restoreIPContext(stream, mode, gp) ) != CIO_OK ) {
267  THROW_CIOERR(iores);
268  }
269  }
270
271  return CIO_OK;
272 }
273
274
275 int
277  MaterialMode matMode)
278 {
279  intdomain = mode;
280
281  switch ( mode ) {
282  case _Line:
283  return this->SetUpPointsOnLine(nPoints, matMode);
284
285  case _Triangle:
286  return this->SetUpPointsOnTriangle(nPoints, matMode);
287
288  case _Square:
289  return this->SetUpPointsOnSquare(nPoints, matMode);
290
291  case _Cube:
292  return this->SetUpPointsOnCube(nPoints, matMode);
293
294  case _Tetrahedra:
295  return this->SetUpPointsOnTetrahedra(nPoints, matMode);
296
297  case _Wedge:
298  // Limited wrapper for now;
299  if ( nPoints == 6 ) {
300  return this->SetUpPointsOnWedge(3, 2, matMode);
301  } else {
302  return this->SetUpPointsOnWedge(3, 3, matMode);
303  }
304
305  default:
306  OOFEM_ERROR("unknown mode (%d)", mode);
307  }
308
309  return 0;
310 }
311
312 int
314  MaterialMode matMode)
315 {
316  intdomain = mode;
317
318  switch ( mode ) {
319  case _3dDegShell:
320  return this->SetUpPointsOn3dDegShell(nPointsXY, nPointsZ, matMode);
321
322  default:
323  OOFEM_ERROR("Unknown mode (%d)", mode);
324  }
325
326  return 0;
327 }
328
329 int
331  const std :: vector< FloatArray > &coords)
332 {
333  intdomain = mode;
334
335  switch ( mode ) {
336  case _Embedded2dLine:
337  if ( coords.size() != 2 ) {
338  OOFEM_ERROR("Exactly 2 coordinates are required for 2D embedded lines!");
339  }
340  return this->SetUpPointsOn2DEmbeddedLine(nPoints, matMode, coords[0], coords[1]);
341
342  default:
343  OOFEM_ERROR("unknown mode");
344  }
345
346  return 0;
347 }
348
349
351 {
352  this->gaussPoints.resize(1);
353  this->gaussPoints [ 0 ] = new GaussPoint(this, 1, FloatArray(0), 1.0, mode);
354  this->intdomain = _Point;
355  return 1;
356 }
357 } // end namespace oofem
MaterialMode giveMaterialMode()
Returns corresponding material mode of receiver.
Definition: gausspoint.h:191
integrationDomain
Used by integrator class to supply integration points for proper domain to be integrated (Area...
int SetUpPoint(MaterialMode mode)
Trivial implementation, only creates a single point.
bool isDynamic
Flag indicating that rule is dynamic, ie, its gauss points (their number, coordinates, weights) can change during computation.
Element * elem
Element which integration rule is coupled to.
integrationDomain intdomain
Integration domain.
int firstLocalStrainIndx
firstLocalStrainIndx and lastLocalStrainIndx indexes describe range of components (strains for exampl...
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 SetUpPointsOnWedge(int nPointsTri, int nPointsDepth, MaterialMode mode)
Sets up receiver&#39;s integration points on a wedge integration domain.
int setUpEmbeddedIntegrationPoints(integrationDomain intdomain, int nPoints, MaterialMode matMode, const std::vector< FloatArray > &coords)
void clear()
Clears the receiver, ie deallocates all integration points.
virtual int SetUpPointsOnTetrahedra(int, MaterialMode mode)
Sets up receiver&#39;s integration points on tetrahedra (volume coords) integration domain.
virtual int SetUpPointsOn3dDegShell(int nPointsXY, int nPointsZ, MaterialMode mode)
Sets up receiver&#39;s integration points on shell integration domain.
virtual contextIOResultType restoreIPContext(DataStream &stream, ContextMode mode, GaussPoint *gp)
Reads integration point state to output stream.
Definition: crosssection.C:128
Abstract base class for all finite elements.
Definition: element.h:145
General IO error.
virtual int SetUpPointsOnSquare(int, MaterialMode mode)
Sets up receiver&#39;s integration points on unit square integration domain.
MaterialMode
Type representing material mode of integration point.
Definition: materialmode.h:89
virtual int read(int *data, int count)=0
Reads count integer values into array pointed by data.
#define THROW_CIOERR(e)
Definition: contextioerr.h:61
virtual int write(const int *data, int count)=0
Writes count integer values from array pointed by data.
IntegrationRule(int n, Element *e, int startIndx, int endIndx, bool dynamic)
Constructor.
virtual void updateYourself(TimeStep *tStep)
Definition: gausspoint.C:141
#define OOFEM_ERROR(...)
Definition: error.h:61
virtual double giveWeight()
Definition: gausspoint.h:181
virtual int SetUpPointsOn2DEmbeddedLine(int nPoints, MaterialMode mode, const FloatArray &coord0, const FloatArray &coord1)
Sets up integration points on 2D embedded line inside 2D volume (the list of local coordinates should...
CrossSection * giveCrossSection()
Returns reference to cross section associated to related element of receiver.
Definition: gausspoint.h:200
contextIOResultType restoreYourself(DataStream &stream)
Definition: floatarray.C:895
virtual int SetUpPointsOnCube(int, MaterialMode mode)
Sets up receiver&#39;s integration points on unit cube integration domain.
double distance_square(const FloatArray &iP1, const FloatArray &iP2, double &oXi, double &oXiUnbounded) const
Definition: floatarray.C:499
virtual contextIOResultType saveIPContext(DataStream &stream, ContextMode mode, GaussPoint *gp)
Stores integration point state to output stream.
Definition: crosssection.C:110
Class representing vector of real numbers.
Definition: floatarray.h:82
GaussPoint * getIntegrationPoint(int n)
Access particular integration point of receiver.
GaussPoint * findIntegrationPointClosestTo(const FloatArray &lcoord)
Scans through the integration points and finds the one closest to the given (local) coordinate...
virtual void printOutputAt(FILE *file, TimeStep *tStep)
Prints output of receiver to file.
Definition: gausspoint.C:78
int giveNumberOfIntegrationPoints() const
Returns number of integration points of receiver.
virtual int SetUpPointsOnLine(int, MaterialMode mode)
Sets up receiver&#39;s integration points on unit line integration domain.
long ContextMode
Definition: contextmode.h:43
virtual int SetUpPointsOnTriangle(int, MaterialMode mode)
Sets up receiver&#39;s integration points on triangular (area coords) integration domain.
#define CM_Definition
Definition: contextmode.h:47
int setUpIntegrationPoints(integrationDomain intdomain, int nPoints, MaterialMode matMode)
virtual void printOutputAt(FILE *file, TimeStep *tStep)
Prints receiver&#39;s output to given stream.
virtual ~IntegrationRule()
Destructor.
the oofem namespace is to define a context or scope in which all oofem names are defined.
virtual contextIOResultType restoreContext(DataStream &stream, ContextMode mode, void *obj)
virtual contextIOResultType saveContext(DataStream &stream, ContextMode mode, void *obj)
Class representing integration point in finite element program.
Definition: gausspoint.h:93
void updateYourself(TimeStep *tStep)