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 "intarray.h"
37 #include "floatarray.h"
38 #include "floatmatrix.h"
39 #include "gaussintegrationrule.h"
40
41 namespace oofem {
42 void
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.125 * ( 1.0 - u ) * ( 1.0 - v ) * ( 1.0 + w ) * ( -u - v + w - 2.0 );
53  answer.at(2) = 0.125 * ( 1.0 - u ) * ( 1.0 + v ) * ( 1.0 + w ) * ( -u + v + w - 2.0 );
54  answer.at(3) = 0.125 * ( 1.0 + u ) * ( 1.0 + v ) * ( 1.0 + w ) * ( u + v + w - 2.0 );
55  answer.at(4) = 0.125 * ( 1.0 + u ) * ( 1.0 - v ) * ( 1.0 + w ) * ( u - v + w - 2.0 );
56  answer.at(5) = 0.125 * ( 1.0 - u ) * ( 1.0 - v ) * ( 1.0 - w ) * ( -u - v - w - 2.0 );
57  answer.at(6) = 0.125 * ( 1.0 - u ) * ( 1.0 + v ) * ( 1.0 - w ) * ( -u + v - w - 2.0 );
58  answer.at(7) = 0.125 * ( 1.0 + u ) * ( 1.0 + v ) * ( 1.0 - w ) * ( u + v - w - 2.0 );
59  answer.at(8) = 0.125 * ( 1.0 + u ) * ( 1.0 - v ) * ( 1.0 - w ) * ( u - v - w - 2.0 );
60
61  answer.at(9) = 0.25 * ( 1.0 - v * v ) * ( 1.0 - u ) * ( 1.0 + w );
62  answer.at(10) = 0.25 * ( 1.0 - u * u ) * ( 1.0 + v ) * ( 1.0 + w );
63  answer.at(11) = 0.25 * ( 1.0 - v * v ) * ( 1.0 + u ) * ( 1.0 + w );
64  answer.at(12) = 0.25 * ( 1.0 - u * u ) * ( 1.0 - v ) * ( 1.0 + w );
65
66  answer.at(13) = 0.25 * ( 1.0 - v * v ) * ( 1.0 - u ) * ( 1.0 - w );
67  answer.at(14) = 0.25 * ( 1.0 - u * u ) * ( 1.0 + v ) * ( 1.0 - w );
68  answer.at(15) = 0.25 * ( 1.0 - v * v ) * ( 1.0 + u ) * ( 1.0 - w );
69  answer.at(16) = 0.25 * ( 1.0 - u * u ) * ( 1.0 - v ) * ( 1.0 - w );
70
71  answer.at(17) = 0.25 * ( 1.0 - u ) * ( 1.0 - v ) * ( 1.0 - w * w );
72  answer.at(18) = 0.25 * ( 1.0 - u ) * ( 1.0 + v ) * ( 1.0 - w * w );
73  answer.at(19) = 0.25 * ( 1.0 + u ) * ( 1.0 + v ) * ( 1.0 - w * w );
74  answer.at(20) = 0.25 * ( 1.0 + u ) * ( 1.0 - v ) * ( 1.0 - w * w );
75 }
76
77 double
79 {
80  FloatMatrix jacobianMatrix, inv, dNduvw, coords;
81  this->giveLocalDerivative(dNduvw, lcoords);
82  coords.resize( 3, dNduvw.giveNumberOfRows() );
83  for ( int i = 1; i <= dNduvw.giveNumberOfRows(); i++ ) {
84  coords.setColumn(* cellgeo.giveVertexCoordinates(i), i);
85  }
86  jacobianMatrix.beProductOf(coords, dNduvw);
87  inv.beInverseOf(jacobianMatrix);
88
90  return jacobianMatrix.giveDeterminant();
91 }
92
93 void
95 {
96  FloatArray n;
97
98  this->evalN(n, lcoords, cellgeo);
100  for ( int i = 1; i <= n.giveSize(); i++ ) {
102  }
103 }
104
106 {
107  const FloatArray *n1 = cellgeo.giveVertexCoordinates(1);
108  const FloatArray *n2 = cellgeo.giveVertexCoordinates(7);
110  return n1->distance(n2);
111 }
112
113
114 #define POINT_TOL 1.e-3
115
116 int
118 {
119  FloatArray res, delta, guess;
120  FloatMatrix jac;
121  double convergence_limit, error = 0.0;
122
123  // find a suitable convergence limit
124  convergence_limit = 1e-6 * this->giveCharacteristicLength(cellgeo);
125
126  // setup initial guess
129
130  // apply Newton-Raphson to solve the problem
131  for ( int nite = 0; nite < 40; nite++ ) {
132  // compute the residual
134  res.beDifferenceOf(gcoords, guess);
135
136  // check for convergence
137  error = res.computeNorm();
138  if ( error < convergence_limit ) {
139  break;
140  }
141
142  // compute the corrections
144  jac.solveForRhs(res, delta);
145
146  // update guess
148  }
149  if ( error > convergence_limit ) { // Imperfect, could give false negatives.
150  //OOFEM_ERROR("no convergence after 10 iterations");
152  return false;
153  }
154
155  // check limits for each local coordinate [-1,1] for quadrilaterals. (different for other elements, typically [0,1]).
156  bool inside = true;
157  for ( int i = 1; i <= answer.giveSize(); i++ ) {
158  if ( answer.at(i) < ( -1. - POINT_TOL ) ) {
160  inside = false;
161  } else if ( answer.at(i) > ( 1. + POINT_TOL ) ) {
163  inside = false;
164  }
165  }
166
167  return inside;
168 }
169
170 void FEI3dHexaQuad :: edgeEvalN(FloatArray &answer, int iedge, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
171 {
172  double u = lcoords.at(1);
174  answer.at(1) = 0.5 * ( u - 1. ) * u;
175  answer.at(2) = 0.5 * ( u + 1. ) * u;
176  answer.at(3) = 1. - u * u;
177 }
178
179 void FEI3dHexaQuad :: edgeLocal2global(FloatArray &answer, int iedge, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
180 {
181  IntArray eNodes;
182  double u = lcoords.at(1);
183  this->computeLocalEdgeMapping(eNodes, iedge);
185  answer.add( 0.5 * ( u - 1. ) * u, * cellgeo.giveVertexCoordinates( eNodes.at(1) ) );
186  answer.add( 0.5 * ( u - 1. ) * u, * cellgeo.giveVertexCoordinates( eNodes.at(2) ) );
187  answer.add( 1. - u * u, * cellgeo.giveVertexCoordinates( eNodes.at(3) ) );
188 }
189
190 void FEI3dHexaQuad :: edgeEvaldNdx(FloatMatrix &answer, int iedge, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
191 {
193  IntArray eNodes;
194  FloatArray dNdu;
195  double u = lcoords.at(1);
196  this->computeLocalEdgeMapping(eNodes, iedge);
197  dNdu.add( u - 0.5, * cellgeo.giveVertexCoordinates( eNodes.at(1) ) );
198  dNdu.add( u + 0.5, * cellgeo.giveVertexCoordinates( eNodes.at(2) ) );
199  dNdu.add( -2. * u, * cellgeo.giveVertexCoordinates( eNodes.at(3) ) );
200  // Why matrix output?
203 }
204
205 double FEI3dHexaQuad :: edgeGiveTransformationJacobian(int iedge, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
206 {
207  IntArray eNodes;
208  FloatArray dNdu;
209  double u = lcoords.at(1);
210  this->computeLocalEdgeMapping(eNodes, iedge);
211  dNdu.add( u - 0.5, * cellgeo.giveVertexCoordinates( eNodes.at(1) ) );
212  dNdu.add( u + 0.5, * cellgeo.giveVertexCoordinates( eNodes.at(2) ) );
213  dNdu.add( -2. * u, * cellgeo.giveVertexCoordinates( eNodes.at(3) ) );
214  return dNdu.computeNorm();
215 }
216
217 void
219 {
220  if ( iedge == 1 ) {
221  edgeNodes = { 1, 2, 9};
222  } else if ( iedge == 2 ) {
223  edgeNodes = { 2, 3, 10};
224  } else if ( iedge == 3 ) {
225  edgeNodes = { 3, 4, 11};
226  } else if ( iedge == 4 ) {
227  edgeNodes = { 4, 1, 12};
228  } else if ( iedge == 5 ) {
229  edgeNodes = { 5, 6, 13};
230  } else if ( iedge == 6 ) {
231  edgeNodes = { 6, 7, 14};
232  } else if ( iedge == 7 ) {
233  edgeNodes = { 7, 8, 15};
234  } else if ( iedge == 8 ) {
235  edgeNodes = { 8, 5, 16};
236  } else if ( iedge == 9 ) {
237  edgeNodes = { 1, 5, 17};
238  } else if ( iedge == 10 ) {
239  edgeNodes = { 2, 6, 18};
240  } else if ( iedge == 11 ) {
241  edgeNodes = { 3, 7, 19};
242  } else if ( iedge == 12 ) {
243  edgeNodes = { 4, 8, 20};
244  } else {
245  OOFEM_ERROR("wrong edge number (%d)", iedge);
246  }
247 }
248
249 void
250 FEI3dHexaQuad :: surfaceEvalN(FloatArray &answer, int isurf, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
251 {
252  double ksi, eta;
254
255  ksi = lcoords.at(1);
256  eta = lcoords.at(2);
257
258  answer.at(1) = ( 1. + ksi ) * ( 1. + eta ) * 0.25 * ( ksi + eta - 1. );
259  answer.at(2) = ( 1. - ksi ) * ( 1. + eta ) * 0.25 * ( -ksi + eta - 1. );
260  answer.at(3) = ( 1. - ksi ) * ( 1. - eta ) * 0.25 * ( -ksi - eta - 1. );
261  answer.at(4) = ( 1. + ksi ) * ( 1. - eta ) * 0.25 * ( ksi - eta - 1. );
262  answer.at(5) = 0.5 * ( 1. - ksi * ksi ) * ( 1. + eta );
263  answer.at(6) = 0.5 * ( 1. - ksi ) * ( 1. - eta * eta );
264  answer.at(7) = 0.5 * ( 1. - ksi * ksi ) * ( 1. - eta );
265  answer.at(8) = 0.5 * ( 1. + ksi ) * ( 1. - eta * eta );
266 }
267
268 void
269 FEI3dHexaQuad :: surfaceEvaldNdx(FloatMatrix &answer, int isurf, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
270 {
271  // Note, this must be in correct order, not just the correct nodes, therefore we must use snodes;
272  IntArray snodes;
273  this->computeLocalSurfaceMapping(snodes, isurf);
274
275  FloatArray lcoords_hex;
276
278 #if 1
279  if ( isurf == 1 ) { // surface 1 - nodes 1 4 3 2
280  lcoords_hex = {-lcoords.at(1), -lcoords.at(2), 1};
281  } else if ( isurf == 2 ) { // surface 2 - nodes 5 6 7 8
282  lcoords_hex = {-lcoords.at(2), -lcoords.at(1), -1};
283  } else if ( isurf == 3 ) { // surface 3 - nodes 1 2 6 5
284  lcoords_hex = {-1, -lcoords.at(1), lcoords.at(2)};
285  } else if ( isurf == 4 ) { // surface 4 - nodes 2 3 7 6
286  lcoords_hex = {-lcoords.at(1), 1, lcoords.at(2)};
287  } else if ( isurf == 5 ) { // surface 5 - nodes 3 4 8 7
288  lcoords_hex = {1, lcoords.at(1), lcoords.at(2)};
289  } else if ( isurf == 6 ) { // surface 6 - nodes 4 1 5 8
290  lcoords_hex = {lcoords.at(1), -1, lcoords.at(2)};
291  } else {
292  OOFEM_ERROR("wrong surface number (%d)", isurf);
293  }
294 #else
295  if ( isurf == 1 ) { // surface 1 - nodes 3 4 8 7
297  lcoords_hex = {-1, lcoords.at(1), lcoords.at(2)};
298  } else if ( isurf == 2 ) { // surface 2 - nodes 2 1 5 6
299  lcoords_hex = {1, lcoords.at(1), lcoords.at(2)};
300  } else if ( isurf == 3 ) { // surface 3 - nodes 3 7 6 2
301  lcoords_hex = {lcoords.at(1), -1, lcoords.at(2)};
302  } else if ( isurf == 4 ) { // surface 4 - nodes 4 8 5 1
303  lcoords_hex = {lcoords.at(1), 1, lcoords.at(2)};
304  } else if ( isurf == 5 ) { // surface 5 - nodes 3 2 1 4
305  lcoords_hex = {lcoords.at(1), lcoords.at(2), -1};
306  } else if ( isurf == 6 ) { // surface 6 - nodes 7 6 5 8
307  lcoords_hex = {lcoords.at(1), lcoords.at(2), 1};
308  } else {
309  OOFEM_ERROR("wrong surface number (%d)", isurf);
310  }
311 #endif
312
313  FloatMatrix fullB;
314  this->evaldNdx(fullB, lcoords_hex, cellgeo);
316  for ( int i = 1; i <= snodes.giveSize(); ++i ) {
317  for ( int j = 1; j <= 3; ++j ) {
318  answer.at(i, j) = fullB.at(snodes.at(i), j);
319  }
320  }
321 }
322
323
324 double
325 FEI3dHexaQuad :: surfaceEvalNormal(FloatArray &answer, int isurf, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
326 {
327  FloatArray a, b, dNdksi(8), dNdeta(8);
328  double ksi, eta;
329  IntArray snodes;
330
331  this->computeLocalSurfaceMapping(snodes, isurf);
332
333  ksi = lcoords.at(1);
334  eta = lcoords.at(2);
335
336  dNdksi.at(1) = 0.25 * ( 1. + eta ) * ( 2.0 * ksi + eta );
337  dNdksi.at(2) = -0.25 * ( 1. + eta ) * ( -2.0 * ksi + eta );
338  dNdksi.at(3) = -0.25 * ( 1. - eta ) * ( -2.0 * ksi - eta );
339  dNdksi.at(4) = 0.25 * ( 1. - eta ) * ( 2.0 * ksi - eta );
340  dNdksi.at(5) = -ksi * ( 1. + eta );
341  dNdksi.at(6) = -0.5 * ( 1. - eta * eta );
342  dNdksi.at(7) = -ksi * ( 1. - eta );
343  dNdksi.at(8) = 0.5 * ( 1. - eta * eta );
344
345  dNdeta.at(1) = 0.25 * ( 1. + ksi ) * ( 2.0 * eta + ksi );
346  dNdeta.at(2) = 0.25 * ( 1. - ksi ) * ( 2.0 * eta - ksi );
347  dNdeta.at(3) = -0.25 * ( 1. - ksi ) * ( -2.0 * eta - ksi );
348  dNdeta.at(4) = -0.25 * ( 1. + ksi ) * ( -2.0 * eta + ksi );
349  dNdeta.at(5) = 0.5 * ( 1. - ksi * ksi );
350  dNdeta.at(6) = -eta * ( 1. - ksi );
351  dNdeta.at(7) = -0.5 * ( 1. - ksi * ksi );
352  dNdeta.at(8) = -eta * ( 1. + ksi );
353
354  for ( int i = 1; i <= 8; ++i ) {
355  a.add( dNdksi.at(i), * cellgeo.giveVertexCoordinates( snodes.at(i) ) );
356  b.add( dNdeta.at(i), * cellgeo.giveVertexCoordinates( snodes.at(i) ) );
357  }
358
361 }
362
363 void
365  const FloatArray &lcoords, const FEICellGeometry &cellgeo)
366 {
367  IntArray nodes;
368  FloatArray n;
369
370  this->computeLocalSurfaceMapping(nodes, isurf);
371
372  this->surfaceEvalN(n, isurf, lcoords, cellgeo);
373
375  for ( int i = 1; i <= n.giveSize(); ++i ) {
377  }
378 }
379
380
381 double
383  const FEICellGeometry &cellgeo)
384 {
385  FloatArray normal;
386  return this->surfaceEvalNormal(normal, isurf, lcoords, cellgeo);
387 }
388
389
390 void
392 {
393  nodes.resize(8);
394
395  // the actual numbering has a positive normal pointing outwards from the element - (LSpace compatible)
396
397  if ( isurf == 1 ) {
398  nodes = { 3, 2, 1, 4, 10, 9, 12, 11};
399  } else if ( isurf == 2 ) {
400  nodes = { 7, 8, 5, 6, 15, 16, 13, 14};
401  } else if ( isurf == 3 ) {
402  nodes = { 2, 6, 5, 1, 18, 13, 17, 9};
403  } else if ( isurf == 4 ) {
404  nodes = { 3, 7, 6, 2, 19, 14, 18, 10};
405  } else if ( isurf == 5 ) {
406  nodes = { 3, 4, 8, 7, 11, 20, 15, 19};
407  } else if ( isurf == 6 ) {
408  nodes = { 4, 1, 5, 8, 12, 17, 16, 20};
409  } else {
410  OOFEM_ERROR("wrong surface number (%d)", isurf);
411  }
412
413 #if 0
414  // this commented numbering is symmetrical with respect to local coordinate axes
415
416  if ( isurf == 1 ) { // surface 1 - nodes 3 2 1 4 10 9 12 11
417  nodes.at(1) = 3;
418  nodes.at(2) = 2;
419  nodes.at(3) = 1;
420  nodes.at(4) = 4;
421  nodes.at(5) = 10;
422  nodes.at(6) = 9;
423  nodes.at(7) = 12;
424  nodes.at(8) = 11;
425  } else if ( iSurf == 2 ) { // surface 2 - nodes 7 6 5 8 14 13 16 15
426  nodes.at(1) = 7;
427  nodes.at(2) = 6;
428  nodes.at(3) = 5;
429  nodes.at(4) = 8;
430  nodes.at(5) = 14;
431  nodes.at(6) = 13;
432  nodes.at(7) = 16;
433  nodes.at(8) = 15;
434  } else if ( iSurf == 3 ) { // surface 3 - nodes 2 1 5 6 9 17 13 18
435  nodes.at(1) = 2;
436  nodes.at(2) = 1;
437  nodes.at(3) = 5;
438  nodes.at(4) = 6;
439  nodes.at(5) = 9;
440  nodes.at(6) = 17;
441  nodes.at(7) = 13;
442  nodes.at(8) = 18;
443  } else if ( isurf == 4 ) { // surface 4 - nodes 3 7 6 2 19 14 18 10
444  nodes.at(1) = 3;
445  nodes.at(2) = 7;
446  nodes.at(3) = 6;
447  nodes.at(4) = 2;
448  nodes.at(5) = 19;
449  nodes.at(6) = 14;
450  nodes.at(7) = 18;
451  nodes.at(8) = 10;
452  } else if ( isurf == 5 ) { // surface 5 - nodes 3 4 8 7 11 20 15 19
453  nodes.at(1) = 3;
454  nodes.at(2) = 4;
455  nodes.at(3) = 8;
456  nodes.at(4) = 7;
457  nodes.at(5) = 11;
458  nodes.at(6) = 20;
459  nodes.at(7) = 15;
460  nodes.at(8) = 19;
461  } else if ( iSurf == 6 ) { // surface 6 - nodes 4 8 5 1 20 16 17 12
462  nodes.at(1) = 4;
463  nodes.at(2) = 8;
464  nodes.at(3) = 5;
465  nodes.at(4) = 1;
466  nodes.at(5) = 20;
467  nodes.at(6) = 16;
468  nodes.at(7) = 17;
469  nodes.at(8) = 12;
470  } else {
471  OOFEM_ERROR("wrong surface number (%d)", isurf);
472  }
473 #endif
474 }
475
476
477 void
478 FEI3dHexaQuad :: giveJacobianMatrixAt(FloatMatrix &jacobianMatrix, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
479 // Returns the jacobian matrix J (x,y,z)/(ksi,eta,dzeta) of the receiver.
480 {
481  FloatMatrix dNduvw, coords;
482  this->giveLocalDerivative(dNduvw, lcoords);
483  coords.resize( 3, dNduvw.giveNumberOfRows() );
484  for ( int i = 1; i <= dNduvw.giveNumberOfRows(); i++ ) {
485  coords.setColumn(* cellgeo.giveVertexCoordinates(i), i);
486  }
487  jacobianMatrix.beProductOf(coords, dNduvw);
488 }
489
490
491 void
493 {
494  double u, v, w;
495  u = lcoords.at(1);
496  v = lcoords.at(2);
497  w = lcoords.at(3);
498
499  dN.resize(20, 3);
500
501  dN.at(1, 1) = 0.125 * ( 1.0 - v ) * ( 1.0 + w ) * ( 2.0 * u + v - w + 1.0 );
502  dN.at(2, 1) = 0.125 * ( 1.0 + v ) * ( 1.0 + w ) * ( 2.0 * u - v - w + 1.0 );
503  dN.at(3, 1) = 0.125 * ( 1.0 + v ) * ( 1.0 + w ) * ( 2.0 * u + v + w - 1.0 );
504  dN.at(4, 1) = 0.125 * ( 1.0 - v ) * ( 1.0 + w ) * ( 2.0 * u - v + w - 1.0 );
505  dN.at(5, 1) = 0.125 * ( 1.0 - v ) * ( 1.0 - w ) * ( 2.0 * u + v + w + 1.0 );
506  dN.at(6, 1) = 0.125 * ( 1.0 + v ) * ( 1.0 - w ) * ( 2.0 * u - v + w + 1.0 );
507  dN.at(7, 1) = 0.125 * ( 1.0 + v ) * ( 1.0 - w ) * ( 2.0 * u + v - w - 1.0 );
508  dN.at(8, 1) = 0.125 * ( 1.0 - v ) * ( 1.0 - w ) * ( 2.0 * u - v - w - 1.0 );
509  dN.at(9, 1) = -0.25 * ( 1.0 - v * v ) * ( 1.0 + w );
510  dN.at(10, 1) = -0.5 * u * ( 1.0 + v ) * ( 1.0 + w );
511  dN.at(11, 1) = 0.25 * ( 1.0 - v * v ) * ( 1.0 + w );
512  dN.at(12, 1) = -0.5 * u * ( 1.0 - v ) * ( 1.0 + w );
513  dN.at(13, 1) = -0.25 * ( 1.0 - v * v ) * ( 1.0 - w );
514  dN.at(14, 1) = -0.5 * u * ( 1.0 + v ) * ( 1.0 - w );
515  dN.at(15, 1) = 0.25 * ( 1.0 - v * v ) * ( 1.0 - w );
516  dN.at(16, 1) = -0.5 * u * ( 1.0 - v ) * ( 1.0 - w );
517  dN.at(17, 1) = -0.25 * ( 1.0 - v ) * ( 1.0 - w * w );
518  dN.at(18, 1) = -0.25 * ( 1.0 + v ) * ( 1.0 - w * w );
519  dN.at(19, 1) = 0.25 * ( 1.0 + v ) * ( 1.0 - w * w );
520  dN.at(20, 1) = 0.25 * ( 1.0 - v ) * ( 1.0 - w * w );
521
522  dN.at(1, 2) = 0.125 * ( 1.0 - u ) * ( 1.0 + w ) * ( 2.0 * v + u - w + 1.0 );
523  dN.at(2, 2) = 0.125 * ( 1.0 - u ) * ( 1.0 + w ) * ( 2.0 * v - u + w - 1.0 );
524  dN.at(3, 2) = 0.125 * ( 1.0 + u ) * ( 1.0 + w ) * ( 2.0 * v + u + w - 1.0 );
525  dN.at(4, 2) = 0.125 * ( 1.0 + u ) * ( 1.0 + w ) * ( 2.0 * v - u - w + 1.0 );
526  dN.at(5, 2) = 0.125 * ( 1.0 - u ) * ( 1.0 - w ) * ( 2.0 * v + u + w + 1.0 );
527  dN.at(6, 2) = 0.125 * ( 1.0 - u ) * ( 1.0 - w ) * ( 2.0 * v - u - w - 1.0 );
528  dN.at(7, 2) = 0.125 * ( 1.0 + u ) * ( 1.0 - w ) * ( 2.0 * v + u - w - 1.0 );
529  dN.at(8, 2) = 0.125 * ( 1.0 + u ) * ( 1.0 - w ) * ( 2.0 * v - u + w + 1.0 );
530  dN.at(9, 2) = -0.5 * v * ( 1.0 - u ) * ( 1.0 + w );
531  dN.at(10, 2) = 0.25 * ( 1.0 - u * u ) * ( 1.0 + w );
532  dN.at(11, 2) = -0.5 * v * ( 1.0 + u ) * ( 1.0 + w );
533  dN.at(12, 2) = -0.25 * ( 1.0 - u * u ) * ( 1.0 + w );
534  dN.at(13, 2) = -0.5 * v * ( 1.0 - u ) * ( 1.0 - w );
535  dN.at(14, 2) = 0.25 * ( 1.0 - u * u ) * ( 1.0 - w );
536  dN.at(15, 2) = -0.5 * v * ( 1.0 + u ) * ( 1.0 - w );
537  dN.at(16, 2) = -0.25 * ( 1.0 - u * u ) * ( 1.0 - w );
538  dN.at(17, 2) = -0.25 * ( 1.0 - u ) * ( 1.0 - w * w );
539  dN.at(18, 2) = 0.25 * ( 1.0 - u ) * ( 1.0 - w * w );
540  dN.at(19, 2) = 0.25 * ( 1.0 + u ) * ( 1.0 - w * w );
541  dN.at(20, 2) = -0.25 * ( 1.0 + u ) * ( 1.0 - w * w );
542
543  dN.at(1, 3) = 0.125 * ( 1.0 - u ) * ( 1.0 - v ) * ( 2.0 * w - u - v - 1.0 );
544  dN.at(2, 3) = 0.125 * ( 1.0 - u ) * ( 1.0 + v ) * ( 2.0 * w - u + v - 1.0 );
545  dN.at(3, 3) = 0.125 * ( 1.0 + u ) * ( 1.0 + v ) * ( 2.0 * w + u + v - 1.0 );
546  dN.at(4, 3) = 0.125 * ( 1.0 + u ) * ( 1.0 - v ) * ( 2.0 * w + u - v - 1.0 );
547  dN.at(5, 3) = 0.125 * ( 1.0 - u ) * ( 1.0 - v ) * ( 2.0 * w + u + v + 1.0 );
548  dN.at(6, 3) = 0.125 * ( 1.0 - u ) * ( 1.0 + v ) * ( 2.0 * w + u - v + 1.0 );
549  dN.at(7, 3) = 0.125 * ( 1.0 + u ) * ( 1.0 + v ) * ( 2.0 * w - u - v + 1.0 );
550  dN.at(8, 3) = 0.125 * ( 1.0 + u ) * ( 1.0 - v ) * ( 2.0 * w - u + v + 1.0 );
551  dN.at(9, 3) = 0.25 * ( 1.0 - v * v ) * ( 1.0 - u );
552  dN.at(10, 3) = 0.25 * ( 1.0 - u * u ) * ( 1.0 + v );
553  dN.at(11, 3) = 0.25 * ( 1.0 - v * v ) * ( 1.0 + u );
554  dN.at(12, 3) = 0.25 * ( 1.0 - u * u ) * ( 1.0 - v );
555  dN.at(13, 3) = -0.25 * ( 1.0 - v * v ) * ( 1.0 - u );
556  dN.at(14, 3) = -0.25 * ( 1.0 - u * u ) * ( 1.0 + v );
557  dN.at(15, 3) = -0.25 * ( 1.0 - v * v ) * ( 1.0 + u );
558  dN.at(16, 3) = -0.25 * ( 1.0 - u * u ) * ( 1.0 - v );
559  dN.at(17, 3) = -0.5 * ( 1.0 - u ) * ( 1.0 - v ) * w;
560  dN.at(18, 3) = -0.5 * ( 1.0 - u ) * ( 1.0 + v ) * w;
561  dN.at(19, 3) = -0.5 * ( 1.0 + u ) * ( 1.0 + v ) * w;
562  dN.at(20, 3) = -0.5 * ( 1.0 + u ) * ( 1.0 - v ) * w;
563 }
564
565
566 double
568 {
569  IntArray fNodes;
570  this->computeLocalSurfaceMapping(fNodes, iEdge);
571
572  const FloatArray &c1 = * cellgeo.giveVertexCoordinates( fNodes.at(1) );
573  const FloatArray &c2 = * cellgeo.giveVertexCoordinates( fNodes.at(2) );
574  const FloatArray &c3 = * cellgeo.giveVertexCoordinates( fNodes.at(3) );
575  const FloatArray &c4 = * cellgeo.giveVertexCoordinates( fNodes.at(4) );
576  const FloatArray &c5 = * cellgeo.giveVertexCoordinates( fNodes.at(5) );
577  const FloatArray &c6 = * cellgeo.giveVertexCoordinates( fNodes.at(6) );
578  const FloatArray &c7 = * cellgeo.giveVertexCoordinates( fNodes.at(7) );
579  const FloatArray &c8 = * cellgeo.giveVertexCoordinates( fNodes.at(8) );
580
581  // Generated with Mathematica (rather unwieldy expression, tried to simplify it as good as possible, but it could probably be better)
582  return (
583  c1(2) * ( c2(1) * ( -3 * c3(0) - 3 * c4(0) - 12 * c5(0) + 14 * c6(0) + 14 * c8(0) ) +
584  c3(1) * ( 3 * c2(0) - 3 * c4(0) - 6 * c5(0) - 6 * c6(0) + 6 * c7(0) + 6 * c8(0) ) +
585  c4(1) * ( 3 * c2(0) + 3 * c3(0) - 14 * c5(0) - 14 * c7(0) + 12 * c8(0) ) +
586  c5(1) * ( 12 * c2(0) + 6 * c3(0) + 14 * c4(0) - 4 * c6(0) - 8 * c7(0) - 60 * c8(0) ) +
587  c6(1) * ( -14 * c2(0) + 6 * c3(0) + 4 * c5(0) + 12 * c7(0) - 8 * c8(0) ) +
588  c7(1) * ( -6 * c3(0) + 14 * c4(0) + 8 * c5(0) - 12 * c6(0) - 4 * c8(0) ) +
589  c8(1) * ( -14 * c2(0) - 6 * c3(0) - 12 * c4(0) + 60 * c5(0) + 8 * c6(0) + 4 * c7(0) ) ) +
590  c2(2) * ( c1(1) * ( 3 * c3(0) + 3 * c4(0) + 12 * c5(0) - 14 * c6(0) - 14 * c8(0) ) +
591  c3(1) * ( -3 * c1(0) - 3 * c4(0) + 14 * c5(0) - 12 * c6(0) + 14 * c7(0) ) +
592  c4(1) * ( -3 * c1(0) + 3 * c3(0) + 6 * c5(0) - 6 * c6(0) - 6 * c7(0) + 6 * c8(0) ) +
593  c5(1) * ( -12 * c1(0) - 14 * c3(0) - 6 * c4(0) + 60 * c6(0) + 8 * c7(0) + 4 * c8(0) ) +
594  c6(1) * ( 14 * c1(0) + 12 * c3(0) + 6 * c4(0) - 60 * c5(0) - 4 * c7(0) - 8 * c8(0) ) +
595  c7(1) * ( -14 * c3(0) + 6 * c4(0) - 8 * c5(0) + 4 * c6(0) + 12 * c8(0) ) +
596  c8(1) * ( 14 * c1(0) - 6 * c4(0) - 4 * c5(0) + 8 * c6(0) - 12 * c7(0) ) ) +
597  c3(2) * ( c1(1) * ( -3 * c2(0) + 3 * c4(0) + 6 * c5(0) + 6 * c6(0) - 6 * c7(0) - 6 * c8(0) ) +
598  c2(1) * ( 3 * c1(0) + 3 * c4(0) - 14 * c5(0) + 12 * c6(0) - 14 * c7(0) ) +
599  c4(1) * ( -3 * c1(0) - 3 * c2(0) + 14 * c6(0) - 12 * c7(0) + 14 * c8(0) ) +
600  c5(1) * ( -6 * c1(0) + 14 * c2(0) - 4 * c6(0) + 8 * c7(0) - 12 * c8(0) ) +
601  c6(1) * ( -6 * c1(0) - 12 * c2(0) - 14 * c4(0) + 4 * c5(0) + 60 * c7(0) + 8 * c8(0) ) +
602  c7(1) * ( 6 * c1(0) + 14 * c2(0) + 12 * c4(0) - 8 * c5(0) - 60 * c6(0) - 4 * c8(0) ) +
603  c8(1) * ( 6 * c1(0) - 14 * c4(0) + 12 * c5(0) - 8 * c6(0) + 4 * c7(0) ) ) +
604  c4(2) * ( c1(1) * ( -3 * c2(0) - 3 * c3(0) + 14 * c5(0) + 14 * c7(0) - 12 * c8(0) ) +
605  c2(1) * ( 3 * c1(0) - 3 * c3(0) - 6 * c5(0) + 6 * c6(0) + 6 * c7(0) - 6 * c8(0) ) +
606  c3(1) * ( 3 * c1(0) + 3 * c2(0) - 14 * c6(0) + 12 * c7(0) - 14 * c8(0) ) +
607  c5(1) * ( -14 * c1(0) + 6 * c2(0) + 12 * c6(0) - 8 * c7(0) + 4 * c8(0) ) +
608  c6(1) * ( -6 * c2(0) + 14 * c3(0) - 12 * c5(0) - 4 * c7(0) + 8 * c8(0) ) +
609  c7(1) * ( -14 * c1(0) - 6 * c2(0) - 12 * c3(0) + 8 * c5(0) + 4 * c6(0) + 60 * c8(0) ) +
610  c8(1) * ( 12 * c1(0) + 6 * c2(0) + 14 * c3(0) - 4 * c5(0) - 8 * c6(0) - 60 * c7(0) ) ) +
611  c5(2) * ( c1(1) * ( -12 * c2(0) - 6 * c3(0) - 14 * c4(0) + 4 * c6(0) + 8 * c7(0) + 60 * c8(0) ) +
612  c2(1) * ( 12 * c1(0) + 14 * c3(0) + 6 * c4(0) - 60 * c6(0) - 8 * c7(0) - 4 * c8(0) ) +
613  c3(1) * ( 6 * c1(0) - 14 * c2(0) + 4 * c6(0) - 8 * c7(0) + 12 * c8(0) ) +
614  c4(1) * ( 14 * c1(0) - 6 * c2(0) - 12 * c6(0) + 8 * c7(0) - 4 * c8(0) ) +
615  c6(1) * ( -4 * c1(0) + 60 * c2(0) - 4 * c3(0) + 12 * c4(0) - 32 * c7(0) - 32 * c8(0) ) +
616  c7(1) * ( -8 * c1(0) + 8 * c2(0) + 8 * c3(0) - 8 * c4(0) + 32 * c6(0) - 32 * c8(0) ) +
617  c8(1) * ( -60 * c1(0) + 4 * c2(0) - 12 * c3(0) + 4 * c4(0) + 32 * c6(0) + 32 * c7(0) ) ) +
618  c6(2) * ( c1(1) * ( 14 * c2(0) - 6 * c3(0) - 4 * c5(0) - 12 * c7(0) + 8 * c8(0) ) +
619  c2(1) * ( -14 * c1(0) - 12 * c3(0) - 6 * c4(0) + 60 * c5(0) + 4 * c7(0) + 8 * c8(0) ) +
620  c3(1) * ( 6 * c1(0) + 12 * c2(0) + 14 * c4(0) - 4 * c5(0) - 60 * c7(0) - 8 * c8(0) ) +
621  c4(1) * ( 6 * c2(0) - 14 * c3(0) + 12 * c5(0) + 4 * c7(0) - 8 * c8(0) ) +
622  c5(1) * ( 4 * c1(0) - 60 * c2(0) + 4 * c3(0) - 12 * c4(0) + 32 * c7(0) + 32 * c8(0) ) +
623  c7(1) * ( 12 * c1(0) - 4 * c2(0) + 60 * c3(0) - 4 * c4(0) - 32 * c5(0) - 32 * c8(0) ) +
624  c8(1) * ( -8 * c1(0) - 8 * c2(0) + 8 * c3(0) + 8 * c4(0) - 32 * c5(0) + 32 * c7(0) ) ) +
625  c7(2) * ( c1(1) * ( 6 * c3(0) - 14 * c4(0) - 8 * c5(0) + 12 * c6(0) + 4 * c8(0) ) +
626  c2(1) * ( 14 * c3(0) - 6 * c4(0) + 8 * c5(0) - 4 * c6(0) - 12 * c8(0) ) +
627  c3(1) * ( -6 * c1(0) - 14 * c2(0) - 12 * c4(0) + 8 * c5(0) + 60 * c6(0) + 4 * c8(0) ) +
628  c4(1) * ( 14 * c1(0) + 6 * c2(0) + 12 * c3(0) - 8 * c5(0) - 4 * c6(0) - 60 * c8(0) ) +
629  c5(1) * ( 8 * c1(0) - 8 * c2(0) - 8 * c3(0) + 8 * c4(0) - 32 * c6(0) + 32 * c8(0) ) +
630  c6(1) * ( -12 * c1(0) + 4 * c2(0) - 60 * c3(0) + 4 * c4(0) + 32 * c5(0) + 32 * c8(0) ) +
631  c8(1) * ( -4 * c1(0) + 12 * c2(0) - 4 * c3(0) + 60 * c4(0) - 32 * c5(0) - 32 * c6(0) ) ) +
632  c8(2) * ( c1(1) * ( 14 * c2(0) + 6 * c3(0) + 12 * c4(0) - 60 * c5(0) - 8 * c6(0) - 4 * c7(0) ) +
633  c2(1) * ( -14 * c1(0) + 6 * c4(0) + 4 * c5(0) - 8 * c6(0) + 12 * c7(0) ) +
634  c3(1) * ( -6 * c1(0) + 14 * c4(0) - 12 * c5(0) + 8 * c6(0) - 4 * c7(0) ) +
635  c4(1) * ( -12 * c1(0) - 6 * c2(0) - 14 * c3(0) + 4 * c5(0) + 8 * c6(0) + 60 * c7(0) ) +
636  c5(1) * ( 60 * c1(0) - 4 * c2(0) + 12 * c3(0) - 4 * c4(0) - 32 * c6(0) - 32 * c7(0) ) +
637  c6(1) * ( 8 * c1(0) + 8 * c2(0) - 8 * c3(0) - 8 * c4(0) + 32 * c5(0) - 32 * c7(0) ) +
638  c7(1) * ( 4 * c1(0) - 12 * c2(0) + 4 * c3(0) - 60 * c4(0) + 32 * c5(0) + 32 * c6(0) ) )
639  ) / 60.0;
640 }
641
642
645 {
646  IntegrationRule *iRule = new GaussIntegrationRule(1, NULL);
647  int points = iRule->getRequiredNumberOfIntegrationPoints(_Cube, order + 9);
648  iRule->SetUpPointsOnCube(points, _Unknown);
649  return iRule;
650 }
651
654 {
655  IntegrationRule *iRule = new GaussIntegrationRule(1, NULL);
656  int points = iRule->getRequiredNumberOfIntegrationPoints(_Square, order + 4);
657  iRule->SetUpPointsOnSquare(points, _Unknown);
658  return iRule;
659 }
660 } // end namespace oofem
double giveDeterminant() const
Returns the trace (sum of diagonal components) of the receiver.
Definition: floatmatrix.C:1408
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 POINT_TOL
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
bool solveForRhs(const FloatArray &b, FloatArray &answer, bool transpose=false)
Solves the system of linear equations .
Definition: floatmatrix.C:1112
double & at(int i)
Coefficient access function.
Definition: floatarray.h:131
virtual const FloatArray * giveVertexCoordinates(int i) const =0
virtual double evalNXIntegral(int iEdge, const FEICellGeometry &cellgeo)
Computes the integral .
Class representing a general abstraction for cell geometry.
Definition: feinterpol.h:62
void clear()
Definition: floatarray.h:206
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...
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
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.
Abstract base class representing integration rule.
void beDifferenceOf(const FloatArray &a, const FloatArray &b)
Sets receiver to be a - b.
Definition: floatarray.C:341
virtual void giveLocalDerivative(FloatMatrix &dN, const FloatArray &lcoords)
double distance(const FloatArray &x) const
Computes the distance between position represented by receiver and position given as parameter...
Definition: floatarray.C:489
virtual void edgeLocal2global(FloatArray &answer, int iedge, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates edge global coordinates from given local ones.
virtual IntegrationRule * giveBoundaryIntegrationRule(int order, int boundary)
Sets up a suitable integration rule for integrating over the requested boundary.
#define OOFEM_ERROR(...)
Definition: error.h:61
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
void resize(int n)
Checks size of receiver towards requested bounds.
Definition: intarray.C:124
virtual double giveCharacteristicLength(const FEICellGeometry &cellgeo) const
virtual double edgeGiveTransformationJacobian(int iedge, 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
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 giveJacobianMatrixAt(FloatMatrix &jacobianMatrix, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Gives the jacobian matrix at the local coordinates.
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.
virtual void evalN(FloatArray &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates the array of interpolation functions (shape functions) at given point.
double computeNorm() const
Computes the norm (or length) of the vector.
Definition: floatarray.C:840
virtual double surfaceEvalNormal(FloatArray &answer, int isurf, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates the normal out of the surface at given point.
void resize(int rows, int cols)
Checks size of receiver towards requested bounds.
Definition: floatmatrix.C:1358
virtual void computeLocalSurfaceMapping(IntArray &nodes, int iSurf)
virtual int getRequiredNumberOfIntegrationPoints(integrationDomain dType, int approxOrder)
Abstract service.
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 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 surfaceLocal2global(FloatArray &answer, int isurf, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates edge global coordinates from given local ones.
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
int giveSize() const
Definition: floatarray.h:218
virtual void computeLocalEdgeMapping(IntArray &edgeNodes, int iedge)
the oofem namespace is to define a context or scope in which all oofem names are defined.
virtual int global2local(FloatArray &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates local coordinates from given global ones.
void beInverseOf(const FloatMatrix &src)
Modifies receiver to become inverse of given parameter.
Definition: floatmatrix.C:835
double normalize()
Definition: floatarray.C:828
int giveNumberOfRows() const
Returns number of rows of receiver.
Definition: floatmatrix.h:156
virtual double surfaceGiveTransformationJacobian(int isurf, 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.