src/MEDLoader/SauvMedConvertor.cxx

   1 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19 // File      : SauvMedConvertor.cxx
  20 // Created   : Tue Aug 16 14:43:20 2011
  21 // Author    : Edward AGAPOV (eap)
  22 //
  23
  24 #include "SauvMedConvertor.hxx"
  25
  26 #include "CellModel.hxx"
  27 #include "MEDFileMesh.hxx"
  28 #include "MEDFileField.hxx"
  29 #include "MEDFileData.hxx"
  30 #include "MEDCouplingFieldDouble.hxx"
  31
  32 #include <iostream>
  33 #include <cassert>
  34 #include <cmath>
  35 #include <queue>
  36 #include <limits>
  37
  38 #include <cstdlib>
  39 #include <cstring>
  40 #include <fcntl.h>
  41
  42 #ifdef WNT
  43 #include <io.h>
  44 #endif
  45
  46 #ifndef WNT
  47 #define HAS_XDR
  48 #include <unistd.h>
  49 #endif
  50
  51 #ifdef HAS_XDR
  52 #include <rpc/xdr.h>
  53 #endif
  54
  55 using namespace SauvUtilities;
  56 using namespace ParaMEDMEM;
  57 using namespace std;
  58
  59 namespace
  60 {
  61   // for ASCII file reader
  62   const int GIBI_MaxOutputLen = 150; // max length of a line in the sauve file
  63   const int GIBI_BufferSize   = 16184; // for buffered reading
  64
  65   using namespace INTERP_KERNEL;
  66
  67   const size_t MaxMedCellType = NORM_ERROR;
  68   const size_t NbGibiCellTypes = 47;
  69   const TCellType GibiTypeToMed[NbGibiCellTypes] =
  70     {
  71       /*1 */ NORM_POINT1 ,/*2 */ NORM_SEG2   ,/*3 */ NORM_SEG3   ,/*4 */ NORM_TRI3   ,/*5 */ NORM_ERROR  ,
  72       /*6 */ NORM_TRI6   ,/*7 */ NORM_ERROR  ,/*8 */ NORM_QUAD4  ,/*9 */ NORM_ERROR  ,/*10*/ NORM_QUAD8  ,
  73       /*11*/ NORM_ERROR  ,/*12*/ NORM_ERROR  ,/*13*/ NORM_ERROR  ,/*14*/ NORM_HEXA8  ,/*15*/ NORM_HEXA20 ,
  74       /*16*/ NORM_PENTA6 ,/*17*/ NORM_PENTA15,/*18*/ NORM_ERROR  ,/*19*/ NORM_ERROR  ,/*20*/ NORM_ERROR  ,
  75       /*21*/ NORM_ERROR  ,/*22*/ NORM_ERROR  ,/*23*/ NORM_TETRA4 ,/*24*/ NORM_TETRA10,/*25*/ NORM_PYRA5  ,
  76       /*26*/ NORM_PYRA13 ,/*27*/ NORM_ERROR  ,/*28*/ NORM_ERROR  ,/*29*/ NORM_ERROR  ,/*30*/ NORM_ERROR  ,
  77       /*31*/ NORM_ERROR  ,/*32*/ NORM_ERROR  ,/*33*/ NORM_ERROR  ,/*34*/ NORM_ERROR  ,/*35*/ NORM_ERROR  ,
  78       /*36*/ NORM_ERROR  ,/*37*/ NORM_ERROR  ,/*38*/ NORM_ERROR  ,/*39*/ NORM_ERROR  ,/*40*/ NORM_ERROR  ,
  79       /*41*/ NORM_ERROR  ,/*42*/ NORM_ERROR  ,/*43*/ NORM_ERROR  ,/*44*/ NORM_ERROR  ,/*45*/ NORM_ERROR  ,
  80       /*46*/ NORM_ERROR  ,/*47*/ NORM_ERROR
  81     };
  82
  83   //================================================================================
  84   /*!
  85    * \brief Return dimension of a group
  86    */
  87   //================================================================================
  88
  89   unsigned getDim( const Group* grp )
  90   {
  91     return SauvUtilities::getDimension( grp->_groups.empty() ? grp->_cellType : grp->_groups[0]->_cellType );
  92   }
  93
  94   //================================================================================
  95   /*!
  96    * \brief Converts connectivity of quadratic elements
  97    */
  98   //================================================================================
  99
 100   inline void ConvertQuadratic( const INTERP_KERNEL::NormalizedCellType type,
 101                                 const Cell &                            aCell )
 102   {
 103     if ( const int * conn = getGibi2MedQuadraticInterlace( type ))
 104       {
 105         Cell* ma = const_cast<Cell*>(&aCell);
 106         vector< Node* > new_nodes( ma->_nodes.size() );
 107         for ( size_t i = 0; i < new_nodes.size(); ++i )
 108           new_nodes[ i ] = ma->_nodes[ conn[ i ]];
 109         ma->_nodes.swap( new_nodes );
 110       }
 111   }
 112
 113   //================================================================================
 114   /*!
 115    * \brief Returns a vector of pairs of node indices to inverse a med volume element
 116    */
 117   //================================================================================
 118
 119   void getReverseVector (const INTERP_KERNEL::NormalizedCellType type,
 120                          vector<pair<int,int> > &                swapVec )
 121   {
 122     swapVec.clear();
 123
 124     switch ( type )
 125       {
 126       case NORM_TETRA4:
 127         swapVec.resize(1);
 128         swapVec[0] = make_pair( 1, 2 );
 129         break;
 130       case NORM_PYRA5:
 131         swapVec.resize(1);
 132         swapVec[0] = make_pair( 1, 3 );
 133         break;
 134       case NORM_PENTA6:
 135         swapVec.resize(2);
 136         swapVec[0] = make_pair( 1, 2 );
 137         swapVec[1] = make_pair( 4, 5 );
 138         break;
 139       case NORM_HEXA8:
 140         swapVec.resize(2);
 141         swapVec[0] = make_pair( 1, 3 );
 142         swapVec[1] = make_pair( 5, 7 );
 143         break;
 144       case NORM_TETRA10:
 145         swapVec.resize(3);
 146         swapVec[0] = make_pair( 1, 2 );
 147         swapVec[1] = make_pair( 4, 6 );
 148         swapVec[2] = make_pair( 8, 9 );
 149         break;
 150       case NORM_PYRA13:
 151         swapVec.resize(4);
 152         swapVec[0] = make_pair( 1, 3 );
 153         swapVec[1] = make_pair( 5, 8 );
 154         swapVec[2] = make_pair( 6, 7 );
 155         swapVec[3] = make_pair( 10, 12 );
 156         break;
 157       case NORM_PENTA15:
 158         swapVec.resize(4);
 159         swapVec[0] = make_pair( 1, 2 );
 160         swapVec[1] = make_pair( 4, 5 );
 161         swapVec[2] = make_pair( 6, 8 );
 162         swapVec[3] = make_pair( 9, 11 );
 163         break;
 164       case NORM_HEXA20:
 165         swapVec.resize(7);
 166         swapVec[0] = make_pair( 1, 3 );
 167         swapVec[1] = make_pair( 5, 7 );
 168         swapVec[2] = make_pair( 8, 11 );
 169         swapVec[3] = make_pair( 9, 10 );
 170         swapVec[4] = make_pair( 12, 15 );
 171         swapVec[5] = make_pair( 13, 14 );
 172         swapVec[6] = make_pair( 17, 19 );
 173         break;
 174         //   case NORM_SEG3: no need to reverse edges
 175         //     swapVec.resize(1);
 176         //     swapVec[0] = make_pair( 1, 2 );
 177         //     break;
 178       case NORM_TRI6:
 179         swapVec.resize(2);
 180         swapVec[0] = make_pair( 1, 2 );
 181         swapVec[1] = make_pair( 3, 5 );
 182         break;
 183       case NORM_QUAD8:
 184         swapVec.resize(3);
 185         swapVec[0] = make_pair( 1, 3 );
 186         swapVec[1] = make_pair( 4, 7 );
 187         swapVec[2] = make_pair( 5, 6 );
 188         break;
 189       default:;
 190       }
 191   }
 192
 193   //================================================================================
 194   /*!
 195    * \brief Inverses element orientation using vector of indices to swap
 196    */
 197   //================================================================================
 198
 199   inline void reverse(const Cell & aCell, const vector<pair<int,int> > & swapVec )
 200   {
 201     Cell* ma = const_cast<Cell*>(&aCell);
 202     for ( unsigned i = 0; i < swapVec.size(); ++i )
 203       std::swap( ma->_nodes[ swapVec[i].first ],
 204                  ma->_nodes[ swapVec[i].second ]);
 205     if ( swapVec.empty() )
 206       ma->_reverse = true;
 207     else
 208       ma->_reverse = false;
 209   }
 210   //================================================================================
 211   /*!
 212    * \brief Comparator of cells by number used for ordering cells within a med group
 213    */
 214   struct TCellByIDCompare
 215   {
 216     bool operator () (const Cell* i1, const Cell* i2)
 217     {
 218       return i1->_number < i2->_number;
 219     }
 220   };
 221   typedef map< const Cell*, unsigned, TCellByIDCompare > TCellToOrderMap;
 222
 223   //================================================================================
 224   /*!
 225    * \brief Fill Group::_relocTable if necessary
 226    */
 227   //================================================================================
 228
 229   void setRelocationTable( Group* grp, TCellToOrderMap& cell2order )
 230   {
 231     if ( !grp->_isProfile ) return;
 232
 233     // check if relocation table is necessary
 234     bool isRelocated = false;
 235     unsigned newOrder = 0;
 236     TCellToOrderMap::iterator c2oIt = cell2order.begin(), c2oEnd = cell2order.end();
 237     for ( ; !isRelocated && c2oIt != c2oEnd; ++c2oIt, ++newOrder )
 238       isRelocated = ( c2oIt->second != newOrder );
 239
 240     if ( isRelocated )
 241     {
 242       grp->_relocTable.resize( cell2order.size() );
 243       for ( newOrder = 0, c2oIt = cell2order.begin(); c2oIt != c2oEnd; ++c2oIt, ++newOrder )
 244         grp->_relocTable[ c2oIt->second ] = newOrder;
 245     }
 246   }
 247 }
 248
 249 //================================================================================
 250 /*!
 251  * \brief Return dimension for the given cell type
 252  */
 253 //================================================================================
 254
 255 unsigned SauvUtilities::getDimension( INTERP_KERNEL::NormalizedCellType type )
 256 {
 257   return type == NORM_ERROR ? -1 : INTERP_KERNEL::CellModel::GetCellModel( type ).getDimension();
 258 }
 259
 260 //================================================================================
 261 /*!
 262  * \brief Returns interlace array to transform a quadratic GIBI element to a MED one
 263  */
 264 //================================================================================
 265
 266 const int * SauvUtilities::getGibi2MedQuadraticInterlace( INTERP_KERNEL::NormalizedCellType type )
 267 {
 268   static vector<const int*> conn;
 269   static const int hexa20 [] = {0,6,4,2, 12,18,16,14, 7,5,3,1, 19,17,15,13, 8,11,10,9};
 270   static const int penta15[] = {0,2,4, 9,11,13, 1,3,5, 10,12,14, 6,7,3};
 271   static const int pyra13 [] = {0,2,4,6, 12, 1,3,5,7, 8,9,10,11};
 272   static const int tetra10[] = {0,2,4, 9, 1,3,5, 6,7,8};
 273   static const int quad8  [] = {0,2,4,6, 1,3,5,7};
 274   static const int tria6  [] = {0,2,4, 1,3,5};
 275   static const int seg3   [] = {0,2,1};
 276   if ( conn.empty() )
 277   {
 278     conn.resize( MaxMedCellType + 1, 0 );
 279     conn[ NORM_HEXA20 ] = hexa20;
 280     conn[ NORM_PENTA15] = penta15;
 281     conn[ NORM_PYRA13 ] = pyra13;
 282     conn[ NORM_TETRA10] = tetra10;
 283     conn[ NORM_SEG3   ] = seg3;
 284     conn[ NORM_TRI6   ] = tria6;
 285     conn[ NORM_QUAD8  ] = quad8;
 286   }
 287   return conn[ type ];
 288 }
 289
 290 //================================================================================
 291 /*!
 292  * \brief Avoid coping sortedNodeIDs
 293  */
 294 //================================================================================
 295
 296 Cell::Cell(const Cell& ma)
 297   : _nodes(ma._nodes), _reverse(ma._reverse), _sortedNodeIDs(0), _number(ma._number)
 298 {
 299   if ( ma._sortedNodeIDs )
 300     {
 301       _sortedNodeIDs = new int[ _nodes.size() ];
 302       std::copy( ma._sortedNodeIDs, ma._sortedNodeIDs + _nodes.size(), _sortedNodeIDs );
 303     }
 304 }
 305
 306 //================================================================================
 307 /*!
 308  * \brief Rerturn the i-th link of face
 309  */
 310 //================================================================================
 311
 312 SauvUtilities::Link Cell::link(int i) const
 313 {
 314   int i2 = ( i + 1 ) % _nodes.size();
 315   if ( _reverse )
 316     return make_pair( _nodes[i2]->_number, _nodes[i]->_number );
 317   else
 318     return make_pair( _nodes[i]->_number, _nodes[i2]->_number );
 319 }
 320
 321 //================================================================================
 322 /*!
 323  * \brief creates if needed and return _sortedNodeIDs
 324  */
 325 //================================================================================
 326
 327 const TID* Cell::getSortedNodes() const
 328 {
 329   if ( !_sortedNodeIDs )
 330   {
 331     size_t l=_nodes.size();
 332     _sortedNodeIDs = new int[ l ];
 333
 334     for (size_t i=0; i!=l; ++i)
 335       _sortedNodeIDs[i]=_nodes[i]->_number;
 336     std::sort( _sortedNodeIDs, _sortedNodeIDs + l );
 337   }
 338   return _sortedNodeIDs;
 339 }
 340
 341 //================================================================================
 342 /*!
 343  * \brief Compare sorted ids of cell nodes
 344  */
 345 //================================================================================
 346
 347 bool Cell::operator< (const Cell& ma) const
 348 {
 349   if ( _nodes.size() == 1 )
 350     return _nodes[0] < ma._nodes[0];
 351
 352   const int* v1 = getSortedNodes();
 353   const int* v2 = ma.getSortedNodes();
 354   for ( const int* vEnd = v1 + _nodes.size(); v1 < vEnd; ++v1, ++v2 )
 355     if(*v1 != *v2)
 356       return *v1 < *v2;
 357   return false;
 358 }
 359
 360 //================================================================================
 361 /*!
 362  * \brief Dump a Cell
 363  */
 364 //================================================================================
 365
 366 std::ostream& SauvUtilities::operator<< (std::ostream& os, const SauvUtilities::Cell& ma)
 367 {
 368   os << "cell " << ma._number << " (" << ma._nodes.size() << " nodes) : < " << ma._nodes[0]->_number;
 369   for( size_t i=1; i!=ma._nodes.size(); ++i)
 370     os << ", " << ma._nodes[i]->_number;
 371 #ifdef _DEBUG_
 372   os << " > sortedNodes: ";
 373   if ( ma._sortedNodeIDs ) {
 374     os << "< ";
 375     for( size_t i=0; i!=ma._nodes.size(); ++i)
 376       os << ( i ? ", " : "" ) << ma._sortedNodeIDs[i];
 377     os << " >";
 378   }
 379   else {
 380     os << "NULL";
 381   }
 382 #endif
 383   return os;
 384 }
 385
 386 //================================================================================
 387 /*!
 388  * \brief Return nb of elements in the group
 389  */
 390 //================================================================================
 391
 392 int Group::size() const
 393 {
 394   int sizze = 0;
 395   if ( !_relocTable.empty() )
 396     sizze =  _relocTable.size();
 397   else if ( _medGroup )
 398     sizze = _medGroup->getNumberOfTuples();
 399   else if ( !_cells.empty() )
 400     sizze = _cells.size();
 401   else
 402     for ( size_t i = 0; i < _groups.size(); ++i )
 403       sizze += _groups[i]->size();
 404   return sizze;
 405 }
 406
 407 //================================================================================
 408 /*!
 409  * \brief Conver gibi element type to med one
 410  */
 411 //================================================================================
 412
 413 INTERP_KERNEL::NormalizedCellType SauvUtilities::gibi2medGeom( size_t gibiType )
 414 {
 415   if ( gibiType < 1 || gibiType > NbGibiCellTypes )
 416     return NORM_ERROR;
 417
 418   return GibiTypeToMed[ gibiType - 1 ];
 419 }
 420
 421 //================================================================================
 422 /*!
 423  * \brief Conver med element type to gibi one
 424  */
 425 //================================================================================
 426
 427 int SauvUtilities::med2gibiGeom( INTERP_KERNEL::NormalizedCellType medGeomType )
 428 {
 429   for ( unsigned int i = 0; i < NbGibiCellTypes; i++ )
 430     if ( GibiTypeToMed[ i ] == medGeomType )
 431       return i + 1;
 432
 433   return -1;
 434 }
 435
 436 //================================================================================
 437 /*!
 438  * \brief Remember the file name
 439  */
 440 //================================================================================
 441
 442 FileReader::FileReader(const char* fileName):_fileName(fileName),_iRead(0),_nbToRead(0)
 443 {
 444 }
 445
 446 //================================================================================
 447 /*!
 448  * \brief Constructor of ASCII sauve file reader
 449  */
 450 //================================================================================
 451
 452 ASCIIReader::ASCIIReader(const char* fileName)
 453   :FileReader(fileName),
 454    _file(-1)
 455 {
 456 }
 457
 458 //================================================================================
 459 /*!
 460  * \brief Return true
 461  */
 462 //================================================================================
 463
 464 bool ASCIIReader::isASCII() const
 465 {
 466   return true;
 467 }
 468
 469 //================================================================================
 470 /*!
 471  * \brief Try to open an ASCII file
 472  */
 473 //================================================================================
 474
 475 bool ASCIIReader::open()
 476 {
 477 #ifdef WNT
 478   _file = ::_open (_fileName.c_str(), _O_RDONLY|_O_BINARY);
 479 #else
 480   _file = ::open (_fileName.c_str(), O_RDONLY);
 481 #endif
 482   if (_file >= 0)
 483     {
 484       _start  = new char [GIBI_BufferSize]; // working buffer beginning
 485       //_tmpBuf = new char [GIBI_MaxOutputLen];
 486       _ptr    = _start;
 487       _eptr   = _start;
 488       _lineNb = 0;
 489     }
 490   else
 491     {
 492       //THROW_IK_EXCEPTION("Can't open file "<<_fileName << " fd: " << _file);
 493     }
 494   return (_file >= 0);
 495 }
 496
 497 //================================================================================
 498 /*!
 499  * \brief Close the file
 500  */
 501 //================================================================================
 502
 503 ASCIIReader::~ASCIIReader()
 504 {
 505   if (_file >= 0)
 506     {
 507       ::close (_file);
 508       if (_start != 0L)
 509         {
 510           delete [] _start;
 511           //delete [] _tmpBuf;
 512           _start = 0;
 513         }
 514       _file = -1;
 515     }
 516 }
 517
 518 //================================================================================
 519 /*!
 520  * \brief Return a next line of the file
 521  */
 522 //================================================================================
 523
 524 bool ASCIIReader::getNextLine (char* & line, bool raiseOEF /*= true*/ )
 525 {
 526   if ( getLine( line )) return true;
 527   if ( raiseOEF )
 528     THROW_IK_EXCEPTION("Unexpected EOF on ln "<<_lineNb);
 529   return false;
 530 }
 531
 532 //================================================================================
 533 /*!
 534  * \brief Read a next line of the file if necessary
 535  */
 536 //================================================================================
 537
 538 bool ASCIIReader::getLine(char* & line)
 539 {
 540   bool aResult = true;
 541   // Check the state of the buffer;
 542   // if there is too little left, read the next portion of data
 543   int nBytesRest = _eptr - _ptr;
 544   if (nBytesRest < GIBI_MaxOutputLen)
 545     {
 546       if (nBytesRest > 0)
 547         {
 548           // move the remaining portion to the buffer beginning
 549           for ( int i = 0; i < nBytesRest; ++i )
 550             _start[i] = _ptr[i];
 551           //memcpy (_tmpBuf, _ptr, nBytesRest);
 552           //memcpy (_start, _tmpBuf, nBytesRest);
 553         }
 554       else
 555         {
 556           nBytesRest = 0;
 557         }
 558       _ptr = _start;
 559       const int nBytesRead = ::read (_file,
 560                                      &_start [nBytesRest],
 561                                      GIBI_BufferSize - nBytesRest);
 562       nBytesRest += nBytesRead;
 563       _eptr = &_start [nBytesRest];
 564     }
 565   // Check the buffer for the end-of-line
 566   char * ptr = _ptr;
 567   while (true)
 568     {
 569       // Check for end-of-the-buffer, the ultimate criterion for termination
 570       if (ptr >= _eptr)
 571         {
 572           if (nBytesRest <= 0)
 573             aResult = false;
 574           else
 575             _eptr[-1] = '\0';
 576           break;
 577         }
 578       // seek the line-feed character
 579       if (ptr[0] == '\n')
 580         {
 581           if (ptr[-1] == '\r')
 582             ptr[-1] = '\0';
 583           ptr[0] = '\0';
 584           ++ptr;
 585           break;
 586         }
 587       ++ptr;
 588     }
 589   // Output the result
 590   line = _ptr;
 591   _ptr = ptr;
 592   _lineNb++;
 593
 594   return aResult;
 595 }
 596
 597 //================================================================================
 598 /*!
 599  * \brief Prepare for iterating over given nb of values
 600  *  \param nbToRead - nb of fields to read
 601  *  \param nbPosInLine - nb of fields in one line
 602  *  \param width - field width
 603  *  \param shift - shift from line beginning to the field start
 604  */
 605 //================================================================================
 606
 607 void ASCIIReader::init( int nbToRead, int nbPosInLine, int width, int shift /*= 0*/ )
 608 {
 609   _nbToRead    = nbToRead;
 610   _nbPosInLine = nbPosInLine;
 611   _width       = width;
 612   _shift       = shift;
 613   _iPos = _iRead = 0;
 614   if ( _nbToRead )
 615     {
 616       getNextLine( _curPos );
 617       _curPos = _curPos + _shift;
 618     }
 619   else
 620     {
 621       _curPos = 0;
 622     }
 623   _curLocale.clear();
 624 }
 625
 626 //================================================================================
 627 /*!
 628  * \brief Prepare for iterating over given nb of string values
 629  */
 630 //================================================================================
 631
 632 void ASCIIReader::initNameReading(int nbValues, int width /*= 8*/)
 633 {
 634   init( nbValues, 72 / ( width + 1 ), width, 1 );
 635 }
 636
 637 //================================================================================
 638 /*!
 639  * \brief Prepare for iterating over given nb of integer values
 640  */
 641 //================================================================================
 642
 643 void ASCIIReader::initIntReading(int nbValues)
 644 {
 645   init( nbValues, 10, 8 );
 646 }
 647
 648 //================================================================================
 649 /*!
 650  * \brief Prepare for iterating over given nb of real values
 651  */
 652 //================================================================================
 653
 654 void ASCIIReader::initDoubleReading(int nbValues)
 655 {
 656   init( nbValues, 3, 22 );
 657
 658   // Correction 2 of getDouble(): set "C" numeric locale to read numbers
 659   // with dot decimal point separator, as it is in SAUVE files
 660   _curLocale = setlocale(LC_NUMERIC, "C");
 661 }
 662
 663 //================================================================================
 664 /*!
 665  * \brief Return true if not all values have been read
 666  */
 667 //================================================================================
 668
 669 bool ASCIIReader::more() const
 670 {
 671   bool result = false;
 672   if ( _iRead < _nbToRead)
 673     {
 674       if ( _curPos ) result = true;
 675     }
 676   return result;
 677 }
 678
 679 //================================================================================
 680 /*!
 681  * \brief Go to the nex value
 682  */
 683 //================================================================================
 684
 685 void ASCIIReader::next()
 686 {
 687   if ( !more() )
 688     THROW_IK_EXCEPTION("SauvUtilities::ASCIIReader::next(): no more() values to read");
 689   ++_iRead;
 690   ++_iPos;
 691   if ( _iRead < _nbToRead )
 692     {
 693       if ( _iPos >= _nbPosInLine )
 694         {
 695           getNextLine( _curPos );
 696           _curPos = _curPos + _shift;
 697           _iPos = 0;
 698         }
 699       else
 700         {
 701           _curPos = _curPos + _width + _shift;
 702         }
 703     }
 704   else
 705     {
 706       _curPos = 0;
 707       if ( !_curLocale.empty() )
 708         {
 709           setlocale(LC_NUMERIC, _curLocale.c_str());
 710           _curLocale.clear();
 711         }
 712     }
 713 }
 714
 715 //================================================================================
 716 /*!
 717  * \brief Return the current integer value
 718  */
 719 //================================================================================
 720
 721 int ASCIIReader::getInt() const
 722 {
 723   // fix for two glued ints (issue 0021009):
 724   // Line nb    |   File contents
 725   // ------------------------------------------------------------------------------------
 726   // 53619905   |       1       2       6       8
 727   // 53619906   |                                                                SCALAIRE
 728   // 53619907   |    -63312600499       1       0       0       0      -2       0       2
 729   //   where -63312600499 is actualy -633 and 12600499
 730   char hold=_curPos[_width];
 731   _curPos[_width] = '\0';
 732   int result = atoi( _curPos );
 733   _curPos[_width] = hold;
 734   return result;
 735   //return atoi(str());
 736 }
 737
 738 //================================================================================
 739 /*!
 740  * \brief Return the current float value
 741  */
 742 //================================================================================
 743
 744 float ASCIIReader::getFloat() const
 745 {
 746   return getDouble();
 747 }
 748
 749 //================================================================================
 750 /*!
 751  * \brief Return the current double value
 752  */
 753 //================================================================================
 754
 755 double ASCIIReader::getDouble() const
 756 {
 757   //std::string aStr (_curPos);
 758
 759   // Correction: add missing 'E' specifier
 760   // int aPosStart = aStr.find_first_not_of(" \t");
 761   // if (aPosStart < (int)aStr.length()) {
 762   //   int aPosSign = aStr.find_first_of("+-", aPosStart + 1); // pass the leading symbol, as it can be a sign
 763   //   if (aPosSign < (int)aStr.length()) {
 764   //     if (aStr[aPosSign - 1] != 'e' && aStr[aPosSign - 1] != 'E')
 765   //       aStr.insert(aPosSign, "E", 1);
 766   //   }
 767   // }
 768
 769   // Different Correction (more optimal)
 770   // Sample:
 771   //  0.00000000000000E+00 -2.37822406690632E+01  6.03062748797469E+01
 772   //  7.70000000000000-100  7.70000000000000+100  7.70000000000000+100
 773   //0123456789012345678901234567890123456789012345678901234567890123456789
 774   const size_t posE = 18;
 775   if ( _curPos[posE] != 'E' && _curPos[posE] != 'e' )
 776     {
 777       std::string aStr (_curPos);
 778       if ( aStr.size() < posE+1 )
 779         THROW_IK_EXCEPTION("No more doubles (line #" << lineNb() << ")");
 780       aStr.insert( posE, "E", 1 );
 781       return atof(aStr.c_str());
 782     }
 783   return atof( _curPos );
 784 }
 785
 786 //================================================================================
 787 /*!
 788  * \brief Return the current string value
 789  */
 790 //================================================================================
 791
 792 string ASCIIReader::getName() const
 793 {
 794   int len = _width;
 795   while (( _curPos[len-1] == ' ' || _curPos[len-1] == 0) && len > 0 )
 796     len--;
 797   return string( _curPos, len );
 798 }
 799
 800 //================================================================================
 801 /*!
 802  * \brief Constructor of a binary sauve file reader
 803  */
 804 //================================================================================
 805
 806 XDRReader::XDRReader(const char* fileName) :FileReader(fileName), _xdrs_file(NULL)
 807 {
 808 }
 809
 810 //================================================================================
 811 /*!
 812  * \brief Close the XDR sauve file
 813  */
 814 //================================================================================
 815
 816 XDRReader::~XDRReader()
 817 {
 818 #ifdef HAS_XDR
 819   if ( _xdrs_file )
 820     {
 821       xdr_destroy((XDR*)_xdrs);
 822       free((XDR*)_xdrs);
 823       ::fclose(_xdrs_file);
 824       _xdrs_file = NULL;
 825     }
 826 #endif
 827 }
 828
 829 //================================================================================
 830 /*!
 831  * \brief Return false
 832  */
 833 //================================================================================
 834
 835 bool XDRReader::isASCII() const
 836 {
 837   return false;
 838 }
 839
 840 //================================================================================
 841 /*!
 842  * \brief Try to open an XRD file
 843  */
 844 //================================================================================
 845
 846 bool XDRReader::open()
 847 {
 848   bool xdr_ok = false;
 849 #ifdef HAS_XDR
 850   if ((_xdrs_file = ::fopen(_fileName.c_str(), "r")))
 851     {
 852       _xdrs = (XDR *)malloc(sizeof(XDR));
 853       xdrstdio_create((XDR*)_xdrs, _xdrs_file, XDR_DECODE);
 854
 855       const int maxsize = 10;
 856       char icha[maxsize+1];
 857       char* icha2 = icha;
 858       if (( xdr_ok = xdr_string((XDR*)_xdrs, &icha2, maxsize)))
 859         {
 860           icha[maxsize] = '\0';
 861           xdr_ok = (strcmp(icha, "CASTEM XDR") == 0);
 862         }
 863       if ( !xdr_ok )
 864         {
 865           xdr_destroy((XDR*)_xdrs);
 866           free((XDR*)_xdrs);
 867           fclose(_xdrs_file);
 868           _xdrs_file = NULL;
 869         }
 870     }
 871 #endif
 872   return xdr_ok;
 873 }
 874
 875 //================================================================================
 876 /*!
 877  * \brief A stub
 878  */
 879 //================================================================================
 880
 881 bool XDRReader::getNextLine (char* &, bool )
 882 {
 883   return true;
 884 }
 885
 886 //================================================================================
 887 /*!
 888  * \brief Prepare for iterating over given nb of values
 889  *  \param nbToRead - nb of fields to read
 890  *  \param width - field width
 891  */
 892 //================================================================================
 893
 894 void XDRReader::init( int nbToRead, int width/*=0*/ )
 895 {
 896   if(_iRead < _nbToRead)
 897     {
 898       cout << "_iRead, _nbToRead : " << _iRead << " " << _nbToRead << endl;
 899       cout << "Unfinished iteration before new one !" << endl;
 900       THROW_IK_EXCEPTION("SauvUtilities::XDRReader::init(): Unfinished iteration before new one !");
 901     }
 902   _iRead    = 0;
 903   _nbToRead = nbToRead;
 904   _width    = width;
 905 }
 906
 907 //================================================================================
 908 /*!
 909  * \brief Prepare for iterating over given nb of string values
 910  */
 911 //================================================================================
 912
 913 void XDRReader::initNameReading(int nbValues, int width)
 914 {
 915   init( nbValues, width );
 916   _xdr_kind = _xdr_kind_char;
 917   if(nbValues*width)
 918     {
 919       unsigned int nels = nbValues*width;
 920       _xdr_cvals = (char*)malloc((nels+1)*sizeof(char));
 921 #ifdef HAS_XDR
 922       xdr_string((XDR*)_xdrs, &_xdr_cvals, nels);
 923 #endif
 924       _xdr_cvals[nels] = '\0';
 925     }
 926 }
 927
 928 //================================================================================
 929 /*!
 930  * \brief Prepare for iterating over given nb of integer values
 931  */
 932 //================================================================================
 933
 934 void XDRReader::initIntReading(int nbValues)
 935 {
 936   init( nbValues );
 937   _xdr_kind = _xdr_kind_int;
 938   if(nbValues)
 939     {
 940 #ifdef HAS_XDR
 941       unsigned int nels = nbValues;
 942       unsigned int actual_nels;
 943       _xdr_ivals = (int*)malloc(nels*sizeof(int));
 944       xdr_array((XDR*)_xdrs, (char **)&_xdr_ivals, &actual_nels, nels, sizeof(int), (xdrproc_t)xdr_int);
 945 #endif
 946     }
 947 }
 948
 949 //================================================================================
 950 /*!
 951  * \brief Prepare for iterating over given nb of real values
 952  */
 953 //================================================================================
 954
 955 void XDRReader::initDoubleReading(int nbValues)
 956 {
 957   init( nbValues );
 958   _xdr_kind = _xdr_kind_double;
 959   if(nbValues)
 960     {
 961 #ifdef HAS_XDR
 962       unsigned int nels = nbValues;
 963       unsigned int actual_nels;
 964       _xdr_dvals = (double*)malloc(nels*sizeof(double));
 965       xdr_array((XDR*)_xdrs, (char **)&_xdr_dvals, &actual_nels, nels, sizeof(double), (xdrproc_t)xdr_double);
 966 #endif
 967     }
 968 }
 969
 970 //================================================================================
 971 /*!
 972  * \brief Return true if not all values have been read
 973  */
 974 //================================================================================
 975
 976 bool XDRReader::more() const
 977 {
 978   return _iRead < _nbToRead;
 979 }
 980
 981 //================================================================================
 982 /*!
 983  * \brief Go to the nex value
 984  */
 985 //================================================================================
 986
 987 void XDRReader::next()
 988 {
 989   if ( !more() )
 990     THROW_IK_EXCEPTION("SauvUtilities::XDRReader::next(): no more() values to read");
 991
 992   ++_iRead;
 993   if ( _iRead < _nbToRead )
 994     {
 995     }
 996   else
 997     {
 998       if(_xdr_kind == _xdr_kind_char) free(_xdr_cvals);
 999       if(_xdr_kind == _xdr_kind_int) free(_xdr_ivals);
1000       if(_xdr_kind == _xdr_kind_double) free(_xdr_dvals);
1001       _xdr_kind = _xdr_kind_null;
1002     }
1003 }
1004
1005 //================================================================================
1006 /*!
1007  * \brief Return the current integer value
1008  */
1009 //================================================================================
1010
1011 int XDRReader::getInt() const
1012 {
1013   if(_iRead < _nbToRead)
1014     {
1015       return _xdr_ivals[_iRead];
1016     }
1017   else
1018     {
1019       int result = 0;
1020 #ifdef HAS_XDR
1021       xdr_int((XDR*)_xdrs, &result);
1022 #endif
1023       return result;
1024     }
1025 }
1026
1027 //================================================================================
1028 /*!
1029  * \brief Return the current float value
1030  */
1031 //================================================================================
1032
1033 float  XDRReader::getFloat() const
1034 {
1035   float result = 0;
1036 #ifdef HAS_XDR
1037   xdr_float((XDR*)_xdrs, &result);
1038 #endif
1039   return result;
1040 }
1041
1042 //================================================================================
1043 /*!
1044  * \brief Return the current double value
1045  */
1046 //================================================================================
1047
1048 double XDRReader::getDouble() const
1049 {
1050   if(_iRead < _nbToRead)
1051     {
1052       return _xdr_dvals[_iRead];
1053     }
1054   else
1055     {
1056       double result = 0;
1057 #ifdef HAS_XDR
1058       xdr_double((XDR*)_xdrs, &result);
1059 #endif
1060       return result;
1061     }
1062 }
1063
1064 //================================================================================
1065 /*!
1066  * \brief Return the current string value
1067  */
1068 //================================================================================
1069
1070 std::string XDRReader::getName() const
1071 {
1072   int len = _width;
1073   char* s = _xdr_cvals + _iRead*_width;
1074   while (( s[len-1] == ' ' || s[len-1] == 0) && len > 0 )
1075     len--;
1076   return string( s, len );
1077 }
1078
1079 //================================================================================
1080 /*!
1081  * \brief Throw an exception if not all needed data is present
1082  */
1083 //================================================================================
1084
1085 void IntermediateMED::checkDataAvailability() const throw(INTERP_KERNEL::Exception)
1086 {
1087   if ( _spaceDim == 0 )
1088     THROW_IK_EXCEPTION("Wrong file format"); // it is the first record in the sauve file
1089
1090   if ( _groups.empty() )
1091     THROW_IK_EXCEPTION("No elements have been read");
1092
1093   if ( _points.empty() || _nbNodes == 0 )
1094     THROW_IK_EXCEPTION("Nodes of elements are not filled");
1095
1096   if ( _coords.empty() )
1097     THROW_IK_EXCEPTION("Node coordinates are missing");
1098
1099   if ( _coords.size() < _nbNodes * _spaceDim )
1100     THROW_IK_EXCEPTION("Nodes and coordinates mismatch");
1101 }
1102
1103 //================================================================================
1104 /*!
1105  * \brief Makes ParaMEDMEM::MEDFileData from self
1106  */
1107 //================================================================================
1108
1109 ParaMEDMEM::MEDFileData* IntermediateMED::convertInMEDFileDS()
1110 {
1111   MEDCouplingAutoRefCountObjectPtr< MEDFileUMesh >  mesh   = makeMEDFileMesh();
1112   MEDCouplingAutoRefCountObjectPtr< MEDFileFields > fields = makeMEDFileFields(mesh);
1113
1114   MEDCouplingAutoRefCountObjectPtr< MEDFileMeshes > meshes = MEDFileMeshes::New();
1115   MEDCouplingAutoRefCountObjectPtr< MEDFileData >  medData = MEDFileData::New();
1116   meshes->pushMesh( mesh );
1117   medData->setMeshes( meshes );
1118   if ( fields ) medData->setFields( fields );
1119
1120   medData->incrRef();
1121   return medData;
1122 }
1123
1124 //================================================================================
1125 /*!
1126  * \brief Creates ParaMEDMEM::MEDFileUMesh from its data
1127  */
1128 //================================================================================
1129
1130 ParaMEDMEM::MEDFileUMesh* IntermediateMED::makeMEDFileMesh()
1131 {
1132   // check if all needed piles are present
1133   checkDataAvailability();
1134
1135   // set long names
1136   setGroupLongNames();
1137
1138   // fix element orientation
1139   if ( _spaceDim == 2 || _spaceDim == 1 )
1140     orientElements2D();
1141   else if ( _spaceDim == 3 )
1142     orientElements3D();
1143
1144   // process groups
1145   decreaseHierarchicalDepthOfSubgroups();
1146   eraseUselessGroups();
1147   //detectMixDimGroups();
1148
1149   // assign IDs
1150   _points.numberNodes();
1151   numberElements();
1152
1153   // make the med mesh
1154
1155   MEDFileUMesh* mesh = MEDFileUMesh::New();
1156
1157   DataArrayDouble *coords = getCoords();
1158   setConnectivity( mesh, coords );
1159   setGroups( mesh );
1160
1161   coords->decrRef();
1162
1163   if ( !mesh->getName() || strlen( mesh->getName() ) == 0 )
1164     mesh->setName( "MESH" );
1165
1166   return mesh;
1167 }
1168
1169 //================================================================================
1170 /*!
1171  * \brief Set long names to groups
1172  */
1173 //================================================================================
1174
1175 void IntermediateMED::setGroupLongNames()
1176 {
1177   // IMP 0020434: mapping GIBI names to MED names
1178   // set med names to objects (mesh, fields, support, group or other)
1179
1180   set<int> treatedGroups;
1181
1182   list<nameGIBItoMED>::iterator itGIBItoMED = _listGIBItoMED_mail.begin();
1183   for (; itGIBItoMED != _listGIBItoMED_mail.end(); itGIBItoMED++)
1184     {
1185       if ( (int)_groups.size() < itGIBItoMED->gibi_id ) continue;
1186
1187       SauvUtilities::Group & grp = _groups[itGIBItoMED->gibi_id - 1];
1188
1189       // if there are several names for grp then the 1st name is the name
1190       // of grp and the rest ones are names of groups referring grp (issue 0021311)
1191       const bool isRefName = !treatedGroups.insert( itGIBItoMED->gibi_id ).second;
1192       if ( !isRefName )
1193         {
1194           grp._name = _mapStrings[ itGIBItoMED->med_id ];
1195         }
1196       else if ( !grp._refNames.empty() && grp._refNames.back().empty() )
1197         {
1198           for ( unsigned i = 0; i < grp._refNames.size(); ++i )
1199             if ( grp._refNames[i].empty() )
1200               grp._refNames[i] = _mapStrings[ (*itGIBItoMED).med_id ];
1201         }
1202       else
1203         {
1204           grp._refNames.push_back( _mapStrings[ (*itGIBItoMED).med_id ]);
1205         }
1206     }
1207 }
1208
1209 //================================================================================
1210 /*!
1211  * \brief Set long names to fields
1212  */
1213 //================================================================================
1214
1215 void IntermediateMED::setFieldLongNames(set< string >& usedNames)
1216 {
1217   list<nameGIBItoMED>::iterator itGIBItoMED = _listGIBItoMED_cham.begin();
1218   for (; itGIBItoMED != _listGIBItoMED_cham.end(); itGIBItoMED++)
1219     {
1220       if (itGIBItoMED->gibi_pile == PILE_FIELD)
1221         {
1222           _cellFields[itGIBItoMED->gibi_id - 1]->_name = _mapStrings[itGIBItoMED->med_id];
1223         }
1224       else if (itGIBItoMED->gibi_pile == PILE_NODES_FIELD)
1225         {
1226           _nodeFields[itGIBItoMED->gibi_id - 1]->_name = _mapStrings[itGIBItoMED->med_id];
1227         }
1228     } // iterate on _listGIBItoMED_cham
1229
1230   for (itGIBItoMED =_listGIBItoMED_comp.begin(); itGIBItoMED != _listGIBItoMED_comp.end(); itGIBItoMED++)
1231     {
1232       string medName  = _mapStrings[itGIBItoMED->med_id];
1233       string gibiName = _mapStrings[itGIBItoMED->gibi_id];
1234
1235       bool name_found = false;
1236       for ( int isNodal = 0; isNodal < 2 && !name_found; ++isNodal )
1237         {
1238           vector<DoubleField* > & fields = isNodal ? _nodeFields : _cellFields;
1239           for ( size_t ifi = 0; ifi < fields.size() && !name_found; ifi++)
1240             {
1241               if (medName.find( fields[ifi]->_name + "." ) == 0 )
1242                 {
1243                   vector<DoubleField::_Sub_data>& aSubDs = fields[ifi]->_sub;
1244                   int nbSub = aSubDs.size();
1245                   for (int isu = 0; isu < nbSub; isu++)
1246                     for (int ico = 0; ico < aSubDs[isu].nbComponents(); ico++)
1247                       {
1248                         if (aSubDs[isu].compName(ico) == gibiName)
1249                           {
1250                             string medNameCompo = medName.substr( fields[ifi]->_name.size() + 1 );
1251                             fields[ifi]->_sub[isu].compName(ico) = medNameCompo;
1252                           }
1253                       }
1254                 }
1255             }
1256         }
1257     } // iterate on _listGIBItoMED_comp
1258
1259   for ( size_t i = 0; i < _nodeFields.size() ; i++)
1260     usedNames.insert( _nodeFields[i]->_name );
1261   for ( size_t i = 0; i < _cellFields.size() ; i++)
1262     usedNames.insert( _cellFields[i]->_name );
1263 }
1264
1265 //================================================================================
1266 /*!
1267  * \brief Decrease hierarchical depth of subgroups
1268  */
1269 //================================================================================
1270
1271 void IntermediateMED::decreaseHierarchicalDepthOfSubgroups()
1272 {
1273   for (size_t i=0; i!=_groups.size(); ++i)
1274   {
1275     Group& grp = _groups[i];
1276     for (size_t j = 0; j < grp._groups.size(); ++j )
1277     {
1278       Group & sub_grp = *grp._groups[j];
1279       if ( !sub_grp._groups.empty() )
1280       {
1281         // replace j with its 1st subgroup
1282         grp._groups[j] = sub_grp._groups[0];
1283         // push back the rest subs
1284         grp._groups.insert( grp._groups.end(), ++sub_grp._groups.begin(), sub_grp._groups.end() );
1285       }
1286     }
1287     // remove empty sub-_groups
1288     vector< Group* > newSubGroups;
1289     newSubGroups.reserve( grp._groups.size() );
1290     for (size_t j = 0; j < grp._groups.size(); ++j )
1291       if ( !grp._groups[j]->empty() )
1292         newSubGroups.push_back( grp._groups[j] );
1293     if ( newSubGroups.size() < grp._groups.size() )
1294       grp._groups.swap( newSubGroups );
1295   }
1296 }
1297
1298 //================================================================================
1299 /*!
1300  * \brief Erase _groups that won't be converted
1301  */
1302 //================================================================================
1303
1304 void IntermediateMED::eraseUselessGroups()
1305 {
1306   // propagate _isProfile=true to sub-groups of composite groups
1307   // for (size_t int i=0; i!=_groups.size(); ++i)
1308   // {
1309   //   Group* grp = _groups[i];
1310   //   if ( grp->_isProfile && !grp->_groups.empty() )
1311   //     for (size_t j = 0; j < grp->_groups.size(); ++j )
1312   //       grp->_groups[j]->_isProfile=true;
1313   // }
1314   std::set<Group*> groups2convert;
1315   // keep not named sub-groups of field supports
1316   for (size_t i=0; i!=_groups.size(); ++i)
1317   {
1318     Group& grp = _groups[i];
1319     if ( grp._isProfile && !grp._groups.empty() )
1320       groups2convert.insert( grp._groups.begin(), grp._groups.end() );
1321   }
1322
1323   // keep named groups and their subgroups
1324   for (size_t i=0; i!=_groups.size(); ++i)
1325   {
1326     Group& grp = _groups[i];
1327     if ( !grp._name.empty() && !grp.empty() )
1328     {
1329       groups2convert.insert( &grp );
1330       groups2convert.insert( grp._groups.begin(), grp._groups.end() );
1331     }
1332   }
1333   // erase groups that are not in groups2convert and not _isProfile
1334   for (size_t i=0; i!=_groups.size(); ++i)
1335   {
1336     Group* grp = &_groups[i];
1337     if ( !grp->_isProfile && !groups2convert.count( grp ) )
1338     {
1339       grp->_cells.clear();
1340       grp->_groups.clear();
1341     }
1342   }
1343 }
1344
1345 //================================================================================
1346 /*!
1347  * \brief Detect _groups of mixed dimension
1348  */
1349 //================================================================================
1350
1351 void IntermediateMED::detectMixDimGroups()
1352 {
1353   //hasMixedCells = false;
1354   for ( size_t i=0; i < _groups.size(); ++i )
1355   {
1356     Group& grp = _groups[i];
1357     if ( grp._groups.size() < 2 )
1358       continue;
1359
1360     // check if sub-groups have different dimension
1361     unsigned dim1 = getDim( &grp );
1362     for ( size_t j = 1; j  < grp._groups.size(); ++j )
1363     {
1364       unsigned dim2 = getDim( grp._groups[j] );
1365       if ( dim1 != dim2 )
1366       {
1367         grp._cells.clear();
1368         grp._groups.clear();
1369         if ( !grp._name.empty() )
1370           cout << "Erase a group with elements of different dim |" << grp._name << "|"<< endl;
1371         break;
1372       }
1373     }
1374   }
1375 }
1376
1377 //================================================================================
1378 /*!
1379  * \brief Fix connectivity of elements in 2D space
1380  */
1381 //================================================================================
1382
1383 void IntermediateMED::orientElements2D()
1384 {
1385   set<Cell>::const_iterator elemIt, elemEnd;
1386   vector< pair<int,int> > swapVec;
1387
1388   // ------------------------------------
1389   // fix connectivity of quadratic edges
1390   // ------------------------------------
1391   set<Cell>& quadEdges = _cellsByType[ INTERP_KERNEL::NORM_SEG3 ];
1392   if ( !quadEdges.empty() )
1393     {
1394       elemIt = quadEdges.begin(), elemEnd = quadEdges.end();
1395       for ( ; elemIt != elemEnd; ++elemIt )
1396         ConvertQuadratic( INTERP_KERNEL::NORM_SEG3, *elemIt );
1397     }
1398
1399   CellsByDimIterator faceIt( *this, 2 );
1400   while ( const set<Cell > * faces = faceIt.nextType() )
1401     {
1402       TCellType cellType = faceIt.type();
1403       bool isQuadratic = getGibi2MedQuadraticInterlace( cellType );
1404
1405       getReverseVector( cellType, swapVec );
1406
1407       // ------------------------------------
1408       // fix connectivity of quadratic faces
1409       // ------------------------------------
1410       if ( isQuadratic )
1411         for ( elemIt = faces->begin(), elemEnd = faces->end(); elemIt != elemEnd; elemIt++ )
1412           ConvertQuadratic( cellType, *elemIt );
1413
1414       // --------------------------
1415       // orient faces clockwise
1416       // --------------------------
1417       int iQuad = isQuadratic ? 2 : 1;
1418       for ( elemIt = faces->begin(), elemEnd = faces->end(); elemIt != elemEnd; elemIt++ )
1419         {
1420           // look for index of the most left node
1421           int iLeft = 0, iNode, nbNodes = elemIt->_nodes.size() / iQuad;
1422           double x, minX = nodeCoords( elemIt->_nodes[0] )[0];
1423           for ( iNode = 1; iNode < nbNodes; ++iNode )
1424             if (( x = nodeCoords( elemIt->_nodes[ iNode ])[ 0 ]) < minX )
1425               minX = x, iLeft = iNode;
1426
1427           // indeces of the nodes neighboring the most left one
1428           int iPrev = ( iLeft - 1 < 0 ) ? nbNodes - 1 : iLeft - 1;
1429           int iNext = ( iLeft + 1 == nbNodes ) ? 0 : iLeft + 1;
1430           // find components of prev-left and left-next vectors
1431           double xP = nodeCoords( elemIt->_nodes[ iPrev ])[ 0 ];
1432           double yP = nodeCoords( elemIt->_nodes[ iPrev ])[ 1 ];
1433           double xN = nodeCoords( elemIt->_nodes[ iNext ])[ 0 ];
1434           double yN = nodeCoords( elemIt->_nodes[ iNext ])[ 1 ];
1435           double xL = nodeCoords( elemIt->_nodes[ iLeft ])[ 0 ];
1436           double yL = nodeCoords( elemIt->_nodes[ iLeft ])[ 1 ];
1437           double xPL = xL - xP, yPL = yL - yP; // components of prev-left vector
1438           double xLN = xN - xL, yLN = yN - yL; // components of left-next vector
1439           // normalise y of the vectors
1440           double modPL = sqrt ( xPL * xPL + yPL * yPL );
1441           double modLN = sqrt ( xLN * xLN + yLN * yLN );
1442           if ( modLN > std::numeric_limits<double>::min() &&
1443                modPL > std::numeric_limits<double>::min() )
1444             {
1445               yPL /= modPL;
1446               yLN /= modLN;
1447               // summary direction of neighboring links must be positive
1448               bool clockwise = ( yPL + yLN > 0 );
1449               if ( !clockwise )
1450                 reverse( *elemIt, swapVec );
1451             }
1452         }
1453     }
1454 }
1455
1456 //================================================================================
1457 /*!
1458  * \brief Fix connectivity of elements in 3D space
1459  */
1460 //================================================================================
1461
1462 void IntermediateMED::orientElements3D()
1463 {
1464   // set _reverse flags of faces
1465   orientFaces3D();
1466
1467   // -----------------
1468   // fix connectivity
1469   // -----------------
1470
1471   set<Cell>::const_iterator elemIt, elemEnd;
1472   vector< pair<int,int> > swapVec;
1473
1474   for ( int dim = 1; dim <= 3; ++dim )
1475   {
1476     CellsByDimIterator cellsIt( *this, dim );
1477     while ( const set<Cell > * elems = cellsIt.nextType() )
1478     {
1479       TCellType cellType = cellsIt.type();
1480       bool isQuadratic = getGibi2MedQuadraticInterlace( cellType );
1481       getReverseVector( cellType, swapVec );
1482
1483       elemIt = elems->begin(), elemEnd = elems->end();
1484       for ( ; elemIt != elemEnd; elemIt++ )
1485       {
1486         // GIBI connectivity -> MED one
1487         if( isQuadratic )
1488           ConvertQuadratic( cellType, *elemIt );
1489
1490         // reverse faces
1491         if ( elemIt->_reverse )
1492           reverse ( *elemIt, swapVec );
1493       }
1494     }
1495   }
1496
1497   orientVolumes();
1498 }
1499
1500 //================================================================================
1501 /*!
1502  * \brief Orient equally (by setting _reverse flag) all connected faces in 3D space
1503  */
1504 //================================================================================
1505
1506 void IntermediateMED::orientFaces3D()
1507 {
1508   // fill map of links and their faces
1509   set<const Cell*> faces;
1510   map<const Cell*, Group*> fgm;
1511   map<Link, list<const Cell*> > linkFacesMap;
1512   map<Link, list<const Cell*> >::iterator lfIt, lfIt2;
1513
1514   for (size_t i=0; i!=_groups.size(); ++i)
1515     {
1516       Group& grp = _groups[i];
1517       if ( !grp._cells.empty() && getDimension( grp._cellType ) == 2 )
1518         for ( size_t j = 0; j < grp._cells.size(); ++j )
1519           if ( faces.insert( grp._cells[j] ).second )
1520             {
1521               for ( size_t k = 0; k < grp._cells[j]->_nodes.size(); ++k )
1522                 linkFacesMap[ grp._cells[j]->link( k ) ].push_back( grp._cells[j] );
1523               fgm.insert( make_pair( grp._cells[j], &grp ));
1524             }
1525     }
1526   // dump linkFacesMap
1527   //     for ( lfIt = linkFacesMap.begin(); lfIt!=linkFacesMap.end(); lfIt++) {
1528   //       cout<< "LINK: " << lfIt->first.first << "-" << lfIt->first.second << endl;
1529   //       list<const Cell*> & fList = lfIt->second;
1530   //       list<const Cell*>::iterator fIt = fList.begin();
1531   //       for ( ; fIt != fList.end(); fIt++ )
1532   //         cout << "\t" << **fIt << fgm[*fIt]->nom << endl;
1533   //     }
1534
1535   // Each oriented link must appear in one face only, else a face is reversed.
1536
1537   queue<const Cell*> faceQueue; /* the queue contains well oriented faces
1538                                      whose neighbors orientation is to be checked */
1539   bool manifold = true;
1540   while ( !linkFacesMap.empty() )
1541     {
1542       if ( faceQueue.empty() )
1543         {
1544           assert( !linkFacesMap.begin()->second.empty() );
1545           faceQueue.push( linkFacesMap.begin()->second.front() );
1546         }
1547       while ( !faceQueue.empty() )
1548         {
1549           const Cell* face = faceQueue.front();
1550           faceQueue.pop();
1551
1552           // loop on links of <face>
1553           for ( int i = 0; i < (int)face->_nodes.size(); ++i )
1554             {
1555               Link link = face->link( i );
1556               // find the neighbor faces
1557               lfIt = linkFacesMap.find( link );
1558               int nbFaceByLink = 0;
1559               list< const Cell* > ml;
1560               if ( lfIt != linkFacesMap.end() )
1561                 {
1562                   list<const Cell*> & fList = lfIt->second;
1563                   list<const Cell*>::iterator fIt = fList.begin();
1564                   assert( fIt != fList.end() );
1565                   for ( ; fIt != fList.end(); fIt++, nbFaceByLink++ )
1566                     {
1567                       ml.push_back( *fIt );
1568                       if ( *fIt != face ) // wrongly oriented neighbor face
1569                         {
1570                           const Cell* badFace = *fIt;
1571                           // reverse and remove badFace from linkFacesMap
1572                           for ( int j = 0; j < (int)badFace->_nodes.size(); ++j )
1573                             {
1574                               Link badlink = badFace->link( j );
1575                               if ( badlink == link ) continue;
1576                               lfIt2 = linkFacesMap.find( badlink );
1577                               if ( lfIt2 != linkFacesMap.end() )
1578                                 {
1579                                   list<const Cell*> & ff = lfIt2->second;
1580                                   ff.erase( find( ff.begin(), ff.end(), badFace ));
1581                                   if ( ff.empty() )
1582                                     linkFacesMap.erase( lfIt2 );
1583                                 }
1584                             }
1585                           badFace->_reverse = true; // reverse
1586                           //INFOS_MED( "REVERSE " << *badFace );
1587                           faceQueue.push( badFace );
1588                         }
1589                     }
1590                   linkFacesMap.erase( lfIt );
1591                 }
1592               // add good neighbors to the queue
1593               Link revLink( link.second, link.first );
1594               lfIt = linkFacesMap.find( revLink );
1595               if ( lfIt != linkFacesMap.end() )
1596                 {
1597                   list<const Cell*> & fList = lfIt->second;
1598                   list<const Cell*>::iterator fIt = fList.begin();
1599                   for ( ; fIt != fList.end(); fIt++, nbFaceByLink++ )
1600                     {
1601                       ml.push_back( *fIt );
1602                       if ( *fIt != face )
1603                         faceQueue.push( *fIt );
1604                     }
1605                   linkFacesMap.erase( lfIt );
1606                 }
1607               if ( nbFaceByLink > 2 )
1608                 {
1609                   if ( manifold )
1610                     {
1611                       list<const Cell*>::iterator ii = ml.begin();
1612                       cout << nbFaceByLink << " faces by 1 link:" << endl;
1613                       for( ; ii!= ml.end(); ii++ )
1614                         cout << "in sub-mesh <" << fgm[ *ii ]->_name << "> " << **ii << endl;
1615                     }
1616                   manifold = false;
1617                 }
1618             } // loop on links of the being checked face
1619         } // loop on the face queue
1620     } // while ( !linkFacesMap.empty() )
1621
1622   if ( !manifold )
1623     cout << " -> Non manifold mesh, faces orientation may be incorrect" << endl;
1624 }
1625
1626 //================================================================================
1627 /*!
1628  * \brief Orient volumes according to MED conventions:
1629  * normal of a bottom (first) face should be outside
1630  */
1631 //================================================================================
1632
1633 void IntermediateMED::orientVolumes()
1634 {
1635   set<Cell>::const_iterator elemIt, elemEnd;
1636   vector< pair<int,int> > swapVec;
1637
1638   CellsByDimIterator cellsIt( *this, 3 );
1639   while ( const set<Cell > * elems = cellsIt.nextType() )
1640     {
1641       TCellType cellType = cellsIt.type();
1642       elemIt = elems->begin(), elemEnd = elems->end();
1643       int nbBottomNodes = 0;
1644       switch ( cellType )
1645         {
1646         case NORM_TETRA4:
1647         case NORM_TETRA10:
1648         case NORM_PENTA6:
1649         case NORM_PENTA15:
1650           nbBottomNodes = 3; break;
1651         case NORM_PYRA5:
1652         case NORM_PYRA13:
1653         case NORM_HEXA8:
1654         case NORM_HEXA20:
1655           nbBottomNodes = 4; break;
1656         default: continue;
1657         }
1658       getReverseVector( cellType, swapVec );
1659
1660       for ( ; elemIt != elemEnd; elemIt++ )
1661         {
1662           // find a normal to the bottom face
1663           const double* n[4];
1664           n[0] = nodeCoords( elemIt->_nodes[0]); // 3 bottom nodes
1665           n[1] = nodeCoords( elemIt->_nodes[1]);
1666           n[2] = nodeCoords( elemIt->_nodes[2]);
1667           n[3] = nodeCoords( elemIt->_nodes[nbBottomNodes]); // a top node
1668           double vec01[3]; // vector n[0]-n[1]
1669           vec01[0] = n[1][0] - n[0][0];
1670           vec01[1] = n[1][1] - n[0][1];
1671           vec01[2] = n[1][2] - n[0][2];
1672           double vec02 [3]; // vector n[0]-n[2]
1673           vec02[0] = n[2][0] - n[0][0];
1674           vec02[1] = n[2][1] - n[0][1];
1675           vec02[2] = n[2][2] - n[0][2];
1676           double normal [3]; // vec01 ^ vec02
1677           normal[0] = vec01[1] * vec02[2] - vec01[2] * vec02[1];
1678           normal[1] = vec01[2] * vec02[0] - vec01[0] * vec02[2];
1679           normal[2] = vec01[0] * vec02[1] - vec01[1] * vec02[0];
1680           // check if the 102 angle is convex
1681           if ( nbBottomNodes > 3 )
1682             {
1683               const double* n3 = nodeCoords( elemIt->_nodes[nbBottomNodes-1] );// last bottom node
1684               double vec03 [3];  // vector n[0]-n3
1685               vec03[0] = n3[0] - n[0][0];
1686               vec03[1] = n3[1] - n[0][1];
1687               vec03[2] = n3[2] - n[0][2];
1688               if ( fabs( normal[0]+normal[1]+normal[2] ) <= numeric_limits<double>::max() ) // vec01 || vec02
1689                 {
1690                   normal[0] = vec01[1] * vec03[2] - vec01[2] * vec03[1]; // vec01 ^ vec03
1691                   normal[1] = vec01[2] * vec03[0] - vec01[0] * vec03[2];
1692                   normal[2] = vec01[0] * vec03[1] - vec01[1] * vec03[0];
1693                 }
1694               else
1695                 {
1696                   double vec [3]; // normal ^ vec01
1697                   vec[0] = normal[1] * vec01[2] - normal[2] * vec01[1];
1698                   vec[1] = normal[2] * vec01[0] - normal[0] * vec01[2];
1699                   vec[2] = normal[0] * vec01[1] - normal[1] * vec01[0];
1700                   double dot2 = vec[0]*vec03[0] + vec[1]*vec03[1] + vec[2]*vec03[2]; // vec*vec03
1701                   if ( dot2 < 0 ) // concave -> reverse normal
1702                     {
1703                       normal[0] *= -1;
1704                       normal[1] *= -1;
1705                       normal[2] *= -1;
1706                     }
1707                 }
1708             }
1709           // direction from top to bottom
1710           double tbDir[3];
1711           tbDir[0] = n[0][0] - n[3][0];
1712           tbDir[1] = n[0][1] - n[3][1];
1713           tbDir[2] = n[0][2] - n[3][2];
1714
1715           // compare 2 directions: normal and top-bottom
1716           double dot = normal[0]*tbDir[0] + normal[1]*tbDir[1] + normal[2]*tbDir[2];
1717           if ( dot < 0. ) // need reverse
1718             reverse( *elemIt, swapVec );
1719
1720         } // loop on volumes of one geometry
1721     } // loop on 3D geometry types
1722
1723 }
1724
1725 //================================================================================
1726 /*!
1727  * \brief Assign new IDs to nodes by skipping not used nodes and return their number
1728  */
1729 //================================================================================
1730
1731 int NodeContainer::numberNodes()
1732 {
1733   int id = 1;
1734   for ( size_t i = 0; i < _nodes.size(); ++i )
1735     for ( size_t j = 0; j < _nodes[i].size(); ++j )
1736       if ( _nodes[i][j].isUsed() )
1737         _nodes[i][j]._number = id++;
1738   return id-1;
1739 }
1740
1741
1742 //================================================================================
1743 /*!
1744  * \brief Assign new IDs to elements
1745  */
1746 //================================================================================
1747
1748 void IntermediateMED::numberElements()
1749 {
1750   set<Cell>::const_iterator elemIt, elemEnd;
1751
1752   // numbering _cells of type NORM_POINT1 by node number
1753   {
1754     const set<Cell>& points = _cellsByType[ INTERP_KERNEL::NORM_POINT1 ];
1755     elemIt = points.begin(), elemEnd = points.end();
1756     for ( ; elemIt != elemEnd; ++elemIt )
1757       elemIt->_number = elemIt->_nodes[0]->_number;
1758   }
1759
1760   // numbering 1D-3D _cells
1761   for ( int dim = 1; dim <= 3; ++dim )
1762     {
1763       // check if re-numeration is needed (to try to keep elem oreder as in sauve file )
1764       bool ok = true, renumEntity = false;
1765       CellsByDimIterator cellsIt( *this, dim );
1766       int prevNbElems = 0;
1767       while ( const set<Cell> * typeCells = cellsIt.nextType() )
1768         {
1769           TID minNumber = std::numeric_limits<TID>::max(), maxNumber = 0;
1770           for ( elemIt = typeCells->begin(), elemEnd = typeCells->end(); elemIt!=elemEnd; ++elemIt)
1771             {
1772               if ( elemIt->_number < minNumber ) minNumber = elemIt->_number;
1773               if ( elemIt->_number > maxNumber ) maxNumber = elemIt->_number;
1774             }
1775           TID typeSize = typeCells->size();
1776           if ( typeSize != maxNumber - minNumber + 1 )
1777             ok = false;
1778           if ( prevNbElems != 0 ) {
1779             if ( minNumber == 1 )
1780               renumEntity = true;
1781             else if ( prevNbElems+1 != (int)minNumber )
1782               ok = false;
1783           }
1784           prevNbElems += typeSize;
1785         }
1786
1787       if ( ok && renumEntity ) // each geom type was numerated separately
1788         {
1789           cellsIt.init( dim );
1790           prevNbElems = cellsIt.nextType()->size(); // no need to renumber the first type
1791           while ( const set<Cell> * typeCells = cellsIt.nextType() )
1792             {
1793               for ( elemIt = typeCells->begin(), elemEnd = typeCells->end(); elemIt!=elemEnd; ++elemIt)
1794                 elemIt->_number += prevNbElems;
1795               prevNbElems += typeCells->size();
1796             }
1797         }
1798       if ( !ok )
1799         {
1800           int cellID=1;
1801           cellsIt.init( dim );
1802           while ( const set<Cell> * typeCells = cellsIt.nextType() )
1803             for ( elemIt = typeCells->begin(), elemEnd = typeCells->end(); elemIt!=elemEnd; ++elemIt)
1804               elemIt->_number = cellID++;
1805         }
1806     }
1807 }
1808
1809 //================================================================================
1810 /*!
1811  * \brief Creates coord array
1812  */
1813 //================================================================================
1814
1815 ParaMEDMEM::DataArrayDouble * IntermediateMED::getCoords()
1816 {
1817   DataArrayDouble* coordArray = DataArrayDouble::New();
1818   coordArray->alloc( _nbNodes, _spaceDim );
1819   double * coordPrt = coordArray->getPointer();
1820   for ( int i = 0, nb = _points.size(); i < nb; ++i )
1821     {
1822       Node* n = getNode( i+1 );
1823       if ( n->isUsed() )
1824         {
1825           const double* nCoords = nodeCoords( n );
1826           std::copy( nCoords, nCoords+_spaceDim, coordPrt );
1827           coordPrt += _spaceDim;
1828         }
1829     }
1830   return coordArray;
1831 }
1832
1833 //================================================================================
1834 /*!
1835  * \brief Sets connectivity of elements to the mesh
1836  *  \param mesh - mesh to fill in
1837  *  \param coords - coordinates that must be shared by all meshes of different dim
1838  */
1839 //================================================================================
1840
1841 void IntermediateMED::setConnectivity( ParaMEDMEM::MEDFileUMesh*    mesh,
1842                                        ParaMEDMEM::DataArrayDouble* coords )
1843 {
1844   int meshDim = 0;
1845
1846   mesh->setCoords( coords );
1847
1848   set<Cell>::const_iterator elemIt, elemEnd;
1849   for ( int dim = 3; dim > 0; --dim )
1850   {
1851     CellsByDimIterator dimCells( *this, dim );
1852
1853     int nbOfCells = 0;
1854     while ( const std::set<Cell > * cells = dimCells.nextType() )
1855       nbOfCells += cells->size();
1856     if ( nbOfCells == 0 )
1857       continue;
1858
1859     if ( !meshDim ) meshDim = dim;
1860
1861     MEDCouplingUMesh* dimMesh = MEDCouplingUMesh::New();
1862     dimMesh->setCoords( coords );
1863     dimMesh->setMeshDimension( dim );
1864     dimMesh->allocateCells( nbOfCells );
1865
1866     int prevNbCells = 0;
1867     dimCells.init( dim );
1868     while ( const std::set<Cell > * cells = dimCells.nextType() )
1869       {
1870         // fill connectivity array to take into account order of elements in the sauv file
1871         const int nbCellNodes = cells->begin()->_nodes.size();
1872         vector< TID > connectivity( cells->size() * nbCellNodes );
1873         int * nodalConnOfCell;
1874         for ( elemIt = cells->begin(), elemEnd = cells->end(); elemIt != elemEnd; ++elemIt )
1875           {
1876             const Cell& cell = *elemIt;
1877             const int index = cell._number - 1 - prevNbCells;
1878             nodalConnOfCell = &connectivity[ index * nbCellNodes ];
1879             if ( cell._reverse )
1880               for ( int i = nbCellNodes-1; i >= 0; --i )
1881                 *nodalConnOfCell++ = cell._nodes[i]->_number - 1;
1882             else
1883               for ( int i = 0; i < nbCellNodes; ++i )
1884                 *nodalConnOfCell++ = cell._nodes[i]->_number - 1;
1885           }
1886         prevNbCells += cells->size();
1887
1888         // fill dimMesh
1889         TCellType cellType = dimCells.type();
1890         nodalConnOfCell = &connectivity[0];
1891         for ( size_t i = 0; i < cells->size(); ++i, nodalConnOfCell += nbCellNodes )
1892           dimMesh->insertNextCell( cellType, nbCellNodes, nodalConnOfCell );
1893       }
1894     dimMesh->finishInsertingCells();
1895     mesh->setMeshAtLevel( dim - meshDim, dimMesh );
1896     dimMesh->decrRef();
1897   }
1898 }
1899
1900 //================================================================================
1901 /*!
1902  * \brief Fill in the mesh with groups
1903  *  \param mesh - mesh to fill in
1904  */
1905 //================================================================================
1906
1907 void IntermediateMED::setGroups( ParaMEDMEM::MEDFileUMesh* mesh )
1908 {
1909   bool isMeshNameSet = false;
1910   const int meshDim = mesh->getMeshDimension();
1911   for ( int dim = 0; dim <= meshDim; ++dim )
1912     {
1913       const int meshDimRelToMaxExt = ( dim == 0 ? 1 : dim - meshDim );
1914
1915       vector<const DataArrayInt *> medGroups;
1916       vector<MEDCouplingAutoRefCountObjectPtr<DataArrayInt> > refGroups;
1917       for ( size_t i = 0; i < _groups.size(); ++i )
1918         {
1919           Group& grp = _groups[i];
1920           if ( (int)getDim( &grp ) != dim &&
1921                grp._groups.empty() ) // to allow groups on diff dims
1922             continue;
1923           // convert only named groups or field supports
1924           if ( grp.empty() || (grp._name.empty() && !grp._isProfile ))
1925             continue;
1926           //if ( grp._medGroup ) continue; // already converted
1927
1928           // sort cells by ID and remember their initial order in the group
1929           TCellToOrderMap cell2order;
1930           unsigned orderInGroup = 0;
1931           vector< Group* > groupVec;
1932           if ( grp._groups.empty() ) groupVec.push_back( & grp );
1933           else                       groupVec = grp._groups;
1934           for ( size_t iG = 0; iG < groupVec.size(); ++iG )
1935             {
1936               Group* aG = groupVec[ iG ];
1937               if ( (int)getDim( aG ) != dim )
1938                 continue;
1939               for ( size_t iC = 0; iC < aG->_cells.size(); ++iC )
1940                 cell2order.insert( cell2order.end(), make_pair( aG->_cells[iC], orderInGroup++ ));
1941             }
1942           if ( cell2order.empty() )
1943             continue;
1944           bool isSelfIntersect = ( orderInGroup != cell2order.size() );
1945           if ( isSelfIntersect ) // self intersecting group
1946             {
1947               ostringstream msg;
1948               msg << "Self intersecting sub-mesh: id = " << i+1
1949                   << ", name = |" << grp._name << "|" << endl
1950                   << " nb unique elements = " << cell2order.size() << endl
1951                   << " total nb elements  = " << orderInGroup;
1952               if ( grp._isProfile )
1953                 {
1954                   THROW_IK_EXCEPTION( msg.str() );
1955                 }
1956               else
1957                 {
1958                   cout << msg.str() << endl;
1959                 }
1960             }
1961           // create a med group
1962           grp._medGroup = DataArrayInt::New();
1963           grp._medGroup->setName( grp._name.c_str() );
1964           grp._medGroup->alloc( cell2order.size(), /*nbOfCompo=*/1 );
1965           int * idsPrt = grp._medGroup->getPointer();
1966           TCellToOrderMap::iterator cell2orderIt, cell2orderEnd = cell2order.end();
1967           for ( cell2orderIt = cell2order.begin(); cell2orderIt != cell2orderEnd; ++cell2orderIt )
1968             *idsPrt++ = (*cell2orderIt).first->_number - 1;
1969
1970           // try to set the mesh name
1971           if ( !isMeshNameSet &&
1972                dim == meshDim &&
1973                !grp._name.empty() &&
1974                grp.size() == mesh->getSizeAtLevel( meshDimRelToMaxExt ))
1975             {
1976               mesh->setName( grp._name.c_str() );
1977               isMeshNameSet = true;
1978             }
1979           if ( !grp._name.empty() )
1980             {
1981               medGroups.push_back( grp._medGroup );
1982             }
1983           // set relocation table
1984           setRelocationTable( &grp, cell2order );
1985
1986           // Issue 0021311. Use case: a gibi group has references (recorded in pile 1)
1987           // and several names (pile 27) refer (pile 10) to this group.
1988           // We create a copy of this group per each named reference
1989           for ( unsigned iRef = 0 ; iRef < grp._refNames.size(); ++iRef )
1990             if ( !grp._refNames[ iRef ].empty() )
1991               {
1992                 refGroups.push_back( grp._medGroup->deepCpy() );
1993                 refGroups.back()->setName( grp._refNames[ iRef ].c_str() );
1994                 medGroups.push_back( refGroups.back() );
1995               }
1996         }
1997       mesh->setGroupsAtLevel( meshDimRelToMaxExt, medGroups );
1998     }
1999 }
2000
2001 //================================================================================
2002 /*!
2003  * \brief Return true if the group is on all elements and return its relative dimension
2004  */
2005 //================================================================================
2006
2007 bool IntermediateMED::isOnAll( const Group* grp, int & dimRel ) const
2008 {
2009   int dim = getDim( grp );
2010
2011   int nbElems = 0;
2012   CellsByDimIterator dimCells( *this, dim );
2013   while ( const set<Cell > * cells = dimCells.nextType() )
2014     nbElems += cells->size();
2015
2016   const bool onAll = ( nbElems == grp->size() );
2017
2018   if ( dim == 0 )
2019     dimRel = 0;
2020   else
2021     {
2022       int meshDim = 3;
2023       for ( ; meshDim > 0; --meshDim )
2024         {
2025           dimCells.init( meshDim );
2026           if ( dimCells.nextType() )
2027             break;
2028         }
2029       dimRel = dim - meshDim;
2030     }
2031   return onAll;
2032 }
2033
2034 //================================================================================
2035 /*!
2036  * \brief Makes fields from own data
2037  */
2038 //================================================================================
2039
2040 ParaMEDMEM::MEDFileFields * IntermediateMED::makeMEDFileFields(ParaMEDMEM::MEDFileUMesh* mesh)
2041 {
2042   if ( _nodeFields.empty() && _cellFields.empty() ) return 0;
2043
2044   // set long names
2045   set< string > usedFieldNames;
2046   setFieldLongNames(usedFieldNames);
2047
2048   MEDFileFields* fields = MEDFileFields::New();
2049
2050   for ( size_t i = 0; i < _nodeFields.size(); ++i )
2051     setFields( _nodeFields[i], fields, mesh, i+1, usedFieldNames );
2052
2053   for ( size_t i = 0; i < _cellFields.size(); ++i )
2054     setFields( _cellFields[i], fields, mesh, i+1, usedFieldNames );
2055
2056   return fields;
2057 }
2058
2059 //================================================================================
2060 /*!
2061  * \brief Make med fields from a SauvUtilities::DoubleField
2062  */
2063 //================================================================================
2064
2065 void IntermediateMED::setFields( SauvUtilities::DoubleField* fld,
2066                                  ParaMEDMEM::MEDFileFields*  medFields,
2067                                  ParaMEDMEM::MEDFileUMesh*   mesh,
2068                                  const TID                   castemID,
2069                                  set< string >&              usedFieldNames)
2070 {
2071   if ( !fld || !fld->isMedCompatible() ) return;
2072
2073   // if ( !fld->hasCommonSupport() ):
2074   //     each sub makes MEDFileFieldMultiTS
2075   // else:
2076   //     unite several subs into a MEDCouplingFieldDouble
2077
2078   const bool uniteSubs = fld->hasCommonSupport();
2079   if ( !uniteSubs )
2080     cout << "Castem field #" << castemID << " " << fld->_name
2081          << " is incompatible with MED format, so we split it into several fields" << endl;
2082
2083   for ( size_t iSub = 0; iSub < fld->_sub.size(); )
2084     {
2085       // set field name
2086       if ( !uniteSubs || fld->_name.empty() )
2087         makeFieldNewName( usedFieldNames, fld );
2088
2089       // allocate values
2090       DataArrayDouble * values = DataArrayDouble::New();
2091       values->alloc( fld->getNbTuples(iSub), fld->_sub[iSub].nbComponents() );
2092
2093       // set values
2094       double * valPtr = values->getPointer();
2095       if ( uniteSubs )
2096         {
2097           int nbElems = fld->_group->size();
2098           for ( int elemShift = 0; elemShift < nbElems; )
2099             elemShift += fld->setValues( valPtr, iSub++, elemShift );
2100           setTS( fld, values, medFields, mesh );
2101         }
2102       else
2103         {
2104           fld->setValues( valPtr, iSub++ );
2105           setTS( fld, values, medFields, mesh, iSub );
2106         }
2107     }
2108 }
2109
2110 //================================================================================
2111 /*!
2112  * \brief Store value array of a field into med fields
2113  */
2114 //================================================================================
2115
2116 void IntermediateMED::setTS( SauvUtilities::DoubleField*  fld,
2117                              ParaMEDMEM::DataArrayDouble* values,
2118                              ParaMEDMEM::MEDFileFields*   medFields,
2119                              ParaMEDMEM::MEDFileUMesh*    mesh,
2120                              const int                    iSub)
2121 {
2122   // analyze a field support
2123   const Group* support = fld->getSupport();
2124   int dimRel;
2125   const bool onAll = isOnAll( support, dimRel );
2126   if ( !onAll && support->_name.empty() )
2127     {
2128       const_cast<Group*>(support)->_name += "PFL_" + fld->_name;
2129       support->_medGroup->setName( support->_name.c_str() );
2130     }
2131
2132   // make a time-stamp
2133   MEDCouplingFieldDouble * timeStamp = MEDCouplingFieldDouble::New( fld->getMedType(),
2134                                                                     fld->getMedTimeDisc() );
2135   timeStamp->setName( fld->_name.c_str() );
2136   timeStamp->setDescription( fld->_description.c_str() );
2137   MEDCouplingAutoRefCountObjectPtr< MEDCouplingUMesh > dimMesh = mesh->getMeshAtLevel( dimRel );
2138   timeStamp->setMesh( dimMesh );
2139   for ( size_t i = 0; i < (size_t)fld->_sub[iSub].nbComponents(); ++i )
2140     values->setInfoOnComponent( i, fld->_sub[iSub]._comp_names[ i ].c_str() );
2141   timeStamp->setArray( values );
2142   values->decrRef();
2143
2144   // get a field to add the time-stamp
2145   bool isNewMedField = false;
2146   if ( !fld->_curMedField || fld->_name != fld->_curMedField->getName() )
2147     {
2148       fld->_curMedField = MEDFileFieldMultiTS::New();
2149       isNewMedField = true;
2150     }
2151
2152   // set an order
2153   timeStamp->setOrder( fld->_curMedField->getNumberOfTS() );
2154
2155   // add the time-stamp
2156   if ( onAll )
2157     fld->_curMedField->appendFieldNoProfileSBT( timeStamp );
2158   else
2159     fld->_curMedField->appendFieldProfile( timeStamp, mesh, dimRel, support->_medGroup );
2160   timeStamp->decrRef();
2161
2162   if ( isNewMedField ) // timeStamp must be added before this
2163     {
2164       medFields->pushField( fld->_curMedField );
2165     }
2166 }
2167
2168 //================================================================================
2169 /*!
2170  * \brief Make a new unique name for a field
2171  */
2172 //================================================================================
2173
2174 void IntermediateMED::makeFieldNewName(std::set< std::string >&    usedNames,
2175                                        SauvUtilities::DoubleField* fld )
2176 {
2177   string base = fld->_name;
2178   if ( base.empty() )
2179     {
2180       base = "F_";
2181     }
2182   else
2183     {
2184       string::size_type pos = base.rfind('_');
2185       if ( pos != string::npos )
2186         base = base.substr( 0, pos+1 );
2187       else
2188         base += '_';
2189     }
2190
2191   int i = 1;
2192   do
2193     {
2194       fld->_name = base + SauvUtilities::toString( i++ );
2195     }
2196   while( !usedNames.insert( fld->_name ).second );
2197 }
2198
2199 //================================================================================
2200 /*!
2201  * \brief Return a vector ready to fill in
2202  */
2203 //================================================================================
2204
2205 std::vector< double >& DoubleField::addComponent( int nb_values )
2206 {
2207   _comp_values.push_back( std::vector< double >() );
2208   std::vector< double >& res = _comp_values.back();
2209   res.resize( nb_values );
2210   return res;
2211 }
2212
2213 DoubleField::~DoubleField()
2214 {
2215   if(_curMedField)
2216     _curMedField->decrRef();
2217 }
2218
2219 //================================================================================
2220 /*!
2221  * \brief Returns a supporting group
2222  */
2223 //================================================================================
2224
2225 const Group* DoubleField::getSupport( const int iSub ) const
2226 {
2227   return _group ? _group : _sub[iSub]._support;
2228 }
2229
2230 //================================================================================
2231 /*!
2232  * \brief Return true if each sub-component is a time stamp
2233  */
2234 //================================================================================
2235
2236 bool DoubleField::isMultiTimeStamps() const
2237 {
2238   if ( _sub.size() < 2 )
2239     return false;
2240   bool sameSupports = true;
2241   Group* grp1 = _sub[0]._support;
2242   for ( size_t i = 1; i < _sub.size() && sameSupports; ++i )
2243     sameSupports = ( grp1 == _sub[i]._support );
2244
2245   return sameSupports;
2246 }
2247
2248 //================================================================================
2249 /*!
2250  * \brief True if the field can be converted into the med field
2251  */
2252 //================================================================================
2253
2254 bool DoubleField::isMedCompatible() const
2255 {
2256   for ( size_t iSub = 0; iSub < _sub.size(); ++iSub )
2257     {
2258       if ( !getSupport(iSub) || !getSupport(iSub)->_medGroup )
2259         THROW_IK_EXCEPTION("SauvReader INTERNAL ERROR: NULL field support");
2260
2261       if ( !_sub[iSub].isValidNbGauss() )
2262         {
2263           cout << "Skip field <" << _name << "> : different nb of gauss points in components" <<endl;
2264           return false;
2265         }
2266     }
2267   // check if there are no gauss or nbGauss() == nbCellNodes,
2268   // else we lack info on gauss point localization
2269   // TODO?
2270   return true;
2271 }
2272
2273 //================================================================================
2274 /*!
2275  * \brief return true if all sub-components has same components and same nbGauss
2276  */
2277 //================================================================================
2278
2279 bool DoubleField::hasSameComponentsBySupport() const
2280 {
2281   vector< _Sub_data >::const_iterator sub_data = _sub.begin();
2282   const _Sub_data& first_sub_data = *sub_data;
2283   for ( ++sub_data ; sub_data != _sub.end(); ++sub_data )
2284   {
2285     if ( first_sub_data._comp_names != sub_data->_comp_names )
2286       return false; // diff names of components
2287
2288     if ( first_sub_data._nb_gauss != sub_data->_nb_gauss &&
2289          first_sub_data._support->_cellType == sub_data->_support->_cellType)
2290       return false; // diff nb of gauss points on same cell type
2291   }
2292   return true;
2293 }
2294
2295 //================================================================================
2296 /*!
2297  * \brief Return type of MEDCouplingFieldDouble
2298  */
2299 //================================================================================
2300
2301 ParaMEDMEM::TypeOfField DoubleField::getMedType( const int iSub ) const
2302 {
2303   using namespace INTERP_KERNEL;
2304
2305   const Group* grp = hasCommonSupport() ? _group : _sub[iSub]._support;
2306   if ( _sub[iSub].nbGauss() > 1 )
2307     {
2308       const CellModel& cm = CellModel::GetCellModel( _sub[iSub]._support->_cellType );
2309       return (int) cm.getNumberOfNodes() == _sub[iSub].nbGauss() ? ON_GAUSS_NE : ON_GAUSS_PT;
2310     }
2311   else
2312     {
2313       return getDim( grp ) == 0 ? ON_NODES : ON_CELLS;
2314     }
2315 }
2316
2317 //================================================================================
2318 /*!
2319  * \brief Return TypeOfTimeDiscretization
2320  */
2321 //================================================================================
2322
2323 ParaMEDMEM::TypeOfTimeDiscretization DoubleField::getMedTimeDisc() const
2324 {
2325   return ONE_TIME;
2326   // NO_TIME = 4,
2327   // ONE_TIME = 5,
2328   // LINEAR_TIME = 6,
2329   // CONST_ON_TIME_INTERVAL = 7
2330 }
2331
2332 //================================================================================
2333 /*!
2334  * \brief Return nb tuples to be used to allocate DataArrayDouble
2335  */
2336 //================================================================================
2337
2338 int DoubleField::getNbTuples( const int iSub ) const
2339 {
2340   int nb = 0;
2341   if ( hasCommonSupport() && !_group->_groups.empty() )
2342     for ( size_t i = 0; i < _group->_groups.size(); ++i )
2343       nb += _sub[i].nbGauss() * _sub[i]._support->size();
2344   else
2345     nb = _sub[iSub].nbGauss() * getSupport(iSub)->size();
2346   return nb;
2347 }
2348
2349 //================================================================================
2350 /*!
2351  * \brief Store values of a sub-component and return nb of elements in the iSub
2352  */
2353 //================================================================================
2354
2355 int DoubleField::setValues( double * valPtr, const int iSub, const int elemShift ) const
2356 {
2357   // find values for iSub
2358   int iComp = 0;
2359   for ( int iS = 0; iS < iSub; ++iS )
2360     iComp += _sub[iS].nbComponents();
2361   const vector< double > * compValues = &_comp_values[ iComp ];
2362
2363   const vector< unsigned >& relocTable = getSupport( iSub )->_relocTable;
2364
2365   // Set values
2366
2367   const int nbElems      = _sub[iSub]._support->size();
2368   const int nbGauss      = _sub[iSub].nbGauss();
2369   const int nbComponents = _sub[iSub].nbComponents();
2370   const int nbValsByElem = nbComponents * nbGauss;
2371   // check nb values
2372   int nbVals = 0;
2373   for ( iComp = 0; iComp < nbComponents; ++iComp )
2374     nbVals += compValues[iComp].size();
2375   if ( nbVals != nbElems * nbValsByElem )
2376     THROW_IK_EXCEPTION("SauvMedConvertor.cxx: support size mismatches field size");
2377   // compute nb values in previous subs
2378   int valsShift = 0;
2379   for ( int iS = iSub-1, shift = elemShift; shift > 0; )
2380   {
2381     int nbE = _sub[iS]._support->size();
2382     shift -= nbE;
2383     valsShift += nbE * _sub[iS].nbComponents() * _sub[iS].nbGauss();
2384   }
2385   for ( int iE = 0; iE < nbElems; ++iE )
2386     {
2387       int iMed = valsShift + nbValsByElem * ( relocTable.empty() ? iE : relocTable[iE+elemShift]-elemShift );
2388       for ( iComp = 0; iComp < nbComponents; ++iComp )
2389         for ( int iG = 0; iG < nbGauss; ++iG )
2390           valPtr[ iMed + iG * nbComponents + iComp ] = compValues[iComp][ iE * nbGauss + iG ];
2391     }
2392   return nbElems;
2393 }
2394
2395 //================================================================================
2396 /*!
2397  * \brief Destructor of IntermediateMED
2398  */
2399 //================================================================================
2400
2401 IntermediateMED::~IntermediateMED()
2402 {
2403   for ( size_t i = 0; i < _nodeFields.size(); ++i )
2404     if ( _nodeFields[i] )
2405       delete _nodeFields[i];
2406   _nodeFields.clear();
2407
2408   for ( size_t i = 0; i < _cellFields.size(); ++i )
2409     if ( _cellFields[i] )
2410       delete _cellFields[i];
2411   _cellFields.clear();
2412
2413   for ( size_t i = 0; i < _groups.size(); ++i )
2414     if ( _groups[i]._medGroup )
2415       _groups[i]._medGroup->decrRef();
2416 }
2417
2418 //================================================================================
2419 /*!
2420  * \brief CellsByDimIterator constructor
2421  */
2422 CellsByDimIterator::CellsByDimIterator( const IntermediateMED & medi, int dimm)
2423 {
2424   myImed = & medi;
2425   init( dimm );
2426 }
2427 /*!
2428  * \brief Initialize iteration on cells of given dimention
2429  */
2430 void CellsByDimIterator::init(const int  dimm)
2431 {
2432   myCurType = -1;
2433   myTypeEnd = INTERP_KERNEL::NORM_HEXA20 + 1;
2434   myDim = dimm;
2435 }
2436 /*!
2437  * \brief return next set of Cell's of required dimension
2438  */
2439 const std::set< Cell > * CellsByDimIterator::nextType()
2440 {
2441   while ( ++myCurType < myTypeEnd )
2442     if ( !myImed->_cellsByType[myCurType].empty() && ( myDim < 0 || dim(false) == myDim ))
2443       return & myImed->_cellsByType[myCurType];
2444   return 0;
2445 }
2446 /*!
2447  * \brief return dimension of cells returned by the last or further next()
2448  */
2449 int CellsByDimIterator::dim(const bool last) const
2450 {
2451   int typp = myCurType;
2452   if ( !last )
2453     while ( typp < myTypeEnd && myImed->_cellsByType[typp].empty() )
2454       ++typp;
2455   return typp < myTypeEnd ? getDimension( TCellType( typp )) : 4;
2456 }
2457 // END CellsByDimIterator ========================================================
2458