OOFEM  2.4
OOFEM.org - Object Oriented Finite Element Solver
gaussintegrationrule.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
22  * License as published by the Free Software Foundation; either
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 "gaussintegrationrule.h"
36 #include "gausspoint.h"
37 #include "element.h"
38 #include "floatarray.h"
39 #include "mathfem.h"
40 
41 namespace oofem {
42 // initialize class member
43 
44 GaussIntegrationRule :: GaussIntegrationRule(int n, Element *e,
45  int startIndx, int endIndx, bool dynamic) :
46  IntegrationRule(n, e, startIndx, endIndx, dynamic) { }
47 
48 GaussIntegrationRule :: GaussIntegrationRule(int n, Element *e) :
49  IntegrationRule(n, e) { }
50 
51 GaussIntegrationRule :: ~GaussIntegrationRule()
52 { }
53 
54 
55 int
56 GaussIntegrationRule :: SetUpPointsOnLine(int nPoints, MaterialMode mode)
57 {
58  FloatArray coords_xi, weights;
59  this->giveLineCoordsAndWeights(nPoints, coords_xi, weights);
60  this->gaussPoints.resize( nPoints );
61 
62  for ( int i = 1; i <= nPoints; i++ ) {
63  this->gaussPoints [ i - 1 ] = new GaussPoint(this, i, {coords_xi.at(i)}, weights.at ( i ), mode);
64  }
65 
66  this->intdomain = _Line;
67  return this->giveNumberOfIntegrationPoints();
68 }
69 
70 
71 int
72 GaussIntegrationRule :: SetUpPointsOn2DEmbeddedLine(int nPoints, MaterialMode mode, const FloatArray &coord0, const FloatArray &coord1)
73 {
74  FloatArray coords_xi, weights;
75  this->giveLineCoordsAndWeights(nPoints, coords_xi, weights);
76  this->gaussPoints.resize( nPoints );
77 
78  for ( int i = 1; i <= nPoints; i++ ) {
79  double x = ( coords_xi.at(i) + 1.0 ) * 0.5;
80  FloatArray subpatchCoord = {x};
81 
82  this->gaussPoints [ i - 1 ] = new GaussPoint(this, i, weights.at ( i ), mode);
83 
84  this->gaussPoints [ i - 1 ]->setSubPatchCoordinates(subpatchCoord);
85 
86  FloatArray globalCoord = {
87  ( 1. - x ) * coord0.at(1) + x * coord1.at(1),
88  ( 1. - x ) * coord0.at(2) + x * coord1.at(2) };
89  this->gaussPoints [ i - 1 ]->setGlobalCoordinates(globalCoord);
90 
91  FloatArray naturalCoord;
92  this->giveElement()->computeLocalCoordinates(naturalCoord, globalCoord);
93  this->gaussPoints [ i - 1 ]->setNaturalCoordinates(naturalCoord);
94  }
95 
96  this->intdomain = _Embedded2dLine;
97  return this->giveNumberOfIntegrationPoints();
98 }
99 
100 
101 int
102 GaussIntegrationRule :: SetUpPointsOnSquare(int nPoints, MaterialMode mode)
103 //GaussIntegrationRule :: SetUpPointsOnSquare(int nPoints_xi1, int nPoints_xi2, MaterialMode mode)
104 {
105  int nPoints_xi1 = ( int ) floor( sqrt( double ( nPoints ) ) );
106  int nPoints_xi2 = nPoints_xi1;
107  FloatArray coords_xi1, weights1, coords_xi2, weights2;
108  this->giveLineCoordsAndWeights(nPoints_xi1, coords_xi1, weights1);
109  this->giveLineCoordsAndWeights(nPoints_xi2, coords_xi2, weights2);
110  this->gaussPoints.resize( nPoints_xi1 * nPoints_xi2 );
111  int count = 0;
112  for ( int i = 1; i <= nPoints_xi1; i++ ) {
113  for ( int j = 1; j <= nPoints_xi2; j++ ) {
114  this->gaussPoints [ count ] = new GaussPoint(this, count + 1, {coords_xi1.at(i), coords_xi2.at(j)},
115  weights1.at ( i ) *weights2.at ( j ), mode);
116  count++;
117  }
118  }
119 
120  this->intdomain = _Square;
121  return this->giveNumberOfIntegrationPoints();
122 }
123 
124 int
125 GaussIntegrationRule :: SetUpPointsOn3dDegShell(int nPointsXY, int nPointsZ, MaterialMode mode)
126 //GaussIntegrationRule :: SetUpPointsOnCube(int nPoints_xi1, int nPoints_xi2, int nPoints_xi3, MaterialMode mode)
127 {
128  int nPoints_xi1 = ( int ) floor(sqrt( double ( nPointsXY ) ) + 0.1 );
129  int nPoints_xi2 = nPoints_xi1;
130  int nPoints_xi3 = nPointsZ;
131  FloatArray coords_xi1, weights1, coords_xi2, weights2, coords_xi3, weights3;
132  this->giveLineCoordsAndWeights(nPoints_xi1, coords_xi1, weights1);
133  this->giveLineCoordsAndWeights(nPoints_xi2, coords_xi2, weights2);
134  this->giveLineCoordsAndWeights(nPoints_xi3, coords_xi3, weights3);
135  this->gaussPoints.resize( nPoints_xi1 * nPoints_xi2 * nPoints_xi3 );
136  int count = 0;
137  // the order of cycles cannot be changed due to MITC4Shell::giveMomentsAt(...)
138  for ( int i = 1; i <= nPoints_xi1; i++ ) {
139  for ( int j = 1; j <= nPoints_xi2; j++ ) {
140  for ( int k = 1; k <= nPoints_xi3; k++ ) {
141  this->gaussPoints [ count ] = new GaussPoint(this, count + 1, {coords_xi1.at(i), coords_xi2.at(j), coords_xi3.at(k)},
142  weights1.at ( i ) *weights2.at ( j ) *weights3.at ( k ), mode);
143  count++;
144  }
145  }
146  }
147 
148  this->intdomain = _3dDegShell;
149  return this->giveNumberOfIntegrationPoints();
150 }
151 
152 int
153 GaussIntegrationRule :: SetUpPointsOn3dDegShellLayers(int nPointsXY, int nPointsZ, MaterialMode mode, const FloatArray &layerThickness)
154 {
155  int nPoints_xi1 = ( int ) floor(sqrt( double ( nPointsXY ) ) + 0.1 );
156  int nPoints_xi2 = nPoints_xi1;
157  int nPoints_xi3 = nPointsZ;
158  FloatArray coords_xi1, weights1, coords_xi2, weights2, coords_xi3, weights3;
159  this->giveLineCoordsAndWeights(nPoints_xi1, coords_xi1, weights1);
160  this->giveLineCoordsAndWeights(nPoints_xi2, coords_xi2, weights2);
161  this->giveLineCoordsAndWeights(nPoints_xi3, coords_xi3, weights3);
162  int pointsPerLayer = nPoints_xi1 * nPoints_xi2 * nPoints_xi3;
163  this->gaussPoints.resize( pointsPerLayer * layerThickness.giveSize() );
164  int count = 0;
165  double totalThickness = layerThickness.sum();
166  // bottom is scaled so that it goes from -1 to 1.
167  double bottom = -1.;
168  double scaledThickness;
169 
170  for ( int t = 1; t <= layerThickness.giveSize(); t++ ) {
171  scaledThickness = layerThickness.at(t) / totalThickness;
172  // the order of cycles cannot be changed due to MITC4Shell::giveMomentsAt(...)
173  for ( int i = 1; i <= nPoints_xi1; i++ ) {
174  for ( int j = 1; j <= nPoints_xi2; j++ ) {
175  for ( int k = 1; k <= nPoints_xi3; k++ ) {
176  this->gaussPoints [ count ] = new GaussPoint(this, count + 1, {coords_xi1.at(i), coords_xi2.at(j), ( coords_xi3.at(k) + 1. ) * scaledThickness + bottom},
177  weights1.at(i)*weights2.at(j)*weights3.at(k)*scaledThickness, mode);
178  count++;
179  }
180  }
181  }
182  bottom += 2.0 * scaledThickness;
183  }
184  this->intdomain = _3dDegShell;
185  return this->giveNumberOfIntegrationPoints();
186 }
187 
188 
189 int
190 GaussIntegrationRule :: SetUpPointsOnCube(int nPoints, MaterialMode mode)
191 //GaussIntegrationRule :: SetUpPointsOnCube(int nPoints_xi1, int nPoints_xi2, int nPoints_xi3, MaterialMode mode)
192 {
193  int nPoints_xi1 = ( int ) floor(cbrt( double ( nPoints ) ) + 0.5);
194  int nPoints_xi2 = nPoints_xi1;
195  int nPoints_xi3 = nPoints_xi1;
196  FloatArray coords_xi1, weights1, coords_xi2, weights2, coords_xi3, weights3;
197  this->giveLineCoordsAndWeights(nPoints_xi1, coords_xi1, weights1);
198  this->giveLineCoordsAndWeights(nPoints_xi2, coords_xi2, weights2);
199  this->giveLineCoordsAndWeights(nPoints_xi3, coords_xi3, weights3);
200  this->gaussPoints.resize( nPoints_xi1 * nPoints_xi2 * nPoints_xi3 );
201  int count = 0;
202  for ( int i = 1; i <= nPoints_xi1; i++ ) {
203  for ( int j = 1; j <= nPoints_xi2; j++ ) {
204  for ( int k = 1; k <= nPoints_xi3; k++ ) {
205  this->gaussPoints [ count ] = new GaussPoint(this, count + 1, {coords_xi1.at(i), coords_xi2.at(j), coords_xi3.at(k)},
206  weights1.at ( i ) *weights2.at ( j ) *weights3.at ( k ), mode);
207  count++;
208  }
209  }
210  }
211 
212  this->intdomain = _Cube;
213  return this->giveNumberOfIntegrationPoints();
214 }
215 
216 
217 int GaussIntegrationRule :: SetUpPointsOnCubeLayers(int nPoints1, int nPoints2, int nPointsDepth, MaterialMode mode, const FloatArray &layerThickness)
218 {
219  FloatArray coords_xi1, weights1, coords_xi2, weights2, coords_xi3, weights3;
220  this->giveLineCoordsAndWeights(nPoints1, coords_xi1, weights1);
221  this->giveLineCoordsAndWeights(nPoints2, coords_xi2, weights2);
222  this->giveLineCoordsAndWeights(nPointsDepth, coords_xi3, weights3);
223  int pointsPerLayer = nPoints1 * nPoints2 * nPointsDepth;
224  this->gaussPoints.resize( pointsPerLayer * layerThickness.giveSize() );
225 
226  int count = 0;
227  double totalThickness = layerThickness.sum();
228 
229  // bottom is scaled so that it goes from -1 to 1.
230  double bottom = -1.;
231  double scaledThickness;
232  for ( int t = 1; t <= layerThickness.giveSize(); t++ ) {
233  scaledThickness = layerThickness.at(t) / totalThickness;
234  for ( int i = 1; i <= nPoints1; i++ ) {
235  for ( int j = 1; j <= nPoints2; j++ ) {
236  for ( int k = 1; k <= nPointsDepth; k++ ) {
237  this->gaussPoints [ count ] = new GaussPoint(this, count + 1,
238  {coords_xi1.at(i), coords_xi2.at(j), ( coords_xi3.at(k) + 1. ) * scaledThickness + bottom},
239  weights1.at(i)*weights2.at(j)*weights3.at(k)*scaledThickness, mode);
240  count++;
241  }
242  }
243  }
244  bottom += 2.0 * scaledThickness;
245  }
246 
247  this->intdomain = _Cube;
248  return this->giveNumberOfIntegrationPoints();
249 }
250 
251 
252 int
253 GaussIntegrationRule :: SetUpPointsOnTriangle(int nPoints, MaterialMode mode)
254 {
255  FloatArray coords_xi1, coords_xi2, weights;
256  this->giveTriCoordsAndWeights(nPoints, coords_xi1, coords_xi2, weights);
257  this->gaussPoints.resize( nPoints );
258 
259  for ( int i = 1; i <= nPoints; i++ ) {
260  this->gaussPoints [ i - 1 ] = new GaussPoint(this, i, {coords_xi1.at(i), coords_xi2.at(i)}, weights.at ( i ), mode);
261  }
262 
263  this->intdomain = _Triangle;
264  return this->giveNumberOfIntegrationPoints();
265 }
266 
267 
268 int
269 GaussIntegrationRule :: SetUpPointsOnTetrahedra(int nPoints, MaterialMode mode)
270 {
271  FloatArray coords_xi1, coords_xi2, coords_xi3, weights;
272  this->giveTetCoordsAndWeights(nPoints, coords_xi1, coords_xi2, coords_xi3, weights);
273  this->gaussPoints.resize( nPoints );
274 
275  for ( int i = 1; i <= nPoints; i++ ) {
276  this->gaussPoints [ i - 1 ] = new GaussPoint(this, i,
277  {coords_xi1.at(i), coords_xi2.at(i), coords_xi3.at(i)}, weights.at ( i ), mode);
278  }
279 
280  this->intdomain = _Tetrahedra;
281  return this->giveNumberOfIntegrationPoints();
282 }
283 
284 
285 int
286 GaussIntegrationRule :: SetUpPointsOnWedge(int nPointsTri, int nPointsDepth, MaterialMode mode)
287 {
288  FloatArray coords_xi1, coords_xi2, coords_xi3, weightsTri, weightsDepth;
289  this->giveTriCoordsAndWeights(nPointsTri, coords_xi1, coords_xi2, weightsTri);
290  this->giveLineCoordsAndWeights(nPointsDepth, coords_xi3, weightsDepth);
291  this->gaussPoints.resize( nPointsTri * nPointsDepth );
292 
293  int count = 0;
294  for ( int i = 1; i <= nPointsTri; i++ ) {
295  for ( int j = 1; j <= nPointsDepth; j++ ) {
296  this->gaussPoints [ count ] = new GaussPoint(this, count + 1, {coords_xi1.at(i), coords_xi2.at(i), coords_xi3.at(j)},
297  weightsTri.at ( i ) *weightsDepth.at ( j ), mode);
298  count++;
299  }
300  }
301 
302  this->intdomain = _Wedge;
303  return this->giveNumberOfIntegrationPoints();
304 }
305 
306 int
307 GaussIntegrationRule :: SetUpPointsOnWedgeLayers(int nPointsTri, int nPointsDepth, MaterialMode mode, const FloatArray &layerThickness)
308 {
309  FloatArray coords_xi1, coords_xi2, coords_xi3, weightsTri, weightsDepth;
310  this->giveTriCoordsAndWeights(nPointsTri, coords_xi1, coords_xi2, weightsTri);
311  this->giveLineCoordsAndWeights(nPointsDepth, coords_xi3, weightsDepth);
312  int pointsPerLayer = nPointsTri * nPointsDepth;
313  this->gaussPoints.resize( pointsPerLayer * layerThickness.giveSize() );
314 
315  int count = 0;
316  double totalThickness = layerThickness.sum();
317 
318  // bottom is scaled so that it goes from -1 to 1.
319  double bottom = -1.;
320  double scaledThickness;
321  for ( int k = 1; k <= layerThickness.giveSize(); k++ ) {
322  scaledThickness = layerThickness.at(k) / totalThickness;
323  for ( int i = 1; i <= nPointsTri; i++ ) {
324  for ( int j = 1; j <= nPointsDepth; j++ ) {
325  this->gaussPoints [ count ] = new GaussPoint(this, count + 1,
326  {coords_xi1.at(i), coords_xi2.at(i), ( coords_xi3.at(j) + 1. ) * scaledThickness + bottom},
327  weightsTri.at ( i ) * ( weightsDepth.at ( j ) *scaledThickness ), mode);
328  count++;
329  }
330  }
331  bottom += 2.0 * scaledThickness;
332  }
333 
334  this->intdomain = _Wedge;
335  return this->giveNumberOfIntegrationPoints();
336 }
337 
338 int
339 GaussIntegrationRule :: getRequiredNumberOfIntegrationPoints(integrationDomain dType,
340  int approxOrder)
341 {
342  // Returns the required number of gp's to exactly integrate a polynomial of order 'approxOrder'
343  int requiredNIP;
344  if ( approxOrder < 0 ) {
345  return 0;
346  }
347 
348  switch ( dType ) {
349  case _Line:
350  requiredNIP = (int)ceil(( approxOrder + 1 ) / 2);
351 
352  if ( requiredNIP <= 8 ) {
353  return requiredNIP;
354  }
355 
356  if ( requiredNIP <= 16 ) {
357  return 16;
358  }
359 
360  if ( requiredNIP <= 24 ) {
361  return 24;
362  }
363 
364  if ( requiredNIP <= 32 ) {
365  return 32;
366  }
367 
368  if ( requiredNIP <= 64 ) {
369  return 64;
370  }
371 
372  return -1;
373 
374  case _Triangle:
375  if ( approxOrder <= 1 ) {
376  return 1;
377  }
378 
379  if ( approxOrder <= 2 ) {
380  return 3;
381  }
382 
383  if ( approxOrder <= 3 ) {
384  return 4;
385  }
386 
387  if ( approxOrder <= 4 ) {
388  return 6;
389  }
390 
391  if ( approxOrder <= 5 ) {
392  return 7;
393  }
394 
395  if ( approxOrder <= 6 ) {
396  return 12;
397  }
398 
399  if ( approxOrder <= 7 ) {
400  return 13;
401  }
402 
403  if ( approxOrder <= 8 ) {
404  return 16;
405  }
406 
407  if ( approxOrder <= 9 ) {
408  return 19;
409  }
410 
411  if ( approxOrder <= 10 ) {
412  return 25;
413  }
414 
415  return -1;
416 
417  case _Square:
418  requiredNIP = max( (int)ceil( (double)( approxOrder + 1 ) / (double)2), 2 );
419  requiredNIP *= requiredNIP;
420  if ( requiredNIP > 64 * 64 ) {
421  return -1;
422  }
423 
424  return requiredNIP;
425 
426  case _Cube:
427  requiredNIP = max( (int)ceil( (double)( approxOrder + 1 ) / (double)2), 2 );
428  requiredNIP *= requiredNIP * requiredNIP;
429  if ( requiredNIP > 64 * 64 * 64 ) { // = 262 144 gp's maybe overkill :)
430  return -1;
431  }
432 
433  return requiredNIP;
434 
435  case _Tetrahedra:
436  if ( approxOrder <= 1 ) {
437  return 1;
438  }
439 
440  if ( approxOrder <= 2 ) {
441  return 4;
442  }
443 
444  if ( approxOrder <= 3 ) {
445  return 5;
446  }
447 
448  if ( approxOrder <= 4 ) {
449  return 11; // unsure
450  }
451 
452  if ( approxOrder <= 5 ) {
453  return 15;
454  }
455 
456  if ( approxOrder <= 6 ) {
457  return 24;
458  }
459 
460  if ( approxOrder <= 8 ) {
461  return 45;
462  }
463 
464  break;
465 
467  case _Wedge:
468  if ( approxOrder <= 1 ) { // (lin tri) x (lin line) = 1 x 1 = 1
469  return 1;
470  }
471 
472  if ( approxOrder <= 2 ) { // (quad tri) x (quad line) = 3 x 2 = 6
473  return 6;
474  }
475 
476  if ( approxOrder <= 3 ) { // (cubic tri) x (cubic line) = 4 x 2 = 8
477  return 8;
478  }
479 
480  if ( approxOrder <= 4 ) { // (quartic tri) x (quartic line) = 6 x 3 = 18
481  return 18;
482  }
483  default:
484  OOFEM_ERROR("unknown integrationDomain");
485  }
486 
487  return -1;
488 }
489 
490 
491 void
492 GaussIntegrationRule :: giveTetCoordsAndWeights(int nPoints, FloatArray &coords_xi1, FloatArray &coords_xi2, FloatArray &coords_xi3, FloatArray &weights)
493 {
494  double a, b, c, w;
495  coords_xi1.resize(nPoints);
496  coords_xi2.resize(nPoints);
497  coords_xi3.resize(nPoints);
498  weights.resize(nPoints);
499  switch ( nPoints ) {
500  case 1:
501  a = 1. / 4.;
502  w = 1. / 6.;
503  coords_xi1(0) = a;
504  coords_xi2(0) = a;
505  coords_xi3(0) = a;
506  weights(0) = w;
507  break;
508 
509  case 4: // quadratic formulae
510  a = ( 5. + 3. * sqrt(5.) ) / 20.;
511  b = ( 5. - sqrt(5.) ) / 20.;
512  w = 1. / 24.;
513  coords_xi1(0) = a;
514  coords_xi2(0) = b;
515  coords_xi3(0) = b;
516  weights(0) = w;
517  coords_xi1(1) = b;
518  coords_xi2(1) = a;
519  coords_xi3(1) = b;
520  weights(1) = w;
521  coords_xi1(2) = b;
522  coords_xi2(2) = b;
523  coords_xi3(2) = a;
524  weights(2) = w;
525  coords_xi1(3) = b;
526  coords_xi2(3) = b;
527  coords_xi3(3) = b;
528  weights(3) = w;
529  break;
530 
531  case 5: // cubic formulae
532  a = 1. / 4.;
533  w = -2. / 15.;
534  coords_xi1(0) = a;
535  coords_xi2(0) = a;
536  coords_xi3(0) = a;
537  weights(0) = w;
538  a = 1. / 2.;
539  b = 1. / 6.;
540  w = 3. / 40.;
541  coords_xi1(1) = a;
542  coords_xi2(1) = b;
543  coords_xi3(1) = b;
544  weights(1) = w;
545  coords_xi1(2) = b;
546  coords_xi2(2) = a;
547  coords_xi3(2) = b;
548  weights(2) = w;
549  coords_xi1(3) = b;
550  coords_xi2(3) = b;
551  coords_xi3(3) = a;
552  weights(3) = w;
553  coords_xi1(4) = b;
554  coords_xi2(4) = b;
555  coords_xi3(4) = b;
556  weights(4) = w;
557  break;
558 
559  case 11: // exact x^4
560  a = 1. / 4.;
561  w = -74. / 5625.;
562  coords_xi1(0) = a;
563  coords_xi2(0) = a;
564  coords_xi3(0) = a;
565  weights(0) = w;
566  a = 5. / 70.;
567  b = 11. / 14.;
568  w = 343. / 45000.;
569  coords_xi1(1) = b;
570  coords_xi2(1) = a;
571  coords_xi3(1) = a;
572  weights(1) = w;
573  coords_xi1(2) = a;
574  coords_xi2(2) = b;
575  coords_xi3(2) = a;
576  weights(2) = w;
577  coords_xi1(3) = a;
578  coords_xi2(3) = a;
579  coords_xi3(3) = b;
580  weights(3) = w;
581  coords_xi1(4) = a;
582  coords_xi2(4) = a;
583  coords_xi3(4) = a;
584  weights(4) = w;
585  a = ( 1. + sqrt(5. / 14.) ) / 4.;
586  b = ( 1. - sqrt(5. / 14.) ) / 4.;
587  w = 28. / 1125.;
588  coords_xi1(5) = a;
589  coords_xi2(5) = a;
590  coords_xi3(5) = b;
591  weights(5) = w;
592  coords_xi1(6) = a;
593  coords_xi2(6) = b;
594  coords_xi3(6) = a;
595  weights(6) = w;
596  coords_xi1(7) = a;
597  coords_xi2(7) = b;
598  coords_xi3(7) = b;
599  weights(7) = w;
600  coords_xi1(8) = b;
601  coords_xi2(8) = a;
602  coords_xi3(8) = a;
603  weights(8) = w;
604  coords_xi1(9) = b;
605  coords_xi2(9) = a;
606  coords_xi3(9) = b;
607  weights(9) = w;
608  coords_xi1(10) = b;
609  coords_xi2(10) = b;
610  coords_xi3(10) = a;
611  weights(10) = w;
612  break;
613 
614  case 15: // exact for x^5
615  // Derived by Patrick Keast, MODERATE-DEGREE TETRAHEDRAL QUADRATURE FORMULAS
616  a = 1. / 4.;
617  w = 0.302836780970891856e-1;
618  coords_xi1(0) = a;
619  coords_xi2(0) = a;
620  coords_xi3(0) = a;
621  weights(0) = w;
622  a = 0.0;
623  b = 1. / 3.;
624  w = 0.602678571428571597e-2;
625  coords_xi1(1) = a;
626  coords_xi2(1) = b;
627  coords_xi3(1) = b;
628  weights(1) = w;
629  coords_xi1(2) = b;
630  coords_xi2(2) = a;
631  coords_xi3(2) = b;
632  weights(2) = w;
633  coords_xi1(3) = b;
634  coords_xi2(3) = b;
635  coords_xi3(3) = a;
636  weights(3) = w;
637  coords_xi1(4) = b;
638  coords_xi2(4) = b;
639  coords_xi3(4) = b;
640  weights(4) = w;
641  a = 8. / 11.;
642  b = 1. / 11.;
643  w = 0.116452490860289742e-1;
644  coords_xi1(5) = a;
645  coords_xi2(5) = b;
646  coords_xi3(5) = b;
647  weights(5) = w;
648  coords_xi1(6) = b;
649  coords_xi2(6) = a;
650  coords_xi3(6) = b;
651  weights(6) = w;
652  coords_xi1(7) = b;
653  coords_xi2(7) = b;
654  coords_xi3(7) = a;
655  weights(7) = w;
656  coords_xi1(8) = b;
657  coords_xi2(8) = b;
658  coords_xi3(8) = b;
659  weights(8) = w;
660  a = 0.665501535736642813e-1;
661  b = 0.433449846426335728;
662  w = 0.109491415613864534e-1;
663  coords_xi1(9) = a;
664  coords_xi2(9) = a;
665  coords_xi3(9) = b;
666  weights(9) = w;
667  coords_xi1(10) = a;
668  coords_xi2(10) = b;
669  coords_xi3(10) = a;
670  weights(10) = w;
671  coords_xi1(11) = a;
672  coords_xi2(11) = b;
673  coords_xi3(11) = b;
674  weights(11) = w;
675  coords_xi1(12) = b;
676  coords_xi2(12) = a;
677  coords_xi3(12) = a;
678  weights(12) = w;
679  coords_xi1(13) = b;
680  coords_xi2(13) = a;
681  coords_xi3(13) = b;
682  weights(13) = w;
683  coords_xi1(14) = b;
684  coords_xi2(14) = b;
685  coords_xi3(14) = a;
686  weights(14) = w;
687  break;
688 
689  case 24: // Exact for x^6
690  // Derived by Patrick Keast, MODERATE-DEGREE TETRAHEDRAL QUADRATURE FORMULAS
691  // See also: Monomial cubature rules since “Stroud”: a compilation
692  a = 0.356191386222544953;
693  b = 0.214602871259151684;
694  w = 0.665379170969464506e-2;
695  coords_xi1(0) = a;
696  coords_xi2(0) = b;
697  coords_xi3(0) = b;
698  weights(0) = w;
699  coords_xi1(1) = b;
700  coords_xi2(1) = a;
701  coords_xi3(1) = b;
702  weights(1) = w;
703  coords_xi1(2) = b;
704  coords_xi2(2) = b;
705  coords_xi3(2) = a;
706  weights(2) = w;
707  coords_xi1(3) = b;
708  coords_xi2(3) = b;
709  coords_xi3(3) = b;
710  weights(3) = w;
711  a = 0.877978124396165982;
712  b = 0.406739585346113397e-1;
713  w = 0.167953517588677620e-2;
714  coords_xi1(4) = a;
715  coords_xi2(4) = b;
716  coords_xi3(4) = b;
717  weights(4) = w;
718  coords_xi1(5) = b;
719  coords_xi2(5) = a;
720  coords_xi3(5) = b;
721  weights(5) = w;
722  coords_xi1(6) = b;
723  coords_xi2(6) = b;
724  coords_xi3(6) = a;
725  weights(6) = w;
726  coords_xi1(7) = b;
727  coords_xi2(7) = b;
728  coords_xi3(7) = b;
729  weights(7) = w;
730  a = 0.329863295731730594e-1;
731  b = 0.322337890142275646;
732  w = 0.922619692394239843e-2;
733  coords_xi1(8) = a;
734  coords_xi2(8) = b;
735  coords_xi3(8) = b;
736  weights(8) = w;
737  coords_xi1(9) = b;
738  coords_xi2(9) = a;
739  coords_xi3(9) = b;
740  weights(9) = w;
741  coords_xi1(10) = b;
742  coords_xi2(10) = b;
743  coords_xi3(10) = a;
744  weights(10) = w;
745  coords_xi1(11) = b;
746  coords_xi2(11) = b;
747  coords_xi3(11) = b;
748  weights(11) = w;
749  a = 0.603005664791649076;
750  b = 0.269672331458315867;
751  c = 0.636610018750175299e-1;
752  w = 0.803571428571428248e-2;
753  // 12 permutations follows here, a, b, c
754  coords_xi1(12) = a;
755  coords_xi2(12) = b;
756  coords_xi3(12) = c;
757  weights(12) = w;
758  coords_xi1(13) = b;
759  coords_xi2(13) = a;
760  coords_xi3(13) = c;
761  weights(13) = w;
762  coords_xi1(14) = a;
763  coords_xi2(14) = c;
764  coords_xi3(14) = b;
765  weights(14) = w;
766  coords_xi1(15) = b;
767  coords_xi2(15) = c;
768  coords_xi3(15) = a;
769  weights(15) = w;
770  coords_xi1(16) = c;
771  coords_xi2(16) = a;
772  coords_xi3(16) = b;
773  weights(16) = w;
774  coords_xi1(17) = c;
775  coords_xi2(17) = b;
776  coords_xi3(17) = a;
777  weights(17) = w;
778  coords_xi1(18) = a;
779  coords_xi2(18) = c;
780  coords_xi3(18) = c;
781  weights(18) = w;
782  coords_xi1(19) = b;
783  coords_xi2(19) = c;
784  coords_xi3(19) = c;
785  weights(19) = w;
786  coords_xi1(20) = c;
787  coords_xi2(20) = a;
788  coords_xi3(20) = c;
789  weights(20) = w;
790  coords_xi1(21) = c;
791  coords_xi2(21) = b;
792  coords_xi3(21) = c;
793  weights(21) = w;
794  coords_xi1(22) = c;
795  coords_xi2(22) = c;
796  coords_xi3(22) = a;
797  weights(22) = w;
798  coords_xi1(23) = c;
799  coords_xi2(23) = c;
800  coords_xi3(23) = b;
801  weights(23) = w;
802  break;
803 
804  case 45: // Exact for x^8
805  // Derived by Patrick Keast, MODERATE-DEGREE TETRAHEDRAL QUADRATURE FORMULAS
806  // See also: Monomial cubature rules since “Stroud”: a compilation
807  a = 1. / 4.;
808  w = -0.393270066412926145e-1;
809  coords_xi1(0) = a;
810  coords_xi2(0) = a;
811  coords_xi3(0) = a;
812  weights(0) = w;
813  a = 0.617587190300082967;
814  b = 0.127470936566639015;
815  w = 0.408131605934270525e-2;
816  coords_xi1(1) = a;
817  coords_xi2(1) = b;
818  coords_xi3(1) = b;
819  weights(1) = w;
820  coords_xi1(2) = b;
821  coords_xi2(2) = a;
822  coords_xi3(2) = b;
823  weights(2) = w;
824  coords_xi1(3) = b;
825  coords_xi2(3) = b;
826  coords_xi3(3) = a;
827  weights(3) = w;
828  coords_xi1(4) = b;
829  coords_xi2(4) = b;
830  coords_xi3(4) = b;
831  weights(4) = w;
832  a = 0.903763508822103123;
833  b = 0.320788303926322960e-1;
834  w = 0.658086773304341943e-3;
835  coords_xi1(5) = a;
836  coords_xi2(5) = b;
837  coords_xi3(5) = b;
838  weights(5) = w;
839  coords_xi1(6) = b;
840  coords_xi2(6) = a;
841  coords_xi3(6) = b;
842  weights(6) = w;
843  coords_xi1(7) = b;
844  coords_xi2(7) = b;
845  coords_xi3(7) = a;
846  weights(7) = w;
847  coords_xi1(8) = b;
848  coords_xi2(8) = b;
849  coords_xi3(8) = b;
850  weights(8) = w;
851  a = 0.450222904356718978;
852  b = 0.497770956432810185e-1;
853  w = 0.438425882512284693e-2;
854  coords_xi1(9) = a;
855  coords_xi2(9) = a;
856  coords_xi3(9) = b;
857  weights(9) = w;
858  coords_xi1(10) = a;
859  coords_xi2(10) = b;
860  coords_xi3(10) = a;
861  weights(10) = w;
862  coords_xi1(11) = a;
863  coords_xi2(11) = b;
864  coords_xi3(11) = b;
865  weights(11) = w;
866  coords_xi1(12) = b;
867  coords_xi2(12) = a;
868  coords_xi3(12) = a;
869  weights(12) = w;
870  coords_xi1(13) = b;
871  coords_xi2(13) = a;
872  coords_xi3(13) = b;
873  weights(13) = w;
874  coords_xi1(14) = b;
875  coords_xi2(14) = b;
876  coords_xi3(14) = a;
877  weights(14) = w;
878  a = 0.316269552601450060;
879  b = 0.183730447398549945;
880  w = 0.138300638425098166e-1;
881  coords_xi1(15) = a;
882  coords_xi2(15) = a;
883  coords_xi3(15) = b;
884  weights(15) = w;
885  coords_xi1(16) = a;
886  coords_xi2(16) = b;
887  coords_xi3(16) = a;
888  weights(16) = w;
889  coords_xi1(17) = a;
890  coords_xi2(17) = b;
891  coords_xi3(17) = b;
892  weights(17) = w;
893  coords_xi1(18) = b;
894  coords_xi2(18) = a;
895  coords_xi3(18) = a;
896  weights(18) = w;
897  coords_xi1(19) = b;
898  coords_xi2(19) = a;
899  coords_xi3(19) = b;
900  weights(19) = w;
901  coords_xi1(20) = b;
902  coords_xi2(20) = b;
903  coords_xi3(20) = a;
904  weights(20) = w;
905  a = 0.513280033360881072;
906  b = 0.229177878448171174e-1;
907  c = 0.231901089397150906;
908  w = 0.424043742468372453e-2;
909  coords_xi1(21) = a;
910  coords_xi2(21) = b;
911  coords_xi3(21) = c;
912  weights(21) = w;
913  coords_xi1(22) = b;
914  coords_xi2(22) = a;
915  coords_xi3(22) = c;
916  weights(22) = w;
917  coords_xi1(23) = a;
918  coords_xi2(23) = c;
919  coords_xi3(23) = b;
920  weights(23) = w;
921  coords_xi1(24) = b;
922  coords_xi2(24) = c;
923  coords_xi3(24) = a;
924  weights(24) = w;
925  coords_xi1(25) = c;
926  coords_xi2(25) = a;
927  coords_xi3(25) = b;
928  weights(25) = w;
929  coords_xi1(26) = c;
930  coords_xi2(26) = b;
931  coords_xi3(26) = a;
932  weights(26) = w;
933  coords_xi1(27) = a;
934  coords_xi2(27) = c;
935  coords_xi3(27) = c;
936  weights(27) = w;
937  coords_xi1(28) = b;
938  coords_xi2(28) = c;
939  coords_xi3(28) = c;
940  weights(28) = w;
941  coords_xi1(29) = c;
942  coords_xi2(29) = a;
943  coords_xi3(29) = c;
944  weights(29) = w;
945  coords_xi1(30) = c;
946  coords_xi2(30) = b;
947  coords_xi3(30) = c;
948  weights(30) = w;
949  coords_xi1(31) = c;
950  coords_xi2(31) = c;
951  coords_xi3(31) = a;
952  weights(31) = w;
953  coords_xi1(32) = c;
954  coords_xi2(32) = c;
955  coords_xi3(32) = b;
956  weights(32) = w;
957  a = 0.193746475248804382;
958  b = 0.730313427807538396;
959  c = 0.379700484718286102e-1;
960  w = 0.223873973961420164e-2;
961  coords_xi1(33) = a;
962  coords_xi2(33) = b;
963  coords_xi3(33) = c;
964  weights(33) = w;
965  coords_xi1(34) = b;
966  coords_xi2(34) = a;
967  coords_xi3(34) = c;
968  weights(34) = w;
969  coords_xi1(35) = a;
970  coords_xi2(35) = c;
971  coords_xi3(35) = b;
972  weights(35) = w;
973  coords_xi1(36) = b;
974  coords_xi2(36) = c;
975  coords_xi3(36) = a;
976  weights(36) = w;
977  coords_xi1(37) = c;
978  coords_xi2(37) = a;
979  coords_xi3(37) = b;
980  weights(37) = w;
981  coords_xi1(38) = c;
982  coords_xi2(38) = b;
983  coords_xi3(38) = a;
984  weights(38) = w;
985  coords_xi1(39) = a;
986  coords_xi2(39) = c;
987  coords_xi3(39) = c;
988  weights(39) = w;
989  coords_xi1(40) = b;
990  coords_xi2(40) = c;
991  coords_xi3(40) = c;
992  weights(40) = w;
993  coords_xi1(41) = c;
994  coords_xi2(41) = a;
995  coords_xi3(41) = c;
996  weights(41) = w;
997  coords_xi1(42) = c;
998  coords_xi2(42) = b;
999  coords_xi3(42) = c;
1000  weights(42) = w;
1001  coords_xi1(43) = c;
1002  coords_xi2(43) = c;
1003  coords_xi3(43) = a;
1004  weights(43) = w;
1005  coords_xi1(44) = c;
1006  coords_xi2(44) = c;
1007  coords_xi3(44) = b;
1008  weights(44) = w;
1009 
1010  break;
1011 
1012  default:
1013  OOFEM_SERROR("unsupported number of IPs (%d)", nPoints);
1014  }
1015 }
1016 
1017 void
1018 GaussIntegrationRule :: giveTriCoordsAndWeights(int nPoints, FloatArray &coords_xi1, FloatArray &coords_xi2, FloatArray &weights)
1019 // Create arrays of coordinates and weights for Gauss Integration Points of a triangle with 'nPoints' integrationpoints
1020 // Coordinates xi1 and xi2 are the two first area coordinates of the triangle.
1021 // Dunavant quadrature rules for triangles in area coordinates. Taken from http://www.mems.rice.edu/~akin/Elsevier/Chap_10.pdf
1022 {
1023  switch ( nPoints ) {
1024  case 1:
1025  coords_xi1 = FloatArray{0.333333333333};
1026  coords_xi2 = FloatArray{0.333333333333};
1027  weights = FloatArray{0.5};
1028  break;
1029 
1030  case 3:
1031  coords_xi1 = {
1032  0.166666666666667,
1033  0.666666666666667,
1034  0.166666666666667
1035  };
1036  coords_xi2 = {
1037  0.166666666666667,
1038  0.166666666666667,
1039  0.666666666666667
1040  };
1041  weights = {
1042  0.166666666666666,
1043  0.166666666666666,
1044  0.166666666666666
1045  };
1046  break;
1047 
1048  case 4:
1049 
1050  coords_xi1 = {
1051  0.333333333333333,
1052  0.200000000000000,
1053  0.200000000000000,
1054  0.600000000000000
1055  };
1056  coords_xi2 = {
1057  0.333333333333333,
1058  0.600000000000000,
1059  0.200000000000000,
1060  0.200000000000000
1061  };
1062  weights = {
1063  -0.281250000000000,
1064  0.260416666666667,
1065  0.260416666666667,
1066  0.260416666666667
1067  };
1068  break;
1069 
1070  case 6:
1071  coords_xi1 = {
1072  0.445948490915965,
1073  0.445948490915965,
1074  0.108103018168070,
1075  0.091576213509771,
1076  0.091576213509771,
1077  0.816847572980459
1078  };
1079  coords_xi2 = {
1080  0.108103018168070,
1081  0.445948490915965,
1082  0.445948490915965,
1083  0.816847572980459,
1084  0.091576213509771,
1085  0.091576213509771
1086  };
1087  weights = {
1088  0.111690794839006,
1089  0.111690794839006,
1090  0.111690794839006,
1091  0.054975871827661,
1092  0.054975871827661,
1093  0.054975871827661
1094  };
1095  break;
1096 
1097  case 7:
1098  coords_xi1 = {
1099  0.333333333333333,
1100  0.470142064105115,
1101  0.470142064105115,
1102  0.059715871789770,
1103  0.101286507323456,
1104  0.101286507323456,
1105  0.797426985353087
1106  };
1107  coords_xi2 = {
1108  0.333333333333333,
1109  0.059715871789770,
1110  0.470142064105115,
1111  0.470142064105115,
1112  0.797426985353087,
1113  0.101286507323456,
1114  0.101286507323456
1115  };
1116  weights = {
1117  0.112500000000000,
1118  0.066197076394253,
1119  0.066197076394253,
1120  0.066197076394253,
1121  0.062969590272414,
1122  0.062969590272414,
1123  0.062969590272414
1124  };
1125  break;
1126 
1127  case 12:
1128  coords_xi1 = {
1129  0.249286745170910,
1130  0.249286745170910,
1131  0.501426509658179,
1132  0.063089014491502,
1133  0.063089014491502,
1134  0.873821971016996,
1135  0.636502499121399,
1136  0.310352451033784,
1137  0.053145049844817,
1138  0.636502499121399,
1139  0.310352451033784,
1140  0.053145049844817
1141  };
1142  coords_xi2 = {
1143  0.501426509658179,
1144  0.249286745170910,
1145  0.249286745170910,
1146  0.873821971016996,
1147  0.063089014491502,
1148  0.063089014491502,
1149  0.053145049844817,
1150  0.636502499121399,
1151  0.310352451033784,
1152  0.053145049844817,
1153  0.636502499121399,
1154  0.310352451033784
1155  };
1156  weights = {
1157  0.058393137863189,
1158  0.058393137863189,
1159  0.058393137863189,
1160  0.025422453185104,
1161  0.025422453185104,
1162  0.025422453185104,
1163  0.041425537809187,
1164  0.041425537809187,
1165  0.041425537809187,
1166  0.041425537809187,
1167  0.041425537809187,
1168  0.041425537809187
1169  };
1170  break;
1171 
1172  case 13:
1173  coords_xi1 = {
1174  0.333333333333333,
1175  0.260345966079040,
1176  0.260345966079040,
1177  0.479308067841920,
1178  0.065130102902216,
1179  0.065130102902216,
1180  0.869739794195568,
1181  0.638444188569810,
1182  0.312865496004874,
1183  0.048690315425316,
1184  0.638444188569810,
1185  0.312865496004874,
1186  0.048690315425316
1187  };
1188  coords_xi2 = {
1189  0.333333333333333,
1190  0.479308067841920,
1191  0.260345966079040,
1192  0.260345966079040,
1193  0.869739794195568,
1194  0.065130102902216,
1195  0.065130102902216,
1196  0.048690315425316,
1197  0.638444188569810,
1198  0.312865496004874,
1199  0.048690315425316,
1200  0.638444188569810,
1201  0.312865496004874
1202  };
1203  weights = {
1204  -0.074785022233841,
1205  0.087807628716604,
1206  0.087807628716604,
1207  0.087807628716604,
1208  0.026673617804419,
1209  0.026673617804419,
1210  0.026673617804419,
1211  0.038556880445128,
1212  0.038556880445128,
1213  0.038556880445128,
1214  0.038556880445128,
1215  0.038556880445128,
1216  0.038556880445128
1217  };
1218  break;
1219 
1220  case 16:
1221  coords_xi1 = {
1222  0.333333333333333,
1223  0.459292588292723,
1224  0.459292588292723,
1225  0.081414823414554,
1226  0.170569307751760,
1227  0.170569307751760,
1228  0.658861384496480,
1229  0.050547228317031,
1230  0.050547228317031,
1231  0.898905543365938,
1232  0.728492392955404,
1233  0.263112829634638,
1234  0.008394777409958,
1235  0.728492392955404,
1236  0.263112829634638,
1237  0.008394777409958
1238  };
1239  coords_xi2 = {
1240  0.333333333333333,
1241  0.081414823414554,
1242  0.459292588292723,
1243  0.459292588292723,
1244  0.658861384496480,
1245  0.170569307751760,
1246  0.170569307751760,
1247  0.898905543365938,
1248  0.050547228317031,
1249  0.050547228317031,
1250  0.008394777409958,
1251  0.728492392955404,
1252  0.263112829634638,
1253  0.008394777409958,
1254  0.728492392955404,
1255  0.263112829634638
1256  };
1257  weights = {
1258  0.072157803838894,
1259  0.047545817133642,
1260  0.047545817133642,
1261  0.047545817133642,
1262  0.051608685267359,
1263  0.051608685267359,
1264  0.051608685267359,
1265  0.016229248811599,
1266  0.016229248811599,
1267  0.016229248811599,
1268  0.013615157087218,
1269  0.013615157087218,
1270  0.013615157087218,
1271  0.013615157087218,
1272  0.013615157087218,
1273  0.013615157087218
1274  };
1275  break;
1276 
1277 
1278  case 19:
1279  coords_xi1 = {
1280  0.333333333333333,
1281  0.489682519198738,
1282  0.489682519198738,
1283  0.020634961602525,
1284  0.437089591492937,
1285  0.437089591492937,
1286  0.125820817014127,
1287  0.188203535619033,
1288  0.188203535619033,
1289  0.623592928761935,
1290  0.044729513394453,
1291  0.044729513394453,
1292  0.910540973211095,
1293  0.741198598784498,
1294  0.221962989160766,
1295  0.036838412054736,
1296  0.741198598784498,
1297  0.221962989160766,
1298  0.036838412054736
1299  };
1300  coords_xi2 = {
1301  0.333333333333333,
1302  0.020634961602525,
1303  0.489682519198738,
1304  0.489682519198738,
1305  0.125820817014127,
1306  0.437089591492937,
1307  0.437089591492937,
1308  0.623592928761935,
1309  0.188203535619033,
1310  0.188203535619033,
1311  0.910540973211095,
1312  0.044729513394453,
1313  0.044729513394453,
1314  0.036838412054736,
1315  0.741198598784498,
1316  0.221962989160766,
1317  0.036838412054736,
1318  0.741198598784498,
1319  0.221962989160766
1320  };
1321  weights = {
1322  0.048567898141400,
1323  0.015667350113570,
1324  0.015667350113570,
1325  0.015667350113570,
1326  0.038913770502387,
1327  0.038913770502387,
1328  0.038913770502387,
1329  0.039823869463605,
1330  0.039823869463605,
1331  0.039823869463605,
1332  0.012788837829349,
1333  0.012788837829349,
1334  0.012788837829349,
1335  0.021641769688645,
1336  0.021641769688645,
1337  0.021641769688645,
1338  0.021641769688645,
1339  0.021641769688645,
1340  0.021641769688645
1341  };
1342  break;
1343 
1344  case 25:
1345  coords_xi1 = {
1346  0.333333333333333,
1347  0.485577633383657,
1348  0.485577633383657,
1349  0.028844733232685,
1350  0.109481575485037,
1351  0.109481575485037,
1352  0.781036849029926,
1353  0.550352941820999,
1354  0.307939838764121,
1355  0.141707219414880,
1356  0.550352941820999,
1357  0.307939838764121,
1358  0.141707219414880,
1359  0.728323904597411,
1360  0.246672560639903,
1361  0.025003534762686,
1362  0.728323904597411,
1363  0.246672560639903,
1364  0.025003534762686,
1365  0.923655933587500,
1366  0.066803251012200,
1367  0.009540815400299,
1368  0.923655933587500,
1369  0.066803251012200,
1370  0.009540815400299
1371  };
1372  coords_xi2 = {
1373  0.333333333333333,
1374  0.028844733232685,
1375  0.485577633383657,
1376  0.485577633383657,
1377  0.781036849029926,
1378  0.109481575485037,
1379  0.109481575485037,
1380  0.141707219414880,
1381  0.550352941820999,
1382  0.307939838764121,
1383  0.141707219414880,
1384  0.550352941820999,
1385  0.307939838764121,
1386  0.025003534762686,
1387  0.728323904597411,
1388  0.246672560639903,
1389  0.025003534762686,
1390  0.728323904597411,
1391  0.246672560639903,
1392  0.009540815400299,
1393  0.923655933587500,
1394  0.066803251012200,
1395  0.009540815400299,
1396  0.923655933587500,
1397  0.066803251012200
1398  };
1399  weights = {
1400  0.045408995191377,
1401  0.018362978878233,
1402  0.018362978878233,
1403  0.018362978878233,
1404  0.022660529717764,
1405  0.022660529717764,
1406  0.022660529717764,
1407  0.036378958422710,
1408  0.036378958422710,
1409  0.036378958422710,
1410  0.036378958422710,
1411  0.036378958422710,
1412  0.036378958422710,
1413  0.014163621265529,
1414  0.014163621265529,
1415  0.014163621265529,
1416  0.014163621265529,
1417  0.014163621265529,
1418  0.014163621265529,
1419  0.004710833481867,
1420  0.004710833481867,
1421  0.004710833481867,
1422  0.004710833481867,
1423  0.004710833481867,
1424  0.004710833481867
1425  };
1426  break;
1427 
1428  default:
1429  OOFEM_SERROR("unsupported number of IPs (%d)", nPoints);
1430  }
1431 }
1432 
1433 
1434 
1435 void
1436 GaussIntegrationRule :: giveLineCoordsAndWeights(int nPoints, FloatArray &coords_xi, FloatArray &weights)
1437 // Create arrays of coordinates and weights for Gauss Integration Points of a line with 'nPoints' integrationpoints
1438 {
1439  switch ( nPoints ) {
1440  case 0:
1441  coords_xi = FloatArray{};
1442  weights = FloatArray{};
1443  break;
1444  case 1:
1445  coords_xi = FloatArray{0.0};
1446  weights = FloatArray{2.0};
1447  break;
1448 
1449  case 2:
1450  coords_xi = {-0.577350269189626, 0.577350269189626};
1451  weights = {1.0, 1.0};
1452  break;
1453 
1454  case 3:
1455  coords_xi = {
1456  -0.774596669241483,
1457  0.0,
1458  0.774596669241483
1459  };
1460  weights = {
1461  0.555555555555555,
1462  0.888888888888888,
1463  0.555555555555555
1464  };
1465  break;
1466 
1467  case 4:
1468  coords_xi = {
1469  -0.861136311594053,
1470  -0.339981043584856,
1471  0.339981043584856,
1472  0.861136311594053
1473  };
1474  weights = {
1475  0.347854845137454,
1476  0.652145154862546,
1477  0.652145154862546,
1478  0.347854845137454
1479  };
1480  break;
1481 
1482  case 5:
1483  coords_xi = {
1484  -0.9061798459386639927976269,
1485  -0.5384693101056830910363144,
1486  0.0,
1487  0.5384693101056830910363144,
1488  0.9061798459386639927976269
1489  };
1490  weights = {
1491  0.2369268850561890875142640,
1492  0.4786286704993664680412915,
1493  0.5688888888888888888888889,
1494  0.4786286704993664680412915,
1495  0.2369268850561890875142640
1496  };
1497  break;
1498 
1499  case 6:
1500  coords_xi = {
1501  -0.2386191860831969086305017,
1502  -0.6612093864662645136613996,
1503  -0.9324695142031520278123016,
1504  0.9324695142031520278123016,
1505  0.6612093864662645136613996,
1506  0.2386191860831969086305017
1507  };
1508  weights = {
1509  0.4679139345726910473898703,
1510  0.3607615730481386075698335,
1511  0.1713244923791703450402961,
1512  0.1713244923791703450402961,
1513  0.3607615730481386075698335,
1514  0.4679139345726910473898703
1515  };
1516  break;
1517 
1518  case 7:
1519  coords_xi = {
1520  -0.9491079123427585245261897,
1521  -0.7415311855993944398638648,
1522  -0.4058451513773971669066064,
1523  0.0,
1524  0.4058451513773971669066064,
1525  0.7415311855993944398638648,
1526  0.9491079123427585245261897
1527  };
1528  weights = {
1529  0.1294849661688696932706114,
1530  0.2797053914892766679014678,
1531  0.3818300505051189449503698,
1532  0.4179591836734693877551020,
1533  0.3818300505051189449503698,
1534  0.2797053914892766679014678,
1535  0.1294849661688696932706114
1536  };
1537  break;
1538 
1539  case 8:
1540  coords_xi = {
1541  -0.960289856497536,
1542  -0.796666477413627,
1543  -0.525532409916329,
1544  -0.183434642495650,
1545  0.183434642495650,
1546  0.525532409916329,
1547  0.796666477413627,
1548  0.960289856497536
1549  };
1550  weights = {
1551  0.101228536290375,
1552  0.222381034453374,
1553  0.313706645877887,
1554  0.362683783378362,
1555  0.362683783378362,
1556  0.313706645877887,
1557  0.222381034453374,
1558  0.101228536290375
1559  };
1560  break;
1561 
1562 
1563  case 16:
1564  coords_xi = {
1565  -0.989400934991650,
1566  -0.944575023073233,
1567  -0.865631202387832,
1568  -0.755404408355003,
1569  -0.617876244402644,
1570  -0.458016777657227,
1571  -0.281603550779259,
1572  -0.095012509837637,
1573  0.095012509837637,
1574  0.281603550779259,
1575  0.458016777657227,
1576  0.617876244402644,
1577  0.755404408355003,
1578  0.865631202387832,
1579  0.944575023073233,
1580  0.989400934991650
1581  };
1582  weights = {
1583  0.027152459411753,
1584  0.062253523938647,
1585  0.095158511682492,
1586  0.124628971255534,
1587  0.149595988816577,
1588  0.169156519395003,
1589  0.182603415044924,
1590  0.189450610455068,
1591  0.189450610455068,
1592  0.182603415044924,
1593  0.169156519395003,
1594  0.149595988816577,
1595  0.124628971255534,
1596  0.095158511682492,
1597  0.062253523938647,
1598  0.027152459411753
1599  };
1600  break;
1601 
1602  case 24:
1603  coords_xi = {
1604  -0.995187219997021,
1605  -0.974728555971309,
1606  -0.938274552002733,
1607  -0.886415527004401,
1608  -0.820001985973903,
1609  -0.740124191578554,
1610  -0.648093651936976,
1611  -0.545421471388840,
1612  -0.433793507626045,
1613  -0.315042679696163,
1614  -0.191118867473616,
1615  -0.064056892862606,
1616  0.064056892862606,
1617  0.191118867473616,
1618  0.315042679696163,
1619  0.433793507626045,
1620  0.545421471388840,
1621  0.648093651936976,
1622  0.740124191578554,
1623  0.820001985973903,
1624  0.886415527004401,
1625  0.938274552002733,
1626  0.974728555971309,
1627  0.995187219997021
1628  };
1629  weights = {
1630  0.012341229799991,
1631  0.028531388628934,
1632  0.044277438817419,
1633  0.059298584915436,
1634  0.073346481411080,
1635  0.086190161531953,
1636  0.097618652104114,
1637  0.107444270115966,
1638  0.115505668053726,
1639  0.121670472927803,
1640  0.125837456346828,
1641  0.127938195346752,
1642  0.127938195346752,
1643  0.125837456346828,
1644  0.121670472927803,
1645  0.115505668053726,
1646  0.107444270115966,
1647  0.097618652104114,
1648  0.086190161531953,
1649  0.073346481411080,
1650  0.059298584915436,
1651  0.044277438817419,
1652  0.028531388628934,
1653  0.012341229799991
1654  };
1655  break;
1656 
1657  case 32:
1658  coords_xi = {
1659  -0.997263861849482,
1660  -0.985611511545268,
1661  -0.964762255587506,
1662  -0.934906075937740,
1663  -0.896321155766052,
1664  -0.849367613732570,
1665  -0.794483795967942,
1666  -0.732182118740290,
1667  -0.663044266930215,
1668  -0.587715757240762,
1669  -0.506899908932229,
1670  -0.421351276130635,
1671  -0.331868602282128,
1672  -0.239287362252137,
1673  -0.144471961582796,
1674  -0.048307665687738,
1675  0.048307665687738,
1676  0.144471961582796,
1677  0.239287362252137,
1678  0.331868602282128,
1679  0.421351276130635,
1680  0.506899908932229,
1681  0.587715757240762,
1682  0.663044266930215,
1683  0.732182118740290,
1684  0.794483795967942,
1685  0.849367613732570,
1686  0.896321155766052,
1687  0.934906075937740,
1688  0.964762255587506,
1689  0.985611511545268,
1690  0.997263861849482
1691  };
1692  weights = {
1693  0.007018610009469,
1694  0.016274394730905,
1695  0.025392065309263,
1696  0.034273862913022,
1697  0.042835898022227,
1698  0.050998059262376,
1699  0.058684093478536,
1700  0.065822222776362,
1701  0.072345794108848,
1702  0.078193895787070,
1703  0.083311924226947,
1704  0.087652093004404,
1705  0.091173878695764,
1706  0.093844399080805,
1707  0.095638720079275,
1708  0.096540088514728,
1709  0.096540088514728,
1710  0.095638720079275,
1711  0.093844399080805,
1712  0.091173878695764,
1713  0.087652093004404,
1714  0.083311924226947,
1715  0.078193895787070,
1716  0.072345794108848,
1717  0.065822222776362,
1718  0.058684093478536,
1719  0.050998059262376,
1720  0.042835898022227,
1721  0.034273862913022,
1722  0.025392065309263,
1723  0.016274394730905,
1724  0.007018610009469
1725  };
1726  break;
1727 
1728  case 64:
1729  coords_xi = {
1730  -0.999305041735772,
1731  -0.996340116771955,
1732  -0.991013371476744,
1733  -0.983336253884626,
1734  -0.973326827789911,
1735  -0.961008799652054,
1736  -0.946411374858403,
1737  -0.929569172131940,
1738  -0.910522137078503,
1739  -0.889315445995114,
1740  -0.865999398154093,
1741  -0.840629296252580,
1742  -0.813265315122798,
1743  -0.783972358943341,
1744  -0.752819907260532,
1745  -0.719881850171611,
1746  -0.685236313054233,
1747  -0.648965471254657,
1748  -0.611155355172393,
1749  -0.571895646202634,
1750  -0.531279464019895,
1751  -0.489403145707053,
1752  -0.446366017253464,
1753  -0.402270157963992,
1754  -0.357220158337668,
1755  -0.311322871990211,
1756  -0.264687162208767,
1757  -0.217423643740007,
1758  -0.169644420423993,
1759  -0.121462819296121,
1760  -0.072993121787799,
1761  -0.024350292663424,
1762  0.024350292663424,
1763  0.072993121787799,
1764  0.121462819296121,
1765  0.169644420423993,
1766  0.217423643740007,
1767  0.264687162208767,
1768  0.311322871990211,
1769  0.357220158337668,
1770  0.402270157963992,
1771  0.446366017253464,
1772  0.489403145707053,
1773  0.531279464019895,
1774  0.571895646202634,
1775  0.611155355172393,
1776  0.648965471254657,
1777  0.685236313054233,
1778  0.719881850171611,
1779  0.752819907260532,
1780  0.783972358943341,
1781  0.813265315122798,
1782  0.840629296252580,
1783  0.865999398154093,
1784  0.889315445995114,
1785  0.910522137078503,
1786  0.929569172131940,
1787  0.946411374858403,
1788  0.961008799652054,
1789  0.973326827789911,
1790  0.983336253884626,
1791  0.991013371476744,
1792  0.996340116771955,
1793  0.999305041735772
1794  };
1795  weights = {
1796  0.001783280721693,
1797  0.004147033260559,
1798  0.006504457968980,
1799  0.008846759826363,
1800  0.011168139460130,
1801  0.013463047896718,
1802  0.015726030476025,
1803  0.017951715775697,
1804  0.020134823153530,
1805  0.022270173808383,
1806  0.024352702568711,
1807  0.026377469715055,
1808  0.028339672614260,
1809  0.030234657072403,
1810  0.032057928354851,
1811  0.033805161837142,
1812  0.035472213256882,
1813  0.037055128540240,
1814  0.038550153178615,
1815  0.039953741132721,
1816  0.041262563242623,
1817  0.042473515123654,
1818  0.043583724529323,
1819  0.044590558163757,
1820  0.045491627927418,
1821  0.046284796581314,
1822  0.046968182816210,
1823  0.047540165714830,
1824  0.047999388596458,
1825  0.048344762234803,
1826  0.048575467441503,
1827  0.048690957009140,
1828  0.048690957009140,
1829  0.048575467441503,
1830  0.048344762234803,
1831  0.047999388596458,
1832  0.047540165714830,
1833  0.046968182816210,
1834  0.046284796581314,
1835  0.045491627927418,
1836  0.044590558163757,
1837  0.043583724529323,
1838  0.042473515123654,
1839  0.041262563242623,
1840  0.039953741132721,
1841  0.038550153178615,
1842  0.037055128540240,
1843  0.035472213256882,
1844  0.033805161837142,
1845  0.032057928354851,
1846  0.030234657072403,
1847  0.028339672614260,
1848  0.026377469715055,
1849  0.024352702568711,
1850  0.022270173808383,
1851  0.020134823153530,
1852  0.017951715775697,
1853  0.015726030476025,
1854  0.013463047896718,
1855  0.011168139460130,
1856  0.008846759826363,
1857  0.006504457968980,
1858  0.004147033260559,
1859  0.001783280721693
1860  };
1861  break;
1862 
1863  default:
1864  OOFEM_SERROR("unsupported number of IPs (%d)", nPoints);
1865  }
1866 }
1867 } // end namespace oofem
oofem::integrationDomain
integrationDomain
Used by integrator class to supply integration points for proper domain to be integrated (Area...
Definition: integrationdomain.h:44
oofem::IntegrationRule::intdomain
integrationDomain intdomain
Integration domain.
Definition: integrationrule.h:102
oofem::GaussIntegrationRule::giveTriCoordsAndWeights
static void giveTriCoordsAndWeights(int nPoints, FloatArray &coords_xi1, FloatArray &coords_xi2, FloatArray &weights)
Definition: gaussintegrationrule.C:1018
oofem::GaussIntegrationRule::giveLineCoordsAndWeights
static void giveLineCoordsAndWeights(int nPoints, FloatArray &coords_xi, FloatArray &weights)
Definition: gaussintegrationrule.C:1436
OOFEM_SERROR
#define OOFEM_SERROR(...)
Definition: error.h:63
oofem::FloatArray::at
double & at(int i)
Coefficient access function.
Definition: floatarray.h:131
oofem::max
int max(int i, int j)
Returns bigger value form two given decimals.
Definition: mathfem.h:71
oofem::GaussIntegrationRule::SetUpPointsOn3dDegShellLayers
virtual int SetUpPointsOn3dDegShellLayers(int nPointsXY, int nPointsZ, MaterialMode mode, const FloatArray &layerThickness)
Sets up receiver&#39;s integration points on shell integration domain wih layers.
Definition: gaussintegrationrule.C:153
oofem::IntegrationRule::giveElement
Element * giveElement()
Returns reference to element containing receiver.
Definition: integrationrule.h:203
oofem::Element
Abstract base class for all finite elements.
Definition: element.h:145
oofem::MaterialMode
MaterialMode
Type representing material mode of integration point.
Definition: materialmode.h:89
oofem::IntegrationRule
Abstract base class representing integration rule.
Definition: integrationrule.h:94
oofem::GaussIntegrationRule::~GaussIntegrationRule
virtual ~GaussIntegrationRule()
Destructor.
Definition: gaussintegrationrule.C:51
oofem::GaussIntegrationRule::SetUpPointsOnSquare
virtual int SetUpPointsOnSquare(int nPoints, MaterialMode mode)
Sets up receiver&#39;s integration points on unit square integration domain.
Definition: gaussintegrationrule.C:102
oofem::FloatArray::sum
double sum() const
Computes the sum of receiver values.
Definition: floatarray.C:852
OOFEM_ERROR
#define OOFEM_ERROR(...)
Definition: error.h:61
oofem::GaussIntegrationRule::SetUpPointsOnTetrahedra
virtual int SetUpPointsOnTetrahedra(int nPoints, MaterialMode mode)
Sets up receiver&#39;s integration points on tetrahedra (volume coords) integration domain.
Definition: gaussintegrationrule.C:269
oofem::GaussIntegrationRule::SetUpPointsOnCube
virtual int SetUpPointsOnCube(int nPoints, MaterialMode mode)
Sets up receiver&#39;s integration points on unit cube integration domain.
Definition: gaussintegrationrule.C:190
oofem::GaussIntegrationRule::getRequiredNumberOfIntegrationPoints
virtual int getRequiredNumberOfIntegrationPoints(integrationDomain dType, int approxOrder)
Abstract service.
Definition: gaussintegrationrule.C:339
oofem::FloatArray
Class representing vector of real numbers.
Definition: floatarray.h:82
oofem::cbrt
double cbrt(double x)
Returns the cubic root of x.
Definition: mathfem.h:109
oofem::GaussIntegrationRule::SetUpPointsOnWedgeLayers
virtual int SetUpPointsOnWedgeLayers(int nPointsTri, int nPointsDepth, MaterialMode mode, const FloatArray &layerThickness)
Sets up receiver&#39;s integration points on a wedge integration domain divided into layers in the zeta-d...
Definition: gaussintegrationrule.C:307
oofem::IntegrationRule::giveNumberOfIntegrationPoints
int giveNumberOfIntegrationPoints() const
Returns number of integration points of receiver.
Definition: integrationrule.h:143
oofem::GaussIntegrationRule::SetUpPointsOn2DEmbeddedLine
virtual int SetUpPointsOn2DEmbeddedLine(int nPoints, MaterialMode mode, const FloatArray &coord0, const FloatArray &coord1)
Sets up integration points on 2D embedded line inside 2D volume (the list of local coordinates should...
Definition: gaussintegrationrule.C:72
oofem::Element::computeLocalCoordinates
virtual bool computeLocalCoordinates(FloatArray &answer, const FloatArray &gcoords)
Computes the element local coordinates from given global coordinates.
Definition: element.C:1222
oofem::FloatArray::giveSize
int giveSize() const
Returns the size of receiver.
Definition: floatarray.h:218
oofem
the oofem namespace is to define a context or scope in which all oofem names are defined.
oofem::GaussIntegrationRule::SetUpPointsOnWedge
virtual int SetUpPointsOnWedge(int nPointsTri, int nPointsDepth, MaterialMode mode)
Sets up receiver&#39;s integration points on a wedge integration domain.
Definition: gaussintegrationrule.C:286
oofem::GaussIntegrationRule::SetUpPointsOnLine
virtual int SetUpPointsOnLine(int nPoints, MaterialMode mode)
Sets up receiver&#39;s integration points on unit line integration domain.
Definition: gaussintegrationrule.C:56
oofem::GaussIntegrationRule::SetUpPointsOn3dDegShell
virtual int SetUpPointsOn3dDegShell(int nPointsXY, int nPointsZ, MaterialMode mode)
Sets up receiver&#39;s integration points on shell integration domain.
Definition: gaussintegrationrule.C:125
oofem::GaussPoint
Class representing integration point in finite element program.
Definition: gausspoint.h:93
oofem::IntegrationRule::gaussPoints
std::vector< GaussPoint * > gaussPoints
Array containing integration points.
Definition: integrationrule.h:105
oofem::GaussIntegrationRule::SetUpPointsOnCubeLayers
virtual int SetUpPointsOnCubeLayers(int nPoints1, int nPoints2, int nPointsDepth, MaterialMode mode, const FloatArray &layerThickness)
Sets up receiver&#39;s integration points on unit cube integration domain divided into layers in the zeta...
Definition: gaussintegrationrule.C:217
oofem::GaussIntegrationRule::SetUpPointsOnTriangle
virtual int SetUpPointsOnTriangle(int nPoints, MaterialMode mode)
Sets up receiver&#39;s integration points on triangular (area coords) integration domain.
Definition: gaussintegrationrule.C:253
oofem::GaussIntegrationRule::giveTetCoordsAndWeights
static void giveTetCoordsAndWeights(int nPoints, FloatArray &coords_xi1, FloatArray &coords_xi2, FloatArray &coords_xi3, FloatArray &weights)
Definition: gaussintegrationrule.C:492
oofem::GaussIntegrationRule::GaussIntegrationRule
GaussIntegrationRule(int n, Element *e, int startIndx, int endIndx, bool dynamic=false)
Constructor.
Definition: gaussintegrationrule.C:44
oofem::FloatArray::resize
void resize(int s)
Resizes receiver towards requested size.
Definition: floatarray.C:631
This page is part of the OOFEM documentation. Copyright (c) 2011 Borek Patzak
Project e-mail: info@oofem.org
Generated at Tue Jan 2 2018 20:07:28 for OOFEM by doxygen 1.8.11 written by Dimitri van Heesch, © 1997-2011