src/MEDWrapper/MED_GaussDef.cxx

   1 // Copyright (C) 2007-2016  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, or (at your option) any later version.
   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
  20 #include "MED_GaussDef.hxx"
  21 #include "MED_Utilities.hxx"
  22 #include "MED_GaussUtils.hxx"
  23
  24 namespace MED
  25 {
  26   using namespace MED;
  27   //---------------------------------------------------------------
  28
  29   void TGaussDef::add(const double x, const double weight)
  30   {
  31     if ( dim() != 1 )
  32       EXCEPTION( std::logic_error,"dim() != 1");
  33     if ( myWeights.capacity() == myWeights.size() )
  34       EXCEPTION( std::logic_error,"Extra gauss point");
  35     myCoords.push_back( x );
  36     myWeights.push_back( weight );
  37   }
  38   void TGaussDef::add(const double x, const double y, const double weight)
  39   {
  40     if ( dim() != 2 )
  41       EXCEPTION( std::logic_error,"dim() != 2");
  42     if ( myWeights.capacity() == myWeights.size() )
  43       EXCEPTION( std::logic_error,"Extra gauss point");
  44     myCoords.push_back( x );
  45     myCoords.push_back( y );
  46     myWeights.push_back( weight );
  47   }
  48   void TGaussDef::add(const double x, const double y, const double z, const double weight)
  49   {
  50     if ( dim() != 3 )
  51       EXCEPTION( std::logic_error,"dim() != 3");
  52     if ( myWeights.capacity() == myWeights.size() )
  53       EXCEPTION( std::logic_error,"Extra gauss point");
  54     myCoords.push_back( x );
  55     myCoords.push_back( y );
  56     myCoords.push_back( z );
  57     myWeights.push_back( weight );
  58   }
  59   void TGaussDef::setRefCoords(const TShapeFun& aShapeFun)
  60   {
  61     myRefCoords.reserve( aShapeFun.myRefCoord.size() );
  62     myRefCoords.assign( aShapeFun.myRefCoord.begin(),
  63                         aShapeFun.myRefCoord.end() );
  64   }
  65
  66
  67   //---------------------------------------------------------------
  68   /*!
  69    * \brief Fill definition of gauss points family
  70    */
  71   //---------------------------------------------------------------
  72
  73   TGaussDef::TGaussDef(const int geom, const int nbGauss, const int variant)
  74   {
  75     myType = geom;
  76     myCoords .reserve( nbGauss * dim() );
  77     myWeights.reserve( nbGauss );
  78
  79     switch ( geom ) {
  80
  81     case eSEG2:
  82     case eSEG3:
  83       if (geom == eSEG2) setRefCoords( TSeg2a() );
  84       else               setRefCoords( TSeg3a() );
  85       switch ( nbGauss ) {
  86       case 1: {
  87         add( 0.0, 2.0 ); break;
  88       }
  89       case 2: {
  90         const double a = 0.577350269189626;
  91         add(  a,  1.0 );
  92         add( -a,  1.0 ); break;
  93       }
  94       case 3: {
  95         const double a = 0.774596669241;
  96         const double P1 = 1./1.8;
  97         const double P2 = 1./1.125;
  98         add( -a,  P1 );
  99         add(  0,  P2 );
 100         add(  a,  P1 ); break;
 101       }
 102       case 4: {
 103         const double a  = 0.339981043584856, b  = 0.861136311594053;
 104         const double P1 = 0.652145154862546, P2 = 0.347854845137454 ;
 105         add(  a,  P1 );
 106         add( -a,  P1 );
 107         add(  b,  P2 );
 108         add( -b,  P2 ); break;
 109       }
 110       default:
 111         EXCEPTION( std::logic_error,"Invalid nb of gauss points for SEG"<<nbGauss);
 112       }
 113       break;
 114
 115     case eTRIA3:
 116     case eTRIA6:
 117       if ( variant == 1 ) {
 118         if (geom == eTRIA3) setRefCoords( TTria3b() );
 119         else                setRefCoords( TTria6b() );
 120         switch ( nbGauss ) {
 121         case 1: { // FPG1
 122           add( 1/3., 1/3., 1/2. ); break;
 123         }
 124         case 3: { // FPG3
 125           // what about COT3 ???
 126           add( 1/6., 1/6., 1/6. );
 127           add( 2/3., 1/6., 1/6. );
 128           add( 1/6., 2/3., 1/6. ); break;
 129         }
 130         case 4: { // FPG4
 131           add( 1/5., 1/5.,  25/(24*4.) );
 132           add( 3/5., 1/5.,  25/(24*4.) );
 133           add( 1/5., 3/5.,  25/(24*4.) );
 134           add( 1/3., 1/3., -27/(24*4.) ); break;
 135         }
 136         case 6: { // FPG6
 137           const double P1 = 0.11169079483905, P2 = 0.0549758718227661;
 138           const double a  = 0.445948490915965, b = 0.091576213509771;
 139           add(     b,     b, P2 );
 140           add( 1-2*b,     b, P2 );
 141           add(     b, 1-2*b, P2 );
 142           add(     a, 1-2*a, P1 );
 143           add(     a,     a, P1 );
 144           add( 1-2*a,     a, P1 ); break;
 145         }
 146         case 7: { // FPG7
 147           const double A  = 0.470142064105115;
 148           const double B  = 0.101286507323456;
 149           const double P1 = 0.066197076394253;
 150           const double P2 = 0.062969590272413;
 151           add(  1/3.,  1/3., 9/80. );
 152           add(     A,     A, P1 );
 153           add( 1-2*A,     A, P1 );
 154           add(     A, 1-2*A, P1 );
 155           add(     B,     B, P2 );
 156           add( 1-2*B,     B, P2 );
 157           add(     B, 1-2*B, P2 ); break;
 158         }
 159         case 12: { // FPG12
 160           const double A  = 0.063089014491502;
 161           const double B  = 0.249286745170910;
 162           const double C  = 0.310352451033785;
 163           const double D  = 0.053145049844816;
 164           const double P1 = 0.025422453185103;
 165           const double P2 = 0.058393137863189;
 166           const double P3 = 0.041425537809187;
 167           add(     A,     A, P1 );
 168           add( 1-2*A,     A, P1 );
 169           add(     A, 1-2*A, P1 );
 170           add(     B,     B, P2 );
 171           add( 1-2*B,     B, P2 );
 172           add(     B, 1-2*B, P2 );
 173           add(     C,     D, P3 );
 174           add(     D,     C, P3 );
 175           add( 1-C-D,     C, P3 );
 176           add( 1-C-D,     D, P3 );
 177           add(     C, 1-C-D, P3 );
 178           add(     D, 1-C-D, P3 ); break;
 179         }
 180         default:
 181           EXCEPTION( std::logic_error,"Invalid nb of gauss points for TRIA, variant 1: "
 182                      <<nbGauss);
 183         }
 184       }
 185       else if ( variant == 2 ) {
 186         if (geom == eTRIA3) setRefCoords( TTria3a() );
 187         else                setRefCoords( TTria6a() );
 188         switch ( nbGauss ) {
 189         case 1: {
 190           add( -1/3., -1/3., 2. ); break;
 191         }
 192         case 3: {
 193           add( -2/3.,  1/3., 2/3. );
 194           add( -2/3., -2/3., 2/3. );
 195           add(  1/3., -2/3., 2/3. ); break;
 196         }
 197         case 6: {
 198           const double P1 = 0.11169079483905, P2 = 0.0549758718227661;
 199           const double A  = 0.445948490915965, B = 0.091576213509771;
 200           add( 2*B-1, 1-4*B, 4*P2 );
 201           add( 2*B-1, 2*B-1, 4*P2 );
 202           add( 1-4*B, 2*B-1, 4*P2 );
 203           add( 1-4*A, 2*A-1, 4*P1 );
 204           add( 2*A-1, 1-4*A, 4*P1 );
 205           add( 2*A-1, 2*A-1, 4*P1 ); break;
 206         }
 207         default:
 208           EXCEPTION( std::logic_error,"Invalid nb of gauss points for TRIA, variant 2: "
 209                      <<nbGauss);
 210         }
 211       }
 212       else if ( variant == 3 ) {
 213         if (geom == eTRIA3) setRefCoords( TTria3b() );
 214         else                setRefCoords( TTria6b() );
 215         switch ( nbGauss ) {
 216         case 4: {
 217           add( 1/3., 1/3., -27/96 );
 218           add( 0.2 , 0.2 ,  25/96 );
 219           add( 0.6 , 0.2 ,  25/96 );
 220           add( 0.2 , 0.6 ,  25/96 ); break;
 221         }
 222         default:
 223           EXCEPTION( std::logic_error,"Invalid nb of gauss points for TRIA, variant 3: "
 224                      <<nbGauss);
 225         }
 226       }
 227       break;
 228
 229     case eQUAD4:
 230     case eQUAD8:
 231       if ( variant == 1 ) {
 232         if (geom == eQUAD4) setRefCoords( TQuad4b() );
 233         else                setRefCoords( TQuad8b() );
 234         switch ( nbGauss ) {
 235         case 1: { // FPG1
 236           add(  0,  0,  4 ); break;
 237         }
 238         case 4: { // FPG4
 239           const double a = 1/sqrt(3.);
 240           add( -a, -a,  1 );
 241           add(  a, -a,  1 );
 242           add(  a,  a,  1 );
 243           add( -a,  a,  1 ); break;
 244         }
 245         case 9: { // FPG9
 246           const double a = 0.774596669241483;
 247           add( -a, -a,  25/81. );
 248           add(  a, -a,  25/81. );
 249           add(  a,  a,  25/81. );
 250           add( -a,  a,  25/81. );
 251           add( 0., -a,  40/81. );
 252           add(  a, 0.,  40/81. );
 253           add( 0.,  a,  40/81. );
 254           add( -a, 0.,  40/81. );
 255           add( 0., 0.,  64/81. ); break;
 256         }
 257         default:
 258           EXCEPTION( std::logic_error,"Invalid nb of gauss points for QUAD, variant 1: "
 259                      <<nbGauss);
 260         }
 261       }
 262       else if ( variant == 2 ) {
 263         if (geom == eQUAD4) setRefCoords( TQuad4a() );
 264         else                setRefCoords( TQuad8a() );
 265         switch ( nbGauss ) {
 266         case 4: {
 267           const double a = 1/sqrt(3.);
 268           add( -a,  a,  1 );
 269           add( -a, -a,  1 );
 270           add(  a, -a,  1 );
 271           add(  a,  a,  1 ); break;
 272         }
 273         case 9: {
 274           const double a = 0.774596669241483;
 275           add( -a,  a,  25/81. );
 276           add( -a, -a,  25/81. );
 277           add(  a, -a,  25/81. );
 278           add(  a,  a,  25/81. );
 279           add( -a, 0.,  40/81. );
 280           add( 0., -a,  40/81. );
 281           add(  a, 0.,  40/81. );
 282           add( 0.,  a,  40/81. );
 283           add( 0., 0.,  64/81. ); break;
 284         }
 285         default:
 286           EXCEPTION( std::logic_error,"Invalid nb of gauss points for QUAD, variant 1: "
 287                      <<nbGauss);
 288         }
 289       }
 290       else if ( variant == 3 ) {
 291         if (geom == eQUAD4) setRefCoords( TQuad4b() );
 292         else                setRefCoords( TQuad8b() );
 293         switch ( nbGauss ) {
 294         case 4: {
 295           const double a = 3/sqrt(3.);
 296           add( -a, -a,  1 );
 297           add( -a,  a,  1 );
 298           add(  a, -a,  1 );
 299           add(  a,  a,  1 ); break;
 300         }
 301         case 9: {
 302           const double a = sqrt(3/5.), c1 = 5/9., c2 = 8/9.;
 303           const double c12 = c1*c2, c22 = c2*c2, c1c2 = c1*c2;
 304           add( -a, -a,  c12  );
 305           add( -a, 0.,  c1c2 );
 306           add( -a,  a,  c12  );
 307           add( 0., -a,  c1c2 );
 308           add( 0., 0.,  c22  );
 309           add( 0.,  a,  c1c2 );
 310           add(  a, -a,  c12  );
 311           add(  a, 0.,  c1c2 );
 312           add(  a,  a,  c12  ); break;
 313         }
 314         default:
 315           EXCEPTION( std::logic_error,"Invalid nb of gauss points for QUAD, variant 3: "
 316                      <<nbGauss);
 317         }
 318       }
 319       break;
 320
 321     case eTETRA4:
 322     case eTETRA10:
 323       if (geom == eTETRA4) setRefCoords( TTetra4a() );
 324       else                 setRefCoords( TTetra10a() );
 325       switch ( nbGauss ) {
 326       case 4: { // FPG4
 327         const double a = (5 - sqrt(5.))/20., b = (5 + 3*sqrt(5.))/20.;
 328         add(  a,  a,  a,  1/24. );
 329         add(  a,  a,  b,  1/24. );
 330         add(  a,  b,  a,  1/24. );
 331         add(  b,  a,  a,  1/24. ); break;
 332       }
 333       case 5: { // FPG5
 334         const double a = 0.25, b = 1/6., c = 0.5;
 335         add(  a,  a,  a, -2/15. );
 336         add(  b,  b,  b,  3/40. );
 337         add(  b,  b,  c,  3/40. );
 338         add(  b,  c,  b,  3/40. );
 339         add(  c,  b,  b,  3/40. ); break;
 340       }
 341       case 15: { // FPG15
 342         const double a = 0.25;
 343         const double b1 = (7 + sqrt(15.))/34., c1 = (13 + 3*sqrt(15.))/34., d = (5 - sqrt(15.))/20.;
 344         const double b2 = (7 - sqrt(15.))/34., c2 = (13 - 3*sqrt(15.))/34., e = (5 + sqrt(15.))/20.;
 345         const double P1 = (2665 - 14*sqrt(15.))/226800.;
 346         const double P2 = (2665 + 14*sqrt(15.))/226800.;
 347         add(  a,  a,  a,  8/405.);//_____
 348         add( b1, b1, b1,  P1    );
 349         add( b1, b1, c1,  P1    );
 350         add( b1, c1, b1,  P1    );
 351         add( c1, b1, b1,  P1    );//_____
 352         add( b2, b2, b2,  P2    );
 353         add( b2, b2, c2,  P2    );
 354         add( b2, c2, b2,  P2    );
 355         add( c2, b2, b2,  P2    );//_____
 356         add(  d,  d,  e,  5/567.);
 357         add(  d,  e,  d,  5/567.);
 358         add(  e,  d,  d,  5/567.);
 359         add(  d,  e,  e,  5/567.);
 360         add(  e,  d,  e,  5/567.);
 361         add(  e,  e,  d,  5/567.);
 362         break;
 363       }
 364       default:
 365         EXCEPTION( std::logic_error,"Invalid nb of gauss points for TETRA: "<<nbGauss);
 366       }
 367       break;
 368
 369     case ePYRA5:
 370     case ePYRA13:
 371       if (geom == ePYRA5) setRefCoords( TPyra5a() );
 372       else                setRefCoords( TPyra13a() );
 373       switch ( nbGauss ) {
 374       case 5: { // FPG5
 375         const double h1 = 0.1531754163448146;
 376         const double h2 = 0.6372983346207416;
 377         add(  .5,  0.,  h1,  2/15. );
 378         add(  0.,  .5,  h1,  2/15. );
 379         add( -.5,  0.,  h1,  2/15. );
 380         add(  0., -.5,  h1,  2/15. );
 381         add(  0.,  0.,  h2,  2/15. ); break;
 382       }
 383       case 6: { // FPG6
 384         const double p1 = 0.1024890634400000 ;
 385         const double p2 = 0.1100000000000000 ;
 386         const double p3 = 0.1467104129066667 ;
 387         const double a  = 0.5702963741068025 ;
 388         const double h1 = 0.1666666666666666 ;
 389         const double h2 = 0.08063183038464675;
 390         const double h3 = 0.6098484849057127 ;
 391         add(  a, 0.,  h1,  p1 );
 392         add( 0.,  a,  h1,  p1 );
 393         add( -a, 0.,  h1,  p1 );
 394         add( 0., -a,  h1,  p1 );
 395         add( 0., 0.,  h2,  p2 );
 396         add( 0., 0.,  h3,  p3 ); break;
 397       }
 398       case 27: { // FPG27
 399         const double a1  = 0.788073483;
 400         const double b6  = 0.499369002;
 401         const double b1  = 0.848418011;
 402         const double c8  = 0.478508449;
 403         const double c1  = 0.652816472;
 404         const double d12 = 0.032303742;
 405         const double d1  = 1.106412899;
 406         double z = 1/2., fz = b1/2*(1 - z);
 407         add(  0.,  0.,   z,  a1 ); // 1
 408         add(  fz,  fz,   z,  b6 ); // 2
 409         add( -fz,  fz,   z,  b6 ); // 3
 410         add( -fz, -fz,   z,  b6 ); // 4
 411         add(  fz, -fz,   z,  b6 ); // 5
 412         z = (1 - b1)/2.;
 413         add(  0.,  0.,   z,  b6 ); // 6
 414         z = (1 + b1)/2.;
 415         add(  0.,  0.,   z,  b6 ); // 7
 416         z = (1 - c1)/2.; fz = c1*(1 - z);
 417         add(  fz,  0.,   z,  c8 ); // 8
 418         add(  0.,  fz,   z,  c8 ); // 9
 419         add( -fz,  0.,   z,  c8 ); // 10
 420         add(  0., -fz,   z,  c8 ); // 11
 421         z = (1 + c1)/2.; fz = c1*(1 - z);
 422         add(  fz,  0.,   z,  c8 ); // 12
 423         add(  0.,  fz,   z,  c8 ); // 13
 424         add( -fz,  0.,   z,  c8 ); // 14
 425         add(  0., -fz,   z,  c8 ); // 15
 426         z = (1 - d1)/2., fz = d1/2*(1 - z);
 427         add(  fz,  fz,   z,  d12); // 16
 428         add( -fz,  fz,   z,  d12); // 17
 429         add( -fz, -fz,   z,  d12); // 18
 430         add(  fz, -fz,   z,  d12); // 19
 431         z = 1/2.; fz = d1*(1 - z);
 432         add(  fz,  0.,   z,  d12); // 20
 433         add(  0.,  fz,   z,  d12); // 21
 434         add( -fz,  0.,   z,  d12); // 22
 435         add(  0., -fz,   z,  d12); // 23
 436         z = (1 + d1)/2., fz = d1/2*(1 - z);
 437         add(  fz,  fz,   z,  d12); // 24
 438         add( -fz,  fz,   z,  d12); // 25
 439         add( -fz, -fz,   z,  d12); // 26
 440         add(  fz, -fz,   z,  d12); // 27
 441         break;
 442       }
 443       default:
 444         EXCEPTION( std::logic_error,"Invalid nb of gauss points for PYRA: "<<nbGauss);
 445       }
 446       break;
 447     case ePENTA6:
 448     case ePENTA15:
 449       if (geom == ePENTA6) setRefCoords( TPenta6a() );
 450       else                 setRefCoords( TPenta15a() );
 451       switch ( nbGauss ) {
 452       case 6: { // FPG6
 453         const double a = sqrt(3.)/3.;
 454         add( -a, .5, .5,  1/6. );
 455         add( -a, 0., .5,  1/6. );
 456         add( -a, .5, 0.,  1/6. );
 457         add(  a, .5, .5,  1/6. );
 458         add(  a, 0., .5,  1/6. );
 459         add(  a, .5, 0.,  1/6. ); break;
 460       }
 461       case 8: { // FPG8
 462         const double a = 0.577350269189626;
 463         add( -a, 1/3., 1/3., -27/96. );
 464         add( -a,  0.6,  0.2,  25/96. );
 465         add( -a,  0.2,  0.6,  25/96. );
 466         add( -a,  0.2,  0.2,  25/96. );
 467         add( +a, 1/3., 1/3., -27/96. );
 468         add( +a,  0.6,  0.2,  25/96. );
 469         add( +a,  0.2,  0.6,  25/96. );
 470         add( +a,  0.2,  0.2,  25/96. ); break;
 471       }
 472       case 21: { // FPG21
 473         const double d = sqrt(3/5.), c1 = 5/9., c2 = 8/9.; // d <=> alfa
 474         const double a = (6 + sqrt(15.))/21.;
 475         const double b = (6 - sqrt(15.))/21.;
 476         const double P1 = (155 + sqrt(15.))/2400.;
 477         const double P2 = (155 - sqrt(15.))/2400.;  //___
 478         add( -d,  1/3.,  1/3., c1*9/80. );//___
 479         add( -d,     a,     a, c1*P1    );
 480         add( -d, 1-2*a,     a, c1*P1    );
 481         add( -d,     a, 1-2*a, c1*P1    );//___
 482         add( -d,     b,     b, c1*P2    );
 483         add( -d, 1-2*b,     b, c1*P2    );
 484         add( -d,     b, 1-2*b, c1*P2    );//___
 485         add( 0.,  1/3.,  1/3., c2*9/80. );//___
 486         add( 0.,     a,     a, c2*P1    );
 487         add( 0., 1-2*a,     a, c2*P1    );
 488         add( 0.,     a, 1-2*a, c2*P1    );//___
 489         add( 0.,     b,     b, c2*P2    );
 490         add( 0., 1-2*b,     b, c2*P2    );
 491         add( 0.,     b, 1-2*b, c2*P2    );//___
 492         add(  d,  1/3.,  1/3., c1*9/80. );//___
 493         add(  d,     a,     a, c1*P1    );
 494         add(  d, 1-2*a,     a, c1*P1    );
 495         add(  d,     a, 1-2*a, c1*P1    );//___
 496         add(  d,     b,     b, c1*P2    );
 497         add(  d, 1-2*b,     b, c1*P2    );
 498         add(  d,     b, 1-2*b, c1*P2    );//___
 499         break;
 500       }
 501       default:
 502         EXCEPTION( std::logic_error,"Invalid nb of gauss points for PENTA: " <<nbGauss);
 503       }
 504       break;
 505
 506     case eHEXA8:
 507     case eHEXA20:
 508       if (geom == eHEXA8) setRefCoords( THexa8a() );
 509       else                setRefCoords( THexa20a() );
 510       switch ( nbGauss ) {
 511       case 8: { // FPG8
 512         const double a = sqrt(3.)/3.;
 513         add( -a, -a, -a,  1. );
 514         add( -a, -a,  a,  1. );
 515         add( -a,  a, -a,  1. );
 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. ); break;
 521       }
 522       case 27: { // FPG27
 523         const double a = sqrt(3/5.), c1 = 5/9., c2 = 8/9.;
 524         const double c12 = c1*c1, c13 = c1*c1*c1;
 525         const double c22 = c2*c2, c23 = c2*c2*c2;
 526         add( -a, -a, -a,   c13  ); // 1
 527         add( -a, -a, 0., c12*c2 ); // 2
 528         add( -a, -a,  a,   c13  ); // 3
 529         add( -a, 0., -a, c12*c2 ); // 4
 530         add( -a, 0., 0., c1*c22 ); // 5
 531         add( -a, 0.,  a, c12*c2 ); // 6
 532         add( -a,  a, -a,   c13  ); // 7
 533         add( -a,  a, 0., c12*c2 ); // 8
 534         add( -a,  a,  a,   c13  ); // 9
 535         add( 0., -a, -a, c12*c2 ); // 10
 536         add( 0., -a, 0., c1*c22 ); // 11
 537         add( 0., -a,  a, c12*c2 ); // 12
 538         add( 0., 0., -a, c1*c22 ); // 13
 539         add( 0., 0., 0.,   c23  ); // 14
 540         add( 0., 0.,  a, c1*c22 ); // 15
 541         add( 0.,  a, -a, c12*c2 ); // 16
 542         add( 0.,  a, 0., c1*c22 ); // 17
 543         add( 0.,  a,  a, c12*c2 ); // 18
 544         add(  a, -a, -a,   c13  ); // 19
 545         add(  a, -a, 0., c12*c2 ); // 20
 546         add(  a, -a,  a,   c13  ); // 21
 547         add(  a, 0., -a, c12*c2 ); // 22
 548         add(  a, 0., 0., c1*c22 ); // 23
 549         add(  a, 0.,  a, c12*c2 ); // 24
 550         add(  a,  a, -a,   c13  ); // 25
 551         add(  a,  a, 0., c12*c2 ); // 26
 552         add(  a,  a,  a,   c13  ); // 27
 553         break;
 554       }
 555       default:
 556         EXCEPTION( std::logic_error,"Invalid nb of gauss points for PENTA: " <<nbGauss);
 557       }
 558       break;
 559
 560     default:
 561       EXCEPTION( std::logic_error,"unexpected EGeometrieElement: "<< geom);
 562     }
 563
 564     if ( myWeights.capacity() != myWeights.size() )
 565       EXCEPTION( std::logic_error,"Not all gauss points defined");
 566   }
 567 }