OOFEM  2.4 OOFEM.org - Object Oriented Finite Element Solver
fei3dhexalin.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 "fei3dhexalin.h"
36 #include "mathfem.h"
37 #include "floatmatrix.h"
38 #include "floatarray.h"
39 #include "gaussintegrationrule.h"
40
41 namespace oofem {
42 void
43 FEI3dHexaLin :: evalN(FloatArray &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
44 {
45  double x, y, z;
47
48  x = lcoords.at(1);
49  y = lcoords.at(2);
50  z = lcoords.at(3);
51
52  answer.at(1) = 0.125 * ( 1. - x ) * ( 1. - y ) * ( 1. + z );
53  answer.at(2) = 0.125 * ( 1. - x ) * ( 1. + y ) * ( 1. + z );
54  answer.at(3) = 0.125 * ( 1. + x ) * ( 1. + y ) * ( 1. + z );
55  answer.at(4) = 0.125 * ( 1. + x ) * ( 1. - y ) * ( 1. + z );
56  answer.at(5) = 0.125 * ( 1. - x ) * ( 1. - y ) * ( 1. - z );
57  answer.at(6) = 0.125 * ( 1. - x ) * ( 1. + y ) * ( 1. - z );
58  answer.at(7) = 0.125 * ( 1. + x ) * ( 1. + y ) * ( 1. - z );
59  answer.at(8) = 0.125 * ( 1. + x ) * ( 1. - y ) * ( 1. - z );
60 }
61
62 double
63 FEI3dHexaLin :: evaldNdx(FloatMatrix &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
64 {
65  FloatMatrix jacobianMatrix, inv, dNduvw, coords;
66
67  this->giveLocalDerivative(dNduvw, lcoords);
68  coords.resize(3, 8);
69  for ( int i = 1; i <= 8; i++ ) {
70  coords.setColumn(* cellgeo.giveVertexCoordinates(i), i);
71  }
72  jacobianMatrix.beProductOf(coords, dNduvw);
73  inv.beInverseOf(jacobianMatrix);
74
76  return jacobianMatrix.giveDeterminant();
77 }
78
79 void
80 FEI3dHexaLin :: local2global(FloatArray &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
81 {
82  FloatArray n;
83  this->evalN(n, lcoords, cellgeo);
84
86  for ( int i = 1; i <= 8; i++ ) {
88  }
89 }
90
91 #define POINT_TOL 1.e-3
92
93 int
94 FEI3dHexaLin :: global2local(FloatArray &answer, const FloatArray &coords, const FEICellGeometry &cellgeo)
95 {
96  double x1, x2, x3, x4, x5, x6, x7, x8, a1, a2, a3, a4, a5, a6, a7, a8;
97  double y1, y2, y3, y4, y5, y6, y7, y8, b1, b2, b3, b4, b5, b6, b7, b8;
98  double z1, z2, z3, z4, z5, z6, z7, z8, c1, c2, c3, c4, c5, c6, c7, c8;
99  double xp, yp, zp, u, v, w;
100  FloatMatrix p(3, 3);
101  FloatArray r(3), delta;
102  int nite = 0;
103
104  x1 = cellgeo.giveVertexCoordinates(1)->at(1);
105  x2 = cellgeo.giveVertexCoordinates(2)->at(1);
106  x3 = cellgeo.giveVertexCoordinates(3)->at(1);
107  x4 = cellgeo.giveVertexCoordinates(4)->at(1);
108  x5 = cellgeo.giveVertexCoordinates(5)->at(1);
109  x6 = cellgeo.giveVertexCoordinates(6)->at(1);
110  x7 = cellgeo.giveVertexCoordinates(7)->at(1);
111  x8 = cellgeo.giveVertexCoordinates(8)->at(1);
112
113  y1 = cellgeo.giveVertexCoordinates(1)->at(2);
114  y2 = cellgeo.giveVertexCoordinates(2)->at(2);
115  y3 = cellgeo.giveVertexCoordinates(3)->at(2);
116  y4 = cellgeo.giveVertexCoordinates(4)->at(2);
117  y5 = cellgeo.giveVertexCoordinates(5)->at(2);
118  y6 = cellgeo.giveVertexCoordinates(6)->at(2);
119  y7 = cellgeo.giveVertexCoordinates(7)->at(2);
120  y8 = cellgeo.giveVertexCoordinates(8)->at(2);
121
122  z1 = cellgeo.giveVertexCoordinates(1)->at(3);
123  z2 = cellgeo.giveVertexCoordinates(2)->at(3);
124  z3 = cellgeo.giveVertexCoordinates(3)->at(3);
125  z4 = cellgeo.giveVertexCoordinates(4)->at(3);
126  z5 = cellgeo.giveVertexCoordinates(5)->at(3);
127  z6 = cellgeo.giveVertexCoordinates(6)->at(3);
128  z7 = cellgeo.giveVertexCoordinates(7)->at(3);
129  z8 = cellgeo.giveVertexCoordinates(8)->at(3);
130
131  xp = coords.at(1);
132  yp = coords.at(2);
133  zp = coords.at(3);
134
135  a1 = x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8;
136  a2 = -x1 - x2 + x3 + x4 - x5 - x6 + x7 + x8;
137  a3 = -x1 + x2 + x3 - x4 - x5 + x6 + x7 - x8;
138  a4 = x1 + x2 + x3 + x4 - x5 - x6 - x7 - x8;
139  a5 = x1 - x2 + x3 - x4 + x5 - x6 + x7 - x8;
140  a6 = -x1 - x2 + x3 + x4 + x5 + x6 - x7 - x8;
141  a7 = -x1 + x2 + x3 - x4 + x5 - x6 - x7 + x8;
142  a8 = x1 - x2 + x3 - x4 - x5 + x6 - x7 + x8;
143
144  b1 = y1 + y2 + y3 + y4 + y5 + y6 + y7 + y8;
145  b2 = -y1 - y2 + y3 + y4 - y5 - y6 + y7 + y8;
146  b3 = -y1 + y2 + y3 - y4 - y5 + y6 + y7 - y8;
147  b4 = y1 + y2 + y3 + y4 - y5 - y6 - y7 - y8;
148  b5 = y1 - y2 + y3 - y4 + y5 - y6 + y7 - y8;
149  b6 = -y1 - y2 + y3 + y4 + y5 + y6 - y7 - y8;
150  b7 = -y1 + y2 + y3 - y4 + y5 - y6 - y7 + y8;
151  b8 = y1 - y2 + y3 - y4 - y5 + y6 - y7 + y8;
152
153  c1 = z1 + z2 + z3 + z4 + z5 + z6 + z7 + z8;
154  c2 = -z1 - z2 + z3 + z4 - z5 - z6 + z7 + z8;
155  c3 = -z1 + z2 + z3 - z4 - z5 + z6 + z7 - z8;
156  c4 = z1 + z2 + z3 + z4 - z5 - z6 - z7 - z8;
157  c5 = z1 - z2 + z3 - z4 + z5 - z6 + z7 - z8;
158  c6 = -z1 - z2 + z3 + z4 + z5 + z6 - z7 - z8;
159  c7 = -z1 + z2 + z3 - z4 + z5 - z6 - z7 + z8;
160  c8 = z1 - z2 + z3 - z4 - z5 + z6 - z7 + z8;
161
162  // setup initial guess
165
166  // apply Newton-Raphson to solve the problem
167  for ( ;; ) {
168  if ( ( ++nite ) > 10 ) {
170  return false;
171  }
172
176
177  // compute the residual
178  r.at(1) = a1 + u * a2 + v * a3 + w * a4 + u * v * a5 + u * w * a6 + v * w * a7 + u * v * w * a8 - 8.0 * xp;
179  r.at(2) = b1 + u * b2 + v * b3 + w * b4 + u * v * b5 + u * w * b6 + v * w * b7 + u * v * w * b8 - 8.0 * yp;
180  r.at(3) = c1 + u * c2 + v * c3 + w * c4 + u * v * c5 + u * w * c6 + v * w * c7 + u * v * w * c8 - 8.0 * zp;
181
182  // check for convergence
183  if ( r.computeSquaredNorm() < 1.e-20 ) {
184  break; // sqrt(1.e-20) = 1.e-10
185  }
186
187  p.at(1, 1) = a2 + v * a5 + w * a6 + v * w * a8;
188  p.at(1, 2) = a3 + u * a5 + w * a7 + u * w * a8;
189  p.at(1, 3) = a4 + u * a6 + v * a7 + u * v * a8;
190
191  p.at(2, 1) = b2 + v * b5 + w * b6 + v * w * b8;
192  p.at(2, 2) = b3 + u * b5 + w * b7 + u * w * b8;
193  p.at(2, 3) = b4 + u * b6 + v * b7 + u * v * b8;
194
195  p.at(3, 1) = c2 + v * c5 + w * c6 + v * w * c8;
196  p.at(3, 2) = c3 + u * c5 + w * c7 + u * w * c8;
197  p.at(3, 3) = c4 + u * c6 + v * c7 + u * v * c8;
198
199  // solve for corrections
200  p.solveForRhs(r, delta);
201
202  // update guess
204  }
205
206  // test if inside
207  bool inside = true;
208  for ( int i = 1; i <= 3; i++ ) {
209  if ( answer.at(i) < ( -1. - POINT_TOL ) ) {
211  inside = false;
212  } else if ( answer.at(i) > ( 1. + POINT_TOL ) ) {
214  inside = false;
215  }
216  }
217
218  return inside;
219 }
220
221 void
222 FEI3dHexaLin :: edgeEvalN(FloatArray &answer, int iedge, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
223 {
224  double ksi = lcoords.at(1);
226
227  answer.at(1) = ( 1. - ksi ) * 0.5;
228  answer.at(2) = ( 1. + ksi ) * 0.5;
229 }
230
231 void
233  const FloatArray &lcoords, const FEICellGeometry &cellgeo)
234 {
235  double l;
236  IntArray edgeNodes;
237  this->computeLocalEdgeMapping(edgeNodes, iedge);
238  l = this->edgeComputeLength(edgeNodes, cellgeo);
239
241  answer.at(1, 1) = -1.0 / l;
242  answer.at(2, 1) = 1.0 / l;
243 }
244
245
246 void
247 FEI3dHexaLin :: edgeEvaldNdxi(FloatArray &answer, int iedge, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
248 {
252 }
253
254
255 void
257  const FloatArray &lcoords, const FEICellGeometry &cellgeo)
258 {
259  IntArray edgeNodes;
260  FloatArray n;
261  this->computeLocalEdgeMapping(edgeNodes, iedge);
262  this->edgeEvalN(n, iedge, lcoords, cellgeo);
263
265  answer.at(1) = ( n.at(1) * cellgeo.giveVertexCoordinates( edgeNodes.at(1) )->at(1) +
266  n.at(2) * cellgeo.giveVertexCoordinates( edgeNodes.at(2) )->at(1) );
267  answer.at(2) = ( n.at(1) * cellgeo.giveVertexCoordinates( edgeNodes.at(1) )->at(2) +
268  n.at(2) * cellgeo.giveVertexCoordinates( edgeNodes.at(2) )->at(2) );
269  answer.at(3) = ( n.at(1) * cellgeo.giveVertexCoordinates( edgeNodes.at(1) )->at(3) +
270  n.at(2) * cellgeo.giveVertexCoordinates( edgeNodes.at(2) )->at(3) );
271 }
272
273
274 double
276 {
277  IntArray edgeNodes;
278  this->computeLocalEdgeMapping(edgeNodes, iedge);
279  return 0.5 * this->edgeComputeLength(edgeNodes, cellgeo);
280 }
281
282
283 void
285 {
286  int aNode = 0, bNode = 0;
287
288  if ( iedge == 1 ) { // edge between nodes 1 2
289  aNode = 1;
290  bNode = 2;
291  } else if ( iedge == 2 ) { // edge between nodes 2 3
292  aNode = 2;
293  bNode = 3;
294  } else if ( iedge == 3 ) { // edge between nodes 3 4
295  aNode = 3;
296  bNode = 4;
297  } else if ( iedge == 4 ) { // edge between nodes 4 1
298  aNode = 4;
299  bNode = 1;
300  } else if ( iedge == 5 ) { // edge between nodes 1 5
301  aNode = 1;
302  bNode = 5;
303  } else if ( iedge == 6 ) { // edge between nodes 2 6
304  aNode = 2;
305  bNode = 6;
306  } else if ( iedge == 7 ) { // edge between nodes 3 7
307  aNode = 3;
308  bNode = 7;
309  } else if ( iedge == 8 ) { // edge between nodes 4 8
310  aNode = 4;
311  bNode = 8;
312  } else if ( iedge == 9 ) { // edge between nodes 5 6
313  aNode = 5;
314  bNode = 6;
315  } else if ( iedge == 10 ) { // edge between nodes 6 7
316  aNode = 6;
317  bNode = 7;
318  } else if ( iedge == 11 ) { // edge between nodes 7 8
319  aNode = 7;
320  bNode = 8;
321  } else if ( iedge == 12 ) { // edge between nodes 8 5
322  aNode = 8;
323  bNode = 5;
324  } else {
325  OOFEM_ERROR("wrong egde number (%d)", iedge);
326  }
327
328  edgeNodes = {aNode, bNode};
329 }
330
331 double
333 {
334  return cellgeo.giveVertexCoordinates( edgeNodes.at(2) )->distance( cellgeo.giveVertexCoordinates( edgeNodes.at(1) ) );
335 }
336
337 void
338 FEI3dHexaLin :: surfaceEvalN(FloatArray &answer, int isurf, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
339 {
340  double ksi, eta;
342
343  ksi = lcoords.at(1);
344  eta = lcoords.at(2);
345
346  answer.at(1) = ( 1. + ksi ) * ( 1. + eta ) * 0.25;
347  answer.at(2) = ( 1. - ksi ) * ( 1. + eta ) * 0.25;
348  answer.at(3) = ( 1. - ksi ) * ( 1. - eta ) * 0.25;
349  answer.at(4) = ( 1. + ksi ) * ( 1. - eta ) * 0.25;
350 }
351
352 void
353 FEI3dHexaLin :: surfaceEvaldNdx(FloatMatrix &answer, int isurf, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
354 {
355  // Note, this must be in correct order, not just the correct nodes, therefore we must use snodes;
356  IntArray snodes;
357  this->computeLocalSurfaceMapping(snodes, isurf);
358
359  FloatArray lcoords_hex;
360
362 #if 1
363  if ( isurf == 1 ) { // surface 1 - nodes 1 4 3 2
364  lcoords_hex = {-lcoords.at(1), -lcoords.at(2), 1};
365  } else if ( isurf == 2 ) { // surface 2 - nodes 5 6 7 8
366  lcoords_hex = {-lcoords.at(2), -lcoords.at(1), -1};
367  } else if ( isurf == 3 ) { // surface 3 - nodes 1 2 6 5
368  lcoords_hex = {-1, -lcoords.at(1), lcoords.at(2)};
369  } else if ( isurf == 4 ) { // surface 4 - nodes 2 3 7 6
370  lcoords_hex = {-lcoords.at(1), 1, lcoords.at(2)};
371  } else if ( isurf == 5 ) { // surface 5 - nodes 3 4 8 7
372  lcoords_hex = {1, lcoords.at(1), lcoords.at(2)};
373  } else if ( isurf == 6 ) { // surface 6 - nodes 4 1 5 8
374  lcoords_hex = {lcoords.at(1), -1, lcoords.at(2)};
375  } else {
376  OOFEM_ERROR("wrong surface number (%d)", isurf);
377  }
378 #else
379  if ( isurf == 1 ) { // surface 1 - nodes 3 4 8 7
381  lcoords_hex = {-1, lcoords.at(1), lcoords.at(2)};
382  } else if ( isurf == 2 ) { // surface 2 - nodes 2 1 5 6
383  lcoords_hex = {1, lcoords.at(1), lcoords.at(2)};
384  } else if ( isurf == 3 ) { // surface 3 - nodes 3 7 6 2
385  lcoords_hex = {lcoords.at(1), -1, lcoords.at(2)};
386  } else if ( isurf == 4 ) { // surface 4 - nodes 4 8 5 1
387  lcoords_hex = {lcoords.at(1), 1, lcoords.at(2)};
388  } else if ( isurf == 5 ) { // surface 5 - nodes 3 2 1 4
389  lcoords_hex = {lcoords.at(1), lcoords.at(2), -1};
390  } else if ( isurf == 6 ) { // surface 6 - nodes 7 6 5 8
391  lcoords_hex = {lcoords.at(1), lcoords.at(2), 1};
392  } else {
393  OOFEM_ERROR("wrong surface number (%d)", isurf);
394  }
395 #endif
396
397  FloatMatrix fullB;
398  this->evaldNdx(fullB, lcoords_hex, cellgeo);
400  for ( int i = 1; i <= snodes.giveSize(); ++i ) {
401  for ( int j = 1; j <= 3; ++j ) {
402  answer.at(i, j) = fullB.at(snodes.at(i), j);
403  }
404  }
405 }
406
407 void
409  const FloatArray &lcoords, const FEICellGeometry &cellgeo)
410 {
411  IntArray nodes;
412  FloatArray n;
413
414  this->computeLocalSurfaceMapping(nodes, iedge);
415  this->surfaceEvalN(n, iedge, lcoords, cellgeo);
416
418  answer.at(1) = n.at(1) * cellgeo.giveVertexCoordinates( nodes.at(1) )->at(1) + n.at(2) * cellgeo.giveVertexCoordinates( nodes.at(2) )->at(1) +
419  n.at(3) * cellgeo.giveVertexCoordinates( nodes.at(3) )->at(1) + n.at(4) * cellgeo.giveVertexCoordinates( nodes.at(4) )->at(1);
420  answer.at(2) = n.at(1) * cellgeo.giveVertexCoordinates( nodes.at(1) )->at(2) + n.at(2) * cellgeo.giveVertexCoordinates( nodes.at(2) )->at(2) +
421  n.at(3) * cellgeo.giveVertexCoordinates( nodes.at(3) )->at(2) + n.at(4) * cellgeo.giveVertexCoordinates( nodes.at(4) )->at(2);
422  answer.at(3) = n.at(1) * cellgeo.giveVertexCoordinates( nodes.at(1) )->at(3) + n.at(2) * cellgeo.giveVertexCoordinates( nodes.at(2) )->at(3) +
423  n.at(3) * cellgeo.giveVertexCoordinates( nodes.at(3) )->at(3) + n.at(4) * cellgeo.giveVertexCoordinates( nodes.at(4) )->at(3);
424 }
425
426 double
427 FEI3dHexaLin :: surfaceEvalNormal(FloatArray &answer, int isurf, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
428 {
429  FloatArray a, b, dNdksi(4), dNdeta(4);
430  double ksi, eta;
431  IntArray snodes;
432
433  this->computeLocalSurfaceMapping(snodes, isurf);
434
435  ksi = lcoords.at(1);
436  eta = lcoords.at(2);
437
438  // No need to divide by 1/4, we'll normalize anyway;
439  dNdksi.at(1) = ( 1. + eta );
440  dNdksi.at(2) = -( 1. + eta );
441  dNdksi.at(3) = -( 1. - eta );
442  dNdksi.at(4) = ( 1. - eta );
443
444  dNdeta.at(1) = ( 1. + ksi );
445  dNdeta.at(2) = ( 1. - ksi );
446  dNdeta.at(3) = -( 1. - ksi );
447  dNdeta.at(4) = -( 1. + ksi );
448
449  for ( int i = 1; i <= 4; ++i ) {
450  a.add( dNdksi.at(i), * cellgeo.giveVertexCoordinates( snodes.at(i) ) );
451  b.add( dNdeta.at(i), * cellgeo.giveVertexCoordinates( snodes.at(i) ) );
452  }
453
456 }
457
458 double
460  const FEICellGeometry &cellgeo)
461 {
462  FloatArray normal;
463  return this->surfaceEvalNormal(normal, isurf, lcoords, cellgeo);
464 }
465
466 void
468 {
469  if ( isurf == 1 ) { // surface 1 - nodes 1 4 3 2
470  surfNodes = {1, 4, 3, 2};
471  } else if ( isurf == 2 ) { // surface 2 - nodes 5 6 7 8
472  surfNodes = {5, 6, 7, 8};
473  } else if ( isurf == 3 ) { // surface 3 - nodes 1 2 6 5
474  surfNodes = {1, 2, 6, 5};
475  } else if ( isurf == 4 ) { // surface 4 - nodes 2 3 7 6
476  surfNodes = {2, 3, 7, 6};
477  } else if ( isurf == 5 ) { // surface 5 - nodes 3 4 8 7
478  surfNodes = {3, 4, 8, 7};
479  } else if ( isurf == 6 ) { // surface 6 - nodes 4 1 5 8
480  surfNodes = {4, 1, 5, 8};
481  } else {
482  OOFEM_ERROR("wrong surface number (%d)", isurf);
483  }
484 }
485
486 void
487 FEI3dHexaLin :: giveJacobianMatrixAt(FloatMatrix &jacobianMatrix, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
488 // Returns the jacobian matrix J (x,y,z)/(ksi,eta,dzeta) of the receiver.
489 {
490  FloatMatrix dNduvw, coords;
491  this->giveLocalDerivative(dNduvw, lcoords);
492  coords.resize(3, 8);
493  for ( int i = 1; i <= 8; i++ ) {
494  coords.setColumn(* cellgeo.giveVertexCoordinates(i), i);
495  }
496  jacobianMatrix.beProductOf(coords, dNduvw);
497 }
498
499 void
501 {
502  double u, v, w;
503  u = lcoords.at(1);
504  v = lcoords.at(2);
505  w = lcoords.at(3);
506
507  dN.resize(8, 3);
508
509  dN.at(1, 1) = -0.125 * ( 1. - v ) * ( 1. + w );
510  dN.at(2, 1) = -0.125 * ( 1. + v ) * ( 1. + w );
511  dN.at(3, 1) = 0.125 * ( 1. + v ) * ( 1. + w );
512  dN.at(4, 1) = 0.125 * ( 1. - v ) * ( 1. + w );
513  dN.at(5, 1) = -0.125 * ( 1. - v ) * ( 1. - w );
514  dN.at(6, 1) = -0.125 * ( 1. + v ) * ( 1. - w );
515  dN.at(7, 1) = 0.125 * ( 1. + v ) * ( 1. - w );
516  dN.at(8, 1) = 0.125 * ( 1. - v ) * ( 1. - w );
517
518  dN.at(1, 2) = -0.125 * ( 1. - u ) * ( 1. + w );
519  dN.at(2, 2) = 0.125 * ( 1. - u ) * ( 1. + w );
520  dN.at(3, 2) = 0.125 * ( 1. + u ) * ( 1. + w );
521  dN.at(4, 2) = -0.125 * ( 1. + u ) * ( 1. + w );
522  dN.at(5, 2) = -0.125 * ( 1. - u ) * ( 1. - w );
523  dN.at(6, 2) = 0.125 * ( 1. - u ) * ( 1. - w );
524  dN.at(7, 2) = 0.125 * ( 1. + u ) * ( 1. - w );
525  dN.at(8, 2) = -0.125 * ( 1. + u ) * ( 1. - w );
526
527  dN.at(1, 3) = 0.125 * ( 1. - u ) * ( 1. - v );
528  dN.at(2, 3) = 0.125 * ( 1. - u ) * ( 1. + v );
529  dN.at(3, 3) = 0.125 * ( 1. + u ) * ( 1. + v );
530  dN.at(4, 3) = 0.125 * ( 1. + u ) * ( 1. - v );
531  dN.at(5, 3) = -0.125 * ( 1. - u ) * ( 1. - v );
532  dN.at(6, 3) = -0.125 * ( 1. - u ) * ( 1. + v );
533  dN.at(7, 3) = -0.125 * ( 1. + u ) * ( 1. + v );
534  dN.at(8, 3) = -0.125 * ( 1. + u ) * ( 1. - v );
535 }
536
537 double
539 {
540  IntArray fNodes;
541  this->computeLocalSurfaceMapping(fNodes, iEdge);
542
543  const FloatArray &c1 = * cellgeo.giveVertexCoordinates( fNodes.at(1) );
544  const FloatArray &c2 = * cellgeo.giveVertexCoordinates( fNodes.at(2) );
545  const FloatArray &c3 = * cellgeo.giveVertexCoordinates( fNodes.at(3) );
546  const FloatArray &c4 = * cellgeo.giveVertexCoordinates( fNodes.at(4) );
547
548  return (
549  c4(2) * ( c1(1) * ( -c2(0) - c3(0) ) + c2(1) * ( c1(0) - c3(0) ) + c3(1) * ( c1(0) + c2(0) ) ) +
550  c3(2) * ( c1(1) * ( -c2(0) + c4(0) ) + c2(1) * ( c1(0) + c4(0) ) + c4(1) * ( -c1(0) - c2(0) ) ) +
551  c2(2) * ( c1(1) * ( c3(0) + c4(0) ) + c3(1) * ( -c1(0) - c4(0) ) + c4(1) * ( -c1(0) + c3(0) ) ) +
552  c1(2) * ( c2(1) * ( -c3(0) - c4(0) ) + c3(1) * ( c2(0) - c4(0) ) + c4(1) * ( c2(0) + c3(0) ) ) ) * 0.25;
553 }
554
557 {
558  IntegrationRule *iRule = new GaussIntegrationRule(1, NULL);
559  int points = iRule->getRequiredNumberOfIntegrationPoints(_Cube, order + 6);
560  iRule->SetUpPointsOnCube(points, _Unknown);
561  return iRule;
562 }
563
566 {
567  IntegrationRule *iRule = new GaussIntegrationRule(1, NULL);
568  int points = iRule->getRequiredNumberOfIntegrationPoints(_Square, order + 2);
569  iRule->SetUpPointsOnSquare(points, _Unknown);
570  return iRule;
571 }
572 } // end namespace oofem
double giveDeterminant() const
Returns the trace (sum of diagonal components) of the receiver.
Definition: floatmatrix.C:1408
virtual void computeLocalEdgeMapping(IntArray &edgeNodes, int iedge)
Definition: fei3dhexalin.C:284
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...
Definition: fei3dhexalin.C:232
void subtract(const FloatArray &src)
Definition: floatarray.C:258
void beVectorProductOf(const FloatArray &v1, const FloatArray &v2)
Computes vector product (or cross product) of vectors given as parameters, , and stores the result in...
Definition: floatarray.C:415
virtual int global2local(FloatArray &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates local coordinates from given global ones.
Definition: fei3dhexalin.C:94
bool solveForRhs(const FloatArray &b, FloatArray &answer, bool transpose=false)
Solves the system of linear equations .
Definition: floatmatrix.C:1112
virtual void edgeLocal2global(FloatArray &answer, int iedge, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates edge global coordinates from given local ones.
Definition: fei3dhexalin.C:256
virtual IntegrationRule * giveBoundaryIntegrationRule(int order, int boundary)
Sets up a suitable integration rule for integrating over the requested boundary.
Definition: fei3dhexalin.C:565
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.
Definition: fei3dhexalin.C:338
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
void clear()
Definition: floatarray.h:206
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: fei3dhexalin.C:63
virtual void surfaceEvaldNdx(FloatMatrix &answer, int isurf, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates the matrix of derivatives of edge interpolation functions (shape functions) at given point...
Definition: fei3dhexalin.C:353
virtual int SetUpPointsOnSquare(int, MaterialMode mode)
Sets up receiver&#39;s integration points on unit square integration domain.
virtual double evalNXIntegral(int iEdge, const FEICellGeometry &cellgeo)
Computes the integral .
Definition: fei3dhexalin.C:538
Class implementing an array of integers.
Definition: intarray.h:61
int & at(int i)
Coefficient access function.
Definition: intarray.h:103
virtual void local2global(FloatArray &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates global coordinates from given local ones.
Definition: fei3dhexalin.C:80
virtual void evalN(FloatArray &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates the array of interpolation functions (shape functions) at given point.
Definition: fei3dhexalin.C:43
Abstract base class representing integration rule.
virtual void surfaceLocal2global(FloatArray &answer, int isurf, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates edge global coordinates from given local ones.
Definition: fei3dhexalin.C:408
#define OOFEM_ERROR(...)
Definition: error.h:61
#define POINT_TOL
Definition: fei3dhexalin.C:91
void giveLocalDerivative(FloatMatrix &dN, const FloatArray &lcoords)
Definition: fei3dhexalin.C:500
virtual int SetUpPointsOnCube(int, MaterialMode mode)
Sets up receiver&#39;s integration points on unit cube integration domain.
double at(int i, int j) const
Coefficient access function.
Definition: floatmatrix.h:176
Class representing vector of real numbers.
Definition: floatarray.h:82
Implementation of matrix containing floating point numbers.
Definition: floatmatrix.h:94
virtual IntegrationRule * giveIntegrationRule(int order)
Sets up a suitable integration rule for numerical integrating over volume.
Definition: fei3dhexalin.C:556
virtual double edgeGiveTransformationJacobian(int iedge, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates the edge jacobian of transformation between local and global coordinates.
Definition: fei3dhexalin.C:275
virtual void edgeEvaldNdxi(FloatArray &answer, int iedge, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates the matrix of derivatives of edge interpolation functions (shape functions) at given point...
Definition: fei3dhexalin.C:247
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.
double edgeComputeLength(IntArray &edgeNodes, const FEICellGeometry &cellgeo)
Definition: fei3dhexalin.C:332
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 void computeLocalSurfaceMapping(IntArray &edgeNodes, int iedge)
Definition: fei3dhexalin.C:467
void beProductOf(const FloatMatrix &a, const FloatMatrix &b)
Assigns to the receiver product of .
Definition: floatmatrix.C:337
int giveSize() const
Definition: intarray.h:203
virtual double surfaceEvalNormal(FloatArray &answer, int isurf, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates the normal out of the surface at given point.
Definition: fei3dhexalin.C:427
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 surfaceGiveTransformationJacobian(int isurf, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates the edge jacobian of transformation between local and global coordinates.
Definition: fei3dhexalin.C:459
Definition: floatarray.C:156
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.
Definition: fei3dhexalin.C:222
virtual void giveJacobianMatrixAt(FloatMatrix &jacobianMatrix, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Gives the jacobian matrix at the local coordinates.
Definition: fei3dhexalin.C:487