src/MEDWrapper/MED_GaussDef.cxx

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