OOFEM 3.0
Loading...
Searching...
No Matches
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 - 2025 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
41namespace oofem {
42// initialize class member
43
44GaussIntegrationRule :: GaussIntegrationRule(int n, Element *e,
45 int startIndx, int endIndx, bool dynamic) :
46 IntegrationRule(n, e, startIndx, endIndx, dynamic) { }
47
48GaussIntegrationRule :: GaussIntegrationRule(int n, Element *e) :
49 IntegrationRule(n, e) { }
50
51GaussIntegrationRule :: ~GaussIntegrationRule()
52{ }
53
54
55int
56GaussIntegrationRule :: 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, Vec1(coords_xi.at(i)), weights.at ( i ), mode);
64 }
65
66 this->intdomain = _Line;
67 return this->giveNumberOfIntegrationPoints();
68}
69
70
71int
72GaussIntegrationRule :: 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 = Vec1(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 = Vec2(
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
101int
102GaussIntegrationRule :: 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, Vec2(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
124int
125GaussIntegrationRule :: 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, Vec3(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
152int
153GaussIntegrationRule :: 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, Vec3(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
189int
190GaussIntegrationRule :: 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, Vec3(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
217int 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 Vec3(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
252int
253GaussIntegrationRule :: 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, Vec2(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
268int
269GaussIntegrationRule :: 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 Vec3(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
285int
286GaussIntegrationRule :: 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, Vec3(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
306int
307GaussIntegrationRule :: 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 Vec3(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
338int
339GaussIntegrationRule :: 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
491void
492GaussIntegrationRule :: 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
1017void
1018GaussIntegrationRule :: 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 = Vec1(0.333333333333);
1026 coords_xi2 = Vec1(0.333333333333);
1027 weights = Vec1(0.5);
1028 break;
1029
1030 case 3:
1031 coords_xi1 = Vec3(
1032 0.166666666666667,
1033 0.666666666666667,
1034 0.166666666666667
1035 );
1036 coords_xi2 = Vec3(
1037 0.166666666666667,
1038 0.166666666666667,
1039 0.666666666666667
1040 );
1041 weights = Vec3(
1042 0.166666666666666,
1043 0.166666666666666,
1044 0.166666666666666
1045 );
1046 break;
1047
1048 case 4:
1049
1050 coords_xi1 = Vec4(
1051 0.333333333333333,
1052 0.200000000000000,
1053 0.200000000000000,
1054 0.600000000000000
1055 );
1056 coords_xi2 = Vec4(
1057 0.333333333333333,
1058 0.600000000000000,
1059 0.200000000000000,
1060 0.200000000000000
1061 );
1062 weights = Vec4(
1063 -0.281250000000000,
1064 0.260416666666667,
1065 0.260416666666667,
1066 0.260416666666667
1067 );
1068 break;
1069
1070 case 6:
1071 coords_xi1 = Vec6(
1072 0.445948490915965,
1073 0.445948490915965,
1074 0.108103018168070,
1075 0.091576213509771,
1076 0.091576213509771,
1077 0.816847572980459
1078 );
1079 coords_xi2 = Vec6(
1080 0.108103018168070,
1081 0.445948490915965,
1082 0.445948490915965,
1083 0.816847572980459,
1084 0.091576213509771,
1085 0.091576213509771
1086 );
1087 weights = Vec6(
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 = Vec7(
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 = Vec7(
1108 0.333333333333333,
1109 0.059715871789770,
1110 0.470142064105115,
1111 0.470142064105115,
1112 0.797426985353087,
1113 0.101286507323456,
1114 0.101286507323456
1115 );
1116 weights = Vec7(
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 = VecX({
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 = VecX({
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 = VecX({
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 = VecX({
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 = VecX({
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 = VecX({
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 = VecX({
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 = VecX({
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 = VecX({
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 = VecX({
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 = VecX({
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 = VecX({
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 = VecX({
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 = VecX({
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 = VecX({
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
1435void
1436GaussIntegrationRule :: 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 = Vec1(0.0);
1446 weights = Vec1(2.0);
1447 break;
1448
1449 case 2:
1450 coords_xi = Vec2(-0.577350269189626, 0.577350269189626);
1451 weights = Vec2(1.0, 1.0);
1452 break;
1453
1454 case 3:
1455 coords_xi = Vec3(
1456 -0.774596669241483,
1457 0.0,
1458 0.774596669241483
1459 );
1460 weights = Vec3(
1461 0.555555555555555,
1462 0.888888888888888,
1463 0.555555555555555
1464 );
1465 break;
1466
1467 case 4:
1468 coords_xi = Vec4(
1469 -0.861136311594053,
1470 -0.339981043584856,
1471 0.339981043584856,
1472 0.861136311594053
1473 );
1474 weights = Vec4(
1475 0.347854845137454,
1476 0.652145154862546,
1477 0.652145154862546,
1478 0.347854845137454
1479 );
1480 break;
1481
1482 case 5:
1483 coords_xi = Vec5(
1484 -0.9061798459386639927976269,
1485 -0.5384693101056830910363144,
1486 0.0,
1487 0.5384693101056830910363144,
1488 0.9061798459386639927976269
1489 );
1490 weights = Vec5(
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 = Vec6(
1501 -0.2386191860831969086305017,
1502 -0.6612093864662645136613996,
1503 -0.9324695142031520278123016,
1504 0.9324695142031520278123016,
1505 0.6612093864662645136613996,
1506 0.2386191860831969086305017
1507 );
1508 weights = Vec6(
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 = Vec7(
1520 -0.9491079123427585245261897,
1521 -0.7415311855993944398638648,
1522 -0.4058451513773971669066064,
1523 0.0,
1524 0.4058451513773971669066064,
1525 0.7415311855993944398638648,
1526 0.9491079123427585245261897
1527 );
1528 weights = Vec7(
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 = Vec8(
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 = Vec8(
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 = VecX({
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 = VecX({
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 = VecX({
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 = VecX({
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 = VecX({
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 = VecX({
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 = VecX({
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 = VecX({
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::FloatArray::resize
void resize(Index s)
Definition floatarray.C:94
oofem::FloatArray::at
double & at(Index i)
Definition floatarray.h:202
oofem::FloatArray::giveSize
Index giveSize() const
Returns the size of receiver.
Definition floatarray.h:261
oofem::FloatArray::sum
double sum() const
Definition floatarray.C:874
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::GaussIntegrationRule::giveTetCoordsAndWeights
static void giveTetCoordsAndWeights(int nPoints, FloatArray &coords_xi1, FloatArray &coords_xi2, FloatArray &coords_xi3, FloatArray &weights)
Definition gaussintegrationrule.C:492
oofem::IntegrationRule::giveNumberOfIntegrationPoints
int giveNumberOfIntegrationPoints() const
Definition integrationrule.h:143
oofem::IntegrationRule::gaussPoints
std::vector< GaussPoint * > gaussPoints
Array containing integration points.
Definition integrationrule.h:105
oofem::IntegrationRule::intdomain
integrationDomain intdomain
Integration domain.
Definition integrationrule.h:102
oofem::IntegrationRule::IntegrationRule
IntegrationRule(int n, Element *e, int startIndx, int endIndx, bool dynamic)
Definition integrationrule.C:44
OOFEM_SERROR
#define OOFEM_SERROR(...)
Definition error.h:81
OOFEM_ERROR
#define OOFEM_ERROR(...)
Definition error.h:79
oofem::Vec2
static FloatArray Vec2(const double &a, const double &b)
Definition floatarray.h:606
oofem::VecX
static FloatArray VecX(std::initializer_list< double > ini)
Definition floatarray.h:614
oofem::cbrt
double cbrt(double x)
Returns the cubic root of x.
Definition mathfem.h:122
oofem::Vec6
static FloatArray Vec6(const double &a, const double &b, const double &c, const double &d, const double &e, const double &f)
Definition floatarray.h:610
oofem::max
FloatArrayF< N > max(const FloatArrayF< N > &a, const FloatArrayF< N > &b)
Definition floatarrayf.h:512
oofem::Vec5
static FloatArray Vec5(const double &a, const double &b, const double &c, const double &d, const double &e)
Definition floatarray.h:609
oofem::Vec8
static FloatArray Vec8(const double &a, const double &b, const double &c, const double &d, const double &e, const double &f, const double &g, const double &h)
Definition floatarray.h:612
oofem::Vec4
static FloatArray Vec4(const double &a, const double &b, const double &c, const double &d)
Definition floatarray.h:608
oofem::Vec7
static FloatArray Vec7(const double &a, const double &b, const double &c, const double &d, const double &e, const double &f, const double &g)
Definition floatarray.h:611
oofem::Vec1
static FloatArray Vec1(const double &a)
Definition floatarray.h:605
oofem::Vec3
static FloatArray Vec3(const double &a, const double &b, const double &c)
Definition floatarray.h:607
This page is part of the OOFEM-3.0 documentation. Copyright Copyright (C) 1994-2025 Borek Patzak Bořek Patzák
Project e-mail: oofem@fsv.cvut.cz
Generated at for OOFEM by doxygen 1.15.0 written by Dimitri van Heesch, © 1997-2011