src/MEDLoader/SauvMedConvertor.cxx

   1 // Copyright (C) 2007-2021  CEA/DEN, EDF R&D
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License, or (at your option) any later version.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19 // File      : SauvMedConvertor.cxx
  20 // Created   : Tue Aug 16 14:43:20 2011
  21 // Author    : Edward AGAPOV (eap)
  22 //
  23
  24 #include "SauvMedConvertor.hxx"
  25
  26 #include "CellModel.hxx"
  27 #include "MEDFileMesh.hxx"
  28 #include "MEDFileField.hxx"
  29 #include "MEDFileData.hxx"
  30 #include "MEDCouplingFieldDouble.hxx"
  31
  32 #include <iostream>
  33 #include <cassert>
  34 #include <cmath>
  35 #include <queue>
  36 #include <limits>
  37
  38 #include <cstdlib>
  39 #include <cstring>
  40 #include <fcntl.h>
  41
  42 #ifdef WIN32
  43 #include <io.h>
  44 #else
  45 #include <unistd.h>
  46 #endif
  47
  48 #ifdef HAS_XDR
  49 #include <rpc/types.h>
  50 #include <rpc/xdr.h>
  51 #endif
  52
  53 using namespace SauvUtilities;
  54 using namespace MEDCoupling;
  55
  56 namespace
  57 {
  58   // for ASCII file reader
  59   const int GIBI_MaxOutputLen = 150; // max length of a line in the sauve file
  60   const int GIBI_BufferSize   = 16184; // for buffered reading
  61
  62   using namespace INTERP_KERNEL;
  63
  64   const size_t MaxMedCellType = NORM_ERROR;
  65   const size_t NbGibiCellTypes = 47;
  66   const TCellType GibiTypeToMed[NbGibiCellTypes] =
  67     {
  68       /*1 */ NORM_POINT1 ,/*2 */ NORM_SEG2   ,/*3 */ NORM_SEG3   ,/*4 */ NORM_TRI3   ,/*5 */ NORM_ERROR  ,
  69       /*6 */ NORM_TRI6   ,/*7 */ NORM_ERROR  ,/*8 */ NORM_QUAD4  ,/*9 */ NORM_ERROR  ,/*10*/ NORM_QUAD8  ,
  70       /*11*/ NORM_ERROR  ,/*12*/ NORM_ERROR  ,/*13*/ NORM_ERROR  ,/*14*/ NORM_HEXA8  ,/*15*/ NORM_HEXA20 ,
  71       /*16*/ NORM_PENTA6 ,/*17*/ NORM_PENTA15,/*18*/ NORM_ERROR  ,/*19*/ NORM_ERROR  ,/*20*/ NORM_ERROR  ,
  72       /*21*/ NORM_ERROR  ,/*22*/ NORM_ERROR  ,/*23*/ NORM_TETRA4 ,/*24*/ NORM_TETRA10,/*25*/ NORM_PYRA5  ,
  73       /*26*/ NORM_PYRA13 ,/*27*/ NORM_ERROR  ,/*28*/ NORM_ERROR  ,/*29*/ NORM_ERROR  ,/*30*/ NORM_ERROR  ,
  74       /*31*/ NORM_ERROR  ,/*32*/ NORM_ERROR  ,/*33*/ NORM_ERROR  ,/*34*/ NORM_ERROR  ,/*35*/ NORM_ERROR  ,
  75       /*36*/ NORM_ERROR  ,/*37*/ NORM_ERROR  ,/*38*/ NORM_ERROR  ,/*39*/ NORM_ERROR  ,/*40*/ NORM_ERROR  ,
  76       /*41*/ NORM_ERROR  ,/*42*/ NORM_ERROR  ,/*43*/ NORM_ERROR  ,/*44*/ NORM_ERROR  ,/*45*/ NORM_ERROR  ,
  77       /*46*/ NORM_ERROR  ,/*47*/ NORM_ERROR
  78     };
  79
  80   //================================================================================
  81   /*!
  82    * \brief Return dimension of a group
  83    */
  84   //================================================================================
  85
  86   unsigned getDim( const Group* grp )
  87   {
  88     return SauvUtilities::getDimension( grp->_groups.empty() ? grp->_cellType : grp->_groups[0]->_cellType );
  89   }
  90
  91   //================================================================================
  92   /*!
  93    * \brief Converts connectivity of quadratic elements
  94    */
  95   //================================================================================
  96
  97   inline void ConvertQuadratic( const INTERP_KERNEL::NormalizedCellType type,
  98                                 const Cell &                            aCell )
  99   {
 100     if ( const int * conn = getGibi2MedQuadraticInterlace( type ))
 101       {
 102         Cell* ma = const_cast<Cell*>(&aCell);
 103         std::vector< Node* > new_nodes( ma->_nodes.size() );
 104         for (std:: size_t i = 0; i < new_nodes.size(); ++i )
 105           new_nodes[ i ] = ma->_nodes[ conn[ i ]];
 106         ma->_nodes.swap( new_nodes );
 107       }
 108   }
 109
 110   //================================================================================
 111   /*!
 112    * \brief Returns a vector of pairs of node indices to inverse a med volume element
 113    */
 114   //================================================================================
 115
 116   void getReverseVector (const INTERP_KERNEL::NormalizedCellType type,
 117                          std::vector<std::pair<int,int> > &                swapVec )
 118   {
 119     swapVec.clear();
 120
 121     switch ( type )
 122       {
 123       case NORM_TETRA4:
 124         swapVec.resize(1);
 125         swapVec[0] = std::make_pair( 1, 2 );
 126         break;
 127       case NORM_PYRA5:
 128         swapVec.resize(1);
 129         swapVec[0] = std::make_pair( 1, 3 );
 130         break;
 131       case NORM_PENTA6:
 132         swapVec.resize(2);
 133         swapVec[0] = std::make_pair( 1, 2 );
 134         swapVec[1] = std::make_pair( 4, 5 );
 135         break;
 136       case NORM_HEXA8:
 137         swapVec.resize(2);
 138         swapVec[0] = std::make_pair( 1, 3 );
 139         swapVec[1] = std::make_pair( 5, 7 );
 140         break;
 141       case NORM_TETRA10:
 142         swapVec.resize(3);
 143         swapVec[0] = std::make_pair( 1, 2 );
 144         swapVec[1] = std::make_pair( 4, 6 );
 145         swapVec[2] = std::make_pair( 8, 9 );
 146         break;
 147       case NORM_PYRA13:
 148         swapVec.resize(4);
 149         swapVec[0] = std::make_pair( 1, 3 );
 150         swapVec[1] = std::make_pair( 5, 8 );
 151         swapVec[2] = std::make_pair( 6, 7 );
 152         swapVec[3] = std::make_pair( 10, 12 );
 153         break;
 154       case NORM_PENTA15:
 155         swapVec.resize(4);
 156         swapVec[0] = std::make_pair( 1, 2 );
 157         swapVec[1] = std::make_pair( 4, 5 );
 158         swapVec[2] = std::make_pair( 6, 8 );
 159         swapVec[3] = std::make_pair( 9, 11 );
 160         break;
 161       case NORM_HEXA20:
 162         swapVec.resize(7);
 163         swapVec[0] = std::make_pair( 1, 3 );
 164         swapVec[1] = std::make_pair( 5, 7 );
 165         swapVec[2] = std::make_pair( 8, 11 );
 166         swapVec[3] = std::make_pair( 9, 10 );
 167         swapVec[4] = std::make_pair( 12, 15 );
 168         swapVec[5] = std::make_pair( 13, 14 );
 169         swapVec[6] = std::make_pair( 17, 19 );
 170         break;
 171         //   case NORM_SEG3: no need to reverse edges
 172         //     swapVec.resize(1);
 173         //     swapVec[0] = std::make_pair( 1, 2 );
 174         //     break;
 175       case NORM_TRI6:
 176         swapVec.resize(2);
 177         swapVec[0] = std::make_pair( 1, 2 );
 178         swapVec[1] = std::make_pair( 3, 5 );
 179         break;
 180       case NORM_QUAD8:
 181         swapVec.resize(3);
 182         swapVec[0] = std::make_pair( 1, 3 );
 183         swapVec[1] = std::make_pair( 4, 7 );
 184         swapVec[2] = std::make_pair( 5, 6 );
 185         break;
 186       default:;
 187       }
 188   }
 189
 190   //================================================================================
 191   /*!
 192    * \brief Inverses element orientation using vector of indices to swap
 193    */
 194   //================================================================================
 195
 196   inline void reverse(const Cell & aCell, const std::vector<std::pair<int,int> > & swapVec )
 197   {
 198     Cell* ma = const_cast<Cell*>(&aCell);
 199     for ( unsigned i = 0; i < swapVec.size(); ++i )
 200       std::swap( ma->_nodes[ swapVec[i].first ],
 201                  ma->_nodes[ swapVec[i].second ]);
 202     if ( swapVec.empty() )
 203       ma->_reverse = true;
 204     else
 205       ma->_reverse = false;
 206   }
 207   //================================================================================
 208   /*!
 209    * \brief Comparator of cells by number used for ordering cells within a med group
 210    */
 211   struct TCellByIDCompare
 212   {
 213     bool operator () (const Cell* i1, const Cell* i2) const
 214     {
 215       return i1->_number < i2->_number;
 216     }
 217   };
 218   typedef std::map< const Cell*, unsigned, TCellByIDCompare > TCellToOrderMap;
 219
 220   //================================================================================
 221   /*!
 222    * \brief Fill Group::_relocTable if necessary
 223    */
 224   //================================================================================
 225
 226   void setRelocationTable( Group* grp, TCellToOrderMap& cell2order )
 227   {
 228     if ( !grp->_isProfile ) return;
 229
 230     // check if relocation table is necessary
 231     bool isRelocated = false;
 232     unsigned newOrder = 0;
 233     TCellToOrderMap::iterator c2oIt = cell2order.begin(), c2oEnd = cell2order.end();
 234     for ( ; !isRelocated && c2oIt != c2oEnd; ++c2oIt, ++newOrder )
 235       isRelocated = ( c2oIt->second != newOrder );
 236
 237     if ( isRelocated )
 238       {
 239         grp->_relocTable.resize( cell2order.size() );
 240         for ( newOrder = 0, c2oIt = cell2order.begin(); c2oIt != c2oEnd; ++c2oIt, ++newOrder )
 241           grp->_relocTable[ c2oIt->second ] = newOrder;
 242       }
 243   }
 244 }
 245
 246 namespace // define default GAUSS points
 247 {
 248   typedef std::vector<double> TDoubleVector;
 249   typedef double*             TCoordSlice;
 250   typedef int                 TInt;
 251   //---------------------------------------------------------------
 252   //! Shape function definitions
 253   //---------------------------------------------------------------
 254   struct TShapeFun
 255   {
 256     TInt myDim;
 257     TInt myNbRef;
 258     TDoubleVector myRefCoord;
 259
 260     TShapeFun(TInt theDim = 0, TInt theNbRef = 0)
 261       :myDim(theDim),myNbRef(theNbRef),myRefCoord(theNbRef*theDim) {}
 262
 263     TInt GetNbRef() const { return myNbRef; }
 264
 265     TCoordSlice GetCoord(std::size_t theRefId) { return &myRefCoord[0] + theRefId * myDim; }
 266   };
 267   //---------------------------------------------------------------
 268   /*!
 269    * \brief Description of family of integration points
 270    */
 271   //---------------------------------------------------------------
 272   struct TGaussDef
 273   {
 274     int           myType;      //!< element geometry (EGeometrieElement or med_geometrie_element)
 275     TDoubleVector myRefCoords; //!< description of reference points
 276     TDoubleVector myCoords;    //!< coordinates of Gauss points
 277     TDoubleVector myWeights;   //!< weights, len(weights)==<nb of gauss points>
 278
 279     /*!
 280      * \brief Creates definition of gauss points family
 281      *  \param geomType - element geometry (EGeometrieElement or med_geometrie_element)
 282      *  \param nbPoints - nb gauss point
 283      *  \param variant - [1-3] to choose the variant of definition
 284      *
 285      * Throws in case of invalid parameters
 286      * variant == 1 refers to "Fonctions de forme et points d'integration
 287      *              des elements finis" v7.4 by J. PELLET, X. DESROCHES, 15/09/05
 288      * variant == 2 refers to the same doc v6.4 by J.P. LEFEBVRE, X. DESROCHES, 03/07/03
 289      * variant == 3 refers to the same doc v6.4, second variant for 2D elements
 290      */
 291     TGaussDef(const int geomType, const int nbPoints, const int variant=1);
 292
 293     int dim() const { return SauvUtilities::getDimension( NormalizedCellType( myType )); }
 294     std::size_t nbPoints() const { return myWeights.capacity(); }
 295
 296   private:
 297     void add(const double x, const double weight);
 298     void add(const double x, const double y, const double weight);
 299     void add(const double x, const double y, const double z, const double weight);
 300     void setRefCoords(const TShapeFun& aShapeFun) { myRefCoords = aShapeFun.myRefCoord; }
 301   };
 302   struct TSeg2a: TShapeFun {
 303     TSeg2a();
 304   };
 305   struct TSeg3a: TShapeFun {
 306     TSeg3a();
 307   };
 308   struct TTria3a: TShapeFun {
 309     TTria3a();
 310   };
 311   struct TTria6a: TShapeFun {
 312     TTria6a();
 313   };
 314   struct TTria3b: TShapeFun {
 315     TTria3b();
 316   };
 317   struct TTria6b: TShapeFun {
 318     TTria6b();
 319   };
 320   struct TQuad4a: TShapeFun {
 321     TQuad4a();
 322   };
 323   struct TQuad8a: TShapeFun {
 324     TQuad8a();
 325   };
 326   struct TQuad4b: TShapeFun {
 327     TQuad4b();
 328   };
 329   struct TQuad8b: TShapeFun {
 330     TQuad8b();
 331   };
 332   struct TTetra4a: TShapeFun {
 333     TTetra4a();
 334   };
 335   struct TTetra10a: TShapeFun {
 336     TTetra10a();
 337   };
 338   struct TTetra4b: TShapeFun {
 339     TTetra4b();
 340   };
 341   struct TTetra10b: TShapeFun {
 342     TTetra10b();
 343   };
 344   struct THexa8a: TShapeFun {
 345     THexa8a();
 346   };
 347   struct THexa20a: TShapeFun {
 348     THexa20a(TInt theDim = 3, TInt theNbRef = 20);
 349   };
 350   struct THexa27a: THexa20a {
 351     THexa27a();
 352   };
 353   struct THexa8b: TShapeFun {
 354     THexa8b();
 355   };
 356   struct THexa20b: TShapeFun {
 357     THexa20b(TInt theDim = 3, TInt theNbRef = 20);
 358   };
 359   struct TPenta6a: TShapeFun {
 360     TPenta6a();
 361   };
 362   struct TPenta6b: TShapeFun {
 363     TPenta6b();
 364   };
 365   struct TPenta15a: TShapeFun {
 366     TPenta15a();
 367   };
 368   struct TPenta15b: TShapeFun {
 369     TPenta15b();
 370   };
 371   struct TPyra5a: TShapeFun {
 372     TPyra5a();
 373   };
 374   struct TPyra5b: TShapeFun {
 375     TPyra5b();
 376   };
 377   struct TPyra13a: TShapeFun {
 378     TPyra13a();
 379   };
 380   struct TPyra13b: TShapeFun {
 381     TPyra13b();
 382   };
 383
 384   void TGaussDef::add(const double x, const double weight)
 385   {
 386     if ( dim() != 1 )
 387       THROW_IK_EXCEPTION("TGaussDef: dim() != 1");
 388     if ( myWeights.capacity() == myWeights.size() )
 389       THROW_IK_EXCEPTION("TGaussDef: Extra gauss point");
 390     myCoords.push_back( x );
 391     myWeights.push_back( weight );
 392   }
 393   void TGaussDef::add(const double x, const double y, const double weight)
 394   {
 395     if ( dim() != 2 )
 396       THROW_IK_EXCEPTION("TGaussDef: dim() != 2");
 397     if ( myWeights.capacity() == myWeights.size() )
 398       THROW_IK_EXCEPTION("TGaussDef: Extra gauss point");
 399     myCoords.push_back( x );
 400     myCoords.push_back( y );
 401     myWeights.push_back( weight );
 402   }
 403   void TGaussDef::add(const double x, const double y, const double z, const double weight)
 404   {
 405     if ( dim() != 3 )
 406       THROW_IK_EXCEPTION("TGaussDef: dim() != 3");
 407     if ( myWeights.capacity() == myWeights.size() )
 408       THROW_IK_EXCEPTION("TGaussDef: Extra gauss point");
 409     myCoords.push_back( x );
 410     myCoords.push_back( y );
 411     myCoords.push_back( z );
 412     myWeights.push_back( weight );
 413   }
 414
 415   //---------------------------------------------------------------
 416   TSeg2a::TSeg2a():TShapeFun(1,2)
 417   {
 418     TInt aNbRef = GetNbRef();
 419     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 420       TCoordSlice aCoord = GetCoord(aRefId);
 421       switch(aRefId){
 422       case  0: aCoord[0] = -1.0; break;
 423       case  1: aCoord[0] =  1.0; break;
 424       }
 425     }
 426   }
 427   //---------------------------------------------------------------
 428   TSeg3a::TSeg3a():TShapeFun(1,3)
 429   {
 430     TInt aNbRef = GetNbRef();
 431     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 432       TCoordSlice aCoord = GetCoord(aRefId);
 433       switch(aRefId){
 434       case  0: aCoord[0] = -1.0; break;
 435       case  1: aCoord[0] =  1.0; break;
 436       case  2: aCoord[0] =  0.0; break;
 437       }
 438     }
 439   }
 440   //---------------------------------------------------------------
 441   TTria3a::TTria3a():
 442     TShapeFun(2,3)
 443   {
 444     TInt aNbRef = GetNbRef();
 445     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 446       TCoordSlice aCoord = GetCoord(aRefId);
 447       switch(aRefId){
 448       case  0: aCoord[0] = -1.0;  aCoord[1] =  1.0; break;
 449       case  1: aCoord[0] = -1.0;  aCoord[1] = -1.0; break;
 450       case  2: aCoord[0] =  1.0;  aCoord[1] = -1.0; break;
 451       }
 452     }
 453   }
 454   //---------------------------------------------------------------
 455   TTria6a::TTria6a():TShapeFun(2,6)
 456   {
 457     TInt aNbRef = GetNbRef();
 458     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 459       TCoordSlice aCoord = GetCoord(aRefId);
 460       switch(aRefId){
 461       case  0: aCoord[0] = -1.0;  aCoord[1] =  1.0; break;
 462       case  1: aCoord[0] = -1.0;  aCoord[1] = -1.0; break;
 463       case  2: aCoord[0] =  1.0;  aCoord[1] = -1.0; break;
 464
 465       case  3: aCoord[0] = -1.0;  aCoord[1] =  1.0; break;
 466       case  4: aCoord[0] =  0.0;  aCoord[1] = -1.0; break;
 467       case  5: aCoord[0] =  0.0;  aCoord[1] =  0.0; break;
 468       }
 469     }
 470   }
 471   //---------------------------------------------------------------
 472   TTria3b::TTria3b():
 473     TShapeFun(2,3)
 474   {
 475     TInt aNbRef = GetNbRef();
 476     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 477       TCoordSlice aCoord = GetCoord(aRefId);
 478       switch(aRefId){
 479       case  0: aCoord[0] =  0.0;  aCoord[1] =  0.0; break;
 480       case  1: aCoord[0] =  1.0;  aCoord[1] =  0.0; break;
 481       case  2: aCoord[0] =  0.0;  aCoord[1] =  1.0; break;
 482       }
 483     }
 484   }
 485   //---------------------------------------------------------------
 486   TTria6b::TTria6b():
 487     TShapeFun(2,6)
 488   {
 489     TInt aNbRef = GetNbRef();
 490     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 491       TCoordSlice aCoord = GetCoord(aRefId);
 492       switch(aRefId){
 493       case  0: aCoord[0] =  0.0;  aCoord[1] =  0.0; break;
 494       case  1: aCoord[0] =  1.0;  aCoord[1] =  0.0; break;
 495       case  2: aCoord[0] =  0.0;  aCoord[1] =  1.0; break;
 496
 497       case  3: aCoord[0] =  0.5;  aCoord[1] =  0.0; break;
 498       case  4: aCoord[0] =  0.5;  aCoord[1] =  0.5; break;
 499       case  5: aCoord[0] =  0.0;  aCoord[1] =  0.5; break;
 500       }
 501     }
 502   }
 503   //---------------------------------------------------------------
 504   TQuad4a::TQuad4a():
 505     TShapeFun(2,4)
 506   {
 507     TInt aNbRef = GetNbRef();
 508     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 509       TCoordSlice aCoord = GetCoord(aRefId);
 510       switch(aRefId){
 511       case  0: aCoord[0] = -1.0;  aCoord[1] =  1.0; break;
 512       case  1: aCoord[0] = -1.0;  aCoord[1] = -1.0; break;
 513       case  2: aCoord[0] =  1.0;  aCoord[1] = -1.0; break;
 514       case  3: aCoord[0] =  1.0;  aCoord[1] =  1.0; break;
 515       }
 516     }
 517   }
 518   //---------------------------------------------------------------
 519   TQuad8a::TQuad8a():
 520     TShapeFun(2,8)
 521   {
 522     TInt aNbRef = GetNbRef();
 523     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 524       TCoordSlice aCoord = GetCoord(aRefId);
 525       switch(aRefId){
 526       case  0: aCoord[0] = -1.0;  aCoord[1] =  1.0; break;
 527       case  1: aCoord[0] = -1.0;  aCoord[1] = -1.0; break;
 528       case  2: aCoord[0] =  1.0;  aCoord[1] = -1.0; break;
 529       case  3: aCoord[0] =  1.0;  aCoord[1] =  1.0; break;
 530
 531       case  4: aCoord[0] = -1.0;  aCoord[1] =  0.0; break;
 532       case  5: aCoord[0] =  0.0;  aCoord[1] = -1.0; break;
 533       case  6: aCoord[0] =  1.0;  aCoord[1] =  0.0; break;
 534       case  7: aCoord[0] =  0.0;  aCoord[1] =  1.0; break;
 535       }
 536     }
 537   }
 538   //---------------------------------------------------------------
 539   TQuad4b::TQuad4b():
 540     TShapeFun(2,4)
 541   {
 542     TInt aNbRef = GetNbRef();
 543     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 544       TCoordSlice aCoord = GetCoord(aRefId);
 545       switch(aRefId){
 546       case  0: aCoord[0] = -1.0;  aCoord[1] = -1.0; break;
 547       case  1: aCoord[0] =  1.0;  aCoord[1] = -1.0; break;
 548       case  2: aCoord[0] =  1.0;  aCoord[1] =  1.0; break;
 549       case  3: aCoord[0] = -1.0;  aCoord[1] =  1.0; break;
 550       }
 551     }
 552   }
 553   //---------------------------------------------------------------
 554   TQuad8b::TQuad8b():
 555     TShapeFun(2,8)
 556   {
 557     TInt aNbRef = GetNbRef();
 558     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 559       TCoordSlice aCoord = GetCoord(aRefId);
 560       switch(aRefId){
 561       case  0: aCoord[0] = -1.0;  aCoord[1] = -1.0; break;
 562       case  1: aCoord[0] =  1.0;  aCoord[1] = -1.0; break;
 563       case  2: aCoord[0] =  1.0;  aCoord[1] =  1.0; break;
 564       case  3: aCoord[0] = -1.0;  aCoord[1] =  1.0; break;
 565
 566       case  4: aCoord[0] =  0.0;  aCoord[1] = -1.0; break;
 567       case  5: aCoord[0] =  1.0;  aCoord[1] =  0.0; break;
 568       case  6: aCoord[0] =  0.0;  aCoord[1] =  1.0; break;
 569       case  7: aCoord[0] = -1.0;  aCoord[1] =  0.0; break;
 570       }
 571     }
 572   }
 573   //---------------------------------------------------------------
 574   TTetra4a::TTetra4a():
 575     TShapeFun(3,4)
 576   {
 577     TInt aNbRef = GetNbRef();
 578     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 579       TCoordSlice aCoord = GetCoord(aRefId);
 580       switch(aRefId){
 581       case  0: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
 582       case  1: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
 583       case  2: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 584       case  3: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 585       }
 586     }
 587   }
 588   //---------------------------------------------------------------
 589   TTetra10a::TTetra10a():
 590     TShapeFun(3,10)
 591   {
 592     TInt aNbRef = GetNbRef();
 593     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 594       TCoordSlice aCoord = GetCoord(aRefId);
 595       switch(aRefId){
 596       case  0: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
 597       case  1: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
 598       case  2: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 599       case  3: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 600
 601       case  4: aCoord[0] =  0.0;  aCoord[1] =  0.5;  aCoord[2] =  0.5; break;
 602       case  5: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
 603       case  6: aCoord[0] =  0.0;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
 604
 605       case  7: aCoord[0] =  0.5;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
 606       case  8: aCoord[0] =  0.5;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
 607       case  9: aCoord[0] =  0.5;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 608       }
 609     }
 610   }
 611   //---------------------------------------------------------------
 612   TTetra4b::TTetra4b():
 613     TShapeFun(3,4)
 614   {
 615     TInt aNbRef = GetNbRef();
 616     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 617       TCoordSlice aCoord = GetCoord(aRefId);
 618       switch(aRefId){
 619       case  0: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
 620       case  2: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
 621       case  1: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 622       case  3: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 623       }
 624     }
 625   }
 626   //---------------------------------------------------------------
 627   TTetra10b::TTetra10b():
 628     TShapeFun(3,10)
 629   {
 630     TInt aNbRef = GetNbRef();
 631     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 632       TCoordSlice aCoord = GetCoord(aRefId);
 633       switch(aRefId){
 634       case  0: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
 635       case  2: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
 636       case  1: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 637       case  3: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 638
 639       case  6: aCoord[0] =  0.0;  aCoord[1] =  0.5;  aCoord[2] =  0.5; break;
 640       case  5: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
 641       case  4: aCoord[0] =  0.0;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
 642
 643       case  7: aCoord[0] =  0.5;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
 644       case  9: aCoord[0] =  0.5;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
 645       case  8: aCoord[0] =  0.5;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 646       }
 647     }
 648   }
 649   //---------------------------------------------------------------
 650   THexa8a::THexa8a():
 651     TShapeFun(3,8)
 652   {
 653     TInt aNbRef = GetNbRef();
 654     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 655       TCoordSlice aCoord = GetCoord(aRefId);
 656       switch(aRefId){
 657       case  0: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
 658       case  1: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
 659       case  2: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
 660       case  3: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
 661       case  4: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
 662       case  5: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
 663       case  6: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
 664       case  7: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
 665       }
 666     }
 667   }
 668   //---------------------------------------------------------------
 669   THexa20a::THexa20a(TInt theDim, TInt theNbRef):
 670     TShapeFun(theDim,theNbRef)
 671   {
 672     std::size_t aNbRef = myRefCoord.size();
 673     for(std::size_t aRefId = 0; aRefId < aNbRef; aRefId++){
 674       TCoordSlice aCoord = GetCoord(aRefId);
 675       switch(aRefId){
 676       case  0: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
 677       case  1: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
 678       case  2: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
 679       case  3: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
 680       case  4: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
 681       case  5: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
 682       case  6: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
 683       case  7: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
 684
 685       case  8: aCoord[0] =  0.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
 686       case  9: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] = -1.0; break;
 687       case 10: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
 688       case 11: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] = -1.0; break;
 689       case 12: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] =  0.0; break;
 690       case 13: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] =  0.0; break;
 691       case 14: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
 692       case 15: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
 693       case 16: aCoord[0] =  0.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
 694       case 17: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
 695       case 18: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
 696       case 19: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
 697       }
 698     }
 699   }
 700   //---------------------------------------------------------------
 701   THexa27a::THexa27a():
 702     THexa20a(3,27)
 703   {
 704     std::size_t aNbRef = myRefCoord.size();
 705     for(std::size_t aRefId = 0; aRefId < aNbRef; aRefId++){
 706       TCoordSlice aCoord = GetCoord(aRefId);
 707       switch(aRefId){
 708       case 20: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] = -1.0; break;
 709       case 21: aCoord[0] =  0.0;  aCoord[1] = -1.0;  aCoord[2] =  0.0; break;
 710       case 22: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 711       case 23: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
 712       case 24: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 713       case 25: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
 714       case 26: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 715       }
 716     }
 717   }
 718   //---------------------------------------------------------------
 719   THexa8b::THexa8b():
 720     TShapeFun(3,8)
 721   {
 722     TInt aNbRef = GetNbRef();
 723     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 724       TCoordSlice aCoord = GetCoord(aRefId);
 725       switch(aRefId){
 726       case  0: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
 727       case  3: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
 728       case  2: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
 729       case  1: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
 730       case  4: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
 731       case  7: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
 732       case  6: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
 733       case  5: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
 734       }
 735     }
 736   }
 737   //---------------------------------------------------------------
 738   THexa20b::THexa20b(TInt theDim, TInt theNbRef):
 739     TShapeFun(theDim,theNbRef)
 740   {
 741     std::size_t aNbRef = myRefCoord.size();
 742     for(std::size_t aRefId = 0; aRefId < aNbRef; aRefId++){
 743       TCoordSlice aCoord = GetCoord(aRefId);
 744       switch(aRefId){
 745       case  0: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
 746       case  3: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
 747       case  2: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
 748       case  1: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
 749       case  4: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
 750       case  7: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
 751       case  6: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
 752       case  5: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
 753
 754       case 11: aCoord[0] =  0.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
 755       case 10: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] = -1.0; break;
 756       case  9: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
 757       case  8: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] = -1.0; break;
 758       case 16: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] =  0.0; break;
 759       case 19: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] =  0.0; break;
 760       case 18: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
 761       case 17: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
 762       case 15: aCoord[0] =  0.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
 763       case 14: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
 764       case 13: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
 765       case 12: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
 766       }
 767     }
 768   }
 769   //---------------------------------------------------------------
 770   TPenta6a::TPenta6a():
 771     TShapeFun(3,6)
 772   {
 773     std::size_t aNbRef = myRefCoord.size();
 774     for(std::size_t aRefId = 0; aRefId < aNbRef; aRefId++){
 775       TCoordSlice aCoord = GetCoord(aRefId);
 776       switch(aRefId){
 777       case  0: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
 778       case  1: aCoord[0] = -1.0;  aCoord[1] = -0.0;  aCoord[2] =  1.0; break;
 779       case  2: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 780       case  3: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
 781       case  4: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
 782       case  5: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 783       }
 784     }
 785   }
 786   //---------------------------------------------------------------
 787   TPenta6b::TPenta6b():
 788     TShapeFun(3,6)
 789   {
 790     std::size_t aNbRef = myRefCoord.size();
 791     for(std::size_t aRefId = 0; aRefId < aNbRef; aRefId++){
 792       TCoordSlice aCoord = GetCoord(aRefId);
 793       switch(aRefId){
 794       case  0: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
 795       case  2: aCoord[0] = -1.0;  aCoord[1] = -0.0;  aCoord[2] =  1.0; break;
 796       case  1: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 797       case  3: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
 798       case  5: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
 799       case  4: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 800       }
 801     }
 802   }
 803   //---------------------------------------------------------------
 804   TPenta15a::TPenta15a():
 805     TShapeFun(3,15)
 806   {
 807     std::size_t aNbRef = myRefCoord.size();
 808     for(std::size_t aRefId = 0; aRefId < aNbRef; aRefId++){
 809       TCoordSlice aCoord = GetCoord(aRefId);
 810       switch(aRefId){
 811       case  0: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
 812       case  1: aCoord[0] = -1.0;  aCoord[1] = -0.0;  aCoord[2] =  1.0; break;
 813       case  2: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 814       case  3: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
 815       case  4: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
 816       case  5: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 817
 818       case  6: aCoord[0] = -1.0;  aCoord[1] =  0.5;  aCoord[2] =  0.5; break;
 819       case  7: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
 820       case  8: aCoord[0] = -1.0;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
 821       case  9: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
 822       case 10: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
 823       case 11: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 824       case 12: aCoord[0] =  1.0;  aCoord[1] =  0.5;  aCoord[2] =  0.5; break;
 825       case 13: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
 826       case 14: aCoord[0] =  1.0;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
 827       }
 828     }
 829   }
 830   //---------------------------------------------------------------
 831   TPenta15b::TPenta15b():
 832     TShapeFun(3,15)
 833   {
 834     std::size_t aNbRef = myRefCoord.size();
 835     for(std::size_t aRefId = 0; aRefId < aNbRef; aRefId++){
 836       TCoordSlice aCoord = GetCoord(aRefId);
 837       switch(aRefId){
 838       case  0: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
 839       case  2: aCoord[0] = -1.0;  aCoord[1] = -0.0;  aCoord[2] =  1.0; break;
 840       case  1: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 841       case  3: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
 842       case  5: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
 843       case  4: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 844
 845       case  8: aCoord[0] = -1.0;  aCoord[1] =  0.5;  aCoord[2] =  0.5; break;
 846       case  7: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
 847       case  6: aCoord[0] = -1.0;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
 848       case 12: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
 849       case 14: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
 850       case 13: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 851       case 11: aCoord[0] =  1.0;  aCoord[1] =  0.5;  aCoord[2] =  0.5; break;
 852       case 10: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
 853       case  9: aCoord[0] =  1.0;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
 854       }
 855     }
 856   }
 857   //---------------------------------------------------------------
 858   TPyra5a::TPyra5a():
 859     TShapeFun(3,5)
 860   {
 861     std::size_t aNbRef = myRefCoord.size();
 862     for(std::size_t aRefId = 0; aRefId < aNbRef; aRefId++){
 863       TCoordSlice aCoord = GetCoord(aRefId);
 864       switch(aRefId){
 865       case  0: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 866       case  1: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
 867       case  2: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 868       case  3: aCoord[0] =  0.0;  aCoord[1] = -1.0;  aCoord[2] =  0.0; break;
 869       case  4: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
 870       }
 871     }
 872   }
 873   //---------------------------------------------------------------
 874   TPyra5b::TPyra5b():
 875     TShapeFun(3,5)
 876   {
 877     std::size_t aNbRef = myRefCoord.size();
 878     for(std::size_t aRefId = 0; aRefId < aNbRef; aRefId++){
 879       TCoordSlice aCoord = GetCoord(aRefId);
 880       switch(aRefId){
 881       case  0: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 882       case  3: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
 883       case  2: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 884       case  1: aCoord[0] =  0.0;  aCoord[1] = -1.0;  aCoord[2] =  0.0; break;
 885       case  4: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
 886       }
 887     }
 888   }
 889   //---------------------------------------------------------------
 890   TPyra13a::TPyra13a():
 891     TShapeFun(3,13)
 892   {
 893     std::size_t aNbRef = myRefCoord.size();
 894     for(std::size_t aRefId = 0; aRefId < aNbRef; aRefId++){
 895       TCoordSlice aCoord = GetCoord(aRefId);
 896       switch(aRefId){
 897       case  0: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 898       case  1: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
 899       case  2: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 900       case  3: aCoord[0] =  0.0;  aCoord[1] = -1.0;  aCoord[2] =  0.0; break;
 901       case  4: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
 902
 903       case  5: aCoord[0] =  0.5;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
 904       case  6: aCoord[0] = -0.5;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
 905       case  7: aCoord[0] = -0.5;  aCoord[1] = -0.5;  aCoord[2] =  0.0; break;
 906       case  8: aCoord[0] =  0.5;  aCoord[1] = -0.5;  aCoord[2] =  0.0; break;
 907       case  9: aCoord[0] =  0.5;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
 908       case 10: aCoord[0] =  0.0;  aCoord[1] =  0.5;  aCoord[2] =  0.5; break;
 909       case 11: aCoord[0] = -0.5;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
 910       case 12: aCoord[0] =  0.0;  aCoord[1] = -0.5;  aCoord[2] =  0.5; break;
 911       }
 912     }
 913   }
 914   //---------------------------------------------------------------
 915   TPyra13b::TPyra13b():
 916     TShapeFun(3,13)
 917   {
 918     std::size_t aNbRef = myRefCoord.size();
 919     for(std::size_t aRefId = 0; aRefId < aNbRef; aRefId++){
 920       TCoordSlice aCoord = GetCoord(aRefId);
 921       switch(aRefId){
 922       case  0: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 923       case  3: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
 924       case  2: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 925       case  1: aCoord[0] =  0.0;  aCoord[1] = -1.0;  aCoord[2] =  0.0; break;
 926       case  4: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
 927
 928       case  8: aCoord[0] =  0.5;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
 929       case  7: aCoord[0] = -0.5;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
 930       case  6: aCoord[0] = -0.5;  aCoord[1] = -0.5;  aCoord[2] =  0.0; break;
 931       case  5: aCoord[0] =  0.5;  aCoord[1] = -0.5;  aCoord[2] =  0.0; break;
 932       case  9: aCoord[0] =  0.5;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
 933       case 12: aCoord[0] =  0.0;  aCoord[1] =  0.5;  aCoord[2] =  0.5; break;
 934       case 11: aCoord[0] = -0.5;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
 935       case 10: aCoord[0] =  0.0;  aCoord[1] = -0.5;  aCoord[2] =  0.5; break;
 936       }
 937     }
 938   }
 939   /*!
 940    * \brief Fill definition of gauss points family
 941    */
 942
 943   TGaussDef::TGaussDef(const int geom, const int nbGauss, const int variant)
 944   {
 945     myType = geom;
 946     myCoords .reserve( nbGauss * dim() );
 947     myWeights.reserve( nbGauss );
 948
 949     switch ( geom ) {
 950
 951     case NORM_SEG2:
 952     case NORM_SEG3:
 953       if (geom == NORM_SEG2) setRefCoords( TSeg2a() );
 954       else                   setRefCoords( TSeg3a() );
 955       switch ( nbGauss ) {
 956       case 1: {
 957         add( 0.0, 2.0 ); break;
 958       }
 959       case 2: {
 960         const double a = 0.577350269189626;
 961         add(  a,  1.0 );
 962         add( -a,  1.0 ); break;
 963       }
 964       case 3: {
 965         const double a = 0.774596669241;
 966         const double P1 = 1./1.8;
 967         const double P2 = 1./1.125;
 968         add( -a,  P1 );
 969         add(  0,  P2 );
 970         add(  a,  P1 ); break;
 971       }
 972       case 4: {
 973         const double a  = 0.339981043584856, b  = 0.861136311594053;
 974         const double P1 = 0.652145154862546, P2 = 0.347854845137454 ;
 975         add(  a,  P1 );
 976         add( -a,  P1 );
 977         add(  b,  P2 );
 978         add( -b,  P2 ); break;
 979       }
 980       default:
 981         THROW_IK_EXCEPTION("TGaussDef: Invalid nb of gauss points for SEG"<<nbGauss);
 982       }
 983       break;
 984
 985     case NORM_TRI3:
 986     case NORM_TRI6:
 987       if ( variant == 1 ) {
 988         if (geom == NORM_TRI3) setRefCoords( TTria3b() );
 989         else                   setRefCoords( TTria6b() );
 990         switch ( nbGauss ) {
 991         case 1: { // FPG1
 992           add( 1/3., 1/3., 1/2. ); break;
 993         }
 994         case 3: { // FPG3
 995           // what about COT3 ???
 996           add( 1/6., 1/6., 1/6. );
 997           add( 2/3., 1/6., 1/6. );
 998           add( 1/6., 2/3., 1/6. ); break;
 999         }
1000         case 4: { // FPG4
1001           add( 1/5., 1/5.,  25/(24*4.) );
1002           add( 3/5., 1/5.,  25/(24*4.) );
1003           add( 1/5., 3/5.,  25/(24*4.) );
1004           add( 1/3., 1/3., -27/(24*4.) ); break;
1005         }
1006         case 6: { // FPG6
1007           const double P1 = 0.11169079483905, P2 = 0.0549758718227661;
1008           const double a  = 0.445948490915965, b = 0.091576213509771;
1009           add(     b,     b, P2 );
1010           add( 1-2*b,     b, P2 );
1011           add(     b, 1-2*b, P2 );
1012           add(     a, 1-2*a, P1 );
1013           add(     a,     a, P1 );
1014           add( 1-2*a,     a, P1 ); break;
1015         }
1016         case 7: { // FPG7
1017           const double A  = 0.470142064105115;
1018           const double B  = 0.101286507323456;
1019           const double P1 = 0.066197076394253;
1020           const double P2 = 0.062969590272413;
1021           add(  1/3.,  1/3., 9/80. );
1022           add(     A,     A, P1 );
1023           add( 1-2*A,     A, P1 );
1024           add(     A, 1-2*A, P1 );
1025           add(     B,     B, P2 );
1026           add( 1-2*B,     B, P2 );
1027           add(     B, 1-2*B, P2 ); break;
1028         }
1029         case 12: { // FPG12
1030           const double A  = 0.063089014491502;
1031           const double B  = 0.249286745170910;
1032           const double C  = 0.310352451033785;
1033           const double D  = 0.053145049844816;
1034           const double P1 = 0.025422453185103;
1035           const double P2 = 0.058393137863189;
1036           const double P3 = 0.041425537809187;
1037           add(     A,     A, P1 );
1038           add( 1-2*A,     A, P1 );
1039           add(     A, 1-2*A, P1 );
1040           add(     B,     B, P2 );
1041           add( 1-2*B,     B, P2 );
1042           add(     B, 1-2*B, P2 );
1043           add(     C,     D, P3 );
1044           add(     D,     C, P3 );
1045           add( 1-C-D,     C, P3 );
1046           add( 1-C-D,     D, P3 );
1047           add(     C, 1-C-D, P3 );
1048           add(     D, 1-C-D, P3 ); break;
1049         }
1050         default:
1051           THROW_IK_EXCEPTION("TGaussDef: Invalid nb of gauss points for TRIA, variant 1: "
1052                              <<nbGauss);
1053         }
1054       }
1055       else if ( variant == 2 ) {
1056         if (geom == NORM_TRI3) setRefCoords( TTria3a() );
1057         else                   setRefCoords( TTria6a() );
1058         switch ( nbGauss ) {
1059         case 1: {
1060           add( -1/3., -1/3., 2. ); break;
1061         }
1062         case 3: {
1063           add( -2/3.,  1/3., 2/3. );
1064           add( -2/3., -2/3., 2/3. );
1065           add(  1/3., -2/3., 2/3. ); break;
1066         }
1067         case 6: {
1068           const double P1 = 0.11169079483905, P2 = 0.0549758718227661;
1069           const double A  = 0.445948490915965, B = 0.091576213509771;
1070           add( 2*B-1, 1-4*B, 4*P2 );
1071           add( 2*B-1, 2*B-1, 4*P2 );
1072           add( 1-4*B, 2*B-1, 4*P2 );
1073           add( 1-4*A, 2*A-1, 4*P1 );
1074           add( 2*A-1, 1-4*A, 4*P1 );
1075           add( 2*A-1, 2*A-1, 4*P1 ); break;
1076         }
1077         default:
1078           THROW_IK_EXCEPTION("TGaussDef: Invalid nb of gauss points for TRIA, variant 2: "
1079                              <<nbGauss);
1080         }
1081       }
1082       else if ( variant == 3 ) {
1083         if (geom == NORM_TRI3) setRefCoords( TTria3b() );
1084         else                   setRefCoords( TTria6b() );
1085         switch ( nbGauss ) {
1086         case 4: {
1087           add( 1/3., 1/3., -27/96 );
1088           add( 0.2 , 0.2 ,  25/96 );
1089           add( 0.6 , 0.2 ,  25/96 );
1090           add( 0.2 , 0.6 ,  25/96 ); break;
1091         }
1092         default:
1093           THROW_IK_EXCEPTION("TGaussDef: Invalid nb of gauss points for TRIA, variant 3: "
1094                              <<nbGauss);
1095         }
1096       }
1097       break;
1098
1099     case NORM_QUAD4:
1100     case NORM_QUAD8:
1101       if ( variant == 1 ) {
1102         if (geom == NORM_QUAD4) setRefCoords( TQuad4b() );
1103         else                    setRefCoords( TQuad8b() );
1104         switch ( nbGauss ) {
1105         case 1: { // FPG1
1106           add(  0,  0,  4 ); break;
1107         }
1108         case 4: { // FPG4
1109           const double a = 1/sqrt(3.);
1110           add( -a, -a,  1 );
1111           add(  a, -a,  1 );
1112           add(  a,  a,  1 );
1113           add( -a,  a,  1 ); break;
1114         }
1115         case 5: { // out from the 3 specs
1116           const double a = 0.774596669241483;
1117           add( -a, -a,  0.5 );
1118           add(  a, -a,  0.5 );
1119           add(  a,  a,  0.5 );
1120           add( -a,  a,  0.5 );
1121           add(  0,  0,  2.0 ); break;
1122         }
1123         case 9: { // FPG9
1124           const double a = 0.774596669241483;
1125           add( -a, -a,  25/81. );
1126           add(  a, -a,  25/81. );
1127           add(  a,  a,  25/81. );
1128           add( -a,  a,  25/81. );
1129           add( 0., -a,  40/81. );
1130           add(  a, 0.,  40/81. );
1131           add( 0.,  a,  40/81. );
1132           add( -a, 0.,  40/81. );
1133           add( 0., 0.,  64/81. ); break;
1134         }
1135         default:
1136           THROW_IK_EXCEPTION("TGaussDef: Invalid nb of gauss points for QUAD, variant 1: "
1137                              <<nbGauss);
1138         }
1139       }
1140       else if ( variant == 2 ) {
1141         if (geom == NORM_QUAD4) setRefCoords( TQuad4a() );
1142         else                    setRefCoords( TQuad8a() );
1143         switch ( nbGauss ) {
1144         case 4: {
1145           const double a = 1/sqrt(3.);
1146           add( -a,  a,  1 );
1147           add( -a, -a,  1 );
1148           add(  a, -a,  1 );
1149           add(  a,  a,  1 ); break;
1150         }
1151         case 9: {
1152           const double a = 0.774596669241483;
1153           add( -a,  a,  25/81. );
1154           add( -a, -a,  25/81. );
1155           add(  a, -a,  25/81. );
1156           add(  a,  a,  25/81. );
1157           add( -a, 0.,  40/81. );
1158           add( 0., -a,  40/81. );
1159           add(  a, 0.,  40/81. );
1160           add( 0.,  a,  40/81. );
1161           add( 0., 0.,  64/81. ); break;
1162         }
1163         default:
1164           THROW_IK_EXCEPTION("TGaussDef: Invalid nb of gauss points for QUAD, variant 1: "
1165                              <<nbGauss);
1166         }
1167       }
1168       else if ( variant == 3 ) {
1169         if (geom == NORM_QUAD4) setRefCoords( TQuad4b() );
1170         else                    setRefCoords( TQuad8b() );
1171         switch ( nbGauss ) {
1172         case 4: {
1173           const double a = 3/sqrt(3.);
1174           add( -a, -a,  1 );
1175           add( -a,  a,  1 );
1176           add(  a, -a,  1 );
1177           add(  a,  a,  1 ); break;
1178         }
1179         case 9: {
1180           const double a = sqrt(3/5.), c1 = 5/9., c2 = 8/9.;
1181           const double c12 = c1*c2, c22 = c2*c2, c1c2 = c1*c2;
1182           add( -a, -a,  c12  );
1183           add( -a, 0.,  c1c2 );
1184           add( -a,  a,  c12  );
1185           add( 0., -a,  c1c2 );
1186           add( 0., 0.,  c22  );
1187           add( 0.,  a,  c1c2 );
1188           add(  a, -a,  c12  );
1189           add(  a, 0.,  c1c2 );
1190           add(  a,  a,  c12  ); break;
1191         }
1192         default:
1193           THROW_IK_EXCEPTION("TGaussDef: Invalid nb of gauss points for QUAD, variant 3: "
1194                              <<nbGauss);
1195         }
1196       }
1197       break;
1198
1199     case NORM_TETRA4:
1200     case NORM_TETRA10:
1201       if (geom == NORM_TETRA4) setRefCoords( TTetra4a() );
1202       else                 setRefCoords( TTetra10a() );
1203       switch ( nbGauss ) {
1204       case 4: { // FPG4
1205         const double a = (5 - sqrt(5.))/20., b = (5 + 3*sqrt(5.))/20.;
1206         add(  a,  a,  a,  1/24. );
1207         add(  a,  a,  b,  1/24. );
1208         add(  a,  b,  a,  1/24. );
1209         add(  b,  a,  a,  1/24. ); break;
1210       }
1211       case 5: { // FPG5
1212         const double a = 0.25, b = 1/6., c = 0.5;
1213         add(  a,  a,  a, -2/15. );
1214         add(  b,  b,  b,  3/40. );
1215         add(  b,  b,  c,  3/40. );
1216         add(  b,  c,  b,  3/40. );
1217         add(  c,  b,  b,  3/40. ); break;
1218       }
1219       case 15: { // FPG15
1220         const double a = 0.25;
1221         const double b1 = (7 + sqrt(15.))/34., c1 = (13 + 3*sqrt(15.))/34., d = (5 - sqrt(15.))/20.;
1222         const double b2 = (7 - sqrt(15.))/34., c2 = (13 - 3*sqrt(15.))/34., e = (5 + sqrt(15.))/20.;
1223         const double P1 = (2665 - 14*sqrt(15.))/226800.;
1224         const double P2 = (2665 + 14*sqrt(15.))/226800.;
1225         add(  a,  a,  a,  8/405.);//_____
1226         add( b1, b1, b1,  P1    );
1227         add( b1, b1, c1,  P1    );
1228         add( b1, c1, b1,  P1    );
1229         add( c1, b1, b1,  P1    );//_____
1230         add( b2, b2, b2,  P2    );
1231         add( b2, b2, c2,  P2    );
1232         add( b2, c2, b2,  P2    );
1233         add( c2, b2, b2,  P2    );//_____
1234         add(  d,  d,  e,  5/567.);
1235         add(  d,  e,  d,  5/567.);
1236         add(  e,  d,  d,  5/567.);
1237         add(  d,  e,  e,  5/567.);
1238         add(  e,  d,  e,  5/567.);
1239         add(  e,  e,  d,  5/567.);
1240         break;
1241       }
1242       default:
1243         THROW_IK_EXCEPTION("TGaussDef: Invalid nb of gauss points for TETRA: "<<nbGauss);
1244       }
1245       break;
1246
1247     case NORM_PYRA5:
1248     case NORM_PYRA13:
1249       if (geom == NORM_PYRA5) setRefCoords( TPyra5a() );
1250       else                setRefCoords( TPyra13a() );
1251       switch ( nbGauss ) {
1252       case 5: { // FPG5
1253         const double h1 = 0.1531754163448146;
1254         const double h2 = 0.6372983346207416;
1255         add(  .5,  0.,  h1,  2/15. );
1256         add(  0.,  .5,  h1,  2/15. );
1257         add( -.5,  0.,  h1,  2/15. );
1258         add(  0., -.5,  h1,  2/15. );
1259         add(  0.,  0.,  h2,  2/15. ); break;
1260       }
1261       case 6: { // FPG6
1262         const double p1 = 0.1024890634400000 ;
1263         const double p2 = 0.1100000000000000 ;
1264         const double p3 = 0.1467104129066667 ;
1265         const double a  = 0.5702963741068025 ;
1266         const double h1 = 0.1666666666666666 ;
1267         const double h2 = 0.08063183038464675;
1268         const double h3 = 0.6098484849057127 ;
1269         add(  a, 0.,  h1,  p1 );
1270         add( 0.,  a,  h1,  p1 );
1271         add( -a, 0.,  h1,  p1 );
1272         add( 0., -a,  h1,  p1 );
1273         add( 0., 0.,  h2,  p2 );
1274         add( 0., 0.,  h3,  p3 ); break;
1275       }
1276       case 27: { // FPG27
1277         const double a1  = 0.788073483;
1278         const double b6  = 0.499369002;
1279         const double b1  = 0.848418011;
1280         const double c8  = 0.478508449;
1281         const double c1  = 0.652816472;
1282         const double d12 = 0.032303742;
1283         const double d1  = 1.106412899;
1284         double z = 1/2., fz = b1/2*(1 - z);
1285         add(  0.,  0.,   z,  a1 ); // 1
1286         add(  fz,  fz,   z,  b6 ); // 2
1287         add( -fz,  fz,   z,  b6 ); // 3
1288         add( -fz, -fz,   z,  b6 ); // 4
1289         add(  fz, -fz,   z,  b6 ); // 5
1290         z = (1 - b1)/2.;
1291         add(  0.,  0.,   z,  b6 ); // 6
1292         z = (1 + b1)/2.;
1293         add(  0.,  0.,   z,  b6 ); // 7
1294         z = (1 - c1)/2.; fz = c1*(1 - z);
1295         add(  fz,  0.,   z,  c8 ); // 8
1296         add(  0.,  fz,   z,  c8 ); // 9
1297         add( -fz,  0.,   z,  c8 ); // 10
1298         add(  0., -fz,   z,  c8 ); // 11
1299         z = (1 + c1)/2.; fz = c1*(1 - z);
1300         add(  fz,  0.,   z,  c8 ); // 12
1301         add(  0.,  fz,   z,  c8 ); // 13
1302         add( -fz,  0.,   z,  c8 ); // 14
1303         add(  0., -fz,   z,  c8 ); // 15
1304         z = (1 - d1)/2., fz = d1/2*(1 - z);
1305         add(  fz,  fz,   z,  d12); // 16
1306         add( -fz,  fz,   z,  d12); // 17
1307         add( -fz, -fz,   z,  d12); // 18
1308         add(  fz, -fz,   z,  d12); // 19
1309         z = 1/2.; fz = d1*(1 - z);
1310         add(  fz,  0.,   z,  d12); // 20
1311         add(  0.,  fz,   z,  d12); // 21
1312         add( -fz,  0.,   z,  d12); // 22
1313         add(  0., -fz,   z,  d12); // 23
1314         z = (1 + d1)/2., fz = d1/2*(1 - z);
1315         add(  fz,  fz,   z,  d12); // 24
1316         add( -fz,  fz,   z,  d12); // 25
1317         add( -fz, -fz,   z,  d12); // 26
1318         add(  fz, -fz,   z,  d12); // 27
1319         break;
1320       }
1321       default:
1322         THROW_IK_EXCEPTION("TGaussDef: Invalid nb of gauss points for PYRA: "<<nbGauss);
1323       }
1324       break;
1325     case NORM_PENTA6:
1326     case NORM_PENTA15:
1327       if (geom == NORM_PENTA6) setRefCoords( TPenta6a() );
1328       else                 setRefCoords( TPenta15a() );
1329       switch ( nbGauss ) {
1330       case 6: { // FPG6
1331         const double a = sqrt(3.)/3.;
1332         add( -a, .5, .5,  1/6. );
1333         add( -a, 0., .5,  1/6. );
1334         add( -a, .5, 0.,  1/6. );
1335         add(  a, .5, .5,  1/6. );
1336         add(  a, 0., .5,  1/6. );
1337         add(  a, .5, 0.,  1/6. ); break;
1338       }
1339       case 8: { // FPG8
1340         const double a = 0.577350269189626;
1341         add( -a, 1/3., 1/3., -27/96. );
1342         add( -a,  0.6,  0.2,  25/96. );
1343         add( -a,  0.2,  0.6,  25/96. );
1344         add( -a,  0.2,  0.2,  25/96. );
1345         add( +a, 1/3., 1/3., -27/96. );
1346         add( +a,  0.6,  0.2,  25/96. );
1347         add( +a,  0.2,  0.6,  25/96. );
1348         add( +a,  0.2,  0.2,  25/96. ); break;
1349       }
1350       case 21: { // FPG21
1351         const double d = sqrt(3/5.), c1 = 5/9., c2 = 8/9.; // d <=> alfa
1352         const double a = (6 + sqrt(15.))/21.;
1353         const double b = (6 - sqrt(15.))/21.;
1354         const double P1 = (155 + sqrt(15.))/2400.;
1355         const double P2 = (155 - sqrt(15.))/2400.;  //___
1356         add( -d,  1/3.,  1/3., c1*9/80. );//___
1357         add( -d,     a,     a, c1*P1    );
1358         add( -d, 1-2*a,     a, c1*P1    );
1359         add( -d,     a, 1-2*a, c1*P1    );//___
1360         add( -d,     b,     b, c1*P2    );
1361         add( -d, 1-2*b,     b, c1*P2    );
1362         add( -d,     b, 1-2*b, c1*P2    );//___
1363         add( 0.,  1/3.,  1/3., c2*9/80. );//___
1364         add( 0.,     a,     a, c2*P1    );
1365         add( 0., 1-2*a,     a, c2*P1    );
1366         add( 0.,     a, 1-2*a, c2*P1    );//___
1367         add( 0.,     b,     b, c2*P2    );
1368         add( 0., 1-2*b,     b, c2*P2    );
1369         add( 0.,     b, 1-2*b, c2*P2    );//___
1370         add(  d,  1/3.,  1/3., c1*9/80. );//___
1371         add(  d,     a,     a, c1*P1    );
1372         add(  d, 1-2*a,     a, c1*P1    );
1373         add(  d,     a, 1-2*a, c1*P1    );//___
1374         add(  d,     b,     b, c1*P2    );
1375         add(  d, 1-2*b,     b, c1*P2    );
1376         add(  d,     b, 1-2*b, c1*P2    );//___
1377         break;
1378       }
1379       default:
1380         THROW_IK_EXCEPTION("TGaussDef: Invalid nb of gauss points for PENTA: " <<nbGauss);
1381       }
1382       break;
1383
1384     case NORM_HEXA8:
1385     case NORM_HEXA20:
1386       if (geom == NORM_HEXA8) setRefCoords( THexa8a() );
1387       else                    setRefCoords( THexa20a() );
1388       switch ( nbGauss ) {
1389       case 8: { // FPG8
1390         const double a = sqrt(3.)/3.;
1391         add( -a, -a, -a,  1. );
1392         add( -a, -a,  a,  1. );
1393         add( -a,  a, -a,  1. );
1394         add( -a,  a,  a,  1. );
1395         add(  a, -a, -a,  1. );
1396         add(  a, -a,  a,  1. );
1397         add(  a,  a, -a,  1. );
1398         add(  a,  a,  a,  1. ); break;
1399       }
1400       case 27: { // FPG27
1401         const double a = sqrt(3/5.), c1 = 5/9., c2 = 8/9.;
1402         const double c12 = c1*c1, c13 = c1*c1*c1;
1403         const double c22 = c2*c2, c23 = c2*c2*c2;
1404         add( -a, -a, -a,   c13  ); // 1
1405         add( -a, -a, 0., c12*c2 ); // 2
1406         add( -a, -a,  a,   c13  ); // 3
1407         add( -a, 0., -a, c12*c2 ); // 4
1408         add( -a, 0., 0., c1*c22 ); // 5
1409         add( -a, 0.,  a, c12*c2 ); // 6
1410         add( -a,  a, -a,   c13  ); // 7
1411         add( -a,  a, 0., c12*c2 ); // 8
1412         add( -a,  a,  a,   c13  ); // 9
1413         add( 0., -a, -a, c12*c2 ); // 10
1414         add( 0., -a, 0., c1*c22 ); // 11
1415         add( 0., -a,  a, c12*c2 ); // 12
1416         add( 0., 0., -a, c1*c22 ); // 13
1417         add( 0., 0., 0.,   c23  ); // 14
1418         add( 0., 0.,  a, c1*c22 ); // 15
1419         add( 0.,  a, -a, c12*c2 ); // 16
1420         add( 0.,  a, 0., c1*c22 ); // 17
1421         add( 0.,  a,  a, c12*c2 ); // 18
1422         add(  a, -a, -a,   c13  ); // 19
1423         add(  a, -a, 0., c12*c2 ); // 20
1424         add(  a, -a,  a,   c13  ); // 21
1425         add(  a, 0., -a, c12*c2 ); // 22
1426         add(  a, 0., 0., c1*c22 ); // 23
1427         add(  a, 0.,  a, c12*c2 ); // 24
1428         add(  a,  a, -a,   c13  ); // 25
1429         add(  a,  a, 0., c12*c2 ); // 26
1430         add(  a,  a,  a,   c13  ); // 27
1431         break;
1432       }
1433       default:
1434         THROW_IK_EXCEPTION("TGaussDef: Invalid nb of gauss points for PENTA: " <<nbGauss);
1435       }
1436       break;
1437
1438     default:
1439       THROW_IK_EXCEPTION("TGaussDef: unexpected EGeometrieElement: "<< geom);
1440     }
1441
1442     if ( myWeights.capacity() != myWeights.size() )
1443       THROW_IK_EXCEPTION("TGaussDef: Not all gauss points defined");
1444   }
1445 }
1446
1447 //================================================================================
1448 /*!
1449  * \brief Return dimension for the given cell type
1450  */
1451 //================================================================================
1452
1453 unsigned SauvUtilities::getDimension( INTERP_KERNEL::NormalizedCellType type )
1454 {
1455   return type == NORM_ERROR ? -1 : INTERP_KERNEL::CellModel::GetCellModel( type ).getDimension();
1456 }
1457
1458 //================================================================================
1459 /*!
1460  * \brief Returns interlace array to transform a quadratic GIBI element to a MED one.
1461  *        i-th array item gives node index in GIBI connectivity for i-th MED node
1462  */
1463 //================================================================================
1464
1465 const int * SauvUtilities::getGibi2MedQuadraticInterlace( INTERP_KERNEL::NormalizedCellType type )
1466 {
1467   static std::vector<const int*> conn;
1468   static const int hexa20 [] = {0,6,4,2, 12,18,16,14, 7,5,3,1, 19,17,15,13, 8,11,10,9};
1469   static const int penta15[] = {0,2,4, 9,11,13, 1,3,5, 10,12,14, 6,8,7};
1470   static const int pyra13 [] = {0,2,4,6, 12, 1,3,5,7, 8,9,10,11};
1471   static const int tetra10[] = {0,2,4, 9, 1,3,5, 6,7,8};
1472   static const int quad8  [] = {0,2,4,6, 1,3,5,7};
1473   static const int tria6  [] = {0,2,4, 1,3,5};
1474   static const int seg3   [] = {0,2,1};
1475   if ( conn.empty() )
1476     {
1477       conn.resize( MaxMedCellType + 1, 0 );
1478       conn[ NORM_HEXA20 ] = hexa20;
1479       conn[ NORM_PENTA15] = penta15;
1480       conn[ NORM_PYRA13 ] = pyra13;
1481       conn[ NORM_TETRA10] = tetra10;
1482       conn[ NORM_SEG3   ] = seg3;
1483       conn[ NORM_TRI6   ] = tria6;
1484       conn[ NORM_QUAD8  ] = quad8;
1485     }
1486   return conn[ type ];
1487 }
1488
1489 //================================================================================
1490 /*!
1491  * \brief Avoid coping sortedNodeIDs
1492  */
1493 //================================================================================
1494
1495 Cell::Cell(const Cell& ma)
1496   : _nodes(ma._nodes), _reverse(ma._reverse), _sortedNodeIDs(0), _number(ma._number)
1497 {
1498   if ( ma._sortedNodeIDs )
1499     {
1500       _sortedNodeIDs = new TID[ _nodes.size() ];
1501       std::copy( ma._sortedNodeIDs, ma._sortedNodeIDs + _nodes.size(), _sortedNodeIDs );
1502     }
1503 }
1504
1505 //================================================================================
1506 /*!
1507  * \brief Rerturn the i-th link of face
1508  */
1509 //================================================================================
1510
1511 SauvUtilities::Link Cell::link(int i) const
1512 {
1513   std::size_t i2 = ( i + 1 ) % _nodes.size();
1514   if ( _reverse )
1515     return std::make_pair( _nodes[i2]->_number, _nodes[i]->_number );
1516   else
1517     return std::make_pair( _nodes[i]->_number, _nodes[i2]->_number );
1518 }
1519
1520 //================================================================================
1521 /*!
1522  * \brief creates if needed and return _sortedNodeIDs
1523  */
1524 //================================================================================
1525
1526 const TID* Cell::getSortedNodes() const
1527 {
1528   if ( !_sortedNodeIDs )
1529     {
1530       size_t l=_nodes.size();
1531       _sortedNodeIDs = new TID[ l ];
1532
1533       for (size_t i=0; i!=l; ++i)
1534         _sortedNodeIDs[i]=_nodes[i]->_number;
1535       std::sort( _sortedNodeIDs, _sortedNodeIDs + l );
1536     }
1537   return _sortedNodeIDs;
1538 }
1539
1540 //================================================================================
1541 /*!
1542  * \brief Compare sorted ids of cell nodes
1543  */
1544 //================================================================================
1545
1546 bool Cell::operator< (const Cell& ma) const
1547 {
1548   if ( _nodes.size() == 1 )
1549     return _nodes[0] < ma._nodes[0];
1550
1551   const TID* v1 = getSortedNodes();
1552   const TID* v2 = ma.getSortedNodes();
1553   for ( const TID* vEnd = v1 + _nodes.size(); v1 < vEnd; ++v1, ++v2 )
1554     if(*v1 != *v2)
1555       return *v1 < *v2;
1556   return false;
1557 }
1558
1559 //================================================================================
1560 /*!
1561  * \brief Dump a Cell
1562  */
1563 //================================================================================
1564
1565 std::ostream& SauvUtilities::operator<< (std::ostream& os, const SauvUtilities::Cell& ma)
1566 {
1567   os << "cell " << ma._number << " (" << ma._nodes.size() << " nodes) : < " << ma._nodes[0]->_number;
1568   for( size_t i=1; i!=ma._nodes.size(); ++i)
1569     os << ", " << ma._nodes[i]->_number;
1570 #ifdef _DEBUG_
1571   os << " > sortedNodes: ";
1572   if ( ma._sortedNodeIDs ) {
1573     os << "< ";
1574     for( size_t i=0; i!=ma._nodes.size(); ++i)
1575       os << ( i ? ", " : "" ) << ma._sortedNodeIDs[i];
1576     os << " >";
1577   }
1578   else {
1579     os << "NULL";
1580   }
1581 #endif
1582   return os;
1583 }
1584
1585 //================================================================================
1586 /*!
1587  * \brief Return nb of elements in the group
1588  */
1589 //================================================================================
1590
1591 mcIdType Group::size() const
1592 {
1593   std::size_t sizze = 0;
1594   if ( !_relocTable.empty() )
1595     sizze =  _relocTable.size();
1596   else if ( _medGroup )
1597     sizze = _medGroup->getNumberOfTuples();
1598   else if ( !_cells.empty() )
1599     sizze = _cells.size();
1600   else
1601     for ( size_t i = 0; i < _groups.size(); ++i )
1602       sizze += _groups[i]->size();
1603   return ToIdType( sizze );
1604 }
1605
1606 //================================================================================
1607 /*!
1608  * \brief Conver gibi element type to med one
1609  */
1610 //================================================================================
1611
1612 INTERP_KERNEL::NormalizedCellType SauvUtilities::gibi2medGeom( size_t gibiType )
1613 {
1614   if ( gibiType < 1 || gibiType > NbGibiCellTypes )
1615     return NORM_ERROR;
1616
1617   return GibiTypeToMed[ gibiType - 1 ];
1618 }
1619
1620 //================================================================================
1621 /*!
1622  * \brief Conver med element type to gibi one
1623  */
1624 //================================================================================
1625
1626 int SauvUtilities::med2gibiGeom( INTERP_KERNEL::NormalizedCellType medGeomType )
1627 {
1628   for ( unsigned int i = 0; i < NbGibiCellTypes; i++ )
1629     if ( GibiTypeToMed[ i ] == medGeomType )
1630       return i + 1;
1631
1632   return -1;
1633 }
1634
1635 //================================================================================
1636 /*!
1637  * \brief Remember the file name
1638  */
1639 //================================================================================
1640
1641 FileReader::FileReader(const char* fileName):_fileName(fileName),_iRead(0),_nbToRead(0)
1642 {
1643 }
1644
1645 //================================================================================
1646 /*!
1647  * \brief Constructor of ASCII sauve file reader
1648  */
1649 //================================================================================
1650
1651 ASCIIReader::ASCIIReader(const char* fileName)
1652   :FileReader(fileName),
1653    _file(-1)
1654 {
1655 }
1656
1657 //================================================================================
1658 /*!
1659  * \brief Return true
1660  */
1661 //================================================================================
1662
1663 bool ASCIIReader::isASCII() const
1664 {
1665   return true;
1666 }
1667
1668 //================================================================================
1669 /*!
1670  * \brief Try to open an ASCII file
1671  */
1672 //================================================================================
1673
1674 bool ASCIIReader::open()
1675 {
1676 #ifdef WIN32
1677   _file = ::_open (_fileName.c_str(), _O_RDONLY|_O_BINARY);
1678 #else
1679   _file = ::open (_fileName.c_str(), O_RDONLY);
1680 #endif
1681   if (_file >= 0)
1682     {
1683       _start  = new char [GIBI_BufferSize]; // working buffer beginning
1684       //_tmpBuf = new char [GIBI_MaxOutputLen];
1685       _ptr    = _start;
1686       _eptr   = _start;
1687       _lineNb = 0;
1688     }
1689   else
1690     {
1691       //THROW_IK_EXCEPTION("Can't open file "<<_fileName << " fd: " << _file);
1692     }
1693   return (_file >= 0);
1694 }
1695
1696 //================================================================================
1697 /*!
1698  * \brief Close the file
1699  */
1700 //================================================================================
1701
1702 ASCIIReader::~ASCIIReader()
1703 {
1704   if (_file >= 0)
1705     {
1706       ::close (_file);
1707       if (_start != 0L)
1708         {
1709           delete [] _start;
1710           //delete [] _tmpBuf;
1711           _start = 0;
1712         }
1713       _file = -1;
1714     }
1715 }
1716
1717 //================================================================================
1718 /*!
1719  * \brief Return a next line of the file
1720  */
1721 //================================================================================
1722
1723 bool ASCIIReader::getNextLine (char* & line, bool raiseOEF /*= true*/ )
1724 {
1725   if ( getLine( line )) return true;
1726   if ( raiseOEF )
1727     THROW_IK_EXCEPTION("Unexpected EOF on ln "<<_lineNb);
1728   return false;
1729 }
1730
1731 //================================================================================
1732 /*!
1733  * \brief Read a next line of the file if necessary
1734  */
1735 //================================================================================
1736
1737 bool ASCIIReader::getLine(char* & line)
1738 {
1739   bool aResult = true;
1740   // Check the state of the buffer;
1741   // if there is too little left, read the next portion of data
1742   std::size_t nBytesRest = _eptr - _ptr;
1743   if (nBytesRest < GIBI_MaxOutputLen)
1744     {
1745       if (nBytesRest > 0)
1746         {
1747           // move the remaining portion to the buffer beginning
1748           for ( std::size_t i = 0; i < nBytesRest; ++i )
1749             _start[i] = _ptr[i];
1750           //memcpy (_tmpBuf, _ptr, nBytesRest);
1751           //memcpy (_start, _tmpBuf, nBytesRest);
1752         }
1753       else
1754         {
1755           nBytesRest = 0;
1756         }
1757       _ptr = _start;
1758       const std::size_t nBytesRead = ::read (_file,
1759                                              &_start [nBytesRest],
1760                                              GIBI_BufferSize - nBytesRest);
1761       nBytesRest += nBytesRead;
1762       _eptr = &_start [nBytesRest];
1763     }
1764   // Check the buffer for the end-of-line
1765   char * ptr = _ptr;
1766   while (true)
1767     {
1768       // Check for end-of-the-buffer, the ultimate criterion for termination
1769       if (ptr >= _eptr)
1770         {
1771           if (nBytesRest <= 0)
1772             aResult = false;
1773           else
1774             _eptr[-1] = '\0';
1775           break;
1776         }
1777       // seek the line-feed character
1778       if (ptr[0] == '\n')
1779         {
1780           if (ptr[-1] == '\r')
1781             ptr[-1] = '\0';
1782           ptr[0] = '\0';
1783           ++ptr;
1784           break;
1785         }
1786       ++ptr;
1787     }
1788   // Output the result
1789   line = _ptr;
1790   _ptr = ptr;
1791   _lineNb++;
1792
1793   return aResult;
1794 }
1795
1796 //================================================================================
1797 /*!
1798  * \brief Prepare for iterating over given nb of values
1799  *  \param nbToRead - nb of fields to read
1800  *  \param nbPosInLine - nb of fields in one line
1801  *  \param width - field width
1802  *  \param shift - shift from line beginning to the field start
1803  */
1804 //================================================================================
1805
1806 void ASCIIReader::init( int nbToRead, int nbPosInLine, int width, int shift /*= 0*/ )
1807 {
1808   _nbToRead    = nbToRead;
1809   _nbPosInLine = nbPosInLine;
1810   _width       = width;
1811   _shift       = shift;
1812   _iPos = _iRead = 0;
1813   if ( _nbToRead )
1814     {
1815       getNextLine( _curPos );
1816       _curPos = _curPos + _shift;
1817     }
1818   else
1819     {
1820       _curPos = 0;
1821     }
1822 }
1823
1824 //================================================================================
1825 /*!
1826  * \brief Prepare for iterating over given nb of string values
1827  */
1828 //================================================================================
1829
1830 void ASCIIReader::initNameReading(int nbValues, int width /*= 8*/)
1831 {
1832   init( nbValues, 72 / ( width + 1 ), width, 1 );
1833 }
1834
1835 //================================================================================
1836 /*!
1837  * \brief Prepare for iterating over given nb of integer values
1838  */
1839 //================================================================================
1840
1841 void ASCIIReader::initIntReading(int nbValues)
1842 {
1843   init( nbValues, 10, 8 );
1844 }
1845
1846 //================================================================================
1847 /*!
1848  * \brief Prepare for iterating over given nb of real values
1849  */
1850 //================================================================================
1851
1852 void ASCIIReader::initDoubleReading(int nbValues)
1853 {
1854   init( nbValues, 3, 22 );
1855 }
1856
1857 //================================================================================
1858 /*!
1859  * \brief Return true if not all values have been read
1860  */
1861 //================================================================================
1862
1863 bool ASCIIReader::more() const
1864 {
1865   bool result = false;
1866   if ( _iRead < _nbToRead)
1867     {
1868       if ( _curPos ) result = true;
1869     }
1870   return result;
1871 }
1872
1873 //================================================================================
1874 /*!
1875  * \brief Go to the nex value
1876  */
1877 //================================================================================
1878
1879 void ASCIIReader::next()
1880 {
1881   if ( !more() )
1882     THROW_IK_EXCEPTION("SauvUtilities::ASCIIReader::next(): no more() values to read");
1883   ++_iRead;
1884   ++_iPos;
1885   if ( _iRead < _nbToRead )
1886     {
1887       if ( _iPos >= _nbPosInLine )
1888         {
1889           getNextLine( _curPos );
1890           _curPos = _curPos + _shift;
1891           _iPos = 0;
1892         }
1893       else
1894         {
1895           _curPos = _curPos + _width + _shift;
1896         }
1897     }
1898   else
1899     {
1900       _curPos = 0;
1901     }
1902 }
1903
1904 //================================================================================
1905 /*!
1906  * \brief Return the current integer value
1907  */
1908 //================================================================================
1909
1910 int ASCIIReader::getInt() const
1911 {
1912   // fix for two glued ints (issue 0021009):
1913   // Line nb    |   File contents
1914   // ------------------------------------------------------------------------------------
1915   // 53619905   |       1       2       6       8
1916   // 53619906   |                                                                SCALAIRE
1917   // 53619907   |    -63312600499       1       0       0       0      -2       0       2
1918   //   where -63312600499 is actually -633 and 12600499
1919   char hold=_curPos[_width];
1920   _curPos[_width] = '\0';
1921   int result = atoi( _curPos );
1922   _curPos[_width] = hold;
1923   return result;
1924   //return atoi(str());
1925 }
1926
1927 //================================================================================
1928 /*!
1929  * \brief Return the current float value
1930  */
1931 //================================================================================
1932
1933 float ASCIIReader::getFloat() const
1934 {
1935   return (float)getDouble();
1936 }
1937
1938 //================================================================================
1939 /*!
1940  * \brief Return the current double value
1941  */
1942 //================================================================================
1943
1944 double ASCIIReader::getDouble() const
1945 {
1946   //std::string aStr (_curPos);
1947
1948   // Correction: add missing 'E' specifier
1949   // int aPosStart = aStr.find_first_not_of(" \t");
1950   // if (aPosStart < (int)aStr.length()) {
1951   //   int aPosSign = aStr.find_first_of("+-", aPosStart + 1); // pass the leading symbol, as it can be a sign
1952   //   if (aPosSign < (int)aStr.length()) {
1953   //     if (aStr[aPosSign - 1] != 'e' && aStr[aPosSign - 1] != 'E')
1954   //       aStr.insert(aPosSign, "E", 1);
1955   //   }
1956   // }
1957
1958   // Different Correction (more optimal)
1959   // Sample:
1960   //  0.00000000000000E+00 -2.37822406690632E+01  6.03062748797469E+01
1961   //  7.70000000000000-100  7.70000000000000+100  7.70000000000000+100
1962   //0123456789012345678901234567890123456789012345678901234567890123456789
1963   const size_t posE = 18;
1964   std::string aStr (_curPos);
1965   if ( aStr.find('E') == std::string::npos && aStr.find('e') == std::string::npos )
1966     {
1967       if ( aStr.size() < posE+1 )
1968         THROW_IK_EXCEPTION("No more doubles (line #" << lineNb() << ")");
1969       aStr.insert( posE, "E", 1 );
1970       return atof(aStr.c_str());
1971     }
1972   return atof( _curPos );
1973 }
1974
1975 //================================================================================
1976 /*!
1977  * \brief Return the current string value
1978  */
1979 //================================================================================
1980
1981 std::string ASCIIReader::getName() const
1982 {
1983   int len = _width;
1984   while (( _curPos[len-1] == ' ' || _curPos[len-1] == 0) && len > 0 )
1985     len--;
1986   return std::string( _curPos, len );
1987 }
1988
1989 //================================================================================
1990 /*!
1991  * \brief Constructor of a binary sauve file reader
1992  */
1993 //================================================================================
1994
1995 XDRReader::XDRReader(const char* fileName) :FileReader(fileName), _xdrs_file(NULL)
1996 {
1997 }
1998
1999 //================================================================================
2000 /*!
2001  * \brief Close the XDR sauve file
2002  */
2003 //================================================================================
2004
2005 XDRReader::~XDRReader()
2006 {
2007 #ifdef HAS_XDR
2008   if ( _xdrs_file )
2009     {
2010       xdr_destroy((XDR*)_xdrs);
2011       free((XDR*)_xdrs);
2012       ::fclose(_xdrs_file);
2013       _xdrs_file = NULL;
2014     }
2015 #endif
2016 }
2017
2018 //================================================================================
2019 /*!
2020  * \brief Return false
2021  */
2022 //================================================================================
2023
2024 bool XDRReader::isASCII() const
2025 {
2026   return false;
2027 }
2028
2029 //================================================================================
2030 /*!
2031  * \brief Try to open an XRD file
2032  */
2033 //================================================================================
2034
2035 bool XDRReader::open()
2036 {
2037   bool xdr_ok = false;
2038 #ifdef HAS_XDR
2039 #ifdef WIN32
2040   if ((_xdrs_file = ::fopen(_fileName.c_str(), "rb")))
2041 #else
2042     if ((_xdrs_file = ::fopen(_fileName.c_str(), "r")))
2043 #endif
2044       {
2045         _xdrs = (XDR *)malloc(sizeof(XDR));
2046         xdrstdio_create((XDR*)_xdrs, _xdrs_file, XDR_DECODE);
2047
2048         const int maxsize = 10;
2049         char icha[maxsize+1];
2050         char* icha2 = icha;
2051         if (( xdr_ok = xdr_string((XDR*)_xdrs, &icha2, maxsize)))
2052           {
2053             icha[maxsize] = '\0';
2054             xdr_ok = (strcmp(icha, "CASTEM XDR") == 0);
2055           }
2056         if ( !xdr_ok )
2057           {
2058             xdr_destroy((XDR*)_xdrs);
2059             free((XDR*)_xdrs);
2060             fclose(_xdrs_file);
2061             _xdrs_file = NULL;
2062           }
2063       }
2064 #endif
2065   return xdr_ok;
2066 }
2067
2068 //================================================================================
2069 /*!
2070  * \brief A stub
2071  */
2072 //================================================================================
2073
2074 bool XDRReader::getNextLine (char* &, bool )
2075 {
2076   return true;
2077 }
2078
2079 //================================================================================
2080 /*!
2081  * \brief Prepare for iterating over given nb of values
2082  *  \param nbToRead - nb of fields to read
2083  *  \param width - field width
2084  */
2085 //================================================================================
2086
2087 void XDRReader::init( int nbToRead, int width/*=0*/ )
2088 {
2089   if(_iRead < _nbToRead)
2090     {
2091       std::cout << "_iRead, _nbToRead : " << _iRead << " " << _nbToRead << std::endl;
2092       std::cout << "Unfinished iteration before new one !" << std::endl;
2093       THROW_IK_EXCEPTION("SauvUtilities::XDRReader::init(): Unfinished iteration before new one !");
2094     }
2095   _iRead    = 0;
2096   _nbToRead = nbToRead;
2097   _width    = width;
2098 }
2099
2100 //================================================================================
2101 /*!
2102  * \brief Prepare for iterating over given nb of string values
2103  */
2104 //================================================================================
2105
2106 void XDRReader::initNameReading(int nbValues, int width)
2107 {
2108   init( nbValues, width );
2109   _xdr_kind = _xdr_kind_char;
2110   if(nbValues*width)
2111     {
2112       unsigned int nels = nbValues*width;
2113       _xdr_cvals = (char*)malloc((nels+1)*sizeof(char));
2114 #ifdef HAS_XDR
2115       xdr_string((XDR*)_xdrs, &_xdr_cvals, nels);
2116 #endif
2117       _xdr_cvals[nels] = '\0';
2118     }
2119 }
2120
2121 //================================================================================
2122 /*!
2123  * \brief Prepare for iterating over given nb of integer values
2124  */
2125 //================================================================================
2126
2127 void XDRReader::initIntReading(int nbValues)
2128 {
2129   init( nbValues );
2130   _xdr_kind = _xdr_kind_int;
2131   if(nbValues)
2132     {
2133 #ifdef HAS_XDR
2134       unsigned int nels = nbValues;
2135       unsigned int actual_nels;
2136       _xdr_ivals = (int*)malloc(nels*sizeof(int));
2137       xdr_array((XDR*)_xdrs, (char **)&_xdr_ivals, &actual_nels, nels, sizeof(int), (xdrproc_t)xdr_int);
2138 #endif
2139     }
2140 }
2141
2142 //================================================================================
2143 /*!
2144  * \brief Prepare for iterating over given nb of real values
2145  */
2146 //================================================================================
2147
2148 void XDRReader::initDoubleReading(int nbValues)
2149 {
2150   init( nbValues );
2151   _xdr_kind = _xdr_kind_double;
2152   if(nbValues)
2153     {
2154 #ifdef HAS_XDR
2155       unsigned int nels = nbValues;
2156       unsigned int actual_nels;
2157       _xdr_dvals = (double*)malloc(nels*sizeof(double));
2158       xdr_array((XDR*)_xdrs, (char **)&_xdr_dvals, &actual_nels, nels, sizeof(double), (xdrproc_t)xdr_double);
2159 #endif
2160     }
2161 }
2162
2163 //================================================================================
2164 /*!
2165  * \brief Return true if not all values have been read
2166  */
2167 //================================================================================
2168
2169 bool XDRReader::more() const
2170 {
2171   return _iRead < _nbToRead;
2172 }
2173
2174 //================================================================================
2175 /*!
2176  * \brief Go to the nex value
2177  */
2178 //================================================================================
2179
2180 void XDRReader::next()
2181 {
2182   if ( !more() )
2183     THROW_IK_EXCEPTION("SauvUtilities::XDRReader::next(): no more() values to read");
2184
2185   ++_iRead;
2186   if ( _iRead < _nbToRead )
2187     {
2188     }
2189   else
2190     {
2191       if(_xdr_kind == _xdr_kind_char) free(_xdr_cvals);
2192       if(_xdr_kind == _xdr_kind_int) free(_xdr_ivals);
2193       if(_xdr_kind == _xdr_kind_double) free(_xdr_dvals);
2194       _xdr_kind = _xdr_kind_null;
2195     }
2196 }
2197
2198 //================================================================================
2199 /*!
2200  * \brief Return the current integer value
2201  */
2202 //================================================================================
2203
2204 int XDRReader::getInt() const
2205 {
2206   if(_iRead < _nbToRead)
2207     {
2208       return _xdr_ivals[_iRead];
2209     }
2210   else
2211     {
2212       int result = 0;
2213 #ifdef HAS_XDR
2214       xdr_int((XDR*)_xdrs, &result);
2215 #endif
2216       return result;
2217     }
2218 }
2219
2220 //================================================================================
2221 /*!
2222  * \brief Return the current float value
2223  */
2224 //================================================================================
2225
2226 float  XDRReader::getFloat() const
2227 {
2228   float result = 0;
2229 #ifdef HAS_XDR
2230   xdr_float((XDR*)_xdrs, &result);
2231 #endif
2232   return result;
2233 }
2234
2235 //================================================================================
2236 /*!
2237  * \brief Return the current double value
2238  */
2239 //================================================================================
2240
2241 double XDRReader::getDouble() const
2242 {
2243   if(_iRead < _nbToRead)
2244     {
2245       return _xdr_dvals[_iRead];
2246     }
2247   else
2248     {
2249       double result = 0;
2250 #ifdef HAS_XDR
2251       xdr_double((XDR*)_xdrs, &result);
2252 #endif
2253       return result;
2254     }
2255 }
2256
2257 //================================================================================
2258 /*!
2259  * \brief Return the current string value
2260  */
2261 //================================================================================
2262
2263 std::string XDRReader::getName() const
2264 {
2265   int len = _width;
2266   char* s = _xdr_cvals + _iRead*_width;
2267   while (( s[len-1] == ' ' || s[len-1] == 0) && len > 0 )
2268     len--;
2269   return std::string( s, len );
2270 }
2271
2272 //================================================================================
2273 /*!
2274  * \brief Throw an exception if not all needed data is present
2275  */
2276 //================================================================================
2277
2278 void IntermediateMED::checkDataAvailability() const
2279 {
2280   if ( _spaceDim == 0 )
2281     THROW_IK_EXCEPTION("Wrong file format"); // it is the first record in the sauve file
2282
2283   if ( _groups.empty() )
2284     THROW_IK_EXCEPTION("No elements have been read");
2285
2286   if ( _points.empty() || _nbNodes == 0 )
2287     THROW_IK_EXCEPTION("Nodes of elements are not filled");
2288
2289   if ( _coords.empty() )
2290     THROW_IK_EXCEPTION("Node coordinates are missing");
2291
2292   if ( _coords.size() < _nbNodes * _spaceDim )
2293     THROW_IK_EXCEPTION("Nodes and coordinates mismatch");
2294 }
2295
2296 //================================================================================
2297 /*!
2298  * \brief Safely adds a new Group
2299  */
2300 //================================================================================
2301
2302 Group* IntermediateMED::addNewGroup(std::vector<SauvUtilities::Group*>* groupsToFix)
2303 {
2304   if ( _groups.size() == _groups.capacity() ) // re-allocation would occur
2305     {
2306       std::vector<Group> newGroups( _groups.size() );
2307       newGroups.push_back( Group() );
2308
2309       for ( size_t i = 0; i < _groups.size(); ++i )
2310         {
2311           // avoid copying _cells
2312           std::vector<const Cell*> cells;
2313           cells.swap( _groups[i]._cells );
2314           newGroups[i] = _groups[i];
2315           newGroups[i]._cells.swap( cells );
2316
2317           // correct pointers to sub-groups
2318           for ( size_t j = 0; j < _groups[i]._groups.size(); ++j )
2319             {
2320               std::size_t iG = _groups[i]._groups[j] - &_groups[0];
2321               newGroups[i]._groups[j] = & newGroups[ iG ];
2322             }
2323         }
2324
2325       // fix given groups
2326       if ( groupsToFix )
2327         for ( size_t i = 0; i < groupsToFix->size(); ++i )
2328           if ( (*groupsToFix)[i] )
2329             {
2330               std::size_t iG = (*groupsToFix)[i] - &_groups[0];
2331               (*groupsToFix)[i] = & newGroups[ iG ];
2332             }
2333
2334       // fix field supports
2335       for ( int isNode = 0; isNode < 2; ++isNode )
2336         {
2337           std::vector<DoubleField* >& fields = isNode ? _nodeFields : _cellFields;
2338           for ( size_t i = 0; i < fields.size(); ++i )
2339             {
2340               if ( !fields[i] ) continue;
2341               for ( size_t j = 0; j < fields[i]->_sub.size(); ++j )
2342                 if ( fields[i]->_sub[j]._support )
2343                   {
2344                     std::size_t iG = fields[i]->_sub[j]._support - &_groups[0];
2345                     fields[i]->_sub[j]._support = & newGroups[ iG ];
2346                   }
2347               if ( fields[i]->_group )
2348                 {
2349                   std::size_t iG = fields[i]->_group - &_groups[0];
2350                   fields[i]->_group = & newGroups[ iG ];
2351                 }
2352             }
2353         }
2354
2355       _groups.swap( newGroups );
2356     }
2357   else
2358     {
2359       _groups.push_back( Group() );
2360     }
2361   return &_groups.back();
2362 }
2363
2364 //================================================================================
2365 /*!
2366  * \brief Makes MEDCoupling::MEDFileData from self
2367  */
2368 //================================================================================
2369
2370 MEDCoupling::MEDFileData* IntermediateMED::convertInMEDFileDS()
2371 {
2372   MCAuto< MEDFileUMesh >  mesh   = makeMEDFileMesh();
2373   MCAuto< MEDFileFields > fields = makeMEDFileFields(mesh);
2374
2375   MCAuto< MEDFileMeshes > meshes = MEDFileMeshes::New();
2376   MCAuto< MEDFileData >  medData = MEDFileData::New();
2377   meshes->pushMesh( mesh );
2378   medData->setMeshes( meshes );
2379   if ( fields ) medData->setFields( fields );
2380
2381   return medData.retn();
2382 }
2383
2384 //================================================================================
2385 /*!
2386  * \brief Creates MEDCoupling::MEDFileUMesh from its data
2387  */
2388 //================================================================================
2389
2390 MEDCoupling::MEDFileUMesh* IntermediateMED::makeMEDFileMesh()
2391 {
2392   // check if all needed piles are present
2393   checkDataAvailability();
2394
2395   decreaseHierarchicalDepthOfSubgroups();
2396
2397   // set long names (before orienting!)
2398   setGroupLongNames();
2399
2400   // fix element orientation
2401   if ( _spaceDim == 2 || _spaceDim == 1 )
2402     orientElements2D();
2403   else if ( _spaceDim == 3 )
2404     orientElements3D();
2405
2406   // process groups
2407   eraseUselessGroups();
2408   //detectMixDimGroups();
2409
2410   // assign IDs
2411   _points.numberNodes();
2412   numberElements();
2413
2414   // make the med mesh
2415
2416   MEDFileUMesh* mesh = MEDFileUMesh::New();
2417
2418   DataArrayDouble *coords = getCoords();
2419   setConnectivity( mesh, coords );
2420   setGroups( mesh );
2421
2422   coords->decrRef();
2423
2424   if ( !mesh->getName().c_str() || strlen( mesh->getName().c_str() ) == 0 )
2425     mesh->setName( "MESH" );
2426
2427   return mesh;
2428 }
2429
2430 //================================================================================
2431 /*!
2432  * \brief Set long names to groups
2433  */
2434 //================================================================================
2435
2436 void IntermediateMED::setGroupLongNames()
2437 {
2438   if ( _listGIBItoMED_mail.empty() )
2439     return;
2440
2441   // IMP 0023285: only keep the meshes named in the table MED_MAIL
2442   // clear all group names
2443   for ( size_t i = 0; i < _groups.size(); ++i )
2444     if ( !_groups[i]._isProfile )
2445       _groups[i]._name.clear();
2446
2447
2448   // IMP 0020434: mapping GIBI names to MED names
2449   // set med names to objects (mesh, fields, support, group or other)
2450
2451   std::set<int> treatedGroups;
2452
2453   std::list<nameGIBItoMED>::iterator itGIBItoMED = _listGIBItoMED_mail.begin();
2454   for (; itGIBItoMED != _listGIBItoMED_mail.end(); itGIBItoMED++)
2455     {
2456       if ( (int)_groups.size() < itGIBItoMED->gibi_id ) continue;
2457
2458       SauvUtilities::Group & grp = _groups[itGIBItoMED->gibi_id - 1];
2459
2460       // if there are several names for grp then the 1st name is the name
2461       // of grp and the rest ones are names of groups referring grp (issue 0021311)
2462       const bool isRefName = !treatedGroups.insert( itGIBItoMED->gibi_id ).second;
2463       if ( !isRefName )
2464         {
2465           grp._name = _mapStrings[ itGIBItoMED->med_id ];
2466         }
2467       else if ( !grp._refNames.empty() && grp._refNames.back().empty() )
2468         {
2469           for ( unsigned i = 0; i < grp._refNames.size(); ++i )
2470             if ( grp._refNames[i].empty() )
2471               grp._refNames[i] = _mapStrings[ (*itGIBItoMED).med_id ];
2472         }
2473       else
2474         {
2475           grp._refNames.push_back( _mapStrings[ (*itGIBItoMED).med_id ]);
2476         }
2477     }
2478
2479   // IMP 0023285: only keep the meshes named in the table MED_MAIL
2480   // remove all cells belonging to non-named groups only
2481
2482   // use Cell::_reverse to mark cells to keep
2483   for ( size_t i = 0; i < _groups.size(); ++i )
2484     {
2485       SauvUtilities::Group & grp = _groups[i];
2486       if ( grp._isProfile || !grp._name.empty() )
2487         {
2488           for ( size_t iC = 0; iC < grp._cells.size(); ++iC )
2489             grp._cells[iC]->_reverse = true;
2490
2491           for (size_t j = 0; j < grp._groups.size(); ++j )
2492             for ( size_t iC = 0; iC < grp._groups[j]->_cells.size(); ++iC )
2493               grp._groups[j]->_cells[iC]->_reverse = true;
2494         }
2495     }
2496   // remove non-marked cells (with _reverse == false)
2497   CellsByDimIterator cellsIt( *this );
2498   while ( cellsIt.nextType() )
2499     {
2500       std::set<Cell> & cells = _cellsByType[ cellsIt.type() ];
2501       std::set<Cell>::iterator cIt = cells.begin();
2502       while ( cIt != cells.end() )
2503         if ( cIt->_reverse )
2504           {
2505             cIt->_reverse = false;
2506             ++cIt;
2507           }
2508         else
2509           {
2510             cells.erase( cIt++ );
2511           }
2512     }
2513 }
2514
2515 //================================================================================
2516 /*!
2517  * \brief Set long names to fields
2518  */
2519 //================================================================================
2520
2521 void IntermediateMED::setFieldLongNames(std::set< std::string >& usedNames)
2522 {
2523   std::list<nameGIBItoMED>::iterator itGIBItoMED = _listGIBItoMED_cham.begin();
2524   for (; itGIBItoMED != _listGIBItoMED_cham.end(); itGIBItoMED++)
2525     {
2526       if (itGIBItoMED->gibi_pile == PILE_FIELD)
2527         {
2528           _cellFields[itGIBItoMED->gibi_id - 1]->_name = _mapStrings[itGIBItoMED->med_id];
2529         }
2530       else if (itGIBItoMED->gibi_pile == PILE_NODES_FIELD)
2531         {
2532           _nodeFields[itGIBItoMED->gibi_id - 1]->_name = _mapStrings[itGIBItoMED->med_id];
2533         }
2534     } // iterate on _listGIBItoMED_cham
2535
2536   for (itGIBItoMED =_listGIBItoMED_comp.begin(); itGIBItoMED != _listGIBItoMED_comp.end(); itGIBItoMED++)
2537     {
2538       std::string medName  = _mapStrings[itGIBItoMED->med_id];
2539       std::string gibiName = _mapStrings[itGIBItoMED->gibi_id];
2540
2541       bool name_found = false;
2542       for ( int isNodal = 0; isNodal < 2 && !name_found; ++isNodal )
2543         {
2544           std::vector<DoubleField* > & fields = isNodal ? _nodeFields : _cellFields;
2545           for ( size_t ifi = 0; ifi < fields.size() && !name_found; ifi++)
2546             {
2547               if (medName.find( fields[ifi]->_name + "." ) == 0 )
2548                 {
2549                   std::vector<DoubleField::_Sub_data>& aSubDs = fields[ifi]->_sub;
2550                   std::size_t nbSub = aSubDs.size();
2551                   for (std::size_t isu = 0; isu < nbSub; isu++)
2552                     for (int ico = 0; ico < aSubDs[isu].nbComponents(); ico++)
2553                       {
2554                         if (aSubDs[isu].compName(ico) == gibiName)
2555                           {
2556                             std::string medNameCompo = medName.substr( fields[ifi]->_name.size() + 1 );
2557                             fields[ifi]->_sub[isu].compName(ico) = medNameCompo;
2558                           }
2559                       }
2560                 }
2561             }
2562         }
2563     } // iterate on _listGIBItoMED_comp
2564
2565   for ( size_t i = 0; i < _nodeFields.size() ; i++)
2566     usedNames.insert( _nodeFields[i]->_name );
2567   for ( size_t i = 0; i < _cellFields.size() ; i++)
2568     usedNames.insert( _cellFields[i]->_name );
2569 }
2570
2571 //================================================================================
2572 /*!
2573  * \brief Decrease hierarchical depth of subgroups
2574  */
2575 //================================================================================
2576
2577 void IntermediateMED::decreaseHierarchicalDepthOfSubgroups()
2578 {
2579   for (size_t i=0; i!=_groups.size(); ++i)
2580     {
2581       Group& grp = _groups[i];
2582       for (size_t j = 0; j < grp._groups.size(); ++j )
2583         {
2584           Group & sub_grp = *grp._groups[j];
2585           if ( !sub_grp._groups.empty() )
2586             {
2587               // replace j with its 1st subgroup
2588               grp._groups[j] = sub_grp._groups[0];
2589               // push back the rest subs
2590               grp._groups.insert( grp._groups.end(), ++sub_grp._groups.begin(), sub_grp._groups.end() );
2591             }
2592         }
2593       // remove empty sub-_groups
2594       std::vector< Group* > newSubGroups;
2595       newSubGroups.reserve( grp._groups.size() );
2596       for (size_t j = 0; j < grp._groups.size(); ++j )
2597         if ( !grp._groups[j]->empty() )
2598           newSubGroups.push_back( grp._groups[j] );
2599       if ( newSubGroups.size() < grp._groups.size() )
2600         grp._groups.swap( newSubGroups );
2601     }
2602 }
2603
2604 //================================================================================
2605 /*!
2606  * \brief Erase _groups that won't be converted
2607  */
2608 //================================================================================
2609
2610 void IntermediateMED::eraseUselessGroups()
2611 {
2612   // propagate _isProfile=true to sub-groups of composite groups
2613   // for (size_t int i=0; i!=_groups.size(); ++i)
2614   // {
2615   //   Group* grp = _groups[i];
2616   //   if ( grp->_isProfile && !grp->_groups.empty() )
2617   //     for (size_t j = 0; j < grp->_groups.size(); ++j )
2618   //       grp->_groups[j]->_isProfile=true;
2619   // }
2620   std::set<Group*> groups2convert;
2621   // keep not named sub-groups of field supports
2622   for (size_t i=0; i!=_groups.size(); ++i)
2623     {
2624       Group& grp = _groups[i];
2625       if ( grp._isProfile && !grp._groups.empty() )
2626         groups2convert.insert( grp._groups.begin(), grp._groups.end() );
2627     }
2628
2629   // keep named groups and their subgroups
2630   for (size_t i=0; i!=_groups.size(); ++i)
2631     {
2632       Group& grp = _groups[i];
2633       if ( !grp._name.empty() && !grp.empty() )
2634         {
2635           groups2convert.insert( &grp );
2636           groups2convert.insert( grp._groups.begin(), grp._groups.end() );
2637         }
2638     }
2639   // erase groups that are not in groups2convert and not _isProfile
2640   for (size_t i=0; i!=_groups.size(); ++i)
2641     {
2642       Group* grp = &_groups[i];
2643       if ( !grp->_isProfile && !groups2convert.count( grp ) )
2644         {
2645           grp->_cells.clear();
2646           grp->_groups.clear();
2647         }
2648     }
2649 }
2650
2651 //================================================================================
2652 /*!
2653  * \brief Detect _groups of mixed dimension
2654  */
2655 //================================================================================
2656
2657 void IntermediateMED::detectMixDimGroups()
2658 {
2659   //hasMixedCells = false;
2660   for ( size_t i=0; i < _groups.size(); ++i )
2661     {
2662       Group& grp = _groups[i];
2663       if ( grp._groups.size() < 2 )
2664         continue;
2665
2666       // check if sub-groups have different dimension
2667       unsigned dim1 = getDim( &grp );
2668       for ( size_t j = 1; j  < grp._groups.size(); ++j )
2669         {
2670           unsigned dim2 = getDim( grp._groups[j] );
2671           if ( dim1 != dim2 )
2672             {
2673               grp._cells.clear();
2674               grp._groups.clear();
2675               if ( !grp._name.empty() )
2676                 std::cout << "Erase a group with elements of different dim |" << grp._name << "|"<< std::endl;
2677               break;
2678             }
2679         }
2680     }
2681 }
2682
2683 //================================================================================
2684 /*!
2685  * \brief Fix connectivity of elements in 2D space
2686  */
2687 //================================================================================
2688
2689 void IntermediateMED::orientElements2D()
2690 {
2691   std::set<Cell>::const_iterator elemIt, elemEnd;
2692   std::vector< std::pair<int,int> > swapVec;
2693
2694   // ------------------------------------
2695   // fix connectivity of quadratic edges
2696   // ------------------------------------
2697   std::set<Cell>& quadEdges = _cellsByType[ INTERP_KERNEL::NORM_SEG3 ];
2698   if ( !quadEdges.empty() )
2699     {
2700       elemIt = quadEdges.begin(), elemEnd = quadEdges.end();
2701       for ( ; elemIt != elemEnd; ++elemIt )
2702         ConvertQuadratic( INTERP_KERNEL::NORM_SEG3, *elemIt );
2703     }
2704
2705   CellsByDimIterator faceIt( *this, 2 );
2706   while ( const std::set<Cell > * faces = faceIt.nextType() )
2707     {
2708       TCellType cellType = faceIt.type();
2709       bool isQuadratic = getGibi2MedQuadraticInterlace( cellType );
2710
2711       getReverseVector( cellType, swapVec );
2712
2713       // ------------------------------------
2714       // fix connectivity of quadratic faces
2715       // ------------------------------------
2716       if ( isQuadratic )
2717         for ( elemIt = faces->begin(), elemEnd = faces->end(); elemIt != elemEnd; elemIt++ )
2718           ConvertQuadratic( cellType, *elemIt );
2719
2720       // --------------------------
2721       // orient faces clockwise
2722       // --------------------------
2723       // COMMENTED for issue 0022612 note 17739
2724       // int iQuad = isQuadratic ? 2 : 1;
2725       // for ( elemIt = faces->begin(), elemEnd = faces->end(); elemIt != elemEnd; elemIt++ )
2726       //   {
2727       //     // look for index of the most left node
2728       //     int iLeft = 0, iNode, nbNodes = elemIt->_nodes.size() / iQuad;
2729       //     double x, minX = nodeCoords( elemIt->_nodes[0] )[0];
2730       //     for ( iNode = 1; iNode < nbNodes; ++iNode )
2731       //       if (( x = nodeCoords( elemIt->_nodes[ iNode ])[ 0 ]) < minX )
2732       //         minX = x, iLeft = iNode;
2733
2734       //     // indices of the nodes neighboring the most left one
2735       //     int iPrev = ( iLeft - 1 < 0 ) ? nbNodes - 1 : iLeft - 1;
2736       //     int iNext = ( iLeft + 1 == nbNodes ) ? 0 : iLeft + 1;
2737       //     // find components of prev-left and left-next vectors
2738       //     double xP = nodeCoords( elemIt->_nodes[ iPrev ])[ 0 ];
2739       //     double yP = nodeCoords( elemIt->_nodes[ iPrev ])[ 1 ];
2740       //     double xN = nodeCoords( elemIt->_nodes[ iNext ])[ 0 ];
2741       //     double yN = nodeCoords( elemIt->_nodes[ iNext ])[ 1 ];
2742       //     double xL = nodeCoords( elemIt->_nodes[ iLeft ])[ 0 ];
2743       //     double yL = nodeCoords( elemIt->_nodes[ iLeft ])[ 1 ];
2744       //     double xPL = xL - xP, yPL = yL - yP; // components of prev-left vector
2745       //     double xLN = xN - xL, yLN = yN - yL; // components of left-next vector
2746       //     // normalise y of the vectors
2747       //     double modPL = sqrt ( xPL * xPL + yPL * yPL );
2748       //     double modLN = sqrt ( xLN * xLN + yLN * yLN );
2749       //     if ( modLN > std::numeric_limits<double>::min() &&
2750       //          modPL > std::numeric_limits<double>::min() )
2751       //       {
2752       //         yPL /= modPL;
2753       //         yLN /= modLN;
2754       //         // summary direction of neighboring links must be positive
2755       //         bool clockwise = ( yPL + yLN > 0 );
2756       //         if ( !clockwise )
2757       //           reverse( *elemIt, swapVec );
2758       //       }
2759       //   }
2760     }
2761 }
2762
2763 //================================================================================
2764 /*!
2765  * \brief Fix connectivity of elements in 3D space
2766  */
2767 //================================================================================
2768
2769 void IntermediateMED::orientElements3D()
2770 {
2771   // set _reverse flags of faces
2772   // COMMENTED for issue 0022612 note 17739
2773   //orientFaces3D();
2774
2775   // -----------------
2776   // fix connectivity
2777   // -----------------
2778
2779   std::set<Cell>::const_iterator elemIt, elemEnd;
2780   std::vector< std::pair<int,int> > swapVec;
2781
2782   for ( int dim = 1; dim <= 3; ++dim )
2783     {
2784       CellsByDimIterator cellsIt( *this, dim );
2785       while ( const std::set<Cell > * elems = cellsIt.nextType() )
2786         {
2787           TCellType cellType = cellsIt.type();
2788           bool isQuadratic = getGibi2MedQuadraticInterlace( cellType );
2789           getReverseVector( cellType, swapVec );
2790
2791           elemIt = elems->begin(), elemEnd = elems->end();
2792           for ( ; elemIt != elemEnd; elemIt++ )
2793             {
2794               // GIBI connectivity -> MED one
2795               if( isQuadratic )
2796                 ConvertQuadratic( cellType, *elemIt );
2797
2798               // reverse faces
2799               if ( elemIt->_reverse )
2800                 reverse ( *elemIt, swapVec );
2801             }
2802         }
2803     }
2804
2805   // COMMENTED for issue 0022612 note 17739
2806   //orientVolumes();
2807 }
2808
2809 //================================================================================
2810 /*!
2811  * \brief Orient equally (by setting _reverse flag) all connected faces in 3D space
2812  */
2813 //================================================================================
2814
2815 void IntermediateMED::orientFaces3D()
2816 {
2817   // fill map of links and their faces
2818   std::set<const Cell*> faces;
2819   std::map<const Cell*, Group*> fgm;
2820   std::map<Link, std::list<const Cell*> > linkFacesMap;
2821   std::map<Link, std::list<const Cell*> >::iterator lfIt, lfIt2;
2822
2823   for (size_t i=0; i!=_groups.size(); ++i)
2824     {
2825       Group& grp = _groups[i];
2826       if ( !grp._cells.empty() && getDimension( grp._cellType ) == 2 )
2827         for ( size_t j = 0; j < grp._cells.size(); ++j )
2828           if ( faces.insert( grp._cells[j] ).second )
2829             {
2830               for ( unsigned int k = 0; k < grp._cells[j]->_nodes.size(); ++k )
2831                 linkFacesMap[ grp._cells[j]->link( k ) ].push_back( grp._cells[j] );
2832               fgm.insert( std::make_pair( grp._cells[j], &grp ));
2833             }
2834     }
2835   // dump linkFacesMap
2836   //     for ( lfIt = linkFacesMap.begin(); lfIt!=linkFacesMap.end(); lfIt++) {
2837   //       cout<< "LINK: " << lfIt->first.first << "-" << lfIt->first.second << std::endl;
2838   //       std::list<const Cell*> & fList = lfIt->second;
2839   //       std::list<const Cell*>::iterator fIt = fList.begin();
2840   //       for ( ; fIt != fList.end(); fIt++ )
2841   //         cout << "\t" << **fIt << fgm[*fIt]->nom << std::endl;
2842   //     }
2843
2844   // Each oriented link must appear in one face only, else a face is reversed.
2845
2846   std::queue<const Cell*> faceQueue; /* the queue contains well oriented faces
2847                                         whose neighbors orientation is to be checked */
2848   bool manifold = true;
2849   while ( !linkFacesMap.empty() )
2850     {
2851       if ( faceQueue.empty() )
2852         {
2853           assert( !linkFacesMap.begin()->second.empty() );
2854           faceQueue.push( linkFacesMap.begin()->second.front() );
2855         }
2856       while ( !faceQueue.empty() )
2857         {
2858           const Cell* face = faceQueue.front();
2859           faceQueue.pop();
2860
2861           // loop on links of <face>
2862           for ( int i = 0; i < (int)face->_nodes.size(); ++i )
2863             {
2864               Link link = face->link( i );
2865               // find the neighbor faces
2866               lfIt = linkFacesMap.find( link );
2867               int nbFaceByLink = 0;
2868               std::list< const Cell* > ml;
2869               if ( lfIt != linkFacesMap.end() )
2870                 {
2871                   std::list<const Cell*> & fList = lfIt->second;
2872                   std::list<const Cell*>::iterator fIt = fList.begin();
2873                   assert( fIt != fList.end() );
2874                   for ( ; fIt != fList.end(); fIt++, nbFaceByLink++ )
2875                     {
2876                       ml.push_back( *fIt );
2877                       if ( *fIt != face ) // wrongly oriented neighbor face
2878                         {
2879                           const Cell* badFace = *fIt;
2880                           // reverse and remove badFace from linkFacesMap
2881                           for ( int j = 0; j < (int)badFace->_nodes.size(); ++j )
2882                             {
2883                               Link badlink = badFace->link( j );
2884                               if ( badlink == link ) continue;
2885                               lfIt2 = linkFacesMap.find( badlink );
2886                               if ( lfIt2 != linkFacesMap.end() )
2887                                 {
2888                                   std::list<const Cell*> & ff = lfIt2->second;
2889                                   std::list<const Cell*>::iterator lfIt3 = find( ff.begin(), ff.end(), badFace );
2890                                   // check if badFace has been found,
2891                                   // else we can't erase it
2892                                   // case of degenerated face in edge
2893                                   if (lfIt3 != ff.end())
2894                                     {
2895                                       ff.erase( lfIt3 );
2896                                       if ( ff.empty() )
2897                                         linkFacesMap.erase( lfIt2 );
2898                                     }
2899                                 }
2900                             }
2901                           badFace->_reverse = true; // reverse
2902                           //INFOS_MED( "REVERSE " << *badFace );
2903                           faceQueue.push( badFace );
2904                         }
2905                     }
2906                   linkFacesMap.erase( lfIt );
2907                 }
2908               // add good neighbors to the queue
2909               Link revLink( link.second, link.first );
2910               lfIt = linkFacesMap.find( revLink );
2911               if ( lfIt != linkFacesMap.end() )
2912                 {
2913                   std::list<const Cell*> & fList = lfIt->second;
2914                   std::list<const Cell*>::iterator fIt = fList.begin();
2915                   for ( ; fIt != fList.end(); fIt++, nbFaceByLink++ )
2916                     {
2917                       ml.push_back( *fIt );
2918                       if ( *fIt != face )
2919                         faceQueue.push( *fIt );
2920                     }
2921                   linkFacesMap.erase( lfIt );
2922                 }
2923               if ( nbFaceByLink > 2 )
2924                 {
2925                   if ( manifold )
2926                     {
2927                       std::list<const Cell*>::iterator ii = ml.begin();
2928                       std::cout << nbFaceByLink << " faces by 1 link:" << std::endl;
2929                       for( ; ii!= ml.end(); ii++ )
2930                         std::cout << "in sub-mesh <" << fgm[ *ii ]->_name << "> " << **ii << std::endl;
2931                     }
2932                   manifold = false;
2933                 }
2934             } // loop on links of the being checked face
2935         } // loop on the face queue
2936     } // while ( !linkFacesMap.empty() )
2937
2938   if ( !manifold )
2939     std::cout << " -> Non manifold mesh, faces orientation may be incorrect" << std::endl;
2940 }
2941
2942 //================================================================================
2943 /*!
2944  * \brief Orient volumes according to MED conventions:
2945  * normal of a bottom (first) face should be outside
2946  */
2947 //================================================================================
2948
2949 void IntermediateMED::orientVolumes()
2950 {
2951   std::set<Cell>::const_iterator elemIt, elemEnd;
2952   std::vector< std::pair<int,int> > swapVec;
2953
2954   CellsByDimIterator cellsIt( *this, 3 );
2955   while ( const std::set<Cell > * elems = cellsIt.nextType() )
2956     {
2957       TCellType cellType = cellsIt.type();
2958       elemIt = elems->begin(), elemEnd = elems->end();
2959       int nbBottomNodes = 0;
2960       switch ( cellType )
2961         {
2962         case NORM_TETRA4:
2963         case NORM_TETRA10:
2964         case NORM_PENTA6:
2965         case NORM_PENTA15:
2966           nbBottomNodes = 3; break;
2967         case NORM_PYRA5:
2968         case NORM_PYRA13:
2969         case NORM_HEXA8:
2970         case NORM_HEXA20:
2971           nbBottomNodes = 4; break;
2972         default: continue;
2973         }
2974       getReverseVector( cellType, swapVec );
2975
2976       for ( ; elemIt != elemEnd; elemIt++ )
2977         {
2978           // find a normal to the bottom face
2979           const double* n[4];
2980           n[0] = nodeCoords( elemIt->_nodes[0]); // 3 bottom nodes
2981           n[1] = nodeCoords( elemIt->_nodes[1]);
2982           n[2] = nodeCoords( elemIt->_nodes[2]);
2983           n[3] = nodeCoords( elemIt->_nodes[nbBottomNodes]); // a top node
2984           double vec01[3]; // vector n[0]-n[1]
2985           vec01[0] = n[1][0] - n[0][0];
2986           vec01[1] = n[1][1] - n[0][1];
2987           vec01[2] = n[1][2] - n[0][2];
2988           double vec02 [3]; // vector n[0]-n[2]
2989           vec02[0] = n[2][0] - n[0][0];
2990           vec02[1] = n[2][1] - n[0][1];
2991           vec02[2] = n[2][2] - n[0][2];
2992           double normal [3]; // vec01 ^ vec02
2993           normal[0] = vec01[1] * vec02[2] - vec01[2] * vec02[1];
2994           normal[1] = vec01[2] * vec02[0] - vec01[0] * vec02[2];
2995           normal[2] = vec01[0] * vec02[1] - vec01[1] * vec02[0];
2996           // check if the 102 angle is convex
2997           if ( nbBottomNodes > 3 )
2998             {
2999               const double* n3 = nodeCoords( elemIt->_nodes[nbBottomNodes-1] );// last bottom node
3000               double vec03 [3];  // vector n[0]-n3
3001               vec03[0] = n3[0] - n[0][0];
3002               vec03[1] = n3[1] - n[0][1];
3003               vec03[2] = n3[2] - n[0][2];
3004               if ( fabs( normal[0]+normal[1]+normal[2] ) <= std::numeric_limits<double>::max() ) // vec01 || vec02
3005                 {
3006                   normal[0] = vec01[1] * vec03[2] - vec01[2] * vec03[1]; // vec01 ^ vec03
3007                   normal[1] = vec01[2] * vec03[0] - vec01[0] * vec03[2];
3008                   normal[2] = vec01[0] * vec03[1] - vec01[1] * vec03[0];
3009                 }
3010               else
3011                 {
3012                   double vec [3]; // normal ^ vec01
3013                   vec[0] = normal[1] * vec01[2] - normal[2] * vec01[1];
3014                   vec[1] = normal[2] * vec01[0] - normal[0] * vec01[2];
3015                   vec[2] = normal[0] * vec01[1] - normal[1] * vec01[0];
3016                   double dot2 = vec[0]*vec03[0] + vec[1]*vec03[1] + vec[2]*vec03[2]; // vec*vec03
3017                   if ( dot2 < 0 ) // concave -> reverse normal
3018                     {
3019                       normal[0] *= -1;
3020                       normal[1] *= -1;
3021                       normal[2] *= -1;
3022                     }
3023                 }
3024             }
3025           // direction from top to bottom
3026           double tbDir[3];
3027           tbDir[0] = n[0][0] - n[3][0];
3028           tbDir[1] = n[0][1] - n[3][1];
3029           tbDir[2] = n[0][2] - n[3][2];
3030
3031           // compare 2 directions: normal and top-bottom
3032           double dot = normal[0]*tbDir[0] + normal[1]*tbDir[1] + normal[2]*tbDir[2];
3033           if ( dot < 0. ) // need reverse
3034             reverse( *elemIt, swapVec );
3035
3036         } // loop on volumes of one geometry
3037     } // loop on 3D geometry types
3038
3039 }
3040
3041 //================================================================================
3042 /*!
3043  * \brief Assign new IDs to nodes by skipping not used nodes and return their number
3044  */
3045 //================================================================================
3046
3047 int NodeContainer::numberNodes()
3048 {
3049   int id = 1;
3050   for ( size_t i = 0; i < _nodes.size(); ++i )
3051     for ( size_t j = 0; j < _nodes[i].size(); ++j )
3052       if ( _nodes[i][j].isUsed() )
3053         _nodes[i][j]._number = id++;
3054   return id-1;
3055 }
3056
3057
3058 //================================================================================
3059 /*!
3060  * \brief Assign new IDs to elements
3061  */
3062 //================================================================================
3063
3064 void IntermediateMED::numberElements()
3065 {
3066   std::set<Cell>::const_iterator elemIt, elemEnd;
3067
3068   // numbering _cells of type NORM_POINT1 by node number
3069   {
3070     const std::set<Cell>& points = _cellsByType[ INTERP_KERNEL::NORM_POINT1 ];
3071     elemIt = points.begin(), elemEnd = points.end();
3072     for ( ; elemIt != elemEnd; ++elemIt )
3073       elemIt->_number = elemIt->_nodes[0]->_number;
3074   }
3075
3076   // numbering 1D-3D _cells
3077   for ( int dim = 1; dim <= 3; ++dim )
3078     {
3079       // check if re-numeration is needed (to try to keep elem oreder as in sauve file )
3080       bool ok = true, renumEntity = false;
3081       CellsByDimIterator cellsIt( *this, dim );
3082       mcIdType prevNbElems = 0;
3083       while ( const std::set<Cell> * typeCells = cellsIt.nextType() )
3084         {
3085           TID minNumber = std::numeric_limits<TID>::max(), maxNumber = 0;
3086           for ( elemIt = typeCells->begin(), elemEnd = typeCells->end(); elemIt!=elemEnd; ++elemIt)
3087             {
3088               if ( elemIt->_number < minNumber ) minNumber = elemIt->_number;
3089               if ( elemIt->_number > maxNumber ) maxNumber = elemIt->_number;
3090             }
3091           mcIdType typeSize = ToIdType( typeCells->size() );
3092           if ( typeSize != maxNumber - minNumber + 1 )
3093             ok = false;
3094           if ( prevNbElems+1 != (int)minNumber )
3095             ok = false;
3096           if ( prevNbElems != 0 && minNumber == 1 )
3097             renumEntity = true;
3098
3099           prevNbElems += typeSize;
3100         }
3101
3102       if ( ok && renumEntity ) // each geom type was numerated separately
3103         {
3104           cellsIt.init( dim );
3105           prevNbElems = ToIdType( cellsIt.nextType()->size()); // no need to renumber the first type
3106           while ( const std::set<Cell> * typeCells = cellsIt.nextType() )
3107             {
3108               for ( elemIt = typeCells->begin(), elemEnd = typeCells->end(); elemIt!=elemEnd; ++elemIt)
3109                 elemIt->_number += prevNbElems;
3110               prevNbElems += ToIdType( typeCells->size() );
3111             }
3112         }
3113       if ( !ok )
3114         {
3115           int cellID=1;
3116           cellsIt.init( dim );
3117           while ( const std::set<Cell> * typeCells = cellsIt.nextType() )
3118             for ( elemIt = typeCells->begin(), elemEnd = typeCells->end(); elemIt!=elemEnd; ++elemIt)
3119               elemIt->_number = cellID++;
3120         }
3121     }
3122 }
3123
3124 //================================================================================
3125 /*!
3126  * \brief Creates coord array
3127  */
3128 //================================================================================
3129
3130 MEDCoupling::DataArrayDouble * IntermediateMED::getCoords()
3131 {
3132   DataArrayDouble* coordArray = DataArrayDouble::New();
3133   coordArray->alloc( _nbNodes, _spaceDim );
3134   double * coordPrt = coordArray->getPointer();
3135   for ( unsigned int i = 0; i < _points.size(); ++i )
3136     {
3137       Node* n = getNode( i+1 );
3138       if ( n->isUsed() )
3139         {
3140           const double* nCoords = nodeCoords( n );
3141           std::copy( nCoords, nCoords+_spaceDim, coordPrt );
3142           coordPrt += _spaceDim;
3143         }
3144     }
3145   return coordArray;
3146 }
3147
3148 //================================================================================
3149 /*!
3150  * \brief Sets connectivity of elements to the mesh
3151  *  \param mesh - mesh to fill in
3152  *  \param coords - coordinates that must be shared by all meshes of different dim
3153  */
3154 //================================================================================
3155
3156 void IntermediateMED::setConnectivity( MEDCoupling::MEDFileUMesh*    mesh,
3157                                        MEDCoupling::DataArrayDouble* coords )
3158 {
3159   int meshDim = 0;
3160
3161   mesh->setCoords( coords );
3162
3163   std::set<Cell>::const_iterator elemIt, elemEnd;
3164   for ( int dim = 3; dim > 0; --dim )
3165     {
3166       CellsByDimIterator dimCells( *this, dim );
3167
3168       mcIdType nbOfCells = 0;
3169       while ( const std::set<Cell > * cells = dimCells.nextType() )
3170         nbOfCells += ToIdType( cells->size() );
3171       if ( nbOfCells == 0 )
3172         continue;
3173
3174       if ( !meshDim ) meshDim = dim;
3175
3176       MEDCouplingUMesh* dimMesh = MEDCouplingUMesh::New();
3177       dimMesh->setCoords( coords );
3178       dimMesh->setMeshDimension( dim );
3179       dimMesh->allocateCells( nbOfCells );
3180
3181       mcIdType prevNbCells = 0;
3182       dimCells.init( dim );
3183       while ( const std::set<Cell > * cells = dimCells.nextType() )
3184         {
3185           // fill connectivity array to take into account order of elements in the sauv file
3186           const mcIdType nbCellNodes = ToIdType( cells->begin()->_nodes.size() );
3187           std::vector< TID > connectivity( cells->size() * nbCellNodes );
3188           TID * nodalConnOfCell;
3189           for ( elemIt = cells->begin(), elemEnd = cells->end(); elemIt != elemEnd; ++elemIt )
3190             {
3191               const Cell& cell = *elemIt;
3192               const TID index = cell._number - 1 - prevNbCells;
3193               nodalConnOfCell = &connectivity[ index * nbCellNodes ];
3194               if ( cell._reverse )
3195                 for ( mcIdType i = nbCellNodes-1; i >= 0; --i )
3196                   *nodalConnOfCell++ = cell._nodes[i]->_number - 1;
3197               else
3198                 for ( mcIdType i = 0; i < nbCellNodes; ++i )
3199                   *nodalConnOfCell++ = cell._nodes[i]->_number - 1;
3200             }
3201           prevNbCells += ToIdType( cells->size() );
3202
3203           // fill dimMesh
3204           TCellType cellType = dimCells.type();
3205           nodalConnOfCell = &connectivity[0];
3206           for ( size_t i = 0; i < cells->size(); ++i, nodalConnOfCell += nbCellNodes )
3207             dimMesh->insertNextCell( cellType, nbCellNodes, nodalConnOfCell );
3208         }
3209       dimMesh->finishInsertingCells();
3210       mesh->setMeshAtLevel( dim - meshDim, dimMesh );
3211       dimMesh->decrRef();
3212     }
3213 }
3214
3215 //================================================================================
3216 /*!
3217  * \brief Fill in the mesh with groups
3218  *  \param mesh - mesh to fill in
3219  */
3220 //================================================================================
3221
3222 void IntermediateMED::setGroups( MEDCoupling::MEDFileUMesh* mesh )
3223 {
3224   bool isMeshNameSet = false;
3225   const int meshDim = mesh->getMeshDimension();
3226   for ( int dim = 0; dim <= meshDim; ++dim )
3227     {
3228       const int meshDimRelToMaxExt = ( dim == 0 ? 1 : dim - meshDim );
3229
3230       std::vector<const DataArrayIdType *> medGroups;
3231       std::vector<MCAuto<DataArrayIdType> > refGroups;
3232       for ( size_t i = 0; i < _groups.size(); ++i )
3233         {
3234           Group& grp = _groups[i];
3235           if ( (int)getDim( &grp ) != dim &&
3236                grp._groups.empty() ) // to allow groups on diff dims
3237             continue;
3238           // convert only named groups or field supports
3239           if ( grp.empty() || (grp._name.empty() && !grp._isProfile ))
3240             continue;
3241           //if ( grp._medGroup ) continue; // already converted
3242
3243           // sort cells by ID and remember their initial order in the group
3244           TCellToOrderMap cell2order;
3245           unsigned orderInGroup = 0;
3246           std::vector< Group* > groupVec;
3247           if ( grp._groups.empty() ) groupVec.push_back( & grp );
3248           else                       groupVec = grp._groups;
3249           for ( size_t iG = 0; iG < groupVec.size(); ++iG )
3250             {
3251               Group* aG = groupVec[ iG ];
3252               if ( (int)getDim( aG ) != dim )
3253                 continue;
3254               for ( size_t iC = 0; iC < aG->_cells.size(); ++iC )
3255                 cell2order.insert( cell2order.end(), std::make_pair( aG->_cells[iC], orderInGroup++ ));
3256             }
3257           if ( cell2order.empty() )
3258             continue;
3259           bool isSelfIntersect = ( orderInGroup != cell2order.size() );
3260           if ( isSelfIntersect ) // self intersecting group
3261             {
3262               std::ostringstream msg;
3263               msg << "Self intersecting sub-mesh: id = " << i+1
3264                   << ", name = |" << grp._name << "|" << std::endl
3265                   << " nb unique elements = " << cell2order.size() << std::endl
3266                   << " total nb elements  = " << orderInGroup;
3267               if ( grp._isProfile )
3268                 {
3269                   THROW_IK_EXCEPTION( msg.str() );
3270                 }
3271               else
3272                 {
3273                   std::cout << msg.str() << std::endl;
3274                 }
3275             }
3276           // create a med group
3277           grp._medGroup = DataArrayIdType::New();
3278           grp._medGroup->setName( grp._name.c_str() );
3279           grp._medGroup->alloc( cell2order.size(), /*nbOfCompo=*/1 );
3280           TID * idsPtr = grp._medGroup->getPointer();
3281           TCellToOrderMap::iterator cell2orderIt, cell2orderEnd = cell2order.end();
3282           for ( cell2orderIt = cell2order.begin(); cell2orderIt != cell2orderEnd; ++cell2orderIt )
3283             *idsPtr++ = (*cell2orderIt).first->_number - 1;
3284
3285           // try to set the mesh name
3286           if ( !isMeshNameSet &&
3287                dim == meshDim &&
3288                !grp._name.empty() &&
3289                grp.size() == mesh->getSizeAtLevel( meshDimRelToMaxExt ))
3290             {
3291               mesh->setName( grp._name.c_str() );
3292               isMeshNameSet = true;
3293             }
3294           if ( !grp._name.empty() )
3295             {
3296               medGroups.push_back( grp._medGroup );
3297             }
3298           // set relocation table
3299           setRelocationTable( &grp, cell2order );
3300
3301           // Issue 0021311. Use case: a gibi group has references (recorded in pile 1)
3302           // and several names (pile 27) refer (pile 10) to this group.
3303           // We create a copy of this group per each named reference
3304           std::set<std::string> uniqueNames;
3305           uniqueNames.insert( grp._name );
3306           for ( unsigned iRef = 0 ; iRef < grp._refNames.size(); ++iRef )
3307             if ( !grp._refNames[ iRef ].empty() &&
3308                  uniqueNames.insert( grp._refNames[ iRef ]).second ) // for name uniqueness (23155)
3309               {
3310                 refGroups.push_back( grp._medGroup->deepCopy() );
3311                 refGroups.back()->setName( grp._refNames[ iRef ].c_str() );
3312                 medGroups.push_back( refGroups.back() );
3313               }
3314         }
3315       mesh->setGroupsAtLevel( meshDimRelToMaxExt, medGroups );
3316     }
3317 }
3318
3319 //================================================================================
3320 /*!
3321  * \brief Return true if the group is on all elements and return its relative dimension
3322  */
3323 //================================================================================
3324
3325 bool IntermediateMED::isOnAll( const Group* grp, int & dimRel ) const
3326 {
3327   int dim = getDim( grp );
3328
3329   mcIdType nbElems = 0;
3330   if ( dim == 0 )
3331     {
3332       nbElems = _nbNodes;
3333       dimRel  = 0;
3334     }
3335   else
3336     {
3337       CellsByDimIterator dimCells( *this, dim );
3338       while ( const std::set<Cell > * cells = dimCells.nextType() )
3339         nbElems += ToIdType( cells->size() );
3340
3341       int meshDim = 3;
3342       for ( ; meshDim > 0; --meshDim )
3343         {
3344           dimCells.init( meshDim );
3345           if ( dimCells.nextType() )
3346             break;
3347         }
3348       dimRel = dim - meshDim;
3349     }
3350
3351   bool onAll = ( nbElems == grp->size() );
3352   return onAll;
3353 }
3354
3355 //================================================================================
3356 /*!
3357  * \brief Makes fields from own data
3358  */
3359 //================================================================================
3360
3361 MEDCoupling::MEDFileFields * IntermediateMED::makeMEDFileFields(MEDCoupling::MEDFileUMesh* mesh)
3362 {
3363   if ( _nodeFields.empty() && _cellFields.empty() ) return 0;
3364
3365   // set long names
3366   std::set< std::string > usedFieldNames;
3367   setFieldLongNames(usedFieldNames);
3368
3369   MEDFileFields* fields = MEDFileFields::New();
3370
3371   for ( unsigned int i = 0; i < _nodeFields.size(); ++i )
3372     setFields( _nodeFields[i], fields, mesh, i+1, usedFieldNames );
3373
3374   for ( unsigned int i = 0; i < _cellFields.size(); ++i )
3375     setFields( _cellFields[i], fields, mesh, i+1, usedFieldNames );
3376
3377   return fields;
3378 }
3379
3380 //================================================================================
3381 /*!
3382  * \brief Make med fields from a SauvUtilities::DoubleField
3383  */
3384 //================================================================================
3385
3386 void IntermediateMED::setFields( SauvUtilities::DoubleField* fld,
3387                                  MEDCoupling::MEDFileFields*  medFields,
3388                                  MEDCoupling::MEDFileUMesh*   mesh,
3389                                  const TID                   castemID,
3390                                  std::set< std::string >&    usedFieldNames)
3391 {
3392   bool sameNbGauss = true;
3393   if ( !fld || !fld->isMedCompatible( sameNbGauss )) return;
3394
3395   if ( !sameNbGauss )
3396     fld->splitSubWithDiffNbGauss();
3397
3398   // if ( !fld->hasCommonSupport() ):
3399   //     each sub makes MEDFileFieldMultiTS
3400   // else:
3401   //     unite several subs into a MEDCouplingFieldDouble
3402
3403   const bool uniteSubs = fld->hasCommonSupport() && sameNbGauss;
3404   if ( !uniteSubs )
3405     std::cout << "Castem field #" << castemID << " <" << fld->_name
3406               << "> is incompatible with MED format, so we split it into several fields:" << std::endl;
3407
3408   for ( unsigned int iSub = 0; iSub < fld->_sub.size(); )
3409     {
3410       // set field name
3411       if ( !uniteSubs || fld->_name.empty() )
3412         makeFieldNewName( usedFieldNames, fld );
3413
3414       // allocate values
3415       DataArrayDouble * values = DataArrayDouble::New();
3416       values->alloc( fld->getNbTuples(iSub), fld->_sub[iSub].nbComponents() );
3417
3418       // set values
3419       double * valPtr = values->getPointer();
3420       if ( uniteSubs )
3421         {
3422           mcIdType nbElems = fld->_group->size();
3423           for ( mcIdType elemShift = 0; elemShift < nbElems && iSub < fld->_sub.size(); )
3424             elemShift += fld->setValues( valPtr, iSub++, elemShift );
3425           setTS( fld, values, medFields, mesh );
3426         }
3427       else
3428         {
3429           fld->setValues( valPtr, iSub );
3430           setTS( fld, values, medFields, mesh, iSub++ );
3431
3432           std::cout << fld->_name << " with components";
3433           for ( size_t i = 0; i < (size_t)fld->_sub[iSub-1].nbComponents(); ++i )
3434             std::cout << " " << fld->_sub[iSub-1]._comp_names[ i ];
3435           std::cout << std::endl;
3436         }
3437     }
3438 }
3439
3440 //================================================================================
3441 /*!
3442  * \brief Store value array of a field into med fields
3443  */
3444 //================================================================================
3445
3446 void IntermediateMED::setTS( SauvUtilities::DoubleField*  fld,
3447                              MEDCoupling::DataArrayDouble* values,
3448                              MEDCoupling::MEDFileFields*   medFields,
3449                              MEDCoupling::MEDFileUMesh*    mesh,
3450                              const int                    iSub)
3451 {
3452   // treat a field support
3453   const Group* support = fld->getSupport( iSub );
3454   int dimRel;
3455   const bool onAll = isOnAll( support, dimRel );
3456   if ( !onAll && support->_name.empty() )
3457     {
3458       const_cast<Group*>(support)->_name += "PFL_" + fld->_name;
3459       support->_medGroup->setName( support->_name.c_str() );
3460     }
3461
3462   // make and fill a time-stamp
3463
3464   MEDCouplingFieldDouble * timeStamp = MEDCouplingFieldDouble::New( fld->getMedType( iSub ),
3465                                                                     fld->getMedTimeDisc() );
3466   timeStamp->setName( fld->_name.c_str() );
3467   timeStamp->setDescription( fld->_description.c_str() );
3468   // set the mesh
3469   if ( onAll )
3470     {
3471       MCAuto
3472         < MEDCouplingUMesh > dimMesh = mesh->getMeshAtLevel( dimRel );
3473       timeStamp->setMesh( dimMesh );
3474     }
3475   else if ( timeStamp->getTypeOfField() == MEDCoupling::ON_NODES )
3476     {
3477       DataArrayDouble * coo = mesh->getCoords();
3478       MCAuto
3479         <DataArrayDouble> subCoo = coo->selectByTupleId(support->_medGroup->begin(),
3480                                                         support->_medGroup->end());
3481       MCAuto< MEDCouplingUMesh > nodeSubMesh =
3482         MEDCouplingUMesh::Build0DMeshFromCoords( subCoo );
3483       timeStamp->setMesh( nodeSubMesh );
3484     }
3485   else
3486     {
3487       MCAuto
3488         < MEDCouplingUMesh > dimMesh = mesh->getMeshAtLevel( dimRel );
3489       MCAuto
3490         <MEDCouplingMesh> subMesh = dimMesh->buildPart(support->_medGroup->begin(),
3491                                                        support->_medGroup->end());
3492       timeStamp->setMesh( subMesh);
3493     }
3494   // set values
3495   for ( size_t i = 0; i < (size_t)fld->_sub[iSub].nbComponents(); ++i )
3496     values->setInfoOnComponent( i, fld->_sub[iSub]._comp_names[ i ].c_str() );
3497   timeStamp->setArray( values );
3498   values->decrRef();
3499   // set gauss points
3500   if ( timeStamp->getTypeOfField() == MEDCoupling::ON_GAUSS_PT )
3501     {
3502       TGaussDef gaussDef( fld->_sub[iSub]._support->_cellType,
3503                           fld->_sub[iSub].nbGauss() );
3504       timeStamp->setGaussLocalizationOnType( fld->_sub[iSub]._support->_cellType,
3505                                              gaussDef.myRefCoords,
3506                                              gaussDef.myCoords,
3507                                              gaussDef.myWeights );
3508     }
3509   // get a field to add the time-stamp
3510   bool isNewMedField = false;
3511   if ( !fld->_curMedField || fld->_name != fld->_curMedField->getName() )
3512     {
3513       fld->_curMedField = MEDFileFieldMultiTS::New();
3514       isNewMedField = true;
3515     }
3516
3517   // set an order
3518   const int nbTS = fld->_curMedField->getNumberOfTS();
3519   if ( nbTS > 0 )
3520     timeStamp->setOrder( nbTS );
3521
3522   // add the time-stamp
3523   timeStamp->checkConsistencyLight();
3524   if ( onAll )
3525     fld->_curMedField->appendFieldNoProfileSBT( timeStamp );
3526   else
3527     fld->_curMedField->appendFieldProfile( timeStamp, mesh, dimRel, support->_medGroup );
3528   timeStamp->decrRef();
3529
3530   if ( isNewMedField ) // timeStamp must be added before this
3531     {
3532       medFields->pushField( fld->_curMedField );
3533     }
3534 }
3535
3536 //================================================================================
3537 /*!
3538  * \brief Make a new unique name for a field
3539  */
3540 //================================================================================
3541
3542 void IntermediateMED::makeFieldNewName(std::set< std::string >&    usedNames,
3543                                        SauvUtilities::DoubleField* fld )
3544 {
3545   std::string base = fld->_name;
3546   if ( base.empty() )
3547     {
3548       base = "F_";
3549     }
3550   else
3551     {
3552       std::string::size_type pos = base.rfind('_');
3553       if ( pos != std::string::npos )
3554         base = base.substr( 0, pos+1 );
3555       else
3556         base += '_';
3557     }
3558
3559   int i = 1;
3560   do
3561     {
3562       fld->_name = base + SauvUtilities::toString( i++ );
3563     }
3564   while( !usedNames.insert( fld->_name ).second );
3565 }
3566
3567 //================================================================================
3568 /*!
3569  * \brief Split sub-components with different nb of gauss points into several sub-components
3570  */
3571 //================================================================================
3572
3573 void DoubleField::splitSubWithDiffNbGauss()
3574 {
3575   for ( size_t iSub = 0; iSub < _sub.size(); ++iSub )
3576     {
3577       if ( _sub[iSub].isSameNbGauss() ) continue;
3578
3579       _sub.insert( _sub.begin() + iSub + 1, 1, _Sub_data() );
3580       _Sub_data & subToSplit = _sub[iSub];
3581       _Sub_data & subNew     = _sub[iSub+1];
3582       size_t iDiff = 1;
3583       while ( subToSplit._nb_gauss[ 0 ] == subToSplit._nb_gauss[ iDiff ] )
3584         ++iDiff;
3585       subNew._support = subToSplit._support;
3586       subNew._comp_names.assign( subToSplit._comp_names.begin() + iDiff,
3587                                  subToSplit._comp_names.end() );
3588       subNew._nb_gauss.assign  ( subToSplit._nb_gauss.begin() + iDiff,
3589                                  subToSplit._nb_gauss.end() );
3590       subToSplit._comp_names.resize( iDiff );
3591       subToSplit._nb_gauss.resize  ( iDiff );
3592     }
3593 }
3594
3595 //================================================================================
3596 /*!
3597  * \brief Return a vector ready to fill in
3598  */
3599 //================================================================================
3600
3601 std::vector< double >& DoubleField::addComponent( int nb_values )
3602 {
3603   _comp_values.push_back( std::vector< double >() );
3604   std::vector< double >& res = _comp_values.back();
3605   res.resize( nb_values );
3606   return res;
3607 }
3608
3609 DoubleField::~DoubleField()
3610 {
3611   if(_curMedField)
3612     _curMedField->decrRef();
3613 }
3614
3615 //================================================================================
3616 /*!
3617  * \brief Returns a supporting group
3618  */
3619 //================================================================================
3620
3621 const Group* DoubleField::getSupport( const int iSub ) const
3622 {
3623   return _group ? _group : _sub[iSub]._support;
3624 }
3625
3626 //================================================================================
3627 /*!
3628  * \brief Return true if each sub-component is a time stamp
3629  */
3630 //================================================================================
3631
3632 bool DoubleField::isMultiTimeStamps() const
3633 {
3634   if ( _sub.size() < 2 )
3635     return false;
3636   bool sameSupports = true;
3637   Group* grpp1 = _sub[0]._support;// grpp NOT grp because XDR under Windows defines grp...
3638   for ( size_t i = 1; i < _sub.size() && sameSupports; ++i )
3639     sameSupports = ( grpp1 == _sub[i]._support );
3640
3641   return sameSupports;
3642 }
3643
3644 //================================================================================
3645 /*!
3646  * \brief True if the field can be converted into the med field
3647  */
3648 //================================================================================
3649
3650 bool DoubleField::isMedCompatible(bool& sameNbGauss) const
3651 {
3652   for ( unsigned int iSub = 0; iSub < _sub.size(); ++iSub )
3653     {
3654       if ( !getSupport(iSub) || !getSupport(iSub)->_medGroup )
3655         THROW_IK_EXCEPTION("SauvReader INTERNAL ERROR: NULL field support");
3656
3657       sameNbGauss = true;
3658       if ( !_sub[iSub].isSameNbGauss() )
3659         {
3660           std::cout << "Field <" << _name << "> : different nb of gauss points in components" << std::endl;
3661           sameNbGauss = false;
3662           //return false;
3663         }
3664     }
3665   return true;
3666 }
3667
3668 //================================================================================
3669 /*!
3670  * \brief return true if all sub-components has same components and same nbGauss
3671  */
3672 //================================================================================
3673
3674 bool DoubleField::hasSameComponentsBySupport() const
3675 {
3676   std::vector< _Sub_data >::const_iterator sub_data = _sub.begin();
3677   const _Sub_data& first_sub_data = *sub_data;
3678   for ( ++sub_data ; sub_data != _sub.end(); ++sub_data )
3679     {
3680       if ( first_sub_data._comp_names != sub_data->_comp_names )
3681         return false; // diff names of components
3682
3683       if ( first_sub_data._nb_gauss != sub_data->_nb_gauss &&
3684            first_sub_data._support->_cellType == sub_data->_support->_cellType)
3685         return false; // diff nb of gauss points on same cell type
3686     }
3687   return true;
3688 }
3689
3690 //================================================================================
3691 /*!
3692  * \brief Return type of MEDCouplingFieldDouble
3693  */
3694 //================================================================================
3695
3696 MEDCoupling::TypeOfField DoubleField::getMedType( const int iSub ) const
3697 {
3698   using namespace INTERP_KERNEL;
3699
3700   const Group* grp = hasCommonSupport() ? _group : _sub[iSub]._support;
3701   if ( _sub[iSub].nbGauss() > 1 )
3702     {
3703       const CellModel& cm = CellModel::GetCellModel( _sub[iSub]._support->_cellType );
3704       return (int) cm.getNumberOfNodes() == _sub[iSub].nbGauss() ? ON_GAUSS_NE : ON_GAUSS_PT;
3705     }
3706   else
3707     {
3708       return getDim( grp ) == 0 ? ON_NODES : ON_CELLS;
3709     }
3710 }
3711
3712 //================================================================================
3713 /*!
3714  * \brief Return TypeOfTimeDiscretization
3715  */
3716 //================================================================================
3717
3718 MEDCoupling::TypeOfTimeDiscretization DoubleField::getMedTimeDisc() const
3719 {
3720   return ONE_TIME;
3721   // NO_TIME = 4,
3722   // ONE_TIME = 5,
3723   // LINEAR_TIME = 6,
3724   // CONST_ON_TIME_INTERVAL = 7
3725 }
3726
3727 //================================================================================
3728 /*!
3729  * \brief Return nb tuples to be used to allocate DataArrayDouble
3730  */
3731 //================================================================================
3732
3733 mcIdType DoubleField::getNbTuples( const int iSub ) const
3734 {
3735   mcIdType nb = 0;
3736   if ( hasCommonSupport() && !_group->_groups.empty() )
3737     for ( size_t i = 0; i < _group->_groups.size(); ++i )
3738       nb += _sub[i].nbGauss() * _sub[i]._support->size();
3739   else
3740     nb = _sub[iSub].nbGauss() * getSupport(iSub)->size();
3741   return nb;
3742 }
3743
3744 //================================================================================
3745 /*!
3746  * \brief Store values of a sub-component and return nb of elements in the iSub
3747  */
3748 //================================================================================
3749
3750 mcIdType DoubleField::setValues( double * valPtr, const int iSub, const mcIdType elemShift ) const
3751 {
3752   // find values for iSub
3753   int iComp = 0;
3754   for ( int iS = 0; iS < iSub; ++iS )
3755     iComp += _sub[iS].nbComponents();
3756   const std::vector< double > * compValues = &_comp_values[ iComp ];
3757
3758   // Set values
3759
3760   const std::vector< mcIdType >& relocTable = getSupport( iSub )->_relocTable;
3761
3762   const mcIdType nbElems = _sub[iSub]._support->size();
3763   const int nbGauss      = _sub[iSub].nbGauss();
3764   const int nbComponents = _sub[iSub].nbComponents();
3765   const int nbValsByElem = nbComponents * nbGauss;
3766
3767   // check nb values
3768   mcIdType nbVals = 0;
3769   for ( iComp = 0; iComp < nbComponents; ++iComp )
3770     nbVals += ToIdType( compValues[iComp].size() );
3771   const bool isConstField = ( nbVals == nbComponents ); // one value per component (issue 22321)
3772   if ( !isConstField && nbVals != nbElems * nbValsByElem )
3773     THROW_IK_EXCEPTION("SauvMedConvertor.cxx: support size mismatches field size");
3774
3775   // compute nb values in previous subs
3776   mcIdType valsShift = 0;
3777   for ( mcIdType iS = iSub-1, shift = elemShift; shift > 0; --iS)
3778     {
3779       mcIdType nbE = _sub[iS]._support->size();
3780       shift -= nbE;
3781       valsShift += nbE * _sub[iS].nbComponents() * _sub[iS].nbGauss();
3782     }
3783
3784   if ( isConstField )
3785     for ( mcIdType iE = 0; iE < nbElems; ++iE )
3786       {
3787         mcIdType iMed = valsShift + nbValsByElem * ( relocTable.empty() ? iE : relocTable[iE+elemShift]-elemShift );
3788         for ( iComp = 0; iComp < nbComponents; ++iComp )
3789           valPtr[ iMed + iComp ] = compValues[iComp][ 0 ];
3790       }
3791   else
3792     for ( mcIdType iE = 0; iE < nbElems; ++iE )
3793       {
3794         mcIdType iMed = valsShift + nbValsByElem * ( relocTable.empty() ? iE : relocTable[iE+elemShift]-elemShift );
3795         for ( iComp = 0; iComp < nbComponents; ++iComp )
3796           for ( int iG = 0; iG < nbGauss; ++iG )
3797             valPtr[ iMed + iG * nbComponents + iComp ] = compValues[iComp][ iE * nbGauss + iG ];
3798       }
3799   return nbElems;
3800 }
3801
3802 //================================================================================
3803 /*!
3804  * \brief Destructor of IntermediateMED
3805  */
3806 //================================================================================
3807
3808 IntermediateMED::~IntermediateMED()
3809 {
3810   for ( size_t i = 0; i < _nodeFields.size(); ++i )
3811     if ( _nodeFields[i] )
3812       delete _nodeFields[i];
3813   _nodeFields.clear();
3814
3815   for ( size_t i = 0; i < _cellFields.size(); ++i )
3816     if ( _cellFields[i] )
3817       delete _cellFields[i];
3818   _cellFields.clear();
3819
3820   for ( size_t i = 0; i < _groups.size(); ++i )
3821     if ( _groups[i]._medGroup )
3822       _groups[i]._medGroup->decrRef();
3823 }
3824
3825 //================================================================================
3826 /*!
3827  * \brief CellsByDimIterator constructor
3828  */
3829 CellsByDimIterator::CellsByDimIterator( const IntermediateMED & medi, int dimm)
3830 {
3831   myImed = & medi;
3832   init( dimm );
3833 }
3834 /*!
3835  * \brief Initialize iteration on cells of given dimension
3836  */
3837 void CellsByDimIterator::init(const int  dimm)
3838 {
3839   myCurType = -1;
3840   myTypeEnd = INTERP_KERNEL::NORM_HEXA20 + 1;
3841   myDim = dimm;
3842 }
3843 /*!
3844  * \brief return next set of Cell's of required dimension
3845  */
3846 const std::set< Cell > * CellsByDimIterator::nextType()
3847 {
3848   while ( ++myCurType < myTypeEnd )
3849     if ( !myImed->_cellsByType[myCurType].empty() && ( myDim < 0 || dim(false) == myDim ))
3850       return & myImed->_cellsByType[myCurType];
3851   return 0;
3852 }
3853 /*!
3854  * \brief return dimension of cells returned by the last or further next()
3855  */
3856 int CellsByDimIterator::dim(const bool last) const
3857 {
3858   int typp = myCurType;
3859   if ( !last )
3860     while ( typp < myTypeEnd && myImed->_cellsByType[typp].empty() )
3861       ++typp;
3862   return typp < myTypeEnd ? getDimension( TCellType( typp )) : 4;
3863 }
3864 // END CellsByDimIterator ========================================================