OOFEM  2.4 OOFEM.org - Object Oriented Finite Element Solver
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
36 #include "mathfem.h"
37 #include "floatmatrix.h"
38 #include "floatarray.h"
39 #include "gaussintegrationrule.h"
40
41 namespace oofem {
42 void
44 {
45  double l1 = lcoords.at(1);
46  double l2 = lcoords.at(2);
47  double l3 = 1. - l1 - l2;
48
50  ( 2. * l1 - 1. ) * l1,
51  ( 2. * l2 - 1. ) * l2,
52  ( 2. * l3 - 1. ) * l3,
53  4. * l1 * l2,
54  4. * l2 * l3,
55  4. * l3 * l1
56  };
57 }
58
59 double
61 {
62  FloatMatrix jacobianMatrix(2, 2), inv, dn;
63
64  this->evaldNdxi(dn, lcoords, cellgeo);
65  for ( int i = 1; i <= dn.giveNumberOfRows(); i++ ) {
66  double x = cellgeo.giveVertexCoordinates(i)->at(xind);
67  double y = cellgeo.giveVertexCoordinates(i)->at(yind);
68
69  jacobianMatrix.at(1, 1) += dn.at(i, 1) * x;
70  jacobianMatrix.at(1, 2) += dn.at(i, 1) * y;
71  jacobianMatrix.at(2, 1) += dn.at(i, 2) * x;
72  jacobianMatrix.at(2, 2) += dn.at(i, 2) * y;
73  }
74  inv.beInverseOf(jacobianMatrix);
75
77  return jacobianMatrix.giveDeterminant();
78 }
79
80 void
82 {
83  double x1, x2, x3, y1, y2, y3, y23, x32, y31, x13, area;
84
86
87  x1 = cellgeo.giveVertexCoordinates(1)->at(xind);
88  x2 = cellgeo.giveVertexCoordinates(2)->at(xind);
89  x3 = cellgeo.giveVertexCoordinates(3)->at(xind);
90
91  y1 = cellgeo.giveVertexCoordinates(1)->at(yind);
92  y2 = cellgeo.giveVertexCoordinates(2)->at(yind);
93  y3 = cellgeo.giveVertexCoordinates(3)->at(yind);
94
95  area = 0.5 * ( x2 * y3 + x1 * y2 + y1 * x3 - x2 * y1 - x3 * y2 - x1 * y3 );
96
97  y23 = ( y2 - y3 ) / ( 2. * area );
98  x32 = ( x3 - x2 ) / ( 2. * area );
99
100  y31 = ( y3 - y1 ) / ( 2. * area );
101  x13 = ( x1 - x3 ) / ( 2. * area );
102
103  answer.at(1, 1) = 4 * y23 * y23;
104  answer.at(1, 2) = 4 * x32 * x32;
105  answer.at(1, 3) = 4 * y23 * x32;
106
107  answer.at(2, 1) = 4 * y31 * y31;
108  answer.at(2, 2) = 4 * x13 * x13;
109  answer.at(2, 3) = 4 * y31 * x13;
110
111  answer.at(3, 1) = 4 * y23 * y23 + 8 * y31 * y23 + 4 * y31 * y31;
112  answer.at(3, 2) = 4 * x32 * x32 + 8 * x13 * x32 + 4 * x13 * x13;
113  answer.at(3, 3) = 4 * y23 * x32 + 4 * y31 * x32 + 4 * y23 * x13 + 4 * y31 * x13;
114
115  answer.at(4, 1) = 8 * y31 * y23;
116  answer.at(4, 2) = 8 * x13 * x32;
117  answer.at(4, 3) = 4 * y31 * x32 + 4 * y23 * x13;
118
119  answer.at(5, 1) = ( -8 ) * y31 * y23 + ( -8 ) * y31 * y31;
120  answer.at(5, 2) = ( -8 ) * x13 * x32 + ( -8 ) * x13 * x13;
121  answer.at(5, 3) = ( -4 ) * y31 * x32 + ( -4 ) * y23 * x13 + ( -8 ) * y31 * x13;
122
123  answer.at(6, 1) = ( -8 ) * y23 * y23 + ( -8 ) * y31 * y23;
124  answer.at(6, 2) = ( -8 ) * x32 * x32 + ( -8 ) * x13 * x32;
125  answer.at(6, 3) = ( -8 ) * y23 * x32 + ( -4 ) * y31 * x32 + ( -4 ) * y23 * x13;
126 }
127
128
129
130 void
132 {
133  FloatArray n;
134  this->evalN(n, lcoords, cellgeo);
135
138  for ( int i = 1; i <= 6; i++ ) {
139  answer.at(1) += n.at(i) * cellgeo.giveVertexCoordinates(i)->at(xind);
140  answer.at(2) += n.at(i) * cellgeo.giveVertexCoordinates(i)->at(yind);
141  }
142 }
143
144
145 void
146 FEI2dTrQuad :: edgeEvalN(FloatArray &answer, int iedge, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
147 {
148  double n3, ksi = lcoords.at(1);
149  n3 = 1. - ksi * ksi;
150
151  answer = { ( 1. - ksi - n3 ) * 0.5, ( 1. + ksi - n3 ) * 0.5, n3 };
152 }
153
154 void
156  const FloatArray &lcoords, const FEICellGeometry &cellgeo)
157 {
158  // I think it at least should return both dNds and J. Both are almost always needed.
159  // In fact, dxdxi is also needed sometimes (surface tension)
160 #if 0
161  IntArray edgeNodes;
162  FloatArray dNdxi(3);
163  FloatArray dxdxi(2);
164  double xi = lcoords.at(1);
165  this->computeLocalEdgeMapping(edgeNodes, iedge);
166  dNdxi.at(1) = xi - 0.5;
167  dNdxi.at(2) = xi + 0.5;
168  dNdxi.at(3) = -2 * xi;
169
170  dxdxi.at(1) = dNdxi.at(1) * cellgeo.giveVertexCoordinates( edgeNodes.at(1) )->at(xind) +
171  dNdxi.at(2) * cellgeo.giveVertexCoordinates( edgeNodes.at(2) )->at(xind) +
172  dNdxi.at(3) * cellgeo.giveVertexCoordinates( edgeNodes.at(3) )->at(xind);
173  dxdxi.at(2) = dNdxi.at(1) * cellgeo.giveVertexCoordinates( edgeNodes.at(1) )->at(yind) +
174  dNdxi.at(2) * cellgeo.giveVertexCoordinates( edgeNodes.at(2) )->at(yind) +
175  dNdxi.at(3) * cellgeo.giveVertexCoordinates( edgeNodes.at(3) )->at(yind);
176
177  double J = dxdxi.computeNorm();
180  return J;
181
182 #endif
183  double xi = lcoords.at(1);
184  double J = edgeGiveTransformationJacobian(iedge, lcoords, cellgeo);
186  ( xi - 0.5 ) / J,
187  ( xi + 0.5 ) / J,
188  -2 * xi / J
189  };
190 }
191
192 double FEI2dTrQuad :: edgeEvalNormal(FloatArray &normal, int iedge, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
193 {
194  IntArray edgeNodes;
195  this->computeLocalEdgeMapping(edgeNodes, iedge);
196  double xi = lcoords(0);
197  double dN1dxi = -0.5 + xi;
198  double dN2dxi = 0.5 + xi;
199  double dN3dxi = -2.0 * xi;
200
201  normal.resize(2);
202
203  normal.at(1) = dN1dxi * cellgeo.giveVertexCoordinates( edgeNodes.at(1) )->at(yind) +
204  dN2dxi *cellgeo.giveVertexCoordinates( edgeNodes.at(2) )->at(yind) +
205  dN3dxi *cellgeo.giveVertexCoordinates( edgeNodes.at(3) )->at(yind);
206
207  normal.at(2) = -dN1dxi *cellgeo.giveVertexCoordinates( edgeNodes.at(1) )->at(xind) +
208  - dN2dxi *cellgeo.giveVertexCoordinates( edgeNodes.at(2) )->at(xind) +
209  - dN3dxi *cellgeo.giveVertexCoordinates( edgeNodes.at(3) )->at(xind);
210
211  return normal.normalize();
212 }
213
214 void
216  const FloatArray &lcoords, const FEICellGeometry &cellgeo)
217 {
218  IntArray edgeNodes;
219  FloatArray n;
220  this->computeLocalEdgeMapping(edgeNodes, iedge);
221  this->edgeEvalN(n, iedge, lcoords, cellgeo);
222
224  answer.at(1) = ( n.at(1) * cellgeo.giveVertexCoordinates( edgeNodes.at(1) )->at(xind) +
225  n.at(2) * cellgeo.giveVertexCoordinates( edgeNodes.at(2) )->at(xind) +
226  n.at(3) * cellgeo.giveVertexCoordinates( edgeNodes.at(3) )->at(xind) );
227  answer.at(2) = ( n.at(1) * cellgeo.giveVertexCoordinates( edgeNodes.at(1) )->at(yind) +
228  n.at(2) * cellgeo.giveVertexCoordinates( edgeNodes.at(2) )->at(yind) +
229  n.at(3) * cellgeo.giveVertexCoordinates( edgeNodes.at(3) )->at(yind) );
230 }
231
232
233 void
235 {
236  int aNode = 0, bNode = 0, cNode = 0;
237  edgeNodes.resize(3);
238
239  if ( iedge == 1 ) { // edge between nodes 1 2
240  aNode = 1;
241  bNode = 2;
242  cNode = 4;
243  } else if ( iedge == 2 ) { // edge between nodes 2 3
244  aNode = 2;
245  bNode = 3;
246  cNode = 5;
247  } else if ( iedge == 3 ) { // edge between nodes 2 3
248  aNode = 3;
249  bNode = 1;
250  cNode = 6;
251  } else {
252  OOFEM_ERROR("wrong egde number (%d)", iedge);
253  }
254
255  edgeNodes.at(1) = aNode;
256  edgeNodes.at(2) = bNode;
257  edgeNodes.at(3) = cNode;
258 }
259
260
261
262 void FEI2dTrQuad :: evaldNdxi(FloatMatrix &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
263 {
264  double l1, l2, l3;
265  l1 = lcoords.at(1);
266  l2 = lcoords.at(2);
267  l3 = 1.0 - l1 - l2;
268
270
271  answer.at(1, 1) = 4.0 * l1 - 1.0;
273  answer.at(3, 1) = -1.0 * ( 4.0 * l3 - 1.0 );
274  answer.at(4, 1) = 4.0 * l2;
275  answer.at(5, 1) = -4.0 * l2;
276  answer.at(6, 1) = 4.0 * l3 - 4.0 * l1;
277
279  answer.at(2, 2) = 4.0 * l2 - 1.0;
280  answer.at(3, 2) = -1.0 * ( 4.0 * l3 - 1.0 );
281  answer.at(4, 2) = 4.0 * l1;
282  answer.at(5, 2) = 4.0 * l3 - 4.0 * l2;
283  answer.at(6, 2) = -4.0 * l1;
284 }
285
286
287 double
289 {
290  const FloatArray *p;
291  double x1, x2, x3, x4, x5, x6, y1, y2, y3, y4, y5, y6;
292
293  p = cellgeo.giveVertexCoordinates(1);
294  x1 = p->at(1);
295  y1 = p->at(2);
296  p = cellgeo.giveVertexCoordinates(2);
297  x2 = p->at(1);
298  y2 = p->at(2);
299  p = cellgeo.giveVertexCoordinates(3);
300  x3 = p->at(1);
301  y3 = p->at(2);
302  p = cellgeo.giveVertexCoordinates(4);
303  x4 = p->at(1);
304  y4 = p->at(2);
305  p = cellgeo.giveVertexCoordinates(5);
306  x5 = p->at(1);
307  y5 = p->at(2);
308  p = cellgeo.giveVertexCoordinates(6);
309  x6 = p->at(1);
310  y6 = p->at(2);
311
312  return fabs( ( 4 * ( -( x4 * y1 ) + x6 * y1 + x4 * y2 - x5 * y2 + x5 * y3 - x6 * y3 ) + x2 * ( y1 - y3 - 4 * y4 + 4 * y5 ) +
313  x1 * ( -y2 + y3 + 4 * y4 - 4 * y6 ) + x3 * ( -y1 + y2 - 4 * y5 + 4 * y6 ) ) / 6 );
314 }
315
316 bool FEI2dTrQuad :: inside(const FloatArray &lcoords) const
317 {
318  const double point_tol = 1.0e-3;
319  bool inside = true;
320  for ( int i = 1; i <= 2; i++ ) {
321  if ( lcoords.at(i) < - point_tol ) {
322  inside = false;
323  } else if ( lcoords.at(i) > ( 1. + point_tol ) ) {
324  inside = false;
325  }
326  }
327
328  if ( 1. - lcoords.at(1) - lcoords.at(2) < - point_tol ) {
329  inside = false;
330  } else if ( 1. - lcoords.at(1) - lcoords.at(2) > ( 1. + point_tol ) ) {
331  inside = false;
332  }
333
334  return inside;
335 }
336
337 double
339 {
340  IntArray eNodes;
341  const FloatArray *node;
342  double x1, x2, x3, y1, y2, y3;
343
344  this->computeLocalEdgeMapping(eNodes, iEdge);
345
346  node = cellgeo.giveVertexCoordinates( eNodes.at(1) );
347  x1 = node->at(xind);
348  y1 = node->at(yind);
349
350  node = cellgeo.giveVertexCoordinates( eNodes.at(2) );
351  x2 = node->at(xind);
352  y2 = node->at(yind);
353
354  node = cellgeo.giveVertexCoordinates( eNodes.at(3) );
355  x3 = node->at(xind);
356  y3 = node->at(yind);
357
358  return -( x1 * y2 - x2 * y1 + 4 * ( x3 * ( y1 - y2 ) + y3 * ( x2 - x1 ) ) ) / 3.0;
359 }
360
363 {
364  IntegrationRule *iRule = new GaussIntegrationRule(1, NULL);
365  int points = iRule->getRequiredNumberOfIntegrationPoints(_Triangle, order + 2);
366  iRule->SetUpPointsOnTriangle(points, _Unknown);
367  return iRule;
368 }
369 } // end namespace oofem
virtual void computeLocalEdgeMapping(IntArray &edgeNodes, int iedge)
double & at(int i)
Coefficient access function.
Definition: floatarray.h:131
virtual const FloatArray * giveVertexCoordinates(int i) const =0
Class representing a general abstraction for cell geometry.
Definition: feinterpol.h:62
virtual double giveArea(const FEICellGeometry &cellgeo) const
Computes the exact area.
virtual double edgeGiveTransformationJacobian(int iedge, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates the edge Jacobian of transformation between local and global coordinates.
Definition: feinterpol2d.C:175
virtual void evald2Ndx2(FloatMatrix &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates the matrix of second derivatives of interpolation functions (shape functions) at given poin...
virtual void edgeEvaldNds(FloatArray &answer, int iedge, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates the matrix of derivatives of edge interpolation functions (shape functions) at given point...
Class implementing an array of integers.
Definition: intarray.h:61
int & at(int i)
Coefficient access function.
Definition: intarray.h:103
virtual double evalNXIntegral(int iEdge, const FEICellGeometry &cellgeo)
Computes the integral .
Abstract base class representing integration rule.
virtual void local2global(FloatArray &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates global coordinates from given local ones.
virtual void evalN(FloatArray &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates the array of interpolation functions (shape functions) at given point.
virtual IntegrationRule * giveIntegrationRule(int order)
Sets up a suitable integration rule for numerical integrating over volume.
#define OOFEM_ERROR(...)
Definition: error.h:61
virtual bool inside(const FloatArray &lcoords) const
virtual double evaldNdx(FloatMatrix &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates the matrix of derivatives of interpolation functions (shape functions) at given point...
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
Class representing vector of real numbers.
Definition: floatarray.h:82
Implementation of matrix containing floating point numbers.
Definition: floatmatrix.h:94
virtual void edgeLocal2global(FloatArray &answer, int iedge, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates edge global coordinates from given local ones.
double computeNorm() const
Computes the norm (or length) of the vector.
Definition: floatarray.C:840
void resize(int rows, int cols)
Checks size of receiver towards requested bounds.
Definition: floatmatrix.C:1358
virtual int getRequiredNumberOfIntegrationPoints(integrationDomain dType, int approxOrder)
Abstract service.
void zero()
Definition: floatarray.C:658
void times(double s)
Definition: floatarray.C:818
virtual int SetUpPointsOnTriangle(int, MaterialMode mode)
Sets up receiver&#39;s integration points on triangular (area coords) integration domain.
void beProductTOf(const FloatMatrix &a, const FloatMatrix &b)
Assigns to the receiver product of .
Definition: floatmatrix.C:397
virtual double edgeEvalNormal(FloatArray &normal, int iedge, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates the normal on the given edge.
the oofem namespace is to define a context or scope in which all oofem names are defined.
void beInverseOf(const FloatMatrix &src)
Modifies receiver to become inverse of given parameter.
Definition: floatmatrix.C:835
double normalize()
Definition: floatarray.C:828
int giveNumberOfRows() const
Returns number of rows of receiver.
Definition: floatmatrix.h:156
virtual void evaldNdxi(FloatMatrix &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates the matrix of derivatives of interpolation functions (shape functions) at given point...