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 FEI2dLineQuad :: evalN(FloatArray &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
43 {
44  double xi = lcoords(0);
46  answer(0) = 0.5 * ( xi - 1.0 ) * xi;
47  answer(1) = 0.5 * ( xi + 1.0 ) * xi;
48  answer(2) = 1.0 - xi * xi;
49 }
50
51 double FEI2dLineQuad :: evaldNdx(FloatMatrix &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
52 {
53  // Not meaningful to return anything.
55  return 0.;
56 }
57
58 void
60 {
61  double xi = lcoords(0);
63  answer.at(1, 1) = xi - 0.5;
64  answer.at(2, 1) = xi + 0.5;
65  answer.at(3, 1) = -2.0 * xi;
66 }
67
68 void FEI2dLineQuad :: local2global(FloatArray &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
69 {
70  FloatArray n;
71  this->evalN(n, lcoords, cellgeo);
73  answer.at(1) = ( n(0) * cellgeo.giveVertexCoordinates(1)->at(xind) +
74  n(1) * cellgeo.giveVertexCoordinates(2)->at(xind) +
75  n(2) * cellgeo.giveVertexCoordinates(3)->at(xind) );
76  answer.at(2) = ( n(0) * cellgeo.giveVertexCoordinates(1)->at(yind) +
77  n(1) * cellgeo.giveVertexCoordinates(2)->at(yind) +
78  n(2) * cellgeo.giveVertexCoordinates(3)->at(yind) );
79 }
80
81 int FEI2dLineQuad :: global2local(FloatArray &answer, const FloatArray &gcoords, const FEICellGeometry &cellgeo)
82 {
83  double x1_x2, y1_y2, px_x3, py_y3, x3_x2_x1, y3_y2_y1;
84  double b0, b1, b2, b3;
85
86  x1_x2 = cellgeo.giveVertexCoordinates(1)->at(xind) - cellgeo.giveVertexCoordinates(2)->at(xind);
87  y1_y2 = cellgeo.giveVertexCoordinates(1)->at(yind) - cellgeo.giveVertexCoordinates(2)->at(yind);
88  px_x3 = gcoords.at(1) - cellgeo.giveVertexCoordinates(3)->at(xind);
89  py_y3 = gcoords.at(2) - cellgeo.giveVertexCoordinates(3)->at(yind);
90  x3_x2_x1 = 2 * cellgeo.giveVertexCoordinates(3)->at(xind) - cellgeo.giveVertexCoordinates(2)->at(xind) - cellgeo.giveVertexCoordinates(1)->at(xind);
91  y3_y2_y1 = 2 * cellgeo.giveVertexCoordinates(3)->at(yind) - cellgeo.giveVertexCoordinates(2)->at(yind) - cellgeo.giveVertexCoordinates(1)->at(yind);
92
93  b0 = 0.50 * ( x1_x2 * px_x3 + y1_y2 * py_y3 );
94  b1 = 0.25 * ( x1_x2 * x1_x2 + y1_y2 * y1_y2 ) + x3_x2_x1 * px_x3 + y3_y2_y1 * py_y3;
95  b2 = 0.75 * ( x1_x2 * x3_x2_x1 + y1_y2 * y3_y2_y1 );
96  b3 = 0.50 * ( x3_x2_x1 * x3_x2_x1 + y3_y2_y1 * y3_y2_y1 );
97
98  double r [ 3 ];
99  int roots;
100  cubic(b3, b2, b1, b0, & r [ 0 ], & r [ 1 ], & r [ 2 ], & roots);
101
102  // Copy them over, along with boundary cases.
103  int points = 2;
104  double p [ 5 ] = {
105  -1.0, 1.0
106  };
107  for ( int i = 0; i < roots; i++ ) {
108  if ( r [ i ] > -1.0 && r [ i ] < 1.0 ) {
109  // The cubic solver has pretty bad accuracy, performing a single newton iteration
110  // which typically improves the solution by many many orders of magnitude.
112  r [ i ] -= ( b0 + b1 * r [ i ] + b2 * r [ i ] * r [ i ] + b3 * r [ i ] * r [ i ] * r [ i ] ) / ( b1 + 2 * b2 * r [ i ] + 3 * b3 * r [ i ] * r [ i ] );
113  p [ points ] = r [ i ];
114  points++;
115  }
116  }
117
118  double min_distance2 = 0.0, min_xi = 0, distance2;
119  FloatArray f(2);
121
122  for ( int i = 0; i < points; i++ ) {
123  answer(0) = p [ i ];
125  distance2 = f.distance_square(gcoords);
126  if ( i == 0 || distance2 < min_distance2 ) {
127  min_distance2 = distance2;
129  }
130  }
131
132  answer(0) = clamp(min_xi, -1., 1.); // Make sure to only give coordinates within the geometry.
133
134  return false; // It will never land exactly on the geometry.
135 }
136
138 {
139  edgeNodes = { 1, 2, 3};
140 }
141
142 void FEI2dLineQuad :: edgeEvalN(FloatArray &answer, int iedge, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
143 {
145 }
146
148  const FloatArray &lcoords, const FEICellGeometry &cellgeo)
149 {
150  double xi = lcoords(0);
152  answer(0) = -0.5 + xi;
153  answer(1) = 0.5 + xi;
154  answer(2) = -2.0 * xi;
155
156  double es1 = answer(0) * cellgeo.giveVertexCoordinates(1)->at(xind) +
159
160  double es2 = answer(0) * cellgeo.giveVertexCoordinates(1)->at(yind) +
163
164  double J = sqrt(es1 * es1 + es2 * es2);
166  //return J;
167 }
168
169 double FEI2dLineQuad :: edgeEvalNormal(FloatArray &normal, int iedge, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
170 {
171  double xi = lcoords(0);
172  double dN1dxi = -0.5 + xi;
173  double dN2dxi = 0.5 + xi;
174  double dN3dxi = -2.0 * xi;
175
176  normal.resize(2);
177
178  normal.at(1) = -dN1dxi *cellgeo.giveVertexCoordinates(1)->at(yind) +
179  - dN2dxi *cellgeo.giveVertexCoordinates(2)->at(yind) +
180  - dN3dxi *cellgeo.giveVertexCoordinates(3)->at(yind);
181
182  normal.at(2) = dN1dxi * cellgeo.giveVertexCoordinates(1)->at(xind) +
183  dN2dxi *cellgeo.giveVertexCoordinates(2)->at(xind) +
184  dN3dxi *cellgeo.giveVertexCoordinates(3)->at(xind);
185
186  return normal.normalize();
187 }
188
190  const FloatArray &lcoords, const FEICellGeometry &cellgeo)
191 {
193 }
194
196 {
197  double xi = lcoords(0);
198  double a1 = -0.5 + xi;
199  double a2 = 0.5 + xi;
200  double a3 = -2.0 * xi;
201
202  double es1 = a1 * cellgeo.giveVertexCoordinates(1)->at(xind) +
203  a2 *cellgeo.giveVertexCoordinates(2)->at(xind) +
204  a3 *cellgeo.giveVertexCoordinates(3)->at(xind);
205  double es2 = a1 * cellgeo.giveVertexCoordinates(1)->at(yind) +
206  a2 *cellgeo.giveVertexCoordinates(2)->at(yind) +
207  a3 *cellgeo.giveVertexCoordinates(3)->at(yind);
208
209  return sqrt(es1 * es1 + es2 * es2);
210 }
211
212 void FEI2dLineQuad :: giveJacobianMatrixAt(FloatMatrix &jacobianMatrix, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
213 {
214  double xi = lcoords(0);
215  double dN1dxi = -0.5 + xi;
216  double dN2dxi = 0.5 + xi;
217  double dN3dxi = -2.0 * xi;
218
219  double es1 = dN1dxi * cellgeo.giveVertexCoordinates(1)->at(xind) +
220  dN2dxi *cellgeo.giveVertexCoordinates(2)->at(xind) +
221  dN3dxi *cellgeo.giveVertexCoordinates(3)->at(xind);
222  double es2 = dN1dxi * cellgeo.giveVertexCoordinates(1)->at(yind) +
223  dN2dxi *cellgeo.giveVertexCoordinates(2)->at(yind) +
224  dN3dxi *cellgeo.giveVertexCoordinates(3)->at(yind);
225
226  double J = sqrt(es1 * es1 + es2 * es2);
227
228  // Only used for determined deformation of element, not sure about the transpose here;
229  jacobianMatrix.resize(2, 2);
230  jacobianMatrix(0, 0) = es1 / J;
231  jacobianMatrix(0, 1) = es2 / J;
232  jacobianMatrix(1, 0) = -es2 / J;
233  jacobianMatrix(1, 1) = es1 / J;
234 }
235
237 {
238  OOFEM_ERROR("Not implemented");
239  return 0.0;
240 }
241
242 double FEI2dLineQuad :: evalNXIntegral(int iEdge, const FEICellGeometry &cellgeo)
243 {
244  const FloatArray *node;
245  double x1, x2, x3, y1, y2, y3;
246
247  node = cellgeo.giveVertexCoordinates(1);
248  x1 = node->at(xind);
249  y1 = node->at(yind);
250
251  node = cellgeo.giveVertexCoordinates(2);
252  x2 = node->at(xind);
253  y2 = node->at(yind);
254
255  node = cellgeo.giveVertexCoordinates(3);
256  x3 = node->at(xind);
257  y3 = node->at(yind);
258
259  return ( x1 * y2 - x2 * y1 + 4 * ( x3 * ( y1 - y2 ) + y3 * ( x2 - x1 ) ) ) / 3.0;
260 }
261
263 {
264  IntegrationRule *iRule = new GaussIntegrationRule(1, NULL);
265  int points = iRule->getRequiredNumberOfIntegrationPoints(_Line, order + 1);
266  iRule->SetUpPointsOnLine(points, _Unknown);
267  return iRule;
268 }
269 } // end namespace oofem
virtual int global2local(FloatArray &answer, const FloatArray &gcoords, const FEICellGeometry &cellgeo)
Default implementation using Newton&#39;s method to find the local coordinates.
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 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
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 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 evalNXIntegral(int iEdge, const FEICellGeometry &cellgeo)
Computes the integral .
virtual double giveTransformationJacobian(const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates the determinant of the transformation.
virtual double evaldNdx(FloatMatrix &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates the matrix of derivatives of interpolation functions (shape functions) at given point...
#define OOFEM_ERROR(...)
Definition: error.h:61
virtual void evalN(FloatArray &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates the array of interpolation functions (shape functions) at given point.
double at(int i, int j) const
Coefficient access function.
Definition: floatmatrix.h:176
double distance_square(const FloatArray &iP1, const FloatArray &iP2, double &oXi, double &oXiUnbounded) const
Definition: floatarray.C:499
Class representing vector of real numbers.
Definition: floatarray.h:82
Implementation of matrix containing floating point numbers.
Definition: floatmatrix.h:94
double clamp(int a, int lower, int upper)
Returns the clamped value of a between upper and lower.
Definition: mathfem.h:75
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 void local2global(FloatArray &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates global coordinates from given local ones.
virtual void edgeLocal2global(FloatArray &answer, int iedge, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates edge global coordinates from given local ones.
void giveJacobianMatrixAt(FloatMatrix &jacobianMatrix, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Gives the jacobian matrix at the local coordinates.
virtual int SetUpPointsOnLine(int, MaterialMode mode)
Sets up receiver&#39;s integration points on unit line integration domain.
virtual void computeLocalEdgeMapping(IntArray &edgeNodes, int iedge)
void times(double s)
Definition: floatarray.C:818
virtual double edgeEvalNormal(FloatArray &normal, int iedge, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates the normal on the given edge.
virtual IntegrationRule * giveIntegrationRule(int order)
Sets up a suitable integration rule for numerical integrating over volume.
the oofem namespace is to define a context or scope in which all oofem names are defined.
void clear()
Sets size of receiver to be an empty matrix. It will have zero rows and zero columns size...
Definition: floatmatrix.h:516
double edgeComputeLength(IntArray &edgeNodes, const FEICellGeometry &cellgeo)
double normalize()