src/MEDMEM/MEDMEM_EnsightMeshDriver.cxx

   1 // Copyright (C) 2007-2013  CEA/DEN, EDF R&D, OPEN CASCADE
   2 //
   3 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   5 //
   6 // This library is free software; you can redistribute it and/or
   7 // modify it under the terms of the GNU Lesser General Public
   8 // License as published by the Free Software Foundation; either
   9 // version 2.1 of the License.
  10 //
  11 // This library is distributed in the hope that it will be useful,
  12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 // Lesser General Public License for more details.
  15 //
  16 // You should have received a copy of the GNU Lesser General Public
  17 // License along with this library; if not, write to the Free Software
  18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  19 //
  20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  21 //
  22
  23 #include "MEDMEM_EnsightMeshDriver.hxx"
  24
  25 #include <sstream>
  26 #include <iomanip>
  27 #include <numeric>
  28
  29 #include "MEDMEM_define.hxx"
  30 #include "MEDMEM_Family.hxx"
  31 #include "MEDMEM_Group.hxx"
  32 #include "MEDMEM_Coordinate.hxx"
  33 #include "MEDMEM_Connectivity.hxx"
  34 #include "MEDMEM_Mesh.hxx"
  35 #include "MEDMEM_CellModel.hxx"
  36 #include "MEDMEM_Grid.hxx"
  37
  38 #include "MEDMEM_MedMeshDriver.hxx"
  39
  40 using namespace std;
  41 using namespace MEDMEM;
  42 using namespace MED_EN;
  43 using namespace MEDMEM_ENSIGHT;
  44
  45 #define TStrTool _ASCIIFileReader
  46 #define TOLERANCE 1e-15
  47
  48 //#define ELEMENT_ID_GIVEN
  49
  50 namespace {
  51
  52   // ---------------------------------------------------------------------
  53   /*!
  54    * \brief The beginning of mesh description used to distinguish files
  55    * generated by ENSIGHT_MESH_WRONLY_DRIVER from others
  56    */
  57   const char* theDescriptionPrefix = "Meshing from MedMemory. ";
  58
  59   // ---------------------------------------------------------------------
  60   /*!
  61    * \brief Default name of a mesh read from EnSight
  62    */
  63   const char* theDefaultMeshName = "EnsightMesh";
  64
  65   // ---------------------------------------------------------------------
  66   /*!
  67    * \brief Max number of types in EnSight part
  68    */
  69   const int   theMaxNbTypes = 20;
  70
  71   // ---------------------------------------------------------------------
  72   /*!
  73    * \brief Make array with elements == index[i+1]-index[i]
  74    *  \param size - result array size
  75    */
  76   int* getNumbersByIndex( const int* index, int size, const int* elemNumbers=0)
  77   {
  78     int* numbers = new int[size];
  79     int* n = numbers;
  80     if ( elemNumbers ) {
  81       const int *elem = elemNumbers-1, *elemEnd = elemNumbers + size;
  82       while ( ++elem < elemEnd )
  83         *n++ = index[*elem] - index[*elem-1];
  84     }
  85     else {
  86       const int *ind = index, *indEnd = index + size + 1;
  87       while ( ++ind < indEnd )
  88         *n++ = ind[0] - ind[-1];
  89     }
  90     return numbers;
  91   }
  92   // ---------------------------------------------------------------------
  93   /*!
  94    * \brief Type used to delete numbers returned by getNumbersByIndex()
  95    */
  96   typedef _ValueOwner<int> TNumbers;
  97
  98 } // namespace
  99
 100
 101 ENSIGHT_MESH_DRIVER::ENSIGHT_MESH_DRIVER(): _CaseFileDriver_User(), _ptrMesh((MESH *)NULL)
 102 {
 103 }
 104
 105 ENSIGHT_MESH_DRIVER::ENSIGHT_MESH_DRIVER(const string & fileName,
 106                                          GMESH *        ptrMesh)
 107   :_CaseFileDriver_User(fileName,MED_EN::RDWR), _ptrMesh(ptrMesh)
 108 {
 109   if ( !_ptrMesh )
 110     throw MEDEXCEPTION("ENSIGHT_MESH_DRIVER(fileName, ptrMesh) : mesh is NULL");
 111   _meshName = _ptrMesh->getName();
 112 }
 113
 114 ENSIGHT_MESH_DRIVER::ENSIGHT_MESH_DRIVER(const string & fileName,
 115                                          GMESH *        ptrMesh,
 116                                          med_mode_acces accessMode)
 117   :_CaseFileDriver_User(fileName,accessMode), _ptrMesh(ptrMesh)
 118 {
 119   if ( !_ptrMesh )
 120     throw MEDEXCEPTION("ENSIGHT_MESH_DRIVER(fileName, ptrMesh) : mesh is NULL");
 121   _meshName = _ptrMesh->getName();
 122 }
 123
 124 ENSIGHT_MESH_DRIVER::ENSIGHT_MESH_DRIVER(const ENSIGHT_MESH_DRIVER & driver)
 125   :_CaseFileDriver_User(driver), _ptrMesh(driver._ptrMesh), _meshName(driver._meshName)
 126 {
 127 }
 128
 129 ENSIGHT_MESH_DRIVER::~ENSIGHT_MESH_DRIVER()
 130 {
 131   MESSAGE_MED("ENSIGHT_MESH_DRIVER::~ENSIGHT_MESH_DRIVER() has been destroyed");
 132 }
 133
 134 void    ENSIGHT_MESH_DRIVER::setMeshName(const string & meshName) { _meshName = meshName; }
 135
 136 string  ENSIGHT_MESH_DRIVER::getMeshName() const { return _meshName; }
 137
 138 void ENSIGHT_MESH_DRIVER::openConst(bool checkDataFile) const
 139 {
 140   const char * LOC ="ENSIGHT_MESH_DRIVER::open() : ";
 141   BEGIN_OF_MED(LOC);
 142
 143   if ( checkDataFile )
 144   {
 145     if ( getDataFileName().empty() )
 146       throw MED_EXCEPTION
 147         ( LOCALIZED( STRING(LOC) << "Internal error, geometry file name is empty"));
 148
 149     if (!canOpenFile( getDataFileName(), getAccessMode() ))
 150       throw MED_EXCEPTION
 151         ( LOCALIZED( STRING(LOC) << "Can not open Ensight Geometry file " << getDataFileName()
 152                      << " in access mode " << getAccessMode()));
 153   }
 154   else
 155   {
 156     if ( getCaseFileName().empty() )
 157       throw MED_EXCEPTION
 158         ( LOCALIZED( STRING(LOC) << "Case file name is empty, "
 159                      "please set a correct fileName before calling open()"));
 160
 161     if ( !canOpenFile( getCaseFileName(), getAccessMode() ))
 162       throw MED_EXCEPTION
 163         ( LOCALIZED( STRING(LOC) << "Can not open Ensight Case file " << getCaseFileName()
 164                      << " in access mode " << getAccessMode()));
 165   }
 166
 167   END_OF_MED(LOC);
 168 }
 169
 170 void ENSIGHT_MESH_DRIVER::open() {
 171   openConst() ;
 172 }
 173
 174 void ENSIGHT_MESH_DRIVER::close() {
 175 }
 176
 177 // ================================================================================
 178 // WRONLY
 179 // ================================================================================
 180
 181 ENSIGHT_MESH_WRONLY_DRIVER::ENSIGHT_MESH_WRONLY_DRIVER() :
 182   ENSIGHT_MESH_DRIVER(), _append(0)
 183 {
 184 }
 185
 186 ENSIGHT_MESH_WRONLY_DRIVER::ENSIGHT_MESH_WRONLY_DRIVER(const string & fileName,
 187                                                        const GMESH *  ptrMesh,
 188                                                        bool           append)
 189   : ENSIGHT_MESH_DRIVER( fileName, (GMESH*)ptrMesh, WRONLY ), _append(append)
 190 {
 191 }
 192
 193 ENSIGHT_MESH_WRONLY_DRIVER::ENSIGHT_MESH_WRONLY_DRIVER(const ENSIGHT_MESH_WRONLY_DRIVER & driver)
 194   : ENSIGHT_MESH_DRIVER(driver),_append(driver._append)
 195 {
 196 }
 197
 198 ENSIGHT_MESH_WRONLY_DRIVER::~ENSIGHT_MESH_WRONLY_DRIVER()
 199 {
 200 }
 201
 202 GENDRIVER * ENSIGHT_MESH_WRONLY_DRIVER::copy() const
 203 {
 204   return new ENSIGHT_MESH_WRONLY_DRIVER(*this) ;
 205 }
 206
 207 void ENSIGHT_MESH_WRONLY_DRIVER::read() throw (MEDEXCEPTION) {
 208   throw MEDEXCEPTION("ENSIGHT_MESH_WRONLY_DRIVER::read : Can't read with a WRONLY driver !");
 209 }
 210
 211 //================================================================================
 212 /*!
 213  * \brief writing
 214  */
 215 //================================================================================
 216
 217 void ENSIGHT_MESH_WRONLY_DRIVER::write() const throw (MEDEXCEPTION)
 218 {
 219   const char * LOC = "ENSIGHT_MESH_WRONLY_DRIVER::write() : ";
 220   BEGIN_OF_MED(LOC);
 221
 222   openConst(false) ; // check if can write to case file
 223
 224   // Ensight case organization requires a main file (filename.case) which defines organization
 225
 226   _CaseFileDriver caseFile( getCaseFileName(), this);
 227   if ( _append )
 228     caseFile.read();
 229   caseFile.addMesh( this );
 230   caseFile.write(); // all fields of _CaseFileDriver_User are set by this method
 231
 232   openConst(true) ; // check if can write to data file
 233
 234   cout << "-> creating the Ensight geometry file " << getDataFileName() << endl ;
 235
 236   // Store mesh description and a special mark in the first two description lines each
 237   // of 79 chars length maximum, while MED mesh description is up to 200 chars
 238   const char* line1 = theDescriptionPrefix;
 239   string      line2 = _ptrMesh->getDescription();
 240   for ( unsigned i = 0; i < line2.size(); ++i ) { // protect from gabage
 241     if ( !line2[ i ] || !isascii( line2[ i ])) {
 242       line2.resize( i );
 243       break;
 244     }
 245   }
 246   if ((int) line2.size() >= MAX_LINE_LENGTH )
 247     line2.resize( MAX_LINE_LENGTH );
 248
 249   // EnSight will assign node/element visible numbers it-self
 250   const char* line3 = "node id assign";
 251 #ifdef ELEMENT_ID_GIVEN
 252   const char* line4 = "element id given";
 253 #else
 254   const char* line4 = "element id assign";
 255 #endif
 256
 257   MESH* mesh = const_cast<MESH*>( _ptrMesh->convertInMESH() ) ; // we write unstructured only
 258
 259   if ( isBinaryEnSightFormatForWriting() )
 260   {
 261     // ======================================================
 262     //                          Binary
 263     // ======================================================
 264
 265     _BinaryFileWriter ensightGeomFile( getDataFileName() );
 266
 267     ensightGeomFile.addString("C Binary");
 268     ensightGeomFile.addString(line1);
 269     ensightGeomFile.addString(line2);
 270     ensightGeomFile.addString(line3);
 271     ensightGeomFile.addString(line4);
 272
 273     // function to write a support as a part
 274     typedef void (ENSIGHT_MESH_WRONLY_DRIVER::* TWritePart) (_BinaryFileWriter&, const SUPPORT*) const;
 275     TWritePart writePart;
 276     if ( isGoldFormat() )
 277     {
 278       // GOLD
 279       writePart = & ENSIGHT_MESH_WRONLY_DRIVER::writePartGoldBinary;
 280     }
 281     else
 282     {
 283       // ENSIGHT 6. Write addionally global nodes
 284       writePart = & ENSIGHT_MESH_WRONLY_DRIVER::writePart6Binary;
 285
 286       // All point are in 3D, so if we are in 1D or 2D, we complete by zero !
 287       int SpaceDimension = mesh->getSpaceDimension() ;
 288       int NumberOfNodes  = mesh->getNumberOfNodes() ;
 289       ensightGeomFile.addString("coordinates");
 290       ensightGeomFile.addInt( NumberOfNodes );
 291       const double *coordinate = mesh->getCoordinates(MED_FULL_INTERLACE) ;
 292       if ( SpaceDimension == 3 ) {
 293         ensightGeomFile.addReal(coordinate, NumberOfNodes * SpaceDimension );
 294       }
 295       else {
 296         typedef _ValueIterator< double > TComponentIt;
 297         vector< TComponentIt > coordCompIt( 3 );
 298         for (int j=0; j<3; j++, coordinate++)
 299           if ( j < SpaceDimension )
 300             coordCompIt[ j ] = TComponentIt( coordinate, SpaceDimension );
 301           else
 302             coordCompIt[ j ] = TComponentIt(); // to iterate on zeros
 303         ensightGeomFile.addReal( coordCompIt, NumberOfNodes, MED_FULL_INTERLACE );
 304       }
 305     }
 306
 307     // We put connectivity
 308
 309     if ( isToWriteEntity( MED_CELL, mesh ))
 310     {
 311       SUPPORT *allCells = const_cast<SUPPORT*>( _ptrMesh->getSupportOnAll( MED_CELL ));
 312       string oldName = allCells->getName();
 313       allCells->setName( _ptrMesh->getName() );
 314       try
 315       {
 316         (this->*writePart)( ensightGeomFile, allCells );
 317       }
 318       catch (MED_EXCEPTION& ex)
 319       {
 320         allCells->setName( oldName );
 321         throw ex;
 322       }
 323       allCells->setName( oldName );
 324     }
 325     // And meshdim-1 connectivity
 326     if ( isToWriteEntity( MED_FACE, mesh ))
 327     {
 328       const SUPPORT *allFaces = _ptrMesh->getSupportOnAll( MED_FACE );
 329       (this->*writePart)( ensightGeomFile, allFaces);
 330     }
 331     else if ( isToWriteEntity(MED_EDGE, mesh))
 332     {
 333       const SUPPORT *allEdges = _ptrMesh->getSupportOnAll( MED_EDGE );
 334       (this->*writePart)( ensightGeomFile, allEdges);
 335     }
 336
 337     // Write all groups as parts
 338
 339     for ( int ent = MED_CELL; ent < MED_ALL_ENTITIES; ++ent )
 340     {
 341       medEntityMesh entity = (medEntityMesh) ent;
 342       int nbGroups = mesh->getNumberOfGroups(entity);
 343       for ( int i=1; i<=nbGroups; i++)
 344       {
 345         const GROUP* group = mesh->getGroup( entity, i );
 346         (this->*writePart)( ensightGeomFile, group );
 347       }
 348     }
 349
 350   }
 351   else
 352   {
 353     // ======================================================
 354     //                           ASCII
 355     // ======================================================
 356     ofstream ensightGeomFile( getDataFileName().c_str(), ios::out);
 357     ensightGeomFile.setf(ios::scientific);
 358     ensightGeomFile.precision(5);
 359
 360     ensightGeomFile << line1 << endl
 361                     << line2 << endl
 362                     << line3 << endl
 363                     << line4 << endl;
 364
 365     // function to write a support as a part
 366     typedef void (ENSIGHT_MESH_WRONLY_DRIVER::* TWritePart) (ofstream&, const SUPPORT*) const;
 367     TWritePart writePart;
 368     if ( isGoldFormat() )
 369     {
 370       // GOLD
 371       writePart = & ENSIGHT_MESH_WRONLY_DRIVER::writePartGoldASCII;
 372     }
 373     else
 374     {
 375       // ENSIGHT 6. Write addionally global nodes
 376       writePart = & ENSIGHT_MESH_WRONLY_DRIVER::writePart6ASCII;
 377
 378       // Put points (all point are in 3D, so if we are in 1D or 2D, we complete by zero !
 379       int SpaceDimension = mesh->getSpaceDimension() ;
 380       int NumberOfNodes  = mesh->getNumberOfNodes() ;
 381       string zeros;
 382       if (SpaceDimension==2) zeros = " 0.00000e+00";
 383       if (SpaceDimension==1) zeros = " 0.00000e+00 0.00000e+00";
 384       ensightGeomFile << "coordinates" << endl
 385                       << setw(8) << NumberOfNodes << endl ;
 386       const double *coordinate = mesh->getCoordinates(MED_FULL_INTERLACE) ;
 387       for (int i=0; i<NumberOfNodes; i++)
 388       {
 389         //ensightGeomFile << setw(8) << i+1 ; // node id
 390         for (int j=0; j<SpaceDimension; j++, coordinate++)
 391           ensightGeomFile << setw(12) << *coordinate;
 392         ensightGeomFile << zeros << endl ;
 393       }
 394     }
 395
 396     // We put connectivity
 397
 398     if ( isToWriteEntity( MED_CELL, mesh ))
 399     {
 400       SUPPORT *allCells = const_cast<SUPPORT*>( _ptrMesh->getSupportOnAll( MED_CELL ));
 401       string oldName = allCells->getName();
 402       allCells->setName( _ptrMesh->getName() );
 403       try
 404       {
 405         (this->*writePart)( ensightGeomFile, allCells );
 406       }
 407       catch (MED_EXCEPTION& ex)
 408       {
 409         allCells->setName( oldName );
 410         throw ex;
 411       }
 412       allCells->setName( oldName );
 413     }
 414     // And meshdim-1 connectivity
 415     if ( isToWriteEntity( MED_FACE, mesh ))
 416     {
 417       const SUPPORT *allFaces = _ptrMesh->getSupportOnAll( MED_FACE );
 418       (this->*writePart)( ensightGeomFile, allFaces);
 419     }
 420     else if ( isToWriteEntity(MED_EDGE, mesh))
 421     {
 422       const SUPPORT *allEdges = _ptrMesh->getSupportOnAll( MED_EDGE );
 423       (this->*writePart)( ensightGeomFile, allEdges);
 424     }
 425
 426     // Write all groups as parts
 427
 428     for ( int ent = MED_CELL; ent < MED_ALL_ENTITIES; ++ent )
 429     {
 430       medEntityMesh entity = (medEntityMesh) ent;
 431       int nbGroups = mesh->getNumberOfGroups(entity);
 432       for ( int i=1; i<=nbGroups; i++)
 433       {
 434         const GROUP* group = mesh->getGroup( entity, i );
 435         (this->*writePart)( ensightGeomFile, group );
 436       }
 437     }
 438
 439     ensightGeomFile.close();
 440
 441     mesh->removeReference();
 442
 443   } // end ASCII format
 444
 445 } // ENSIGHT_MESH_WRONLY_DRIVER::write()
 446
 447 //================================================================================
 448 /*!
 449  * \brief Write support as an EnSight Gold part
 450  */
 451 //================================================================================
 452
 453 void ENSIGHT_MESH_WRONLY_DRIVER::writePartGoldBinary(_BinaryFileWriter& ensightGeomFile,
 454                                                      const SUPPORT*     support) const
 455 {
 456   // part number
 457   int partNum = getPartNumber( support );
 458   if ( !partNum )
 459     throw MED_EXCEPTION ( LOCALIZED( STRING("Internal error: invalid part number")));
 460   ensightGeomFile.addString( "part" );
 461   ensightGeomFile.addInt( partNum );
 462
 463   // group/mesh name
 464   ensightGeomFile.addString( support->getName() );
 465
 466   // get geom types
 467   medEntityMesh entity = support->getEntity();
 468   int nbTypes = support->getNumberOfTypes();
 469   const medGeometryElement* geoType = support->getTypes();
 470
 471   const int * connectivity = 0;
 472   const int * elemConnectivity = 0;
 473   const int * index = 0;
 474   int j;
 475
 476   const MESH* mesh = support->getMesh()->convertInMESH() ; // we write unstructured only
 477
 478   // COORDINATES                                                             Gold binary
 479   // ===================================================================================
 480   // In Ensight, coordinates of nodes of support elements are in MED_NO_INTERLACE mode.
 481   // We are to write only nodes belonging to elements of the support and
 482   // nodal connectivity should refer to these nodes.
 483   map<int, int> med2ensIds;
 484   map<int, int>::iterator medEnsIt;
 485   int SpaceDimension = _ptrMesh->getSpaceDimension() ;
 486   int NumberOfNodes  = _ptrMesh->getNumberOfNodes() ;
 487   // -------------------------------------------------
 488   if ( support->isOnAllElements() )
 489   {
 490     // nb of nodes
 491     ensightGeomFile.addString( "coordinates" );
 492     ensightGeomFile.addInt( NumberOfNodes );
 493
 494     // coordinates
 495     const double *coordinate = mesh->getCoordinates(MED_FULL_INTERLACE);
 496     typedef _ValueIterator< double > TComponentIt;
 497     vector< TComponentIt > coordCompIt( 1 );
 498     for (int j=0; j<SPACE_DIM; j++, coordinate++) { // loop on dimensions
 499       if ( j < SpaceDimension )
 500         coordCompIt[ 0 ] = TComponentIt( coordinate, SpaceDimension );
 501       else
 502         coordCompIt[ 0 ] = TComponentIt(); // to iterate on zeros
 503       ensightGeomFile.addReal( coordCompIt, NumberOfNodes, MED_NO_INTERLACE );
 504     }
 505   }
 506   // -------------------------------------------------
 507   else // support is not on all elements
 508   {
 509     // nb of nodes
 510     getSupportNodes( support, med2ensIds );
 511     NumberOfNodes = med2ensIds.size();
 512     ensightGeomFile.addString( "coordinates" );
 513     ensightGeomFile.addInt( NumberOfNodes );
 514
 515     // coordinates
 516     vector<float> floatCoords( NumberOfNodes );
 517     for ( j=0; j < SPACE_DIM; j++) { // loop on dimensions
 518       medEnsIt = med2ensIds.begin();
 519       if ( j < SpaceDimension ) {
 520         const double *coordinate = mesh->getCoordinates(MED_FULL_INTERLACE) + j;
 521         for (int i=0; i<NumberOfNodes; i++, ++medEnsIt )
 522           floatCoords[ i ] = (float) coordinate[ (medEnsIt->first-1) * SpaceDimension];
 523       }
 524       else if ( j-1 < SpaceDimension ) {
 525         for (int i=0; i<NumberOfNodes; i++)
 526           floatCoords[ i ] = 0.;
 527       }
 528       ensightGeomFile.addReal( &floatCoords[0], NumberOfNodes );
 529     }
 530     // assign local node ids
 531     for ( medEnsIt = med2ensIds.begin(), j=1; j<=NumberOfNodes; j++, ++medEnsIt )
 532       medEnsIt->second = j;
 533   }
 534
 535   // CONNECTIVITY                                                            Gold binary
 536   // ===================================================================================
 537   // loop on types
 538   for (int i=0; i<nbTypes; i++)
 539   {
 540     const medGeometryElement    medType = geoType[i];
 541     const TEnSightElemType& ensightType = getEnSightType(medType);
 542     const int numberOfCell              = support->getNumberOfElements(medType);
 543     int nbCellNodes                     = ensightType._medIndex.size();
 544
 545     // type name and nb cells
 546     ensightGeomFile.addString( ensightType._name );
 547     ensightGeomFile.addInt(  numberOfCell );
 548
 549     vector<int> nodeIds;
 550
 551     // -------------------------------------------------
 552     if ( support->isOnAllElements() )
 553     {
 554 #ifdef ELEMENT_ID_GIVEN
 555       // elem numbers
 556       nodeIds.resize( numberOfCell );
 557       for ( j = 1; j <= numberOfCell; j++)
 558         nodeIds[ j-1 ] = j;
 559       ensightGeomFile.addInt( nodeIds );
 560
 561       if ( entity != MED_NODE ) nodeIds.clear();
 562 #endif
 563
 564       if ( nbCellNodes > 1 ) // STANDARD ELEMENTS connectivity
 565       {
 566         connectivity = mesh->getConnectivity( MED_NODAL, entity, medType);
 567         nodeIds.reserve( numberOfCell * nbCellNodes);
 568         for (j = 0 ; j < numberOfCell; j++, connectivity += nbCellNodes)
 569           for (int k=0; k<nbCellNodes; k++)
 570             nodeIds.push_back( connectivity[ ensightType._medIndex[k] ]);
 571         ensightGeomFile.addInt( nodeIds );
 572       }
 573       else if ( entity == MED_NODE ) // NODES connectivity
 574       {
 575 #if !defined(ELEMENT_ID_GIVEN)
 576         nodeIds.resize( numberOfCell );
 577         for ( j = 1; j <= numberOfCell; j++)
 578           nodeIds[ j-1 ] = j;
 579 #endif
 580         ensightGeomFile.addInt( nodeIds );
 581       }
 582       else if ( medType == MED_POLYGON ) // POLYGONs connectivity
 583       {
 584         int nbStdCells = mesh->getGlobalNumberingIndex(entity)[i]-1;
 585         connectivity   = mesh->getConnectivity( MED_NODAL, entity, medType);
 586         int connLength = mesh->getConnectivityLength( MED_NODAL, entity, medType);
 587         index          = mesh->getConnectivityIndex(MED_NODAL, entity);
 588
 589         // number of nodes in each element
 590         {
 591           TIntOwner nbNodesInPoly( getNumbersByIndex( index+nbStdCells, numberOfCell ));
 592           ensightGeomFile.addInt( nbNodesInPoly.myValues, numberOfCell );
 593         } // nbNodesInPoly is deleted here
 594
 595         // connectivity
 596         ensightGeomFile.addInt( connectivity, connLength );
 597       }
 598       else // POLYHEDRA connectivity
 599       {
 600         connectivity   = mesh->getConnectivity( MED_NODAL, entity, medType);
 601         int nbStdCells = mesh->getGlobalNumberingIndex(entity)[i]-1;
 602         index          = mesh->getConnectivityIndex(MED_NODAL, entity) + nbStdCells;
 603
 604         vector<int> nbFacesInPolyhedron, nbNodesInFace, faceConn;
 605         for ( int j = 0; j < numberOfCell; ++j )
 606         {
 607           int nbFaces = 0, nbNodes = 0;
 608           for ( int k = index[j]; k < index[j+1]; ++k )
 609             if ( connectivity[k-1] == -1 )
 610             {
 611               nbNodesInFace.push_back( nbNodes );
 612               nbNodes = 0;
 613               ++nbFaces;
 614             }
 615             else
 616             {
 617               faceConn.push_back( connectivity[k-1] );
 618               ++nbNodes;
 619             }
 620           nbNodesInFace.push_back( nbNodes );
 621           nbFacesInPolyhedron.push_back( nbFaces+1 );
 622         }
 623
 624         // nb of faces in each polyhedron
 625         ensightGeomFile.addInt( nbFacesInPolyhedron );
 626         // number of nodes in each face
 627         ensightGeomFile.addInt( nbNodesInFace );
 628         // connectivity
 629         ensightGeomFile.addInt( faceConn );
 630       }
 631     }
 632     // -------------------------------------------------
 633     else // support is not on all elements
 634     {
 635       const int *number = support->getNumber(medType);
 636
 637 #ifdef ELEMENT_ID_GIVEN
 638       ensightGeomFile.addInt( number, numberOfCell );
 639 #endif
 640       if ( nbCellNodes > 1 ) // STANDARD ELEMENTS connectivity
 641       {
 642         connectivity = mesh->getConnectivity( MED_NODAL, entity, MED_ALL_ELEMENTS);
 643         index = mesh->getConnectivityIndex(MED_NODAL, entity);
 644
 645         nodeIds.reserve( numberOfCell * nbCellNodes);
 646         for (j=0; j<numberOfCell; j++) {
 647           int elem = number[j];
 648           elemConnectivity = connectivity + index[elem-1]-1;
 649           for (int k=0; k<nbCellNodes; k++)
 650           {
 651             int node = elemConnectivity[ ensightType._medIndex[k] ];
 652             nodeIds.push_back( med2ensIds[ node ]);
 653           }
 654         }
 655         ensightGeomFile.addInt( nodeIds );
 656       }
 657       else if ( entity == MED_NODE )  // NODES connectivity
 658       {
 659         nodeIds.resize( numberOfCell );
 660         for ( j = 1; j <= numberOfCell; j++)
 661           nodeIds[ j-1 ] = j;
 662         ensightGeomFile.addInt( nodeIds );
 663       }
 664       else if ( medType == MED_POLYGON ) // POLYGONs connectivity
 665       {
 666         connectivity = mesh->getConnectivity( MED_NODAL, entity, MED_ALL_ELEMENTS);
 667         index = mesh->getConnectivityIndex(MED_NODAL, entity);
 668
 669         // number of nodes in each element
 670         {
 671           TIntOwner nbNodesInPoly( getNumbersByIndex( index, numberOfCell, number ));
 672           ensightGeomFile.addInt( nbNodesInPoly.myValues, numberOfCell );
 673         } // nbNodesInPoly is deleted here
 674
 675         // connectivity
 676         for ( j = 0; j < numberOfCell; ++j )
 677         {
 678           int elem = number[ j ];
 679           elemConnectivity   = connectivity + index[ elem-1 ]-1;
 680           const int* connEnd = connectivity + index[ elem   ]-1;
 681           while ( elemConnectivity < connEnd )
 682             nodeIds.push_back( med2ensIds[ *elemConnectivity++ ]);
 683         }
 684         ensightGeomFile.addInt( nodeIds );
 685       }
 686       else // POLYHEDRA connectivity
 687       {
 688         connectivity = mesh->getConnectivity( MED_NODAL, entity, MED_ALL_ELEMENTS);
 689         index = mesh->getConnectivityIndex(MED_NODAL, entity);
 690         vector<int> nbFacesInPolyhedron, nbNodesInFace, faceConn;
 691         for ( int j = 0; j < numberOfCell; ++j )
 692         {
 693           int elem    = number[ j ];
 694           int nbFaces = 0, nbNodes = 0;
 695           for ( int k = index[elem]; k < index[elem+1]; ++k )
 696             if ( connectivity[k-1] == -1 )
 697             {
 698               nbNodesInFace.push_back( nbNodes );
 699               nbNodes = 0;
 700               ++nbFaces;
 701             }
 702             else
 703             {
 704               faceConn.push_back( connectivity[k-1] );
 705               ++nbNodes;
 706             }
 707           nbNodesInFace.push_back( nbNodes );
 708           nbFacesInPolyhedron.push_back( nbFaces+1 );
 709         }
 710
 711         // nb of faces in each polyhedron
 712         ensightGeomFile.addInt( nbFacesInPolyhedron );
 713         // number of nodes in each face
 714         ensightGeomFile.addInt( nbNodesInFace );
 715         // connectivity
 716         ensightGeomFile.addInt( faceConn );
 717       }
 718     }
 719   }
 720
 721   mesh->removeReference();
 722
 723 } // writePartGoldBinary()
 724
 725 //================================================================================
 726 /*!
 727  * \brief Write support as an EnSight Gold part
 728  */
 729 //================================================================================
 730
 731 void ENSIGHT_MESH_WRONLY_DRIVER::writePartGoldASCII(ofstream&      ensightGeomFile,
 732                                                     const SUPPORT* support) const
 733 {
 734   const int iw = 10;
 735
 736   // part number
 737   int partNum = getPartNumber( support );
 738   ensightGeomFile << "part" << endl
 739                   << setw(iw) << partNum << endl;
 740   if ( !partNum )
 741     throw MED_EXCEPTION ( LOCALIZED( STRING("Internal error: invalid part number")));
 742
 743   // group/mesh name
 744   ensightGeomFile << support->getName() << endl;
 745
 746   // get geom types
 747   medEntityMesh entity = support->getEntity();
 748   int nbTypes = support->getNumberOfTypes();
 749   const medGeometryElement* geoType = support->getTypes();
 750
 751   const int * connectivity = 0;
 752   const int * elemConnectivity = 0;
 753   const int * index = 0;
 754   int j;
 755
 756   const MESH* mesh = support->getMesh()->convertInMESH() ; // we write unstructured only
 757
 758   // COORDINATES                                                              Gold ASCII
 759   // ===================================================================================
 760   // In Ensight, coordinates of nodes of support elements are in MED_NO_INTERLACE mode.
 761   // We are to write only nodes belonging to elements of the support and
 762   // nodal connectivity should refer to these nodes.
 763   map<int, int> med2ensIds;
 764   map<int, int>::iterator medEnsIt;
 765   int SpaceDimension = mesh->getSpaceDimension() ;
 766   int NumberOfNodes  = mesh->getNumberOfNodes() ;
 767   string zeroStr = " 0.00000e+00";
 768   // -----------------------------------
 769   if ( support->isOnAllElements() )
 770   {
 771     // nb of nodes
 772     ensightGeomFile << "coordinates" << endl
 773                     << setw(iw) << NumberOfNodes << endl ;
 774
 775     // coordinates
 776     for (j=0; j<SPACE_DIM; j++) { // loop on dimensions
 777       if ( j < SpaceDimension ) {
 778         const double *coordinate = mesh->getCoordinates(MED_FULL_INTERLACE) + j;
 779         for (int i=0; i<NumberOfNodes; i++, coordinate += SpaceDimension)
 780           ensightGeomFile << setw(12) << (float) *coordinate << endl;
 781       }
 782       else {
 783         for (int i=0; i<NumberOfNodes; i++)
 784           ensightGeomFile << zeroStr << endl;
 785       }
 786     }
 787   }
 788   // -----------------------------------
 789   else // support is not on all elements
 790   {
 791     // nb of nodes
 792     getSupportNodes( support, med2ensIds );
 793     NumberOfNodes = med2ensIds.size();
 794     ensightGeomFile << "coordinates" << endl
 795                     << setw(iw) << NumberOfNodes << endl ;
 796
 797     // coordinates
 798     for ( j=0; j<SPACE_DIM; j++) { // loop on dimensions
 799       medEnsIt = med2ensIds.begin();
 800       if ( j < SpaceDimension ) {
 801         const double *coordinate = mesh->getCoordinates(MED_FULL_INTERLACE) + j;
 802         for (int i=0; i<NumberOfNodes; i++, ++medEnsIt )
 803           ensightGeomFile << setw(12)
 804                           << (float) coordinate[ (medEnsIt->first-1) * SpaceDimension] << endl;
 805       }
 806       else {
 807         for (int i=0; i<NumberOfNodes; i++)
 808           ensightGeomFile << zeroStr << endl;
 809       }
 810     }
 811     // assign local node ids
 812     for ( medEnsIt = med2ensIds.begin(), j=1; j<=NumberOfNodes; j++, ++medEnsIt )
 813       medEnsIt->second = j;
 814   }
 815
 816   // CONNECTIVITY                                                             Gold ASCII
 817   // ===================================================================================
 818   // loop on types
 819   for (int i=0; i<nbTypes; i++)
 820   {
 821     const medGeometryElement    medType = geoType[i];
 822     const TEnSightElemType& ensightType = getEnSightType(medType);
 823     const int numberOfCell              = support->getNumberOfElements(medType);
 824     int nbCellNodes                     = ensightType._medIndex.size();
 825
 826     // type name and nb cells
 827     ensightGeomFile << ensightType._name        << endl
 828                     << setw(iw) << numberOfCell << endl;
 829
 830     // -----------------------------------
 831     if ( support->isOnAllElements() )
 832     {
 833 #ifdef ELEMENT_ID_GIVEN
 834       for ( j = 1; j <= numberOfCell; j++)
 835         ensightGeomFile << setw(iw) << j << endl;
 836 #endif
 837
 838       if ( nbCellNodes > 1 ) // STANDARD ELEMENTS connectivity
 839       {
 840         connectivity = mesh->getConnectivity( MED_NODAL, entity, medType);
 841         for (j = 0 ; j < numberOfCell; j++, connectivity += nbCellNodes) {
 842           for (int k=0; k<nbCellNodes; k++)
 843             ensightGeomFile << setw(iw) << connectivity[ ensightType._medIndex[k] ];
 844           ensightGeomFile << endl ;
 845         }
 846       }
 847       else if ( entity == MED_NODE ) // NODES connectivity
 848       {
 849         for ( j = 1; j <= numberOfCell; j++)
 850           ensightGeomFile << setw(iw) << j << endl;
 851       }
 852       else if ( medType == MED_POLYGON ) // POLYGONs connectivity
 853       {
 854         int nbStdCells = mesh->getGlobalNumberingIndex(entity)[i]-1;
 855         connectivity   = mesh->getConnectivity( MED_NODAL, entity, MED_ALL_ELEMENTS);
 856         index          = mesh->getConnectivityIndex(MED_NODAL, entity) + nbStdCells;
 857         // number of nodes in each element
 858         const int* ind = index;
 859         for (j = 0 ; j < numberOfCell; j++, ++ind)
 860           ensightGeomFile << setw(iw) << ( ind[1] - ind[0] ) << endl;
 861
 862         // connectivity
 863         for (j = 0; j < numberOfCell; j++, ++index) {
 864           nbCellNodes = index[1] - index[0];
 865           for (int k=0; k<nbCellNodes; k++, ++connectivity)
 866             ensightGeomFile << setw(iw) << *connectivity;
 867           ensightGeomFile << endl;
 868         }
 869       }
 870       else // POLYHEDRA connectivity
 871       {
 872         int nbStdCells = mesh->getGlobalNumberingIndex(entity)[i]-1;
 873         connectivity   = mesh->getConnectivity( MED_NODAL, entity, medType);
 874         index          = mesh->getConnectivityIndex(MED_NODAL, entity) + nbStdCells;
 875         ostringstream nbNodesInFace, faceConn;
 876         for ( j = 0; j < numberOfCell; ++j )
 877         {
 878           int nbFaces = 0, nbNodes = 0;
 879           for ( int k = index[j]; k < index[j+1]; ++k )
 880             if ( connectivity[k-1] == -1 )
 881             {
 882               faceConn << endl;
 883               nbNodesInFace << setw(iw) << nbNodes << endl;
 884               nbNodes = 0;
 885               ++nbFaces;
 886             }
 887             else
 888             {
 889               faceConn << setw(iw) << connectivity[k-1] ;
 890               ++nbNodes;
 891             }
 892           faceConn << endl;
 893           nbNodesInFace << setw(iw) << nbNodes << endl;
 894           ensightGeomFile << setw(iw) << nbFaces+1 << endl ;// nb of faces in each polyhedron
 895         }
 896         ensightGeomFile << nbNodesInFace.str();// number of nodes in each face
 897         ensightGeomFile << faceConn.str();// connectivity of each face
 898       }
 899     }
 900     // -----------------------------------
 901     else // support is not on all elements
 902     {
 903       const int *number = support->getNumber(medType);
 904
 905 #ifdef ELEMENT_ID_GIVEN
 906       for ( j = 0; j < numberOfCell; j++)
 907         ensightGeomFile << setw(iw) << number[j] << endl;
 908 #endif
 909       if ( nbCellNodes > 1 ) // STANDARD ELEMENTS connectivity
 910       {
 911         connectivity = mesh->getConnectivity( MED_NODAL, entity, MED_ALL_ELEMENTS);
 912         index = mesh->getConnectivityIndex(MED_NODAL, entity);
 913
 914         for (j=0; j<numberOfCell; j++) {
 915           int elem = number[j];
 916           elemConnectivity = connectivity + index[elem-1]-1;
 917           for (int k=0; k<nbCellNodes; k++) {
 918             int node = elemConnectivity[ ensightType._medIndex[k] ];
 919             ensightGeomFile << setw(iw) << med2ensIds[ node ];
 920           }
 921           ensightGeomFile << endl;
 922         }
 923       }
 924       else if ( entity == MED_NODE )  // NODES connectivity
 925       {
 926         for (j=0; j<numberOfCell; j++) {
 927           int node = med2ensIds[ number[j] ];
 928           ensightGeomFile << setw(iw) << node << endl ;
 929         }
 930       }
 931       else if ( medType == MED_POLYGON ) // POLYGONs connectivity
 932       {
 933         connectivity   = mesh->getConnectivity( MED_NODAL, entity, MED_ALL_ELEMENTS);
 934         index          = mesh->getConnectivityIndex(MED_NODAL, entity);
 935         // number of nodes in each element
 936         for (j = 0 ; j < numberOfCell; j++) {
 937           int elem = number[j];
 938           ensightGeomFile << setw(iw) << ( index[elem] - index[elem-1] ) << endl;
 939         }
 940         // connectivity
 941         for ( j = 0; j < numberOfCell; ++j ) {
 942           int elem = number[ j ];
 943           elemConnectivity   = connectivity + index[ elem-1 ]-1;
 944           const int* connEnd = connectivity + index[ elem   ]-1;
 945           while ( elemConnectivity < connEnd )
 946             ensightGeomFile << setw(iw) << med2ensIds[ *elemConnectivity++ ];
 947           ensightGeomFile << endl;
 948         }
 949       }
 950       else // POLYHEDRA connectivity
 951       {
 952         connectivity   = mesh->getConnectivity( MED_NODAL, entity, medType);
 953         index          = mesh->getConnectivityIndex(MED_NODAL, entity);
 954         ostringstream nbNodesInFace, faceConn;
 955         for ( j = 0; j < numberOfCell; ++j )
 956         {
 957           int elem = number[j];
 958           int nbFaces = 0, nbNodes = 0;
 959           for ( int k = index[elem]; k < index[elem+1]; ++k )
 960             if ( connectivity[k-1] == -1 )
 961             {
 962               faceConn << endl;
 963               nbNodesInFace << setw(iw) << nbNodes << endl;
 964               nbNodes = 0;
 965               ++nbFaces;
 966             }
 967             else
 968             {
 969               faceConn << setw(iw) << connectivity[k-1] ;
 970               ++nbNodes;
 971             }
 972           faceConn << endl;
 973           nbNodesInFace << setw(iw) << nbNodes << endl;
 974           ensightGeomFile << setw(iw) << nbFaces+1 << endl ;// nb of faces in each polyhedron
 975         }
 976         ensightGeomFile << nbNodesInFace.str();// number of nodes in each face
 977         ensightGeomFile << faceConn.str();// connectivity of each face
 978       }
 979     }
 980   }
 981
 982   mesh->removeReference();
 983
 984 }  // writePartGoldASCII()
 985
 986 //================================================================================
 987 /*!
 988  * \brief Write support as an Ensight6 part
 989  */
 990 //================================================================================
 991
 992 void ENSIGHT_MESH_WRONLY_DRIVER::writePart6Binary(_BinaryFileWriter& ensightGeomFile,
 993                                                   const SUPPORT*     support) const
 994 {
 995   // part number
 996   int partNum = getPartNumber( support );
 997   ensightGeomFile.addString( STRING("part ") << partNum );
 998   if ( !partNum )
 999     throw MED_EXCEPTION ( LOCALIZED( STRING("Internal error: invalid part number")));
1000
1001   // group/mesh name
1002   ensightGeomFile.addString( support->getName() );
1003
1004   // get geom types
1005   medEntityMesh entity = support->getEntity();
1006   int nbTypes = support->getNumberOfTypes();
1007   const medGeometryElement* geoType = support->getTypes();
1008
1009   const MESH* mesh = support->getMesh()->convertInMESH() ; // we write unstructured only
1010
1011   int j = 1;
1012   const int * connectivity = 0;
1013   if ( entity != MED_NODE )
1014     connectivity = mesh->getConnectivity( MED_NODAL, entity, MED_ALL_ELEMENTS);
1015   const int * elemConnectivity = connectivity;
1016
1017   // CONNECTIVITY                                                       Ensight 6 binary
1018   // ===================================================================================
1019   // loop on types
1020   for (int i=0; i<nbTypes; i++)
1021   {
1022     const medGeometryElement    medType = geoType[i];
1023     const TEnSightElemType& ensightType = getEnSightType(medType);
1024     int nbCellNodes = ensightType._medIndex.size();
1025     if ( nbCellNodes == 0 )
1026       continue; // poly?
1027
1028     // type name and nb cells
1029     int numberOfCell = support->getNumberOfElements(medType);
1030     ensightGeomFile.addString( ensightType._name );
1031     ensightGeomFile.addInt( numberOfCell );
1032
1033     vector<int> nodeIds;
1034     // -------------------------------------------------
1035     if ( support->isOnAllElements() )
1036     {
1037 #ifdef ELEMENT_ID_GIVEN
1038       nodeIds.resize( numberOfCell );
1039       for ( j = 1; j <= numberOfCell; j++)
1040         nodeIds[ j-1 ] = j;
1041       ensightGeomFile.addInt( nodeIds );
1042 #endif
1043       if ( entity == MED_NODE ) {
1044 #if !defined(ELEMENT_ID_GIVEN)
1045         nodeIds.resize( numberOfCell * nbCellNodes);
1046         for ( j = 1; j <= numberOfCell; j++)
1047           nodeIds[ j-1 ] = j;
1048 #endif
1049       }
1050       else {
1051         nodeIds.clear();
1052         nodeIds.reserve( numberOfCell * nbCellNodes );
1053         for (j = 0 ; j < numberOfCell; j++, elemConnectivity += nbCellNodes)
1054           for (int k=0; k<nbCellNodes; k++)
1055             nodeIds.push_back( elemConnectivity[ ensightType._medIndex[k] ]);
1056       }
1057       ensightGeomFile.addInt( nodeIds );
1058     }
1059     // -------------------------------------------------
1060     else // support is not on all elements
1061     {
1062       const int *number = support->getNumber(medType);
1063
1064 #ifdef ELEMENT_ID_GIVEN
1065       ensightGeomFile.addInt( number, numberOfCell );
1066 #endif
1067       if ( entity == MED_NODE ) {
1068         ensightGeomFile.addInt( number, numberOfCell );
1069       }
1070       else {
1071         const int* index = mesh->getConnectivityIndex(MED_NODAL, entity);
1072
1073         nodeIds.reserve( numberOfCell * nbCellNodes);
1074         for (j=0; j<numberOfCell; j++) {
1075           int elem = number[j];
1076           elemConnectivity = connectivity + index[elem-1]-1;
1077           for (int k=0; k<nbCellNodes; k++)
1078             nodeIds.push_back( elemConnectivity[ ensightType._medIndex[k] ]);
1079         }
1080         ensightGeomFile.addInt( nodeIds );
1081       }
1082     }
1083   } // loop on types
1084
1085   mesh->removeReference();
1086
1087
1088 } // writePart6Binary()
1089
1090 //================================================================================
1091 /*!
1092  * \brief Write support as an Ensight6 part
1093  */
1094 //================================================================================
1095
1096 void ENSIGHT_MESH_WRONLY_DRIVER::writePart6ASCII(ofstream&      ensightGeomFile,
1097                                                  const SUPPORT* support) const
1098 {
1099   const int iw = 8;
1100
1101   // part number
1102   int partNum = getPartNumber( support );
1103   ensightGeomFile << "part " << partNum << endl;
1104   if ( !partNum )
1105     throw MED_EXCEPTION ( LOCALIZED( STRING("Internal error: invalid part number")));
1106
1107   // group/mesh name
1108   ensightGeomFile << support->getName() << endl;
1109
1110   // get geom types
1111   medEntityMesh entity = support->getEntity();
1112   int nbTypes = support->getNumberOfTypes();
1113   const medGeometryElement* geoType = support->getTypes();
1114
1115   const MESH* mesh = support->getMesh()->convertInMESH() ; // we write unstructured only
1116
1117   int j = 1;
1118   const int * connectivity = 0;
1119   if ( entity != MED_NODE )
1120     connectivity = mesh->getConnectivity( MED_NODAL,entity, MED_ALL_ELEMENTS);
1121   const int * elemConnectivity = connectivity;
1122
1123   // CONNECTIVITY                                                        Ensight 6 ASCII
1124   // ===================================================================================
1125   // loop on types
1126   for (int i=0; i<nbTypes; i++)
1127   {
1128     const medGeometryElement    medType = geoType[i];
1129     const TEnSightElemType& ensightType = getEnSightType(medType);
1130     int nbCellNodes = ensightType._medIndex.size();
1131     if ( nbCellNodes == 0 )
1132       continue; // poly?
1133
1134     // type name and nb cells
1135     int numberOfCell = support->getNumberOfElements(medType);
1136     ensightGeomFile << ensightType._name       << endl
1137                     << setw(iw) << numberOfCell << endl;
1138
1139     // -------------------------------------------------
1140     if ( support->isOnAllElements() )
1141     {
1142       if ( entity == MED_NODE ) {
1143         for ( j = 1; j <= numberOfCell; j++) {
1144 #ifdef ELEMENT_ID_GIVEN
1145           ensightGeomFile << setw(iw) << j;
1146 #endif
1147           ensightGeomFile << setw(iw) << j << endl;
1148         }
1149       }
1150       else {
1151         for (j = 1 ; j <= numberOfCell; j++, elemConnectivity += nbCellNodes) {
1152 #ifdef ELEMENT_ID_GIVEN
1153           ensightGeomFile << setw(iw) << elem++;
1154 #endif
1155           for (int k=0; k<nbCellNodes; k++)
1156           {
1157             ensightGeomFile << setw(iw) << elemConnectivity[ ensightType._medIndex[k] ];
1158           }
1159           ensightGeomFile << endl ;
1160         }
1161       }
1162     }
1163     // -------------------------------------------------
1164     else  // support is not on all elements
1165     {
1166       const int *number = support->getNumber(medType);
1167       if ( entity == MED_NODE ) {
1168         for (j=0; j<numberOfCell; j++) {
1169           int node = number[j];
1170 #ifdef ELEMENT_ID_GIVEN
1171           ensightGeomFile << setw(iw) << node;
1172 #endif
1173           ensightGeomFile << setw(iw) << node << endl ;
1174         }
1175       }
1176       else {
1177         const int* index = mesh->getConnectivityIndex(MED_NODAL, entity);
1178
1179         for (j=0; j<numberOfCell; j++) {
1180           int elem = number[j];
1181 #ifdef ELEMENT_ID_GIVEN
1182           ensightGeomFile << setw(iw) << elem;
1183 #endif
1184           elemConnectivity = connectivity + index[elem-1]-1;
1185           for (int k=0; k<nbCellNodes; k++)
1186           {
1187             ensightGeomFile << setw(iw) << elemConnectivity[ ensightType._medIndex[k] ];
1188           }
1189           ensightGeomFile << endl ;
1190         }
1191       }
1192     }
1193   } // loop on types
1194
1195   mesh->removeReference();
1196
1197
1198 } // writePart6ASCII()
1199
1200 //================================================================================
1201 /*!
1202  * \brief Return nb of part to write
1203  */
1204 //================================================================================
1205
1206 int ENSIGHT_MESH_WRONLY_DRIVER::nbPartsToWrite() const
1207 {
1208   int nbParts = 0;
1209   nbParts += (int) isToWriteEntity( MED_CELL, _ptrMesh );
1210   nbParts += (int) isToWriteEntity( MED_FACE, _ptrMesh );
1211   nbParts += (int) isToWriteEntity( MED_EDGE, _ptrMesh );
1212
1213   // all groups
1214   for ( int ent = MED_CELL; ent < MED_ALL_ENTITIES; ++ent ) {
1215     int nbGroups = _ptrMesh->getNumberOfGroups(medEntityMesh(ent));
1216     nbParts += nbGroups;
1217   }
1218   return nbParts;
1219 }
1220
1221 // ================================================================================
1222 // RDONLY
1223 // ================================================================================
1224
1225 ENSIGHT_MESH_RDONLY_DRIVER::ENSIGHT_MESH_RDONLY_DRIVER()
1226   : ENSIGHT_MESH_DRIVER(), _indexInCaseFile(1)
1227 {
1228 }
1229
1230 ENSIGHT_MESH_RDONLY_DRIVER::ENSIGHT_MESH_RDONLY_DRIVER(const string & fileName,
1231                                                        MESH *         ptrMesh,
1232                                                        int            index)
1233   : ENSIGHT_MESH_DRIVER(fileName,ptrMesh,RDONLY), _indexInCaseFile( index )
1234 {
1235 }
1236
1237 ENSIGHT_MESH_RDONLY_DRIVER::ENSIGHT_MESH_RDONLY_DRIVER(const ENSIGHT_MESH_RDONLY_DRIVER & driver) : ENSIGHT_MESH_DRIVER(driver), _indexInCaseFile( driver._indexInCaseFile )
1238 {
1239 }
1240
1241 ENSIGHT_MESH_RDONLY_DRIVER::~ENSIGHT_MESH_RDONLY_DRIVER()
1242 {
1243 }
1244
1245 GENDRIVER * ENSIGHT_MESH_RDONLY_DRIVER::copy() const
1246 {
1247   return new ENSIGHT_MESH_RDONLY_DRIVER(*this) ;
1248 }
1249
1250 void ENSIGHT_MESH_RDONLY_DRIVER::write() const throw (MEDEXCEPTION)
1251 {
1252   throw MEDEXCEPTION("ENSIGHT_MESH_RDONLY_DRIVER::write : Can't write with a RDONLY driver !");
1253 }
1254
1255 void ENSIGHT_MESH_RDONLY_DRIVER::merge ( const GENDRIVER& driver )
1256 {
1257   _CaseFileDriver_User::merge( driver );
1258
1259   const ENSIGHT_MESH_RDONLY_DRIVER* other =
1260     dynamic_cast< const ENSIGHT_MESH_RDONLY_DRIVER* >( &driver );
1261   if ( other ) {
1262     if ( _indexInCaseFile < other->_indexInCaseFile )
1263       _indexInCaseFile = other->_indexInCaseFile;
1264   }
1265 }
1266
1267 //================================================================================
1268 /*!
1269  * \brief Read mesh in all supported formats
1270  */
1271 //================================================================================
1272
1273 void ENSIGHT_MESH_RDONLY_DRIVER::read() throw (MEDEXCEPTION)
1274 {
1275   const char * LOC = "ENSIGHT_MESH_RDONLY_DRIVER::read() : " ;
1276   BEGIN_OF_MED(LOC);
1277
1278   openConst(false); // check if can read case file
1279
1280   _CaseFileDriver caseFile( getCaseFileName(), this);
1281   caseFile.read();
1282   caseFile.setDataFileName( _indexInCaseFile, this ); // data from Case File is passed here
1283
1284   openConst(true); // check if can read data file
1285
1286   cout << "-> Entering into the geometry file " << getDataFileName() << endl  ;
1287
1288   MESH* mesh = (MESH*) getMesh();
1289
1290   _InterMed* imed = new _InterMed();
1291   imed->_medMesh = mesh;
1292   imed->_isOwnMedMesh = false;
1293   imed->_needSubParts = ( caseFile.getNbVariables() > 0 );
1294   imed->groupes.reserve(1000);
1295
1296   // to let field drivers know eventual indices of values
1297   if ( imed->_needSubParts )
1298     setInterData( imed );
1299
1300   if ( isBinaryDataFile( getDataFileName() )) // binary
1301   {
1302     if ( isGoldFormat() ) // Gold
1303     {
1304       readGoldBinary( *imed );
1305     }
1306     else // EnSight6
1307     {
1308       read6Binary( *imed );
1309     }
1310   }
1311   else // ASCII
1312   {
1313     if ( isGoldFormat() ) // Gold
1314     {
1315       readGoldASCII( *imed );
1316     }
1317     else // EnSight6
1318     {
1319       read6ASCII( *imed );
1320     }
1321   }
1322
1323   if ( _isMadeByMed && !imed->groupes.empty() ) {
1324     mesh->_name = imed->groupes[0].nom;
1325     imed->groupes[0].nom = "SupportOnAll_";
1326     imed->groupes[0].nom += entNames[MED_CELL];
1327   }
1328   else {
1329     mesh->_name = theDefaultMeshName;
1330   }
1331   mesh->_spaceDimension = SPACE_DIM;
1332   mesh->_numberOfNodes  = imed->points.size() - imed->nbMerged( MED_POINT1 );
1333   mesh->_coordinate     = imed->getCoordinate();
1334
1335   //Construction des groupes
1336   imed->getGroups(mesh->_groupCell,
1337                   mesh->_groupFace,
1338                   mesh->_groupEdge,
1339                   mesh->_groupNode, mesh);
1340
1341   mesh->_connectivity = imed->getConnectivity();
1342
1343   mesh->createFamilies();
1344
1345   // add attributes to families
1346   set<string> famNames;
1347   for (medEntityMesh entity=MED_CELL; entity<MED_ALL_ENTITIES; ++entity)
1348   {
1349     int i, nb = mesh->getNumberOfFamilies(entity);
1350     for ( i = 1; i <= nb; ++i ) {
1351       FAMILY* f = const_cast<FAMILY*>( mesh->getFamily( entity, i ));
1352       f->setNumberOfAttributes( 1 );
1353       int* attIDs = new int[1];
1354       attIDs[0] = 1;
1355       f->setAttributesIdentifiers( attIDs );
1356       int* attVals = new int[1];
1357       attVals[0] = 1;
1358       f->setAttributesValues( attVals );
1359       string* attDescr = new string[1];
1360       attDescr[0] = "med_family";
1361       f->setAttributesDescriptions( attDescr );
1362       delete [] attDescr;
1363       if ( f->getName().length() > 31 ) // limit a name length
1364         f->setName( STRING("FAM_") << f->getIdentifier());
1365       // setAll() for groups
1366       nb = mesh->getNumberOfGroups(entity);
1367       for ( i = 1; i <= nb; ++i ) {
1368         GROUP * g = const_cast<GROUP*>( mesh->getGroup( entity, i ));
1369         if (mesh->getNumberOfElements( entity, MED_ALL_ELEMENTS ) ==
1370             g->getNumberOfElements( MED_ALL_ELEMENTS ))
1371         {
1372           g->setAll( true );
1373           g->update();
1374         }
1375       }
1376     }
1377   }
1378
1379   if ( !imed->_needSubParts )
1380     delete imed;
1381
1382   END_OF_MED(LOC);
1383 }
1384
1385 //================================================================================
1386 /*!
1387  * \brief Read mesh in Gold ASCII format
1388  */
1389 //================================================================================
1390
1391 void ENSIGHT_MESH_RDONLY_DRIVER::readGoldASCII(_InterMed & imed)
1392 {
1393   const char * LOC ="ENSIGHT_MESH_RDONLY_DRIVER::readGoldASCII() : ";
1394   BEGIN_OF_MED(LOC);
1395
1396   _ASCIIFileReader geoFile( getDataFileName() );
1397
1398   if ( isSingleFileMode() ) {
1399     int curTimeStep = 1;
1400     while ( curTimeStep++ < getIndexInDataFile() ) {
1401       while ( !geoFile.isTimeStepEnd())
1402         geoFile.getLine();
1403     }
1404     while ( !geoFile.isTimeStepBeginning() )
1405       geoFile.getLine();
1406   }
1407   // ----------------------
1408   // Read mesh description
1409   // ----------------------
1410   {
1411     string descriptionLine1 = geoFile.getLine();
1412     string descriptionLine2 = geoFile.getLine();
1413
1414     // find out if the file was created by MED driver
1415     int prefixSize = strlen( theDescriptionPrefix );
1416     _isMadeByMed = ( descriptionLine1.substr(0, prefixSize ) == theDescriptionPrefix );
1417
1418     if ( _isMadeByMed )
1419       descriptionLine1 = descriptionLine1.substr( prefixSize );
1420     _ptrMesh->setDescription( descriptionLine1 + descriptionLine2 );
1421   }
1422
1423   // ----------------------------------------
1424   // Find out presence of node/elem numbers
1425   // ----------------------------------------
1426
1427   // EnSight User Manual (for v8) says:
1428   //    You do not have to assign node IDs. If you do, the element connectivities are
1429   //    based on the node numbers. If you let EnSight assign the node IDs, the nodes
1430   //    are considered to be sequential starting at node 1, and element connectivity is
1431   //    done accordingly. If node IDs are set to off, they are numbered internally;
1432   //    however, you will not be able to display or query on them. If you have node
1433   //    IDs in your data, you can have EnSight ignore them by specifying "node id
1434   //    ignore." Using this option may reduce some of the memory taken up by the
1435   //    Client and Server, but display and query on the nodes will not be available.
1436
1437   // read "node|element id <off|given|assign|ignore>"
1438   geoFile.getWord(); geoFile.getWord();
1439   string nodeIds = geoFile.getWord();
1440   geoFile.getWord(); geoFile.getWord();
1441   string elemIds = geoFile.getWord();
1442
1443   bool haveNodeIds = ( nodeIds == "given" || nodeIds == "ignore" );
1444   bool haveElemIds = ( elemIds == "given" || elemIds == "ignore" );
1445
1446   // extents: xmin xmax ymin ymax zmin zmax
1447   vector<double> extents;
1448   geoFile.toNextLine();
1449   if ( strncmp( "extents", geoFile.getCurrentPtr(), 7 ) == 0 ) {
1450     geoFile.skip( /*width =*/ 7, /*nbLines =*/ 1 );
1451     extents.reserve( 6 );
1452     while ( extents.size() < extents.capacity() )
1453       extents.push_back( geoFile.getReal() );
1454   }
1455
1456   typedef map<int,_noeud>::iterator INoeud;
1457   map<int,_noeud> & points = imed.points;
1458   INoeud firstNode;
1459
1460   _groupe * partGroupe = 0;
1461   int partNum = 0, nbParts = 0;
1462
1463   while ( !geoFile.isTimeStepEnd() )
1464   {
1465     string word, restLine, line = geoFile.getLine();
1466     TStrTool::split( line, word, restLine );
1467
1468     const TEnSightElemType & partType = getEnSightType( word );
1469     if ( partType._medType != MED_ALL_ELEMENTS )
1470     {
1471       //  Unstructured element type encounters
1472       // --------------------------------------
1473       int  nbElemNodes = partType._medType % 100;
1474       int      nbElems = geoFile.getInt(); // ne
1475       bool     isGhost = isGhostType( word );
1476       int    nodeShift = points.empty() ? 0 : points.rbegin()->first;
1477
1478       // read element ids
1479       vector<int> elemIds;
1480       if ( haveElemIds ) {
1481         elemIds.reserve( nbElems );
1482         while ((int) elemIds.size() < nbElems )
1483           elemIds.push_back( geoFile.getInt() ); // id_e
1484       }
1485       if ( isGhost ) { // do not store ghost elements (?)
1486         int nbNodes = nbElems * nbElemNodes;
1487         if ( partType._name == "nsided" ) // polygons
1488         {
1489           for ( int i = 0; i < nbElems; ++i )
1490             nbNodes += geoFile.getInt();
1491           geoFile.skip( nbNodes * INT_WIDTH_GOLD, /*nbLines = */nbElems );
1492         }
1493         else if ( partType._name == "nfaced" ) // polyhedrons
1494         {
1495           int nbFaces = 0;
1496           for ( int i = 0; i < nbElems; ++i )
1497             nbFaces += geoFile.getInt();
1498           for ( int f = 0; f < nbFaces; ++f )
1499             nbNodes += geoFile.getInt();
1500           geoFile.skip( nbNodes * INT_WIDTH_GOLD, /*nbLines = */nbFaces );
1501         }
1502         else // standard types
1503         {
1504           geoFile.skip( nbNodes, nbElemNodes, INT_WIDTH_GOLD );
1505         }
1506         continue;
1507       }
1508
1509       // add a group corresponding to subPart (geoType)
1510       imed.groupes.push_back(_groupe());
1511       _groupe & groupe = imed.groupes.back();
1512       groupe.mailles.resize( nbElems );
1513
1514       // find out if "coordinates" has already been encountered
1515       _SubPartDesc coordDesc( partNum , "coordinates");
1516       map< _SubPartDesc, _SubPart >::iterator descPart =
1517         imed._subPartDescribed.find( coordDesc );
1518       bool haveCoords = ( descPart != imed._subPartDescribed.end() );
1519       if ( haveCoords ) {
1520         firstNode = descPart->second.myFirstNode;
1521         nodeShift -= descPart->second.myNbNodes;
1522       }
1523
1524       // read poly element data
1525       bool isPoly = ( !nbElemNodes );
1526       vector<int> nbElemNodesVec( 1, nbElemNodes);
1527       vector<int> nbElemFaces, nbFaceNodes;
1528       if ( partType._name == "nsided" ) // polygons
1529       {
1530         nbElemNodesVec.resize( nbElems );
1531         for ( int i = 0; i < nbElems; ++i )
1532           nbElemNodesVec[ i ] = geoFile.getInt(); // npi
1533       }
1534       else if ( partType._name == "nfaced" ) // polyhedrons
1535       {
1536         nbElemFaces.resize( nbElems );
1537         nbElemNodesVec.resize( nbElems );
1538         int totalNbFaces = 0;
1539         for ( int i = 0; i < nbElems; ++i )
1540           totalNbFaces += ( nbElemFaces[ i ] = geoFile.getInt() ); // nf_ei
1541
1542         nbFaceNodes.resize( totalNbFaces );
1543         vector<int>::iterator nbFN = nbFaceNodes.begin();
1544         for ( int i = 0; i < nbElems; ++i ) {
1545           nbElemNodesVec[ i ] = 0;
1546           for ( int nbFaces = nbElemFaces[ i ]; nbFaces; --nbFaces, ++nbFN )
1547             nbElemNodesVec[ i ] += ( *nbFN = geoFile.getInt() ); // np(f_ei)
1548         }
1549       }
1550       // iterator returning nbElemNodes for standard elems and
1551       // next value from nbElemNodesVec for poly elements
1552       _ValueIterator<int> nbElemNodesIt( & nbElemNodesVec[0], isPoly ? 1 : 0);
1553
1554       // iterator returning values form partType._medIndex for standard elems
1555       // and node index (n) for poly elements
1556       int n;
1557       _ValueIterator<int> medIndexIt( isPoly ? &n : &partType._medIndex[0],
1558                                       isPoly ? 0  : 1);
1559       // read connectivity
1560       _maille ma( partType._medType, nbElemNodes );
1561       ma.sommets.resize( nbElemNodes );
1562       INoeud node;
1563       for ( int i = 0; i < nbElems; ++i ) {
1564         _ValueIterator<int> medIndex = medIndexIt;
1565         nbElemNodes = nbElemNodesIt.next();
1566         if ((int) ma.sommets.size() != nbElemNodes )
1567           ma.sommets.resize( nbElemNodes );
1568         for ( n = 0; n < nbElemNodes; ++n ) {
1569           int nodeID = geoFile.getInt(); // nn_ei
1570           if ( haveCoords )
1571             node = points.find( nodeID + nodeShift );
1572           else
1573             node = points.insert( make_pair( nodeID + nodeShift, _noeud())).first;
1574           ma.sommets[ medIndex.next() ] = node;
1575         }
1576         if ( haveElemIds )
1577           ma.setOrdre( elemIds[ i ] );
1578         groupe.mailles[i] = imed.insert(ma);
1579       }
1580       // store nb nodes in polyhedron faces
1581       if ( !nbFaceNodes.empty() ) {
1582         const int* nbFaceNodesPtr = & nbFaceNodes[0];
1583         for ( int i = 0; i < nbElems; ++i ) {
1584           vector<int> & nbNodesByFace = imed.polyherdalNbFaceNodes[ &groupe.maille( i ) ];
1585           nbNodesByFace.assign( nbFaceNodesPtr, nbFaceNodesPtr + nbElemFaces[ i ] );
1586           nbFaceNodesPtr += nbElemFaces[ i ];
1587         }
1588       }
1589       // create subPart for "coordinates"
1590       if ( !haveCoords ) {
1591         _SubPart & coordSubPart = imed._subPartDescribed[ coordDesc ];
1592         coordSubPart.myFirstNode = points.insert( make_pair( 1 + nodeShift, _noeud())).first;
1593       }
1594       // add subPart group to part group
1595       int groupeIndex = imed.groupes.size();
1596       partGroupe->groupes.push_back( groupeIndex );
1597
1598       // create subPart
1599       _SubPart subPart( partNum, partType._name );
1600       subPart.myNbCells = nbElems;
1601       subPart.myCellGroupIndex = groupeIndex;
1602       imed.addSubPart( subPart );
1603     }
1604     else if ( word == "coordinates" )
1605     {
1606       // Local node coordinates of a part
1607       // ------------------------------------
1608       int nbNodes = geoFile.getInt(); // nn
1609
1610       // read node ids
1611       vector<int> nodeIds;
1612       if ( haveNodeIds ) {
1613         nodeIds.reserve( nbNodes );
1614         while ((int) nodeIds.size() < nbNodes )
1615           nodeIds.push_back( geoFile.getInt() ); // id_n
1616       }
1617
1618       // find out if "coordinates" has already been add at reading connectivity
1619       _SubPartDesc coordDesc( partNum , "coordinates");
1620       map< _SubPartDesc, _SubPart >::iterator descPart =
1621         imed._subPartDescribed.find( coordDesc );
1622       bool haveCoords = ( descPart != imed._subPartDescribed.end() );
1623
1624       if ( haveCoords ) {
1625         // check that all nodes have been added
1626         firstNode = descPart->second.myFirstNode;
1627         descPart->second.myNbNodes = nbNodes;
1628         INoeud inoeud = firstNode, inoEnd = points.end();
1629         int id = inoeud->first, idEnd = id + nbNodes;
1630         for ( ; id < idEnd; ++id ) {
1631           if ( inoeud == inoEnd || inoeud->first > id ) {
1632             INoeud in = points.insert( inoeud, make_pair( id, _noeud() ));
1633             in->second.number = id;
1634             in->second.coord.resize( SPACE_DIM );
1635           } else {
1636             ++inoeud;
1637           }
1638         }
1639       }
1640       else {
1641         // add nodes
1642         int nodeShift = points.empty() ? 0 : points.rbegin()->first;
1643         for ( int iNode = 1; iNode <= nbNodes; ++iNode ) {
1644           INoeud inoeud = points.insert( points.end(), make_pair( iNode + nodeShift, _noeud()));
1645           inoeud->second.number = inoeud->first;
1646           inoeud->second.coord.resize( SPACE_DIM );
1647         }
1648         firstNode = points.find( 1 + nodeShift );
1649         // create "coordinates" subPart
1650         _SubPart & subPart  = imed._subPartDescribed[ coordDesc ];
1651         subPart.myNbNodes   = nbNodes;
1652         subPart.myFirstNode = firstNode;
1653       }
1654
1655       // read coordinates in no interlace mode
1656       INoeud endNode = points.end();
1657       for ( int j = 0; j < SPACE_DIM; ++j ) {
1658         for ( INoeud in = firstNode; in != endNode; ++in ) {
1659           _noeud & node = in->second;
1660           node.coord[ j ] = geoFile.getReal();
1661         }
1662       }
1663     }
1664     else if ( word == "part" )
1665     {
1666       // Another part encounters
1667       // -----------------------
1668       partNum = geoFile.getInt();
1669       nbParts++;
1670       geoFile.toNextLine();
1671
1672       string partName = geoFile.getLine();
1673       if ( partName.empty() )
1674         partName = "Part_" + restLine;
1675
1676       if (int( imed.groupes.capacity() - imed.groupes.size() ) < theMaxNbTypes )
1677         imed.groupes.reserve( size_t( 1.5 * imed.groupes.size() ));
1678       imed.groupes.push_back(_groupe());
1679       partGroupe = & imed.groupes.back();
1680       partGroupe->nom = partName;
1681       partGroupe->groupes.reserve( theMaxNbTypes );
1682     }
1683     else if ( word == "block" )
1684     {
1685       // Structured type
1686       // ------------------
1687       bool rectilinear = ( restLine.find( "rectilinear" ) != restLine.npos );
1688       bool uniform     = ( restLine.find( "uniform" )     != restLine.npos );
1689       bool curvilinear = ( !rectilinear && !uniform );
1690       bool iblanked    = ( restLine.find( "iblanked" )    != restLine.npos );
1691       bool with_ghost  = ( restLine.find( "with_ghost" )  != restLine.npos );
1692       bool range       = ( restLine.find( "range" )       != restLine.npos );
1693
1694       // dimension
1695       int I = geoFile.getInt();
1696       int J = geoFile.getInt();
1697       int K = geoFile.getInt();
1698       int NumberOfNodes = I*J*K;
1699       if ( !NumberOfNodes ) continue;
1700
1701       // range
1702       if ( range ) {
1703         vector<int> ijkRange; // imin imax jmin jmax kmin kmax
1704         ijkRange.reserve(6);
1705         while ( ijkRange.size() < 6 )
1706           ijkRange.push_back( geoFile.getInt() );
1707         I = ijkRange[1]-ijkRange[0]+1;
1708         J = ijkRange[3]-ijkRange[2]+1;
1709         K = ijkRange[5]-ijkRange[4]+1;
1710         NumberOfNodes = I*J*K;
1711       }
1712       // add nodes
1713       int nodeShift = points.empty() ? 0 : points.rbegin()->first;
1714       for ( int iNode = 1; iNode <= NumberOfNodes; ++iNode ) {
1715         INoeud inoeud = points.insert( points.end(), make_pair( iNode + nodeShift, _noeud()));
1716         _noeud & node = inoeud->second;
1717         node.number   = inoeud->first;
1718         node.coord.resize( SPACE_DIM );
1719       }
1720       INoeud firstNode = points.find( nodeShift + 1 );
1721       INoeud endNode   = points.end();
1722
1723       GRID grid; // calculator of unstructured data
1724       grid._iArrayLength   = I;
1725       grid._jArrayLength   = J;
1726       grid._kArrayLength   = K;
1727       grid._spaceDimension = SPACE_DIM;
1728       if ( J < 2 ) { grid._spaceDimension--; grid._jArrayLength = 0; }
1729       if ( K < 2 ) { grid._spaceDimension--; grid._kArrayLength = 0; }
1730       int nbElems = grid.getNumberOfElements(MED_CELL, MED_ALL_ELEMENTS);
1731
1732       if ( curvilinear ) // read coordinates for all nodes
1733       {
1734         for ( int j = 0; j < SPACE_DIM; ++j ) {
1735           for ( INoeud in = firstNode; in != endNode; ++in )
1736             in->second.coord[ j ] = geoFile.getReal();
1737         }
1738         grid._gridType   = MED_BODY_FITTED;
1739         grid._coordinate = new COORDINATE(3, 1, 0); // to avoid exception in convertInMESH()
1740       }
1741       else if ( rectilinear ) // read delta vectors with non-regular spacing
1742       {
1743         grid._iArray = (double*)geoFile.convertReals<double>( I );
1744         grid._jArray = (double*)geoFile.convertReals<double>( J );
1745         grid._kArray = (double*)geoFile.convertReals<double>( K );
1746         grid._gridType = MED_CARTESIAN;
1747       }
1748       else // uniform: read grid origine and delta vectors for regular spacing grid
1749       {
1750         TDblOwner xyzOrigin( (double*)geoFile.convertReals<double>( 3 ));
1751         TDblOwner xyzDelta ( (double*)geoFile.convertReals<double>( 3 ));
1752         // compute full delta vectors
1753         grid._iArray = new double[ I ];
1754         grid._jArray = new double[ J ];
1755         grid._kArray = new double[ K ];
1756         double* coors[SPACE_DIM] = { grid._iArray, grid._jArray, grid._kArray };
1757         int     size [SPACE_DIM] = { I, J, K };
1758         for ( int j = 0; j < SPACE_DIM; ++j ) {
1759           double* coo    = coors[ j ];
1760           double* cooEnd = coo + size[ j ];
1761           coo[0]         = xyzOrigin[ j ];
1762           while ( ++coo < cooEnd )
1763             *coo = coo[-1] + xyzDelta[ j ];
1764         }
1765         grid._gridType = MED_CARTESIAN;
1766       }
1767
1768       // iblanks
1769       if ( iblanked )
1770         geoFile.skip( NumberOfNodes, /*nbPerLine =*/ 1, INT_WIDTH_GOLD);
1771       // ghosts
1772       if ( with_ghost ) {
1773         geoFile.getWord(); // "ghost_flags"
1774         geoFile.skip( nbElems, /*nbPerLine =*/ 1, INT_WIDTH_GOLD);
1775       }
1776       // node ids
1777       if ( haveNodeIds && geoFile.lookAt( "node_ids" )) {
1778         geoFile.getWord(); // "node_ids"
1779         geoFile.skip( NumberOfNodes, /*nbPerLine =*/ 1, INT_WIDTH_GOLD);
1780       }
1781       // element ids
1782       if ( haveElemIds && geoFile.lookAt( "element_ids" ) ) {
1783         geoFile.getWord(); // "element_ids"
1784         geoFile.skip( nbElems, /*nbPerLine =*/ 1, INT_WIDTH_GOLD);
1785       }
1786
1787       // let GRID compute all coordinates and connectivity
1788       const MESH* unstruct = grid.convertInMESH();
1789       if ( !curvilinear )
1790       {
1791         const double * coo = unstruct->getCoordinates(MED_FULL_INTERLACE);
1792         typedef _ValueIterator< double > TCoordIt;
1793         TCoordIt xCoo( coo+0, grid._spaceDimension);
1794         TCoordIt yCoo( coo+1, grid._spaceDimension);
1795         TCoordIt zCoo( coo+2, grid._spaceDimension);
1796         if ( grid._spaceDimension < 3 ) zCoo = TCoordIt( grid._kArray, 0 );
1797         if ( grid._spaceDimension < 2 ) yCoo = TCoordIt( grid._jArray, 0 );
1798         for ( INoeud in = firstNode; in != endNode; ++in ) {
1799           _noeud& node = in->second;
1800           node.coord[ 0 ] = xCoo.next();
1801           node.coord[ 1 ] = yCoo.next();
1802           node.coord[ 2 ] = zCoo.next();
1803         }
1804       }
1805
1806       // store connectivity
1807       const int * conn = unstruct->getConnectivity( MED_NODAL, MED_CELL, MED_ALL_ELEMENTS );
1808       medGeometryElement elemType = grid.getElementType( MED_CELL, 1 );
1809       int  nbElemNodes = elemType % 100;
1810
1811       partGroupe->mailles.resize( nbElems );
1812       _maille ma( elemType, nbElemNodes );
1813       ma.sommets.resize( nbElemNodes );
1814
1815       for ( int i = 0, nIndex = 0; i < nbElems; ++i ) {
1816         for ( int n = 0; n < nbElemNodes; ++n ) {
1817           int nodeID = conn[ nIndex++ ];
1818           ma.sommets[n] = points.find( nodeID + nodeShift );
1819         }
1820         //ma.ordre = ++order;
1821         partGroupe->mailles[i] = imed.insert(ma);
1822       }
1823       _SubPart subPart( partNum, "block" );
1824       subPart.myNbCells    = nbElems;
1825       subPart.myNbNodes    = NumberOfNodes;
1826       subPart.myFirstNode  = firstNode;
1827       subPart.myCellGroupIndex = imed.groupes.size();
1828       imed.addSubPart( subPart );
1829
1830       unstruct->removeReference();
1831     }
1832     else
1833     {
1834       throw MEDEXCEPTION
1835         ( LOCALIZED( STRING(LOC) << "Unexpected word: " << word <<
1836                      " in " << getDataFileName()));
1837     }
1838   } // while ( !geoFile.eof() )
1839
1840   if ( nbParts > 1 )
1841     imed.mergeNodesAndElements(TOLERANCE);
1842
1843   END_OF_MED(LOC);
1844 }
1845
1846 //================================================================================
1847 /*!
1848  * \brief Read mesh in Gold Binary format
1849  */
1850 //================================================================================
1851
1852 void ENSIGHT_MESH_RDONLY_DRIVER::readGoldBinary(_InterMed & imed)
1853 {
1854   const char * LOC ="ENSIGHT_MESH_RDONLY_DRIVER::readGoldBinary() : ";
1855   BEGIN_OF_MED(LOC);
1856
1857   _BinaryFileReader geoFile( getDataFileName() );
1858
1859   // check if swapping bytes needed
1860   try {
1861     countPartsBinary( geoFile, isSingleFileMode());
1862   }
1863   catch (...) {
1864     geoFile.swapBytes();
1865     geoFile.rewind();
1866   }
1867   if ( getIndexInDataFile() <= 1 )
1868     geoFile.rewind();
1869   if ( geoFile.getPosition() == 0 ) {
1870     TStrOwner format( geoFile.getLine() ); // "C|Fortran Binary"
1871     if ( !contains( "C Binary", format )) {
1872       if ( contains( "Fortran Binary", format ))
1873         throw(MEDEXCEPTION(STRING(LOC) << "Fortran Binary format not supported"));
1874       else
1875         throw(MEDEXCEPTION(STRING(LOC) << "unexpected line in " << getDataFileName()
1876                            << "\n" << format.myValues));
1877     }
1878   }
1879   if ( isSingleFileMode() ) {
1880     // one time step may be skipped by countPartsBinary
1881     int curTimeStep = geoFile.getPosition() ? 2 : 1 ;
1882     while ( curTimeStep < getIndexInDataFile() ) {
1883       countPartsBinary( geoFile, true ); // skip time step
1884       curTimeStep++;
1885     }
1886     while (1) {
1887       TStrOwner line( geoFile.getLine() );
1888       if ( isTimeStepBeginning( line.myValues ))
1889         break;
1890     }
1891   }
1892   // ----------------------
1893   // Read mesh description
1894   // ----------------------
1895   {
1896     TStrOwner descriptionLine1 ( geoFile.getLine() );
1897     TStrOwner descriptionLine2 ( geoFile.getLine() );
1898
1899     // find out if the file was created by MED driver
1900     _isMadeByMed = contains( theDescriptionPrefix, descriptionLine1 );
1901
1902     if ( _isMadeByMed )
1903       _ptrMesh->setDescription( descriptionLine2.myValues );
1904     else
1905       _ptrMesh->setDescription( string(descriptionLine1) + descriptionLine2.myValues );
1906   }
1907
1908   // ----------------------------------------
1909   // Find out presence of node/elem numbers
1910   // ----------------------------------------
1911
1912   // EnSight User Manual (for v8) says:
1913   //    You do not have to assign node IDs. If you do, the element connectivities are
1914   //    based on the node numbers. If you let EnSight assign the node IDs, the nodes
1915   //    are considered to be sequential starting at node 1, and element connectivity is
1916   //    done accordingly. If node IDs are set to off, they are numbered internally;
1917   //    however, you will not be able to display or query on them. If you have node
1918   //    IDs in your data, you can have EnSight ignore them by specifying "node id
1919   //    ignore." Using this option may reduce some of the memory taken up by the
1920   //    Client and Server, but display and query on the nodes will not be available.
1921
1922   // read "node|element id <off|given|assign|ignore>"
1923   bool haveNodeIds, haveElemIds;
1924   {
1925     TStrOwner nodeIds( geoFile.getLine() );
1926     TStrOwner elemIds( geoFile.getLine() );
1927
1928     haveNodeIds = ( contains( "given", nodeIds ) || contains( "ignore", nodeIds ) );
1929     haveElemIds = ( contains( "given", elemIds ) || contains( "ignore", elemIds ) );
1930   }
1931
1932   typedef map<int,_noeud>::iterator INoeud;
1933   map<int,_noeud> & points = imed.points;
1934   INoeud firstNode;
1935
1936   _groupe * partGroupe = 0;
1937   int partNum = 0, nbParts = 0;
1938
1939   TFltOwner extents(0); // extents: xmin xmax ymin ymax zmin zmax
1940
1941   while ( !geoFile.eof() )
1942   {
1943     TStrOwner line( geoFile.getLine() );
1944     if ( isSingleFileMode() && isTimeStepEnd( line.myValues ))
1945       break;
1946     string word, restLine;
1947     TStrTool::split( line.myValues, word, restLine );
1948
1949     const TEnSightElemType & partType = getEnSightType( word );
1950     if ( partType._medType != MED_ALL_ELEMENTS )
1951     {
1952       //  Unstructured element type encounters
1953       // --------------------------------------
1954       int      nbElems = *TIntOwner( geoFile.getInt(1) ); // ne
1955       int  nbElemNodes = partType._medType % 100;
1956       bool     isGhost = isGhostType( word );
1957       int    nodeShift = points.empty() ? 0 : points.rbegin()->first;
1958
1959       // read element ids
1960       TIntOwner elemIds( haveElemIds ? geoFile.getInt( nbElems ): 0 ); // id_e
1961
1962       if ( isGhost ) { // do not store ghost elements (?)
1963         int nbNodes = nbElems * nbElemNodes;
1964         if ( partType._name == "nsided" ) // polygons
1965         {
1966           TIntOwner nbNodesInFace( geoFile.getInt( nbElems ));
1967           nbNodes = std::accumulate( nbNodesInFace.myValues, nbNodesInFace.myValues + nbElems, 0 );
1968         }
1969         else if ( partType._name == "nfaced" ) // polyhedrons
1970         {
1971           TIntOwner nbElemFaces( geoFile.getInt( nbElems ));
1972           int nbFaces = accumulate( nbElemFaces.myValues, nbElemFaces.myValues + nbElems, 0 );
1973           TIntOwner nbNodesInFace( geoFile.getInt( nbFaces ));
1974           nbNodes = std::accumulate( nbNodesInFace.myValues, nbNodesInFace.myValues + nbElems, 0 );
1975         }
1976         geoFile.skip( nbNodes * sizeof(int) );
1977         continue;
1978       }
1979
1980       // add a group corresponding to subPart (geoType)
1981       imed.groupes.push_back(_groupe());
1982       _groupe & groupe = imed.groupes.back();
1983       groupe.mailles.resize( nbElems );
1984
1985       // find out if "coordinates" has already been encountered
1986       _SubPartDesc coordDesc( partNum , "coordinates");
1987       map< _SubPartDesc, _SubPart >::iterator descPart =
1988         imed._subPartDescribed.find( coordDesc );
1989       bool haveCoords = ( descPart != imed._subPartDescribed.end() );
1990       if ( haveCoords ) {
1991         firstNode = descPart->second.myFirstNode;
1992         nodeShift -= descPart->second.myNbNodes;
1993       }
1994
1995       // read poly element data
1996       bool isPoly = ( !nbElemNodes );
1997       int nbNodes = 0;
1998       TIntOwner nbElemNodesVec(0), nbElemFaces(0), nbFaceNodes(0);
1999       if ( partType._name == "nsided" ) // polygons
2000       {
2001         nbElemNodesVec.myValues = geoFile.getInt( nbElems ); // npi
2002         nbNodes = accumulate( nbElemNodesVec.myValues, nbElemNodesVec.myValues + nbElems, 0 );
2003       }
2004       else if ( partType._name == "nfaced" ) // polyhedrons
2005       {
2006         nbElemFaces.myValues = geoFile.getInt( nbElems ); // nf_ei
2007         int totalNbFaces = accumulate( nbElemFaces.myValues, nbElemFaces.myValues + nbElems, 0 );
2008
2009         nbFaceNodes.myValues = geoFile.getInt( totalNbFaces ); // np(f_ei)
2010         // calculate nb of nodes in each polyhedron
2011         int* nbEN = nbElemNodesVec.myValues = new int[ nbElems ];
2012         const int *nbFN = nbFaceNodes, *nbEF = nbElemFaces, *nbEND = nbEN + nbElems;
2013         for ( ; nbEN < nbEND; ++nbEN, ++nbEF ) {
2014           nbNodes += *nbEN = accumulate( nbFN, nbFN + *nbEF, 0 );
2015           nbFN    += *nbEF;
2016         }
2017       }
2018       else // standard types
2019       {
2020         nbElemNodesVec.myValues = new int[ 1 ];
2021         nbElemNodesVec[ 0 ] = nbElemNodes;
2022         nbNodes = nbElems * nbElemNodes;
2023       }
2024       // iterator returning nbElemNodes for standard elems and
2025       // next value from nbElemNodesVec for poly elements
2026       _ValueIterator<int> nbElemNodesIt( nbElemNodesVec, isPoly ? 1 : 0);
2027
2028       // iterator returning values form partType._medIndex for standard elems
2029       // and node index (n) for poly elements
2030       int n;
2031       _ValueIterator<int> medIndexIt( isPoly ? &n : &partType._medIndex[0],
2032                                       isPoly ? 0  : 1);
2033       // read connectivity
2034       _maille ma( partType._medType, nbElemNodes );
2035       ma.sommets.resize( nbElemNodes );
2036       TIntOwner connectivity( geoFile.getInt( nbNodes )); // nn_ei
2037       int* nodeID = connectivity;
2038       INoeud node;
2039       for ( int i = 0; i < nbElems; ++i ) {
2040         _ValueIterator<int> medIndex = medIndexIt;
2041         nbElemNodes = nbElemNodesIt.next();
2042         if ((int) ma.sommets.size() != nbElemNodes )
2043           ma.sommets.resize( nbElemNodes );
2044         for ( n = 0; n < nbElemNodes; ++n, ++nodeID ) {
2045           if ( haveCoords )
2046             node = points.find( *nodeID + nodeShift );
2047           else
2048             node = points.insert( make_pair( *nodeID + nodeShift, _noeud())).first;
2049           ma.sommets[ medIndex.next() ] = node;
2050         }
2051         if ( haveElemIds )
2052           ma.setOrdre( elemIds[ i ] );
2053         groupe.mailles[i] = imed.insert(ma);
2054       }
2055       // store nb nodes in polyhedron faces
2056       if ( nbFaceNodes.myValues ) {
2057         const int* nbFaceNodesPtr = nbFaceNodes.myValues;
2058         for ( int i = 0; i < nbElems; ++i ) {
2059           vector<int> & nbNodesByFace = imed.polyherdalNbFaceNodes[ &groupe.maille( i ) ];
2060           nbNodesByFace.assign( nbFaceNodesPtr, nbFaceNodesPtr + nbElemFaces[ i ] );
2061           nbFaceNodesPtr += nbElemFaces[ i ];
2062         }
2063       }
2064       // create subPart for "coordinates"
2065       if ( !haveCoords ) {
2066         _SubPart & coordSubPart = imed._subPartDescribed[ coordDesc ];
2067         coordSubPart.myFirstNode = points.insert( make_pair( 1 + nodeShift, _noeud())).first;
2068       }
2069       // add subPart group to part group
2070       int groupeIndex = imed.groupes.size();
2071       partGroupe->groupes.push_back( groupeIndex );
2072
2073       // create subPart
2074       _SubPart subPart( partNum, partType._name );
2075       subPart.myNbCells = nbElems;
2076       subPart.myCellGroupIndex = groupeIndex;
2077       imed.addSubPart( subPart );
2078     }
2079     else if ( word == "coordinates" )
2080     {
2081       // Local node coordinates of a part
2082       // ------------------------------------
2083       int nbNodes = *TIntOwner( geoFile.getInt(1) ); // nn
2084
2085       // read node ids
2086       TIntOwner nodeIds(0);
2087       if ( haveNodeIds )
2088         nodeIds.myValues = geoFile.getInt( nbNodes ); // id_n
2089
2090       // find out if "coordinates" has already been add at reading connectivity
2091       _SubPartDesc coordDesc( partNum , "coordinates");
2092       map< _SubPartDesc, _SubPart >::iterator descPart =
2093         imed._subPartDescribed.find( coordDesc );
2094       bool haveCoords = ( descPart != imed._subPartDescribed.end() );
2095
2096       if ( haveCoords ) {
2097         // check that all nodes have been added
2098         firstNode = descPart->second.myFirstNode;
2099         descPart->second.myNbNodes = nbNodes;
2100         INoeud inoeud = firstNode, inoEnd = points.end();
2101         int id = inoeud->first, idEnd = id + nbNodes;
2102         for ( ; id < idEnd; ++id ) {
2103           if ( inoeud == inoEnd || inoeud->first > id ) {
2104             INoeud in = points.insert( inoeud, make_pair( id, _noeud() ));
2105             in->second.number = id;
2106           } else {
2107             ++inoeud;
2108           }
2109         }
2110       }
2111       else {
2112         // add nodes
2113         int nodeShift = points.empty() ? 0 : points.rbegin()->first;
2114         for ( int iNode = 1; iNode <= nbNodes; ++iNode ) {
2115           INoeud inoeud = points.insert( points.end(), make_pair( iNode + nodeShift, _noeud()));
2116           inoeud->second.number = inoeud->first;
2117         }
2118         firstNode = points.find( 1 + nodeShift );
2119         // create "coordinates" subPart
2120         _SubPart & subPart  = imed._subPartDescribed[ coordDesc ];
2121         subPart.myNbNodes   = nbNodes;
2122         subPart.myFirstNode = firstNode;
2123       }
2124
2125       // read coordinates in no interlace mode
2126       TFltOwner noInterlaceCoords( geoFile.getFlt( nbNodes * SPACE_DIM ));
2127       float* x = noInterlaceCoords;
2128       float* y = x + nbNodes;
2129       float* z = y + nbNodes;
2130       INoeud endNode = points.end();
2131       for ( INoeud in = firstNode; in != endNode; ++in ) {
2132         _noeud & node = in->second;
2133         node.coord.resize( SPACE_DIM );
2134         node.coord[ 0 ] = *x++;
2135         node.coord[ 1 ] = *y++;
2136         node.coord[ 2 ] = *z++;
2137       }
2138     }
2139     else if ( word == "part" )
2140     {
2141       // Another part encounters
2142       // -----------------------
2143       partNum = *TIntOwner( geoFile.getInt(1) );
2144       nbParts++;
2145
2146       string partName( TStrOwner( geoFile.getLine() ));
2147       if ( partName.empty() )
2148         partName = "Part_" + restLine;
2149
2150       if (int( imed.groupes.capacity() - imed.groupes.size()) < theMaxNbTypes )
2151         imed.groupes.reserve( size_t( 1.5 * imed.groupes.size() ));
2152       imed.groupes.push_back(_groupe());
2153       partGroupe = & imed.groupes.back();
2154       partGroupe->nom = partName;
2155       partGroupe->groupes.reserve( theMaxNbTypes );
2156     }
2157     else if ( word == "block" )
2158     {
2159       // Structured type
2160       // ------------------
2161       bool rectilinear = ( restLine.find( "rectilinear" ) != restLine.npos );
2162       bool uniform     = ( restLine.find( "uniform" )     != restLine.npos );
2163       bool curvilinear = ( !rectilinear && !uniform );
2164       bool iblanked    = ( restLine.find( "iblanked" )    != restLine.npos );
2165       bool with_ghost  = ( restLine.find( "with_ghost" )  != restLine.npos );
2166       bool range       = ( restLine.find( "range" )       != restLine.npos );
2167
2168       // dimension
2169       TIntOwner ijk( geoFile.getInt(3) );
2170       int I = ijk[0];
2171       int J = ijk[1];
2172       int K = ijk[2];
2173       int NumberOfNodes = I*J*K;
2174       if ( !NumberOfNodes ) continue;
2175
2176       // range
2177       if ( range ) {
2178         TIntOwner ijkRange( geoFile.getInt( 6 ));// imin imax jmin jmax kmin kmax
2179         I = ijkRange[1]-ijkRange[0]+1;
2180         J = ijkRange[3]-ijkRange[2]+1;
2181         K = ijkRange[5]-ijkRange[4]+1;
2182         NumberOfNodes = I*J*K;
2183       }
2184       // add nodes
2185       int nodeShift = points.empty() ? 0 : points.rbegin()->first;
2186       for ( int iNode = 1; iNode <= NumberOfNodes; ++iNode ) {
2187         INoeud inoeud = points.insert( points.end(), make_pair( iNode + nodeShift, _noeud()));
2188         _noeud & node = inoeud->second;
2189         node.number   = inoeud->first;
2190         node.coord.resize( SPACE_DIM );
2191       }
2192       INoeud firstNode = points.find( nodeShift + 1 );
2193       INoeud endNode   = points.end();
2194
2195       GRID grid; // calculator of unstructured data
2196       grid._iArrayLength   = I;
2197       grid._jArrayLength   = J;
2198       grid._kArrayLength   = K;
2199       grid._spaceDimension = SPACE_DIM;
2200       if ( J < 2 ) { grid._spaceDimension--; grid._jArrayLength = 0; }
2201       if ( K < 2 ) { grid._spaceDimension--; grid._kArrayLength = 0; }
2202       int nbElems = grid.getNumberOfElements(MED_CELL, MED_ALL_ELEMENTS);
2203
2204       if ( curvilinear ) // read coordinates for all nodes
2205       {
2206         TFltOwner noInterlaceCoords( geoFile.getFlt( NumberOfNodes * SPACE_DIM ));
2207         float* x = noInterlaceCoords;
2208         float* y = x + NumberOfNodes;
2209         float* z = y + NumberOfNodes;
2210         for ( INoeud in = firstNode; in != endNode; ++in ) {
2211           _noeud & node = in->second;
2212           node.coord.resize( SPACE_DIM );
2213           node.coord[ 0 ] = *x++;
2214           node.coord[ 1 ] = *y++;
2215           node.coord[ 2 ] = *z++;
2216         }
2217         grid._gridType   = MED_BODY_FITTED;
2218         grid._coordinate = new COORDINATE(3, 1, 0); // to avoid exception in convertInMESH()
2219       }
2220       else if ( rectilinear ) // read delta vectors with non-regular spacing
2221       {
2222         grid._iArray = (double*)geoFile.convertReals<double>( I );
2223         grid._jArray = (double*)geoFile.convertReals<double>( J );
2224         grid._kArray = (double*)geoFile.convertReals<double>( K );
2225         grid._gridType = MED_CARTESIAN;
2226       }
2227       else // uniform: read grid origine and delta vectors for regular spacing grid
2228       {
2229         TFltOwner xyzOrigin( geoFile.getFlt( 3 ));
2230         TFltOwner xyzDelta ( geoFile.getFlt( 3 ));
2231         // compute full delta vectors
2232         grid._iArray = new double[ I ];
2233         grid._jArray = new double[ J ];
2234         grid._kArray = new double[ K ];
2235         double* coors[SPACE_DIM] = { grid._iArray, grid._jArray, grid._kArray };
2236         int     size [SPACE_DIM] = { I, J, K };
2237         for ( int j = 0; j < SPACE_DIM; ++j ) {
2238           double* coo    = coors[ j ];
2239           double* cooEnd = coo + size[ j ];
2240           coo[0]         = xyzOrigin[ j ];
2241           while ( ++coo < cooEnd )
2242             *coo = coo[-1] + xyzDelta[ j ];
2243         }
2244         grid._gridType = MED_CARTESIAN;
2245       }
2246
2247       // iblanks
2248       if ( iblanked )
2249         geoFile.skip( NumberOfNodes * sizeof(int) );
2250       // ghosts
2251       if ( with_ghost ) {
2252         TStrOwner( geoFile.getLine() ); // "ghost_flags"
2253         geoFile.skip( nbElems * sizeof(int) );
2254       }
2255       // node ids
2256       if ( haveNodeIds && !geoFile.eof() ) {
2257         TStrOwner nextLine( geoFile.getLine() ); // "node_ids"
2258         if ( contains( "node_ids", nextLine ) )
2259           geoFile.skip( NumberOfNodes * sizeof(int) );
2260         else
2261           geoFile.skip( -MAX_LINE_LENGTH );
2262       }
2263       // element ids
2264       TIntOwner elemIdOwner(0);
2265       _ValueIterator<int> elemIds;
2266       if ( haveElemIds && !geoFile.eof() ) {
2267         TStrOwner nextLine( geoFile.getLine() ); // "element_ids"
2268         if ( contains( "element_ids", nextLine ) ) {
2269           elemIdOwner.myValues = geoFile.getInt( nbElems );
2270           elemIds = _ValueIterator<int>( elemIdOwner, 1);
2271         } else {
2272           geoFile.skip( -MAX_LINE_LENGTH );
2273         }
2274       }
2275
2276       // let GRID compute all coordinates and connectivity
2277       const MESH* unstruct = grid.convertInMESH();
2278       if ( !curvilinear )
2279       {
2280         const double * coo = unstruct->getCoordinates(MED_FULL_INTERLACE);
2281         typedef _ValueIterator< double > TCoordIt;
2282         TCoordIt xCoo( coo+0, grid._spaceDimension);
2283         TCoordIt yCoo( coo+1, grid._spaceDimension);
2284         TCoordIt zCoo( coo+2, grid._spaceDimension);
2285         if ( grid._spaceDimension < 3 ) zCoo = TCoordIt( grid._kArray, 0 );
2286         if ( grid._spaceDimension < 2 ) yCoo = TCoordIt( grid._jArray, 0 );
2287         for ( INoeud in = firstNode; in != endNode; ++in ) {
2288           _noeud& node = in->second;
2289           node.coord[ 0 ] = xCoo.next();
2290           node.coord[ 1 ] = yCoo.next();
2291           node.coord[ 2 ] = zCoo.next();
2292         }
2293       }
2294
2295       // store connectivity
2296       const int * conn = unstruct->getConnectivity( MED_NODAL, MED_CELL, MED_ALL_ELEMENTS );
2297       medGeometryElement elemType = grid.getElementType( MED_CELL, 1 );
2298       int  nbElemNodes = elemType % 100;
2299
2300       partGroupe->mailles.resize( nbElems );
2301       _maille ma( elemType, nbElemNodes );
2302       ma.sommets.resize( nbElemNodes );
2303
2304       for ( int i = 0, nIndex = 0; i < nbElems; ++i ) {
2305         for ( int n = 0; n < nbElemNodes; ++n ) {
2306           int nodeID = conn[ nIndex++ ];
2307           ma.sommets[n] = points.find( nodeID + nodeShift );
2308         }
2309         ma.setOrdre( elemIds.next() );
2310         partGroupe->mailles[i] = imed.insert(ma);
2311       }
2312
2313       _SubPart subPart( partNum, "block" );
2314       subPart.myNbCells    = nbElems;
2315       subPart.myNbNodes    = NumberOfNodes;
2316       subPart.myFirstNode  = firstNode;
2317       subPart.myCellGroupIndex = imed.groupes.size();
2318       imed.addSubPart( subPart );
2319
2320       unstruct->removeReference();
2321     }
2322     else if ( word == "extents" )
2323     {
2324       extents.myValues = geoFile.getFlt( 6 );
2325     }
2326     else
2327     {
2328       throw MED_EXCEPTION
2329         ( LOCALIZED( STRING(LOC) << "Unexpected word: " << word <<
2330                      " in " << getDataFileName()));
2331     }
2332   } // while ( !geoFile.eof() )
2333
2334   if ( nbParts > 1 )
2335     imed.mergeNodesAndElements(TOLERANCE);
2336
2337   END_OF_MED(LOC);
2338 }
2339
2340 //================================================================================
2341 /*!
2342  * \brief Read mesh in Ensight6 ASCII format
2343  */
2344 //================================================================================
2345
2346 void ENSIGHT_MESH_RDONLY_DRIVER::read6ASCII(_InterMed & imed)
2347 {
2348   const char * LOC ="ENSIGHT_MESH_RDONLY_DRIVER::read6ASCII() : ";
2349   BEGIN_OF_MED(LOC);
2350
2351   _ASCIIFileReader geoFile( getDataFileName() );
2352
2353   if ( isSingleFileMode() ) {
2354     int curTimeStep = 1;
2355     while ( curTimeStep < getIndexInDataFile() ) {
2356       while ( !geoFile.isTimeStepEnd())
2357         geoFile.getLine();
2358       curTimeStep++;
2359     }
2360     while ( !geoFile.isTimeStepBeginning() )
2361       geoFile.getLine();
2362   }
2363   // ----------------------
2364   // Read mesh description
2365   // ----------------------
2366   {
2367     string descriptionLine1 = geoFile.getLine();
2368     string descriptionLine2 = geoFile.getLine();
2369
2370     // find out if the file was created by MED driver
2371     int prefixSize = strlen( theDescriptionPrefix );
2372     _isMadeByMed = ( descriptionLine1.substr(0, prefixSize ) == theDescriptionPrefix );
2373
2374     if ( _isMadeByMed )
2375       descriptionLine1 = descriptionLine1.substr( prefixSize );
2376     _ptrMesh->setDescription( descriptionLine1 + descriptionLine2 );
2377   }
2378
2379   // ----------------------------------------
2380   // Find out presence of node/elem numbers
2381   // ----------------------------------------
2382
2383   // EnSight User Manual (for v8) says:
2384   //    You do not have to assign node IDs. If you do, the element connectivities are
2385   //    based on the node numbers. If you let EnSight assign the node IDs, the nodes
2386   //    are considered to be sequential starting at node 1, and element connectivity is
2387   //    done accordingly. If node IDs are set to off, they are numbered internally;
2388   //    however, you will not be able to display or query on them. If you have node
2389   //    IDs in your data, you can have EnSight ignore them by specifying "node id
2390   //    ignore." Using this option may reduce some of the memory taken up by the
2391   //    Client and Server, but display and query on the nodes will not be available.
2392
2393   // read "node|element id <off|given|assign|ignore>"
2394   geoFile.getWord(); geoFile.getWord();
2395   string nodeIds = geoFile.getWord();
2396   geoFile.getWord(); geoFile.getWord();
2397   string elemIds = geoFile.getWord();
2398
2399   bool haveNodeIds = ( nodeIds == "given" || nodeIds == "ignore" );
2400   bool haveElemIds = ( elemIds == "given" || elemIds == "ignore" );
2401
2402   map<int,_noeud> & points = imed.points;
2403   typedef map<int,_noeud>::iterator INoeud;
2404
2405   int haveStructuredParts = 0, haveUnstructuredParts = 0;
2406
2407   _groupe * partGroupe = 0;
2408   int       partNum = 0;
2409
2410   while ( !geoFile.isTimeStepEnd() )
2411   {
2412     string word, restLine, line = geoFile.getLine();
2413     TStrTool::split( line, word, restLine );
2414
2415     const TEnSightElemType & partType = getEnSightType( word );
2416     if ( !partType._medIndex.empty() )
2417     {
2418       //  Unstructured element type encounters
2419       // --------------------------------------
2420       int  nbElemNodes = partType._medType % 100;
2421       int      nbElems = geoFile.getInt();
2422       if ( nbElems > 0 )
2423         haveUnstructuredParts++;
2424
2425       imed.groupes.push_back(_groupe());
2426       _groupe & groupe = imed.groupes.back();
2427       groupe.mailles.resize( nbElems );
2428
2429       // read connectivity
2430       _maille ma( partType._medType, nbElemNodes );
2431       ma.sommets.resize( nbElemNodes );
2432       INoeud node;
2433       for ( int i = 0; i < nbElems; ++i ) {
2434         if ( haveElemIds )
2435           geoFile.getInt();
2436         for ( int n = 0; n < nbElemNodes; ++n ) {
2437           int nodeID = geoFile.getInt();
2438           ma.sommets[ partType._medIndex[n] ] = points.find( nodeID );
2439         }
2440         //ma.ordre = ++order;
2441         groupe.mailles[i] = imed.insert(ma);
2442       }
2443
2444       int groupeIndex = imed.groupes.size();
2445       partGroupe->groupes.push_back( groupeIndex );
2446
2447       _SubPart subPart( partNum, partType._name );
2448       subPart.myNbCells = nbElems;
2449       subPart.myCellGroupIndex = groupeIndex;
2450       imed.addSubPart( subPart );
2451     }
2452     else if ( word == "part" )
2453     {
2454       // Another part encounters
2455       // -----------------------
2456       partNum = atoi( restLine.c_str() );
2457
2458       string partName = geoFile.getLine();
2459       if ( partName.empty() )
2460         partName = "Part_" + restLine;
2461
2462       if (int( imed.groupes.capacity() - imed.groupes.size()) < theMaxNbTypes )
2463         imed.groupes.reserve( size_t( 1.5 * imed.groupes.size() ));
2464       imed.groupes.push_back(_groupe());
2465       partGroupe = & imed.groupes.back();
2466       partGroupe->nom = partName;
2467       partGroupe->groupes.reserve( theMaxNbTypes );
2468     }
2469     else if ( word == "block" )
2470     {
2471       // Structured type
2472       // ------------------
2473       bool iblanked  = ( restLine == "iblanked" );
2474
2475       // dimension
2476       int I = geoFile.getInt();
2477       int J = geoFile.getInt();
2478       int K = geoFile.getInt();
2479       int NumberOfNodes = I*J*K;
2480       if ( !NumberOfNodes ) continue;
2481       haveStructuredParts++;
2482
2483       // add nodes
2484       int nodeShift = points.empty() ? 0 : points.rbegin()->first;
2485       for ( int iNode = 1; iNode <= NumberOfNodes; ++iNode ) {
2486         INoeud inoeud = points.insert( points.end(), make_pair( iNode + nodeShift, _noeud()));
2487         _noeud & node = inoeud->second;
2488         node.number   = inoeud->first;
2489         node.coord.resize( SPACE_DIM );
2490       }
2491       // read coordinates
2492       INoeud firstNode = points.find( nodeShift + 1 );
2493       INoeud endNode   = points.end();
2494       for ( int j = 0; j < SPACE_DIM; ++j ) {
2495         for ( INoeud in = firstNode; in != endNode; ++in ) {
2496           _noeud & node = in->second;
2497           node.coord[ j ] = geoFile.getReal();
2498         }
2499       }
2500       // iblanks
2501       if ( iblanked )
2502         geoFile.skip(NumberOfNodes, /*nbPerLine =*/ 10, INT_WIDTH_6);
2503
2504       // let GRID calculate connectivity
2505       GRID grid;
2506       grid._iArrayLength  = I;
2507       grid._jArrayLength  = J;
2508       grid._kArrayLength  = K;
2509       grid._gridType      = MED_BODY_FITTED;
2510       grid._spaceDimension= SPACE_DIM;
2511       grid._coordinate    = new COORDINATE(3, 1, 0); // to avoid exception in convertInMESH()
2512       if ( J < 2 ) { grid._spaceDimension--; grid._jArrayLength = 0; }
2513       if ( K < 2 ) { grid._spaceDimension--; grid._kArrayLength = 0; }
2514
2515       const MESH* unstruct = grid.convertInMESH();
2516       const int * conn = unstruct->getConnectivity( MED_NODAL, MED_CELL, MED_ALL_ELEMENTS );
2517       medGeometryElement elemType = grid.getElementType( MED_CELL, 1 );
2518       int  nbElemNodes = elemType % 100;
2519       int      nbElems = grid.getNumberOfElements(MED_CELL, MED_ALL_ELEMENTS);
2520
2521       partGroupe->mailles.resize( nbElems );
2522       _maille ma( elemType, nbElemNodes );
2523       ma.sommets.resize( nbElemNodes );
2524
2525       for ( int i = 0, nIndex = 0; i < nbElems; ++i ) {
2526         for ( int n = 0; n < nbElemNodes; ++n ) {
2527           int nodeID = conn[ nIndex++ ];
2528           ma.sommets[n] = points.find( nodeID + nodeShift );
2529         }
2530         //ma.ordre = ++order;
2531         partGroupe->mailles[i] = imed.insert(ma);
2532       }
2533
2534       _SubPart subPart( partNum, "block" );
2535       subPart.myNbCells    = nbElems;
2536       subPart.myNbNodes    = NumberOfNodes;
2537       subPart.myFirstNode  = firstNode;
2538       subPart.myCellGroupIndex = imed.groupes.size();
2539       imed.addSubPart( subPart );
2540
2541       unstruct->removeReference();
2542     }
2543     else if ( word == "coordinates" )
2544     {
2545       // ------------------------------------
2546       // Unstructured global node coordinates
2547       // ------------------------------------
2548       int nbNodes = geoFile.getInt();
2549
2550       cout << "-> loading coordinates of " << nbNodes << " nodes " << endl ;
2551
2552       INoeud inoeud;
2553       for ( int i=0 ; i < nbNodes ; i++ )
2554       {
2555         if ( haveNodeIds ) {
2556           int nodeID = geoFile.getInt();
2557           inoeud = points.insert( make_pair( nodeID, _noeud() )).first;
2558           inoeud->second.number = nodeID;
2559         }
2560         else {
2561           int nodeID = i + 1;
2562           inoeud = points.insert( points.end(), make_pair( nodeID, _noeud()));
2563           inoeud->second.number = nodeID;
2564         }
2565         _noeud & node = inoeud->second;
2566         node.coord.resize( SPACE_DIM );
2567         node.coord[ 0 ] = geoFile.getReal();
2568         node.coord[ 1 ] = geoFile.getReal();
2569         node.coord[ 2 ] = geoFile.getReal();
2570       }
2571
2572       _SubPartDesc cooDesc = _SubPartDesc::globalCoordDesc();
2573       _SubPart subPart( cooDesc.partNumber(), cooDesc.typeName() );
2574       subPart.myNbNodes    = nbNodes;
2575       subPart.myFirstNode  = points.begin();
2576       imed.addSubPart( subPart );
2577     }
2578     else
2579     {
2580       throw MED_EXCEPTION
2581         ( LOCALIZED( STRING(LOC) << "Unexpected word: " << word <<
2582                      " in " << getDataFileName()));
2583     }
2584   } // while ( !geoFile.eof() )
2585
2586   if ( ( haveStructuredParts && haveUnstructuredParts ) || haveStructuredParts > 1 )
2587     imed.mergeNodesAndElements(TOLERANCE);
2588
2589   END_OF_MED(LOC);
2590 }
2591
2592 //================================================================================
2593 /*!
2594  * \brief Read mesh in Ensight6 ASCII format
2595  */
2596 //================================================================================
2597
2598 void ENSIGHT_MESH_RDONLY_DRIVER::read6Binary(_InterMed & imed)
2599 {
2600   const char * LOC ="ENSIGHT_MESH_RDONLY_DRIVER::read6Binary() : ";
2601   BEGIN_OF_MED(LOC);
2602
2603   _BinaryFileReader geoFile( getDataFileName() );
2604
2605   // check if swapping bytes needed
2606   try {
2607     countPartsBinary( geoFile, isSingleFileMode());
2608   }
2609   catch (...) {
2610     geoFile.swapBytes();
2611     geoFile.rewind();
2612   }
2613   if ( getIndexInDataFile() <= 1 )
2614     geoFile.rewind();
2615   if ( geoFile.getPosition() == 0 ) {
2616     TStrOwner format( geoFile.getLine() ); // "C|Fortran Binary"
2617     if ( !contains( "C Binary", format )) {
2618       if ( contains( "Fortran Binary", format ))
2619         throw(MEDEXCEPTION(STRING(LOC) << "Fortran Binary format not supported"));
2620       else
2621         throw(MEDEXCEPTION(STRING(LOC) << "unexpected line in " << getDataFileName()
2622                            << "\n" << format.myValues));
2623     }
2624   }
2625
2626   if ( isSingleFileMode() ) {
2627     // one time step may be skipped by countPartsBinary
2628     int curTimeStep = geoFile.getPosition() ? 2 : 1 ;
2629     while ( curTimeStep < getIndexInDataFile() ) {
2630       countPartsBinary( geoFile, true ); // skip time step
2631       curTimeStep++;
2632     }
2633     while (1) {
2634       TStrOwner line( geoFile.getLine() );
2635       if ( isTimeStepBeginning( line.myValues ))
2636         break;
2637     }
2638   }
2639   // ----------------------
2640   // Read mesh description
2641   // ----------------------
2642   {
2643     TStrOwner descriptionLine1( geoFile.getLine() );
2644     TStrOwner descriptionLine2( geoFile.getLine() );
2645
2646     // find out if the file was created by MED driver
2647     _isMadeByMed = contains( theDescriptionPrefix, descriptionLine1 );
2648
2649     if ( _isMadeByMed )
2650       _ptrMesh->setDescription( descriptionLine2.myValues );
2651     else
2652       _ptrMesh->setDescription( string(descriptionLine1) + descriptionLine2.myValues );
2653   }
2654
2655   // ----------------------------------------
2656   // Find out presence of node/elem numbers
2657   // ----------------------------------------
2658
2659   // EnSight User Manual (for v8) says:
2660   //    You do not have to assign node IDs. If you do, the element connectivities are
2661   //    based on the node numbers. If you let EnSight assign the node IDs, the nodes
2662   //    are considered to be sequential starting at node 1, and element connectivity is
2663   //    done accordingly. If node IDs are set to off, they are numbered internally;
2664   //    however, you will not be able to display or query on them. If you have node
2665   //    IDs in your data, you can have EnSight ignore them by specifying "node id
2666   //    ignore." Using this option may reduce some of the memory taken up by the
2667   //    Client and Server, but display and query on the nodes will not be available.
2668
2669   // read "node|element id <off|given|assign|ignore>"
2670   bool haveNodeIds, haveElemIds;
2671   {
2672     TStrOwner nodeIds( geoFile.getLine() );
2673     TStrOwner elemIds( geoFile.getLine() );
2674
2675     haveNodeIds = ( contains( "given", nodeIds ) || contains( "ignore", nodeIds ) );
2676     haveElemIds = ( contains( "given", elemIds ) || contains( "ignore", elemIds ) );
2677   }
2678   map<int,_noeud> & points = imed.points;
2679   typedef map<int,_noeud>::iterator INoeud;
2680
2681   int haveStructuredParts = 0, haveUnstructuredParts = 0;
2682
2683   _groupe * partGroupe = 0;
2684   int       partNum = 0;
2685
2686   while ( !geoFile.eof() )
2687   {
2688     TStrOwner line( geoFile.getLine() );
2689     if ( isSingleFileMode() && isTimeStepEnd( line.myValues ))
2690       break;
2691     string word, restLine;
2692     TStrTool::split( line.myValues, word, restLine );
2693
2694     const TEnSightElemType & partType = getEnSightType( word );
2695     if ( !partType._medIndex.empty() )
2696     {
2697       //  Unstructured element type encounters
2698       // --------------------------------------
2699       int  nbElemNodes = partType._medType % 100;
2700       int      nbElems = *TIntOwner( geoFile.getInt(1) ); // ne
2701       if ( nbElems > 0 )
2702         haveUnstructuredParts++;
2703
2704       TIntOwner numbers(0);
2705       if ( haveElemIds )
2706         numbers.myValues = geoFile.getInt( nbElems ); // id_e
2707
2708       imed.groupes.push_back(_groupe());
2709       _groupe & groupe = imed.groupes.back();
2710       groupe.mailles.resize( nbElems );
2711
2712       // read connectivity
2713       _maille ma( partType._medType, nbElemNodes );
2714       ma.sommets.resize( nbElemNodes );
2715       TIntOwner connectivity( geoFile.getInt( nbElems * nbElemNodes ));
2716       int* nodeID = connectivity;
2717       INoeud node;
2718       for ( int i = 0; i < nbElems; ++i ) {
2719         for ( int n = 0; n < nbElemNodes; ++n, ++nodeID )
2720           ma.sommets[ partType._medIndex[n] ] = points.find( *nodeID );
2721         //ma.ordre = ++order;
2722         groupe.mailles[i] = imed.insert(ma);
2723       }
2724
2725       int groupeIndex = imed.groupes.size();
2726       partGroupe->groupes.push_back( groupeIndex );
2727
2728       _SubPart subPart( partNum, partType._name );
2729       subPart.myNbCells = nbElems;
2730       subPart.myCellGroupIndex = groupeIndex;
2731       imed.addSubPart( subPart );
2732     }
2733     else if ( word == "part" )
2734     {
2735       // Another part encounters
2736       // -----------------------
2737       partNum = atoi( restLine.c_str() );
2738
2739       string partName( TStrOwner( geoFile.getLine() ));
2740       if ( partName.empty() )
2741         partName = "Part_" + restLine;
2742
2743       if (int( imed.groupes.capacity() - imed.groupes.size()) < theMaxNbTypes )
2744         imed.groupes.reserve( size_t( 1.5 * imed.groupes.size() ));
2745       imed.groupes.push_back(_groupe());
2746       partGroupe = & imed.groupes.back();
2747       partGroupe->nom = partName;
2748       partGroupe->groupes.reserve( theMaxNbTypes );
2749     }
2750     else if ( word == "block" )
2751     {
2752       // Structured type
2753       // ------------------
2754       bool iblanked  = ( restLine == "iblanked" );
2755
2756       // dimension
2757       TIntOwner ijk( geoFile.getInt(3) );
2758       int I = ijk[0];
2759       int J = ijk[1];
2760       int K = ijk[2];
2761       int NumberOfNodes = I*J*K;
2762       if ( !NumberOfNodes ) continue;
2763       haveStructuredParts++;
2764
2765       // read coordinates
2766       int nodeShift = points.empty() ? 0 : points.rbegin()->first;
2767       {
2768         TFltOwner noInterlaceCoords( geoFile.getFlt( NumberOfNodes * SPACE_DIM ));
2769         float* x = noInterlaceCoords;
2770         float* y = x + NumberOfNodes;
2771         float* z = y + NumberOfNodes;
2772         for ( int iNode = 1; iNode <= NumberOfNodes; ++iNode ) {
2773           INoeud inoeud = points.insert( points.end(), make_pair( iNode + nodeShift, _noeud()));
2774           _noeud & node = inoeud->second;
2775           node.number   = inoeud->first;
2776           node.coord.resize( SPACE_DIM );
2777           node.coord[0] = *x++;
2778           node.coord[1] = *y++;
2779           node.coord[2] = *z++;
2780         }
2781       }
2782       // iblanks
2783       if ( iblanked )
2784         geoFile.skip(NumberOfNodes * sizeof(int));
2785
2786       // let GRID calculate connectivity
2787       GRID grid;
2788       grid._iArrayLength  = I;
2789       grid._jArrayLength  = J;
2790       grid._kArrayLength  = K;
2791       grid._gridType      = MED_BODY_FITTED;
2792       grid._spaceDimension= SPACE_DIM;
2793       grid._coordinate    = new COORDINATE(3, 1, 0); // to avoid exception in convertInMESH()
2794       if ( J < 2 ) { grid._spaceDimension--; grid._jArrayLength = 0; }
2795       if ( K < 2 ) { grid._spaceDimension--; grid._kArrayLength = 0; }
2796
2797       const MESH* unstruct = grid.convertInMESH();
2798       const int * conn = unstruct->getConnectivity( MED_NODAL, MED_CELL, MED_ALL_ELEMENTS );
2799       medGeometryElement elemType = grid.getElementType( MED_CELL, 1 );
2800       int  nbElemNodes = elemType % 100;
2801       int      nbElems = grid.getNumberOfElements(MED_CELL, MED_ALL_ELEMENTS);
2802
2803       partGroupe->mailles.resize( nbElems );
2804       _maille ma( elemType, nbElemNodes );
2805       ma.sommets.resize( nbElemNodes );
2806
2807       for ( int i = 0, nIndex = 0; i < nbElems; ++i ) {
2808         for ( int n = 0; n < nbElemNodes; ++n ) {
2809           int nodeID = conn[ nIndex++ ];
2810           ma.sommets[n] = points.find( nodeID + nodeShift );
2811         }
2812         //ma.ordre = ++order;
2813         partGroupe->mailles[i] = imed.insert(ma);
2814       }
2815
2816       _SubPart subPart( partNum, "block" );
2817       subPart.myNbCells    = nbElems;
2818       subPart.myNbNodes    = NumberOfNodes;
2819       subPart.myFirstNode  = points.find( nodeShift + 1 );
2820       subPart.myCellGroupIndex = imed.groupes.size();
2821       imed.addSubPart( subPart );
2822
2823       unstruct->removeReference();
2824     }
2825     else if ( word == "coordinates" )
2826     {
2827       // ------------------------------
2828       // Unstructured node coordinates
2829       // ------------------------------
2830       int nbNodes = *TIntOwner( geoFile.getInt(1) );
2831
2832       TIntOwner numbers(0);
2833       if ( haveNodeIds )
2834         numbers.myValues = geoFile.getInt( nbNodes );
2835
2836       TFltOwner fullInterlaceCoords( geoFile.getFlt( nbNodes * SPACE_DIM ));
2837       float* coord = fullInterlaceCoords;
2838
2839       cout << "-> loading coordinates of " << nbNodes << " nodes " << endl ;
2840
2841       INoeud inoeud;
2842       for ( int i=0 ; i < nbNodes ; i++ )
2843       {
2844         if ( haveNodeIds ) {
2845           int nodeID = numbers[ i ];
2846           inoeud = points.insert( make_pair( nodeID, _noeud() )).first;
2847           inoeud->second.number = nodeID;
2848         }
2849         else {
2850           int nodeID = i + 1;
2851           inoeud = points.insert( points.end(), make_pair( nodeID, _noeud()));
2852           inoeud->second.number = nodeID;
2853         }
2854         _noeud & node = inoeud->second;
2855         node.coord.resize( SPACE_DIM );
2856         node.coord[ 0 ] = *coord++;
2857         node.coord[ 1 ] = *coord++;
2858         node.coord[ 2 ] = *coord++;
2859       }
2860
2861       _SubPartDesc cooDesc = _SubPartDesc::globalCoordDesc();
2862       _SubPart subPart( cooDesc.partNumber(), cooDesc.typeName() );
2863       subPart.myNbNodes    = nbNodes;
2864       subPart.myFirstNode  = points.begin();
2865       imed.addSubPart( subPart );
2866     }
2867     else
2868     {
2869       throw MED_EXCEPTION
2870         ( LOCALIZED( STRING(LOC) << "Unexpected word: " << word <<
2871                      " in " << getDataFileName()));
2872     }
2873   } // while ( !geoFile.eof() )
2874
2875   if ( ( haveStructuredParts && haveUnstructuredParts ) || haveStructuredParts > 1 )
2876     imed.mergeNodesAndElements(TOLERANCE);
2877
2878   END_OF_MED(LOC);
2879 }
2880
2881 //================================================================================
2882 /*!
2883  * \brief count number of parts in EnSight geometry file
2884  */
2885 //================================================================================
2886
2887 int ENSIGHT_MESH_RDONLY_DRIVER::countParts(const string& geomFileName,
2888                                            const bool    isSingleFileMode)
2889 {
2890   const char * LOC ="ENSIGHT_MESH_RDONLY_DRIVER::countParts() : ";
2891
2892   int nbParts = 0;
2893   if ( isBinaryDataFile( geomFileName ))
2894   {
2895     _BinaryFileReader geoFile(geomFileName);
2896     // check if swapping bytes needed
2897     try {
2898       return countPartsBinary( geoFile, isSingleFileMode );
2899     }
2900     catch (...) {
2901     }
2902     geoFile.swapBytes();
2903     geoFile.rewind();
2904     nbParts = countPartsBinary( geoFile, isSingleFileMode );
2905   }
2906   else
2907   {
2908     _ASCIIFileReader geoFile(geomFileName);
2909
2910     if ( isSingleFileMode )
2911       while ( !isTimeStepBeginning( geoFile.getLine() ));
2912
2913     geoFile.getLine(); // description line 1
2914     geoFile.getLine(); // description line 2
2915
2916     // read "node|element id <off|given|assign|ignore>"
2917     geoFile.getWord(); geoFile.getWord();
2918     string nodeIds = geoFile.getWord();
2919     geoFile.getWord(); geoFile.getWord();
2920     string elemIds = geoFile.getWord();
2921     bool haveNodeIds = ( nodeIds == "given" || nodeIds == "ignore" );
2922     bool haveElemIds = ( elemIds == "given" || elemIds == "ignore" );
2923
2924     bool isGold = true;
2925     while ( !geoFile.isTimeStepEnd() )
2926     {
2927       string word, restLine, line = geoFile.getLine();
2928       TStrTool::split( line, word, restLine );
2929
2930       const TEnSightElemType & partType = getEnSightType( word );
2931       if ( partType._medType != MED_ALL_ELEMENTS )
2932       {
2933         //  Unstructured element type encounters
2934         // --------------------------------------
2935         int  nbElemNodes = partType._medType % 100;
2936         int      nbElems = geoFile.getInt(); // ne
2937
2938         // element ids
2939         if ( haveElemIds && isGold )
2940           geoFile.skip( nbElems, /*nbPerLine =*/ 1, INT_WIDTH_GOLD ); // id_e
2941
2942         // skip connectivity
2943         int nbNodes = nbElems * nbElemNodes;
2944         if ( partType._name == "nsided" ) // polygons
2945         {
2946           for ( int i = 0; i < nbElems; ++i )
2947             nbNodes += geoFile.getInt();
2948           geoFile.skip( nbNodes * INT_WIDTH_GOLD, /*nbLines = */nbElems );
2949         }
2950         else if ( partType._name == "nfaced" ) // polyhedrons
2951         {
2952           int nbFaces = 0;
2953           for ( int i = 0; i < nbElems; ++i )
2954             nbFaces += geoFile.getInt();
2955           for ( int f = 0; f < nbFaces; ++f )
2956             nbNodes += geoFile.getInt();
2957           geoFile.skip( nbNodes * INT_WIDTH_GOLD, /*nbLines = */nbFaces );
2958         }
2959         else // standard types
2960         {
2961           if ( isGold )
2962             geoFile.skip( nbNodes, nbElemNodes, INT_WIDTH_GOLD );
2963           else if ( haveElemIds )
2964             geoFile.skip( nbNodes + nbElems, nbElemNodes+1, INT_WIDTH_6 );
2965           else
2966             geoFile.skip( nbNodes, nbElemNodes, INT_WIDTH_6 );
2967         }
2968       }
2969       else if ( word == "coordinates" )
2970       {
2971         isGold = ( nbParts > 0 );
2972         int nbNodes = geoFile.getInt(); // nn
2973
2974         if ( isGold )
2975         {
2976           if ( haveNodeIds )
2977             geoFile.skip( nbNodes, /*nbPerLine =*/ 1, INT_WIDTH_GOLD ); // node ids
2978           geoFile.skip( nbNodes * SPACE_DIM, /*nbPerLine =*/ 1, FLT_WIDTH ); //  coordinates
2979         }
2980         else {
2981           int coordWidth = 3 * FLT_WIDTH;
2982           if ( haveNodeIds )
2983             coordWidth += INT_WIDTH_6;
2984           geoFile.skip(nbNodes, /*nbPerLine =*/ 1, coordWidth);
2985         }
2986       }
2987       else if ( word == "part" )
2988       {
2989         nbParts++;
2990         if ( isGold )
2991           geoFile.skip( 1, /*nbPerLine =*/ 1, INT_WIDTH_GOLD ); //part number
2992         else
2993           geoFile.getWord(); // part number
2994         geoFile.toNextLine();
2995         geoFile.getLine(); // description line
2996       }
2997       else if ( word == "block" )
2998       {
2999         // Structured type
3000         // ------------------
3001         bool rectilinear = ( restLine.find( "rectilinear" ) != restLine.npos );
3002         bool uniform     = ( restLine.find( "uniform" )     != restLine.npos );
3003         bool curvilinear = ( !rectilinear && !uniform );
3004         bool iblanked    = ( restLine.find( "iblanked" )    != restLine.npos );
3005         bool with_ghost  = ( restLine.find( "with_ghost" )  != restLine.npos );
3006         bool range       = ( restLine.find( "range" )       != restLine.npos );
3007
3008         // dimension
3009         int I = geoFile.getInt();
3010         int J = geoFile.getInt();
3011         int K = geoFile.getInt();
3012         int nbNodes = I*J*K;
3013         if ( !nbNodes ) continue;
3014
3015         // range
3016         if ( range ) {
3017           vector<int> ijkRange; // imin imax jmin jmax kmin kmax
3018           ijkRange.reserve(6);
3019           while ( ijkRange.size() < 6 )
3020             ijkRange.push_back( geoFile.getInt() );
3021           I = ijkRange[1]-ijkRange[0]+1;
3022           J = ijkRange[3]-ijkRange[2]+1;
3023           K = ijkRange[5]-ijkRange[4]+1;
3024           nbNodes = I*J*K;
3025         }
3026         int nbElems = (I-1)*(J-1)*(K-1);
3027
3028         if ( curvilinear ) // read coordinates for all nodes
3029         {
3030           if ( isGold )
3031             geoFile.skip( nbNodes, /*nbPerLine =*/ 1, FLT_WIDTH );
3032           else
3033             geoFile.skip( nbNodes * SPACE_DIM, /*nbPerLine =*/ 6, FLT_WIDTH );
3034         }
3035         else if ( rectilinear ) // read delta vectors with non-regular spacing
3036         {
3037           geoFile.skip( I + J + K, /*nbPerLine =*/ 1, FLT_WIDTH );
3038         }
3039         else // uniform: read grid origine and delta vectors for regular spacing grid
3040         {
3041           geoFile.skip( 6, /*nbPerLine =*/ 1, FLT_WIDTH );
3042         }
3043
3044         // iblanks
3045         if ( iblanked ) {
3046           if ( isGold )
3047             geoFile.skip( nbNodes, /*nbPerLine =*/ 1, INT_WIDTH_GOLD);
3048           else
3049             geoFile.skip( nbNodes, /*nbPerLine =*/ 10, INT_WIDTH_6 );
3050         }
3051         // ghosts
3052         if ( with_ghost ) {
3053           geoFile.getWord(); // "ghost_flags"
3054           geoFile.skip( nbElems, /*nbPerLine =*/ 1, INT_WIDTH_GOLD);
3055         }
3056         // node ids
3057         if ( haveNodeIds && geoFile.lookAt( "node_ids" )) {
3058           geoFile.getWord(); // "node_ids"
3059           geoFile.skip( nbNodes, /*nbPerLine =*/ 1, INT_WIDTH_GOLD);
3060         }
3061         // element ids
3062         if ( haveElemIds && geoFile.lookAt( "element_ids" ) ) {
3063           geoFile.getWord(); // "element_ids"
3064           geoFile.skip( nbElems, /*nbPerLine =*/ 1, INT_WIDTH_GOLD);
3065         }
3066       }
3067       else if ( word == "extents" ) {
3068         geoFile.getLine(); geoFile.getLine(); geoFile.getLine();// 3 x 2E12.5
3069       }
3070       else
3071       {
3072         throw MED_EXCEPTION ( LOCALIZED( STRING(LOC) << "Unexpected word: " << word));
3073       }
3074     }
3075   }
3076   return nbParts;
3077 }
3078
3079 //================================================================================
3080 /*!
3081  * \brief count number of parts in EnSight geometry file
3082  */
3083 //================================================================================
3084
3085 int ENSIGHT_MESH_RDONLY_DRIVER::countPartsBinary(_BinaryFileReader& geoFile,
3086                                                  const bool         isSingleFileMode)
3087 {
3088   const char * LOC ="ENSIGHT_MESH_RDONLY_DRIVER::countPartsBinary() : ";
3089
3090   if ( geoFile.getPosition() == 0 ) {
3091     TStrOwner format( geoFile.getLine() ); // "C|Fortran Binary"
3092     if ( !contains( "C Binary", format )) {
3093       if ( contains( "Fortran Binary", format ))
3094         throw(MEDEXCEPTION(STRING(LOC) << "Fortran Binary format not supported"));
3095       else
3096         throw(MEDEXCEPTION(STRING(LOC) << "unexpected line: \n" << format.myValues));
3097     }
3098   }
3099
3100   if ( isSingleFileMode ) {
3101     while (1) {
3102       TStrOwner line( geoFile.getLine() );
3103       if ( isTimeStepBeginning( line.myValues ))
3104         break;
3105     }
3106   }
3107
3108   // 2 description lines
3109   // ----------------------
3110   geoFile.skip( 2 * MAX_LINE_LENGTH );
3111
3112   // read "node|element id <off|given|assign|ignore>"
3113   bool haveNodeIds, haveElemIds;
3114   {
3115     TStrOwner nodeIds( geoFile.getLine() );
3116     TStrOwner elemIds( geoFile.getLine() );
3117
3118     haveNodeIds = ( contains( "given", nodeIds ) || contains( "ignore", nodeIds ) );
3119     haveElemIds = ( contains( "given", elemIds ) || contains( "ignore", elemIds ) );
3120   }
3121
3122   int nbParts = 0; // the result
3123   bool isGold = true;
3124
3125   while ( !geoFile.eof() )
3126   {
3127     TStrOwner line( geoFile.getLine() );
3128     if ( isSingleFileMode && isTimeStepEnd( line.myValues ))
3129       break;
3130     string word, restLine;
3131     TStrTool::split( line.myValues, word, restLine );
3132
3133     const TEnSightElemType & partType = getEnSightType( word );
3134     if ( partType._medType != MED_ALL_ELEMENTS )
3135     {
3136       //  Unstructured element type encounters
3137       // --------------------------------------
3138       int      nbElems = *TIntOwner( geoFile.getInt(1) ); // ne
3139       int  nbElemNodes = partType._medType % 100;
3140
3141       // read element ids
3142       if ( haveElemIds )
3143         geoFile.skip( nbElems * sizeof(int) ); // id_e
3144
3145       int nbNodes = nbElems * nbElemNodes;
3146       if ( partType._name == "nsided" ) // polygons
3147       {
3148         TIntOwner nbNodesInFace( geoFile.getInt( nbElems ));
3149         nbNodes = std::accumulate( nbNodesInFace.myValues, nbNodesInFace.myValues + nbElems, 0 );
3150       }
3151       else if ( partType._name == "nfaced" ) // polyhedrons
3152       {
3153         TIntOwner nbElemFaces( geoFile.getInt( nbElems ));
3154         int nbFaces = accumulate( nbElemFaces.myValues, nbElemFaces.myValues + nbElems, 0 );
3155         TIntOwner nbNodesInFace( geoFile.getInt( nbFaces ));
3156         nbNodes = std::accumulate( nbNodesInFace.myValues, nbNodesInFace.myValues + nbFaces, 0 );
3157       }
3158       geoFile.skip( nbNodes * sizeof(int) );
3159     }
3160     else if ( word == "coordinates" )
3161     {
3162       if ( nbParts == 0 )
3163         isGold = false;
3164       int nbNodes = *TIntOwner( geoFile.getInt(1) ); // nn
3165
3166       // node ids
3167       if ( haveNodeIds )
3168         geoFile.skip( nbNodes * sizeof(int) ); // id_n
3169
3170       // coordinates
3171       geoFile.skip( nbNodes * SPACE_DIM * sizeof(int) );
3172     }
3173     else if ( word == "part" )
3174     {
3175       ++nbParts;
3176       if ( isGold ) geoFile.skip(sizeof(int)); // part #
3177
3178       geoFile.skip(MAX_LINE_LENGTH); // description line
3179     }
3180     else if ( word == "block" )
3181     {
3182       // Structured type
3183       // ------------------
3184       bool rectilinear = ( restLine.find( "rectilinear" ) != restLine.npos );
3185       bool uniform     = ( restLine.find( "uniform" )     != restLine.npos );
3186       bool curvilinear = ( !rectilinear && !uniform );
3187       bool iblanked    = ( restLine.find( "iblanked" )    != restLine.npos );
3188       bool with_ghost  = ( restLine.find( "with_ghost" )  != restLine.npos );
3189       bool range       = ( restLine.find( "range" )       != restLine.npos );
3190
3191       // dimension
3192       TIntOwner ijk( geoFile.getInt(3) );
3193       int I = ijk[0];
3194       int J = ijk[1];
3195       int K = ijk[2];
3196       int NumberOfNodes = I*J*K;
3197       if ( !NumberOfNodes ) {
3198         if ( I != 0 && J != 0 && K != 0 )
3199           throw MEDEXCEPTION( "Need to swap bytes" );
3200         continue;
3201       }
3202
3203       // range
3204       if ( range ) {
3205         TIntOwner ijkRange( geoFile.getInt( 6 ));// imin imax jmin jmax kmin kmax
3206         I = ijkRange[1]-ijkRange[0]+1;
3207         J = ijkRange[3]-ijkRange[2]+1;
3208         K = ijkRange[5]-ijkRange[4]+1;
3209         NumberOfNodes = I*J*K;
3210       }
3211       int nbElems = (I-1)*(J-1)*(K-1);
3212
3213       if ( curvilinear ) // read coordinates for all nodes
3214       {
3215         geoFile.skip( NumberOfNodes * SPACE_DIM * sizeof(float) );
3216       }
3217       else if ( rectilinear ) // read delta vectors with non-regular spacing
3218       {
3219         geoFile.skip( (I+J+K) * sizeof(float) );
3220       }
3221       else // uniform: read grid origine and delta vectors for regular spacing grid
3222       {
3223         geoFile.skip( 6 * sizeof(float) );
3224       }
3225
3226       // iblanks
3227       if ( iblanked )
3228         geoFile.skip( NumberOfNodes * sizeof(int) );
3229       // ghosts
3230       if ( with_ghost ) {
3231         geoFile.skip( MAX_LINE_LENGTH ); // "ghost_flags"
3232         geoFile.skip( nbElems * sizeof(int) );
3233       }
3234       // node ids
3235       if ( haveNodeIds && isGold && !geoFile.eof()  ) {
3236         TStrOwner nextLine( geoFile.getLine() ); // "node_ids"
3237         if ( contains( "node_ids", nextLine ) )
3238           geoFile.skip( NumberOfNodes * sizeof(int) );
3239         else
3240           geoFile.skip( -MAX_LINE_LENGTH );
3241       }
3242       // element ids
3243       if ( haveElemIds && isGold && !geoFile.eof() ) {
3244         TStrOwner nextLine( geoFile.getLine() ); // "element_ids"
3245         if ( contains( "element_ids", nextLine ) )
3246           geoFile.skip( nbElems * sizeof(int) );
3247         else
3248           geoFile.skip( -MAX_LINE_LENGTH );
3249       }
3250     }
3251     else if ( word == "extents" )
3252     {
3253       geoFile.skip( 6 * sizeof(float) );
3254     }
3255     else
3256     {
3257       throw MED_EXCEPTION ( LOCALIZED( STRING(LOC) << "Unexpected word: " << word ));
3258     }
3259   } // while ( !geoFile.eof() )
3260
3261   return nbParts;
3262 }
3263
3264 GROUP* ENSIGHT_MESH_RDONLY_DRIVER::makeGroup( _groupe&     grp,
3265                                               _InterMed & imed)
3266 {
3267   //const char* LOC = "ENSIGHT_MESH_RDONLY_DRIVER::makeGroup(): error";
3268
3269   // prevent creation of other groups but only this one
3270   for (size_t i=0; i < imed.groupes.size(); ++i)
3271     imed.groupes[i].nom.clear();
3272
3273   // let _intermediateMED create a GROUP from grp
3274   grp.medGroup = 0; // the new GROUP should appear in grp.medGroup
3275   grp.nom = "TMP";
3276   vector<GROUP *> tmp;
3277   MESH* mesh = imed._isOwnMedMesh ? (MESH*) 0 : imed._medMesh;
3278   imed.getGroups( tmp, tmp, tmp, tmp, mesh );
3279   if ( !grp.medGroup )
3280     throw MEDEXCEPTION(LOCALIZED("Can't create a GROUP from _groupe"));
3281
3282   grp.medGroup->setName(""); // to let a caller set a proper name
3283   grp.nom = "";
3284
3285   // find pre-existing equal _groupe
3286   _groupe * equalGroupe = 0;
3287   for (unsigned int i=0; i < imed.groupes.size() && !equalGroupe; ++i) {
3288     _groupe& g = imed.groupes[i];
3289     if ( &g != &grp && g.medGroup && g.medGroup->deepCompare( *grp.medGroup ))
3290       equalGroupe = & g;
3291   }
3292   if ( equalGroupe ) {
3293     if ( grp.medGroup->getMesh() )
3294       grp.medGroup->getMesh()->addReference(); // let the mesh survive
3295     grp.medGroup->removeReference();
3296     grp.medGroup = equalGroupe->medGroup;
3297   }
3298   else { // new unique group
3299
3300     if ( grp.medGroup->isOnAllElements() ) // on all elems
3301       grp.medGroup->setName( string("SupportOnAll_")+entNames[ grp.medGroup->getEntity() ] );
3302
3303     // add a new group to mesh
3304     if ( !imed._isOwnMedMesh ) {
3305       vector<GROUP*> * groups = 0;
3306       switch ( grp.medGroup->getEntity() ) {
3307       case MED_CELL: groups = & imed._medMesh->_groupCell; break;
3308       case MED_FACE: groups = & imed._medMesh->_groupFace; break;
3309       case MED_EDGE: groups = & imed._medMesh->_groupEdge; break;
3310       case MED_NODE: groups = & imed._medMesh->_groupNode; break;
3311       default:;
3312       }
3313       if ( groups ) {
3314         groups->resize( groups->size() + 1 );
3315         groups->at( groups->size() - 1) = grp.medGroup;
3316       }
3317     }
3318   }
3319   return grp.medGroup;
3320 }