src/MEDLoader/SauvReader.cxx

   1 // Copyright (C) 2007-2019  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      : SauvReader.cxx
  20 // Created   : Tue Aug 16 13:57:42 2011
  21 // Author    : Edward AGAPOV (eap)
  22 //
  23
  24 #include "SauvReader.hxx"
  25
  26 #include "SauvMedConvertor.hxx"
  27 #include "MCAuto.hxx"
  28 #include "NormalizedUnstructuredMesh.hxx"
  29 #include "MEDCouplingRefCountObject.hxx"
  30
  31 #include <cstring>
  32 #include <sstream>
  33 #include <iostream>
  34
  35 using namespace MEDCoupling;
  36 using namespace SauvUtilities;
  37 using namespace std;
  38
  39 #define GIBI_EQUAL(var_str, stat_str) (strncmp (var_str, stat_str, strlen(stat_str)) == 0)
  40
  41 //================================================================================
  42 /*!
  43  * \brief Creates a reader of a given sauve file
  44  */
  45 //================================================================================
  46
  47 SauvReader* SauvReader::New(const std::string& fileName)
  48 {
  49   if ( fileName.empty() ) THROW_IK_EXCEPTION("Invalid file name");
  50
  51   MEDCoupling::MCAuto< SauvUtilities::FileReader> parser;
  52
  53   // try to open as XRD
  54   parser = new XDRReader( fileName.c_str() );
  55   if ( parser->open() )
  56     {
  57       SauvReader* reader = new SauvReader;
  58       reader->_fileReader = parser.retn();
  59       return reader;
  60     }
  61
  62   // try to open as ASCII
  63   parser = new ASCIIReader( fileName.c_str() );
  64   if ( parser->open() )
  65     {
  66       SauvReader* reader = new SauvReader;
  67       reader->_fileReader = parser.retn();
  68       return reader;
  69     }
  70
  71   THROW_IK_EXCEPTION("Unable to open file |"<< fileName << "|");
  72 }
  73 //================================================================================
  74 /*!
  75  * \brief Destructor
  76  */
  77 //================================================================================
  78
  79 SauvReader::~SauvReader()
  80 {
  81   _fileReader->decrRef();
  82 }
  83
  84 std::size_t SauvReader::getHeapMemorySizeWithoutChildren() const
  85 {
  86   return 0;
  87 }
  88
  89 std::vector<const BigMemoryObject *> SauvReader::getDirectChildrenWithNull() const
  90 {
  91   return std::vector<const BigMemoryObject *>();
  92 }
  93
  94 //================================================================================
  95 /*!
  96  * \brief Return current line of ASCII file to report an error
  97  */
  98 //================================================================================
  99
 100 std::string SauvReader::lineNb() const
 101 {
 102   if ( isASCII() )
 103     return string(" (line #") + SauvUtilities::toString
 104       ( static_cast<SauvUtilities::ASCIIReader*>( _fileReader )->lineNb() ) + ")";
 105
 106   return "";
 107 }
 108
 109 //================================================================================
 110 /*!
 111  * \brief Reads contents of the sauve file and convert it to MEDFileData
 112  */
 113 //================================================================================
 114
 115 MEDCoupling::MEDFileData * SauvReader::loadInMEDFileDS()
 116 {
 117   SauvUtilities::IntermediateMED iMed; // intermadiate DS
 118   _iMed = &iMed;
 119
 120   char* line; // a current line
 121   const char* enregistrement_type=" ENREGISTREMENT DE TYPE";
 122
 123   while ( getNextLine(line, /*raiseOEF=*/false)) // external loop looking for "ENREGISTREMENT DE TYPE"
 124     {
 125       if ( isASCII() && !GIBI_EQUAL( line, enregistrement_type ))
 126         continue; // "ENREGISTREMENT DE TYPE" not found -> read the next line
 127
 128       // read the number of a record
 129       int recordNumber;
 130       if ( isASCII() )
 131         recordNumber = atoi( line + strlen(enregistrement_type) + 1 );
 132       else
 133         recordNumber = getInt();
 134
 135       // read the record
 136       if ( recordNumber == 2 )
 137         readRecord2();
 138       else if (recordNumber == 4 )
 139         readRecord4();
 140       else if (recordNumber == 7 )
 141         readRecord7();
 142       else if (recordNumber == 8 )
 143         readRecord8();
 144       else if (recordNumber == 5 )
 145         break; // stop reading
 146       else
 147         if ( !isASCII() )
 148           THROW_IK_EXCEPTION("XDR : ENREGISTREMENT DE TYPE " << recordNumber << " not implemented!!!");
 149     }
 150
 151   MEDCoupling::MEDFileData* medFileData = iMed.convertInMEDFileDS();
 152
 153   return medFileData;
 154 }
 155
 156 //================================================================================
 157 /*!
 158  * \brief Reads "ENREGISTREMENT DE TYPE 4"
 159  */
 160 //================================================================================
 161
 162 void SauvReader::readRecord4()
 163 {
 164   if ( !isASCII() )
 165     {
 166       getInt(); // skip NIVEAU
 167       getInt(); // skip ERREUR
 168       _iMed->_spaceDim = getInt();
 169       getFloat(); // skip DENSITE
 170     }
 171   else
 172     {
 173       char* line;
 174       getNextLine(line);
 175       const char* s = " NIVEAU  15 NIVEAU ERREUR   0 DIMENSION";
 176       _iMed->_spaceDim = atoi( line + strlen( s ) + 1 );
 177       if ( !GIBI_EQUAL( line, " NIVEAU" ))
 178         THROW_IK_EXCEPTION( "Could not read space dimension" << lineNb() );
 179       }
 180   if ( _iMed->_spaceDim < 1 )
 181     THROW_IK_EXCEPTION( "Invalid space dimension:" << _iMed->_spaceDim );
 182 }
 183
 184 //================================================================================
 185 /*!
 186  * \brief Reads "ENREGISTREMENT DE TYPE 7"
 187  */
 188 //================================================================================
 189
 190 void SauvReader::readRecord7()
 191 {
 192   if ( !isASCII() )
 193     {
 194       getInt(); // skip NOMBRE INFO CASTEM2000
 195       getInt(); // skip IFOUR
 196       getInt(); // skip NIFOUR
 197       getInt(); // skip IFOMOD
 198       getInt(); // skip IECHO
 199       getInt(); // skip IIMPI
 200       getInt(); // skip IOSPI
 201       getInt(); // skip ISOTYP
 202       getInt(); // skip NSDPGE
 203     }
 204   else
 205     {
 206       // skip 3 lines:
 207       // NOMBRE INFO CASTEM2000   8
 208       // IFOUR   2 NIFOUR   0 IFOMOD   2 IECHO   1 IIMPI   0 IOSPI   0 ISOTYP   1
 209       // NSDPGE     0
 210       char* line;
 211       getNextLine(line);
 212       getNextLine(line);
 213       getNextLine(line);
 214     }
 215 }
 216
 217 //================================================================================
 218 /*!
 219  * \brief Reads "ENREGISTREMENT DE TYPE 8"
 220  */
 221 //================================================================================
 222
 223 void SauvReader::readRecord8()
 224 {
 225   // This record is useless (a constant table)
 226   // => we skip it
 227   int info;
 228   int nbIntToSkip;
 229   if ( !isASCII() )
 230     {
 231       getInt();
 232       info = getInt();
 233       int i = 0;
 234       if (info == 3) {
 235         // castem >= 18
 236         // 1 more line
 237         nbIntToSkip = 145;
 238       }
 239       else
 240         nbIntToSkip = 141;
 241
 242       while (i <= nbIntToSkip) {
 243         getInt();
 244         i ++;
 245       }
 246     }
 247 }
 248
 249 //================================================================================
 250 /*!
 251  * \brief Reads the pile number, nb of objects and nb named of objects
 252  */
 253 //================================================================================
 254
 255 int SauvReader::readPileNumber(int& nbNamedObjects, int& nbObjects)
 256 {
 257   // FORMAT(' PILE NUMERO',I4,'NBRE OBJETS NOMMES',I8,'NBRE OBJETS',I8)
 258   int pileNumber;
 259   if ( !isASCII() )
 260     {
 261       initIntReading(3);
 262       pileNumber     = getInt(); next();
 263       nbNamedObjects = getInt(); next();
 264       nbObjects      = getInt(); next();
 265     }
 266   else
 267     {
 268       char* line;
 269       getNextLine(line);
 270       const char *s1 = " PILE NUMERO", *s2 = "NBRE OBJETS NOMMES", *s3 = "NBRE OBJETS";
 271       if ( ! GIBI_EQUAL( line, s1 ) )
 272         THROW_IK_EXCEPTION("Could not read the pile number " << lineNb() );
 273       line           = line + strlen(s1);
 274       pileNumber     = atoi( line );
 275       line           = line + 4 + strlen(s2);
 276       nbNamedObjects = atoi( line );
 277       line           = line + 8 + strlen(s3);
 278       nbObjects      = atoi( line );
 279     }
 280   if ( nbNamedObjects<0 )
 281     THROW_IK_EXCEPTION("Invalid nb of named objects: " << nbNamedObjects  << lineNb() );
 282   if ( nbObjects<0)
 283     THROW_IK_EXCEPTION("Invalid nb of objects: " << nbObjects  << lineNb() );
 284   // It appears to be a valid case
 285   // if ( nbObjects<nbNamedObjects)
 286   //   THROW_IK_EXCEPTION("In PILE " << pileNumber <<
 287   //                      " nb of objects is less than nb of named objects"  << lineNb() );
 288   return pileNumber;
 289 }
 290
 291 //================================================================================
 292 /*!
 293  * \brief Reads "ENREGISTREMENT DE TYPE 2"
 294  */
 295 //================================================================================
 296
 297 void SauvReader::readRecord2()
 298 {
 299   if ( _iMed->_spaceDim == 0 )
 300     THROW_IK_EXCEPTION("Missing ENREGISTREMENT DE TYPE   4");
 301
 302   // read a pile number
 303   int pileNumber, nbNamedObjects, nbObjects;
 304   pileNumber = readPileNumber(nbNamedObjects, nbObjects);
 305
 306   if ( !_encounteredPiles.insert( pileNumber ).second && // piles may repeat
 307        isASCII())
 308     return;
 309
 310   // read object names and their indices
 311   vector<string> objectNames(nbNamedObjects);
 312   for ( initNameReading( nbNamedObjects ); more(); next() )
 313     objectNames[ index() ] = getName();
 314
 315   vector<int> nameIndices(nbNamedObjects);
 316   for ( initIntReading( nbNamedObjects ); more(); next() )
 317     nameIndices[ index() ] = getInt();
 318
 319   switch ( pileNumber )
 320     {
 321     case PILE_SOUS_MAILLAGE:
 322       read_PILE_SOUS_MAILLAGE(nbObjects, objectNames, nameIndices);
 323       break;
 324     case PILE_NODES_FIELD:
 325       read_PILE_NODES_FIELD(nbObjects, objectNames, nameIndices);
 326       break;
 327     case PILE_TABLES:
 328       read_PILE_TABLES(nbObjects, objectNames, nameIndices);
 329       break;
 330     case PILE_LREEL:
 331       read_PILE_LREEL(nbObjects, objectNames, nameIndices);
 332       break;
 333     case PILE_LOGIQUES:
 334       read_PILE_LOGIQUES(nbObjects, objectNames, nameIndices);
 335       break;
 336     case PILE_FLOATS:
 337       read_PILE_FLOATS(nbObjects, objectNames, nameIndices);
 338       break;
 339     case PILE_INTEGERS:
 340       read_PILE_INTEGERS(nbObjects, objectNames, nameIndices);
 341       break;
 342     case PILE_STRINGS:
 343       read_PILE_STRINGS(nbObjects, objectNames, nameIndices);
 344       break;
 345     case PILE_LMOTS:
 346       read_PILE_LMOTS(nbObjects, objectNames, nameIndices);
 347       break;
 348     case PILE_NOEUDS:
 349       read_PILE_NOEUDS(nbObjects, objectNames, nameIndices);
 350       break;
 351     case PILE_COORDONNEES:
 352       read_PILE_COORDONNEES(nbObjects, objectNames, nameIndices);
 353       break;
 354     case PILE_MODL:
 355       read_PILE_MODL(nbObjects, objectNames, nameIndices);
 356       break;
 357     case PILE_FIELD:
 358       read_PILE_FIELD(nbObjects, objectNames, nameIndices);
 359       break;
 360     default:
 361       if ( !isASCII() )
 362         THROW_IK_EXCEPTION("XDR : reading PILE " << pileNumber << " not implemented !!!");
 363     }
 364 }
 365
 366 //================================================================================
 367 /*!
 368  * \brief Reads "PILE NUMERO   1": gibi sub-meshes that are converted into med groups
 369  */
 370 //================================================================================
 371
 372 void SauvReader::read_PILE_SOUS_MAILLAGE(const int                 nbObjects,
 373                                          std::vector<std::string>& objectNames,
 374                                          std::vector<int>&         nameIndices)
 375 {
 376   _iMed->_groups.reserve(nbObjects*2); // fields may add some groups
 377
 378   char* line;
 379   map<int,int> strangeGroupType;
 380   int i;
 381
 382   for (int object=0; object!=nbObjects; ++object) // loop on sub-groups
 383     {
 384       initIntReading( 5 );
 385       int castemCellType = getIntNext();
 386       int nbSubGroups    = getIntNext();
 387       int nbReferences   = getIntNext();
 388       int nbNodesPerElem = getIntNext();
 389       int nbElements     = getIntNext();
 390
 391       _iMed->_groups.push_back(Group());
 392       SauvUtilities::Group & group = _iMed->_groups.back();
 393
 394       // Issue 0021311. Allocate places for names of referring groups
 395       // that will be possibly filled after reading long names from
 396       // PILE_TABLES and PILE_STRINGS
 397       group._refNames.resize( nbReferences );
 398
 399       // castemCellType=0 corresponds to a sub-mesh composed of other sub-meshes
 400       if (castemCellType==0 && nbSubGroups>0)
 401         {
 402           group._groups.resize( nbSubGroups );
 403           for ( initIntReading( nbSubGroups ); more(); next() )
 404             group._groups[ index() ] = & _iMed->_groups[ getInt() - 1 ];
 405           //std::sort( group._groups.begin(), group._groups.end() ); // for _groups comparison in getFieldSupport()
 406         }
 407       // skip references
 408       if ( isASCII() )
 409         for ( i = 0; i < nbReferences; i += 10 ) // FORMAT(10I8)
 410           getNextLine(line);
 411       else
 412         for (initIntReading(nbReferences); more(); next());
 413
 414       // skip colors
 415       if ( isASCII() )
 416         for ( i = 0; i < nbElements; i += 10 )
 417           getNextLine(line);
 418       else
 419         for (initIntReading(nbElements); more(); next());
 420
 421       // not a composite group
 422       if (castemCellType>0 && nbSubGroups==0)
 423         {
 424           group._cellType = SauvUtilities::gibi2medGeom(castemCellType);
 425
 426           initIntReading( nbElements * nbNodesPerElem );
 427           if ( group._cellType == INTERP_KERNEL::NORM_ERROR ) // look for group end
 428             {
 429               for ( ; more();  next());
 430               strangeGroupType.insert( make_pair( object, castemCellType ));
 431             }
 432           else
 433             {
 434               // if ( group._cellType == MED_POINT1 ) group._cellType = NORM_ERROR; // issue 21199
 435
 436               // read connectivity of elements of a group
 437               SauvUtilities::Cell ma( nbNodesPerElem );
 438               SauvUtilities::Node* pNode;
 439               group._cells.resize( nbElements );
 440               for ( i = 0; i < nbElements; ++i )
 441                 {
 442                   ma.init();
 443                   for ( int n = 0; n < nbNodesPerElem; ++n )
 444                     {
 445                       int nodeID = getIntNext();
 446                       pNode = _iMed->getNode( nodeID );
 447                       ma._nodes[n] = pNode;
 448                       _iMed->_nbNodes += ( !pNode->isUsed() );
 449                       pNode->_number = nodeID;
 450                     }
 451                   ma._number = _iMed->getNbCellsOfType( group._cellType ) + 1;
 452                   group._cells[i] = _iMed->insert( group._cellType, ma );
 453                 }
 454             }
 455         }
 456     } // loop on groups
 457
 458   // set group names
 459   for (i=0; i!=(int)objectNames.size(); ++i)
 460     {
 461       int grpID = nameIndices[i];
 462       SauvUtilities::Group & grp = _iMed->_groups[ grpID-1 ];
 463       if ( !grp._name.empty() ) // a group has several names
 464         { // create a group with subgroup grp and named grp.name
 465           SauvUtilities::Group* newGroup = _iMed->addNewGroup();
 466           newGroup->_groups.push_back( &_iMed->_groups[ grpID-1 ]);
 467           newGroup->_name = grp._name;
 468         }
 469       grp._name=objectNames[i];
 470 #ifdef _DEBUG
 471       map<int,int>::iterator it = strangeGroupType.find( grpID - 1 );
 472       if ( it != strangeGroupType.end() )
 473         cout << "Skip " << grp._name << " of not supported CASTEM type: " << it->second << endl;
 474 #endif
 475     }
 476 } // read_PILE_SOUS_MAILLAGE()
 477
 478 //================================================================================
 479 /*!
 480  * \brief Skip "PILE NUMERO  18" of XDR file
 481  */
 482 //================================================================================
 483
 484 void SauvReader::read_PILE_LREEL (const int nbObjects, std::vector<std::string>&, std::vector<int>&)
 485 {
 486   if ( isXRD() )
 487     {
 488       for (int object=0; object!=nbObjects; ++object) // pour chaque Group
 489         {
 490           initIntReading(1);
 491           int nb_vals = getIntNext();
 492           initDoubleReading(nb_vals);
 493           for(int i=0; i<nb_vals; i++) next();
 494         }
 495     }
 496 }
 497
 498 //================================================================================
 499 /*!
 500  * \brief Skip "PILE NUMERO  24" of XDR file
 501  */
 502 //================================================================================
 503
 504 void SauvReader::read_PILE_LOGIQUES (const int, std::vector<std::string>&, std::vector<int>&)
 505 {
 506   if ( isXRD() )
 507     {
 508       initIntReading(1);
 509       int nb_vals = getIntNext();
 510       initIntReading(nb_vals);
 511       for(int i=0; i<nb_vals; i++) next();
 512     }
 513 }
 514
 515 //================================================================================
 516 /*!
 517  * \brief Skip "PILE NUMERO  25" of XDR file
 518  */
 519 //================================================================================
 520
 521 void SauvReader::read_PILE_FLOATS (const int, std::vector<std::string>&, std::vector<int>&)
 522 {
 523   if ( isXRD() )
 524     {
 525       initIntReading(1);
 526       int nb_vals = getIntNext();
 527       initDoubleReading(nb_vals);
 528       for(int i=0; i<nb_vals; i++) next();
 529     }
 530 }
 531
 532 //================================================================================
 533 /*!
 534  * \brief Skip "PILE NUMERO  26" of XDR file
 535  */
 536 //================================================================================
 537
 538 void SauvReader::read_PILE_INTEGERS (const int, std::vector<std::string>&, std::vector<int>&)
 539 {
 540   if ( isXRD() )
 541     {
 542       initIntReading(1);
 543       int nb_vals = getIntNext();
 544       initIntReading(nb_vals);
 545       for(int i=0; i<nb_vals; i++) next();
 546     }
 547 }
 548
 549 //================================================================================
 550 /*!
 551  * \brief Skip "PILE NUMERO  29" of XDR file
 552  */
 553 //================================================================================
 554
 555 void SauvReader::read_PILE_LMOTS (const int nbObjects, std::vector<std::string>&, std::vector<int>&)
 556 {
 557   if ( isXRD() )
 558     {
 559       for (int object=0; object!=nbObjects; ++object) // pour chaque Group
 560         {
 561           initIntReading(2);
 562           int len = getIntNext();
 563           int nb_vals = getIntNext();
 564           int nb_char = len*nb_vals;
 565           int nb_char_tmp = 0;
 566           int fixed_length = 71;
 567           while (nb_char_tmp < nb_char)
 568             {
 569               int remain_len = nb_char - nb_char_tmp;
 570               int width;
 571               if ( remain_len > fixed_length )
 572                 {
 573                   width = fixed_length;
 574                 }
 575               else
 576                 {
 577                   width = remain_len;
 578                 }
 579               initNameReading(1, width);
 580               next();
 581               nb_char_tmp += width;
 582             }
 583         }
 584     }
 585 }
 586
 587 //================================================================================
 588 /*!
 589  * \brief Skip "PILE NUMERO  38" of XDR file
 590  */
 591 //================================================================================
 592
 593 void SauvReader::read_PILE_MODL (const int nbObjects, std::vector<std::string>&, std::vector<int>&)
 594 {
 595   if ( isXRD() )
 596     {
 597       for (int object=0; object!=nbObjects; ++object) // pour chaque Group
 598         {
 599           // see wrmodl.eso
 600           initIntReading(10);
 601           int n1  = getIntNext();
 602           int nm2 = getIntNext();
 603           int nm3 = getIntNext();
 604           int nm4 = getIntNext();
 605           int nm5 = getIntNext();
 606           int n45 = getIntNext();
 607           /*int nm6 =*/ getIntNext();
 608           /*int nm7 =*/ getIntNext();
 609           next();
 610           next();
 611           int nm1 = n1 * n45;
 612           int nm9 = n1 * 16;
 613           for (initIntReading(nm1); more(); next());
 614           for (initIntReading(nm9); more(); next());
 615           for (initNameReading(nm5, 8); more(); next());
 616           for (initNameReading(nm2, 8); more(); next());
 617           for (initNameReading(nm3, 8); more(); next());
 618           for (initIntReading(nm4); more(); next());
 619         }
 620     }
 621 } // Fin case pile 38
 622
 623 //================================================================================
 624 /*!
 625  * \brief Read "PILE NUMERO  32": links to node coordinates
 626  */
 627 //================================================================================
 628
 629 void SauvReader::read_PILE_NOEUDS (const int nbObjects, std::vector<std::string>&, std::vector<int>&)
 630 {
 631   initIntReading(1);
 632   int nb_indices = getIntNext();
 633
 634   if (nb_indices != nbObjects)
 635     THROW_IK_EXCEPTION("Error of reading PILE NUMERO  " << PILE_NOEUDS << lineNb() );
 636
 637   for ( initIntReading( nbObjects ); more(); next() )
 638     {
 639       int coordID = getInt();
 640       _iMed->getNode( index()+1 )->_coordID = coordID;
 641     }
 642 }
 643
 644 //================================================================================
 645 /*!
 646  * \brief Read "PILE NUMERO  33": node coordinates
 647  */
 648 //================================================================================
 649
 650 void SauvReader::read_PILE_COORDONNEES (const int nbObjects, std::vector<std::string>&, std::vector<int>&)
 651 {
 652   initIntReading(1);
 653   int nbReals = getIntNext();
 654
 655   if ( nbReals < (int)(_iMed->_nbNodes*(_iMed->_spaceDim+1)) )
 656     THROW_IK_EXCEPTION("Error of reading PILE NUMERO  " << PILE_COORDONNEES << lineNb() );
 657
 658   // there are coordinates + density for each node
 659   _iMed->_coords.resize( nbReals - nbReals/(_iMed->_spaceDim+1));
 660   double* coordPtr = &_iMed->_coords[0];
 661
 662   initDoubleReading( nbReals );
 663   while ( more() )
 664     {
 665       for (unsigned j = 0; j < _iMed->_spaceDim; ++j, next())
 666         *coordPtr++ = getDouble();
 667       // skip density
 668       getDouble();
 669       next();
 670     }
 671 }
 672
 673 //================================================================================
 674 /*!
 675  * \brief Find or create a Group equal to a given field support
 676  */
 677 //================================================================================
 678
 679 void SauvReader::setFieldSupport(const vector<SauvUtilities::Group*>& supports,
 680                                  SauvUtilities::DoubleField*          field)
 681 {
 682   SauvUtilities::Group* group = NULL;
 683   set<SauvUtilities::Group*> sup_set( supports.begin(), supports.end() );
 684   if ( sup_set.size() == 1 ) // one or equal supports
 685     {
 686       group = supports[0];
 687     }
 688   else
 689     {
 690       // check if sub-components are on cells of different types
 691       map<int,int> nbGaussByCellType;
 692       for ( size_t i = 0; i < supports.size(); ++i )
 693         {
 694           map<int,int>::iterator ct2ng = nbGaussByCellType.find( supports[i]->_cellType );
 695           if ( ct2ng == nbGaussByCellType.end() )
 696             nbGaussByCellType[ supports[i]->_cellType ] = field->_sub[i].nbGauss();
 697           else if ( ct2ng->second != field->_sub[i].nbGauss() )
 698             return;
 699         }
 700       bool isSameCellType = ( nbGaussByCellType.size() == 1 );
 701       // try to find an existing composite group with the same sub-groups
 702       if ( isSameCellType )
 703         for ( size_t i = 0; i < _iMed->_groups.size() && !group; ++i )
 704           {
 705             Group & grp = _iMed->_groups[i];
 706             if (sup_set.size() == grp._groups.size())
 707               {
 708                 bool sameOrder = true;
 709                 for ( size_t j = 0; j < supports.size() && sameOrder; ++j )
 710                   sameOrder = ( supports[j] == grp._groups[ j % grp._groups.size() ]);
 711                 if ( sameOrder )
 712                   group = & _iMed->_groups[i];
 713               }
 714           }
 715       if ( !group ) // no such a group, add a new one
 716         {
 717           vector<SauvUtilities::Group*> newGroups( supports.begin(),
 718                                                    supports.begin() + sup_set.size() );
 719           // check if supports includes newGroups in the same order
 720           bool sameOrder = true;
 721           for ( size_t j = newGroups.size(); j < supports.size() && sameOrder; ++j )
 722             sameOrder = ( supports[j] == newGroups[ j % newGroups.size() ]);
 723           if ( sameOrder )
 724             {
 725               group = _iMed->addNewGroup( & newGroups );
 726               group->_groups.swap( newGroups );
 727             }
 728         }
 729       // sort field sub-components and supports by cell type
 730       if ( group && !isSameCellType )
 731         {
 732           // sort groups
 733           vector<SauvUtilities::Group*>& groups = group->_groups;
 734           bool isModified = false, isSwapped = true;
 735           while ( isSwapped )
 736             {
 737               isSwapped = false;
 738               for ( size_t i = 1; i < groups.size(); ++i )
 739                 {
 740                   int nbN1 = groups[i-1]->empty() ? 0 : groups[i-1]->_cells[0]->_nodes.size();
 741                   int nbN2 = groups[i  ]->empty() ? 0 : groups[i  ]->_cells[0]->_nodes.size();
 742                   if ( nbN1 > nbN2 )
 743                     {
 744                       isSwapped = isModified = true;
 745                       std::swap( groups[i], groups[i-1] );
 746                     }
 747                 }
 748             }
 749           // relocate sub-components according to a new order of groups
 750           if ( isModified )
 751             {
 752               vector< DoubleField::_Sub_data > newSub   ( field->_sub.size() );
 753               vector< vector< double > >       newValues( field->_comp_values.size() );
 754               size_t iFromSub = 0, iNewSub = 0, iNewComp = 0;
 755               for ( ; iFromSub < field->_sub.size(); iFromSub += groups.size() )
 756                 {
 757                   size_t iFromComp = iNewComp;
 758                   for ( size_t iG = 0; iG < groups.size(); ++iG )
 759                     {
 760                       size_t iComp = iFromComp;
 761                       for ( size_t iSub = iFromSub; iSub < field->_sub.size(); ++iSub )
 762                         if ( field->_sub[ iSub ]._support == groups[ iG ] )
 763                           {
 764                             newSub[ iNewSub++ ] = field->_sub[ iSub ];
 765                             int iC = 0, nbC = field->_sub[ iSub ].nbComponents();
 766                             for ( ; iC < nbC; ++iC )
 767                               newValues[ iNewComp++ ].swap( field->_comp_values[ iComp++ ]);
 768                             break;
 769                           }
 770                         else
 771                           {
 772                             iComp += field->_sub[ iSub ].nbComponents();
 773                           }
 774                     }
 775                 }
 776               field->_sub.swap( newSub );
 777               field->_comp_values.swap( newValues );
 778             }
 779         }
 780     }
 781   if ( group )
 782     group->_isProfile = true;
 783
 784   field->_group = group;
 785 }
 786
 787 //================================================================================
 788 /*!
 789  * \brief Set field names
 790  */
 791 //================================================================================
 792
 793 void SauvReader::setFieldNames(const vector<SauvUtilities::DoubleField* >& fields,
 794                                const vector<string>&                       objets_nommes,
 795                                const vector<int>&                          indices_objets_nommes)
 796 {
 797   unsigned i;
 798   for ( i = 0; i < indices_objets_nommes.size(); ++i )
 799     {
 800       int fieldIndex = indices_objets_nommes[ i ];
 801       if ( fields[ fieldIndex - 1 ] )
 802         fields[ fieldIndex - 1 ]->_name = objets_nommes[ i ];
 803     }
 804 }
 805
 806 //================================================================================
 807 /*!
 808  * \brief Read "PILE NUMERO   2": NODE FIELDS
 809  */
 810 //================================================================================
 811
 812 void SauvReader::read_PILE_NODES_FIELD (const int                 nbObjects,
 813                                         std::vector<std::string>& objectNames,
 814                                         std::vector<int>&         nameIndices)
 815 {
 816   _iMed->_nodeFields.resize( nbObjects, (SauvUtilities::DoubleField*) 0 );
 817   for (int object=0; object!=nbObjects; ++object) // loop on fields
 818     {
 819       // EXAMPLE ( with no values )
 820
 821       // (1)       4       7       2       1
 822       // (2)     -88       0       3     -89       0       1     -90       0       2     -91
 823       // (2)       0       1
 824       // (3) FX   FY   FZ   FZ   FX   FY   FLX
 825       // (4)       0       0       0       0       0       0       0
 826       // (5)           cree  par  muc pri
 827       // (6)
 828       // (7)       2
 829
 830       // (1): nb subcomponents, nb components(total), IFOUR, nb attributes
 831       int nb_sub, total_nb_comp, nb_attr;
 832       int i_sub, i_comp;
 833       initIntReading( 4 );
 834       nb_sub        = getIntNext();
 835       total_nb_comp = getIntNext();
 836       next(); // ignore IFOUR
 837       nb_attr       = getIntNext();
 838       if ( nb_sub < 0 || total_nb_comp < 0 || nb_attr < 0 )
 839         THROW_IK_EXCEPTION("Error of field reading " << lineNb());
 840
 841       // (2) loop on subcomponents of a field, for each read
 842       // (a) support, (b) number of values and (c) number of components
 843       vector<Group*> supports( nb_sub );
 844       vector<int> nb_values  ( nb_sub );
 845       vector<int> nb_comps   ( nb_sub );
 846       int total_nb_values = 0;
 847       initIntReading( nb_sub * 3 );
 848       for ( i_sub = 0; i_sub < nb_sub; ++i_sub )
 849         {
 850           int supId = -getIntNext(); // (a) reference to support
 851           if ( supId < 1 || supId > (int)_iMed->_groups.size() )
 852             THROW_IK_EXCEPTION("Wrong mesh reference: "<< supId << lineNb() );
 853           supports[ i_sub ] = &_iMed->_groups[ supId-1 ]; // (a) reference to support
 854
 855           nb_values[ i_sub ] = getIntNext();    // (b) nb points
 856           total_nb_values += nb_values[ i_sub ];
 857           if ( nb_values[ i_sub ] < 0 )
 858             THROW_IK_EXCEPTION(" Wrong nb of points: " << nb_values[ i_sub ]  << lineNb() );
 859           nb_comps[ i_sub ] = getInt(); next();     // (c) nb of components in i_sub
 860         }
 861
 862       // create a field if there are values
 863       SauvUtilities::DoubleField* fdouble = 0;
 864       if ( total_nb_values > 0 )
 865         fdouble = new DoubleField( nb_sub, total_nb_comp );
 866       _iMed->_nodeFields[ object ] = fdouble;
 867
 868       // (3) component names
 869       initNameReading( total_nb_comp, 4 );
 870       for ( i_sub = 0; i_sub < nb_sub; ++i_sub )
 871         {
 872           // store support id and nb components of a sub
 873           if ( fdouble )
 874             fdouble->_sub[ i_sub ].setData( nb_comps[ i_sub ], supports[ i_sub ] );
 875           for ( i_comp = 0; i_comp < nb_comps[ i_sub ]; ++i_comp, next() )
 876             {
 877               // store component name
 878               string compName = getName();
 879               if ( fdouble )
 880                 fdouble->_sub[ i_sub ].compName( i_comp ) = compName;
 881             }
 882         }
 883       // (4) nb harmonics ( ignored )
 884       for ( initIntReading( total_nb_comp ); more(); next() );
 885       // (5) TYPE ( ignored )
 886       for (initNameReading(1, /*length=*/71); more(); next());
 887       // (6) TITRE ( ignored )
 888       for (initNameReading(1, /*length=*/71); more(); next());
 889       // (7) attributes ( ignored )
 890       for ( initIntReading( nb_attr ); more(); next() );
 891
 892       for ( i_sub = 0; i_sub < nb_sub; ++i_sub )
 893         {
 894           // loop on components: read values
 895           initDoubleReading( nb_values[ i_sub ] * nb_comps[ i_sub ] );
 896           for ( i_comp = 0; i_comp < nb_comps[ i_sub ]; ++i_comp )
 897             {
 898               if ( fdouble )
 899                 {
 900                   vector<double>& vals = fdouble->addComponent( nb_values[ i_sub ] );
 901                   for ( int i = 0; more() && i < nb_values[ i_sub ]; next(), ++i )
 902                     vals[ i ] = getDouble();
 903                 }
 904               else
 905                 {
 906                   for ( int i = 0; i < nb_values[ i_sub ]; next(), ++i );
 907                 }
 908             }
 909         } // loop on subcomponents of a field
 910
 911       // set a supporting group including all subs supports but only
 912       // if all subs have the same components
 913       if ( fdouble && fdouble->hasSameComponentsBySupport() )
 914         setFieldSupport( supports, fdouble );
 915       else
 916         for ( i_sub = 0; i_sub < nb_sub; ++i_sub )
 917           fdouble->_sub[ i_sub ]._support->_isProfile = true;
 918
 919     } // end loop on field objects
 920
 921   // set field names
 922   setFieldNames( _iMed->_nodeFields, objectNames, nameIndices );
 923
 924 }  // read_PILE_NODES_FIELD()
 925
 926 //================================================================================
 927 /*!
 928  * \brief Read "PILE NUMERO  39": FIELDS
 929  */
 930 //================================================================================
 931
 932 void SauvReader::read_PILE_FIELD (const int                 nbObjects,
 933                                   std::vector<std::string>& objectNames,
 934                                   std::vector<int>&         nameIndices)
 935 {
 936   // REAL EXAMPLE
 937
 938   // (1)        1       2       6      16
 939   // (2)                                                         CARACTERISTIQUES
 940   // (3)      -15  317773       4       0       0       0      -2       0       3
 941   // (4)             317581
 942   // (5)  0
 943   // (6)   317767  317761  317755  317815
 944   // (7)  YOUN     NU       H        SIGY
 945   // (8)  REAL*8            REAL*8            REAL*8            REAL*8
 946   // (9)        1       1       0       0
 947   // (10)  2.00000000000000E+05
 948   // (11)       1       1       0       0
 949   // (12)  3.30000000000000E-01
 950   // (13)       1       1       0       0
 951   // (14)  1.00000000000000E+04
 952   // (15)       6     706       0       0
 953   // (16)  1.00000000000000E+02  1.00000000000000E+02  1.00000000000000E+02
 954   // (17)  1.00000000000000E+02  1.00000000000000E+02  1.00000000000000E+02
 955   // (18)  ...
 956
 957   _iMed->_cellFields.resize( nbObjects, (SauvUtilities::DoubleField*) 0 );
 958   for (int object=0; object!=nbObjects; ++object) // pour chaque field
 959     {
 960       initIntReading( 4 );
 961       int i_sub, nb_sub = getIntNext(); // (1) <nb_sub> 2 6 <title length>
 962       next(); // skip "2"
 963       next(); // skip "6"
 964       int title_length = getIntNext(); // <title length>
 965       if ( nb_sub < 1 )
 966         THROW_IK_EXCEPTION("Error of field reading: wrong nb of subcomponents " << nb_sub << lineNb() );
 967
 968       string description;
 969       if ( title_length )
 970         {
 971           if ( isXRD() )
 972             {
 973               initNameReading(1, title_length);
 974               description = getName();
 975               next();
 976             }
 977           else
 978             {
 979               char* line;
 980               getNextLine( line ); // (2) title
 981               const int len = 72; // line length
 982               description = string(line + len - title_length, title_length); // title is right justified
 983             }
 984         }
 985       // look for a line starting with '-' : <reference to support>
 986       if ( isXRD() )
 987         {
 988           initIntReading( nb_sub * 9 );
 989         }
 990       else
 991         {
 992           do {
 993             initIntReading( nb_sub * 9 );
 994           } while ( getInt() >= 0 );
 995         }
 996       int total_nb_comp = 0;
 997       vector<Group*> supports( nb_sub );
 998       vector<int>     nb_comp( nb_sub );
 999       for ( i_sub = 0; i_sub < nb_sub; ++i_sub )
1000         {                                    // (3)
1001           int supportId     = -getIntNext(); // <reference to support>
1002           next();                            // ignore <address>
1003           nb_comp [ i_sub ] =  getIntNext(); // <nb of components in the sub>
1004           for ( int i = 0; i < 6; ++i ) next();  // ignore 6 ints, in example "0 0 0 -2 0 3"
1005
1006           if ( supportId < 1 || supportId > (int)_iMed->_groups.size() )
1007             THROW_IK_EXCEPTION("Error of field reading: wrong mesh reference "<< supportId << lineNb() );
1008           if ( nb_comp[ i_sub ] < 0 )
1009             THROW_IK_EXCEPTION("Error of field reading: wrong nb of components " <<nb_comp[ i_sub ] << lineNb() );
1010
1011           supports[ i_sub ] = &_iMed->_groups[ supportId-1 ];
1012           total_nb_comp += nb_comp[ i_sub ];
1013         }
1014       // (4) dummy strings
1015       for ( initNameReading( nb_sub, 17 ); more(); next() );
1016       // (5) dummy strings
1017       for ( initNameReading( nb_sub ); more(); next() );
1018
1019       // loop on subcomponents of a field, each of which refers to
1020       // a certain support and has its own number of components;
1021       // read component values
1022       SauvUtilities::DoubleField* fdouble = 0;
1023       for ( i_sub = 0; i_sub < nb_sub; ++ i_sub )
1024         {
1025           vector<string> comp_names( nb_comp[ i_sub ]), comp_type( nb_comp[ i_sub ]);
1026           // (6) nb_comp addresses of MELVAL structure
1027           for ( initIntReading( nb_comp[ i_sub ] ); more(); next() );
1028           // (7) component names
1029           for ( initNameReading( nb_comp[ i_sub ] ); more(); next() )
1030             comp_names[ index() ] = getName();
1031           // (8) component type
1032           for ( initNameReading( nb_comp[ i_sub ], 17 ); more(); next() ) // 17 is name width
1033             {
1034               comp_type[ index() ] = getName();
1035               // component types must be the same
1036               if ( index() > 0 && comp_type[ index() ] != comp_type[ index() - 1] )
1037                 THROW_IK_EXCEPTION( "Error of field reading: diff component types <"
1038                                     << comp_type[ index() ] << "> != <" << comp_type[ index() - 1 ]
1039                                     << ">" << lineNb() );
1040             }
1041           // now type is known, create a field, one for all subs
1042           bool isReal = (nb_comp[i_sub] > 0) ? (comp_type[0] == "REAL*8") : true;
1043           if ( !fdouble && total_nb_comp )
1044             {
1045               if ( !isReal )
1046                 cout << "Warning: read NOT REAL field, type <" << comp_type[0] << ">" << lineNb() << endl;
1047               _iMed->_cellFields[ object ] = fdouble = new SauvUtilities::DoubleField( nb_sub, total_nb_comp );
1048               fdouble->_description = description;
1049             }
1050           // store support id and nb components of a sub
1051           if ( fdouble )
1052             fdouble->_sub[ i_sub ].setData( nb_comp[ i_sub ], supports[ i_sub ]);
1053           // loop on components: read values
1054           for ( int i_comp = 0; i_comp < nb_comp[ i_sub ]; ++i_comp )
1055             {
1056               // (9) nb of values
1057               initIntReading( 4 );
1058               int nb_val_by_elem = getIntNext();
1059               int nb_values      = getIntNext();
1060               next();
1061               next();
1062               fdouble->_sub[ i_sub ]._nb_gauss[ i_comp ] = nb_val_by_elem;
1063
1064               // (10) values
1065               nb_values *= nb_val_by_elem;
1066               if ( fdouble )
1067                 {
1068                   vector<double> & vals = fdouble->addComponent( nb_values );
1069                   for ( isReal ? initDoubleReading( nb_values ) : initIntReading( nb_values ); more(); next())
1070                     vals[ index() ] = getDouble();
1071                   // store component name
1072                   fdouble->_sub[ i_sub ].compName( i_comp ) = comp_names[ i_comp ];
1073                 }
1074               else
1075                 {
1076                   for ( isReal ? initDoubleReading( nb_values ) : initIntReading( nb_values ); more(); next() ) ;
1077                 }
1078             }
1079         } // loop on subcomponents of a field
1080
1081       // set id of a group including all sub supports but only
1082       // if all subs have the same nb of components
1083       if ( fdouble && fdouble->hasSameComponentsBySupport() )
1084         setFieldSupport( supports, fdouble );
1085       else
1086         for ( i_sub = 0; i_sub < nb_sub; ++i_sub )
1087           fdouble->_sub[ i_sub ]._support->_isProfile = true;
1088
1089     } // end loop on field objects
1090
1091   // set field names
1092   setFieldNames( _iMed->_cellFields, objectNames, nameIndices );
1093
1094 } // read_PILE_FIELD()
1095
1096 //================================================================================
1097 /*!
1098  * \brief Read "PILE NUMERO  10": TABLES
1099  */
1100 //================================================================================
1101
1102 void SauvReader::read_PILE_TABLES (const int                 nbObjects,
1103                                    std::vector<std::string>& objectNames,
1104                                    std::vector<int>&         nameIndices)
1105 {
1106   // IMP 0020434: mapping GIBI names to MED names
1107
1108   string table_med_mail = "MED_MAIL";
1109   string table_med_cham = "MED_CHAM";
1110   string table_med_comp = "MED_COMP";
1111   int table_med_mail_id = -1;
1112   int table_med_cham_id = -1;
1113   int table_med_comp_id = -1;
1114   for (size_t iname = 0; iname < objectNames.size(); iname++)
1115     if      (objectNames[iname] == table_med_mail) table_med_mail_id = nameIndices[iname];
1116     else if (objectNames[iname] == table_med_cham) table_med_cham_id = nameIndices[iname];
1117     else if (objectNames[iname] == table_med_comp) table_med_comp_id = nameIndices[iname];
1118
1119   if ( isASCII() )
1120     if (table_med_mail_id < 0 && table_med_cham_id < 0 && table_med_comp_id < 0)
1121       return;
1122
1123   for (int itable = 1; itable <= nbObjects; itable++)
1124     {
1125       // read tables "MED_MAIL", "MED_CHAM" and "MED_COMP", that keeps correspondence
1126       // between GIBI names (8 symbols if any) and MED names (possibly longer)
1127       initIntReading(1);
1128       int nb_table_vals = getIntNext();
1129       if (nb_table_vals < 0)
1130         THROW_IK_EXCEPTION("Error of reading PILE NUMERO  10" << lineNb() );
1131
1132       int name_i_med_pile;
1133       initIntReading(nb_table_vals);
1134       for (int i = 0; i < nb_table_vals/4; i++)
1135         {
1136           if (itable == table_med_mail_id ||
1137               itable == table_med_cham_id ||
1138               itable == table_med_comp_id)
1139             {
1140               nameGIBItoMED name_i;
1141               name_i_med_pile  = getIntNext();
1142               name_i.med_id    = getIntNext();
1143               name_i.gibi_pile = getIntNext();
1144               name_i.gibi_id   = getIntNext();
1145
1146               if (name_i_med_pile != PILE_STRINGS)
1147                 {
1148                   // error: med(long) names are always kept in PILE_STRINGS
1149                 }
1150               if (itable == table_med_mail_id)
1151                 {
1152                   if (name_i.gibi_pile != PILE_SOUS_MAILLAGE) {
1153                     // error: must be PILE_SOUS_MAILLAGE
1154                   }
1155                   _iMed->_listGIBItoMED_mail.push_back(name_i);
1156                 }
1157               else if (itable == table_med_cham_id)
1158                 {
1159                   if (name_i.gibi_pile != PILE_FIELD &&
1160                       name_i.gibi_pile != PILE_NODES_FIELD)
1161                     {
1162                       // error: must be PILE_FIELD or PILE_NODES_FIELD
1163                     }
1164                   _iMed->_listGIBItoMED_cham.push_back(name_i);
1165                 }
1166               else if (itable == table_med_comp_id)
1167                 {
1168                   if (name_i.gibi_pile != PILE_STRINGS)
1169                     {
1170                       // error: gibi(short) names of components are kept in PILE_STRINGS
1171                     }
1172                   _iMed->_listGIBItoMED_comp.push_back(name_i);
1173                 }
1174             }
1175           else
1176             {
1177               // pass table
1178               for ( int ii = 0; ii < 4; ++ii ) next();
1179             }
1180         }
1181     } // for (int itable = 0; itable < nbObjects; itable++)
1182 }
1183
1184 //================================================================================
1185 /*!
1186  * \brief Read "PILE NUMERO  27"
1187  */
1188 //================================================================================
1189
1190 void SauvReader::read_PILE_STRINGS (const int                 nbObjects,
1191                                     std::vector<std::string>& objectNames,
1192                                     std::vector<int>&         nameIndices)
1193 {
1194   // IMP 0020434: mapping GIBI names to MED names
1195   initIntReading(2);
1196   int stringLen    = getIntNext();
1197   int nbSubStrings = getIntNext();
1198   if (nbSubStrings != nbObjects)
1199     THROW_IK_EXCEPTION("Error of reading PILE NUMERO  27" << lineNb() );
1200
1201   string aWholeString;
1202   if ( isXRD() )
1203     {
1204       const int fixedLength = 71;
1205       while ((int)aWholeString.length() < stringLen)
1206         {
1207           int remainLen = stringLen - aWholeString.length();
1208           int len;
1209           if ( remainLen > fixedLength )
1210             {
1211               len = fixedLength;
1212             }
1213           else
1214             {
1215               len = remainLen;
1216             }
1217           initNameReading(1, len);
1218           aWholeString += getName();
1219           next();
1220         }
1221     }
1222   else
1223     {
1224       char* line;
1225       const int fixedLength = 71;
1226       while ((int)aWholeString.length() < stringLen)
1227         {
1228           getNextLine( line );
1229           int remainLen = stringLen - aWholeString.length();
1230           if ( remainLen > fixedLength )
1231             {
1232               aWholeString += line + 1;
1233             }
1234           else
1235             {
1236               aWholeString += line + ( 72 - remainLen );
1237             }
1238         }
1239     }
1240   int prevOffset = 0;
1241   int currOffset = 0;
1242   initIntReading(nbSubStrings);
1243   for (int istr = 1; istr <= nbSubStrings; istr++, next())
1244     {
1245       currOffset = getInt();
1246       // fill mapStrings
1247       _iMed->_mapStrings[istr] = aWholeString.substr(prevOffset, currOffset - prevOffset);
1248       prevOffset = currOffset;
1249     }
1250 }