OOFEM  2.4
OOFEM.org - Object Oriented Finite Element Solver
structural3delement.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 "feinterpol3d.h"
37 #include "gausspoint.h"
39 #include "gaussintegrationrule.h"
40 #include "mathfem.h"
41 
42 namespace oofem {
44  NLStructuralElement(n, aDomain),
45  matRotation(false)
46 {
47 }
48 
49 
52 {
55 }
56 
57 
58 void
60 // Returns the [ 6 x (nno*3) ] strain-displacement matrix {B} of the receiver, eva-
61 // luated at gp.
62 // B matrix - 6 rows : epsilon-X, epsilon-Y, epsilon-Z, gamma-YZ, gamma-ZX, gamma-XY :
63 {
64  FEInterpolation *interp = this->giveInterpolation();
65  FloatMatrix dNdx;
66  interp->evaldNdx( dNdx, gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(this) );
67 
68  answer.resize(6, dNdx.giveNumberOfRows() * 3);
69  answer.zero();
70 
71  for ( int i = 1; i <= dNdx.giveNumberOfRows(); i++ ) {
72  answer.at(1, 3 * i - 2) = dNdx.at(i, 1);
73  answer.at(2, 3 * i - 1) = dNdx.at(i, 2);
74  answer.at(3, 3 * i - 0) = dNdx.at(i, 3);
75 
76  answer.at(5, 3 * i - 2) = answer.at(4, 3 * i - 1) = dNdx.at(i, 3);
77  answer.at(6, 3 * i - 2) = answer.at(4, 3 * i - 0) = dNdx.at(i, 2);
78  answer.at(6, 3 * i - 1) = answer.at(5, 3 * i - 0) = dNdx.at(i, 1);
79  }
80 }
81 
82 
83 
84 void
86 // Returns the [ 9 x (nno * 3) ] displacement gradient matrix {BH} of the receiver,
87 // evaluated at gp.
88 // BH matrix - 9 rows : du/dx, dv/dy, dw/dz, dv/dz, du/dz, du/dy, dw/dy, dw/dx, dv/dx
89 {
90  FEInterpolation *interp = this->giveInterpolation();
91  FloatMatrix dNdx;
92  interp->evaldNdx( dNdx, gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(this) );
93 
94  answer.resize(9, dNdx.giveNumberOfRows() * 3);
95  answer.zero();
96 
97  for ( int i = 1; i <= dNdx.giveNumberOfRows(); i++ ) {
98  answer.at(1, 3 * i - 2) = dNdx.at(i, 1); // du/dx
99  answer.at(2, 3 * i - 1) = dNdx.at(i, 2); // dv/dy
100  answer.at(3, 3 * i - 0) = dNdx.at(i, 3); // dw/dz
101  answer.at(4, 3 * i - 1) = dNdx.at(i, 3); // dv/dz
102  answer.at(7, 3 * i - 0) = dNdx.at(i, 2); // dw/dy
103  answer.at(5, 3 * i - 2) = dNdx.at(i, 3); // du/dz
104  answer.at(8, 3 * i - 0) = dNdx.at(i, 1); // dw/dx
105  answer.at(6, 3 * i - 2) = dNdx.at(i, 2); // du/dy
106  answer.at(9, 3 * i - 1) = dNdx.at(i, 1); // dv/dx
107  }
108 
109 }
110 
111 
114 {
115  return _3dMat;
116 }
117 
118 
119 void
121 {
122  if ( this->elemLocalCS.isNotEmpty() ) { // User specified orientation
123  x.beColumnOf(this->elemLocalCS, 1);
124  y.beColumnOf(this->elemLocalCS, 2);
125  z.beColumnOf(this->elemLocalCS, 3);
126  } else {
128  FloatMatrix jac;
129  FloatArray help;
130  this->giveInterpolation()->giveJacobianMatrixAt( jac, lcoords, FEIElementGeometryWrapper(this) );
131  x.beColumnOf(jac, 1); // This is {dx/dxi, dy/dxi, dz/dxi}
132  x.normalize();
133  help.beColumnOf(jac, 2);
134  z.beVectorProductOf(x, help); // Normal to the xi-eta plane.
135  z.normalize();
136  y.beVectorProductOf(z, x);
137  }
138 }
139 
140 
141 void
143 {
144  if ( this->matRotation ) {
146  FloatArray x, y, z;
147  FloatArray rotStrain, s;
148 
149  this->giveMaterialOrientationAt( x, y, z, gp->giveNaturalCoordinates() );
150  // Transform from global c.s. to material c.s.
151  rotStrain = {
152  e(0) * x(0) * x(0) + e(5) * x(0) * x(1) + e(1) * x(1) * x(1) + e(4) * x(0) * x(2) + e(3) * x(1) * x(2) + e(2) * x(2) * x(2),
153  e(0) * y(0) * y(0) + e(5) * y(0) * y(1) + e(1) * y(1) * y(1) + e(4) * y(0) * y(2) + e(3) * y(1) * y(2) + e(2) * y(2) * y(2),
154  e(0) * z(0) * z(0) + e(5) * z(0) * z(1) + e(1) * z(1) * z(1) + e(4) * z(0) * z(2) + e(3) * z(1) * z(2) + e(2) * z(2) * z(2),
155  2 * e(0) * y(0) * z(0) + e(4) * y(2) * z(0) + 2 * e(1) * y(1) * z(1) + e(3) * y(2) * z(1) + e(5) * ( y(1) * z(0) + y(0) * z(1) ) + ( e(4) * y(0) + e(3) * y(1) + 2 * e(2) * y(2) ) * z(2),
156  2 * e(0) * x(0) * z(0) + e(4) * x(2) * z(0) + 2 * e(1) * x(1) * z(1) + e(3) * x(2) * z(1) + e(5) * ( x(1) * z(0) + x(0) * z(1) ) + ( e(4) * x(0) + e(3) * x(1) + 2 * e(2) * x(2) ) * z(2),
157  2 * e(0) * x(0) * y(0) + e(4) * x(2) * y(0) + 2 * e(1) * x(1) * y(1) + e(3) * x(2) * y(1) + e(5) * ( x(1) * y(0) + x(0) * y(1) ) + ( e(4) * x(0) + e(3) * x(1) + 2 * e(2) * x(2) ) * y(2)
158  };
159  this->giveStructuralCrossSection()->giveRealStress_3d(s, gp, rotStrain, tStep);
160  answer = {
161  s(0) * x(0) * x(0) + 2 * s(5) * x(0) * y(0) + s(1) * y(0) * y(0) + 2 * ( s(4) * x(0) + s(3) * y(0) ) * z(0) + s(2) * z(0) * z(0),
162  s(0) * x(1) * x(1) + 2 * s(5) * x(1) * y(1) + s(1) * y(1) * y(1) + 2 * ( s(4) * x(1) + s(3) * y(1) ) * z(1) + s(2) * z(1) * z(1),
163  s(0) * x(2) * x(2) + 2 * s(5) * x(2) * y(2) + s(1) * y(2) * y(2) + 2 * ( s(4) * x(2) + s(3) * y(2) ) * z(2) + s(2) * z(2) * z(2),
164  y(2) * ( s(5) * x(1) + s(1) * y(1) + s(3) * z(1) ) + x(2) * ( s(0) * x(1) + s(5) * y(1) + s(4) * z(1) ) + ( s(4) * x(1) + s(3) * y(1) + s(2) * z(1) ) * z(2),
165  y(2) * ( s(5) * x(0) + s(1) * y(0) + s(3) * z(0) ) + x(2) * ( s(0) * x(0) + s(5) * y(0) + s(4) * z(0) ) + ( s(4) * x(0) + s(3) * y(0) + s(2) * z(0) ) * z(2),
166  y(1) * ( s(5) * x(0) + s(1) * y(0) + s(3) * z(0) ) + x(1) * ( s(0) * x(0) + s(5) * y(0) + s(4) * z(0) ) + ( s(4) * x(0) + s(3) * y(0) + s(2) * z(0) ) * z(1)
167  };
168  } else {
169  this->giveStructuralCrossSection()->giveRealStress_3d(answer, gp, e, tStep);
170  }
171 }
172 
173 void
175 {
176  this->giveStructuralCrossSection()->giveStiffnessMatrix_3d(answer, rMode, gp, tStep);
177  if ( this->matRotation ) {
178  FloatArray x, y, z;
179  FloatMatrix Q;
180 
181  this->giveMaterialOrientationAt( x, y, z, gp->giveNaturalCoordinates() );
182 
183  Q = {
184  { x(0) * x(0), x(1) * x(1), x(2) * x(2), x(1) * x(2), x(0) * x(2), x(0) * x(1) },
185  { y(0) * y(0), y(1) * y(1), y(2) * y(2), y(1) * y(2), y(0) * y(2), y(0) * y(1) },
186  { z(0) * z(0), z(1) * z(1), z(2) * z(2), z(1) * z(2), z(0) * z(2), z(0) * z(1) },
187  { 2 * y(0) * z(0), 2 * y(1) * z(1), 2 * y(2) * z(2), y(2) * z(1) + y(1) * z(2), y(2) * z(0) + y(0) * z(2), y(1) * z(0) + y(0) * z(1) },
188  { 2 * x(0) * z(0), 2 * x(1) * z(1), 2 * x(2) * z(2), x(2) * z(1) + x(1) * z(2), x(2) * z(0) + x(0) * z(2), x(1) * z(0) + x(0) * z(1) },
189  { 2 * x(0) * y(0), 2 * x(1) * y(1), 2 * x(2) * y(2), x(2) * y(1) + x(1) * y(2), x(2) * y(0) + x(0) * y(2), x(1) * y(0) + x(0) * y(1) }
190  };
191  answer.rotatedWith(Q, 't');
192  }
193 }
194 
195 
196 void
198 {
199  answer = {D_u, D_v, D_w};
200 }
201 
202 
203 int
205 {
207  IntArray dofIdMask;
208  this->giveDofManDofIDMask(-1, dofIdMask); // ok for standard elements
209  return this->giveInterpolation()->giveNumberOfNodes() * dofIdMask.giveSize();
210 }
211 
212 
214 // Sets up the array containing the four Gauss points of the receiver.
215 {
216  if ( integrationRulesArray.size() == 0 ) {
217  integrationRulesArray.resize( 1 );
218  integrationRulesArray [ 0 ].reset( new GaussIntegrationRule(1, this, 1, 6) );
220  }
221 }
222 
223 
224 
225 double
227 // Returns the portion of the receiver which is attached to gp.
228 {
229  double determinant, weight, volume;
230  determinant = fabs( this->giveInterpolation()->giveTransformationJacobian( gp->giveNaturalCoordinates(),
231  FEIElementGeometryWrapper(this) ) );
232 
233  weight = gp->giveWeight();
234  volume = determinant * weight;
235 
236  return volume;
237 }
238 
239 
240 
241 double
243 {
244  return this->giveLengthInDir(normalToCrackPlane);
245 }
246 
247 
248 
249 
250 
251 // Surface support
252 void
254 {
255  /* Returns the [ 3 x (nno*3) ] shape function matrix {N} of the receiver,
256  * evaluated at the given gp.
257  * {u} = {N}*{a} gives the displacements at the integration point.
258  */
259 
260  // Evaluate the shape functions at the position of the gp.
261  FloatArray N;
262  static_cast< FEInterpolation3d* > ( this->giveInterpolation() )->
263  surfaceEvalN( N, iSurf, sgp->giveNaturalCoordinates(), FEIElementGeometryWrapper(this) );
264  answer.beNMatrixOf(N, 3);
265 }
266 
267 void
269 {
270  IntArray nodes;
271  const int ndofsn = 3;
272 
273  static_cast< FEInterpolation3d* > ( this->giveInterpolation() )->
274  computeLocalSurfaceMapping(nodes, iSurf);
275 
276  answer.resize(nodes.giveSize() *3 );
277 
278  for ( int i = 1; i <= nodes.giveSize(); i++ ) {
279  answer.at(i * ndofsn - 2) = nodes.at(i) * ndofsn - 2;
280  answer.at(i * ndofsn - 1) = nodes.at(i) * ndofsn - 1;
281  answer.at(i * ndofsn) = nodes.at(i) * ndofsn;
282  }
283 }
284 
285 double
287 {
288  double determinant, weight, volume;
289  determinant = fabs( static_cast< FEInterpolation3d* > ( this->giveInterpolation() )->
290  surfaceGiveTransformationJacobian( iSurf, gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(this) ) );
291 
292  weight = gp->giveWeight();
293  volume = determinant * weight;
294 
295  return volume;
296 }
297 
298 
299 
300 int
302 {
303  OOFEM_ERROR("surface local coordinate system not supported");
304  return 1;
305 }
306 
307 
308 
309 
310 
311 
312 // Edge support
313 
314 void
316 {
317  /*
318  * provides dof mapping of local edge dofs (only nonzero are taken into account)
319  * to global element dofs
320  */
321  IntArray eNodes;
322  static_cast< FEInterpolation3d* > ( this->giveInterpolation() )->computeLocalEdgeMapping(eNodes, iEdge);
323 
324  answer.resize( eNodes.giveSize() * 3 );
325  for ( int i = 1; i <= eNodes.giveSize(); i++ ) {
326  answer.at(i * 3 - 2) = eNodes.at(i) * 3 - 2;
327  answer.at(i * 3 - 1) = eNodes.at(i) * 3 - 1;
328  answer.at(i * 3) = eNodes.at(i) * 3;
329  }
330 }
331 
332 
333 
334 double
336 {
337  /* Returns the line element ds associated with the given gp on the specific edge.
338  * Note: The name is misleading since there is no volume to speak of in this case.
339  * The returned value is used for integration of a line integral (external forces).
340  */
341  double detJ = static_cast< FEInterpolation3d* > ( this->giveInterpolation() )->
342  edgeGiveTransformationJacobian( iEdge, gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(this) );
343  return detJ * gp->giveWeight();
344 }
345 
346 
347 int
349 {
350  // returns transformation matrix from
351  // edge local coordinate system
352  // to element local c.s
353  // (same as global c.s in this case)
354  //
355  // i.e. f(element local) = T * f(edge local)
356  //
358  OOFEM_ERROR("egde local coordinate system not supported");
359  return 1;
360 }
361 
362 
363 
364 } // end namespace oofem
CrossSection * giveCrossSection()
Definition: element.C:495
virtual void giveDofManDofIDMask(int inode, IntArray &answer) const
Returns dofmanager dof mask for node.
virtual MaterialMode giveMaterialMode()
Returns material mode for receiver integration points.
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
virtual void computeBHmatrixAt(GaussPoint *gp, FloatMatrix &answer)
Computes a matrix which, multiplied by the column matrix of nodal displacements, gives the displaceme...
Abstract base class for "structural" finite elements with geometrical nonlinearities.
virtual bool hasField(InputFieldType id)=0
Returns true if record contains field identified by idString keyword.
void giveMaterialOrientationAt(FloatArray &x, FloatArray &y, FloatArray &z, const FloatArray &lcoords)
virtual void giveSurfaceDofMapping(IntArray &answer, int) const
Assembles surface dof mapping mask, which provides mapping between surface local DOFs and "global" el...
MaterialMode
Type representing material mode of integration point.
Definition: materialmode.h:89
Class implementing an array of integers.
Definition: intarray.h:61
int & at(int i)
Coefficient access function.
Definition: intarray.h:103
MatResponseMode
Describes the character of characteristic material matrix.
virtual void computeSurfaceNMatrixAt(FloatMatrix &answer, int iSurf, GaussPoint *gp)
virtual int giveNumberOfNodes() const
Returns the number of geometric nodes of the receiver.
Definition: feinterpol.h:486
virtual FEInterpolation * giveInterpolation() const
Definition: element.h:629
virtual double computeEdgeVolumeAround(GaussPoint *gp, int iEdge)
Computes volume related to integration point on local edge.
void beColumnOf(const FloatMatrix &mat, int col)
Reciever will be set to a given column in a matrix.
Definition: floatarray.C:1114
virtual void computeConstitutiveMatrixAt(FloatMatrix &answer, MatResponseMode rMode, GaussPoint *gp, TimeStep *tStep)
Computes constitutive matrix of receiver.
Class representing a general abstraction for finite element interpolation class.
Definition: feinterpol.h:132
bool isNotEmpty() const
Tests for empty matrix.
Definition: floatmatrix.h:162
FloatMatrix elemLocalCS
Transformation material matrix, used in orthotropic and anisotropic materials, global->local transfor...
Definition: element.h:173
virtual void computeGaussPoints()
Initializes the array of integration rules member variable.
#define OOFEM_ERROR(...)
Definition: error.h:61
virtual double giveWeight()
Returns integration weight of receiver.
Definition: gausspoint.h:181
StructuralCrossSection * giveStructuralCrossSection()
Helper function which returns the structural cross-section for the element.
void rotatedWith(const FloatMatrix &r, char mode= 'n')
Returns the receiver &#39;a&#39; transformed using give transformation matrix r.
Definition: floatmatrix.C:1557
Wrapper around element definition to provide FEICellGeometry interface.
Definition: feinterpol.h:95
#define N(p, q)
Definition: mdm.C:367
virtual double giveCharacteristicLength(const FloatArray &normalToCrackPlane)
Returns the size of element in the given direction, in some cases adjusted (e.g.
virtual int computeLoadLEToLRotationMatrix(FloatMatrix &answer, int iEdge, GaussPoint *gp)
Returns transformation matrix from local edge c.s to element local coordinate system of load vector c...
Structural3DElement(int n, Domain *d)
Constructor.
virtual double giveLengthInDir(const FloatArray &normalToCrackPlane)
Default implementation returns length of element projection into specified direction.
Definition: element.C:1143
double at(int i, int j) const
Coefficient access function.
Definition: floatmatrix.h:176
void resize(int n)
Checks size of receiver towards requested bounds.
Definition: intarray.C:124
virtual int setupIntegrationPoints(IntegrationRule &irule, int npoints, Element *element)
Sets up integration rule for the given element.
Definition: crosssection.C:54
int numberOfGaussPoints
Number of integration points as specified by nip.
Definition: element.h:188
void beNMatrixOf(const FloatArray &n, int nsd)
Assigns the receiver to be a repeated diagonal matrix.
Definition: floatmatrix.C:505
virtual void giveRealStress_3d(FloatArray &answer, GaussPoint *gp, const FloatArray &reducedStrain, TimeStep *tStep)=0
virtual int computeLoadLSToLRotationMatrix(FloatMatrix &answer, int, GaussPoint *gp)
Returns transformation matrix from local surface c.s to element local coordinate system of load vecto...
Class representing vector of real numbers.
Definition: floatarray.h:82
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
virtual void giveJacobianMatrixAt(FloatMatrix &jacobianMatrix, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Gives the jacobian matrix at the local coordinates.
Definition: feinterpol.h:232
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
Class representing a general abstraction for surface finite element interpolation class...
Definition: feinterpol3d.h:44
virtual double computeSurfaceVolumeAround(GaussPoint *gp, int)
Computes volume related to integration point on local surface.
virtual void computeBmatrixAt(GaussPoint *gp, FloatMatrix &answer, int lowerIndx=1, int upperIndx=ALL_STRAINS)
Computes the geometrical matrix of receiver in given integration point.
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 IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description stored in input record.
virtual int computeNumberOfDofs()
Computes or simply returns total number of element&#39;s local DOFs.
void zero()
Zeroes all coefficient of receiver.
Definition: floatmatrix.C:1326
virtual void giveEdgeDofMapping(IntArray &answer, int iEdge) const
Assembles edge dof mapping mask, which provides mapping between edge local DOFs and "global" element ...
virtual double computeVolumeAround(GaussPoint *gp)
Returns volume related to given integration point.
int giveSize() const
Definition: intarray.h:203
virtual void giveStiffnessMatrix_3d(FloatMatrix &answer, MatResponseMode mode, GaussPoint *gp, TimeStep *tStep)=0
Method for computing the stiffness matrix.
the oofem namespace is to define a context or scope in which all oofem names are defined.
virtual void computeStressVector(FloatArray &answer, const FloatArray &strain, GaussPoint *gp, TimeStep *tStep)
Computes the stress vector of receiver at given integration point, at time step tStep.
double normalize()
Normalizes receiver.
Definition: floatarray.C:828
virtual IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description stored in input record.
int giveNumberOfRows() const
Returns number of rows of receiver.
Definition: floatmatrix.h:156
Class representing integration point in finite element program.
Definition: gausspoint.h:93
#define _IFT_Structural3DElement_materialCoordinateSystem
[optional] Support for material directions based on element orientation.
Class representing solution step.
Definition: timestep.h:80
virtual double evaldNdx(FloatMatrix &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo)=0
Evaluates the matrix of derivatives of interpolation functions (shape functions) at given point...
const FloatArray & giveNaturalCoordinates()
Returns coordinate array of receiver.
Definition: gausspoint.h:138
Class representing Gaussian-quadrature integration rule.

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