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 double
44 {
45  double x1, x2, x3, x4, y1, y2, y3, y4;
46  double x85, x56, x67, x78, y85, y56, y67, y78;
47
48  const FloatArray *node1 = cellgeo.giveVertexCoordinates(1);
49  const FloatArray *node2 = cellgeo.giveVertexCoordinates(2);
50  const FloatArray *node3 = cellgeo.giveVertexCoordinates(3);
51  const FloatArray *node4 = cellgeo.giveVertexCoordinates(4);
52  const FloatArray *node5 = cellgeo.giveVertexCoordinates(5);
53  const FloatArray *node6 = cellgeo.giveVertexCoordinates(6);
54  const FloatArray *node7 = cellgeo.giveVertexCoordinates(7);
55  const FloatArray *node8 = cellgeo.giveVertexCoordinates(8);
56
57  x1 = node1->at(xind);
58  x2 = node2->at(xind);
59  x3 = node3->at(xind);
60  x4 = node4->at(xind);
61
62  y1 = node1->at(yind);
63  y2 = node2->at(yind);
64  y3 = node3->at(yind);
65  y4 = node4->at(yind);
66
67  x85 = node8->at(xind) - node5->at(xind);
68  x56 = node5->at(xind) - node6->at(xind);
69  x67 = node6->at(xind) - node7->at(xind);
70  x78 = node7->at(xind) - node8->at(xind);
71
72  y85 = node8->at(yind) - node5->at(yind);
73  y56 = node5->at(yind) - node6->at(yind);
74  y67 = node6->at(yind) - node7->at(yind);
75  y78 = node7->at(yind) - node8->at(yind);
76
77  double p1 = ( x2 - x4 ) * ( y1 - y3 ) - ( x1 - x3 ) * ( y2 - y4 );
78  double p2 = y1 * x85 + y2 * x56 + y3 * x67 + y4 * x78 - x1 * y85 - x2 * y56 - x3 * y67 - x4 * y78;
79
80  return fabs(p1 + p2 * 4.0) / 6.; // Expression derived with mathematica, but not verified in any computations
81 }
82
83 void
85 {
86  double ksi, eta;
87
88  ksi = lcoords.at(1);
89  eta = lcoords.at(2);
90
92  ( 1. + ksi ) * ( 1. + eta ) * 0.25 * ( ksi + eta - 1. ),
93  ( 1. - ksi ) * ( 1. + eta ) * 0.25 * ( -ksi + eta - 1. ),
94  ( 1. - ksi ) * ( 1. - eta ) * 0.25 * ( -ksi - eta - 1. ),
95  ( 1. + ksi ) * ( 1. - eta ) * 0.25 * ( ksi - eta - 1. ),
96  0.5 * ( 1. - ksi * ksi ) * ( 1. + eta ),
97  0.5 * ( 1. - ksi ) * ( 1. - eta * eta ),
98  0.5 * ( 1. - ksi * ksi ) * ( 1. - eta ),
99  0.5 * ( 1. + ksi ) * ( 1. - eta * eta )
100  };
101 }
102
103 double
105 {
106  FloatMatrix jacobianMatrix(2, 2), inv, dn;
107
108  this->evaldNdxi(dn, lcoords, cellgeo);
109  for ( int i = 1; i <= dn.giveNumberOfRows(); i++ ) {
110  double x = cellgeo.giveVertexCoordinates(i)->at(xind);
111  double y = cellgeo.giveVertexCoordinates(i)->at(yind);
112
113  jacobianMatrix.at(1, 1) += dn.at(i, 1) * x;
114  jacobianMatrix.at(1, 2) += dn.at(i, 1) * y;
115  jacobianMatrix.at(2, 1) += dn.at(i, 2) * x;
116  jacobianMatrix.at(2, 2) += dn.at(i, 2) * y;
117  }
118  inv.beInverseOf(jacobianMatrix);
119
121  return jacobianMatrix.giveDeterminant();
122 }
123
124 void
126 {
127  FloatArray n;
128
129  this->evalN(n, lcoords, cellgeo);
130
133  for ( int i = 1; i <= n.giveSize(); i++ ) {
134  answer.at(1) += n.at(i) * cellgeo.giveVertexCoordinates(i)->at(xind);
135  answer.at(2) += n.at(i) * cellgeo.giveVertexCoordinates(i)->at(yind);
136  }
137 }
138
140 {
141  const FloatArray *n1 = cellgeo.giveVertexCoordinates(1);
142  const FloatArray *n2 = cellgeo.giveVertexCoordinates(3);
143  return n1->distance(n2);
144 }
145
147 {
148  const double point_tol = 1.0e-3;
149  bool inside = true;
150  for ( int i = 1; i <= 2; i++ ) {
151  if ( lcoords.at(i) < ( -1. - point_tol ) ) {
152  inside = false;
153  } else if ( lcoords.at(i) > ( 1. + point_tol ) ) {
154  inside = false;
155  }
156  }
157
158  return inside;
159 }
160
161
162 void
164 {
165  // 1-------3-------2
166
167  double n3, ksi = lcoords.at(1);
168  n3 = 1. - ksi * ksi;
169
170  answer = { ( 1. - ksi - n3 ) * 0.5, ( 1. + ksi - n3 ) * 0.5, n3 };
171 }
172
173 void
175  const FloatArray &lcoords, const FEICellGeometry &cellgeo)
176 {
177  double ksi = lcoords.at(1);
178  answer = { ksi - 0.5, ksi + 0.5, ksi * 2.0 };
179 }
180
181 void
183  const FloatArray &lcoords, const FEICellGeometry &cellgeo)
184 {
185  IntArray edgeNodes;
186  FloatArray n;
187  this->computeLocalEdgeMapping(edgeNodes, iedge);
188  this->edgeEvalN(n, iedge, lcoords, cellgeo);
189
191  answer.at(1) = ( n.at(1) * cellgeo.giveVertexCoordinates( edgeNodes.at(1) )->at(xind) +
192  n.at(2) * cellgeo.giveVertexCoordinates( edgeNodes.at(2) )->at(xind) +
193  n.at(3) * cellgeo.giveVertexCoordinates( edgeNodes.at(3) )->at(xind) );
194  answer.at(2) = ( n.at(1) * cellgeo.giveVertexCoordinates( edgeNodes.at(1) )->at(yind) +
195  n.at(2) * cellgeo.giveVertexCoordinates( edgeNodes.at(2) )->at(yind) +
196  n.at(3) * cellgeo.giveVertexCoordinates( edgeNodes.at(3) )->at(yind) );
197 }
198
199
200 void
202 {
203  int aNode = 0, bNode = 0, cNode = 0;
204  edgeNodes.resize(3);
205
206  if ( iedge == 1 ) { // edge between nodes 1 2
207  aNode = 1;
208  bNode = 2;
209  cNode = 5;
210  } else if ( iedge == 2 ) { // edge between nodes 2 3
211  aNode = 2;
212  bNode = 3;
213  cNode = 6;
214  } else if ( iedge == 3 ) { // edge between nodes 2 3
215  aNode = 3;
216  bNode = 4;
217  cNode = 7;
218  } else if ( iedge == 4 ) { // edge between nodes 3 4
219  aNode = 4;
220  bNode = 1;
221  cNode = 8;
222  } else {
223  OOFEM_ERROR("wrong edge number (%d)", iedge);
224  }
225
226  edgeNodes.at(1) = aNode;
227  edgeNodes.at(2) = bNode;
228  edgeNodes.at(3) = cNode;
229 }
230
231 double FEI2dQuadQuad :: edgeEvalNormal(FloatArray &normal, int iedge, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
232 {
233  IntArray edgeNodes;
234  this->computeLocalEdgeMapping(edgeNodes, iedge);
235  double xi = lcoords(0);
236  double dN1dxi = -0.5 + xi;
237  double dN2dxi = 0.5 + xi;
238  double dN3dxi = -2.0 * xi;
239
240  normal.resize(2);
241
242  normal.at(1) = dN1dxi * cellgeo.giveVertexCoordinates( edgeNodes.at(1) )->at(yind) +
243  dN2dxi *cellgeo.giveVertexCoordinates( edgeNodes.at(2) )->at(yind) +
244  dN3dxi *cellgeo.giveVertexCoordinates( edgeNodes.at(3) )->at(yind);
245
246  normal.at(2) = -dN1dxi *cellgeo.giveVertexCoordinates( edgeNodes.at(1) )->at(xind) +
247  - dN2dxi *cellgeo.giveVertexCoordinates( edgeNodes.at(2) )->at(xind) +
248  - dN3dxi *cellgeo.giveVertexCoordinates( edgeNodes.at(3) )->at(xind);
249
250  return normal.normalize();
251 }
252
253
255 {
256  double ksi, eta;
257  ksi = lcoords.at(1);
258  eta = lcoords.at(2);
260
261  // dn/dxi
262  answer.at(1, 1) = 0.25 * ( 1. + eta ) * ( 2.0 * ksi + eta );
263  answer.at(2, 1) = -0.25 * ( 1. + eta ) * ( -2.0 * ksi + eta );
264  answer.at(3, 1) = -0.25 * ( 1. - eta ) * ( -2.0 * ksi - eta );
265  answer.at(4, 1) = 0.25 * ( 1. - eta ) * ( 2.0 * ksi - eta );
266  answer.at(5, 1) = -ksi * ( 1. + eta );
267  answer.at(6, 1) = -0.5 * ( 1. - eta * eta );
268  answer.at(7, 1) = -ksi * ( 1. - eta );
269  answer.at(8, 1) = 0.5 * ( 1. - eta * eta );
270
271  answer.at(1, 2) = 0.25 * ( 1. + ksi ) * ( 2.0 * eta + ksi );
272  answer.at(2, 2) = 0.25 * ( 1. - ksi ) * ( 2.0 * eta - ksi );
273  answer.at(3, 2) = -0.25 * ( 1. - ksi ) * ( -2.0 * eta - ksi );
274  answer.at(4, 2) = -0.25 * ( 1. + ksi ) * ( -2.0 * eta + ksi );
275  answer.at(5, 2) = 0.5 * ( 1. - ksi * ksi );
276  answer.at(6, 2) = -eta * ( 1. - ksi );
277  answer.at(7, 2) = -0.5 * ( 1. - ksi * ksi );
278  answer.at(8, 2) = -eta * ( 1. + ksi );
279 }
280
281
282 double
284 {
285  IntArray eNodes;
286  const FloatArray *node;
287  double x1, x2, x3, y1, y2, y3;
288
289  this->computeLocalEdgeMapping(eNodes, iEdge);
290
291  node = cellgeo.giveVertexCoordinates( eNodes.at(1) );
292  x1 = node->at(xind);
293  y1 = node->at(yind);
294
295  node = cellgeo.giveVertexCoordinates( eNodes.at(2) );
296  x2 = node->at(xind);
297  y2 = node->at(yind);
298
299  node = cellgeo.giveVertexCoordinates( eNodes.at(3) );
300  x3 = node->at(xind);
301  y3 = node->at(yind);
302
303  return -( x1 * y2 - x2 * y1 + 4 * ( x3 * ( y1 - y2 ) + y3 * ( x2 - x1 ) ) ) / 3.0;
304 }
305
306
309 {
310  IntegrationRule *iRule = new GaussIntegrationRule(1, NULL);
311  int points = iRule->getRequiredNumberOfIntegrationPoints(_Square, order + 4);
312  iRule->SetUpPointsOnSquare(points, _Unknown);
313  return iRule;
314 }
315
316
317 /*
319  */
320
321 double
323 {
324  FloatArray N;
325  this->evalN( N, lcoords, cellgeo);
326
327  double r = 0.0;
328  for ( int i = 1; i <= 8; i++ ) {
329  double x = cellgeo.giveVertexCoordinates(i)->at(1);
330  r += x * N.at(i);
331  }
332
334 }
335
336 double
338  const FEICellGeometry &cellgeo)
339 {
340  IntArray edgeNodes;
341  FloatArray n;
342  this->computeLocalEdgeMapping(edgeNodes, iedge);
343  this->edgeEvalN(n, iedge, lcoords, cellgeo);
344
345  double r = n.at(1)*cellgeo.giveVertexCoordinates(edgeNodes.at(1))->at(1) +
346  n.at(2)*cellgeo.giveVertexCoordinates(edgeNodes.at(2))->at(1) +
347  n.at(3)*cellgeo.giveVertexCoordinates(edgeNodes.at(3))->at(1);
349
350 }
351
352 double
354 {
355  return this->edgeGiveTransformationJacobian(boundary, lcoords, cellgeo);
356 }
357
358 double
360 {
361  return this->edgeGiveTransformationJacobian(boundary, lcoords, cellgeo);
362 }
363
364
365
366
367 } // end namespace oofem
virtual double boundaryEdgeGiveTransformationJacobian(int boundary, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates the determinant of the transformation Jacobian on the requested boundary.
virtual void edgeEvalN(FloatArray &answer, int iedge, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates the array of edge interpolation functions (shape functions) at given point.
virtual double giveTransformationJacobian(const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates the determinant of the transformation.
Definition: feinterpol.C:43
virtual double giveCharacteristicLength(const FEICellGeometry &cellgeo) const
Returns a characteristic length of the geometry, typically a diagonal or edge length.
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 void edgeLocal2global(FloatArray &answer, int iedge, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates edge global coordinates from given local ones.
virtual void evaldNdxi(FloatMatrix &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates the matrix of derivatives of interpolation functions (shape functions) at given point...
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 int SetUpPointsOnSquare(int, MaterialMode mode)
Sets up receiver&#39;s integration points on unit square integration domain.
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
Abstract base class representing integration rule.
virtual double giveTransformationJacobian(const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates the determinant of the transformation.
double distance(const FloatArray &x) const
Computes the distance between position represented by receiver and position given as parameter...
Definition: floatarray.C:489
#define OOFEM_ERROR(...)
Definition: error.h:61
virtual double edgeEvalNormal(FloatArray &normal, int iedge, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates the normal on the given edge.
virtual double boundaryGiveTransformationJacobian(int boundary, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates the determinant of the transformation Jacobian on the requested boundary.
virtual void computeLocalEdgeMapping(IntArray &edgeNodes, int iedge)
#define N(p, q)
Definition: mdm.C:367
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 IntegrationRule * giveIntegrationRule(int order)
Sets up a suitable integration rule for numerical integrating over volume.
Class representing vector of real numbers.
Definition: floatarray.h:82
Implementation of matrix containing floating point numbers.
Definition: floatmatrix.h:94
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.
virtual double evaldNdx(FloatMatrix &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates the matrix of derivatives of interpolation functions (shape functions) at given point...
virtual void evalN(FloatArray &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates the array of interpolation functions (shape functions) at given point.
void zero()
Definition: floatarray.C:658
virtual double giveArea(const FEICellGeometry &cellgeo) const
Computes the exact area.
void beProductTOf(const FloatMatrix &a, const FloatMatrix &b)
Assigns to the receiver product of .
Definition: floatmatrix.C:397
virtual bool inside(const FloatArray &lcoords) const
int giveSize() const
Definition: floatarray.h:218
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
virtual void local2global(FloatArray &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates global coordinates from given local ones.
int giveNumberOfRows() const
Returns number of rows of receiver.
Definition: floatmatrix.h:156
virtual double evalNXIntegral(int iEdge, const FEICellGeometry &cellgeo)
Computes the integral .
virtual double edgeGiveTransformationJacobian(int iedge, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates the edge Jacobian of transformation between local and global coordinates.