OOFEM  2.4
OOFEM.org - Object Oriented Finite Element Solver
truss2d.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 "Elements/Bars/truss2d.h"
37 #include "node.h"
38 #include "material.h"
39 #include "gausspoint.h"
40 #include "gaussintegrationrule.h"
41 #include "floatmatrix.h"
42 #include "floatarray.h"
43 #include "intarray.h"
44 #include "mathfem.h"
45 #include "classfactory.h"
46 #include "fei2dlinelin.h"
47 
48 #ifdef __OOFEG
49  #include "oofeggraphiccontext.h"
50 #endif
51 
52 namespace oofem {
53 REGISTER_Element(Truss2d);
54 
55  FEI2dLineLin Truss2d :: interp[3] = {FEI2dLineLin(1,3),
56  FEI2dLineLin(1,2),
57  FEI2dLineLin(2,3)};
58 
59 Truss2d :: Truss2d(int n, Domain *aDomain) :
60  NLStructuralElement(n, aDomain)
61 {
62  numberOfDofMans = 2;
63  length = 0.;
64  pitch = 10.; // a dummy value
65  cs_mode = 0;
66 }
67 
68 
69 void
70 Truss2d :: computeBmatrixAt(GaussPoint *gp, FloatMatrix &answer, int li, int ui)
71 //
72 // Returns linear part of geometrical equations of the receiver at gp.
73 // Returns the linear part of the B matrix
74 //
75 {
76  // eps = B*a = [-cos -sin cos sin]/l *[u1 v1 u2 v2]^t
77  // cos = (x2 - x1) / l
78  // sin = (z2 - z1) / l
79 
80  double l, x1, x2, z1, z2;
81  //determine in which plane the truss is defined
82  int c1 = 0, c2 = 0;
83  resolveCoordIndices(c1, c2);
84 
85  x1 = this->giveNode(1)->giveCoordinate(c1);
86  z1 = this->giveNode(1)->giveCoordinate(c2);
87  x2 = this->giveNode(2)->giveCoordinate(c1);
88  z2 = this->giveNode(2)->giveCoordinate(c2);
89 
90  answer.resize(1, 4);
91 
92  answer.at(1, 1) = x1 - x2;
93  answer.at(1, 2) = z1 - z2;
94  answer.at(1, 3) = x2 - x1;
95  answer.at(1, 4) = z2 - z1;
96 
97  l = this->computeLength();
98  answer.times(1.0 / l / l);
99 }
100 
101 
102 void
104 {
105  // Will be the same as the regular B-matrix
106  this->computeBmatrixAt(gp, answer);
107 }
108 
109 
110 
112 // Sets up the array of Gauss Points of the receiver.
113 {
114  if ( integrationRulesArray.size() == 0 ) {
115  integrationRulesArray.resize( 1 );
116  integrationRulesArray [ 0 ].reset( new GaussIntegrationRule(1, this, 1, 2) );
118  }
119 }
120 
121 
122 void
124 // Returns the lumped mass matrix of the receiver. This expression is
125 // valid in both local and global axes.
126 {
127  answer.resize(4, 4);
128  answer.zero();
129  if ( !isActivated(tStep) ) {
130  return;
131  }
132 
133  GaussPoint *gp = integrationRulesArray [ 0 ]->getIntegrationPoint(0);
134  double density = this->giveStructuralCrossSection()->give('d', gp);
135  double halfMass = density * this->giveCrossSection()->give(CS_Area, gp) * this->computeLength() * 0.5;
136  answer.at(1, 1) = halfMass;
137  answer.at(2, 2) = halfMass;
138  answer.at(3, 3) = halfMass;
139  answer.at(4, 4) = halfMass;
140 }
141 
142 
143 void
145 // Returns the displacement interpolation matrix {N} of the receiver, eva-
146 // luated at gp.
147 {
148  double ksi, n1, n2;
149 
150  ksi = iLocCoord.at(1);
151  n1 = ( 1. - ksi ) * 0.5;
152  n2 = ( 1. + ksi ) * 0.5;
153 
154  answer.resize(2, 4);
155  answer.zero();
156  answer.at(1, 1) = n1;
157  answer.at(1, 3) = n2;
158  answer.at(2, 2) = n1;
159  answer.at(2, 4) = n2;
160 }
161 
162 
163 int
165 {
166  //determine in which plane the truss is defined
167  int c1 = 0, c2 = 0;
168  resolveCoordIndices(c1, c2);
169 
170  double ksi, n1, n2;
171 
172  ksi = lcoords.at(1);
173  n1 = ( 1. - ksi ) * 0.5;
174  n2 = ( 1. + ksi ) * 0.5;
175 
176  answer.resize(3);
177  answer.zero();
178  answer.at(c1) = n1 * this->giveNode(1)->giveCoordinate(c1) + n2 *this->giveNode(2)->giveCoordinate(c1);
179  answer.at(c2) = n1 * this->giveNode(1)->giveCoordinate(c2) + n2 *this->giveNode(2)->giveCoordinate(c2);
180 
181  return 1;
182 }
183 
184 
186 // Returns the length of the receiver. This method is valid only if 1
187 // Gauss point is used.
188 {
189  double weight = gp->giveWeight();
190  return 0.5 * this->computeLength() * weight * this->giveCrossSection()->give(CS_Area, gp);
191 }
192 
193 
195 // Returns the length of the receiver.
196 {
197  //determine in which plane the truss is defined
198  int c1 = 0, c2 = 0;
199  resolveCoordIndices(c1, c2);
200 
201  double dx, dz;
202  Node *nodeA, *nodeB;
203 
204  if ( length == 0. ) {
205  nodeA = this->giveNode(1);
206  nodeB = this->giveNode(2);
207  dx = nodeB->giveCoordinate(c1) - nodeA->giveCoordinate(c1);
208  dz = nodeB->giveCoordinate(c2) - nodeA->giveCoordinate(c2);
209  length = sqrt(dx * dx + dz * dz);
210  }
211 
212  return length;
213 }
214 
215 
217 // Returns the pitch of the receiver.
218 {
219  double xA, xB, zA, zB;
220  Node *nodeA, *nodeB;
221  //determine in which plane the truss is defined
222  int c1 = 0, c2 = 0;
223  resolveCoordIndices(c1, c2);
224 
225  if ( pitch == 10. ) { // 10. : dummy initialization value
226  nodeA = this->giveNode(1);
227  nodeB = this->giveNode(2);
228  xA = nodeA->giveCoordinate(c1);
229  xB = nodeB->giveCoordinate(c1);
230  zA = nodeA->giveCoordinate(c2);
231  zB = nodeB->giveCoordinate(c2);
232  pitch = atan2(zB - zA, xB - xA);
233  }
234 
235  return pitch;
236 }
237 
238 
239 int
241 //
242 // returns a unit vectors of local coordinate system at element
243 // stored rowwise (mainly used by some materials with ortho and anisotrophy)
244 //
245 {
246  double sine, cosine;
247 
248  answer.resize(3, 3);
249  answer.zero();
250 
251  sine = sin( this->givePitch() );
252  cosine = cos(pitch);
253 
254  answer.at(1, 1) = cosine;
255  answer.at(1, 2) = sine;
256  answer.at(2, 1) = -sine;
257  answer.at(2, 2) = cosine;
258  answer.at(3, 3) = 1.0;
259 
260  return 1;
261 }
262 
263 
264 void
266 {
267  if ( cs_mode == 0 ) {
268  //xz-plane
269  c1 = 1;
270  c2 = 3;
271  } else if ( cs_mode == 1 ) {
272  //xy-plane
273  c1 = 1;
274  c2 = 2;
275  } else if ( cs_mode == 2 ) {
276  //yz-plane
277  c1 = 2;
278  c2 = 3;
279  } else {
280  OOFEM_ERROR("Unknow cs_mode");
281  }
282 }
283 
286 {
287  IRResultType result; // Required by IR_GIVE_FIELD macro
288 
289  cs_mode = 0;
291 
292  if ( cs_mode != 0 && cs_mode != 1 && cs_mode != 2 ) {
293  OOFEM_WARNING("Unsupported value of cs_mode");
294  return IRRT_BAD_FORMAT;
295  }
296 
298 }
299 
300 
301 void
302 Truss2d :: giveDofManDofIDMask(int inode, IntArray &answer) const
303 {
304  if ( cs_mode == 0 ) {
305  answer = {D_u, D_w};
306  } else if ( cs_mode == 1 ) {
307  answer = {D_u, D_v};
308  } else if ( cs_mode == 2 ) {
309  answer = {D_v, D_w};
310  }
311 }
312 
313 
314 void
316 {
317  this->giveStructuralCrossSection()->giveRealStress_1d(answer, gp, strain, tStep);
318 }
319 
320 
321 void
323 {
324  this->giveStructuralCrossSection()->giveStiffnessMatrix_1d(answer, rMode, gp, tStep);
325 }
326 
327 
328 void
329 Truss2d :: giveEdgeDofMapping(IntArray &answer, int iEdge) const
330 {
331  /*
332  * provides dof mapping of local edge dofs (only nonzero are taken into account)
333  * to global element dofs
334  */
335 
336  if ( iEdge != 1 ) {
337  OOFEM_ERROR("wrong edge number");
338  }
339 
340 
341  answer.resize(4);
342  answer.at(1) = 1;
343  answer.at(2) = 2;
344  answer.at(3) = 3;
345  answer.at(4) = 4;
346 }
347 
348 double
350 {
351  if ( iEdge != 1 ) { // edge between nodes 1 2
352  OOFEM_ERROR("wrong egde number");
353  }
354 
355  double weight = gp->giveWeight();
356  return 0.5 * this->computeLength() * weight;
357 }
358 
359 int
361 {
362  // returns transformation matrix from
363  // edge local coordinate system
364  // to element local c.s
365  // (same as global c.s in this case)
366  //
367  // i.e. f(element local) = T * f(edge local)
368  //
369  //double dx,dy, length ;
370  double sine, cosine;
371 
372  answer.resize(2, 2);
373  answer.zero();
374 
375  sine = sin( this->givePitch() );
376  cosine = cos(pitch);
377 
378  answer.at(1, 1) = cosine;
379  answer.at(1, 2) = -sine;
380  answer.at(2, 1) = sine;
381  answer.at(2, 2) = cosine;
382 
383  return 1;
384 }
385 
386 
388 {
389  // return interpolator
390  return & interp[cs_mode];
391 }
392 
393 
394 #ifdef __OOFEG
396 {
397  int c1, c2;
398  resolveCoordIndices(c1, c2);
399 
400  GraphicObj *go;
401  // if (!go) { // create new one
402  WCRec p [ 2 ]; /* point */
403  if ( !gc.testElementGraphicActivity(this) ) {
404  return;
405  }
406 
407  EASValsSetLineWidth(OOFEG_RAW_GEOMETRY_WIDTH);
408  EASValsSetColor( gc.getElementColor() );
409  EASValsSetLayer(OOFEG_RAW_GEOMETRY_LAYER);
410  if ( cs_mode == 0 ) {
411  p [ 0 ].x = ( FPNum ) this->giveNode(1)->giveCoordinate(c1);
412  p [ 0 ].y = 0.;
413  p [ 0 ].z = ( FPNum ) this->giveNode(1)->giveCoordinate(c2);
414  p [ 1 ].x = ( FPNum ) this->giveNode(2)->giveCoordinate(c1);
415  p [ 1 ].y = 0.;
416  p [ 1 ].z = ( FPNum ) this->giveNode(2)->giveCoordinate(c2);
417  } else if ( cs_mode == 1 ) {
418  p [ 0 ].x = ( FPNum ) this->giveNode(1)->giveCoordinate(c1);
419  p [ 0 ].y = ( FPNum ) this->giveNode(1)->giveCoordinate(c2);
420  p [ 0 ].z = 0.;
421  p [ 1 ].x = ( FPNum ) this->giveNode(2)->giveCoordinate(c1);
422  p [ 1 ].y = ( FPNum ) this->giveNode(2)->giveCoordinate(c2);
423  p [ 1 ].z = 0.;
424  } else if ( cs_mode == 2 ) {
425  p [ 0 ].x = 0.;
426  p [ 0 ].y = ( FPNum ) this->giveNode(1)->giveCoordinate(c1);
427  p [ 0 ].z = ( FPNum ) this->giveNode(1)->giveCoordinate(c2);
428  p [ 1 ].x = 0.;
429  p [ 1 ].y = ( FPNum ) this->giveNode(2)->giveCoordinate(c1);
430  p [ 1 ].z = ( FPNum ) this->giveNode(2)->giveCoordinate(c2);
431  }
432 
433  go = CreateLine3D(p);
434  EGWithMaskChangeAttributes(WIDTH_MASK | COLOR_MASK | LAYER_MASK, go);
435  EGAttachObject(go, ( EObjectP ) this);
436  EMAddGraphicsToModel(ESIModel(), go);
437 }
438 
439 
441 {
442  int c1, c2;
443  resolveCoordIndices(c1, c2);
444 
445  GraphicObj *go;
446  double defScale = gc.getDefScale();
447  // if (!go) { // create new one
448  WCRec p [ 2 ]; /* point */
449  if ( !gc.testElementGraphicActivity(this) ) {
450  return;
451  }
452 
453  EASValsSetLineWidth(OOFEG_DEFORMED_GEOMETRY_WIDTH);
454  EASValsSetColor( gc.getDeformedElementColor() );
455  EASValsSetLayer(OOFEG_DEFORMED_GEOMETRY_LAYER);
456 
457  p [ 0 ].x = ( FPNum ) this->giveNode(1)->giveUpdatedCoordinate(1, tStep, defScale);
458  p [ 0 ].y = ( FPNum ) this->giveNode(1)->giveUpdatedCoordinate(2, tStep, defScale);
459  p [ 0 ].z = ( FPNum ) this->giveNode(1)->giveUpdatedCoordinate(3, tStep, defScale);
460 
461  p [ 1 ].x = ( FPNum ) this->giveNode(2)->giveUpdatedCoordinate(1, tStep, defScale);
462  p [ 1 ].y = ( FPNum ) this->giveNode(2)->giveUpdatedCoordinate(2, tStep, defScale);
463  p [ 1 ].z = ( FPNum ) this->giveNode(2)->giveUpdatedCoordinate(3, tStep, defScale);
464 
465  go = CreateLine3D(p);
466  EGWithMaskChangeAttributes(WIDTH_MASK | COLOR_MASK | LAYER_MASK, go);
467  EMAddGraphicsToModel(ESIModel(), go);
468 }
469 #endif
470 } // end namespace oofem
CrossSection * giveCrossSection()
Definition: element.C:495
int testElementGraphicActivity(Element *)
Test if particular element passed fulfills various filtering criteria for its graphics output...
virtual bool isActivated(TimeStep *tStep)
Definition: element.C:793
Truss2d(int n, Domain *d)
Definition: truss2d.C:59
double givePitch()
Definition: truss2d.C:216
Class and object Domain.
Definition: domain.h:115
virtual FEInterpolation * giveInterpolation() const
Definition: truss2d.C:387
Abstract base class for "structural" finite elements with geometrical nonlinearities.
virtual void computeBmatrixAt(GaussPoint *gp, FloatMatrix &, int=1, int=ALL_STRAINS)
Computes the geometrical matrix of receiver in given integration point.
Definition: truss2d.C:70
virtual void computeConstitutiveMatrixAt(FloatMatrix &answer, MatResponseMode rMode, GaussPoint *gp, TimeStep *tStep)
Computes constitutive matrix of receiver.
Definition: truss2d.C:322
double & at(int i)
Coefficient access function.
Definition: floatarray.h:131
#define OOFEG_RAW_GEOMETRY_LAYER
oofem::oofegGraphicContext gc[OOFEG_LAST_LAYER]
virtual IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description stored in input record.
Definition: truss2d.C:285
virtual int computeLoadLEToLRotationMatrix(FloatMatrix &, int, GaussPoint *gp)
Returns transformation matrix from local edge c.s to element local coordinate system of load vector c...
Definition: truss2d.C:360
void resolveCoordIndices(int &c1, int &c2)
Definition: truss2d.C:265
virtual double computeEdgeVolumeAround(GaussPoint *gp, int)
Computes volume related to integration point on local edge.
Definition: truss2d.C:349
virtual double giveCoordinate(int i)
Definition: node.C:82
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.
#define OOFEG_DEFORMED_GEOMETRY_LAYER
static FEI2dLineLin interp[3]
Definition: truss2d.h:67
Class representing a general abstraction for finite element interpolation class.
Definition: feinterpol.h:132
virtual double computeLength()
Computes the length (zero for all but 1D geometries)
Definition: truss2d.C:194
virtual void drawRawGeometry(oofegGraphicContext &gc, TimeStep *tStep)
Definition: truss2d.C:395
#define OOFEM_ERROR(...)
Definition: error.h:61
REGISTER_Element(LSpace)
#define OOFEG_RAW_GEOMETRY_WIDTH
virtual double giveWeight()
Returns integration weight of receiver.
Definition: gausspoint.h:181
UnknownType
Type representing particular unknown (its physical meaning).
Definition: unknowntype.h:55
StructuralCrossSection * giveStructuralCrossSection()
Helper function which returns the structural cross-section for the element.
void times(double f)
Multiplies receiver by factor f.
Definition: floatmatrix.C:1594
virtual void computeBHmatrixAt(GaussPoint *gp, FloatMatrix &)
Computes a matrix which, multiplied by the column matrix of nodal displacements, gives the displaceme...
Definition: truss2d.C:103
virtual double giveUpdatedCoordinate(int ic, TimeStep *tStep, double scale=1.)
Returns updated ic-th coordinate of receiver.
Definition: node.C:245
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 void giveEdgeDofMapping(IntArray &answer, int) const
Assembles edge dof mapping mask, which provides mapping between edge local DOFs and "global" element ...
Definition: truss2d.C:329
virtual int setupIntegrationPoints(IntegrationRule &irule, int npoints, Element *element)
Sets up integration rule for the given element.
Definition: crosssection.C:54
double length
Definition: truss2d.h:63
virtual void computeGaussPoints()
Initializes the array of integration rules member variable.
Definition: truss2d.C:111
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 double give(CrossSectionProperty a, GaussPoint *gp)
Returns the value of cross section property at given point.
Definition: crosssection.C:151
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
virtual void giveRealStress_1d(FloatArray &answer, GaussPoint *gp, const FloatArray &reducedStrain, TimeStep *tStep)=0
void zero()
Zeroes all coefficients of receiver.
Definition: floatarray.C:658
virtual double computeVolumeAround(GaussPoint *gp)
Returns volume related to given integration point.
Definition: truss2d.C:185
#define OOFEG_DEFORMED_GEOMETRY_WIDTH
#define _IFT_Truss2d_cs
Definition: truss2d.h:43
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 computeGlobalCoordinates(FloatArray &answer, const FloatArray &lcoords)
Computes the global coordinates from given element&#39;s local coordinates.
Definition: truss2d.C:164
void zero()
Zeroes all coefficient of receiver.
Definition: floatmatrix.C:1326
virtual void giveStiffnessMatrix_1d(FloatMatrix &answer, MatResponseMode mode, GaussPoint *gp, TimeStep *tStep)=0
#define IR_GIVE_OPTIONAL_FIELD(__ir, __value, __id)
Macro facilitating the use of input record reading methods.
Definition: inputrecord.h:78
virtual void giveDofManDofIDMask(int inode, IntArray &answer) const
Returns dofmanager dof mask for node.
Definition: truss2d.C:302
double pitch
Definition: truss2d.h:64
the oofem namespace is to define a context or scope in which all oofem names are defined.
virtual void drawDeformedGeometry(oofegGraphicContext &gc, TimeStep *tStep, UnknownType)
Definition: truss2d.C:440
Class implementing node in finite element mesh.
Definition: node.h:87
Node * giveNode(int i) const
Returns reference to the i-th node of element.
Definition: element.h:610
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.
Definition: truss2d.C:315
Class representing integration point in finite element program.
Definition: gausspoint.h:93
#define OOFEM_WARNING(...)
Definition: error.h:62
virtual void computeNmatrixAt(const FloatArray &iLocCoord, FloatMatrix &)
Computes interpolation matrix for element unknowns.
Definition: truss2d.C:144
Class representing solution step.
Definition: timestep.h:80
virtual void computeLumpedMassMatrix(FloatMatrix &answer, TimeStep *tStep)
Computes lumped mass matrix of receiver.
Definition: truss2d.C:123
int numberOfDofMans
Number of dofmanagers.
Definition: element.h:149
Class representing Gaussian-quadrature integration rule.
void resize(int s)
Resizes receiver towards requested size.
Definition: floatarray.C:631
virtual int giveLocalCoordinateSystem(FloatMatrix &answer)
Returns local coordinate system of receiver Required by material models with ortho- and anisotrophy...
Definition: truss2d.C:240

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:32 for OOFEM by doxygen 1.8.11 written by Dimitri van Heesch, © 1997-2011