OOFEM  2.4 OOFEM.org - Object Oriented Finite Element Solver
fei3dwedgelin.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
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 "fei3dwedgelin.h"
36 #include "floatarray.h"
37 #include "floatmatrix.h"
38 #include "intarray.h"
39 #include "gaussintegrationrule.h"
40
41 namespace oofem {
42 void
43 FEI3dWedgeLin :: evalN(FloatArray &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
44 {
45  double u, v, w;
47
48  u = lcoords.at(1);
49  v = lcoords.at(2);
50  w = lcoords.at(3);
51
52  answer.at(1) = 0.5 * ( 1. - w ) * ( 1. - u - v );
53  answer.at(2) = 0.5 * ( 1. - w ) * u;
54  answer.at(3) = 0.5 * ( 1. - w ) * v;
55  answer.at(4) = 0.5 * ( 1. + w ) * ( 1. - u - v );
56  answer.at(5) = 0.5 * ( 1. + w ) * u;
57  answer.at(6) = 0.5 * ( 1. + w ) * v;
58 }
59
60
61 double
62 FEI3dWedgeLin :: evaldNdx(FloatMatrix &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
63 {
64  FloatMatrix jacobianMatrix, inv, dNduvw, coords;
65
66  this->giveLocalDerivative(dNduvw, lcoords);
67  coords.resize(3, 6);
68  for ( int i = 1; i <= 6; i++ ) {
69  coords.setColumn(* cellgeo.giveVertexCoordinates(i), i);
70  }
71  jacobianMatrix.beProductOf(coords, dNduvw);
72  inv.beInverseOf(jacobianMatrix);
73
75  return jacobianMatrix.giveDeterminant();
76 }
77
78
79 void
80 FEI3dWedgeLin :: local2global(FloatArray &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
81 {
82  FloatArray n;
83
84  this->evalN(n, lcoords, cellgeo);
85
87  for ( int i = 1; i <= 6; i++ ) {
89  }
90 }
91
92
94 {
95  const FloatArray *n1 = cellgeo.giveVertexCoordinates(1);
96  const FloatArray *n2 = cellgeo.giveVertexCoordinates(6);
98  return n1->distance(n2);
99 }
100
101 #define POINT_TOL 1.e-3
102
103 int
104 FEI3dWedgeLin :: global2local(FloatArray &answer, const FloatArray &gcoords, const FEICellGeometry &cellgeo)
105 {
106  FloatArray res, delta, guess;
107  FloatMatrix jac;
108  double convergence_limit, error = 0.0;
109
110  // find a suitable convergence limit
111  convergence_limit = 1e-6 * this->giveCharacteristicLength(cellgeo);
112
113  // setup initial guess
116
117  // apply Newton-Raphson to solve the problem
118  for ( int nite = 0; nite < 10; nite++ ) {
119  // compute the residual
121  res.beDifferenceOf(gcoords, guess);
122
123  // check for convergence
124  error = res.computeNorm();
125  if ( error < convergence_limit ) {
126  break;
127  }
128
129  // compute the corrections
131  jac.solveForRhs(res, delta);
132
133  // update guess
135  }
136  if ( error > convergence_limit ) { // Imperfect, could give false negatives.
137  //OOFEM_ERROR("no convergence after 10 iterations");
138  answer = {1. / 3., 1. / 3., 1. / 3.};
139  return false;
140  }
141
142  // check limits for each local coordinate [-1,1] for quadrilaterals. (different for other elements, typically [0,1]).
143  bool inside = true;
144  for ( int i = 1; i <= answer.giveSize(); i++ ) {
145  if ( answer.at(i) < ( 0. - POINT_TOL ) ) {
147  inside = false;
148  } else if ( answer.at(i) > ( 1. + POINT_TOL ) ) {
150  inside = false;
151  }
152  }
153
154  return inside;
155 }
156
157
158 double
160 {
161  FloatMatrix jacobianMatrix;
162
163  this->giveJacobianMatrixAt(jacobianMatrix, lcoords, cellgeo);
164  return jacobianMatrix.giveDeterminant() / 2.;
165 }
166
167
168 void
169 FEI3dWedgeLin :: giveJacobianMatrixAt(FloatMatrix &jacobianMatrix, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
170 // Returns the jacobian matrix J (x,y,z)/(ksi,eta,dzeta) of the receiver.
171 {
172  FloatMatrix dNduvw, coords;
173  this->giveLocalDerivative(dNduvw, lcoords);
174  coords.resize(3, 6);
175  for ( int i = 1; i <= 6; i++ ) {
176  coords.setColumn(* cellgeo.giveVertexCoordinates(i), i);
177  }
178  jacobianMatrix.beProductOf(coords, dNduvw);
179 }
180
181
182 void
184 {
185  double u, v, w;
186  u = lcoords.at(1);
187  v = lcoords.at(2);
188  w = lcoords.at(3);
189
190  dN.resize(6, 3);
191
192  dN.at(1, 1) = -0.5 * ( 1. - w );
193  dN.at(2, 1) = 0.5 * ( 1. - w );
194  dN.at(3, 1) = 0.;
195  dN.at(4, 1) = -0.5 * ( 1. + w );
196  dN.at(5, 1) = 0.5 * ( 1. + w );
197  dN.at(6, 1) = 0.;
198
199  dN.at(1, 2) = -0.5 * ( 1. - w );
200  dN.at(2, 2) = 0.;
201  dN.at(3, 2) = 0.5 * ( 1. - w );
202  dN.at(4, 2) = -0.5 * ( 1. + w );
203  dN.at(5, 2) = 0.;
204  dN.at(6, 2) = 0.5 * ( 1. + w );
205
206  dN.at(1, 3) = -0.5 * ( 1. - u - v );
207  dN.at(2, 3) = -0.5 * u;
208  dN.at(3, 3) = -0.5 * v;
209  dN.at(4, 3) = 0.5 * ( 1. - u - v );
210  dN.at(5, 3) = 0.5 * u;
211  dN.at(6, 3) = 0.5 * v;
212 }
213
214
215 void
216 FEI3dWedgeLin :: edgeEvalN(FloatArray &answer, int iedge, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
217 {
218  double ksi = lcoords.at(1);
220  answer.at(1) = ( 1. - ksi ) * 0.5;
221  answer.at(2) = ( 1. - ksi ) * 0.5;
222 }
223
224
225 void
226 FEI3dWedgeLin :: edgeEvaldNdx(FloatMatrix &answer, int iedge, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
227 {
228  OOFEM_ERROR("not implemented");
229 }
230
231
232 void
233 FEI3dWedgeLin :: edgeLocal2global(FloatArray &answer, int iedge, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
234 {
235  IntArray nodes;
236  FloatArray n;
237
238  this->computeLocalEdgeMapping(nodes, iedge);
239  this->edgeEvalN(n, iedge, lcoords, cellgeo);
240
242  for ( int i = 1; i <= n.giveSize(); ++i ) {
244  }
245 }
246
247
248 void
250 {
251  if ( iedge == 1 ) {
252  edgeNodes = {1, 2};
253  } else if ( iedge == 2 ) {
254  edgeNodes = {2, 3};
255  } else if ( iedge == 3 ) {
256  edgeNodes = {3, 1};
257  } else if ( iedge == 4 ) {
258  edgeNodes = {4, 5};
259  } else if ( iedge == 5 ) {
260  edgeNodes = {5, 6};
261  } else if ( iedge == 6 ) {
262  edgeNodes = {6, 4};
263  } else if ( iedge == 7 ) {
264  edgeNodes = {1, 4};
265  } else if ( iedge == 8 ) {
266  edgeNodes = {2, 5};
267  } else if ( iedge == 9 ) {
268  edgeNodes = {3, 6};
269  } else {
270  OOFEM_ERROR("Edge %d doesn't exist.\n", iedge);
271  }
272 }
273
274
275 double
277 {
278  OOFEM_ERROR("not implemented");
279  return 0.0;
280 }
281
282
283 void
284 FEI3dWedgeLin :: surfaceEvalN(FloatArray &answer, int isurf, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
285 {
286  double ksi = lcoords.at(1);
287  double eta = lcoords.at(2);
288
289  if ( isurf <= 2 ) {
293  answer.at(3) = 1.0 - ksi - eta;
294  } else {
296  answer.at(1) = ( 1. + ksi ) * ( 1. + eta ) * 0.25;
297  answer.at(2) = ( 1. - ksi ) * ( 1. + eta ) * 0.25;
298  answer.at(3) = ( 1. - ksi ) * ( 1. - eta ) * 0.25;
299  answer.at(4) = ( 1. + ksi ) * ( 1. - eta ) * 0.25;
300  }
301 }
302
303
304 void
306  const FloatArray &lcoords, const FEICellGeometry &cellgeo)
307 {
308  IntArray nodes;
309  FloatArray n;
310
311  this->computeLocalSurfaceMapping(nodes, isurf);
312  this->surfaceEvalN(n, isurf, lcoords, cellgeo);
313
315  for ( int i = 1; i <= n.giveSize(); ++i ) {
317  }
318 }
319
320
321 void
323 {
324  if ( isurf == 1 ) {
325  nodes = {1, 2, 3};
326  } else if ( isurf == 2 ) {
327  nodes = {4, 5, 6};
328  } else if ( isurf == 3 ) {
329  nodes = {1, 2, 5, 4};
330  } else if ( isurf == 4 ) {
331  nodes = {2, 3, 6, 5};
332  } else if ( isurf == 5 ) {
333  nodes = {3, 1, 4, 6};
334  } else {
335  OOFEM_ERROR("Surface %d doesn't exist.\n", isurf);
336  }
337 }
338
339
340 double
342  const FEICellGeometry &cellgeo)
343 {
344  OOFEM_ERROR("not implemented");
345  return 0;
346 }
347
348
351 {
352  IntegrationRule *iRule = new GaussIntegrationRule(1, NULL);
354  //int points = iRule->getRequiredNumberOfIntegrationPoints(_Wedge, order);
355  OOFEM_WARNING("Warning.. ignoring 'order' argument: FIXME");
356  int pointsZeta = 1;
357  int pointsTriangle = 1;
358  iRule->SetUpPointsOnWedge(pointsTriangle, pointsZeta, _Unknown);
359  return iRule;
360 }
361
362
365 {
366  IntegrationRule *iRule = new GaussIntegrationRule(1, NULL);
367  if ( boundary <= 2 ) {
368  int points = iRule->getRequiredNumberOfIntegrationPoints(_Triangle, order + 0);
369  iRule->SetUpPointsOnTriangle(points, _Unknown);
370  } else {
371  int points = iRule->getRequiredNumberOfIntegrationPoints(_Square, order + 2);
372  iRule->SetUpPointsOnSquare(points, _Unknown);
373  }
374  return iRule;
375 }
376 } // end namespace oofem
double giveDeterminant() const
Returns the trace (sum of diagonal components) of the receiver.
Definition: floatmatrix.C:1408
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 void surfaceLocal2global(FloatArray &answer, int iedge, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates edge global coordinates from given local ones.
virtual double giveCharacteristicLength(const FEICellGeometry &cellgeo) const
Definition: fei3dwedgelin.C:93
virtual int global2local(FloatArray &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates local coordinates from given global ones.
virtual IntegrationRule * giveIntegrationRule(int order)
Sets up a suitable integration rule for numerical integrating over volume.
bool solveForRhs(const FloatArray &b, FloatArray &answer, bool transpose=false)
Solves the system of linear equations .
Definition: floatmatrix.C:1112
virtual double giveTransformationJacobian(const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates the determinant of the transformation.
double & at(int i)
Coefficient access function.
Definition: floatarray.h:131
virtual const FloatArray * giveVertexCoordinates(int i) const =0
virtual int SetUpPointsOnWedge(int nPointsTri, int nPointsDepth, MaterialMode mode)
Sets up receiver&#39;s integration points on a wedge integration domain.
Class representing a general abstraction for cell geometry.
Definition: feinterpol.h:62
void clear()
Definition: floatarray.h:206
virtual void local2global(FloatArray &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates global coordinates from given local ones.
Definition: fei3dwedgelin.C:80
virtual void computeLocalEdgeMapping(IntArray &edgeNodes, int iedge)
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
virtual double edgeGiveTransformationJacobian(int iedge, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates the edge jacobian of transformation between local and global coordinates.
int & at(int i)
Coefficient access function.
Definition: intarray.h:103
Abstract base class representing integration rule.
void beDifferenceOf(const FloatArray &a, const FloatArray &b)
Sets receiver to be a - b.
Definition: floatarray.C:341
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 void giveJacobianMatrixAt(FloatMatrix &jacobianMatrix, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Gives the jacobian matrix at the local coordinates.
double at(int i, int j) const
Coefficient access function.
Definition: floatmatrix.h:176
virtual double surfaceGiveTransformationJacobian(int isurf, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates the edge jacobian of transformation between local and global coordinates.
Class representing vector of real numbers.
Definition: floatarray.h:82
Implementation of matrix containing floating point numbers.
Definition: floatmatrix.h:94
virtual void edgeEvaldNdx(FloatMatrix &answer, int iedge, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates the matrix of derivatives of edge interpolation functions (shape functions) at given point...
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 giveLocalDerivative(FloatMatrix &dN, const FloatArray &lcoords)
virtual IntegrationRule * giveBoundaryIntegrationRule(int order, int boundary)
Sets up a suitable integration rule for integrating over the requested boundary.
void zero()
Definition: floatarray.C:658
void setColumn(const FloatArray &src, int c)
Sets the values of the matrix in specified column.
Definition: floatmatrix.C:648
virtual int SetUpPointsOnTriangle(int, MaterialMode mode)
Sets up receiver&#39;s integration points on triangular (area coords) integration domain.
virtual void surfaceEvalN(FloatArray &answer, int isurf, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates the array of edge interpolation functions (shape functions) at given point.
void beProductOf(const FloatMatrix &a, const FloatMatrix &b)
Assigns to the receiver product of .
Definition: floatmatrix.C:337
int giveSize() const
Definition: floatarray.h:218
the oofem namespace is to define a context or scope in which all oofem names are defined.
virtual double evaldNdx(FloatMatrix &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates the matrix of derivatives of interpolation functions (shape functions) at given point...
Definition: fei3dwedgelin.C:62
void beInverseOf(const FloatMatrix &src)
Modifies receiver to become inverse of given parameter.
Definition: floatmatrix.C:835
#define OOFEM_WARNING(...)
Definition: error.h:62
virtual void evalN(FloatArray &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates the array of interpolation functions (shape functions) at given point.
Definition: fei3dwedgelin.C:43