OOFEM  2.4
OOFEM.org - Object Oriented Finite Element Solver
linquad3d_planestress.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/Elements/PlaneStress/linquad3d_planestress.h"
36 #include "classfactory.h"
37 #include "../sm/Materials/structuralms.h"
38 #include "../sm/CrossSections/structuralcrosssection.h"
39 #include "gausspoint.h"
40 
41 #ifdef __OOFEG
42  #include "oofeggraphiccontext.h"
43  #include "oofegutils.h"
44  #include "connectivitytable.h"
45  #include "Materials/rcm2.h"
46 #endif
47 
48 namespace oofem {
49 REGISTER_Element(LinQuad3DPlaneStress);
50 
52  PlaneStress2d(n, aDomain)
53  // Constructor.
54 {
55  this->GtoLRotationMatrix = NULL;
56 }
57 
59 // Destructor
60 {
61  delete this->GtoLRotationMatrix;
62 }
63 
64 Interface *
66 {
67  if ( interface == ZZNodalRecoveryModelInterfaceType ) {
68  return static_cast< ZZNodalRecoveryModelInterface * >(this);
69  } else if ( interface == SPRNodalRecoveryModelInterfaceType ) {
70  return static_cast< SPRNodalRecoveryModelInterface * >(this);
71  } else if ( interface == SpatialLocalizerInterfaceType ) {
72  return static_cast< SpatialLocalizerInterface * >(this);
73  }
74 
75  return NULL;
76 }
77 
78 
81 {
82  if (cellGeometryWrapper) {
83  return cellGeometryWrapper;
84  } else {
87  }
88 }
89 
90 
91 void
93 // Returns global coordinates given in global vector
94 // transformed into local coordinate system of the
95 // receiver
96 {
97  // test if previously computed
98  if ( GtoLRotationMatrix == NULL ) {
100  }
101 
102 
103  lxy.resize(4);
104  const FloatArray *nc;
105  for ( int i = 0; i < 4; i++ ) {
106  nc = this->giveNode(i + 1)->giveCoordinates();
107  lxy[i].beProductOf(* GtoLRotationMatrix, *nc);
108  }
109 }
110 
111 void
113 {
114  answer = {D_u, D_v, D_w};
115 }
116 
117 
118 
119 
120 const FloatMatrix *
122 // Returns the rotation matrix of the receiver of the size [3,3]
123 // coords(local) = T * coords(global)
124 //
125 // local coordinate (described by vector triplet e1',e2',e3') is defined as follows:
126 //
127 // e1' : [N2-N1] Ni - means i - th node
128 // help : [N3-N1]
129 // e3' : e1' x help
130 // e2' : e3' x e1'
131 {
132  if ( GtoLRotationMatrix == NULL ) {
133  FloatArray e1, e2, e3, help;
134 
135  // compute e1' = [N2-N1] and help = [N3-N1]
136  e1.beDifferenceOf( * this->giveNode(2)->giveCoordinates(), * this->giveNode(1)->giveCoordinates() );
137  help.beDifferenceOf( * this->giveNode(3)->giveCoordinates(), * this->giveNode(1)->giveCoordinates() );
138 
139  // let us normalize e1'
140  e1.normalize();
141 
142  // compute e3' : vector product of e1' x help
143  e3.beVectorProductOf(e1, help);
144  // let us normalize
145  e3.normalize();
146 
147  // now from e3' x e1' compute e2'
148  e2.beVectorProductOf(e3, e1);
149 
150  //
151  GtoLRotationMatrix = new FloatMatrix(3, 3);
152 
153  for ( int i = 1; i <= 3; i++ ) {
154  GtoLRotationMatrix->at(1, i) = e1.at(i);
155  GtoLRotationMatrix->at(2, i) = e2.at(i);
156  GtoLRotationMatrix->at(3, i) = e3.at(i);
157  }
158  }
159 
160  return GtoLRotationMatrix;
161 }
162 
163 
164 bool
166 // Returns the rotation matrix of the receiver of the size [8,12]
167 // r(local) = T * r(global)
168 // for one node (r written transposed): {u,v} = T * {u,v,w}
169 {
170  // test if pereviously computed
171  if ( GtoLRotationMatrix == NULL ) {
173  }
174 
175  answer.resize(8, 12);
176  answer.zero();
177 
178  for ( int i = 1; i <= 3; i++ ) {
179  answer.at(1, i) = answer.at(3, i + 3) = answer.at(5, i + 6) = answer.at(7, i+9) = GtoLRotationMatrix->at(1, i);
180  answer.at(2, i) = answer.at(4, i + 3) = answer.at(6, i + 6) = answer.at(8, i+9) = GtoLRotationMatrix->at(2, i);
181  }
182 
183  return 1;
184 }
185 int
187 // Returns the rotation matrix of the receiver of the size [6,6]
188 // f(local) = T * f(global)
189 {
190  // test if previously computed
191  if ( GtoLRotationMatrix == NULL ) {
193  }
194 
195  answer.resize(6, 6);
196  answer.zero();
197 
198  for ( int i = 1; i <= 3; i++ ) {
199  answer.at(1, i) = answer.at(4, i + 3) = GtoLRotationMatrix->at(1, i);
200  answer.at(2, i) = answer.at(5, i + 3) = GtoLRotationMatrix->at(2, i);
201  answer.at(3, i) = answer.at(6, i + 3) = GtoLRotationMatrix->at(3, i);
202  }
203 
204  return 1;
205 }
206 
207 void
209 // returns characteristic tensor of the receiver at given gp and tStep
210 // strain vector = (Eps_X, Eps_y, Gamma_xy, Kappa_z)
211 {
212  FloatArray charVect;
214 
215  answer.resize(3, 3);
216  answer.zero();
217 
218  if ( ( type == LocalForceTensor ) || ( type == GlobalForceTensor ) ) {
219  //this->computeStressVector(charVect, gp, tStep);
220  charVect = ms->giveStressVector();
221 
222  answer.at(1, 1) = charVect.at(1);
223  answer.at(2, 2) = charVect.at(2);
224  answer.at(1, 2) = charVect.at(3);
225  answer.at(2, 1) = charVect.at(3);
226  } else if ( ( type == LocalMomentTensor ) || ( type == GlobalMomentTensor ) ) {
227  } else if ( ( type == LocalStrainTensor ) || ( type == GlobalStrainTensor ) ) {
228  //this->computeStrainVector(charVect, gp, tStep);
229  charVect = ms->giveStrainVector();
230 
231  answer.at(1, 1) = charVect.at(1);
232  answer.at(2, 2) = charVect.at(2);
233  answer.at(1, 2) = charVect.at(3) / 2.;
234  answer.at(2, 1) = charVect.at(3) / 2.;
235  } else if ( ( type == LocalCurvatureTensor ) || ( type == GlobalCurvatureTensor ) ) {
236  } else {
237  OOFEM_ERROR("unsupported tensor mode");
238  exit(1);
239  }
240 
241  if ( ( type == GlobalForceTensor ) || ( type == GlobalMomentTensor ) ||
242  ( type == GlobalStrainTensor ) || ( type == GlobalCurvatureTensor ) ) {
245  }
246 }
247 
248 
249 int
251 {
252  FloatMatrix globTensor;
253  CharTensor cht;
254 
255  answer.resize(6);
256 
257  if ( type == IST_ShellForceTensor || type == IST_ShellStrainTensor ) {
258  double c = 1.0;
259  if ( type == IST_ShellForceTensor ) {
260  cht = GlobalForceTensor;
261  } else {
262  cht = GlobalStrainTensor;
263  c = 2.0;
264  }
265 
266  this->giveCharacteristicTensor(globTensor, cht, gp, tStep);
267 
268  answer.at(1) = globTensor.at(1, 1); //xx
269  answer.at(2) = globTensor.at(2, 2); //yy
270  answer.at(3) = globTensor.at(3, 3); //zz
271  answer.at(4) = c * globTensor.at(2, 3); //yz
272  answer.at(5) = c * globTensor.at(1, 3); //xz
273  answer.at(6) = c * globTensor.at(1, 2); //xy
274 
275  return 1;
276  } else if ( type == IST_ShellMomentTensor || type == IST_CurvatureTensor ) {
277  answer.clear();
278  return 1;
279  } else {
280  return PlaneStress2d :: giveIPValue(answer, gp, type, tStep);
281  }
282 }
283 
284 
285 void
287 // Performs end-of-step operations.
288 {
289  FloatArray v;
290 
291  fprintf( file, "element %d (%8d) :\n", this->giveLabel(), this->giveNumber() );
292 
293  for ( int i = 0; i < (int)integrationRulesArray.size(); i++ ) {
294  for ( GaussPoint *gp: *integrationRulesArray [ i ] ) {
295 
296  fprintf( file, " GP %2d.%-2d :", i + 1, gp->giveNumber() );
297 
298  this->giveIPValue(v, gp, IST_ShellStrainTensor, tStep);
299  fprintf(file, " strains ");
300  for ( auto &val : v ) fprintf(file, " %.4e", val);
301 
302  /*
303  this->giveIPValue(v, gp, IST_CurvatureTensor, tStep);
304  fprintf(file, "\n curvatures ");
305  for ( auto &val : v ) fprintf(file, " %.4e", val);
306  */
307  // Forces - Moments
308  this->giveIPValue(v, gp, IST_ShellForceTensor, tStep);
309  fprintf(file, "\n forces ");
310  for ( auto &val : v ) fprintf(file, " %.4e", val);
311  /*
312  this->giveIPValue(v, gp, IST_ShellMomentTensor, tStep);
313  fprintf(file, "\n moments ");
314  for ( auto &val : v ) fprintf(file, " %.4e", val);
315  */
316  fprintf(file, "\n");
317  }
318  }
319 }
320 
321 
322 } // end namespace oofem
InternalStateType
Type representing the physical meaning of element or constitutive model internal variable.
This class implements an isoparametric four-node quadrilateral plane- stress elasticity finite elemen...
Definition: planstrss.h:56
The element interface required by ZZNodalRecoveryModel.
virtual int giveIPValue(FloatArray &answer, GaussPoint *gp, InternalStateType type, TimeStep *tStep)
Returns the integration point corresponding value in full form.
void beVectorProductOf(const FloatArray &v1, const FloatArray &v2)
Computes vector product (or cross product) of vectors given as parameters, , and stores the result in...
Definition: floatarray.C:415
Class and object Domain.
Definition: domain.h:115
Wrapper around cell with vertex coordinates stored in FloatArray**.
Definition: feinterpol.h:115
The element interface required by ZZNodalRecoveryModel.
double & at(int i)
Coefficient access function.
Definition: floatarray.h:131
This class implements a structural material status information.
Definition: structuralms.h:65
Class representing a general abstraction for cell geometry.
Definition: feinterpol.h:62
void clear()
Clears receiver (zero size).
Definition: floatarray.h:206
virtual int giveIPValue(FloatArray &answer, GaussPoint *gp, InternalStateType type, TimeStep *tStep)
Returns the integration point corresponding value in full form.
void computeLocalNodalCoordinates(std::vector< FloatArray > &lxy)
virtual int computeLoadGToLRotationMtrx(FloatMatrix &answer)
Returns transformation matrix from global coordinate system to local element coordinate system for el...
Class implementing an array of integers.
Definition: intarray.h:61
void beDifferenceOf(const FloatArray &a, const FloatArray &b)
Sets receiver to be a - b.
Definition: floatarray.C:341
#define OOFEM_ERROR(...)
Definition: error.h:61
REGISTER_Element(LSpace)
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 void giveDofManDofIDMask(int inode, IntArray &) const
Returns dofmanager dof mask for node.
void beProductOf(const FloatMatrix &aMatrix, const FloatArray &anArray)
Receiver becomes the result of the product of aMatrix and anArray.
Definition: floatarray.C:676
virtual void printOutputAt(FILE *file, TimeStep *tStep)
Prints output of receiver to stream, for given time step.
double at(int i, int j) const
Coefficient access function.
Definition: floatmatrix.h:176
std::vector< FloatArray > lc
Local vertex coordinates.
FloatMatrix * GtoLRotationMatrix
Transformation Matrix form GtoL(3,3) is stored at the element level for computation efficiency...
Class representing vector of real numbers.
Definition: floatarray.h:82
Implementation of matrix containing floating point numbers.
Definition: floatmatrix.h:94
const FloatArray & giveStressVector() const
Returns the const pointer to receiver&#39;s stress vector.
Definition: structuralms.h:107
IntegrationPointStatus * giveMaterialStatus()
Returns reference to associated material status (NULL if not defined).
Definition: gausspoint.h:205
void resize(int rows, int cols)
Checks size of receiver towards requested bounds.
Definition: floatmatrix.C:1358
const FloatMatrix * computeGtoLRotationMatrix()
Class Interface.
Definition: interface.h:82
virtual Interface * giveInterface(InterfaceType it)
Interface requesting service.
The spatial localizer element interface associated to spatial localizer.
std::vector< std::unique_ptr< IntegrationRule > > integrationRulesArray
List of integration rules of receiver (each integration rule contains associated integration points a...
Definition: element.h:170
virtual FloatArray * giveCoordinates()
Definition: node.h:114
void zero()
Zeroes all coefficient of receiver.
Definition: floatmatrix.C:1326
InterfaceType
Enumerative type, used to identify interface type.
Definition: interfacetype.h:43
void giveCharacteristicTensor(FloatMatrix &answer, CharTensor type, GaussPoint *gp, TimeStep *tStep)
FEICellGeometry * cellGeometryWrapper
To facilitate the transformation of 2d elements into 3d, the complexity of transformation from 3d to ...
int giveLabel() const
Definition: element.h:1055
the oofem namespace is to define a context or scope in which all oofem names are defined.
int giveNumber() const
Definition: femcmpnn.h:107
double normalize()
Normalizes receiver.
Definition: floatarray.C:828
Node * giveNode(int i) const
Returns reference to the i-th node of element.
Definition: element.h:610
const FloatArray & giveStrainVector() const
Returns the const pointer to receiver&#39;s strain vector.
Definition: structuralms.h:105
Class representing integration point in finite element program.
Definition: gausspoint.h:93
virtual FEICellGeometry * giveCellGeometryWrapper()
Returns the Cell Geometry Wrapper.
Class representing solution step.
Definition: timestep.h:80
void resize(int s)
Resizes receiver towards requested size.
Definition: floatarray.C:631
LinQuad3DPlaneStress(int n, Domain *d)

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