src/MEDLoader/SauvReader.cxx

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