src/MEDMEM/MEDMEM_EnsightMeshDriver.cxx

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