src/MEDMEM/MEDMEM_GaussLocalization.cxx

   1 // Copyright (C) 2007-2012  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.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19
  20 // File      : MEDMEM_GaussLocalization.cxx
  21 // Created   : Thu Dec 20 12:26:33 2007
  22 // Author    : Edward AGAPOV (eap)
  23 //
  24 #include "MEDMEM_GaussLocalization.hxx"
  25
  26 #include <vector>
  27 #include <math.h>
  28
  29 #define EXCEPTION(type,msg) \
  30   MEDEXCEPTION(LOCALIZED(MEDMEM::STRING(#type)<< msg))
  31
  32 namespace // matters used by GAUSS_LOCALIZATION_::makeDefaultLocalization()
  33 {
  34   typedef std::vector<double> TDoubleVector;
  35   typedef double*             TCoordSlice;
  36   typedef int                 TInt;
  37   //---------------------------------------------------------------
  38   //! Shape function definitions
  39   //---------------------------------------------------------------
  40   struct TShapeFun
  41   {
  42     TInt myDim;
  43     TInt myNbRef;
  44     TDoubleVector myRefCoord;
  45
  46     TShapeFun(TInt theDim = 0, TInt theNbRef = 0)
  47       :myDim(theDim),myNbRef(theNbRef),myRefCoord(theNbRef*theDim) {}
  48
  49     TInt GetNbRef() const { return myNbRef; }
  50
  51     TCoordSlice GetCoord(TInt theRefId) { return &myRefCoord[0] + theRefId * myDim; }
  52   };
  53   //---------------------------------------------------------------
  54   /*!
  55    * \brief Description of family of integration points
  56    */
  57   //---------------------------------------------------------------
  58   struct TGaussDef
  59   {
  60     int           myType;      //!< element geometry (EGeometrieElement or med_geometrie_element)
  61     TDoubleVector myRefCoords; //!< description of reference points
  62     TDoubleVector myCoords;    //!< coordinates of Gauss points
  63     TDoubleVector myWeights;   //!< weights, len(weights)==<nb of gauss points>
  64
  65     /*!
  66      * \brief Creates definition of gauss points family
  67      *  \param geomType - element geometry (EGeometrieElement or med_geometrie_element)
  68      *  \param nbPoints - nb gauss point
  69      *  \param variant - [1-3] to choose the variant of definition
  70      *
  71      * Throws in case of invalid parameters
  72      * variant == 1 refers to "Fonctions de forme et points d'integration
  73      *              des elements finis" v7.4 by J. PELLET, X. DESROCHES, 15/09/05
  74      * variant == 2 refers to the same doc v6.4 by J.P. LEFEBVRE, X. DESROCHES, 03/07/03
  75      * variant == 3 refers to the same doc v6.4, second variant for 2D elements
  76      */
  77     TGaussDef(const int geomType, const int nbPoints, const int variant=1);
  78
  79     int dim() const { return myType/100; }
  80     int nbPoints() const { return myWeights.capacity(); }
  81
  82   private:
  83     void add(const double x, const double weight);
  84     void add(const double x, const double y, const double weight);
  85     void add(const double x, const double y, const double z, const double weight);
  86     void setRefCoords(const TShapeFun& aShapeFun) { myRefCoords = aShapeFun.myRefCoord; }
  87   };
  88   struct TSeg2a: TShapeFun {
  89     TSeg2a();
  90   };
  91   struct TSeg3a: TShapeFun {
  92     TSeg3a();
  93   };
  94   struct TTria3a: TShapeFun {
  95     TTria3a();
  96   };
  97   struct TTria6a: TShapeFun {
  98     TTria6a();
  99   };
 100   struct TTria3b: TShapeFun {
 101     TTria3b();
 102   };
 103   struct TTria6b: TShapeFun {
 104     TTria6b();
 105   };
 106   struct TQuad4a: TShapeFun {
 107     TQuad4a();
 108   };
 109   struct TQuad8a: TShapeFun {
 110     TQuad8a();
 111   };
 112   struct TQuad4b: TShapeFun {
 113     TQuad4b();
 114   };
 115   struct TQuad8b: TShapeFun {
 116     TQuad8b();
 117   };
 118   struct TTetra4a: TShapeFun {
 119     TTetra4a();
 120   };
 121   struct TTetra10a: TShapeFun {
 122     TTetra10a();
 123   };
 124   struct TTetra4b: TShapeFun {
 125     TTetra4b();
 126   };
 127   struct TTetra10b: TShapeFun {
 128     TTetra10b();
 129   };
 130   struct THexa8a: TShapeFun {
 131     THexa8a();
 132   };
 133   struct THexa20a: TShapeFun {
 134     THexa20a(TInt theDim = 3, TInt theNbRef = 20);
 135   };
 136   struct THexa27a: THexa20a {
 137     THexa27a();
 138   };
 139   struct THexa8b: TShapeFun {
 140     THexa8b();
 141   };
 142   struct THexa20b: TShapeFun {
 143     THexa20b(TInt theDim = 3, TInt theNbRef = 20);
 144   };
 145   struct TPenta6a: TShapeFun {
 146     TPenta6a();
 147   };
 148   struct TPenta6b: TShapeFun {
 149     TPenta6b();
 150   };
 151   struct TPenta15a: TShapeFun {
 152     TPenta15a();
 153   };
 154   struct TPenta15b: TShapeFun {
 155     TPenta15b();
 156   };
 157   struct TPyra5a: TShapeFun {
 158     TPyra5a();
 159   };
 160   struct TPyra5b: TShapeFun {
 161     TPyra5b();
 162   };
 163   struct TPyra13a: TShapeFun {
 164     TPyra13a();
 165   };
 166   struct TPyra13b: TShapeFun {
 167     TPyra13b();
 168   };
 169
 170   void TGaussDef::add(const double x, const double weight)
 171   {
 172     if ( dim() != 1 )
 173       EXCEPTION( logic_error,"dim() != 1");
 174     if ( myWeights.capacity() == myWeights.size() )
 175       EXCEPTION( logic_error,"Extra gauss point");
 176     myCoords.push_back( x );
 177     myWeights.push_back( weight );
 178   }
 179   void TGaussDef::add(const double x, const double y, const double weight)
 180   {
 181     if ( dim() != 2 )
 182       EXCEPTION( logic_error,"dim() != 2");
 183     if ( myWeights.capacity() == myWeights.size() )
 184       EXCEPTION( logic_error,"Extra gauss point");
 185     myCoords.push_back( x );
 186     myCoords.push_back( y );
 187     myWeights.push_back( weight );
 188   }
 189   void TGaussDef::add(const double x, const double y, const double z, const double weight)
 190   {
 191     if ( dim() != 3 )
 192       EXCEPTION( logic_error,"dim() != 3");
 193     if ( myWeights.capacity() == myWeights.size() )
 194       EXCEPTION( logic_error,"Extra gauss point");
 195     myCoords.push_back( x );
 196     myCoords.push_back( y );
 197     myCoords.push_back( z );
 198     myWeights.push_back( weight );
 199   }
 200
 201   //---------------------------------------------------------------
 202   TSeg2a::TSeg2a():TShapeFun(1,2)
 203   {
 204     TInt aNbRef = GetNbRef();
 205     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 206       TCoordSlice aCoord = GetCoord(aRefId);
 207       switch(aRefId){
 208       case  0: aCoord[0] = -1.0; break;
 209       case  1: aCoord[0] =  1.0; break;
 210       }
 211     }
 212   }
 213   //---------------------------------------------------------------
 214   TSeg3a::TSeg3a():TShapeFun(1,3)
 215   {
 216     TInt aNbRef = GetNbRef();
 217     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 218       TCoordSlice aCoord = GetCoord(aRefId);
 219       switch(aRefId){
 220       case  0: aCoord[0] = -1.0; break;
 221       case  1: aCoord[0] =  1.0; break;
 222       case  2: aCoord[0] =  0.0; break;
 223       }
 224     }
 225   }
 226   //---------------------------------------------------------------
 227   TTria3a::TTria3a():
 228     TShapeFun(2,3)
 229   {
 230     TInt aNbRef = GetNbRef();
 231     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 232       TCoordSlice aCoord = GetCoord(aRefId);
 233       switch(aRefId){
 234       case  0: aCoord[0] = -1.0;  aCoord[1] =  1.0; break;
 235       case  1: aCoord[0] = -1.0;  aCoord[1] = -1.0; break;
 236       case  2: aCoord[0] =  1.0;  aCoord[1] = -1.0; break;
 237       }
 238     }
 239   }
 240   //---------------------------------------------------------------
 241   TTria6a::TTria6a():TShapeFun(2,6)
 242   {
 243     TInt aNbRef = GetNbRef();
 244     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 245       TCoordSlice aCoord = GetCoord(aRefId);
 246       switch(aRefId){
 247       case  0: aCoord[0] = -1.0;  aCoord[1] =  1.0; break;
 248       case  1: aCoord[0] = -1.0;  aCoord[1] = -1.0; break;
 249       case  2: aCoord[0] =  1.0;  aCoord[1] = -1.0; break;
 250
 251       case  3: aCoord[0] = -1.0;  aCoord[1] =  1.0; break;
 252       case  4: aCoord[0] =  0.0;  aCoord[1] = -1.0; break;
 253       case  5: aCoord[0] =  0.0;  aCoord[1] =  0.0; break;
 254       }
 255     }
 256   }
 257   //---------------------------------------------------------------
 258   TTria3b::TTria3b():
 259     TShapeFun(2,3)
 260   {
 261     TInt aNbRef = GetNbRef();
 262     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 263       TCoordSlice aCoord = GetCoord(aRefId);
 264       switch(aRefId){
 265       case  0: aCoord[0] =  0.0;  aCoord[1] =  0.0; break;
 266       case  1: aCoord[0] =  1.0;  aCoord[1] =  0.0; break;
 267       case  2: aCoord[0] =  0.0;  aCoord[1] =  1.0; break;
 268       }
 269     }
 270   }
 271   //---------------------------------------------------------------
 272   TTria6b::TTria6b():
 273     TShapeFun(2,6)
 274   {
 275     TInt aNbRef = GetNbRef();
 276     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 277       TCoordSlice aCoord = GetCoord(aRefId);
 278       switch(aRefId){
 279       case  0: aCoord[0] =  0.0;  aCoord[1] =  0.0; break;
 280       case  1: aCoord[0] =  1.0;  aCoord[1] =  0.0; break;
 281       case  2: aCoord[0] =  0.0;  aCoord[1] =  1.0; break;
 282
 283       case  3: aCoord[0] =  0.5;  aCoord[1] =  0.0; break;
 284       case  4: aCoord[0] =  0.5;  aCoord[1] =  0.5; break;
 285       case  5: aCoord[0] =  0.0;  aCoord[1] =  0.5; break;
 286       }
 287     }
 288   }
 289   //---------------------------------------------------------------
 290   TQuad4a::TQuad4a():
 291     TShapeFun(2,4)
 292   {
 293     TInt aNbRef = GetNbRef();
 294     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 295       TCoordSlice aCoord = GetCoord(aRefId);
 296       switch(aRefId){
 297       case  0: aCoord[0] = -1.0;  aCoord[1] =  1.0; break;
 298       case  1: aCoord[0] = -1.0;  aCoord[1] = -1.0; break;
 299       case  2: aCoord[0] =  1.0;  aCoord[1] = -1.0; break;
 300       case  3: aCoord[0] =  1.0;  aCoord[1] =  1.0; break;
 301       }
 302     }
 303   }
 304   //---------------------------------------------------------------
 305   TQuad8a::TQuad8a():
 306     TShapeFun(2,8)
 307   {
 308     TInt aNbRef = GetNbRef();
 309     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 310       TCoordSlice aCoord = GetCoord(aRefId);
 311       switch(aRefId){
 312       case  0: aCoord[0] = -1.0;  aCoord[1] =  1.0; break;
 313       case  1: aCoord[0] = -1.0;  aCoord[1] = -1.0; break;
 314       case  2: aCoord[0] =  1.0;  aCoord[1] = -1.0; break;
 315       case  3: aCoord[0] =  1.0;  aCoord[1] =  1.0; break;
 316
 317       case  4: aCoord[0] = -1.0;  aCoord[1] =  0.0; break;
 318       case  5: aCoord[0] =  0.0;  aCoord[1] = -1.0; break;
 319       case  6: aCoord[0] =  1.0;  aCoord[1] =  0.0; break;
 320       case  7: aCoord[0] =  0.0;  aCoord[1] =  1.0; break;
 321       }
 322     }
 323   }
 324   //---------------------------------------------------------------
 325   TQuad4b::TQuad4b():
 326     TShapeFun(2,4)
 327   {
 328     TInt aNbRef = GetNbRef();
 329     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 330       TCoordSlice aCoord = GetCoord(aRefId);
 331       switch(aRefId){
 332       case  0: aCoord[0] = -1.0;  aCoord[1] = -1.0; break;
 333       case  1: aCoord[0] =  1.0;  aCoord[1] = -1.0; break;
 334       case  2: aCoord[0] =  1.0;  aCoord[1] =  1.0; break;
 335       case  3: aCoord[0] = -1.0;  aCoord[1] =  1.0; break;
 336       }
 337     }
 338   }
 339   //---------------------------------------------------------------
 340   TQuad8b::TQuad8b():
 341     TShapeFun(2,8)
 342   {
 343     TInt aNbRef = GetNbRef();
 344     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 345       TCoordSlice aCoord = GetCoord(aRefId);
 346       switch(aRefId){
 347       case  0: aCoord[0] = -1.0;  aCoord[1] = -1.0; break;
 348       case  1: aCoord[0] =  1.0;  aCoord[1] = -1.0; break;
 349       case  2: aCoord[0] =  1.0;  aCoord[1] =  1.0; break;
 350       case  3: aCoord[0] = -1.0;  aCoord[1] =  1.0; break;
 351
 352       case  4: aCoord[0] =  0.0;  aCoord[1] = -1.0; break;
 353       case  5: aCoord[0] =  1.0;  aCoord[1] =  0.0; break;
 354       case  6: aCoord[0] =  0.0;  aCoord[1] =  1.0; break;
 355       case  7: aCoord[0] = -1.0;  aCoord[1] =  0.0; break;
 356       }
 357     }
 358   }
 359   //---------------------------------------------------------------
 360   TTetra4a::TTetra4a():
 361     TShapeFun(3,4)
 362   {
 363     TInt aNbRef = GetNbRef();
 364     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 365       TCoordSlice aCoord = GetCoord(aRefId);
 366       switch(aRefId){
 367       case  0: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
 368       case  1: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
 369       case  2: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 370       case  3: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 371       }
 372     }
 373   }
 374   //---------------------------------------------------------------
 375   TTetra10a::TTetra10a():
 376     TShapeFun(3,10)
 377   {
 378     TInt aNbRef = GetNbRef();
 379     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 380       TCoordSlice aCoord = GetCoord(aRefId);
 381       switch(aRefId){
 382       case  0: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
 383       case  1: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
 384       case  2: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 385       case  3: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 386
 387       case  4: aCoord[0] =  0.0;  aCoord[1] =  0.5;  aCoord[2] =  0.5; break;
 388       case  5: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
 389       case  6: aCoord[0] =  0.0;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
 390
 391       case  7: aCoord[0] =  0.5;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
 392       case  8: aCoord[0] =  0.5;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
 393       case  9: aCoord[0] =  0.5;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 394       }
 395     }
 396   }
 397   //---------------------------------------------------------------
 398   TTetra4b::TTetra4b():
 399     TShapeFun(3,4)
 400   {
 401     TInt aNbRef = GetNbRef();
 402     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 403       TCoordSlice aCoord = GetCoord(aRefId);
 404       switch(aRefId){
 405       case  0: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
 406       case  2: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
 407       case  1: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 408       case  3: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 409       }
 410     }
 411   }
 412   //---------------------------------------------------------------
 413   TTetra10b::TTetra10b():
 414     TShapeFun(3,10)
 415   {
 416     TInt aNbRef = GetNbRef();
 417     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 418       TCoordSlice aCoord = GetCoord(aRefId);
 419       switch(aRefId){
 420       case  0: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
 421       case  2: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
 422       case  1: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 423       case  3: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 424
 425       case  6: aCoord[0] =  0.0;  aCoord[1] =  0.5;  aCoord[2] =  0.5; break;
 426       case  5: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
 427       case  4: aCoord[0] =  0.0;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
 428
 429       case  7: aCoord[0] =  0.5;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
 430       case  9: aCoord[0] =  0.5;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
 431       case  8: aCoord[0] =  0.5;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 432       }
 433     }
 434   }
 435   //---------------------------------------------------------------
 436   THexa8a::THexa8a():
 437     TShapeFun(3,8)
 438   {
 439     TInt aNbRef = GetNbRef();
 440     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 441       TCoordSlice aCoord = GetCoord(aRefId);
 442       switch(aRefId){
 443       case  0: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
 444       case  1: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
 445       case  2: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
 446       case  3: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
 447       case  4: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
 448       case  5: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
 449       case  6: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
 450       case  7: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
 451       }
 452     }
 453   }
 454   //---------------------------------------------------------------
 455   THexa20a::THexa20a(TInt theDim, TInt theNbRef):
 456     TShapeFun(theDim,theNbRef)
 457   {
 458     TInt aNbRef = myRefCoord.size();
 459     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 460       TCoordSlice aCoord = GetCoord(aRefId);
 461       switch(aRefId){
 462       case  0: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
 463       case  1: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
 464       case  2: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
 465       case  3: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
 466       case  4: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
 467       case  5: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
 468       case  6: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
 469       case  7: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
 470
 471       case  8: aCoord[0] =  0.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
 472       case  9: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] = -1.0; break;
 473       case 10: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
 474       case 11: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] = -1.0; break;
 475       case 12: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] =  0.0; break;
 476       case 13: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] =  0.0; break;
 477       case 14: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
 478       case 15: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
 479       case 16: aCoord[0] =  0.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
 480       case 17: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
 481       case 18: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
 482       case 19: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
 483       }
 484     }
 485   }
 486   //---------------------------------------------------------------
 487   THexa27a::THexa27a():
 488     THexa20a(3,27)
 489   {
 490     TInt aNbRef = myRefCoord.size();
 491     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 492       TCoordSlice aCoord = GetCoord(aRefId);
 493       switch(aRefId){
 494       case 20: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] = -1.0; break;
 495       case 21: aCoord[0] =  0.0;  aCoord[1] = -1.0;  aCoord[2] =  0.0; break;
 496       case 22: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 497       case 23: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
 498       case 24: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 499       case 25: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
 500       case 26: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 501       }
 502     }
 503   }
 504   //---------------------------------------------------------------
 505   THexa8b::THexa8b():
 506     TShapeFun(3,8)
 507   {
 508     TInt aNbRef = GetNbRef();
 509     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 510       TCoordSlice aCoord = GetCoord(aRefId);
 511       switch(aRefId){
 512       case  0: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
 513       case  3: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
 514       case  2: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
 515       case  1: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
 516       case  4: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
 517       case  7: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
 518       case  6: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
 519       case  5: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
 520       }
 521     }
 522   }
 523   //---------------------------------------------------------------
 524   THexa20b::THexa20b(TInt theDim, TInt theNbRef):
 525     TShapeFun(theDim,theNbRef)
 526   {
 527     TInt aNbRef = myRefCoord.size();
 528     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 529       TCoordSlice aCoord = GetCoord(aRefId);
 530       switch(aRefId){
 531       case  0: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
 532       case  3: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
 533       case  2: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
 534       case  1: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
 535       case  4: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
 536       case  7: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
 537       case  6: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
 538       case  5: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
 539
 540       case 11: aCoord[0] =  0.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
 541       case 10: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] = -1.0; break;
 542       case  9: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
 543       case  8: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] = -1.0; break;
 544       case 16: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] =  0.0; break;
 545       case 19: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] =  0.0; break;
 546       case 18: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
 547       case 17: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
 548       case 15: aCoord[0] =  0.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
 549       case 14: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
 550       case 13: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
 551       case 12: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
 552       }
 553     }
 554   }
 555   //---------------------------------------------------------------
 556   TPenta6a::TPenta6a():
 557     TShapeFun(3,6)
 558   {
 559     TInt aNbRef = myRefCoord.size();
 560     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 561       TCoordSlice aCoord = GetCoord(aRefId);
 562       switch(aRefId){
 563       case  0: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
 564       case  1: aCoord[0] = -1.0;  aCoord[1] = -0.0;  aCoord[2] =  1.0; break;
 565       case  2: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 566       case  3: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
 567       case  4: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
 568       case  5: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 569       }
 570     }
 571   }
 572   //---------------------------------------------------------------
 573   TPenta6b::TPenta6b():
 574     TShapeFun(3,6)
 575   {
 576     TInt aNbRef = myRefCoord.size();
 577     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 578       TCoordSlice aCoord = GetCoord(aRefId);
 579       switch(aRefId){
 580       case  0: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
 581       case  2: aCoord[0] = -1.0;  aCoord[1] = -0.0;  aCoord[2] =  1.0; break;
 582       case  1: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 583       case  3: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
 584       case  5: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
 585       case  4: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 586       }
 587     }
 588   }
 589   //---------------------------------------------------------------
 590   TPenta15a::TPenta15a():
 591     TShapeFun(3,15)
 592   {
 593     TInt aNbRef = myRefCoord.size();
 594     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 595       TCoordSlice aCoord = GetCoord(aRefId);
 596       switch(aRefId){
 597       case  0: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
 598       case  1: aCoord[0] = -1.0;  aCoord[1] = -0.0;  aCoord[2] =  1.0; break;
 599       case  2: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 600       case  3: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
 601       case  4: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
 602       case  5: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 603
 604       case  6: aCoord[0] = -1.0;  aCoord[1] =  0.5;  aCoord[2] =  0.5; break;
 605       case  7: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
 606       case  8: aCoord[0] = -1.0;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
 607       case  9: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
 608       case 10: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
 609       case 11: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 610       case 12: aCoord[0] =  1.0;  aCoord[1] =  0.5;  aCoord[2] =  0.5; break;
 611       case 13: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
 612       case 14: aCoord[0] =  1.0;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
 613       }
 614     }
 615   }
 616   //---------------------------------------------------------------
 617   TPenta15b::TPenta15b():
 618     TShapeFun(3,15)
 619   {
 620     TInt aNbRef = myRefCoord.size();
 621     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 622       TCoordSlice aCoord = GetCoord(aRefId);
 623       switch(aRefId){
 624       case  0: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
 625       case  2: aCoord[0] = -1.0;  aCoord[1] = -0.0;  aCoord[2] =  1.0; break;
 626       case  1: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 627       case  3: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
 628       case  5: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
 629       case  4: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 630
 631       case  8: aCoord[0] = -1.0;  aCoord[1] =  0.5;  aCoord[2] =  0.5; break;
 632       case  7: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
 633       case  6: aCoord[0] = -1.0;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
 634       case 12: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
 635       case 14: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
 636       case 13: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 637       case 11: aCoord[0] =  1.0;  aCoord[1] =  0.5;  aCoord[2] =  0.5; break;
 638       case 10: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
 639       case  9: aCoord[0] =  1.0;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
 640       }
 641     }
 642   }
 643   //---------------------------------------------------------------
 644   TPyra5a::TPyra5a():
 645     TShapeFun(3,5)
 646   {
 647     TInt aNbRef = myRefCoord.size();
 648     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 649       TCoordSlice aCoord = GetCoord(aRefId);
 650       switch(aRefId){
 651       case  0: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 652       case  1: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
 653       case  2: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 654       case  3: aCoord[0] =  0.0;  aCoord[1] = -1.0;  aCoord[2] =  0.0; break;
 655       case  4: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
 656       }
 657     }
 658   }
 659   //---------------------------------------------------------------
 660   TPyra5b::TPyra5b():
 661     TShapeFun(3,5)
 662   {
 663     TInt aNbRef = myRefCoord.size();
 664     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 665       TCoordSlice aCoord = GetCoord(aRefId);
 666       switch(aRefId){
 667       case  0: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 668       case  3: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
 669       case  2: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 670       case  1: aCoord[0] =  0.0;  aCoord[1] = -1.0;  aCoord[2] =  0.0; break;
 671       case  4: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
 672       }
 673     }
 674   }
 675   //---------------------------------------------------------------
 676   TPyra13a::TPyra13a():
 677     TShapeFun(3,13)
 678   {
 679     TInt aNbRef = myRefCoord.size();
 680     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 681       TCoordSlice aCoord = GetCoord(aRefId);
 682       switch(aRefId){
 683       case  0: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 684       case  1: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
 685       case  2: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 686       case  3: aCoord[0] =  0.0;  aCoord[1] = -1.0;  aCoord[2] =  0.0; break;
 687       case  4: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
 688
 689       case  5: aCoord[0] =  0.5;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
 690       case  6: aCoord[0] = -0.5;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
 691       case  7: aCoord[0] = -0.5;  aCoord[1] = -0.5;  aCoord[2] =  0.0; break;
 692       case  8: aCoord[0] =  0.5;  aCoord[1] = -0.5;  aCoord[2] =  0.0; break;
 693       case  9: aCoord[0] =  0.5;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
 694       case 10: aCoord[0] =  0.0;  aCoord[1] =  0.5;  aCoord[2] =  0.5; break;
 695       case 11: aCoord[0] = -0.5;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
 696       case 12: aCoord[0] =  0.0;  aCoord[1] = -0.5;  aCoord[2] =  0.5; break;
 697       }
 698     }
 699   }
 700   //---------------------------------------------------------------
 701   TPyra13b::TPyra13b():
 702     TShapeFun(3,13)
 703   {
 704     TInt aNbRef = myRefCoord.size();
 705     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 706       TCoordSlice aCoord = GetCoord(aRefId);
 707       switch(aRefId){
 708       case  0: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 709       case  3: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
 710       case  2: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 711       case  1: aCoord[0] =  0.0;  aCoord[1] = -1.0;  aCoord[2] =  0.0; break;
 712       case  4: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
 713
 714       case  8: aCoord[0] =  0.5;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
 715       case  7: aCoord[0] = -0.5;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
 716       case  6: aCoord[0] = -0.5;  aCoord[1] = -0.5;  aCoord[2] =  0.0; break;
 717       case  5: aCoord[0] =  0.5;  aCoord[1] = -0.5;  aCoord[2] =  0.0; break;
 718       case  9: aCoord[0] =  0.5;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
 719       case 12: aCoord[0] =  0.0;  aCoord[1] =  0.5;  aCoord[2] =  0.5; break;
 720       case 11: aCoord[0] = -0.5;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
 721       case 10: aCoord[0] =  0.0;  aCoord[1] = -0.5;  aCoord[2] =  0.5; break;
 722       }
 723     }
 724   }
 725   /*!
 726    * \brief Fill definition of gauss points family
 727    */
 728
 729   TGaussDef::TGaussDef(const int geom, const int nbGauss, const int variant)
 730   {
 731     myType = geom;
 732     myCoords .reserve( nbGauss * dim() );
 733     myWeights.reserve( nbGauss );
 734
 735     switch ( geom ) {
 736
 737     case MED_SEG2:
 738     case MED_SEG3:
 739       if (geom == MED_SEG2) setRefCoords( TSeg2a() );
 740       else               setRefCoords( TSeg3a() );
 741       switch ( nbGauss ) {
 742       case 1: {
 743         add( 0.0, 2.0 ); break;
 744       }
 745       case 2: {
 746         const double a = 0.577350269189626;
 747         add(  a,  1.0 );
 748         add( -a,  1.0 ); break;
 749       }
 750       case 3: {
 751         const double a = 0.774596669241;
 752         const double P1 = 1./1.8;
 753         const double P2 = 1./1.125;
 754         add( -a,  P1 );
 755         add(  0,  P2 );
 756         add(  a,  P1 ); break;
 757       }
 758       case 4: {
 759         const double a  = 0.339981043584856, b  = 0.861136311594053;
 760         const double P1 = 0.652145154862546, P2 = 0.347854845137454 ;
 761         add(  a,  P1 );
 762         add( -a,  P1 );
 763         add(  b,  P2 );
 764         add( -b,  P2 ); break;
 765       }
 766       default:
 767         EXCEPTION( logic_error,"Invalid nb of gauss points for SEG"<<nbGauss);
 768       }
 769       break;
 770
 771     case MED_TRIA3:
 772     case MED_TRIA6:
 773       if ( variant == 1 ) {
 774         if (geom == MED_TRIA3) setRefCoords( TTria3b() );
 775         else                setRefCoords( TTria6b() );
 776         switch ( nbGauss ) {
 777         case 1: { // FPG1
 778           add( 1/3., 1/3., 1/2. ); break;
 779         }
 780         case 3: { // FPG3
 781           // what about COT3 ???
 782           add( 1/6., 1/6., 1/6. );
 783           add( 2/3., 1/6., 1/6. );
 784           add( 1/6., 2/3., 1/6. ); break;
 785         }
 786         case 4: { // FPG4
 787           add( 1/5., 1/5.,  25/(24*4.) );
 788           add( 3/5., 1/5.,  25/(24*4.) );
 789           add( 1/5., 3/5.,  25/(24*4.) );
 790           add( 1/3., 1/3., -27/(24*4.) ); break;
 791         }
 792         case 6: { // FPG6
 793           const double P1 = 0.11169079483905, P2 = 0.0549758718227661;
 794           const double a  = 0.445948490915965, b = 0.091576213509771;
 795           add(     b,     b, P2 );
 796           add( 1-2*b,     b, P2 );
 797           add(     b, 1-2*b, P2 );
 798           add(     a, 1-2*a, P1 );
 799           add(     a,     a, P1 );
 800           add( 1-2*a,     a, P1 ); break;
 801         }
 802         case 7: { // FPG7
 803           const double A  = 0.470142064105115;
 804           const double B  = 0.101286507323456;
 805           const double P1 = 0.066197076394253;
 806           const double P2 = 0.062969590272413;
 807           add(  1/3.,  1/3., 9/80. );
 808           add(     A,     A, P1 );
 809           add( 1-2*A,     A, P1 );
 810           add(     A, 1-2*A, P1 );
 811           add(     B,     B, P2 );
 812           add( 1-2*B,     B, P2 );
 813           add(     B, 1-2*B, P2 ); break;
 814         }
 815         case 12: { // FPG12
 816           const double A  = 0.063089014491502;
 817           const double B  = 0.249286745170910;
 818           const double C  = 0.310352451033785;
 819           const double D  = 0.053145049844816;
 820           const double P1 = 0.025422453185103;
 821           const double P2 = 0.058393137863189;
 822           const double P3 = 0.041425537809187;
 823           add(     A,     A, P1 );
 824           add( 1-2*A,     A, P1 );
 825           add(     A, 1-2*A, P1 );
 826           add(     B,     B, P2 );
 827           add( 1-2*B,     B, P2 );
 828           add(     B, 1-2*B, P2 );
 829           add(     C,     D, P3 );
 830           add(     D,     C, P3 );
 831           add( 1-C-D,     C, P3 );
 832           add( 1-C-D,     D, P3 );
 833           add(     C, 1-C-D, P3 );
 834           add(     D, 1-C-D, P3 ); break;
 835         }
 836         default:
 837           EXCEPTION( logic_error,"Invalid nb of gauss points for TRIA, variant 1: "
 838                      <<nbGauss);
 839         }
 840       }
 841       else if ( variant == 2 ) {
 842         if (geom == MED_TRIA3) setRefCoords( TTria3a() );
 843         else                setRefCoords( TTria6a() );
 844         switch ( nbGauss ) {
 845         case 1: {
 846           add( -1/3., -1/3., 2. ); break;
 847         }
 848         case 3: {
 849           add( -2/3.,  1/3., 2/3. );
 850           add( -2/3., -2/3., 2/3. );
 851           add(  1/3., -2/3., 2/3. ); break;
 852         }
 853         case 6: {
 854           const double P1 = 0.11169079483905, P2 = 0.0549758718227661;
 855           const double A  = 0.445948490915965, B = 0.091576213509771;
 856           add( 2*B-1, 1-4*B, 4*P2 );
 857           add( 2*B-1, 2*B-1, 4*P2 );
 858           add( 1-4*B, 2*B-1, 4*P2 );
 859           add( 1-4*A, 2*A-1, 4*P1 );
 860           add( 2*A-1, 1-4*A, 4*P1 );
 861           add( 2*A-1, 2*A-1, 4*P1 ); break;
 862         }
 863         default:
 864           EXCEPTION( logic_error,"Invalid nb of gauss points for TRIA, variant 2: "
 865                      <<nbGauss);
 866         }
 867       }
 868       else if ( variant == 3 ) {
 869         if (geom == MED_TRIA3) setRefCoords( TTria3b() );
 870         else                setRefCoords( TTria6b() );
 871         switch ( nbGauss ) {
 872         case 4: {
 873           add( 1/3., 1/3., -27/96 );
 874           add( 0.2 , 0.2 ,  25/96 );
 875           add( 0.6 , 0.2 ,  25/96 );
 876           add( 0.2 , 0.6 ,  25/96 ); break;
 877         }
 878         default:
 879           EXCEPTION( logic_error,"Invalid nb of gauss points for TRIA, variant 3: "
 880                      <<nbGauss);
 881         }
 882       }
 883       break;
 884
 885     case MED_QUAD4:
 886     case MED_QUAD8:
 887       if ( variant == 1 ) {
 888         if (geom == MED_QUAD4) setRefCoords( TQuad4b() );
 889         else                setRefCoords( TQuad8b() );
 890         switch ( nbGauss ) {
 891         case 1: { // FPG1
 892           add(  0,  0,  4 ); break;
 893         }
 894         case 4: { // FPG4
 895           const double a = 1/sqrt(3.);
 896           add( -a, -a,  1 );
 897           add(  a, -a,  1 );
 898           add(  a,  a,  1 );
 899           add( -a,  a,  1 ); break;
 900         }
 901         case 9: { // FPG9
 902           const double a = 0.774596669241483;
 903           add( -a, -a,  25/81. );
 904           add(  a, -a,  25/81. );
 905           add(  a,  a,  25/81. );
 906           add( -a,  a,  25/81. );
 907           add( 0., -a,  40/81. );
 908           add(  a, 0.,  40/81. );
 909           add( 0.,  a,  40/81. );
 910           add( -a, 0.,  40/81. );
 911           add( 0., 0.,  64/81. ); break;
 912         }
 913         default:
 914           EXCEPTION( logic_error,"Invalid nb of gauss points for QUAD, variant 1: "
 915                      <<nbGauss);
 916         }
 917       }
 918       else if ( variant == 2 ) {
 919         if (geom == MED_QUAD4) setRefCoords( TQuad4a() );
 920         else                setRefCoords( TQuad8a() );
 921         switch ( nbGauss ) {
 922         case 4: {
 923           const double a = 1/sqrt(3.);
 924           add( -a,  a,  1 );
 925           add( -a, -a,  1 );
 926           add(  a, -a,  1 );
 927           add(  a,  a,  1 ); break;
 928         }
 929         case 9: {
 930           const double a = 0.774596669241483;
 931           add( -a,  a,  25/81. );
 932           add( -a, -a,  25/81. );
 933           add(  a, -a,  25/81. );
 934           add(  a,  a,  25/81. );
 935           add( -a, 0.,  40/81. );
 936           add( 0., -a,  40/81. );
 937           add(  a, 0.,  40/81. );
 938           add( 0.,  a,  40/81. );
 939           add( 0., 0.,  64/81. ); break;
 940         }
 941         default:
 942           EXCEPTION( logic_error,"Invalid nb of gauss points for QUAD, variant 1: "
 943                      <<nbGauss);
 944         }
 945       }
 946       else if ( variant == 3 ) {
 947         if (geom == MED_QUAD4) setRefCoords( TQuad4b() );
 948         else                setRefCoords( TQuad8b() );
 949         switch ( nbGauss ) {
 950         case 4: {
 951           const double a = 3/sqrt(3.);
 952           add( -a, -a,  1 );
 953           add( -a,  a,  1 );
 954           add(  a, -a,  1 );
 955           add(  a,  a,  1 ); break;
 956         }
 957         case 9: {
 958           const double a = sqrt(3/5.), c1 = 5/9., c2 = 8/9.;
 959           const double c12 = c1*c2, c22 = c2*c2, c1c2 = c1*c2;
 960           add( -a, -a,  c12  );
 961           add( -a, 0.,  c1c2 );
 962           add( -a,  a,  c12  );
 963           add( 0., -a,  c1c2 );
 964           add( 0., 0.,  c22  );
 965           add( 0.,  a,  c1c2 );
 966           add(  a, -a,  c12  );
 967           add(  a, 0.,  c1c2 );
 968           add(  a,  a,  c12  ); break;
 969         }
 970         default:
 971           EXCEPTION( logic_error,"Invalid nb of gauss points for QUAD, variant 3: "
 972                      <<nbGauss);
 973         }
 974       }
 975       break;
 976
 977     case MED_TETRA4:
 978     case MED_TETRA10:
 979       if (geom == MED_TETRA4) setRefCoords( TTetra4a() );
 980       else                 setRefCoords( TTetra10a() );
 981       switch ( nbGauss ) {
 982       case 4: { // FPG4
 983         const double a = (5 - sqrt(5.))/20., b = (5 + 3*sqrt(5.))/20.;
 984         add(  a,  a,  a,  1/24. );
 985         add(  a,  a,  b,  1/24. );
 986         add(  a,  b,  a,  1/24. );
 987         add(  b,  a,  a,  1/24. ); break;
 988       }
 989       case 5: { // FPG5
 990         const double a = 0.25, b = 1/6., c = 0.5;
 991         add(  a,  a,  a, -2/15. );
 992         add(  b,  b,  b,  3/40. );
 993         add(  b,  b,  c,  3/40. );
 994         add(  b,  c,  b,  3/40. );
 995         add(  c,  b,  b,  3/40. ); break;
 996       }
 997       case 15: { // FPG15
 998         const double a = 0.25;
 999         const double b1 = (7 + sqrt(15.))/34., c1 = (13 + 3*sqrt(15.))/34., d = (5 - sqrt(15.))/20.;
1000         const double b2 = (7 - sqrt(15.))/34., c2 = (13 - 3*sqrt(15.))/34., e = (5 + sqrt(15.))/20.;
1001         const double P1 = (2665 - 14*sqrt(15.))/226800.;
1002         const double P2 = (2665 + 14*sqrt(15.))/226800.;
1003         add(  a,  a,  a,  8/405.);//_____
1004         add( b1, b1, b1,  P1    );
1005         add( b1, b1, c1,  P1    );
1006         add( b1, c1, b1,  P1    );
1007         add( c1, b1, b1,  P1    );//_____
1008         add( b2, b2, b2,  P2    );
1009         add( b2, b2, c2,  P2    );
1010         add( b2, c2, b2,  P2    );
1011         add( c2, b2, b2,  P2    );//_____
1012         add(  d,  d,  e,  5/567.);
1013         add(  d,  e,  d,  5/567.);
1014         add(  e,  d,  d,  5/567.);
1015         add(  d,  e,  e,  5/567.);
1016         add(  e,  d,  e,  5/567.);
1017         add(  e,  e,  d,  5/567.);
1018         break;
1019       }
1020       default:
1021         EXCEPTION( logic_error,"Invalid nb of gauss points for TETRA: "<<nbGauss);
1022       }
1023       break;
1024
1025     case MED_PYRA5:
1026     case MED_PYRA13:
1027       if (geom == MED_PYRA5) setRefCoords( TPyra5a() );
1028       else                setRefCoords( TPyra13a() );
1029       switch ( nbGauss ) {
1030       case 5: { // FPG5
1031         const double h1 = 0.1531754163448146;
1032         const double h2 = 0.6372983346207416;
1033         add(  .5,  0.,  h1,  2/15. );
1034         add(  0.,  .5,  h1,  2/15. );
1035         add( -.5,  0.,  h1,  2/15. );
1036         add(  0., -.5,  h1,  2/15. );
1037         add(  0.,  0.,  h2,  2/15. ); break;
1038       }
1039       case 6: { // FPG6
1040         const double p1 = 0.1024890634400000 ;
1041         const double p2 = 0.1100000000000000 ;
1042         const double p3 = 0.1467104129066667 ;
1043         const double a  = 0.5702963741068025 ;
1044         const double h1 = 0.1666666666666666 ;
1045         const double h2 = 0.08063183038464675;
1046         const double h3 = 0.6098484849057127 ;
1047         add(  a, 0.,  h1,  p1 );
1048         add( 0.,  a,  h1,  p1 );
1049         add( -a, 0.,  h1,  p1 );
1050         add( 0., -a,  h1,  p1 );
1051         add( 0., 0.,  h2,  p2 );
1052         add( 0., 0.,  h3,  p3 ); break;
1053       }
1054       case 27: { // FPG27
1055         const double a1  = 0.788073483;
1056         const double b6  = 0.499369002;
1057         const double b1  = 0.848418011;
1058         const double c8  = 0.478508449;
1059         const double c1  = 0.652816472;
1060         const double d12 = 0.032303742;
1061         const double d1  = 1.106412899;
1062         double z = 1/2., fz = b1/2*(1 - z);
1063         add(  0.,  0.,   z,  a1 ); // 1
1064         add(  fz,  fz,   z,  b6 ); // 2
1065         add( -fz,  fz,   z,  b6 ); // 3
1066         add( -fz, -fz,   z,  b6 ); // 4
1067         add(  fz, -fz,   z,  b6 ); // 5
1068         z = (1 - b1)/2.;
1069         add(  0.,  0.,   z,  b6 ); // 6
1070         z = (1 + b1)/2.;
1071         add(  0.,  0.,   z,  b6 ); // 7
1072         z = (1 - c1)/2.; fz = c1*(1 - z);
1073         add(  fz,  0.,   z,  c8 ); // 8
1074         add(  0.,  fz,   z,  c8 ); // 9
1075         add( -fz,  0.,   z,  c8 ); // 10
1076         add(  0., -fz,   z,  c8 ); // 11
1077         z = (1 + c1)/2.; fz = c1*(1 - z);
1078         add(  fz,  0.,   z,  c8 ); // 12
1079         add(  0.,  fz,   z,  c8 ); // 13
1080         add( -fz,  0.,   z,  c8 ); // 14
1081         add(  0., -fz,   z,  c8 ); // 15
1082         z = (1 - d1)/2., fz = d1/2*(1 - z);
1083         add(  fz,  fz,   z,  d12); // 16
1084         add( -fz,  fz,   z,  d12); // 17
1085         add( -fz, -fz,   z,  d12); // 18
1086         add(  fz, -fz,   z,  d12); // 19
1087         z = 1/2.; fz = d1*(1 - z);
1088         add(  fz,  0.,   z,  d12); // 20
1089         add(  0.,  fz,   z,  d12); // 21
1090         add( -fz,  0.,   z,  d12); // 22
1091         add(  0., -fz,   z,  d12); // 23
1092         z = (1 + d1)/2., fz = d1/2*(1 - z);
1093         add(  fz,  fz,   z,  d12); // 24
1094         add( -fz,  fz,   z,  d12); // 25
1095         add( -fz, -fz,   z,  d12); // 26
1096         add(  fz, -fz,   z,  d12); // 27
1097         break;
1098       }
1099       default:
1100         EXCEPTION( logic_error,"Invalid nb of gauss points for PYRA: "<<nbGauss);
1101       }
1102       break;
1103     case MED_PENTA6:
1104     case MED_PENTA15:
1105       if (geom == MED_PENTA6) setRefCoords( TPenta6a() );
1106       else                 setRefCoords( TPenta15a() );
1107       switch ( nbGauss ) {
1108       case 6: { // FPG6
1109         const double a = sqrt(3.)/3.;
1110         add( -a, .5, .5,  1/6. );
1111         add( -a, 0., .5,  1/6. );
1112         add( -a, .5, 0.,  1/6. );
1113         add(  a, .5, .5,  1/6. );
1114         add(  a, 0., .5,  1/6. );
1115         add(  a, .5, 0.,  1/6. ); break;
1116       }
1117       case 8: { // FPG8
1118         const double a = 0.577350269189626;
1119         add( -a, 1/3., 1/3., -27/96. );
1120         add( -a,  0.6,  0.2,  25/96. );
1121         add( -a,  0.2,  0.6,  25/96. );
1122         add( -a,  0.2,  0.2,  25/96. );
1123         add( +a, 1/3., 1/3., -27/96. );
1124         add( +a,  0.6,  0.2,  25/96. );
1125         add( +a,  0.2,  0.6,  25/96. );
1126         add( +a,  0.2,  0.2,  25/96. ); break;
1127       }
1128       case 21: { // FPG21
1129         const double d = sqrt(3/5.), c1 = 5/9., c2 = 8/9.; // d <=> alfa
1130         const double a = (6 + sqrt(15.))/21.;
1131         const double b = (6 - sqrt(15.))/21.;
1132         const double P1 = (155 + sqrt(15.))/2400.;
1133         const double P2 = (155 - sqrt(15.))/2400.;  //___
1134         add( -d,  1/3.,  1/3., c1*9/80. );//___
1135         add( -d,     a,     a, c1*P1    );
1136         add( -d, 1-2*a,     a, c1*P1    );
1137         add( -d,     a, 1-2*a, c1*P1    );//___
1138         add( -d,     b,     b, c1*P2    );
1139         add( -d, 1-2*b,     b, c1*P2    );
1140         add( -d,     b, 1-2*b, c1*P2    );//___
1141         add( 0.,  1/3.,  1/3., c2*9/80. );//___
1142         add( 0.,     a,     a, c2*P1    );
1143         add( 0., 1-2*a,     a, c2*P1    );
1144         add( 0.,     a, 1-2*a, c2*P1    );//___
1145         add( 0.,     b,     b, c2*P2    );
1146         add( 0., 1-2*b,     b, c2*P2    );
1147         add( 0.,     b, 1-2*b, c2*P2    );//___
1148         add(  d,  1/3.,  1/3., c1*9/80. );//___
1149         add(  d,     a,     a, c1*P1    );
1150         add(  d, 1-2*a,     a, c1*P1    );
1151         add(  d,     a, 1-2*a, c1*P1    );//___
1152         add(  d,     b,     b, c1*P2    );
1153         add(  d, 1-2*b,     b, c1*P2    );
1154         add(  d,     b, 1-2*b, c1*P2    );//___
1155         break;
1156       }
1157       default:
1158         EXCEPTION( logic_error,"Invalid nb of gauss points for PENTA: " <<nbGauss);
1159       }
1160       break;
1161
1162     case MED_HEXA8:
1163     case MED_HEXA20:
1164       if (geom == MED_HEXA8) setRefCoords( THexa8a() );
1165       else                setRefCoords( THexa20a() );
1166       switch ( nbGauss ) {
1167       case 8: { // FPG8
1168         const double a = sqrt(3.)/3.;
1169         add( -a, -a, -a,  1. );
1170         add( -a, -a,  a,  1. );
1171         add( -a,  a, -a,  1. );
1172         add( -a,  a,  a,  1. );
1173         add(  a, -a, -a,  1. );
1174         add(  a, -a,  a,  1. );
1175         add(  a,  a, -a,  1. );
1176         add(  a,  a,  a,  1. ); break;
1177       }
1178       case 27: { // FPG27
1179         const double a = sqrt(3/5.), c1 = 5/9., c2 = 8/9.;
1180         const double c12 = c1*c1, c13 = c1*c1*c1;
1181         const double c22 = c2*c2, c23 = c2*c2*c2;
1182         add( -a, -a, -a,   c13  ); // 1
1183         add( -a, -a, 0., c12*c2 ); // 2
1184         add( -a, -a,  a,   c13  ); // 3
1185         add( -a, 0., -a, c12*c2 ); // 4
1186         add( -a, 0., 0., c1*c22 ); // 5
1187         add( -a, 0.,  a, c12*c2 ); // 6
1188         add( -a,  a, -a,   c13  ); // 7
1189         add( -a,  a, 0., c12*c2 ); // 8
1190         add( -a,  a,  a,   c13  ); // 9
1191         add( 0., -a, -a, c12*c2 ); // 10
1192         add( 0., -a, 0., c1*c22 ); // 11
1193         add( 0., -a,  a, c12*c2 ); // 12
1194         add( 0., 0., -a, c1*c22 ); // 13
1195         add( 0., 0., 0.,   c23  ); // 14
1196         add( 0., 0.,  a, c1*c22 ); // 15
1197         add( 0.,  a, -a, c12*c2 ); // 16
1198         add( 0.,  a, 0., c1*c22 ); // 17
1199         add( 0.,  a,  a, c12*c2 ); // 18
1200         add(  a, -a, -a,   c13  ); // 19
1201         add(  a, -a, 0., c12*c2 ); // 20
1202         add(  a, -a,  a,   c13  ); // 21
1203         add(  a, 0., -a, c12*c2 ); // 22
1204         add(  a, 0., 0., c1*c22 ); // 23
1205         add(  a, 0.,  a, c12*c2 ); // 24
1206         add(  a,  a, -a,   c13  ); // 25
1207         add(  a,  a, 0., c12*c2 ); // 26
1208         add(  a,  a,  a,   c13  ); // 27
1209         break;
1210       }
1211       default:
1212         EXCEPTION( logic_error,"Invalid nb of gauss points for PENTA: " <<nbGauss);
1213       }
1214       break;
1215
1216     default:
1217       EXCEPTION( logic_error,"unexpected EGeometrieElement: "<< geom);
1218     }
1219
1220     if ( myWeights.capacity() != myWeights.size() )
1221       EXCEPTION( logic_error,"Not all gauss points defined");
1222   }
1223 }
1224
1225 //=======================================================================
1226 /*!
1227  * Creates a localization filled with default values. The caller gets pointer ownership
1228  */
1229 //=======================================================================
1230
1231 namespace MEDMEM {
1232
1233   GAUSS_LOCALIZATION_*
1234   GAUSS_LOCALIZATION_::makeDefaultLocalization(const string &     locName,
1235                                                medGeometryElement typeGeo,
1236                                                int                nGauss) throw (MEDEXCEPTION)
1237   {
1238     TGaussDef gaussDef( typeGeo, nGauss, 1 );
1239
1240     GAUSS_LOCALIZATION_ * gsloc =
1241       new GAUSS_LOCALIZATION<FullInterlace> ( locName.c_str(),
1242                                               typeGeo,
1243                                               nGauss,
1244                                               &gaussDef.myRefCoords[0],
1245                                               &gaussDef.myCoords[0],
1246                                               &gaussDef.myWeights[0] );
1247     return gsloc;
1248   }
1249 }