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  const FloatArray *node1 = cellgeo.giveVertexCoordinates(1);
46  const FloatArray *node2 = cellgeo.giveVertexCoordinates(2);
47  const FloatArray *node3 = cellgeo.giveVertexCoordinates(3);
48  const FloatArray *node4 = cellgeo.giveVertexCoordinates(4);
49
50  double x13 = node1->at(xind) - node3->at(xind);
51  double y13 = node1->at(yind) - node3->at(yind);
52  double x24 = node2->at(xind) - node4->at(xind);
53  double y24 = node2->at(yind) - node4->at(yind);
54
55  return fabs( 0.5 * ( x13 * y24 - x24 * y13 ) );
56 }
57
58 void
60 {
61  double ksi, eta;
62
63  ksi = lcoords.at(1);
64  eta = lcoords.at(2);
65
67  ( 1. + ksi ) * ( 1. + eta ) * 0.25,
68  ( 1. - ksi ) * ( 1. + eta ) * 0.25,
69  ( 1. - ksi ) * ( 1. - eta ) * 0.25,
70  ( 1. + ksi ) * ( 1. - eta ) * 0.25
71  };
72 }
73
74 double
76 {
77  FloatMatrix jacobianMatrix(2, 2), inv, dn;
78
79  this->evaldNdxi(dn, lcoords, cellgeo);
80  for ( int i = 1; i <= dn.giveNumberOfRows(); i++ ) {
81  double x = cellgeo.giveVertexCoordinates(i)->at(xind);
82  double y = cellgeo.giveVertexCoordinates(i)->at(yind);
83
84  jacobianMatrix.at(1, 1) += dn.at(i, 1) * x;
85  jacobianMatrix.at(1, 2) += dn.at(i, 1) * y;
86  jacobianMatrix.at(2, 1) += dn.at(i, 2) * x;
87  jacobianMatrix.at(2, 2) += dn.at(i, 2) * y;
88  }
89  inv.beInverseOf(jacobianMatrix);
90
92  return jacobianMatrix.giveDeterminant();
93 }
94
95 void
97 {
98  double ksi = lcoords.at(1);
99  double eta = lcoords.at(2);
100
101  double n1 = ( 1. + ksi ) * ( 1. + eta ) * 0.25;
102  double n2 = ( 1. - ksi ) * ( 1. + eta ) * 0.25;
103  double n3 = ( 1. - ksi ) * ( 1. - eta ) * 0.25;
104  double n4 = ( 1. + ksi ) * ( 1. - eta ) * 0.25;
105
106  const FloatArray* const p1 = cellgeo.giveVertexCoordinates(1);
107  const FloatArray* const p2 = cellgeo.giveVertexCoordinates(2);
108  const FloatArray* const p3 = cellgeo.giveVertexCoordinates(3);
109  const FloatArray* const p4 = cellgeo.giveVertexCoordinates(4);
110
111  answer = {n1 * p1->at(xind) + n2 * p2->at(xind) +
112  n3 * p3->at(xind) + n4 * p4->at(xind),
113  n1 * p1->at(yind) + n2 * p2->at(yind) +
114  n3 * p3->at(yind) + n4 * p4->at(yind)} ;
115 }
116
117 #define POINT_TOL 1.e-6
118
119 int
121 {
122  double x1, x2, x3, x4, y1, y2, y3, y4, a1, a2, a3, a4, b1, b2, b3, b4;
123  double a, b, c, ksi1, ksi2, ksi3, eta1 = 0.0, eta2 = 0.0, denom;
124  int nroot;
125
127
128  x1 = cellgeo.giveVertexCoordinates(1)->at(xind);
129  x2 = cellgeo.giveVertexCoordinates(2)->at(xind);
130  x3 = cellgeo.giveVertexCoordinates(3)->at(xind);
131  x4 = cellgeo.giveVertexCoordinates(4)->at(xind);
132
133  y1 = cellgeo.giveVertexCoordinates(1)->at(yind);
134  y2 = cellgeo.giveVertexCoordinates(2)->at(yind);
135  y3 = cellgeo.giveVertexCoordinates(3)->at(yind);
136  y4 = cellgeo.giveVertexCoordinates(4)->at(yind);
137
138  a1 = x1 + x2 + x3 + x4;
139  a2 = x1 - x2 - x3 + x4;
140  a3 = x1 + x2 - x3 - x4;
141  a4 = x1 - x2 + x3 - x4;
142
143  b1 = y1 + y2 + y3 + y4;
144  b2 = y1 - y2 - y3 + y4;
145  b3 = y1 + y2 - y3 - y4;
146  b4 = y1 - y2 + y3 - y4;
147
148  a = a2 * b4 - b2 * a4;
149  b = a1 * b4 + a2 * b3 - a3 * b2 - b1 * a4 - b4 * 4.0 * coords.at(xind) + a4 * 4.0 * coords.at(yind);
150  c = a1 * b3 - a3 * b1 - 4.0 * coords.at(xind) * b3 + 4.0 * coords.at(yind) * a3;
151
152  // solve quadratic equation for ksi
153  cubic(0.0, a, b, c, & ksi1, & ksi2, & ksi3, & nroot);
154
155  if ( nroot == 0 ) {
157  return false;
158  }
159
160  if ( nroot ) {
161  denom = ( b3 + ksi1 * b4 );
162  if ( fabs(denom) <= 1.0e-10 ) {
163  eta1 = ( 4.0 * coords.at(xind) - a1 - ksi1 * a2 ) / ( a3 + ksi1 * a4 );
164  } else {
165  eta1 = ( 4.0 * coords.at(yind) - b1 - ksi1 * b2 ) / denom;
166  }
167  }
168
169  if ( nroot > 1 ) {
170  double diff_ksi1, diff_eta1, diff_ksi2, diff_eta2, diff1, diff2;
171
172  denom = b3 + ksi2 * b4;
173  if ( fabs(denom) <= 1.0e-10 ) {
174  eta2 = ( 4.0 * coords.at(xind) - a1 - ksi2 * a2 ) / ( a3 + ksi2 * a4 );
175  } else {
176  eta2 = ( 4.0 * coords.at(yind) - b1 - ksi2 * b2 ) / denom;
177  }
178
179  // choose the one which seems to be closer to the parametric space (square <-1;1>x<-1;1>)
180  diff_ksi1 = 0.0;
181  if ( ksi1 > 1.0 ) {
182  diff_ksi1 = ksi1 - 1.0;
183  }
184
185  if ( ksi1 < -1.0 ) {
186  diff_ksi1 = ksi1 + 1.0;
187  }
188
189  diff_eta1 = 0.0;
190  if ( eta1 > 1.0 ) {
191  diff_eta1 = eta1 - 1.0;
192  }
193
194  if ( eta1 < -1.0 ) {
195  diff_eta1 = eta1 + 1.0;
196  }
197
198  diff_ksi2 = 0.0;
199  if ( ksi2 > 1.0 ) {
200  diff_ksi2 = ksi2 - 1.0;
201  }
202
203  if ( ksi2 < -1.0 ) {
204  diff_ksi2 = ksi2 + 1.0;
205  }
206
207  diff_eta2 = 0.0;
208  if ( eta2 > 1.0 ) {
209  diff_eta2 = eta2 - 1.0;
210  }
211
212  if ( eta2 < -1.0 ) {
213  diff_eta2 = eta2 + 1.0;
214  }
215
216  diff1 = diff_ksi1 * diff_ksi1 + diff_eta1 * diff_eta1;
217  diff2 = diff_ksi2 * diff_ksi2 + diff_eta2 * diff_eta2;
218
219  // ksi2, eta2 seems to be closer
220  if ( diff1 > diff2 ) {
221  ksi1 = ksi2;
222  eta1 = eta2;
223  }
224  }
225
228
229  // test if inside
230  bool inside = true;
231  for ( int i = 1; i <= 2; i++ ) {
232  if ( answer.at(i) < ( -1. - POINT_TOL ) ) {
234  inside = false;
235  } else if ( answer.at(i) > ( 1. + POINT_TOL ) ) {
237  inside = false;
238  }
239  }
240
241  return inside;
242 }
243
244 void
245 FEI2dQuadLin :: edgeEvalN(FloatArray &answer, int iedge, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
246 {
247  double ksi = lcoords.at(1);
248  answer = { ( 1. - ksi ) * 0.5, ( 1. + ksi ) * 0.5 };
249 }
250
251 double
252 FEI2dQuadLin :: edgeEvalNormal(FloatArray &answer, int iedge, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
253 {
254  int nodeA, nodeB;
255  IntArray edgeNodes;
256  this->computeLocalEdgeMapping(edgeNodes, iedge);
257  nodeA = edgeNodes.at(1);
258  nodeB = edgeNodes.at(2);
259
261  cellgeo.giveVertexCoordinates(nodeA)->at(yind) - cellgeo.giveVertexCoordinates(nodeB)->at(yind),
262  cellgeo.giveVertexCoordinates(nodeB)->at(xind) - cellgeo.giveVertexCoordinates(nodeA)->at(xind)
263  };
265 }
266
267 void
269  const FloatArray &lcoords, const FEICellGeometry &cellgeo)
270 {
271  IntArray edgeNodes;
272  this->computeLocalEdgeMapping(edgeNodes, iedge);
273  double l = this->edgeComputeLength(edgeNodes, cellgeo);
274
275  answer = { -1.0 / l, 1.0 / l };
276 }
277
278 void
280  const FloatArray &lcoords, const FEICellGeometry &cellgeo)
281 {
282  IntArray edgeNodes;
283  FloatArray n;
284  this->computeLocalEdgeMapping(edgeNodes, iedge);
285  this->edgeEvalN(n, iedge, lcoords, cellgeo);
286
288  answer.at(1) = ( n.at(1) * cellgeo.giveVertexCoordinates( edgeNodes.at(1) )->at(xind) +
289  n.at(2) * cellgeo.giveVertexCoordinates( edgeNodes.at(2) )->at(xind) );
290  answer.at(2) = ( n.at(1) * cellgeo.giveVertexCoordinates( edgeNodes.at(1) )->at(yind) +
291  n.at(2) * cellgeo.giveVertexCoordinates( edgeNodes.at(2) )->at(yind) );
292 }
293
294 void
296 {
297  int aNode = 0, bNode = 0;
298  edgeNodes.resize(2);
299
300  if ( iedge == 1 ) { // edge between nodes 1 2
301  aNode = 1;
302  bNode = 2;
303  } else if ( iedge == 2 ) { // edge between nodes 2 3
304  aNode = 2;
305  bNode = 3;
306  } else if ( iedge == 3 ) { // edge between nodes 3 4
307  aNode = 3;
308  bNode = 4;
309  } else if ( iedge == 4 ) { // edge between nodes 4 1
310  aNode = 4;
311  bNode = 1;
312  } else {
313  OOFEM_ERROR("wrong egde number (%d)", iedge);
314  }
315
316  edgeNodes.at(1) = aNode;
317  edgeNodes.at(2) = bNode;
318 }
319
320 double
322 {
323  double dx, dy;
324  int nodeA, nodeB;
325
326  nodeA = edgeNodes.at(1);
327  nodeB = edgeNodes.at(2);
328
329  dx = cellgeo.giveVertexCoordinates(nodeB)->at(xind) - cellgeo.giveVertexCoordinates(nodeA)->at(xind);
330  dy = cellgeo.giveVertexCoordinates(nodeB)->at(yind) - cellgeo.giveVertexCoordinates(nodeA)->at(yind);
331  return sqrt(dx * dx + dy * dy);
332 }
333
334
335 bool FEI2dQuadLin :: inside(const FloatArray &lcoords) const
336 {
337  const double point_tol = 1.0e-3;
338  bool inside = true;
339  for ( int i = 1; i <= 2; i++ ) {
340  if ( lcoords.at(i) < ( -1. - point_tol ) ) {
341  inside = false;
342  } else if ( lcoords.at(i) > ( 1. + point_tol ) ) {
343  inside = false;
344  }
345  }
346
347  return inside;
348 }
349
350 void FEI2dQuadLin :: evaldNdxi(FloatMatrix &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
351 {
352  const double &ksi = lcoords[0];
353  const double &eta = lcoords[1];
354
356
357  // dn/dxi
358  answer.at(1, 1) = 0.25 * ( 1. + eta );
359  answer.at(2, 1) = -0.25 * ( 1. + eta );
360  answer.at(3, 1) = -0.25 * ( 1. - eta );
361  answer.at(4, 1) = 0.25 * ( 1. - eta );
362
363  // dn/deta
364  answer.at(1, 2) = 0.25 * ( 1. + ksi );
365  answer.at(2, 2) = 0.25 * ( 1. - ksi );
366  answer.at(3, 2) = -0.25 * ( 1. - ksi );
367  answer.at(4, 2) = -0.25 * ( 1. + ksi );
368 }
369
370 double FEI2dQuadLin :: evalNXIntegral(int iEdge, const FEICellGeometry &cellgeo)
371 {
372  IntArray eNodes;
373  const FloatArray *node;
374  double x1, x2, y1, y2;
375
376  this->computeLocalEdgeMapping(eNodes, iEdge);
377
378  node = cellgeo.giveVertexCoordinates( eNodes.at(1) );
379  x1 = node->at(xind);
380  y1 = node->at(yind);
381
382  node = cellgeo.giveVertexCoordinates( eNodes.at(2) );
383  x2 = node->at(xind);
384  y2 = node->at(yind);
385
386  return -( x2 * y1 - x1 * y2 );
387 }
388
391 {
392  IntegrationRule *iRule = new GaussIntegrationRule(1, NULL);
393  int points = iRule->getRequiredNumberOfIntegrationPoints(_Square, order + 2);
394  iRule->SetUpPointsOnSquare(points, _Unknown);
395  return iRule;
396 }
397
398
399 /*
401  */
402
403 double
405 {
406  FloatArray N;
407  this->evalN( N, lcoords, cellgeo);
408
409  double r = 0.0;
410  for ( int i = 1; i <= 4; i++ ) {
411  double x = cellgeo.giveVertexCoordinates(i)->at(1);
412  r += x * N.at(i);
413  }
414
415  return r * FEI2dQuadLin::giveTransformationJacobian(lcoords, cellgeo);
416 }
417
418 double
420  const FEICellGeometry &cellgeo)
421 {
422  IntArray edgeNodes;
423  FloatArray n;
424  this->computeLocalEdgeMapping(edgeNodes, iedge);
425  this->edgeEvalN(n, iedge, lcoords, cellgeo);
426
427  double r = n.at(1)*cellgeo.giveVertexCoordinates(edgeNodes.at(1))->at(1) + n.at(2)*cellgeo.giveVertexCoordinates(edgeNodes.at(2))->at(1);
428  return r * FEI2dQuadLin::edgeGiveTransformationJacobian(iedge, lcoords, cellgeo);
429
430 }
431
432 double
434 {
435  return this->edgeGiveTransformationJacobian(boundary, lcoords, cellgeo);
436 }
437
438 double
440 {
441  return this->edgeGiveTransformationJacobian(boundary, lcoords, cellgeo);
442 }
443
444
445
446 } // end namespace oofem
virtual int global2local(FloatArray &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Default implementation using Newton&#39;s method to find the local coordinates.
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 edgeGiveTransformationJacobian(int iedge, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates the edge Jacobian of transformation between local and global coordinates.
virtual void local2global(FloatArray &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates global coordinates from given local ones.
double & at(int i)
Coefficient access function.
Definition: floatarray.h:131
virtual const FloatArray * giveVertexCoordinates(int i) const =0
virtual double edgeEvalNormal(FloatArray &answer, int iedge, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates the normal on the given edge.
Class representing a general abstraction for cell geometry.
Definition: feinterpol.h:62
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 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...
virtual int SetUpPointsOnSquare(int, MaterialMode mode)
Sets up receiver&#39;s integration points on unit square integration domain.
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.
void cubic(double a, double b, double c, double d, double *r1, double *r2, double *r3, int *num)
Solves cubic equation for real roots.
Definition: mathfem.C:43
virtual double giveTransformationJacobian(const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates the determinant of the transformation.
virtual void computeLocalEdgeMapping(IntArray &edgeNodes, int iedge)
#define OOFEM_ERROR(...)
Definition: error.h:61
virtual double boundaryEdgeGiveTransformationJacobian(int boundary, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates the determinant of the transformation Jacobian on the requested boundary.
#define N(p, q)
Definition: mdm.C:367
virtual IntegrationRule * giveIntegrationRule(int order)
Sets up a suitable integration rule for numerical integrating over volume.
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
double edgeComputeLength(IntArray &edgeNodes, const FEICellGeometry &cellgeo)
Implementation of matrix containing floating point numbers.
Definition: floatmatrix.h:94
#define POINT_TOL
virtual void edgeLocal2global(FloatArray &answer, int iedge, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates edge global coordinates from given local ones.
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
virtual double evalNXIntegral(int iEdge, const FEICellGeometry &cellgeo)
Computes the integral .
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
virtual void evaldNdxi(FloatMatrix &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates the matrix of derivatives of interpolation functions (shape functions) at given point...
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 double evaldNdx(FloatMatrix &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates the matrix of derivatives of interpolation functions (shape functions) at given point...
int giveNumberOfRows() const
Returns number of rows of receiver.
Definition: floatmatrix.h:156
virtual double giveArea(const FEICellGeometry &cellgeo) const
Computes the exact area.
virtual double boundaryGiveTransformationJacobian(int boundary, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates the determinant of the transformation Jacobian on the requested boundary.