Back to index

salome-med  6.5.0
MED_GaussDef.cxx
Go to the documentation of this file.
00001 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D, OPEN CASCADE
00002 //
00003 // This library is free software; you can redistribute it and/or
00004 // modify it under the terms of the GNU Lesser General Public
00005 // License as published by the Free Software Foundation; either
00006 // version 2.1 of the License.
00007 //
00008 // This library is distributed in the hope that it will be useful,
00009 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00010 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00011 // Lesser General Public License for more details.
00012 //
00013 // You should have received a copy of the GNU Lesser General Public
00014 // License along with this library; if not, write to the Free Software
00015 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
00016 //
00017 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
00018 //
00019 
00020 // File   : MED_GaussDef.hxx
00021 // Module : MED
00022 // Author : Edward AGAPOV (eap)
00023 // $Header:
00024 //
00025 #include "MED_GaussDef.hxx"
00026 #include "MED_Utilities.hxx"
00027 #include "MED_GaussUtils.hxx"
00028 
00029 namespace MED
00030 {
00031   using namespace std;
00032   using namespace MED;
00033   //---------------------------------------------------------------
00034 
00035   void TGaussDef::add(const double x, const double weight)
00036   {
00037     if ( dim() != 1 )
00038       EXCEPTION( logic_error,"dim() != 1");
00039     if ( myWeights.capacity() == myWeights.size() )
00040       EXCEPTION( logic_error,"Extra gauss point");
00041     myCoords.push_back( x );
00042     myWeights.push_back( weight );
00043   }
00044   void TGaussDef::add(const double x, const double y, const double weight)
00045   {
00046     if ( dim() != 2 )
00047       EXCEPTION( logic_error,"dim() != 2");
00048     if ( myWeights.capacity() == myWeights.size() )
00049       EXCEPTION( logic_error,"Extra gauss point");
00050     myCoords.push_back( x );
00051     myCoords.push_back( y );
00052     myWeights.push_back( weight );
00053   }
00054   void TGaussDef::add(const double x, const double y, const double z, const double weight)
00055   {
00056     if ( dim() != 3 )
00057       EXCEPTION( logic_error,"dim() != 3");
00058     if ( myWeights.capacity() == myWeights.size() )
00059       EXCEPTION( logic_error,"Extra gauss point");
00060     myCoords.push_back( x );
00061     myCoords.push_back( y );
00062     myCoords.push_back( z );
00063     myWeights.push_back( weight );
00064   }
00065   void TGaussDef::setRefCoords(const TShapeFun& aShapeFun)
00066   {
00067     myRefCoords.reserve( aShapeFun.myRefCoord.size() );
00068     myRefCoords.assign( aShapeFun.myRefCoord.begin(),
00069                         aShapeFun.myRefCoord.end() );
00070   }
00071 
00072 
00073   //---------------------------------------------------------------
00077   //---------------------------------------------------------------
00078 
00079   TGaussDef::TGaussDef(const int geom, const int nbGauss, const int variant)
00080   {
00081     myType = geom;
00082     myCoords .reserve( nbGauss * dim() );
00083     myWeights.reserve( nbGauss );
00084 
00085     switch ( geom ) {
00086 
00087     case eSEG2:
00088     case eSEG3:
00089       if (geom == eSEG2) setRefCoords( TSeg2a() );
00090       else               setRefCoords( TSeg3a() );
00091       switch ( nbGauss ) {
00092       case 1: {
00093         add( 0.0, 2.0 ); break;
00094       }
00095       case 2: {
00096         const double a = 0.577350269189626;
00097         add(  a,  1.0 );
00098         add( -a,  1.0 ); break;
00099       }
00100       case 3: {
00101         const double a = 0.774596669241;
00102         const double P1 = 1./1.8;
00103         const double P2 = 1./1.125;
00104         add( -a,  P1 );
00105         add(  0,  P2 ); 
00106         add(  a,  P1 ); break;
00107       }
00108       case 4: {
00109         const double a  = 0.339981043584856, b  = 0.861136311594053;
00110         const double P1 = 0.652145154862546, P2 = 0.347854845137454 ;
00111         add(  a,  P1 );
00112         add( -a,  P1 );
00113         add(  b,  P2 ); 
00114         add( -b,  P2 ); break;
00115       }
00116       default:
00117         EXCEPTION( logic_error,"Invalid nb of gauss points for SEG"<<nbGauss);
00118       }
00119       break;
00120 
00121     case eTRIA3:
00122     case eTRIA6:
00123       if ( variant == 1 ) {
00124         if (geom == eTRIA3) setRefCoords( TTria3b() );
00125         else                setRefCoords( TTria6b() );
00126         switch ( nbGauss ) {
00127         case 1: { // FPG1
00128           add( 1/3., 1/3., 1/2. ); break;
00129         }
00130         case 3: { // FPG3
00131           // what about COT3 ???
00132           add( 1/6., 1/6., 1/6. );
00133           add( 2/3., 1/6., 1/6. );
00134           add( 1/6., 2/3., 1/6. ); break;
00135         }
00136         case 4: { // FPG4
00137           add( 1/5., 1/5.,  25/(24*4.) );
00138           add( 3/5., 1/5.,  25/(24*4.) );
00139           add( 1/5., 3/5.,  25/(24*4.) );
00140           add( 1/3., 1/3., -27/(24*4.) ); break;
00141         }
00142         case 6: { // FPG6
00143           const double P1 = 0.11169079483905, P2 = 0.0549758718227661;
00144           const double a  = 0.445948490915965, b = 0.091576213509771;
00145           add(     b,     b, P2 ); 
00146           add( 1-2*b,     b, P2 );
00147           add(     b, 1-2*b, P2 );
00148           add(     a, 1-2*a, P1 );
00149           add(     a,     a, P1 ); 
00150           add( 1-2*a,     a, P1 ); break;
00151         }
00152         case 7: { // FPG7
00153           const double A  = 0.470142064105115;
00154           const double B  = 0.101286507323456;
00155           const double P1 = 0.066197076394253;
00156           const double P2 = 0.062969590272413;
00157           add(  1/3.,  1/3., 9/80. ); 
00158           add(     A,     A, P1 ); 
00159           add( 1-2*A,     A, P1 );
00160           add(     A, 1-2*A, P1 );
00161           add(     B,     B, P2 ); 
00162           add( 1-2*B,     B, P2 );
00163           add(     B, 1-2*B, P2 ); break;
00164         }
00165         case 12: { // FPG12
00166           const double A  = 0.063089014491502;
00167           const double B  = 0.249286745170910;
00168           const double C  = 0.310352451033785;
00169           const double D  = 0.053145049844816;
00170           const double P1 = 0.025422453185103;
00171           const double P2 = 0.058393137863189;
00172           const double P3 = 0.041425537809187;
00173           add(     A,     A, P1 ); 
00174           add( 1-2*A,     A, P1 );
00175           add(     A, 1-2*A, P1 );
00176           add(     B,     B, P2 ); 
00177           add( 1-2*B,     B, P2 );
00178           add(     B, 1-2*B, P2 );
00179           add(     C,     D, P3 );
00180           add(     D,     C, P3 );
00181           add( 1-C-D,     C, P3 );
00182           add( 1-C-D,     D, P3 );
00183           add(     C, 1-C-D, P3 );
00184           add(     D, 1-C-D, P3 ); break;
00185         }
00186         default:
00187           EXCEPTION( logic_error,"Invalid nb of gauss points for TRIA, variant 1: "
00188                      <<nbGauss);
00189         }
00190       }
00191       else if ( variant == 2 ) {
00192         if (geom == eTRIA3) setRefCoords( TTria3a() );
00193         else                setRefCoords( TTria6a() );
00194         switch ( nbGauss ) {
00195         case 1: {
00196           add( -1/3., -1/3., 2. ); break;
00197         }
00198         case 3: {
00199           add( -2/3.,  1/3., 2/3. );
00200           add( -2/3., -2/3., 2/3. );
00201           add(  1/3., -2/3., 2/3. ); break;
00202         }
00203         case 6: {
00204           const double P1 = 0.11169079483905, P2 = 0.0549758718227661;
00205           const double A  = 0.445948490915965, B = 0.091576213509771;
00206           add( 2*B-1, 1-4*B, 4*P2 ); 
00207           add( 2*B-1, 2*B-1, 4*P2 );
00208           add( 1-4*B, 2*B-1, 4*P2 );
00209           add( 1-4*A, 2*A-1, 4*P1 );
00210           add( 2*A-1, 1-4*A, 4*P1 ); 
00211           add( 2*A-1, 2*A-1, 4*P1 ); break;
00212         }
00213         default:
00214           EXCEPTION( logic_error,"Invalid nb of gauss points for TRIA, variant 2: "
00215                      <<nbGauss);
00216         }
00217       }
00218       else if ( variant == 3 ) {
00219         if (geom == eTRIA3) setRefCoords( TTria3b() );
00220         else                setRefCoords( TTria6b() );
00221         switch ( nbGauss ) {
00222         case 4: {
00223           add( 1/3., 1/3., -27/96 );
00224           add( 0.2 , 0.2 ,  25/96 );
00225           add( 0.6 , 0.2 ,  25/96 );
00226           add( 0.2 , 0.6 ,  25/96 ); break;
00227         }
00228         default:
00229           EXCEPTION( logic_error,"Invalid nb of gauss points for TRIA, variant 3: "
00230                      <<nbGauss);
00231         }
00232       }
00233       break;
00234 
00235     case eQUAD4:
00236     case eQUAD8:
00237       if ( variant == 1 ) {
00238         if (geom == eQUAD4) setRefCoords( TQuad4b() );
00239         else                setRefCoords( TQuad8b() );
00240         switch ( nbGauss ) {
00241         case 1: { // FPG1
00242           add(  0,  0,  4 ); break;
00243         }
00244         case 4: { // FPG4
00245           const double a = 1/sqrt(3.);
00246           add( -a, -a,  1 );
00247           add(  a, -a,  1 );
00248           add(  a,  a,  1 );
00249           add( -a,  a,  1 ); break;
00250         }
00251         case 9: { // FPG9
00252           const double a = 0.774596669241483;
00253           add( -a, -a,  25/81. );
00254           add(  a, -a,  25/81. );
00255           add(  a,  a,  25/81. );
00256           add( -a,  a,  25/81. );
00257           add( 0., -a,  40/81. );
00258           add(  a, 0.,  40/81. );
00259           add( 0.,  a,  40/81. );
00260           add( -a, 0.,  40/81. );
00261           add( 0., 0.,  64/81. ); break;
00262         }
00263         default:
00264           EXCEPTION( logic_error,"Invalid nb of gauss points for QUAD, variant 1: "
00265                      <<nbGauss);
00266         }
00267       }
00268       else if ( variant == 2 ) {
00269         if (geom == eQUAD4) setRefCoords( TQuad4a() );
00270         else                setRefCoords( TQuad8a() );
00271         switch ( nbGauss ) {
00272         case 4: {
00273           const double a = 1/sqrt(3.);
00274           add( -a,  a,  1 );
00275           add( -a, -a,  1 );
00276           add(  a, -a,  1 );
00277           add(  a,  a,  1 ); break;
00278         }
00279         case 9: {
00280           const double a = 0.774596669241483;
00281           add( -a,  a,  25/81. );
00282           add( -a, -a,  25/81. );
00283           add(  a, -a,  25/81. );
00284           add(  a,  a,  25/81. );
00285           add( -a, 0.,  40/81. );
00286           add( 0., -a,  40/81. );
00287           add(  a, 0.,  40/81. );
00288           add( 0.,  a,  40/81. );
00289           add( 0., 0.,  64/81. ); break;
00290         }
00291         default:
00292           EXCEPTION( logic_error,"Invalid nb of gauss points for QUAD, variant 1: "
00293                      <<nbGauss);
00294         }
00295       }
00296       else if ( variant == 3 ) {
00297         if (geom == eQUAD4) setRefCoords( TQuad4b() );
00298         else                setRefCoords( TQuad8b() );
00299         switch ( nbGauss ) {
00300         case 4: {
00301           const double a = 3/sqrt(3.);
00302           add( -a, -a,  1 );
00303           add( -a,  a,  1 );
00304           add(  a, -a,  1 );
00305           add(  a,  a,  1 ); break;
00306         }
00307         case 9: {
00308           const double a = sqrt(3/5.), c1 = 5/9., c2 = 8/9.;
00309           const double c12 = c1*c2, c22 = c2*c2, c1c2 = c1*c2;
00310           add( -a, -a,  c12  );
00311           add( -a, 0.,  c1c2 );
00312           add( -a,  a,  c12  );
00313           add( 0., -a,  c1c2 );
00314           add( 0., 0.,  c22  );
00315           add( 0.,  a,  c1c2 );
00316           add(  a, -a,  c12  );
00317           add(  a, 0.,  c1c2 );
00318           add(  a,  a,  c12  ); break;
00319         }
00320         default:
00321           EXCEPTION( logic_error,"Invalid nb of gauss points for QUAD, variant 3: "
00322                      <<nbGauss);
00323         }
00324       }
00325       break;
00326 
00327     case eTETRA4:
00328     case eTETRA10:
00329       if (geom == eTETRA4) setRefCoords( TTetra4a() );
00330       else                 setRefCoords( TTetra10a() );
00331       switch ( nbGauss ) {
00332       case 4: { // FPG4
00333         const double a = (5 - sqrt(5.))/20., b = (5 + 3*sqrt(5.))/20.;
00334         add(  a,  a,  a,  1/24. );
00335         add(  a,  a,  b,  1/24. );
00336         add(  a,  b,  a,  1/24. );
00337         add(  b,  a,  a,  1/24. ); break;
00338       }
00339       case 5: { // FPG5
00340         const double a = 0.25, b = 1/6., c = 0.5;
00341         add(  a,  a,  a, -2/15. );
00342         add(  b,  b,  b,  3/40. );
00343         add(  b,  b,  c,  3/40. );
00344         add(  b,  c,  b,  3/40. );
00345         add(  c,  b,  b,  3/40. ); break;
00346       }
00347       case 15: { // FPG15
00348         const double a = 0.25;
00349         const double b1 = (7 + sqrt(15.))/34., c1 = (13 + 3*sqrt(15.))/34., d = (5 - sqrt(15.))/20.;
00350         const double b2 = (7 - sqrt(15.))/34., c2 = (13 - 3*sqrt(15.))/34., e = (5 + sqrt(15.))/20.;
00351         const double P1 = (2665 - 14*sqrt(15.))/226800.;
00352         const double P2 = (2665 + 14*sqrt(15.))/226800.;
00353         add(  a,  a,  a,  8/405.);//_____
00354         add( b1, b1, b1,  P1    );
00355         add( b1, b1, c1,  P1    );
00356         add( b1, c1, b1,  P1    );
00357         add( c1, b1, b1,  P1    );//_____
00358         add( b2, b2, b2,  P2    );
00359         add( b2, b2, c2,  P2    );
00360         add( b2, c2, b2,  P2    );
00361         add( c2, b2, b2,  P2    );//_____
00362         add(  d,  d,  e,  5/567.);
00363         add(  d,  e,  d,  5/567.);
00364         add(  e,  d,  d,  5/567.);
00365         add(  d,  e,  e,  5/567.);
00366         add(  e,  d,  e,  5/567.);
00367         add(  e,  e,  d,  5/567.);
00368         break;
00369       }
00370       default:
00371         EXCEPTION( logic_error,"Invalid nb of gauss points for TETRA: "<<nbGauss);
00372       }
00373       break;
00374 
00375     case ePYRA5:
00376     case ePYRA13:
00377       if (geom == ePYRA5) setRefCoords( TPyra5a() );
00378       else                setRefCoords( TPyra13a() );
00379       switch ( nbGauss ) {
00380       case 5: { // FPG5
00381         const double h1 = 0.1531754163448146;
00382         const double h2 = 0.6372983346207416;
00383         add(  .5,  0.,  h1,  2/15. );
00384         add(  0.,  .5,  h1,  2/15. );
00385         add( -.5,  0.,  h1,  2/15. );
00386         add(  0., -.5,  h1,  2/15. );
00387         add(  0.,  0.,  h2,  2/15. ); break;
00388       }
00389       case 6: { // FPG6
00390         const double p1 = 0.1024890634400000 ;
00391         const double p2 = 0.1100000000000000 ;
00392         const double p3 = 0.1467104129066667 ;
00393         const double a  = 0.5702963741068025 ;
00394         const double h1 = 0.1666666666666666 ;
00395         const double h2 = 0.08063183038464675;
00396         const double h3 = 0.6098484849057127 ;
00397         add(  a, 0.,  h1,  p1 );
00398         add( 0.,  a,  h1,  p1 );
00399         add( -a, 0.,  h1,  p1 );
00400         add( 0., -a,  h1,  p1 );
00401         add( 0., 0.,  h2,  p2 );
00402         add( 0., 0.,  h3,  p3 ); break;
00403       }
00404       case 27: { // FPG27
00405         const double a1  = 0.788073483; 
00406         const double b6  = 0.499369002; 
00407         const double b1  = 0.848418011; 
00408         const double c8  = 0.478508449; 
00409         const double c1  = 0.652816472; 
00410         const double d12 = 0.032303742; 
00411         const double d1  = 1.106412899;
00412         double z = 1/2., fz = b1/2*(1 - z);
00413         add(  0.,  0.,   z,  a1 ); // 1
00414         add(  fz,  fz,   z,  b6 ); // 2
00415         add( -fz,  fz,   z,  b6 ); // 3
00416         add( -fz, -fz,   z,  b6 ); // 4
00417         add(  fz, -fz,   z,  b6 ); // 5
00418         z = (1 - b1)/2.;
00419         add(  0.,  0.,   z,  b6 ); // 6
00420         z = (1 + b1)/2.;
00421         add(  0.,  0.,   z,  b6 ); // 7
00422         z = (1 - c1)/2.; fz = c1*(1 - z);
00423         add(  fz,  0.,   z,  c8 ); // 8
00424         add(  0.,  fz,   z,  c8 ); // 9
00425         add( -fz,  0.,   z,  c8 ); // 10
00426         add(  0., -fz,   z,  c8 ); // 11
00427         z = (1 + c1)/2.; fz = c1*(1 - z);
00428         add(  fz,  0.,   z,  c8 ); // 12
00429         add(  0.,  fz,   z,  c8 ); // 13
00430         add( -fz,  0.,   z,  c8 ); // 14
00431         add(  0., -fz,   z,  c8 ); // 15
00432         z = (1 - d1)/2., fz = d1/2*(1 - z);
00433         add(  fz,  fz,   z,  d12); // 16
00434         add( -fz,  fz,   z,  d12); // 17
00435         add( -fz, -fz,   z,  d12); // 18
00436         add(  fz, -fz,   z,  d12); // 19
00437         z = 1/2.; fz = d1*(1 - z);
00438         add(  fz,  0.,   z,  d12); // 20
00439         add(  0.,  fz,   z,  d12); // 21
00440         add( -fz,  0.,   z,  d12); // 22
00441         add(  0., -fz,   z,  d12); // 23
00442         z = (1 + d1)/2., fz = d1/2*(1 - z);
00443         add(  fz,  fz,   z,  d12); // 24
00444         add( -fz,  fz,   z,  d12); // 25
00445         add( -fz, -fz,   z,  d12); // 26
00446         add(  fz, -fz,   z,  d12); // 27
00447         break;
00448       }
00449       default:
00450         EXCEPTION( logic_error,"Invalid nb of gauss points for PYRA: "<<nbGauss);
00451       }
00452       break;
00453     case ePENTA6:
00454     case ePENTA15:
00455       if (geom == ePENTA6) setRefCoords( TPenta6a() );
00456       else                 setRefCoords( TPenta15a() );
00457       switch ( nbGauss ) {
00458       case 6: { // FPG6
00459         const double a = sqrt(3.)/3.;
00460         add( -a, .5, .5,  1/6. );
00461         add( -a, 0., .5,  1/6. );
00462         add( -a, .5, 0.,  1/6. );
00463         add(  a, .5, .5,  1/6. );
00464         add(  a, 0., .5,  1/6. );
00465         add(  a, .5, 0.,  1/6. ); break;
00466       }
00467       case 8: { // FPG8
00468         const double a = 0.577350269189626;
00469         add( -a, 1/3., 1/3., -27/96. );
00470         add( -a,  0.6,  0.2,  25/96. );
00471         add( -a,  0.2,  0.6,  25/96. );
00472         add( -a,  0.2,  0.2,  25/96. );
00473         add( +a, 1/3., 1/3., -27/96. );
00474         add( +a,  0.6,  0.2,  25/96. );
00475         add( +a,  0.2,  0.6,  25/96. );
00476         add( +a,  0.2,  0.2,  25/96. ); break;
00477       }
00478       case 21: { // FPG21
00479         const double d = sqrt(3/5.), c1 = 5/9., c2 = 8/9.; // d <=> alfa
00480         const double a = (6 + sqrt(15.))/21.;
00481         const double b = (6 - sqrt(15.))/21.;
00482         const double P1 = (155 + sqrt(15.))/2400.;
00483         const double P2 = (155 - sqrt(15.))/2400.;  //___
00484         add( -d,  1/3.,  1/3., c1*9/80. );//___
00485         add( -d,     a,     a, c1*P1    );
00486         add( -d, 1-2*a,     a, c1*P1    );
00487         add( -d,     a, 1-2*a, c1*P1    );//___
00488         add( -d,     b,     b, c1*P2    );
00489         add( -d, 1-2*b,     b, c1*P2    );
00490         add( -d,     b, 1-2*b, c1*P2    );//___
00491         add( 0.,  1/3.,  1/3., c2*9/80. );//___
00492         add( 0.,     a,     a, c2*P1    );
00493         add( 0., 1-2*a,     a, c2*P1    );
00494         add( 0.,     a, 1-2*a, c2*P1    );//___
00495         add( 0.,     b,     b, c2*P2    );
00496         add( 0., 1-2*b,     b, c2*P2    );
00497         add( 0.,     b, 1-2*b, c2*P2    );//___
00498         add(  d,  1/3.,  1/3., c1*9/80. );//___
00499         add(  d,     a,     a, c1*P1    );
00500         add(  d, 1-2*a,     a, c1*P1    );
00501         add(  d,     a, 1-2*a, c1*P1    );//___
00502         add(  d,     b,     b, c1*P2    );
00503         add(  d, 1-2*b,     b, c1*P2    );
00504         add(  d,     b, 1-2*b, c1*P2    );//___
00505         break;
00506       }
00507       default:
00508         EXCEPTION( logic_error,"Invalid nb of gauss points for PENTA: " <<nbGauss);
00509       }
00510       break;
00511 
00512     case eHEXA8:
00513     case eHEXA20:
00514       if (geom == eHEXA8) setRefCoords( THexa8a() );
00515       else                setRefCoords( THexa20a() );
00516       switch ( nbGauss ) {
00517       case 8: { // FPG8
00518         const double a = sqrt(3.)/3.;
00519         add( -a, -a, -a,  1. );
00520         add( -a, -a,  a,  1. );
00521         add( -a,  a, -a,  1. );
00522         add( -a,  a,  a,  1. );
00523         add(  a, -a, -a,  1. );
00524         add(  a, -a,  a,  1. );
00525         add(  a,  a, -a,  1. );
00526         add(  a,  a,  a,  1. ); break;
00527       }
00528       case 27: { // FPG27
00529         const double a = sqrt(3/5.), c1 = 5/9., c2 = 8/9.;
00530         const double c12 = c1*c1, c13 = c1*c1*c1;
00531         const double c22 = c2*c2, c23 = c2*c2*c2;
00532         add( -a, -a, -a,   c13  ); // 1
00533         add( -a, -a, 0., c12*c2 ); // 2
00534         add( -a, -a,  a,   c13  ); // 3
00535         add( -a, 0., -a, c12*c2 ); // 4
00536         add( -a, 0., 0., c1*c22 ); // 5
00537         add( -a, 0.,  a, c12*c2 ); // 6
00538         add( -a,  a, -a,   c13  ); // 7
00539         add( -a,  a, 0., c12*c2 ); // 8
00540         add( -a,  a,  a,   c13  ); // 9
00541         add( 0., -a, -a, c12*c2 ); // 10
00542         add( 0., -a, 0., c1*c22 ); // 11
00543         add( 0., -a,  a, c12*c2 ); // 12
00544         add( 0., 0., -a, c1*c22 ); // 13
00545         add( 0., 0., 0.,   c23  ); // 14
00546         add( 0., 0.,  a, c1*c22 ); // 15
00547         add( 0.,  a, -a, c12*c2 ); // 16
00548         add( 0.,  a, 0., c1*c22 ); // 17
00549         add( 0.,  a,  a, c12*c2 ); // 18
00550         add(  a, -a, -a,   c13  ); // 19
00551         add(  a, -a, 0., c12*c2 ); // 20
00552         add(  a, -a,  a,   c13  ); // 21
00553         add(  a, 0., -a, c12*c2 ); // 22
00554         add(  a, 0., 0., c1*c22 ); // 23
00555         add(  a, 0.,  a, c12*c2 ); // 24
00556         add(  a,  a, -a,   c13  ); // 25
00557         add(  a,  a, 0., c12*c2 ); // 26
00558         add(  a,  a,  a,   c13  ); // 27
00559         break;
00560       }
00561       default:
00562         EXCEPTION( logic_error,"Invalid nb of gauss points for PENTA: " <<nbGauss);
00563       }
00564       break;
00565 
00566     default:
00567       EXCEPTION( logic_error,"unexpected EGeometrieElement: "<< geom);
00568     }
00569 
00570     if ( myWeights.capacity() != myWeights.size() )
00571       EXCEPTION( logic_error,"Not all gauss points defined");
00572   }
00573 }