src/MEDLoader/SauvMedConvertor.cxx

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