src/MEDLoader/SauvMedConvertor.cxx

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