Back to index

salome-med  6.5.0
MEDMEM_GaussLocalization.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      : MEDMEM_GaussLocalization.cxx
00021 // Created   : Thu Dec 20 12:26:33 2007
00022 // Author    : Edward AGAPOV (eap)
00023 //
00024 #include "MEDMEM_GaussLocalization.hxx"
00025 
00026 #include <vector>
00027 #include <math.h>
00028 
00029 #define EXCEPTION(type,msg) \
00030   MEDEXCEPTION(LOCALIZED(MEDMEM::STRING(#type)<< msg))
00031 
00032 namespace // matters used by GAUSS_LOCALIZATION_::makeDefaultLocalization()
00033 {
00034   typedef std::vector<double> TDoubleVector;
00035   typedef double*             TCoordSlice;
00036   typedef int                 TInt;
00037   //---------------------------------------------------------------
00039   //---------------------------------------------------------------
00040   struct TShapeFun
00041   {
00042     TInt myDim;
00043     TInt myNbRef;
00044     TDoubleVector myRefCoord;
00045 
00046     TShapeFun(TInt theDim = 0, TInt theNbRef = 0)
00047       :myDim(theDim),myNbRef(theNbRef),myRefCoord(theNbRef*theDim) {}
00048 
00049     TInt GetNbRef() const { return myNbRef; }
00050 
00051     TCoordSlice GetCoord(TInt theRefId) { return &myRefCoord[0] + theRefId * myDim; }
00052   };
00053   //---------------------------------------------------------------
00057   //---------------------------------------------------------------
00058   struct TGaussDef
00059   {
00060     int           myType;      
00061     TDoubleVector myRefCoords; 
00062     TDoubleVector myCoords;    
00063     TDoubleVector myWeights;   
00064 
00077     TGaussDef(const int geomType, const int nbPoints, const int variant=1);
00078 
00079     int dim() const { return myType/100; }
00080     int nbPoints() const { return myWeights.capacity(); }
00081 
00082   private:
00083     void add(const double x, const double weight);
00084     void add(const double x, const double y, const double weight);
00085     void add(const double x, const double y, const double z, const double weight);
00086     void setRefCoords(const TShapeFun& aShapeFun) { myRefCoords = aShapeFun.myRefCoord; }
00087   };
00088   struct TSeg2a: TShapeFun {
00089     TSeg2a();
00090   };
00091   struct TSeg3a: TShapeFun {
00092     TSeg3a();
00093   };
00094   struct TTria3a: TShapeFun {
00095     TTria3a();
00096   };
00097   struct TTria6a: TShapeFun {
00098     TTria6a();
00099   };
00100   struct TTria3b: TShapeFun {
00101     TTria3b();
00102   };
00103   struct TTria6b: TShapeFun {
00104     TTria6b();
00105   };
00106   struct TQuad4a: TShapeFun {
00107     TQuad4a();
00108   };
00109   struct TQuad8a: TShapeFun {
00110     TQuad8a();
00111   };
00112   struct TQuad4b: TShapeFun {
00113     TQuad4b();
00114   };
00115   struct TQuad8b: TShapeFun {
00116     TQuad8b();
00117   };
00118   struct TTetra4a: TShapeFun {
00119     TTetra4a();
00120   };
00121   struct TTetra10a: TShapeFun {
00122     TTetra10a();
00123   };
00124   struct TTetra4b: TShapeFun {
00125     TTetra4b();
00126   };
00127   struct TTetra10b: TShapeFun {
00128     TTetra10b();
00129   };
00130   struct THexa8a: TShapeFun {
00131     THexa8a();
00132   };
00133   struct THexa20a: TShapeFun {
00134     THexa20a(TInt theDim = 3, TInt theNbRef = 20);
00135   };
00136   struct THexa27a: THexa20a {
00137     THexa27a();
00138   };
00139   struct THexa8b: TShapeFun {
00140     THexa8b();
00141   };
00142   struct THexa20b: TShapeFun {
00143     THexa20b(TInt theDim = 3, TInt theNbRef = 20);
00144   };
00145   struct TPenta6a: TShapeFun {
00146     TPenta6a();
00147   };
00148   struct TPenta6b: TShapeFun {
00149     TPenta6b();
00150   };
00151   struct TPenta15a: TShapeFun {
00152     TPenta15a();
00153   };
00154   struct TPenta15b: TShapeFun {
00155     TPenta15b();
00156   };
00157   struct TPyra5a: TShapeFun {
00158     TPyra5a();
00159   };
00160   struct TPyra5b: TShapeFun {
00161     TPyra5b();
00162   };
00163   struct TPyra13a: TShapeFun {
00164     TPyra13a();
00165   };
00166   struct TPyra13b: TShapeFun {
00167     TPyra13b();
00168   };
00169 
00170   void TGaussDef::add(const double x, const double weight)
00171   {
00172     if ( dim() != 1 )
00173       EXCEPTION( logic_error,"dim() != 1");
00174     if ( myWeights.capacity() == myWeights.size() )
00175       EXCEPTION( logic_error,"Extra gauss point");
00176     myCoords.push_back( x );
00177     myWeights.push_back( weight );
00178   }
00179   void TGaussDef::add(const double x, const double y, const double weight)
00180   {
00181     if ( dim() != 2 )
00182       EXCEPTION( logic_error,"dim() != 2");
00183     if ( myWeights.capacity() == myWeights.size() )
00184       EXCEPTION( logic_error,"Extra gauss point");
00185     myCoords.push_back( x );
00186     myCoords.push_back( y );
00187     myWeights.push_back( weight );
00188   }
00189   void TGaussDef::add(const double x, const double y, const double z, const double weight)
00190   {
00191     if ( dim() != 3 )
00192       EXCEPTION( logic_error,"dim() != 3");
00193     if ( myWeights.capacity() == myWeights.size() )
00194       EXCEPTION( logic_error,"Extra gauss point");
00195     myCoords.push_back( x );
00196     myCoords.push_back( y );
00197     myCoords.push_back( z );
00198     myWeights.push_back( weight );
00199   }
00200 
00201   //---------------------------------------------------------------
00202   TSeg2a::TSeg2a():TShapeFun(1,2)
00203   {
00204     TInt aNbRef = GetNbRef();
00205     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
00206       TCoordSlice aCoord = GetCoord(aRefId);
00207       switch(aRefId){
00208       case  0: aCoord[0] = -1.0; break;
00209       case  1: aCoord[0] =  1.0; break;
00210       }
00211     }
00212   }
00213   //---------------------------------------------------------------
00214   TSeg3a::TSeg3a():TShapeFun(1,3)
00215   {
00216     TInt aNbRef = GetNbRef();
00217     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
00218       TCoordSlice aCoord = GetCoord(aRefId);
00219       switch(aRefId){
00220       case  0: aCoord[0] = -1.0; break;
00221       case  1: aCoord[0] =  1.0; break;
00222       case  2: aCoord[0] =  0.0; break;
00223       }
00224     }
00225   }
00226   //---------------------------------------------------------------
00227   TTria3a::TTria3a():
00228     TShapeFun(2,3)
00229   {
00230     TInt aNbRef = GetNbRef();
00231     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
00232       TCoordSlice aCoord = GetCoord(aRefId);
00233       switch(aRefId){
00234       case  0: aCoord[0] = -1.0;  aCoord[1] =  1.0; break;
00235       case  1: aCoord[0] = -1.0;  aCoord[1] = -1.0; break;
00236       case  2: aCoord[0] =  1.0;  aCoord[1] = -1.0; break;
00237       }
00238     }
00239   }
00240   //---------------------------------------------------------------
00241   TTria6a::TTria6a():TShapeFun(2,6)
00242   {
00243     TInt aNbRef = GetNbRef();
00244     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
00245       TCoordSlice aCoord = GetCoord(aRefId);
00246       switch(aRefId){
00247       case  0: aCoord[0] = -1.0;  aCoord[1] =  1.0; break;
00248       case  1: aCoord[0] = -1.0;  aCoord[1] = -1.0; break;
00249       case  2: aCoord[0] =  1.0;  aCoord[1] = -1.0; break;
00250 
00251       case  3: aCoord[0] = -1.0;  aCoord[1] =  1.0; break;
00252       case  4: aCoord[0] =  0.0;  aCoord[1] = -1.0; break;
00253       case  5: aCoord[0] =  0.0;  aCoord[1] =  0.0; break;
00254       }
00255     }
00256   }
00257   //---------------------------------------------------------------
00258   TTria3b::TTria3b():
00259     TShapeFun(2,3)
00260   {
00261     TInt aNbRef = GetNbRef();
00262     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
00263       TCoordSlice aCoord = GetCoord(aRefId);
00264       switch(aRefId){
00265       case  0: aCoord[0] =  0.0;  aCoord[1] =  0.0; break;
00266       case  1: aCoord[0] =  1.0;  aCoord[1] =  0.0; break;
00267       case  2: aCoord[0] =  0.0;  aCoord[1] =  1.0; break;
00268       }
00269     }
00270   }
00271   //---------------------------------------------------------------
00272   TTria6b::TTria6b():
00273     TShapeFun(2,6)
00274   {
00275     TInt aNbRef = GetNbRef();
00276     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
00277       TCoordSlice aCoord = GetCoord(aRefId);
00278       switch(aRefId){
00279       case  0: aCoord[0] =  0.0;  aCoord[1] =  0.0; break;
00280       case  1: aCoord[0] =  1.0;  aCoord[1] =  0.0; break;
00281       case  2: aCoord[0] =  0.0;  aCoord[1] =  1.0; break;
00282 
00283       case  3: aCoord[0] =  0.5;  aCoord[1] =  0.0; break;
00284       case  4: aCoord[0] =  0.5;  aCoord[1] =  0.5; break;
00285       case  5: aCoord[0] =  0.0;  aCoord[1] =  0.5; break;
00286       }
00287     }
00288   }
00289   //---------------------------------------------------------------
00290   TQuad4a::TQuad4a():
00291     TShapeFun(2,4)
00292   {
00293     TInt aNbRef = GetNbRef();
00294     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
00295       TCoordSlice aCoord = GetCoord(aRefId);
00296       switch(aRefId){
00297       case  0: aCoord[0] = -1.0;  aCoord[1] =  1.0; break;
00298       case  1: aCoord[0] = -1.0;  aCoord[1] = -1.0; break;
00299       case  2: aCoord[0] =  1.0;  aCoord[1] = -1.0; break;
00300       case  3: aCoord[0] =  1.0;  aCoord[1] =  1.0; break;
00301       }
00302     }
00303   }
00304   //---------------------------------------------------------------
00305   TQuad8a::TQuad8a():
00306     TShapeFun(2,8)
00307   {
00308     TInt aNbRef = GetNbRef();
00309     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
00310       TCoordSlice aCoord = GetCoord(aRefId);
00311       switch(aRefId){
00312       case  0: aCoord[0] = -1.0;  aCoord[1] =  1.0; break;
00313       case  1: aCoord[0] = -1.0;  aCoord[1] = -1.0; break;
00314       case  2: aCoord[0] =  1.0;  aCoord[1] = -1.0; break;
00315       case  3: aCoord[0] =  1.0;  aCoord[1] =  1.0; break;
00316 
00317       case  4: aCoord[0] = -1.0;  aCoord[1] =  0.0; break;
00318       case  5: aCoord[0] =  0.0;  aCoord[1] = -1.0; break;
00319       case  6: aCoord[0] =  1.0;  aCoord[1] =  0.0; break;
00320       case  7: aCoord[0] =  0.0;  aCoord[1] =  1.0; break;
00321       }
00322     }
00323   }
00324   //---------------------------------------------------------------
00325   TQuad4b::TQuad4b():
00326     TShapeFun(2,4)
00327   {
00328     TInt aNbRef = GetNbRef();
00329     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
00330       TCoordSlice aCoord = GetCoord(aRefId);
00331       switch(aRefId){
00332       case  0: aCoord[0] = -1.0;  aCoord[1] = -1.0; break;
00333       case  1: aCoord[0] =  1.0;  aCoord[1] = -1.0; break;
00334       case  2: aCoord[0] =  1.0;  aCoord[1] =  1.0; break;
00335       case  3: aCoord[0] = -1.0;  aCoord[1] =  1.0; break;
00336       }
00337     }
00338   }
00339   //---------------------------------------------------------------
00340   TQuad8b::TQuad8b():
00341     TShapeFun(2,8)
00342   {
00343     TInt aNbRef = GetNbRef();
00344     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
00345       TCoordSlice aCoord = GetCoord(aRefId);
00346       switch(aRefId){
00347       case  0: aCoord[0] = -1.0;  aCoord[1] = -1.0; break;
00348       case  1: aCoord[0] =  1.0;  aCoord[1] = -1.0; break;
00349       case  2: aCoord[0] =  1.0;  aCoord[1] =  1.0; break;
00350       case  3: aCoord[0] = -1.0;  aCoord[1] =  1.0; break;
00351 
00352       case  4: aCoord[0] =  0.0;  aCoord[1] = -1.0; break;
00353       case  5: aCoord[0] =  1.0;  aCoord[1] =  0.0; break;
00354       case  6: aCoord[0] =  0.0;  aCoord[1] =  1.0; break;
00355       case  7: aCoord[0] = -1.0;  aCoord[1] =  0.0; break;
00356       }
00357     }
00358   }
00359   //---------------------------------------------------------------
00360   TTetra4a::TTetra4a():
00361     TShapeFun(3,4)
00362   {
00363     TInt aNbRef = GetNbRef();
00364     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
00365       TCoordSlice aCoord = GetCoord(aRefId);
00366       switch(aRefId){
00367       case  0: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
00368       case  1: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
00369       case  2: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
00370       case  3: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
00371       }
00372     }
00373   }
00374   //---------------------------------------------------------------
00375   TTetra10a::TTetra10a():
00376     TShapeFun(3,10)
00377   {
00378     TInt aNbRef = GetNbRef();
00379     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
00380       TCoordSlice aCoord = GetCoord(aRefId);
00381       switch(aRefId){
00382       case  0: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
00383       case  1: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
00384       case  2: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
00385       case  3: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
00386 
00387       case  4: aCoord[0] =  0.0;  aCoord[1] =  0.5;  aCoord[2] =  0.5; break;
00388       case  5: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
00389       case  6: aCoord[0] =  0.0;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
00390 
00391       case  7: aCoord[0] =  0.5;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
00392       case  8: aCoord[0] =  0.5;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
00393       case  9: aCoord[0] =  0.5;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
00394       }
00395     }
00396   }
00397   //---------------------------------------------------------------
00398   TTetra4b::TTetra4b():
00399     TShapeFun(3,4)
00400   {
00401     TInt aNbRef = GetNbRef();
00402     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
00403       TCoordSlice aCoord = GetCoord(aRefId);
00404       switch(aRefId){
00405       case  0: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
00406       case  2: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
00407       case  1: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
00408       case  3: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
00409       }
00410     }
00411   }
00412   //---------------------------------------------------------------
00413   TTetra10b::TTetra10b():
00414     TShapeFun(3,10)
00415   {
00416     TInt aNbRef = GetNbRef();
00417     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
00418       TCoordSlice aCoord = GetCoord(aRefId);
00419       switch(aRefId){
00420       case  0: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
00421       case  2: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
00422       case  1: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
00423       case  3: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
00424 
00425       case  6: aCoord[0] =  0.0;  aCoord[1] =  0.5;  aCoord[2] =  0.5; break;
00426       case  5: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
00427       case  4: aCoord[0] =  0.0;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
00428 
00429       case  7: aCoord[0] =  0.5;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
00430       case  9: aCoord[0] =  0.5;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
00431       case  8: aCoord[0] =  0.5;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
00432       }
00433     }
00434   }
00435   //---------------------------------------------------------------
00436   THexa8a::THexa8a():
00437     TShapeFun(3,8)
00438   {
00439     TInt aNbRef = GetNbRef();
00440     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
00441       TCoordSlice aCoord = GetCoord(aRefId);
00442       switch(aRefId){
00443       case  0: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
00444       case  1: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
00445       case  2: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
00446       case  3: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
00447       case  4: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
00448       case  5: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
00449       case  6: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
00450       case  7: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
00451       }
00452     }
00453   }
00454   //---------------------------------------------------------------
00455   THexa20a::THexa20a(TInt theDim, TInt theNbRef):
00456     TShapeFun(theDim,theNbRef)
00457   {
00458     TInt aNbRef = myRefCoord.size();
00459     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
00460       TCoordSlice aCoord = GetCoord(aRefId);
00461       switch(aRefId){
00462       case  0: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
00463       case  1: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
00464       case  2: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
00465       case  3: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
00466       case  4: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
00467       case  5: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
00468       case  6: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
00469       case  7: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
00470 
00471       case  8: aCoord[0] =  0.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
00472       case  9: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] = -1.0; break;
00473       case 10: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
00474       case 11: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] = -1.0; break;
00475       case 12: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] =  0.0; break;
00476       case 13: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] =  0.0; break;
00477       case 14: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
00478       case 15: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
00479       case 16: aCoord[0] =  0.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
00480       case 17: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
00481       case 18: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
00482       case 19: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
00483       }
00484     }
00485   }
00486   //---------------------------------------------------------------
00487   THexa27a::THexa27a():
00488     THexa20a(3,27)
00489   {
00490     TInt aNbRef = myRefCoord.size();
00491     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
00492       TCoordSlice aCoord = GetCoord(aRefId);
00493       switch(aRefId){
00494       case 20: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] = -1.0; break;
00495       case 21: aCoord[0] =  0.0;  aCoord[1] = -1.0;  aCoord[2] =  0.0; break;
00496       case 22: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
00497       case 23: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
00498       case 24: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
00499       case 25: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
00500       case 26: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
00501       }
00502     }
00503   }
00504   //---------------------------------------------------------------
00505   THexa8b::THexa8b():
00506     TShapeFun(3,8)
00507   {
00508     TInt aNbRef = GetNbRef();
00509     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
00510       TCoordSlice aCoord = GetCoord(aRefId);
00511       switch(aRefId){
00512       case  0: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
00513       case  3: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
00514       case  2: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
00515       case  1: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
00516       case  4: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
00517       case  7: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
00518       case  6: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
00519       case  5: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
00520       }
00521     }
00522   }
00523   //---------------------------------------------------------------
00524   THexa20b::THexa20b(TInt theDim, TInt theNbRef):
00525     TShapeFun(theDim,theNbRef)
00526   {
00527     TInt aNbRef = myRefCoord.size();
00528     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
00529       TCoordSlice aCoord = GetCoord(aRefId);
00530       switch(aRefId){
00531       case  0: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
00532       case  3: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
00533       case  2: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
00534       case  1: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
00535       case  4: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
00536       case  7: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
00537       case  6: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
00538       case  5: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
00539 
00540       case 11: aCoord[0] =  0.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
00541       case 10: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] = -1.0; break;
00542       case  9: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
00543       case  8: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] = -1.0; break;
00544       case 16: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] =  0.0; break;
00545       case 19: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] =  0.0; break;
00546       case 18: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
00547       case 17: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
00548       case 15: aCoord[0] =  0.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
00549       case 14: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
00550       case 13: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
00551       case 12: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
00552       }
00553     }
00554   }
00555   //---------------------------------------------------------------
00556   TPenta6a::TPenta6a():
00557     TShapeFun(3,6)
00558   {
00559     TInt aNbRef = myRefCoord.size();
00560     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
00561       TCoordSlice aCoord = GetCoord(aRefId);
00562       switch(aRefId){
00563       case  0: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
00564       case  1: aCoord[0] = -1.0;  aCoord[1] = -0.0;  aCoord[2] =  1.0; break;
00565       case  2: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
00566       case  3: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
00567       case  4: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
00568       case  5: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
00569       }
00570     }
00571   }
00572   //---------------------------------------------------------------
00573   TPenta6b::TPenta6b():
00574     TShapeFun(3,6)
00575   {
00576     TInt aNbRef = myRefCoord.size();
00577     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
00578       TCoordSlice aCoord = GetCoord(aRefId);
00579       switch(aRefId){
00580       case  0: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
00581       case  2: aCoord[0] = -1.0;  aCoord[1] = -0.0;  aCoord[2] =  1.0; break;
00582       case  1: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
00583       case  3: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
00584       case  5: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
00585       case  4: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
00586       }
00587     }
00588   }
00589   //---------------------------------------------------------------
00590   TPenta15a::TPenta15a():
00591     TShapeFun(3,15)
00592   {
00593     TInt aNbRef = myRefCoord.size();
00594     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
00595       TCoordSlice aCoord = GetCoord(aRefId);
00596       switch(aRefId){
00597       case  0: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
00598       case  1: aCoord[0] = -1.0;  aCoord[1] = -0.0;  aCoord[2] =  1.0; break;
00599       case  2: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
00600       case  3: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
00601       case  4: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
00602       case  5: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
00603 
00604       case  6: aCoord[0] = -1.0;  aCoord[1] =  0.5;  aCoord[2] =  0.5; break;
00605       case  7: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
00606       case  8: aCoord[0] = -1.0;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
00607       case  9: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
00608       case 10: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
00609       case 11: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
00610       case 12: aCoord[0] =  1.0;  aCoord[1] =  0.5;  aCoord[2] =  0.5; break;
00611       case 13: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
00612       case 14: aCoord[0] =  1.0;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
00613       }
00614     }
00615   }
00616   //---------------------------------------------------------------
00617   TPenta15b::TPenta15b():
00618     TShapeFun(3,15)
00619   {
00620     TInt aNbRef = myRefCoord.size();
00621     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
00622       TCoordSlice aCoord = GetCoord(aRefId);
00623       switch(aRefId){
00624       case  0: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
00625       case  2: aCoord[0] = -1.0;  aCoord[1] = -0.0;  aCoord[2] =  1.0; break;
00626       case  1: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
00627       case  3: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
00628       case  5: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
00629       case  4: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
00630 
00631       case  8: aCoord[0] = -1.0;  aCoord[1] =  0.5;  aCoord[2] =  0.5; break;
00632       case  7: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
00633       case  6: aCoord[0] = -1.0;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
00634       case 12: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
00635       case 14: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
00636       case 13: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
00637       case 11: aCoord[0] =  1.0;  aCoord[1] =  0.5;  aCoord[2] =  0.5; break;
00638       case 10: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
00639       case  9: aCoord[0] =  1.0;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
00640       }
00641     }
00642   }
00643   //---------------------------------------------------------------
00644   TPyra5a::TPyra5a():
00645     TShapeFun(3,5)
00646   {
00647     TInt aNbRef = myRefCoord.size();
00648     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
00649       TCoordSlice aCoord = GetCoord(aRefId);
00650       switch(aRefId){
00651       case  0: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
00652       case  1: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
00653       case  2: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
00654       case  3: aCoord[0] =  0.0;  aCoord[1] = -1.0;  aCoord[2] =  0.0; break;
00655       case  4: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
00656       }
00657     }
00658   }
00659   //---------------------------------------------------------------
00660   TPyra5b::TPyra5b():
00661     TShapeFun(3,5)
00662   {
00663     TInt aNbRef = myRefCoord.size();
00664     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
00665       TCoordSlice aCoord = GetCoord(aRefId);
00666       switch(aRefId){        
00667       case  0: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
00668       case  3: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
00669       case  2: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
00670       case  1: aCoord[0] =  0.0;  aCoord[1] = -1.0;  aCoord[2] =  0.0; break;
00671       case  4: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
00672       }
00673     }
00674   }
00675   //---------------------------------------------------------------
00676   TPyra13a::TPyra13a():
00677     TShapeFun(3,13)
00678   {
00679     TInt aNbRef = myRefCoord.size();
00680     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
00681       TCoordSlice aCoord = GetCoord(aRefId);
00682       switch(aRefId){
00683       case  0: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
00684       case  1: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
00685       case  2: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
00686       case  3: aCoord[0] =  0.0;  aCoord[1] = -1.0;  aCoord[2] =  0.0; break;
00687       case  4: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
00688 
00689       case  5: aCoord[0] =  0.5;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
00690       case  6: aCoord[0] = -0.5;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
00691       case  7: aCoord[0] = -0.5;  aCoord[1] = -0.5;  aCoord[2] =  0.0; break;
00692       case  8: aCoord[0] =  0.5;  aCoord[1] = -0.5;  aCoord[2] =  0.0; break;
00693       case  9: aCoord[0] =  0.5;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
00694       case 10: aCoord[0] =  0.0;  aCoord[1] =  0.5;  aCoord[2] =  0.5; break;
00695       case 11: aCoord[0] = -0.5;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
00696       case 12: aCoord[0] =  0.0;  aCoord[1] = -0.5;  aCoord[2] =  0.5; break;
00697       }
00698     }
00699   }
00700   //---------------------------------------------------------------
00701   TPyra13b::TPyra13b():
00702     TShapeFun(3,13)
00703   {
00704     TInt aNbRef = myRefCoord.size();
00705     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
00706       TCoordSlice aCoord = GetCoord(aRefId);
00707       switch(aRefId){
00708       case  0: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
00709       case  3: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
00710       case  2: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
00711       case  1: aCoord[0] =  0.0;  aCoord[1] = -1.0;  aCoord[2] =  0.0; break;
00712       case  4: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
00713 
00714       case  8: aCoord[0] =  0.5;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
00715       case  7: aCoord[0] = -0.5;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
00716       case  6: aCoord[0] = -0.5;  aCoord[1] = -0.5;  aCoord[2] =  0.0; break;
00717       case  5: aCoord[0] =  0.5;  aCoord[1] = -0.5;  aCoord[2] =  0.0; break;
00718       case  9: aCoord[0] =  0.5;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
00719       case 12: aCoord[0] =  0.0;  aCoord[1] =  0.5;  aCoord[2] =  0.5; break;
00720       case 11: aCoord[0] = -0.5;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
00721       case 10: aCoord[0] =  0.0;  aCoord[1] = -0.5;  aCoord[2] =  0.5; break;
00722       }
00723     }
00724   }
00729   TGaussDef::TGaussDef(const int geom, const int nbGauss, const int variant)
00730   {
00731     myType = geom;
00732     myCoords .reserve( nbGauss * dim() );
00733     myWeights.reserve( nbGauss );
00734 
00735     switch ( geom ) {
00736 
00737     case MED_SEG2:
00738     case MED_SEG3:
00739       if (geom == MED_SEG2) setRefCoords( TSeg2a() );
00740       else               setRefCoords( TSeg3a() );
00741       switch ( nbGauss ) {
00742       case 1: {
00743         add( 0.0, 2.0 ); break;
00744       }
00745       case 2: {
00746         const double a = 0.577350269189626;
00747         add(  a,  1.0 );
00748         add( -a,  1.0 ); break;
00749       }
00750       case 3: {
00751         const double a = 0.774596669241;
00752         const double P1 = 1./1.8;
00753         const double P2 = 1./1.125;
00754         add( -a,  P1 );
00755         add(  0,  P2 ); 
00756         add(  a,  P1 ); break;
00757       }
00758       case 4: {
00759         const double a  = 0.339981043584856, b  = 0.861136311594053;
00760         const double P1 = 0.652145154862546, P2 = 0.347854845137454 ;
00761         add(  a,  P1 );
00762         add( -a,  P1 );
00763         add(  b,  P2 ); 
00764         add( -b,  P2 ); break;
00765       }
00766       default:
00767         EXCEPTION( logic_error,"Invalid nb of gauss points for SEG"<<nbGauss);
00768       }
00769       break;
00770 
00771     case MED_TRIA3:
00772     case MED_TRIA6:
00773       if ( variant == 1 ) {
00774         if (geom == MED_TRIA3) setRefCoords( TTria3b() );
00775         else                setRefCoords( TTria6b() );
00776         switch ( nbGauss ) {
00777         case 1: { // FPG1
00778           add( 1/3., 1/3., 1/2. ); break;
00779         }
00780         case 3: { // FPG3
00781           // what about COT3 ???
00782           add( 1/6., 1/6., 1/6. );
00783           add( 2/3., 1/6., 1/6. );
00784           add( 1/6., 2/3., 1/6. ); break;
00785         }
00786         case 4: { // FPG4
00787           add( 1/5., 1/5.,  25/(24*4.) );
00788           add( 3/5., 1/5.,  25/(24*4.) );
00789           add( 1/5., 3/5.,  25/(24*4.) );
00790           add( 1/3., 1/3., -27/(24*4.) ); break;
00791         }
00792         case 6: { // FPG6
00793           const double P1 = 0.11169079483905, P2 = 0.0549758718227661;
00794           const double a  = 0.445948490915965, b = 0.091576213509771;
00795           add(     b,     b, P2 ); 
00796           add( 1-2*b,     b, P2 );
00797           add(     b, 1-2*b, P2 );
00798           add(     a, 1-2*a, P1 );
00799           add(     a,     a, P1 ); 
00800           add( 1-2*a,     a, P1 ); break;
00801         }
00802         case 7: { // FPG7
00803           const double A  = 0.470142064105115;
00804           const double B  = 0.101286507323456;
00805           const double P1 = 0.066197076394253;
00806           const double P2 = 0.062969590272413;
00807           add(  1/3.,  1/3., 9/80. ); 
00808           add(     A,     A, P1 ); 
00809           add( 1-2*A,     A, P1 );
00810           add(     A, 1-2*A, P1 );
00811           add(     B,     B, P2 ); 
00812           add( 1-2*B,     B, P2 );
00813           add(     B, 1-2*B, P2 ); break;
00814         }
00815         case 12: { // FPG12
00816           const double A  = 0.063089014491502;
00817           const double B  = 0.249286745170910;
00818           const double C  = 0.310352451033785;
00819           const double D  = 0.053145049844816;
00820           const double P1 = 0.025422453185103;
00821           const double P2 = 0.058393137863189;
00822           const double P3 = 0.041425537809187;
00823           add(     A,     A, P1 ); 
00824           add( 1-2*A,     A, P1 );
00825           add(     A, 1-2*A, P1 );
00826           add(     B,     B, P2 ); 
00827           add( 1-2*B,     B, P2 );
00828           add(     B, 1-2*B, P2 );
00829           add(     C,     D, P3 );
00830           add(     D,     C, P3 );
00831           add( 1-C-D,     C, P3 );
00832           add( 1-C-D,     D, P3 );
00833           add(     C, 1-C-D, P3 );
00834           add(     D, 1-C-D, P3 ); break;
00835         }
00836         default:
00837           EXCEPTION( logic_error,"Invalid nb of gauss points for TRIA, variant 1: "
00838                      <<nbGauss);
00839         }
00840       }
00841       else if ( variant == 2 ) {
00842         if (geom == MED_TRIA3) setRefCoords( TTria3a() );
00843         else                setRefCoords( TTria6a() );
00844         switch ( nbGauss ) {
00845         case 1: {
00846           add( -1/3., -1/3., 2. ); break;
00847         }
00848         case 3: {
00849           add( -2/3.,  1/3., 2/3. );
00850           add( -2/3., -2/3., 2/3. );
00851           add(  1/3., -2/3., 2/3. ); break;
00852         }
00853         case 6: {
00854           const double P1 = 0.11169079483905, P2 = 0.0549758718227661;
00855           const double A  = 0.445948490915965, B = 0.091576213509771;
00856           add( 2*B-1, 1-4*B, 4*P2 ); 
00857           add( 2*B-1, 2*B-1, 4*P2 );
00858           add( 1-4*B, 2*B-1, 4*P2 );
00859           add( 1-4*A, 2*A-1, 4*P1 );
00860           add( 2*A-1, 1-4*A, 4*P1 ); 
00861           add( 2*A-1, 2*A-1, 4*P1 ); break;
00862         }
00863         default:
00864           EXCEPTION( logic_error,"Invalid nb of gauss points for TRIA, variant 2: "
00865                      <<nbGauss);
00866         }
00867       }
00868       else if ( variant == 3 ) {
00869         if (geom == MED_TRIA3) setRefCoords( TTria3b() );
00870         else                setRefCoords( TTria6b() );
00871         switch ( nbGauss ) {
00872         case 4: {
00873           add( 1/3., 1/3., -27/96 );
00874           add( 0.2 , 0.2 ,  25/96 );
00875           add( 0.6 , 0.2 ,  25/96 );
00876           add( 0.2 , 0.6 ,  25/96 ); break;
00877         }
00878         default:
00879           EXCEPTION( logic_error,"Invalid nb of gauss points for TRIA, variant 3: "
00880                      <<nbGauss);
00881         }
00882       }
00883       break;
00884 
00885     case MED_QUAD4:
00886     case MED_QUAD8:
00887       if ( variant == 1 ) {
00888         if (geom == MED_QUAD4) setRefCoords( TQuad4b() );
00889         else                setRefCoords( TQuad8b() );
00890         switch ( nbGauss ) {
00891         case 1: { // FPG1
00892           add(  0,  0,  4 ); break;
00893         }
00894         case 4: { // FPG4
00895           const double a = 1/sqrt(3.);
00896           add( -a, -a,  1 );
00897           add(  a, -a,  1 );
00898           add(  a,  a,  1 );
00899           add( -a,  a,  1 ); break;
00900         }
00901         case 9: { // FPG9
00902           const double a = 0.774596669241483;
00903           add( -a, -a,  25/81. );
00904           add(  a, -a,  25/81. );
00905           add(  a,  a,  25/81. );
00906           add( -a,  a,  25/81. );
00907           add( 0., -a,  40/81. );
00908           add(  a, 0.,  40/81. );
00909           add( 0.,  a,  40/81. );
00910           add( -a, 0.,  40/81. );
00911           add( 0., 0.,  64/81. ); break;
00912         }
00913         default:
00914           EXCEPTION( logic_error,"Invalid nb of gauss points for QUAD, variant 1: "
00915                      <<nbGauss);
00916         }
00917       }
00918       else if ( variant == 2 ) {
00919         if (geom == MED_QUAD4) setRefCoords( TQuad4a() );
00920         else                setRefCoords( TQuad8a() );
00921         switch ( nbGauss ) {
00922         case 4: {
00923           const double a = 1/sqrt(3.);
00924           add( -a,  a,  1 );
00925           add( -a, -a,  1 );
00926           add(  a, -a,  1 );
00927           add(  a,  a,  1 ); break;
00928         }
00929         case 9: {
00930           const double a = 0.774596669241483;
00931           add( -a,  a,  25/81. );
00932           add( -a, -a,  25/81. );
00933           add(  a, -a,  25/81. );
00934           add(  a,  a,  25/81. );
00935           add( -a, 0.,  40/81. );
00936           add( 0., -a,  40/81. );
00937           add(  a, 0.,  40/81. );
00938           add( 0.,  a,  40/81. );
00939           add( 0., 0.,  64/81. ); break;
00940         }
00941         default:
00942           EXCEPTION( logic_error,"Invalid nb of gauss points for QUAD, variant 1: "
00943                      <<nbGauss);
00944         }
00945       }
00946       else if ( variant == 3 ) {
00947         if (geom == MED_QUAD4) setRefCoords( TQuad4b() );
00948         else                setRefCoords( TQuad8b() );
00949         switch ( nbGauss ) {
00950         case 4: {
00951           const double a = 3/sqrt(3.);
00952           add( -a, -a,  1 );
00953           add( -a,  a,  1 );
00954           add(  a, -a,  1 );
00955           add(  a,  a,  1 ); break;
00956         }
00957         case 9: {
00958           const double a = sqrt(3/5.), c1 = 5/9., c2 = 8/9.;
00959           const double c12 = c1*c2, c22 = c2*c2, c1c2 = c1*c2;
00960           add( -a, -a,  c12  );
00961           add( -a, 0.,  c1c2 );
00962           add( -a,  a,  c12  );
00963           add( 0., -a,  c1c2 );
00964           add( 0., 0.,  c22  );
00965           add( 0.,  a,  c1c2 );
00966           add(  a, -a,  c12  );
00967           add(  a, 0.,  c1c2 );
00968           add(  a,  a,  c12  ); break;
00969         }
00970         default:
00971           EXCEPTION( logic_error,"Invalid nb of gauss points for QUAD, variant 3: "
00972                      <<nbGauss);
00973         }
00974       }
00975       break;
00976 
00977     case MED_TETRA4:
00978     case MED_TETRA10:
00979       if (geom == MED_TETRA4) setRefCoords( TTetra4a() );
00980       else                 setRefCoords( TTetra10a() );
00981       switch ( nbGauss ) {
00982       case 4: { // FPG4
00983         const double a = (5 - sqrt(5.))/20., b = (5 + 3*sqrt(5.))/20.;
00984         add(  a,  a,  a,  1/24. );
00985         add(  a,  a,  b,  1/24. );
00986         add(  a,  b,  a,  1/24. );
00987         add(  b,  a,  a,  1/24. ); break;
00988       }
00989       case 5: { // FPG5
00990         const double a = 0.25, b = 1/6., c = 0.5;
00991         add(  a,  a,  a, -2/15. );
00992         add(  b,  b,  b,  3/40. );
00993         add(  b,  b,  c,  3/40. );
00994         add(  b,  c,  b,  3/40. );
00995         add(  c,  b,  b,  3/40. ); break;
00996       }
00997       case 15: { // FPG15
00998         const double a = 0.25;
00999         const double b1 = (7 + sqrt(15.))/34., c1 = (13 + 3*sqrt(15.))/34., d = (5 - sqrt(15.))/20.;
01000         const double b2 = (7 - sqrt(15.))/34., c2 = (13 - 3*sqrt(15.))/34., e = (5 + sqrt(15.))/20.;
01001         const double P1 = (2665 - 14*sqrt(15.))/226800.;
01002         const double P2 = (2665 + 14*sqrt(15.))/226800.;
01003         add(  a,  a,  a,  8/405.);//_____
01004         add( b1, b1, b1,  P1    );
01005         add( b1, b1, c1,  P1    );
01006         add( b1, c1, b1,  P1    );
01007         add( c1, b1, b1,  P1    );//_____
01008         add( b2, b2, b2,  P2    );
01009         add( b2, b2, c2,  P2    );
01010         add( b2, c2, b2,  P2    );
01011         add( c2, b2, b2,  P2    );//_____
01012         add(  d,  d,  e,  5/567.);
01013         add(  d,  e,  d,  5/567.);
01014         add(  e,  d,  d,  5/567.);
01015         add(  d,  e,  e,  5/567.);
01016         add(  e,  d,  e,  5/567.);
01017         add(  e,  e,  d,  5/567.);
01018         break;
01019       }
01020       default:
01021         EXCEPTION( logic_error,"Invalid nb of gauss points for TETRA: "<<nbGauss);
01022       }
01023       break;
01024 
01025     case MED_PYRA5:
01026     case MED_PYRA13:
01027       if (geom == MED_PYRA5) setRefCoords( TPyra5a() );
01028       else                setRefCoords( TPyra13a() );
01029       switch ( nbGauss ) {
01030       case 5: { // FPG5
01031         const double h1 = 0.1531754163448146;
01032         const double h2 = 0.6372983346207416;
01033         add(  .5,  0.,  h1,  2/15. );
01034         add(  0.,  .5,  h1,  2/15. );
01035         add( -.5,  0.,  h1,  2/15. );
01036         add(  0., -.5,  h1,  2/15. );
01037         add(  0.,  0.,  h2,  2/15. ); break;
01038       }
01039       case 6: { // FPG6
01040         const double p1 = 0.1024890634400000 ;
01041         const double p2 = 0.1100000000000000 ;
01042         const double p3 = 0.1467104129066667 ;
01043         const double a  = 0.5702963741068025 ;
01044         const double h1 = 0.1666666666666666 ;
01045         const double h2 = 0.08063183038464675;
01046         const double h3 = 0.6098484849057127 ;
01047         add(  a, 0.,  h1,  p1 );
01048         add( 0.,  a,  h1,  p1 );
01049         add( -a, 0.,  h1,  p1 );
01050         add( 0., -a,  h1,  p1 );
01051         add( 0., 0.,  h2,  p2 );
01052         add( 0., 0.,  h3,  p3 ); break;
01053       }
01054       case 27: { // FPG27
01055         const double a1  = 0.788073483; 
01056         const double b6  = 0.499369002; 
01057         const double b1  = 0.848418011; 
01058         const double c8  = 0.478508449; 
01059         const double c1  = 0.652816472; 
01060         const double d12 = 0.032303742; 
01061         const double d1  = 1.106412899;
01062         double z = 1/2., fz = b1/2*(1 - z);
01063         add(  0.,  0.,   z,  a1 ); // 1
01064         add(  fz,  fz,   z,  b6 ); // 2
01065         add( -fz,  fz,   z,  b6 ); // 3
01066         add( -fz, -fz,   z,  b6 ); // 4
01067         add(  fz, -fz,   z,  b6 ); // 5
01068         z = (1 - b1)/2.;
01069         add(  0.,  0.,   z,  b6 ); // 6
01070         z = (1 + b1)/2.;
01071         add(  0.,  0.,   z,  b6 ); // 7
01072         z = (1 - c1)/2.; fz = c1*(1 - z);
01073         add(  fz,  0.,   z,  c8 ); // 8
01074         add(  0.,  fz,   z,  c8 ); // 9
01075         add( -fz,  0.,   z,  c8 ); // 10
01076         add(  0., -fz,   z,  c8 ); // 11
01077         z = (1 + c1)/2.; fz = c1*(1 - z);
01078         add(  fz,  0.,   z,  c8 ); // 12
01079         add(  0.,  fz,   z,  c8 ); // 13
01080         add( -fz,  0.,   z,  c8 ); // 14
01081         add(  0., -fz,   z,  c8 ); // 15
01082         z = (1 - d1)/2., fz = d1/2*(1 - z);
01083         add(  fz,  fz,   z,  d12); // 16
01084         add( -fz,  fz,   z,  d12); // 17
01085         add( -fz, -fz,   z,  d12); // 18
01086         add(  fz, -fz,   z,  d12); // 19
01087         z = 1/2.; fz = d1*(1 - z);
01088         add(  fz,  0.,   z,  d12); // 20
01089         add(  0.,  fz,   z,  d12); // 21
01090         add( -fz,  0.,   z,  d12); // 22
01091         add(  0., -fz,   z,  d12); // 23
01092         z = (1 + d1)/2., fz = d1/2*(1 - z);
01093         add(  fz,  fz,   z,  d12); // 24
01094         add( -fz,  fz,   z,  d12); // 25
01095         add( -fz, -fz,   z,  d12); // 26
01096         add(  fz, -fz,   z,  d12); // 27
01097         break;
01098       }
01099       default:
01100         EXCEPTION( logic_error,"Invalid nb of gauss points for PYRA: "<<nbGauss);
01101       }
01102       break;
01103     case MED_PENTA6:
01104     case MED_PENTA15:
01105       if (geom == MED_PENTA6) setRefCoords( TPenta6a() );
01106       else                 setRefCoords( TPenta15a() );
01107       switch ( nbGauss ) {
01108       case 6: { // FPG6
01109         const double a = sqrt(3.)/3.;
01110         add( -a, .5, .5,  1/6. );
01111         add( -a, 0., .5,  1/6. );
01112         add( -a, .5, 0.,  1/6. );
01113         add(  a, .5, .5,  1/6. );
01114         add(  a, 0., .5,  1/6. );
01115         add(  a, .5, 0.,  1/6. ); break;
01116       }
01117       case 8: { // FPG8
01118         const double a = 0.577350269189626;
01119         add( -a, 1/3., 1/3., -27/96. );
01120         add( -a,  0.6,  0.2,  25/96. );
01121         add( -a,  0.2,  0.6,  25/96. );
01122         add( -a,  0.2,  0.2,  25/96. );
01123         add( +a, 1/3., 1/3., -27/96. );
01124         add( +a,  0.6,  0.2,  25/96. );
01125         add( +a,  0.2,  0.6,  25/96. );
01126         add( +a,  0.2,  0.2,  25/96. ); break;
01127       }
01128       case 21: { // FPG21
01129         const double d = sqrt(3/5.), c1 = 5/9., c2 = 8/9.; // d <=> alfa
01130         const double a = (6 + sqrt(15.))/21.;
01131         const double b = (6 - sqrt(15.))/21.;
01132         const double P1 = (155 + sqrt(15.))/2400.;
01133         const double P2 = (155 - sqrt(15.))/2400.;  //___
01134         add( -d,  1/3.,  1/3., c1*9/80. );//___
01135         add( -d,     a,     a, c1*P1    );
01136         add( -d, 1-2*a,     a, c1*P1    );
01137         add( -d,     a, 1-2*a, c1*P1    );//___
01138         add( -d,     b,     b, c1*P2    );
01139         add( -d, 1-2*b,     b, c1*P2    );
01140         add( -d,     b, 1-2*b, c1*P2    );//___
01141         add( 0.,  1/3.,  1/3., c2*9/80. );//___
01142         add( 0.,     a,     a, c2*P1    );
01143         add( 0., 1-2*a,     a, c2*P1    );
01144         add( 0.,     a, 1-2*a, c2*P1    );//___
01145         add( 0.,     b,     b, c2*P2    );
01146         add( 0., 1-2*b,     b, c2*P2    );
01147         add( 0.,     b, 1-2*b, c2*P2    );//___
01148         add(  d,  1/3.,  1/3., c1*9/80. );//___
01149         add(  d,     a,     a, c1*P1    );
01150         add(  d, 1-2*a,     a, c1*P1    );
01151         add(  d,     a, 1-2*a, c1*P1    );//___
01152         add(  d,     b,     b, c1*P2    );
01153         add(  d, 1-2*b,     b, c1*P2    );
01154         add(  d,     b, 1-2*b, c1*P2    );//___
01155         break;
01156       }
01157       default:
01158         EXCEPTION( logic_error,"Invalid nb of gauss points for PENTA: " <<nbGauss);
01159       }
01160       break;
01161 
01162     case MED_HEXA8:
01163     case MED_HEXA20:
01164       if (geom == MED_HEXA8) setRefCoords( THexa8a() );
01165       else                setRefCoords( THexa20a() );
01166       switch ( nbGauss ) {
01167       case 8: { // FPG8
01168         const double a = sqrt(3.)/3.;
01169         add( -a, -a, -a,  1. );
01170         add( -a, -a,  a,  1. );
01171         add( -a,  a, -a,  1. );
01172         add( -a,  a,  a,  1. );
01173         add(  a, -a, -a,  1. );
01174         add(  a, -a,  a,  1. );
01175         add(  a,  a, -a,  1. );
01176         add(  a,  a,  a,  1. ); break;
01177       }
01178       case 27: { // FPG27
01179         const double a = sqrt(3/5.), c1 = 5/9., c2 = 8/9.;
01180         const double c12 = c1*c1, c13 = c1*c1*c1;
01181         const double c22 = c2*c2, c23 = c2*c2*c2;
01182         add( -a, -a, -a,   c13  ); // 1
01183         add( -a, -a, 0., c12*c2 ); // 2
01184         add( -a, -a,  a,   c13  ); // 3
01185         add( -a, 0., -a, c12*c2 ); // 4
01186         add( -a, 0., 0., c1*c22 ); // 5
01187         add( -a, 0.,  a, c12*c2 ); // 6
01188         add( -a,  a, -a,   c13  ); // 7
01189         add( -a,  a, 0., c12*c2 ); // 8
01190         add( -a,  a,  a,   c13  ); // 9
01191         add( 0., -a, -a, c12*c2 ); // 10
01192         add( 0., -a, 0., c1*c22 ); // 11
01193         add( 0., -a,  a, c12*c2 ); // 12
01194         add( 0., 0., -a, c1*c22 ); // 13
01195         add( 0., 0., 0.,   c23  ); // 14
01196         add( 0., 0.,  a, c1*c22 ); // 15
01197         add( 0.,  a, -a, c12*c2 ); // 16
01198         add( 0.,  a, 0., c1*c22 ); // 17
01199         add( 0.,  a,  a, c12*c2 ); // 18
01200         add(  a, -a, -a,   c13  ); // 19
01201         add(  a, -a, 0., c12*c2 ); // 20
01202         add(  a, -a,  a,   c13  ); // 21
01203         add(  a, 0., -a, c12*c2 ); // 22
01204         add(  a, 0., 0., c1*c22 ); // 23
01205         add(  a, 0.,  a, c12*c2 ); // 24
01206         add(  a,  a, -a,   c13  ); // 25
01207         add(  a,  a, 0., c12*c2 ); // 26
01208         add(  a,  a,  a,   c13  ); // 27
01209         break;
01210       }
01211       default:
01212         EXCEPTION( logic_error,"Invalid nb of gauss points for PENTA: " <<nbGauss);
01213       }
01214       break;
01215 
01216     default:
01217       EXCEPTION( logic_error,"unexpected EGeometrieElement: "<< geom);
01218     }
01219 
01220     if ( myWeights.capacity() != myWeights.size() )
01221       EXCEPTION( logic_error,"Not all gauss points defined");
01222   }
01223 }
01224 
01225 //=======================================================================
01229 //=======================================================================
01230 
01231 namespace MEDMEM {
01232 
01233   GAUSS_LOCALIZATION_*
01234   GAUSS_LOCALIZATION_::makeDefaultLocalization(const string &     locName,
01235                                                medGeometryElement typeGeo,
01236                                                int                nGauss) throw (MEDEXCEPTION)
01237   {
01238     TGaussDef gaussDef( typeGeo, nGauss, 1 );
01239 
01240     GAUSS_LOCALIZATION_ * gsloc = 
01241       new GAUSS_LOCALIZATION<FullInterlace> ( locName.c_str(),
01242                                               typeGeo,
01243                                               nGauss,
01244                                               &gaussDef.myRefCoords[0],
01245                                               &gaussDef.myCoords[0],
01246                                               &gaussDef.myWeights[0] );
01247     return gsloc;
01248   }
01249 }