medtool/src/MEDLoader/SauvReader.cxx

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