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 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::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 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               _iMed->_groups.push_back( SauvUtilities::Group() );
 692               group = & _iMed->_groups.back();
 693               group->_groups.swap( newGroups );
 694             }
 695         }
 696       // sort field sub-components and supports by cell type
 697       if ( group && !isSameCellType )
 698         {
 699           // sort groups
 700           vector<SauvUtilities::Group*>& groups = group->_groups;
 701           bool isModified = false, isSwapped = true;
 702           while ( isSwapped )
 703             {
 704               isSwapped = false;
 705               for ( size_t i = 1; i < groups.size(); ++i )
 706                 {
 707                   int nbN1 = groups[i-1]->empty() ? 0 : groups[i-1]->_cells[0]->_nodes.size();
 708                   int nbN2 = groups[i  ]->empty() ? 0 : groups[i  ]->_cells[0]->_nodes.size();
 709                   if ( nbN1 > nbN2 )
 710                     {
 711                       isSwapped = isModified = true;
 712                       std::swap( groups[i], groups[i-1] );
 713                     }
 714                 }
 715             }
 716           // relocate sub-components according to a new order of groups
 717           if ( isModified )
 718             {
 719               vector< DoubleField::_Sub_data > newSub( field->_sub.size() );
 720               vector< vector< double > >       newValues( field->_comp_values.size() );
 721               size_t iFromSub = 0, iNewSub = 0, iNewComp = 0;
 722               for ( ; iFromSub < field->_sub.size(); iFromSub += groups.size() )
 723                 {
 724                   size_t iFromComp = iNewComp;
 725                   for ( size_t iG = 0; iG < groups.size(); ++iG )
 726                     {
 727                       size_t iComp = iFromComp;
 728                       for ( size_t iSub = iFromSub; iSub < field->_sub.size(); ++iSub )
 729                         if ( field->_sub[ iSub ]._support == groups[ iG ] )
 730                           {
 731                             newSub[ iNewSub++ ] = field->_sub[ iSub ];
 732                             int iC = 0, nbC = field->_sub[ iSub ].nbComponents();
 733                             for ( ; iC < nbC; ++iC )
 734                               newValues[ iNewComp++ ].swap( field->_comp_values[ iComp++ ]);
 735                             break;
 736                           }
 737                         else
 738                           {
 739                             iComp += field->_sub[ iSub ].nbComponents();
 740                           }
 741                     }
 742                 }
 743               field->_sub.swap( newSub );
 744               field->_comp_values.swap( newValues );
 745             }
 746         }
 747     }
 748   if ( group )
 749     group->_isProfile = true;
 750
 751   field->_group = group;
 752 }
 753
 754 //================================================================================
 755 /*!
 756  * \brief Set field names
 757  */
 758 //================================================================================
 759
 760 void SauvReader::setFieldNames(const vector<SauvUtilities::DoubleField* >& fields,
 761                                const vector<string>&                       objets_nommes,
 762                                const vector<int>&                          indices_objets_nommes)
 763 {
 764   unsigned i;
 765   for ( i = 0; i < indices_objets_nommes.size(); ++i )
 766     {
 767       int fieldIndex = indices_objets_nommes[ i ];
 768       if ( fields[ fieldIndex - 1 ] )
 769         fields[ fieldIndex - 1 ]->_name = objets_nommes[ i ];
 770     }
 771 }
 772
 773 //================================================================================
 774 /*!
 775  * \brief Read "PILE NUMERO   2": NODE FIELDS
 776  */
 777 //================================================================================
 778
 779 void SauvReader::read_PILE_NODES_FIELD (const int                 nbObjects,
 780                                         std::vector<std::string>& objectNames,
 781                                         std::vector<int>&         nameIndices)
 782 {
 783   _iMed->_nodeFields.resize( nbObjects, (SauvUtilities::DoubleField*) 0 );
 784   for (int object=0; object!=nbObjects; ++object) // loop on fields
 785     {
 786       // EXAMPLE ( with no values )
 787
 788       // (1)       4       7       2       1
 789       // (2)     -88       0       3     -89       0       1     -90       0       2     -91
 790       // (2)       0       1
 791       // (3) FX   FY   FZ   FZ   FX   FY   FLX
 792       // (4)       0       0       0       0       0       0       0
 793       // (5)           cree  par  muc pri
 794       // (6)
 795       // (7)       2
 796
 797       // (1): nb subcomponents, nb components(total), IFOUR, nb attributes
 798       int nb_sub, total_nb_comp, nb_attr;
 799       int i_sub, i_comp;
 800       initIntReading( 4 );
 801       nb_sub        = getIntNext();
 802       total_nb_comp = getIntNext();
 803       next(); // ignore IFOUR
 804       nb_attr       = getIntNext();
 805       if ( nb_sub < 0 || total_nb_comp < 0 || nb_attr < 0 )
 806         THROW_IK_EXCEPTION("Error of field reading " << lineNb());
 807
 808       // (2) loop on subcomponents of a field, for each read
 809       // (a) support, (b) number of values and (c) number of components
 810       vector<Group*> supports( nb_sub );
 811       vector<int> nb_values  ( nb_sub );
 812       vector<int> nb_comps   ( nb_sub );
 813       int total_nb_values = 0;
 814       initIntReading( nb_sub * 3 );
 815       for ( i_sub = 0; i_sub < nb_sub; ++i_sub )
 816         {
 817           int supId = -getIntNext(); // (a) reference to support
 818           if ( supId < 1 || supId > (int)_iMed->_groups.size() )
 819             THROW_IK_EXCEPTION("Wrong mesh reference: "<< supId << lineNb() );
 820           supports[ i_sub ] = &_iMed->_groups[ supId-1 ]; // (a) reference to support
 821
 822           nb_values[ i_sub ] = getIntNext();    // (b) nb points
 823           total_nb_values += nb_values[ i_sub ];
 824           if ( nb_values[ i_sub ] < 0 )
 825             THROW_IK_EXCEPTION(" Wrong nb of points: " << nb_values[ i_sub ]  << lineNb() );
 826           nb_comps[ i_sub ] = getInt(); next();     // (c) nb of components in i_sub
 827         }
 828
 829       // create a field if there are values
 830       SauvUtilities::DoubleField* fdouble = 0;
 831       if ( total_nb_values > 0 )
 832         fdouble = new DoubleField( nb_sub, total_nb_comp );
 833       _iMed->_nodeFields[ object ] = fdouble;
 834
 835       // (3) component names
 836       initNameReading( total_nb_comp, 4 );
 837       for ( i_sub = 0; i_sub < nb_sub; ++i_sub )
 838         {
 839           // store support id and nb components of a sub
 840           if ( fdouble )
 841             fdouble->_sub[ i_sub ].setData( nb_comps[ i_sub ], supports[ i_sub ] );
 842           for ( i_comp = 0; i_comp < nb_comps[ i_sub ]; ++i_comp, next() )
 843             {
 844               // store component name
 845               string compName = getName();
 846               if ( fdouble )
 847                 fdouble->_sub[ i_sub ].compName( i_comp ) = compName;
 848             }
 849         }
 850       // (4) nb harmonics ( ignored )
 851       for ( initIntReading( total_nb_comp ); more(); next() );
 852       // (5) TYPE ( ignored )
 853       for (initNameReading(1, /*length=*/71); more(); next());
 854       // (6) TITRE ( ignored )
 855       for (initNameReading(1, /*length=*/71); more(); next());
 856       // (7) attributes ( ignored )
 857       for ( initIntReading( nb_attr ); more(); next() );
 858
 859       for ( i_sub = 0; i_sub < nb_sub; ++i_sub )
 860         {
 861           // loop on components: read values
 862           initDoubleReading( nb_values[ i_sub ] * nb_comps[ i_sub ] );
 863           for ( i_comp = 0; i_comp < nb_comps[ i_sub ]; ++i_comp )
 864             {
 865               if ( fdouble )
 866                 {
 867                   vector<double>& vals = fdouble->addComponent( nb_values[ i_sub ] );
 868                   for ( int i = 0; more() && i < nb_values[ i_sub ]; next(), ++i )
 869                     vals[ i ] = getDouble();
 870                 }
 871               else
 872                 {
 873                   for ( int i = 0; i < nb_values[ i_sub ]; next(), ++i );
 874                 }
 875             }
 876         } // loop on subcomponents of a field
 877
 878       // set a supporting group including all subs supports but only
 879       // if all subs have the same components
 880       if ( fdouble && fdouble->hasSameComponentsBySupport() )
 881         setFieldSupport( supports, fdouble );
 882       else
 883         for ( i_sub = 0; i_sub < nb_sub; ++i_sub )
 884           fdouble->_sub[ i_sub ]._support->_isProfile = true;
 885
 886     } // end loop on field objects
 887
 888   // set field names
 889   setFieldNames( _iMed->_nodeFields, objectNames, nameIndices );
 890
 891 }  // read_PILE_NODES_FIELD()
 892
 893 //================================================================================
 894 /*!
 895  * \brief Read "PILE NUMERO  39": FIELDS
 896  */
 897 //================================================================================
 898
 899 void SauvReader::read_PILE_FIELD (const int                 nbObjects,
 900                                   std::vector<std::string>& objectNames,
 901                                   std::vector<int>&         nameIndices)
 902 {
 903   // REAL EXAMPLE
 904
 905   // (1)        1       2       6      16
 906   // (2)                                                         CARACTERISTIQUES
 907   // (3)      -15  317773       4       0       0       0      -2       0       3
 908   // (4)             317581
 909   // (5)  0
 910   // (6)   317767  317761  317755  317815
 911   // (7)  YOUN     NU       H        SIGY
 912   // (8)  REAL*8            REAL*8            REAL*8            REAL*8
 913   // (9)        1       1       0       0
 914   // (10)  2.00000000000000E+05
 915   // (11)       1       1       0       0
 916   // (12)  3.30000000000000E-01
 917   // (13)       1       1       0       0
 918   // (14)  1.00000000000000E+04
 919   // (15)       6     706       0       0
 920   // (16)  1.00000000000000E+02  1.00000000000000E+02  1.00000000000000E+02
 921   // (17)  1.00000000000000E+02  1.00000000000000E+02  1.00000000000000E+02
 922   // (18)  ...
 923
 924   _iMed->_cellFields.resize( nbObjects, (SauvUtilities::DoubleField*) 0 );
 925   for (int object=0; object!=nbObjects; ++object) // pour chaque field
 926     {
 927       initIntReading( 4 );
 928       int i_sub, nb_sub = getIntNext(); // (1) <nb_sub> 2 6 <title length>
 929       next(); // skip "2"
 930       next(); // skip "6"
 931       int title_length = getIntNext(); // <title length>
 932       if ( nb_sub < 1 )
 933         THROW_IK_EXCEPTION("Error of field reading: wrong nb of subcomponents " << nb_sub << lineNb() );
 934
 935       string description;
 936       if ( title_length )
 937         {
 938           if ( isXRD() )
 939             {
 940               initNameReading(1, title_length);
 941               description = getName();
 942               next();
 943             }
 944           else
 945             {
 946               char* line;
 947               getNextLine( line ); // (2) title
 948               const int len = 72; // line length
 949               description = string(line + len - title_length, title_length); // title is right justified
 950             }
 951         }
 952       // look for a line starting with '-' : <reference to support>
 953       if ( isXRD() )
 954         {
 955           initIntReading( nb_sub * 9 );
 956         }
 957       else
 958         {
 959           do {
 960             initIntReading( nb_sub * 9 );
 961           } while ( getInt() >= 0 );
 962         }
 963       int total_nb_comp = 0;
 964       vector<Group*> supports( nb_sub );
 965       vector<int>     nb_comp( nb_sub );
 966       for ( i_sub = 0; i_sub < nb_sub; ++i_sub )
 967         {                                    // (3)
 968           int supportId     = -getIntNext(); // <reference to support>
 969           next();                            // ignore <address>
 970           nb_comp [ i_sub ] =  getIntNext(); // <nb of components in the sub>
 971           for ( int i = 0; i < 6; ++i ) next();  // ignore 6 ints, in example "0 0 0 -2 0 3"
 972
 973           if ( supportId < 1 || supportId > (int)_iMed->_groups.size() )
 974             THROW_IK_EXCEPTION("Error of field reading: wrong mesh reference "<< supportId << lineNb() );
 975           if ( nb_comp[ i_sub ] < 0 )
 976             THROW_IK_EXCEPTION("Error of field reading: wrong nb of components " <<nb_comp[ i_sub ] << lineNb() );
 977
 978           supports[ i_sub ] = &_iMed->_groups[ supportId-1 ];
 979           total_nb_comp += nb_comp[ i_sub ];
 980         }
 981       // (4) dummy strings
 982       for ( initNameReading( nb_sub, 17 ); more(); next() );
 983       // (5) dummy strings
 984       for ( initNameReading( nb_sub ); more(); next() );
 985
 986       // loop on subcomponents of a field, each of which refers to
 987       // a certain support and has its own number of components;
 988       // read component values
 989       SauvUtilities::DoubleField* fdouble = 0;
 990       for ( i_sub = 0; i_sub < nb_sub; ++ i_sub )
 991         {
 992           vector<string> comp_names( nb_comp[ i_sub ]), comp_type( nb_comp[ i_sub ]);
 993           // (6) nb_comp addresses of MELVAL structure
 994           for ( initIntReading( nb_comp[ i_sub ] ); more(); next() );
 995           // (7) component names
 996           for ( initNameReading( nb_comp[ i_sub ] ); more(); next() )
 997             comp_names[ index() ] = getName();
 998           // (8) component type
 999           for ( initNameReading( nb_comp[ i_sub ], 17 ); more(); next() ) // 17 is name width
1000             {
1001               comp_type[ index() ] = getName();
1002               // component types must be the same
1003               if ( index() > 0 && comp_type[ index() ] != comp_type[ index() - 1] )
1004                 THROW_IK_EXCEPTION( "Error of field reading: diff component types <"
1005                                     << comp_type[ index() ] << "> != <" << comp_type[ index() - 1 ]
1006                                     << ">" << lineNb() );
1007             }
1008           // now type is known, create a field, one for all subs
1009           bool isReal = (nb_comp[i_sub] > 0) ? (comp_type[0] == "REAL*8") : true;
1010           if ( !fdouble && total_nb_comp )
1011             {
1012               if ( !isReal )
1013                 cout << "Warning: read NOT REAL field, type <" << comp_type[0] << ">" << lineNb() << endl;
1014               _iMed->_cellFields[ object ] = fdouble = new SauvUtilities::DoubleField( nb_sub, total_nb_comp );
1015               fdouble->_description = description;
1016             }
1017           // store support id and nb components of a sub
1018           if ( fdouble )
1019             fdouble->_sub[ i_sub ].setData( nb_comp[ i_sub ], supports[ i_sub ]);
1020           // loop on components: read values
1021           for ( int i_comp = 0; i_comp < nb_comp[ i_sub ]; ++i_comp )
1022             {
1023               // (9) nb of values
1024               initIntReading( 4 );
1025               int nb_val_by_elem = getIntNext();
1026               int nb_values      = getIntNext();
1027               next();
1028               next();
1029               fdouble->_sub[ i_sub ]._nb_gauss[ i_comp ] = nb_val_by_elem;
1030
1031               // (10) values
1032               nb_values *= nb_val_by_elem;
1033               if ( fdouble )
1034                 {
1035                   vector<double> & vals = fdouble->addComponent( nb_values );
1036                   for ( isReal ? initDoubleReading( nb_values ) : initIntReading( nb_values ); more(); next())
1037                     vals[ index() ] = getDouble();
1038                   // store component name
1039                   fdouble->_sub[ i_sub ].compName( i_comp ) = comp_names[ i_comp ];
1040                 }
1041               else
1042                 {
1043                   for ( isReal ? initDoubleReading( nb_values ) : initIntReading( nb_values ); more(); next() ) ;
1044                 }
1045             }
1046         } // loop on subcomponents of a field
1047
1048       // set id of a group including all sub supports but only
1049       // if all subs have the same nb of components
1050       if ( fdouble && fdouble->hasSameComponentsBySupport() )
1051         setFieldSupport( supports, fdouble );
1052       else
1053         for ( i_sub = 0; i_sub < nb_sub; ++i_sub )
1054           fdouble->_sub[ i_sub ]._support->_isProfile = true;
1055
1056     } // end loop on field objects
1057
1058   // set field names
1059   setFieldNames( _iMed->_cellFields, objectNames, nameIndices );
1060
1061 } // read_PILE_FIELD()
1062
1063 //================================================================================
1064 /*!
1065  * \brief Read "PILE NUMERO  10": TABLES
1066  */
1067 //================================================================================
1068
1069 void SauvReader::read_PILE_TABLES (const int                 nbObjects,
1070                                    std::vector<std::string>& objectNames,
1071                                    std::vector<int>&         nameIndices)
1072 {
1073   // IMP 0020434: mapping GIBI names to MED names
1074
1075   string table_med_mail = "MED_MAIL";
1076   string table_med_cham = "MED_CHAM";
1077   string table_med_comp = "MED_COMP";
1078   int table_med_mail_id = -1;
1079   int table_med_cham_id = -1;
1080   int table_med_comp_id = -1;
1081   for (size_t iname = 0; iname < objectNames.size(); iname++)
1082     if      (objectNames[iname] == table_med_mail) table_med_mail_id = nameIndices[iname];
1083     else if (objectNames[iname] == table_med_cham) table_med_cham_id = nameIndices[iname];
1084     else if (objectNames[iname] == table_med_comp) table_med_comp_id = nameIndices[iname];
1085
1086   if ( isASCII() )
1087     if (table_med_mail_id < 0 && table_med_cham_id < 0 && table_med_comp_id < 0)
1088       return;
1089
1090   for (int itable = 1; itable <= nbObjects; itable++)
1091     {
1092       // read tables "MED_MAIL", "MED_CHAM" and "MED_COMP", that keeps correspondence
1093       // between GIBI names (8 symbols if any) and MED names (possibly longer)
1094       initIntReading(1);
1095       int nb_table_vals = getIntNext();
1096       if (nb_table_vals < 0)
1097         THROW_IK_EXCEPTION("Error of reading PILE NUMERO  10" << lineNb() );
1098
1099       int name_i_med_pile;
1100       initIntReading(nb_table_vals);
1101       for (int i = 0; i < nb_table_vals/4; i++)
1102         {
1103           if (itable == table_med_mail_id ||
1104               itable == table_med_cham_id ||
1105               itable == table_med_comp_id)
1106             {
1107               nameGIBItoMED name_i;
1108               name_i_med_pile  = getIntNext();
1109               name_i.med_id    = getIntNext();
1110               name_i.gibi_pile = getIntNext();
1111               name_i.gibi_id   = getIntNext();
1112
1113               if (name_i_med_pile != PILE_STRINGS)
1114                 {
1115                   // error: med(long) names are always kept in PILE_STRINGS
1116                 }
1117               if (itable == table_med_mail_id)
1118                 {
1119                   if (name_i.gibi_pile != PILE_SOUS_MAILLAGE) {
1120                     // error: must be PILE_SOUS_MAILLAGE
1121                   }
1122                   _iMed->_listGIBItoMED_mail.push_back(name_i);
1123                 }
1124               else if (itable == table_med_cham_id)
1125                 {
1126                   if (name_i.gibi_pile != PILE_FIELD &&
1127                       name_i.gibi_pile != PILE_NODES_FIELD)
1128                     {
1129                       // error: must be PILE_FIELD or PILE_NODES_FIELD
1130                     }
1131                   _iMed->_listGIBItoMED_cham.push_back(name_i);
1132                 }
1133               else if (itable == table_med_comp_id)
1134                 {
1135                   if (name_i.gibi_pile != PILE_STRINGS)
1136                     {
1137                       // error: gibi(short) names of components are kept in PILE_STRINGS
1138                     }
1139                   _iMed->_listGIBItoMED_comp.push_back(name_i);
1140                 }
1141             }
1142           else
1143             {
1144               // pass table
1145               for ( int ii = 0; ii < 4; ++ii ) next();
1146             }
1147         }
1148     } // for (int itable = 0; itable < nbObjects; itable++)
1149 }
1150
1151 //================================================================================
1152 /*!
1153  * \brief Read "PILE NUMERO  27"
1154  */
1155 //================================================================================
1156
1157 void SauvReader::read_PILE_STRINGS (const int                 nbObjects,
1158                                     std::vector<std::string>& objectNames,
1159                                     std::vector<int>&         nameIndices)
1160 {
1161   // IMP 0020434: mapping GIBI names to MED names
1162   initIntReading(2);
1163   int stringLen    = getIntNext();
1164   int nbSubStrings = getIntNext();
1165   if (nbSubStrings != nbObjects)
1166     THROW_IK_EXCEPTION("Error of reading PILE NUMERO  27" << lineNb() );
1167
1168   string aWholeString;
1169   if ( isXRD() )
1170     {
1171       const int fixedLength = 71;
1172       while ((int)aWholeString.length() < stringLen)
1173         {
1174           int remainLen = stringLen - aWholeString.length();
1175           int len;
1176           if ( remainLen > fixedLength )
1177             {
1178               len = fixedLength;
1179             }
1180           else
1181             {
1182               len = remainLen;
1183             }
1184           initNameReading(1, len);
1185           aWholeString += getName();
1186           next();
1187         }
1188     }
1189   else
1190     {
1191       char* line;
1192       const int fixedLength = 71;
1193       while ((int)aWholeString.length() < stringLen)
1194         {
1195           getNextLine( line );
1196           int remainLen = stringLen - aWholeString.length();
1197           if ( remainLen > fixedLength )
1198             {
1199               aWholeString += line + 1;
1200             }
1201           else
1202             {
1203               aWholeString += line + ( 72 - remainLen );
1204             }
1205         }
1206     }
1207   int prevOffset = 0;
1208   int currOffset = 0;
1209   initIntReading(nbSubStrings);
1210   for (int istr = 1; istr <= nbSubStrings; istr++, next())
1211     {
1212       currOffset = getInt();
1213       // fill mapStrings
1214       _iMed->_mapStrings[istr] = aWholeString.substr(prevOffset, currOffset - prevOffset);
1215       prevOffset = currOffset;
1216     }
1217 }