src/MEDLoader/SauvReader.cxx

   1 // Copyright (C) 2007-2013  CEA/DEN, EDF R&D
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19 // File      : 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 char *fileName)
  48 {
  49   if ( !fileName || strlen(fileName) < 1 ) 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 );
  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 );
  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::getDirectChildren() 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 composit 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           _iMed->_groups.push_back(Group());
 432           _iMed->_groups.back()._groups.push_back( &_iMed->_groups[ grpID-1 ]);
 433           _iMed->_groups.back()._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 Finds or create a Group equal to a given field support
 642  */
 643 //================================================================================
 644
 645 SauvUtilities::Group* SauvReader::getFieldSupport(const vector<SauvUtilities::Group*>& supports)
 646 {
 647   SauvUtilities::Group* group = NULL;
 648   set<SauvUtilities::Group*> sup_set( supports.begin(), supports.end() );
 649   if (sup_set.size() == 1 ) // one or equal supports
 650     {
 651       group = supports[0];
 652     }
 653   else
 654     {
 655       // try to find an existing composite group with the same sub-groups
 656       for ( size_t i = 0; i < _iMed->_groups.size() && !group; ++i )
 657         {
 658           Group & grp = _iMed->_groups[i];
 659           if (sup_set.size() == grp._groups.size())
 660             {
 661               bool sameOrder = true;
 662               for ( size_t j = 0; j < supports.size() && sameOrder; ++j )
 663                 sameOrder = ( supports[j] == grp._groups[ j % grp._groups.size() ]);
 664               if ( sameOrder )
 665                 group = & _iMed->_groups[i];
 666             }
 667         }
 668       if ( !group ) // no such a group, add a new one
 669         {
 670           vector<SauvUtilities::Group*> newGroups( supports.begin(),
 671                                                    supports.begin() + sup_set.size() );
 672           // check if supports includes newGroups in the same order
 673           bool sameOrder = true;
 674           for ( size_t j = newGroups.size(); j < supports.size() && sameOrder; ++j )
 675             sameOrder = ( supports[j] == newGroups[ j % newGroups.size() ]);
 676           if ( sameOrder )
 677             {
 678               _iMed->_groups.push_back( SauvUtilities::Group() );
 679               group = & _iMed->_groups.back();
 680               group->_groups.swap( newGroups );
 681             }
 682         }
 683     }
 684   if ( group )
 685     group->_isProfile = true;
 686   return group;
 687 }
 688
 689 //================================================================================
 690 /*!
 691  * \brief set field names
 692  */
 693 //================================================================================
 694
 695 void SauvReader::setFieldNames(const vector<SauvUtilities::DoubleField* >& fields,
 696                                const vector<string>&                       objets_nommes,
 697                                const vector<int>&                          indices_objets_nommes)
 698 {
 699   unsigned i;
 700   for ( i = 0; i < indices_objets_nommes.size(); ++i )
 701     {
 702       int fieldIndex = indices_objets_nommes[ i ];
 703       if ( fields[ fieldIndex - 1 ] )
 704         fields[ fieldIndex - 1 ]->_name = objets_nommes[ i ];
 705     }
 706 }
 707
 708 //================================================================================
 709 /*!
 710  * \brief Read "PILE NUMERO   2": NODE FIELDS
 711  */
 712 //================================================================================
 713
 714 void SauvReader::read_PILE_NODES_FIELD (const int                 nbObjects,
 715                                         std::vector<std::string>& objectNames,
 716                                         std::vector<int>&         nameIndices)
 717 {
 718   _iMed->_nodeFields.resize( nbObjects, (SauvUtilities::DoubleField*) 0 );
 719   for (int object=0; object!=nbObjects; ++object) // loop on fields
 720     {
 721       // EXAMPLE ( with no values )
 722
 723       // (1)       4       7       2       1
 724       // (2)     -88       0       3     -89       0       1     -90       0       2     -91
 725       // (2)       0       1
 726       // (3) FX   FY   FZ   FZ   FX   FY   FLX
 727       // (4)       0       0       0       0       0       0       0
 728       // (5)           cree  par  muc pri
 729       // (6)
 730       // (7)       2
 731
 732       // (1): nb subcomponents, nb components(total), IFOUR, nb attributes
 733       int nb_sub, total_nb_comp, nb_attr;
 734       int i_sub, i_comp;
 735       initIntReading( 4 );
 736       nb_sub        = getIntNext();
 737       total_nb_comp = getIntNext();
 738       next(); // ignore IFOUR
 739       nb_attr       = getIntNext();
 740       if ( nb_sub < 0 || total_nb_comp < 0 || nb_attr < 0 )
 741         THROW_IK_EXCEPTION("Error of field reading " << lineNb());
 742
 743       // (2) loop on subcomponents of a field, for each read
 744       // (a) support, (b) number of values and (c) number of components
 745       vector<Group*> supports( nb_sub );
 746       vector<int> nb_values  ( nb_sub );
 747       vector<int> nb_comps   ( nb_sub );
 748       int total_nb_values = 0;
 749       initIntReading( nb_sub * 3 );
 750       for ( i_sub = 0; i_sub < nb_sub; ++i_sub )
 751         {
 752           int supId = -getIntNext(); // (a) reference to support
 753           if ( supId < 1 || supId > (int)_iMed->_groups.size() )
 754             THROW_IK_EXCEPTION("Wrong mesh reference: "<< supId << lineNb() );
 755           supports[ i_sub ] = &_iMed->_groups[ supId-1 ]; // (a) reference to support
 756
 757           nb_values[ i_sub ] = getIntNext();    // (b) nb points
 758           total_nb_values += nb_values[ i_sub ];
 759           if ( nb_values[ i_sub ] < 0 )
 760             THROW_IK_EXCEPTION(" Wrong nb of points: " << nb_values[ i_sub ]  << lineNb() );
 761           nb_comps[ i_sub ] = getInt(); next();     // (c) nb of components in i_sub
 762         }
 763
 764       // create a field if there are values
 765       SauvUtilities::DoubleField* fdouble = 0;
 766       if ( total_nb_values > 0 )
 767         fdouble = new DoubleField( nb_sub, total_nb_comp );
 768       _iMed->_nodeFields[ object ] = fdouble;
 769
 770       // (3) component names
 771       initNameReading( total_nb_comp, 4 );
 772       for ( i_sub = 0; i_sub < nb_sub; ++i_sub )
 773         {
 774           // store support id and nb components of a sub
 775           if ( fdouble )
 776             fdouble->_sub[ i_sub ].setData( nb_comps[ i_sub ], supports[ i_sub ] );
 777           for ( i_comp = 0; i_comp < nb_comps[ i_sub ]; ++i_comp, next() )
 778             {
 779               // store component name
 780               string compName = getName();
 781               if ( fdouble )
 782                 fdouble->_sub[ i_sub ].compName( i_comp ) = compName;
 783             }
 784         }
 785       // (4) nb harmonics ( ignored )
 786       for ( initIntReading( total_nb_comp ); more(); next() );
 787       // (5) TYPE ( ignored )
 788       for (initNameReading(1, /*length=*/71); more(); next());
 789       // (6) TITRE ( ignored )
 790       for (initNameReading(1, /*length=*/71); more(); next());
 791       // (7) attributes ( ignored )
 792       for ( initIntReading( nb_attr ); more(); next() );
 793
 794       for ( i_sub = 0; i_sub < nb_sub; ++i_sub )
 795         {
 796           // loop on components: read values
 797           initDoubleReading( nb_values[ i_sub ] * nb_comps[ i_sub ] );
 798           for ( i_comp = 0; i_comp < nb_comps[ i_sub ]; ++i_comp )
 799             {
 800               if ( fdouble )
 801                 {
 802                   vector<double>& vals = fdouble->addComponent( nb_values[ i_sub ] );
 803                   for ( int i = 0; more() && i < nb_values[ i_sub ]; next(), ++i )
 804                     vals[ i ] = getDouble();
 805                 }
 806               else
 807                 {
 808                   for ( int i = 0; i < nb_values[ i_sub ]; next(), ++i );
 809                 }
 810             }
 811         } // loop on subcomponents of a field
 812
 813       // set a supporting group including all subs supports but only
 814       // if all subs have the same components
 815       if ( fdouble && fdouble->hasSameComponentsBySupport() )
 816         fdouble->_group = getFieldSupport( supports );
 817       else
 818         for ( i_sub = 0; i_sub < nb_sub; ++i_sub )
 819           fdouble->_sub[ i_sub ]._support->_isProfile = true;
 820
 821     } // end loop on field objects
 822
 823   // set field names
 824   setFieldNames( _iMed->_nodeFields, objectNames, nameIndices );
 825
 826 }  // read_PILE_NODES_FIELD()
 827
 828 //================================================================================
 829 /*!
 830  * \brief Read "PILE NUMERO  39": FIELDS
 831  */
 832 //================================================================================
 833
 834 void SauvReader::read_PILE_FIELD (const int                 nbObjects,
 835                                   std::vector<std::string>& objectNames,
 836                                   std::vector<int>&         nameIndices)
 837 {
 838   // REAL EXAMPLE
 839
 840   // (1)        1       2       6      16
 841   // (2)                                                         CARACTERISTIQUES
 842   // (3)      -15  317773       4       0       0       0      -2       0       3
 843   // (4)             317581
 844   // (5)  0
 845   // (6)   317767  317761  317755  317815
 846   // (7)  YOUN     NU       H        SIGY
 847   // (8)  REAL*8            REAL*8            REAL*8            REAL*8
 848   // (9)        1       1       0       0
 849   // (10)  2.00000000000000E+05
 850   // (11)       1       1       0       0
 851   // (12)  3.30000000000000E-01
 852   // (13)       1       1       0       0
 853   // (14)  1.00000000000000E+04
 854   // (15)       6     706       0       0
 855   // (16)  1.00000000000000E+02  1.00000000000000E+02  1.00000000000000E+02
 856   // (17)  1.00000000000000E+02  1.00000000000000E+02  1.00000000000000E+02
 857   // (18)  ...
 858
 859   _iMed->_cellFields.resize( nbObjects, (SauvUtilities::DoubleField*) 0 );
 860   for (int object=0; object!=nbObjects; ++object) // pour chaque field
 861     {
 862       initIntReading( 4 );
 863       int i_sub, nb_sub = getIntNext(); // (1) <nb_sub> 2 6 <title length>
 864       next(); // skip "2"
 865       next(); // skip "6"
 866       int title_length = getIntNext(); // <title length>
 867       if ( nb_sub < 1 )
 868         THROW_IK_EXCEPTION("Error of field reading: wrong nb of subcomponents " << nb_sub << lineNb() );
 869
 870       string description;
 871       if ( title_length )
 872         {
 873           if ( isXRD() )
 874             {
 875               initNameReading(1, title_length);
 876               description = getName();
 877               next();
 878             }
 879           else
 880             {
 881               char* line;
 882               getNextLine( line ); // (2) title
 883               const int len = 72; // line length
 884               description = string(line + len - title_length, title_length); // title is right justified
 885             }
 886         }
 887       // look for a line starting with '-' : <reference to support>
 888       if ( isXRD() )
 889         {
 890           initIntReading( nb_sub * 9 );
 891         }
 892       else
 893         {
 894           do {
 895             initIntReading( nb_sub * 9 );
 896           } while ( getInt() >= 0 );
 897         }
 898       int total_nb_comp = 0;
 899       vector<Group*> supports( nb_sub );
 900       vector<int>     nb_comp( nb_sub );
 901       for ( i_sub = 0; i_sub < nb_sub; ++i_sub )
 902         {                                    // (3)
 903           int supportId     = -getIntNext(); // <reference to support>
 904           next();                            // ignore <address>
 905           nb_comp [ i_sub ] =  getIntNext(); // <nb of components in the sub>
 906           for ( int i = 0; i < 6; ++i ) next();  // ignore 6 ints, in example "0 0 0 -2 0 3"
 907
 908           if ( supportId < 1 || supportId > (int)_iMed->_groups.size() )
 909             THROW_IK_EXCEPTION("Error of field reading: wrong mesh reference "<< supportId << lineNb() );
 910           if ( nb_comp[ i_sub ] < 0 )
 911             THROW_IK_EXCEPTION("Error of field reading: wrong nb of components " <<nb_comp[ i_sub ] << lineNb() );
 912
 913           supports[ i_sub ] = &_iMed->_groups[ supportId-1 ];
 914           total_nb_comp += nb_comp[ i_sub ];
 915         }
 916       // (4) dummy strings
 917       for ( initNameReading( nb_sub, 17 ); more(); next() );
 918       // (5) dummy strings
 919       for ( initNameReading( nb_sub ); more(); next() );
 920
 921       // loop on subcomponents of a field, each of which refers to
 922       // a certain support and has its own number of components;
 923       // read component values
 924       SauvUtilities::DoubleField* fdouble = 0;
 925       for ( i_sub = 0; i_sub < nb_sub; ++ i_sub )
 926         {
 927           vector<string> comp_names( nb_comp[ i_sub ]), comp_type( nb_comp[ i_sub ]);
 928           // (6) nb_comp addresses of MELVAL structure
 929           for ( initIntReading( nb_comp[ i_sub ] ); more(); next() );
 930           // (7) component names
 931           for ( initNameReading( nb_comp[ i_sub ] ); more(); next() )
 932             comp_names[ index() ] = getName();
 933           // (8) component type
 934           for ( initNameReading( nb_comp[ i_sub ], 17 ); more(); next() ) // 17 is name width
 935             {
 936               comp_type[ index() ] = getName();
 937               // component types must be the same
 938               if ( index() > 0 && comp_type[ index() ] != comp_type[ index() - 1] )
 939                 THROW_IK_EXCEPTION( "Error of field reading: diff component types <"
 940                                     << comp_type[ index() ] << "> != <" << comp_type[ index() - 1 ]
 941                                     << ">" << lineNb() );
 942             }
 943           // now type is known, create a field, one for all subs
 944           bool isReal = (nb_comp[i_sub] > 0) ? (comp_type[0] == "REAL*8") : true;
 945           if ( !fdouble && total_nb_comp )
 946             {
 947               if ( !isReal )
 948                 cout << "Warning: read NOT REAL field, type <" << comp_type[0] << ">" << lineNb() << endl;
 949               _iMed->_cellFields[ object ] = fdouble = new SauvUtilities::DoubleField( nb_sub, total_nb_comp );
 950               fdouble->_description = description;
 951             }
 952           // store support id and nb components of a sub
 953           if ( fdouble )
 954             fdouble->_sub[ i_sub ].setData( nb_comp[ i_sub ], supports[ i_sub ]);
 955           // loop on components: read values
 956           for ( int i_comp = 0; i_comp < nb_comp[ i_sub ]; ++i_comp )
 957             {
 958               // (9) nb of values
 959               initIntReading( 4 );
 960               int nb_val_by_elem = getIntNext();
 961               int nb_values      = getIntNext();
 962               next();
 963               next();
 964               fdouble->_sub[ i_sub ]._nb_gauss[ i_comp ] = nb_val_by_elem;
 965
 966               // (10) values
 967               nb_values *= nb_val_by_elem;
 968               if ( fdouble )
 969                 {
 970                   vector<double> & vals = fdouble->addComponent( nb_values );
 971                   for ( isReal ? initDoubleReading( nb_values ) : initIntReading( nb_values ); more(); next())
 972                     vals[ index() ] = getDouble();
 973                   // store component name
 974                   fdouble->_sub[ i_sub ].compName( i_comp ) = comp_names[ i_comp ];
 975                 }
 976               else
 977                 {
 978                   for ( isReal ? initDoubleReading( nb_values ) : initIntReading( nb_values ); more(); next() ) ;
 979                 }
 980             }
 981         } // loop on subcomponents of a field
 982
 983       // set id of a group including all sub supports but only
 984       // if all subs have the same nb of components
 985       if ( fdouble && fdouble->hasSameComponentsBySupport() )
 986         fdouble->_group = getFieldSupport( supports );
 987       else
 988         for ( i_sub = 0; i_sub < nb_sub; ++i_sub )
 989           fdouble->_sub[ i_sub ]._support->_isProfile = true;
 990
 991     } // end loop on field objects
 992
 993   // set field names
 994   setFieldNames( _iMed->_cellFields, objectNames, nameIndices );
 995
 996 } // read_PILE_FIELD()
 997
 998 //================================================================================
 999 /*!
1000  * \brief Read "PILE NUMERO  10": TABLES
1001  */
1002 //================================================================================
1003
1004 void SauvReader::read_PILE_TABLES (const int                 nbObjects,
1005                                    std::vector<std::string>& objectNames,
1006                                    std::vector<int>&         nameIndices)
1007 {
1008   // IMP 0020434: mapping GIBI names to MED names
1009
1010   string table_med_mail = "MED_MAIL";
1011   string table_med_cham = "MED_CHAM";
1012   string table_med_comp = "MED_COMP";
1013   int table_med_mail_id = -1;
1014   int table_med_cham_id = -1;
1015   int table_med_comp_id = -1;
1016   for (size_t iname = 0; iname < objectNames.size(); iname++)
1017     if      (objectNames[iname] == table_med_mail) table_med_mail_id = nameIndices[iname];
1018     else if (objectNames[iname] == table_med_cham) table_med_cham_id = nameIndices[iname];
1019     else if (objectNames[iname] == table_med_comp) table_med_comp_id = nameIndices[iname];
1020
1021   if ( isASCII() )
1022     if (table_med_mail_id < 0 && table_med_cham_id < 0 && table_med_comp_id < 0)
1023       return;
1024
1025   for (int itable = 1; itable <= nbObjects; itable++)
1026     {
1027       // read tables "MED_MAIL", "MED_CHAM" and "MED_COMP", that keeps correspondence
1028       // between GIBI names (8 symbols if any) and MED names (possibly longer)
1029       initIntReading(1);
1030       int nb_table_vals = getIntNext();
1031       if (nb_table_vals < 0)
1032         THROW_IK_EXCEPTION("Error of reading PILE NUMERO  10" << lineNb() );
1033
1034       int name_i_med_pile;
1035       initIntReading(nb_table_vals);
1036       for (int i = 0; i < nb_table_vals/4; i++)
1037         {
1038           if (itable == table_med_mail_id ||
1039               itable == table_med_cham_id ||
1040               itable == table_med_comp_id)
1041             {
1042               nameGIBItoMED name_i;
1043               name_i_med_pile  = getIntNext();
1044               name_i.med_id    = getIntNext();
1045               name_i.gibi_pile = getIntNext();
1046               name_i.gibi_id   = getIntNext();
1047
1048               if (name_i_med_pile != PILE_STRINGS)
1049                 {
1050                   // error: med(long) names are always kept in PILE_STRINGS
1051                 }
1052               if (itable == table_med_mail_id)
1053                 {
1054                   if (name_i.gibi_pile != PILE_SOUS_MAILLAGE) {
1055                     // error: must be PILE_SOUS_MAILLAGE
1056                   }
1057                   _iMed->_listGIBItoMED_mail.push_back(name_i);
1058                 }
1059               else if (itable == table_med_cham_id)
1060                 {
1061                   if (name_i.gibi_pile != PILE_FIELD &&
1062                       name_i.gibi_pile != PILE_NODES_FIELD)
1063                     {
1064                       // error: must be PILE_FIELD or PILE_NODES_FIELD
1065                     }
1066                   _iMed->_listGIBItoMED_cham.push_back(name_i);
1067                 }
1068               else if (itable == table_med_comp_id)
1069                 {
1070                   if (name_i.gibi_pile != PILE_STRINGS)
1071                     {
1072                       // error: gibi(short) names of components are kept in PILE_STRINGS
1073                     }
1074                   _iMed->_listGIBItoMED_comp.push_back(name_i);
1075                 }
1076             }
1077           else
1078             {
1079               // pass table
1080               for ( int ii = 0; ii < 4; ++ii ) next();
1081             }
1082         }
1083     } // for (int itable = 0; itable < nbObjects; itable++)
1084 }
1085
1086 //================================================================================
1087 /*!
1088  * \brief Read "PILE NUMERO  27"
1089  */
1090 //================================================================================
1091
1092 void SauvReader::read_PILE_STRINGS (const int                 nbObjects,
1093                                     std::vector<std::string>& objectNames,
1094                                     std::vector<int>&         nameIndices)
1095 {
1096   // IMP 0020434: mapping GIBI names to MED names
1097   initIntReading(2);
1098   int stringLen    = getIntNext();
1099   int nbSubStrings = getIntNext();
1100   if (nbSubStrings != nbObjects)
1101     THROW_IK_EXCEPTION("Error of reading PILE NUMERO  27" << lineNb() );
1102
1103   string aWholeString;
1104   if ( isXRD() )
1105     {
1106       const int fixedLength = 71;
1107       while ((int)aWholeString.length() < stringLen)
1108         {
1109           int remainLen = stringLen - aWholeString.length();
1110           int len;
1111           if ( remainLen > fixedLength )
1112             {
1113               len = fixedLength;
1114             }
1115           else
1116             {
1117               len = remainLen;
1118             }
1119           initNameReading(1, len);
1120           aWholeString += getName();
1121           next();
1122         }
1123     }
1124   else
1125     {
1126       char* line;
1127       const int fixedLength = 71;
1128       while ((int)aWholeString.length() < stringLen)
1129         {
1130           getNextLine( line );
1131           int remainLen = stringLen - aWholeString.length();
1132           if ( remainLen > fixedLength )
1133             {
1134               aWholeString += line + 1;
1135             }
1136           else
1137             {
1138               aWholeString += line + ( 72 - remainLen );
1139             }
1140         }
1141     }
1142   int prevOffset = 0;
1143   int currOffset = 0;
1144   initIntReading(nbSubStrings);
1145   for (int istr = 1; istr <= nbSubStrings; istr++, next())
1146     {
1147       currOffset = getInt();
1148       // fill mapStrings
1149       _iMed->_mapStrings[istr] = aWholeString.substr(prevOffset, currOffset - prevOffset);
1150       prevOffset = currOffset;
1151     }
1152 }