src/MEDWrapper/Base/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 // 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 }