]> SALOME platform Git repositories - modules/smesh.git/blob - src/MEDWrapper/Base/MED_GaussDef.cxx
Salome HOME
CMake:
[modules/smesh.git] / src / MEDWrapper / Base / MED_GaussDef.cxx
1 // Copyright (C) 2007-2013  CEA/DEN, EDF R&D, OPEN CASCADE
2 //
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.
7 //
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.
12 //
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
16 //
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
18 //
19 // File   : MED_GaussDef.hxx
20 // Author : Edward AGAPOV (eap)
21 //
22 #include "MED_GaussDef.hxx"
23 #include "MED_Utilities.hxx"
24 #include "MED_GaussUtils.hxx"
25
26 namespace MED
27 {
28   using namespace std;
29   using namespace MED;
30   //---------------------------------------------------------------
31
32   void TGaussDef::add(const double x, const double weight)
33   {
34     if ( dim() != 1 )
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 );
40   }
41   void TGaussDef::add(const double x, const double y, const double weight)
42   {
43     if ( dim() != 2 )
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 );
50   }
51   void TGaussDef::add(const double x, const double y, const double z, const double weight)
52   {
53     if ( dim() != 3 )
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 );
61   }
62   void TGaussDef::setRefCoords(const TShapeFun& aShapeFun)
63   {
64     myRefCoords.reserve( aShapeFun.myRefCoord.size() );
65     myRefCoords.assign( aShapeFun.myRefCoord.begin(),
66                         aShapeFun.myRefCoord.end() );
67   }
68
69
70   //---------------------------------------------------------------
71   /*!
72    * \brief Fill definition of gauss points family
73    */
74   //---------------------------------------------------------------
75
76   TGaussDef::TGaussDef(const int geom, const int nbGauss, const int variant)
77   {
78     myType = geom;
79     myCoords .reserve( nbGauss * dim() );
80     myWeights.reserve( nbGauss );
81
82     switch ( geom ) {
83
84     case eSEG2:
85     case eSEG3:
86       if (geom == eSEG2) setRefCoords( TSeg2a() );
87       else               setRefCoords( TSeg3a() );
88       switch ( nbGauss ) {
89       case 1: {
90         add( 0.0, 2.0 ); break;
91       }
92       case 2: {
93         const double a = 0.577350269189626;
94         add(  a,  1.0 );
95         add( -a,  1.0 ); break;
96       }
97       case 3: {
98         const double a = 0.774596669241;
99         const double P1 = 1./1.8;
100         const double P2 = 1./1.125;
101         add( -a,  P1 );
102         add(  0,  P2 ); 
103         add(  a,  P1 ); break;
104       }
105       case 4: {
106         const double a  = 0.339981043584856, b  = 0.861136311594053;
107         const double P1 = 0.652145154862546, P2 = 0.347854845137454 ;
108         add(  a,  P1 );
109         add( -a,  P1 );
110         add(  b,  P2 ); 
111         add( -b,  P2 ); break;
112       }
113       default:
114         EXCEPTION( logic_error,"Invalid nb of gauss points for SEG"<<nbGauss);
115       }
116       break;
117
118     case eTRIA3:
119     case eTRIA6:
120       if ( variant == 1 ) {
121         if (geom == eTRIA3) setRefCoords( TTria3b() );
122         else                setRefCoords( TTria6b() );
123         switch ( nbGauss ) {
124         case 1: { // FPG1
125           add( 1/3., 1/3., 1/2. ); break;
126         }
127         case 3: { // FPG3
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;
132         }
133         case 4: { // FPG4
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;
138         }
139         case 6: { // FPG6
140           const double P1 = 0.11169079483905, P2 = 0.0549758718227661;
141           const double a  = 0.445948490915965, b = 0.091576213509771;
142           add(     b,     b, P2 ); 
143           add( 1-2*b,     b, P2 );
144           add(     b, 1-2*b, P2 );
145           add(     a, 1-2*a, P1 );
146           add(     a,     a, P1 ); 
147           add( 1-2*a,     a, P1 ); break;
148         }
149         case 7: { // FPG7
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. ); 
155           add(     A,     A, P1 ); 
156           add( 1-2*A,     A, P1 );
157           add(     A, 1-2*A, P1 );
158           add(     B,     B, P2 ); 
159           add( 1-2*B,     B, P2 );
160           add(     B, 1-2*B, P2 ); break;
161         }
162         case 12: { // FPG12
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;
170           add(     A,     A, P1 ); 
171           add( 1-2*A,     A, P1 );
172           add(     A, 1-2*A, P1 );
173           add(     B,     B, P2 ); 
174           add( 1-2*B,     B, P2 );
175           add(     B, 1-2*B, P2 );
176           add(     C,     D, P3 );
177           add(     D,     C, P3 );
178           add( 1-C-D,     C, P3 );
179           add( 1-C-D,     D, P3 );
180           add(     C, 1-C-D, P3 );
181           add(     D, 1-C-D, P3 ); break;
182         }
183         default:
184           EXCEPTION( logic_error,"Invalid nb of gauss points for TRIA, variant 1: "
185                      <<nbGauss);
186         }
187       }
188       else if ( variant == 2 ) {
189         if (geom == eTRIA3) setRefCoords( TTria3a() );
190         else                setRefCoords( TTria6a() );
191         switch ( nbGauss ) {
192         case 1: {
193           add( -1/3., -1/3., 2. ); break;
194         }
195         case 3: {
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;
199         }
200         case 6: {
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;
209         }
210         default:
211           EXCEPTION( logic_error,"Invalid nb of gauss points for TRIA, variant 2: "
212                      <<nbGauss);
213         }
214       }
215       else if ( variant == 3 ) {
216         if (geom == eTRIA3) setRefCoords( TTria3b() );
217         else                setRefCoords( TTria6b() );
218         switch ( nbGauss ) {
219         case 4: {
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;
224         }
225         default:
226           EXCEPTION( logic_error,"Invalid nb of gauss points for TRIA, variant 3: "
227                      <<nbGauss);
228         }
229       }
230       break;
231
232     case eQUAD4:
233     case eQUAD8:
234       if ( variant == 1 ) {
235         if (geom == eQUAD4) setRefCoords( TQuad4b() );
236         else                setRefCoords( TQuad8b() );
237         switch ( nbGauss ) {
238         case 1: { // FPG1
239           add(  0,  0,  4 ); break;
240         }
241         case 4: { // FPG4
242           const double a = 1/sqrt(3.);
243           add( -a, -a,  1 );
244           add(  a, -a,  1 );
245           add(  a,  a,  1 );
246           add( -a,  a,  1 ); break;
247         }
248         case 9: { // FPG9
249           const double a = 0.774596669241483;
250           add( -a, -a,  25/81. );
251           add(  a, -a,  25/81. );
252           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;
259         }
260         default:
261           EXCEPTION( logic_error,"Invalid nb of gauss points for QUAD, variant 1: "
262                      <<nbGauss);
263         }
264       }
265       else if ( variant == 2 ) {
266         if (geom == eQUAD4) setRefCoords( TQuad4a() );
267         else                setRefCoords( TQuad8a() );
268         switch ( nbGauss ) {
269         case 4: {
270           const double a = 1/sqrt(3.);
271           add( -a,  a,  1 );
272           add( -a, -a,  1 );
273           add(  a, -a,  1 );
274           add(  a,  a,  1 ); break;
275         }
276         case 9: {
277           const double a = 0.774596669241483;
278           add( -a,  a,  25/81. );
279           add( -a, -a,  25/81. );
280           add(  a, -a,  25/81. );
281           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;
287         }
288         default:
289           EXCEPTION( logic_error,"Invalid nb of gauss points for QUAD, variant 1: "
290                      <<nbGauss);
291         }
292       }
293       else if ( variant == 3 ) {
294         if (geom == eQUAD4) setRefCoords( TQuad4b() );
295         else                setRefCoords( TQuad8b() );
296         switch ( nbGauss ) {
297         case 4: {
298           const double a = 3/sqrt(3.);
299           add( -a, -a,  1 );
300           add( -a,  a,  1 );
301           add(  a, -a,  1 );
302           add(  a,  a,  1 ); break;
303         }
304         case 9: {
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;
307           add( -a, -a,  c12  );
308           add( -a, 0.,  c1c2 );
309           add( -a,  a,  c12  );
310           add( 0., -a,  c1c2 );
311           add( 0., 0.,  c22  );
312           add( 0.,  a,  c1c2 );
313           add(  a, -a,  c12  );
314           add(  a, 0.,  c1c2 );
315           add(  a,  a,  c12  ); break;
316         }
317         default:
318           EXCEPTION( logic_error,"Invalid nb of gauss points for QUAD, variant 3: "
319                      <<nbGauss);
320         }
321       }
322       break;
323
324     case eTETRA4:
325     case eTETRA10:
326       if (geom == eTETRA4) setRefCoords( TTetra4a() );
327       else                 setRefCoords( TTetra10a() );
328       switch ( nbGauss ) {
329       case 4: { // FPG4
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;
335       }
336       case 5: { // FPG5
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;
343       }
344       case 15: { // FPG15
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.);
365         break;
366       }
367       default:
368         EXCEPTION( logic_error,"Invalid nb of gauss points for TETRA: "<<nbGauss);
369       }
370       break;
371
372     case ePYRA5:
373     case ePYRA13:
374       if (geom == ePYRA5) setRefCoords( TPyra5a() );
375       else                setRefCoords( TPyra13a() );
376       switch ( nbGauss ) {
377       case 5: { // FPG5
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;
385       }
386       case 6: { // FPG6
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;
400       }
401       case 27: { // FPG27
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
415         z = (1 - b1)/2.;
416         add(  0.,  0.,   z,  b6 ); // 6
417         z = (1 + b1)/2.;
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
444         break;
445       }
446       default:
447         EXCEPTION( logic_error,"Invalid nb of gauss points for PYRA: "<<nbGauss);
448       }
449       break;
450     case ePENTA6:
451     case ePENTA15:
452       if (geom == ePENTA6) setRefCoords( TPenta6a() );
453       else                 setRefCoords( TPenta15a() );
454       switch ( nbGauss ) {
455       case 6: { // FPG6
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;
463       }
464       case 8: { // FPG8
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;
474       }
475       case 21: { // FPG21
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    );//___
502         break;
503       }
504       default:
505         EXCEPTION( logic_error,"Invalid nb of gauss points for PENTA: " <<nbGauss);
506       }
507       break;
508
509     case eHEXA8:
510     case eHEXA20:
511       if (geom == eHEXA8) setRefCoords( THexa8a() );
512       else                setRefCoords( THexa20a() );
513       switch ( nbGauss ) {
514       case 8: { // FPG8
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. );
519         add( -a,  a,  a,  1. );
520         add(  a, -a, -a,  1. );
521         add(  a, -a,  a,  1. );
522         add(  a,  a, -a,  1. );
523         add(  a,  a,  a,  1. ); break;
524       }
525       case 27: { // FPG27
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
556         break;
557       }
558       default:
559         EXCEPTION( logic_error,"Invalid nb of gauss points for PENTA: " <<nbGauss);
560       }
561       break;
562
563     default:
564       EXCEPTION( logic_error,"unexpected EGeometrieElement: "<< geom);
565     }
566
567     if ( myWeights.capacity() != myWeights.size() )
568       EXCEPTION( logic_error,"Not all gauss points defined");
569   }
570 }