1 // Copyright (C) 2007-2013 CEA/DEN, EDF R&D, OPEN CASCADE
3 // This library is free software; you can redistribute it and/or
4 // modify it under the terms of the GNU Lesser General Public
5 // License as published by the Free Software Foundation; either
6 // version 2.1 of the License.
8 // This library is distributed in the hope that it will be useful,
9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
11 // Lesser General Public License for more details.
13 // You should have received a copy of the GNU Lesser General Public
14 // License along with this library; if not, write to the Free Software
15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
19 // File : MED_GaussDef.hxx
20 // Author : Edward AGAPOV (eap)
22 #include "MED_GaussDef.hxx"
23 #include "MED_Utilities.hxx"
24 #include "MED_GaussUtils.hxx"
30 //---------------------------------------------------------------
32 void TGaussDef::add(const double x, const double weight)
35 EXCEPTION( logic_error,"dim() != 1");
36 if ( myWeights.capacity() == myWeights.size() )
37 EXCEPTION( logic_error,"Extra gauss point");
38 myCoords.push_back( x );
39 myWeights.push_back( weight );
41 void TGaussDef::add(const double x, const double y, const double weight)
44 EXCEPTION( logic_error,"dim() != 2");
45 if ( myWeights.capacity() == myWeights.size() )
46 EXCEPTION( logic_error,"Extra gauss point");
47 myCoords.push_back( x );
48 myCoords.push_back( y );
49 myWeights.push_back( weight );
51 void TGaussDef::add(const double x, const double y, const double z, const double weight)
54 EXCEPTION( logic_error,"dim() != 3");
55 if ( myWeights.capacity() == myWeights.size() )
56 EXCEPTION( logic_error,"Extra gauss point");
57 myCoords.push_back( x );
58 myCoords.push_back( y );
59 myCoords.push_back( z );
60 myWeights.push_back( weight );
62 void TGaussDef::setRefCoords(const TShapeFun& aShapeFun)
64 myRefCoords.reserve( aShapeFun.myRefCoord.size() );
65 myRefCoords.assign( aShapeFun.myRefCoord.begin(),
66 aShapeFun.myRefCoord.end() );
70 //---------------------------------------------------------------
72 * \brief Fill definition of gauss points family
74 //---------------------------------------------------------------
76 TGaussDef::TGaussDef(const int geom, const int nbGauss, const int variant)
79 myCoords .reserve( nbGauss * dim() );
80 myWeights.reserve( nbGauss );
86 if (geom == eSEG2) setRefCoords( TSeg2a() );
87 else setRefCoords( TSeg3a() );
90 add( 0.0, 2.0 ); break;
93 const double a = 0.577350269189626;
95 add( -a, 1.0 ); break;
98 const double a = 0.774596669241;
99 const double P1 = 1./1.8;
100 const double P2 = 1./1.125;
106 const double a = 0.339981043584856, b = 0.861136311594053;
107 const double P1 = 0.652145154862546, P2 = 0.347854845137454 ;
111 add( -b, P2 ); break;
114 EXCEPTION( logic_error,"Invalid nb of gauss points for SEG"<<nbGauss);
120 if ( variant == 1 ) {
121 if (geom == eTRIA3) setRefCoords( TTria3b() );
122 else setRefCoords( TTria6b() );
125 add( 1/3., 1/3., 1/2. ); break;
128 // what about COT3 ???
129 add( 1/6., 1/6., 1/6. );
130 add( 2/3., 1/6., 1/6. );
131 add( 1/6., 2/3., 1/6. ); break;
134 add( 1/5., 1/5., 25/(24*4.) );
135 add( 3/5., 1/5., 25/(24*4.) );
136 add( 1/5., 3/5., 25/(24*4.) );
137 add( 1/3., 1/3., -27/(24*4.) ); break;
140 const double P1 = 0.11169079483905, P2 = 0.0549758718227661;
141 const double a = 0.445948490915965, b = 0.091576213509771;
147 add( 1-2*a, a, P1 ); break;
150 const double A = 0.470142064105115;
151 const double B = 0.101286507323456;
152 const double P1 = 0.066197076394253;
153 const double P2 = 0.062969590272413;
154 add( 1/3., 1/3., 9/80. );
160 add( B, 1-2*B, P2 ); break;
163 const double A = 0.063089014491502;
164 const double B = 0.249286745170910;
165 const double C = 0.310352451033785;
166 const double D = 0.053145049844816;
167 const double P1 = 0.025422453185103;
168 const double P2 = 0.058393137863189;
169 const double P3 = 0.041425537809187;
181 add( D, 1-C-D, P3 ); break;
184 EXCEPTION( logic_error,"Invalid nb of gauss points for TRIA, variant 1: "
188 else if ( variant == 2 ) {
189 if (geom == eTRIA3) setRefCoords( TTria3a() );
190 else setRefCoords( TTria6a() );
193 add( -1/3., -1/3., 2. ); break;
196 add( -2/3., 1/3., 2/3. );
197 add( -2/3., -2/3., 2/3. );
198 add( 1/3., -2/3., 2/3. ); break;
201 const double P1 = 0.11169079483905, P2 = 0.0549758718227661;
202 const double A = 0.445948490915965, B = 0.091576213509771;
203 add( 2*B-1, 1-4*B, 4*P2 );
204 add( 2*B-1, 2*B-1, 4*P2 );
205 add( 1-4*B, 2*B-1, 4*P2 );
206 add( 1-4*A, 2*A-1, 4*P1 );
207 add( 2*A-1, 1-4*A, 4*P1 );
208 add( 2*A-1, 2*A-1, 4*P1 ); break;
211 EXCEPTION( logic_error,"Invalid nb of gauss points for TRIA, variant 2: "
215 else if ( variant == 3 ) {
216 if (geom == eTRIA3) setRefCoords( TTria3b() );
217 else setRefCoords( TTria6b() );
220 add( 1/3., 1/3., -27/96 );
221 add( 0.2 , 0.2 , 25/96 );
222 add( 0.6 , 0.2 , 25/96 );
223 add( 0.2 , 0.6 , 25/96 ); break;
226 EXCEPTION( logic_error,"Invalid nb of gauss points for TRIA, variant 3: "
234 if ( variant == 1 ) {
235 if (geom == eQUAD4) setRefCoords( TQuad4b() );
236 else setRefCoords( TQuad8b() );
239 add( 0, 0, 4 ); break;
242 const double a = 1/sqrt(3.);
246 add( -a, a, 1 ); break;
249 const double a = 0.774596669241483;
250 add( -a, -a, 25/81. );
251 add( a, -a, 25/81. );
253 add( -a, a, 25/81. );
254 add( 0., -a, 40/81. );
255 add( a, 0., 40/81. );
256 add( 0., a, 40/81. );
257 add( -a, 0., 40/81. );
258 add( 0., 0., 64/81. ); break;
261 EXCEPTION( logic_error,"Invalid nb of gauss points for QUAD, variant 1: "
265 else if ( variant == 2 ) {
266 if (geom == eQUAD4) setRefCoords( TQuad4a() );
267 else setRefCoords( TQuad8a() );
270 const double a = 1/sqrt(3.);
274 add( a, a, 1 ); break;
277 const double a = 0.774596669241483;
278 add( -a, a, 25/81. );
279 add( -a, -a, 25/81. );
280 add( a, -a, 25/81. );
282 add( -a, 0., 40/81. );
283 add( 0., -a, 40/81. );
284 add( a, 0., 40/81. );
285 add( 0., a, 40/81. );
286 add( 0., 0., 64/81. ); break;
289 EXCEPTION( logic_error,"Invalid nb of gauss points for QUAD, variant 1: "
293 else if ( variant == 3 ) {
294 if (geom == eQUAD4) setRefCoords( TQuad4b() );
295 else setRefCoords( TQuad8b() );
298 const double a = 3/sqrt(3.);
302 add( a, a, 1 ); break;
305 const double a = sqrt(3/5.), c1 = 5/9., c2 = 8/9.;
306 const double c12 = c1*c2, c22 = c2*c2, c1c2 = c1*c2;
315 add( a, a, c12 ); break;
318 EXCEPTION( logic_error,"Invalid nb of gauss points for QUAD, variant 3: "
326 if (geom == eTETRA4) setRefCoords( TTetra4a() );
327 else setRefCoords( TTetra10a() );
330 const double a = (5 - sqrt(5.))/20., b = (5 + 3*sqrt(5.))/20.;
331 add( a, a, a, 1/24. );
332 add( a, a, b, 1/24. );
333 add( a, b, a, 1/24. );
334 add( b, a, a, 1/24. ); break;
337 const double a = 0.25, b = 1/6., c = 0.5;
338 add( a, a, a, -2/15. );
339 add( b, b, b, 3/40. );
340 add( b, b, c, 3/40. );
341 add( b, c, b, 3/40. );
342 add( c, b, b, 3/40. ); break;
345 const double a = 0.25;
346 const double b1 = (7 + sqrt(15.))/34., c1 = (13 + 3*sqrt(15.))/34., d = (5 - sqrt(15.))/20.;
347 const double b2 = (7 - sqrt(15.))/34., c2 = (13 - 3*sqrt(15.))/34., e = (5 + sqrt(15.))/20.;
348 const double P1 = (2665 - 14*sqrt(15.))/226800.;
349 const double P2 = (2665 + 14*sqrt(15.))/226800.;
350 add( a, a, a, 8/405.);//_____
351 add( b1, b1, b1, P1 );
352 add( b1, b1, c1, P1 );
353 add( b1, c1, b1, P1 );
354 add( c1, b1, b1, P1 );//_____
355 add( b2, b2, b2, P2 );
356 add( b2, b2, c2, P2 );
357 add( b2, c2, b2, P2 );
358 add( c2, b2, b2, P2 );//_____
359 add( d, d, e, 5/567.);
360 add( d, e, d, 5/567.);
361 add( e, d, d, 5/567.);
362 add( d, e, e, 5/567.);
363 add( e, d, e, 5/567.);
364 add( e, e, d, 5/567.);
368 EXCEPTION( logic_error,"Invalid nb of gauss points for TETRA: "<<nbGauss);
374 if (geom == ePYRA5) setRefCoords( TPyra5a() );
375 else setRefCoords( TPyra13a() );
378 const double h1 = 0.1531754163448146;
379 const double h2 = 0.6372983346207416;
380 add( .5, 0., h1, 2/15. );
381 add( 0., .5, h1, 2/15. );
382 add( -.5, 0., h1, 2/15. );
383 add( 0., -.5, h1, 2/15. );
384 add( 0., 0., h2, 2/15. ); break;
387 const double p1 = 0.1024890634400000 ;
388 const double p2 = 0.1100000000000000 ;
389 const double p3 = 0.1467104129066667 ;
390 const double a = 0.5702963741068025 ;
391 const double h1 = 0.1666666666666666 ;
392 const double h2 = 0.08063183038464675;
393 const double h3 = 0.6098484849057127 ;
394 add( a, 0., h1, p1 );
395 add( 0., a, h1, p1 );
396 add( -a, 0., h1, p1 );
397 add( 0., -a, h1, p1 );
398 add( 0., 0., h2, p2 );
399 add( 0., 0., h3, p3 ); break;
402 const double a1 = 0.788073483;
403 const double b6 = 0.499369002;
404 const double b1 = 0.848418011;
405 const double c8 = 0.478508449;
406 const double c1 = 0.652816472;
407 const double d12 = 0.032303742;
408 const double d1 = 1.106412899;
409 double z = 1/2., fz = b1/2*(1 - z);
410 add( 0., 0., z, a1 ); // 1
411 add( fz, fz, z, b6 ); // 2
412 add( -fz, fz, z, b6 ); // 3
413 add( -fz, -fz, z, b6 ); // 4
414 add( fz, -fz, z, b6 ); // 5
416 add( 0., 0., z, b6 ); // 6
418 add( 0., 0., z, b6 ); // 7
419 z = (1 - c1)/2.; fz = c1*(1 - z);
420 add( fz, 0., z, c8 ); // 8
421 add( 0., fz, z, c8 ); // 9
422 add( -fz, 0., z, c8 ); // 10
423 add( 0., -fz, z, c8 ); // 11
424 z = (1 + c1)/2.; fz = c1*(1 - z);
425 add( fz, 0., z, c8 ); // 12
426 add( 0., fz, z, c8 ); // 13
427 add( -fz, 0., z, c8 ); // 14
428 add( 0., -fz, z, c8 ); // 15
429 z = (1 - d1)/2., fz = d1/2*(1 - z);
430 add( fz, fz, z, d12); // 16
431 add( -fz, fz, z, d12); // 17
432 add( -fz, -fz, z, d12); // 18
433 add( fz, -fz, z, d12); // 19
434 z = 1/2.; fz = d1*(1 - z);
435 add( fz, 0., z, d12); // 20
436 add( 0., fz, z, d12); // 21
437 add( -fz, 0., z, d12); // 22
438 add( 0., -fz, z, d12); // 23
439 z = (1 + d1)/2., fz = d1/2*(1 - z);
440 add( fz, fz, z, d12); // 24
441 add( -fz, fz, z, d12); // 25
442 add( -fz, -fz, z, d12); // 26
443 add( fz, -fz, z, d12); // 27
447 EXCEPTION( logic_error,"Invalid nb of gauss points for PYRA: "<<nbGauss);
452 if (geom == ePENTA6) setRefCoords( TPenta6a() );
453 else setRefCoords( TPenta15a() );
456 const double a = sqrt(3.)/3.;
457 add( -a, .5, .5, 1/6. );
458 add( -a, 0., .5, 1/6. );
459 add( -a, .5, 0., 1/6. );
460 add( a, .5, .5, 1/6. );
461 add( a, 0., .5, 1/6. );
462 add( a, .5, 0., 1/6. ); break;
465 const double a = 0.577350269189626;
466 add( -a, 1/3., 1/3., -27/96. );
467 add( -a, 0.6, 0.2, 25/96. );
468 add( -a, 0.2, 0.6, 25/96. );
469 add( -a, 0.2, 0.2, 25/96. );
470 add( +a, 1/3., 1/3., -27/96. );
471 add( +a, 0.6, 0.2, 25/96. );
472 add( +a, 0.2, 0.6, 25/96. );
473 add( +a, 0.2, 0.2, 25/96. ); break;
476 const double d = sqrt(3/5.), c1 = 5/9., c2 = 8/9.; // d <=> alfa
477 const double a = (6 + sqrt(15.))/21.;
478 const double b = (6 - sqrt(15.))/21.;
479 const double P1 = (155 + sqrt(15.))/2400.;
480 const double P2 = (155 - sqrt(15.))/2400.; //___
481 add( -d, 1/3., 1/3., c1*9/80. );//___
482 add( -d, a, a, c1*P1 );
483 add( -d, 1-2*a, a, c1*P1 );
484 add( -d, a, 1-2*a, c1*P1 );//___
485 add( -d, b, b, c1*P2 );
486 add( -d, 1-2*b, b, c1*P2 );
487 add( -d, b, 1-2*b, c1*P2 );//___
488 add( 0., 1/3., 1/3., c2*9/80. );//___
489 add( 0., a, a, c2*P1 );
490 add( 0., 1-2*a, a, c2*P1 );
491 add( 0., a, 1-2*a, c2*P1 );//___
492 add( 0., b, b, c2*P2 );
493 add( 0., 1-2*b, b, c2*P2 );
494 add( 0., b, 1-2*b, c2*P2 );//___
495 add( d, 1/3., 1/3., c1*9/80. );//___
496 add( d, a, a, c1*P1 );
497 add( d, 1-2*a, a, c1*P1 );
498 add( d, a, 1-2*a, c1*P1 );//___
499 add( d, b, b, c1*P2 );
500 add( d, 1-2*b, b, c1*P2 );
501 add( d, b, 1-2*b, c1*P2 );//___
505 EXCEPTION( logic_error,"Invalid nb of gauss points for PENTA: " <<nbGauss);
511 if (geom == eHEXA8) setRefCoords( THexa8a() );
512 else setRefCoords( THexa20a() );
515 const double a = sqrt(3.)/3.;
516 add( -a, -a, -a, 1. );
517 add( -a, -a, a, 1. );
518 add( -a, a, -a, 1. );
520 add( a, -a, -a, 1. );
523 add( a, a, a, 1. ); break;
526 const double a = sqrt(3/5.), c1 = 5/9., c2 = 8/9.;
527 const double c12 = c1*c1, c13 = c1*c1*c1;
528 const double c22 = c2*c2, c23 = c2*c2*c2;
529 add( -a, -a, -a, c13 ); // 1
530 add( -a, -a, 0., c12*c2 ); // 2
531 add( -a, -a, a, c13 ); // 3
532 add( -a, 0., -a, c12*c2 ); // 4
533 add( -a, 0., 0., c1*c22 ); // 5
534 add( -a, 0., a, c12*c2 ); // 6
535 add( -a, a, -a, c13 ); // 7
536 add( -a, a, 0., c12*c2 ); // 8
537 add( -a, a, a, c13 ); // 9
538 add( 0., -a, -a, c12*c2 ); // 10
539 add( 0., -a, 0., c1*c22 ); // 11
540 add( 0., -a, a, c12*c2 ); // 12
541 add( 0., 0., -a, c1*c22 ); // 13
542 add( 0., 0., 0., c23 ); // 14
543 add( 0., 0., a, c1*c22 ); // 15
544 add( 0., a, -a, c12*c2 ); // 16
545 add( 0., a, 0., c1*c22 ); // 17
546 add( 0., a, a, c12*c2 ); // 18
547 add( a, -a, -a, c13 ); // 19
548 add( a, -a, 0., c12*c2 ); // 20
549 add( a, -a, a, c13 ); // 21
550 add( a, 0., -a, c12*c2 ); // 22
551 add( a, 0., 0., c1*c22 ); // 23
552 add( a, 0., a, c12*c2 ); // 24
553 add( a, a, -a, c13 ); // 25
554 add( a, a, 0., c12*c2 ); // 26
555 add( a, a, a, c13 ); // 27
559 EXCEPTION( logic_error,"Invalid nb of gauss points for PENTA: " <<nbGauss);
564 EXCEPTION( logic_error,"unexpected EGeometrieElement: "<< geom);
567 if ( myWeights.capacity() != myWeights.size() )
568 EXCEPTION( logic_error,"Not all gauss points defined");