Back to index

salome-med  6.5.0
MEDMEM_EnsightMeshDriver.cxx
Go to the documentation of this file.
00001 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D, OPEN CASCADE
00002 //
00003 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
00004 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
00005 //
00006 // This library is free software; you can redistribute it and/or
00007 // modify it under the terms of the GNU Lesser General Public
00008 // License as published by the Free Software Foundation; either
00009 // version 2.1 of the License.
00010 //
00011 // This library is distributed in the hope that it will be useful,
00012 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00013 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00014 // Lesser General Public License for more details.
00015 //
00016 // You should have received a copy of the GNU Lesser General Public
00017 // License along with this library; if not, write to the Free Software
00018 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
00019 //
00020 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
00021 //
00022 
00023 #include "MEDMEM_EnsightMeshDriver.hxx"
00024 
00025 #include <sstream>
00026 #include <iomanip>
00027 #include <numeric>
00028 
00029 #include "MEDMEM_define.hxx"
00030 #include "MEDMEM_Family.hxx"
00031 #include "MEDMEM_Group.hxx"
00032 #include "MEDMEM_Coordinate.hxx"
00033 #include "MEDMEM_Connectivity.hxx"
00034 #include "MEDMEM_Mesh.hxx"
00035 #include "MEDMEM_CellModel.hxx"
00036 #include "MEDMEM_Grid.hxx"
00037 
00038 #include "MEDMEM_MedMeshDriver.hxx"
00039 
00040 using namespace std;
00041 using namespace MEDMEM;
00042 using namespace MED_EN;
00043 using namespace MEDMEM_ENSIGHT;
00044 
00045 #define TStrTool _ASCIIFileReader
00046 #define TOLERANCE 1e-15
00047 
00048 //#define ELEMENT_ID_GIVEN
00049 
00050 namespace {
00051 
00052   // ---------------------------------------------------------------------
00057   const char* theDescriptionPrefix = "Meshing from MedMemory. ";
00058 
00059   // ---------------------------------------------------------------------
00063   const char* theDefaultMeshName = "EnsightMesh";
00064 
00065   // ---------------------------------------------------------------------
00069   const int   theMaxNbTypes = 20;
00070 
00071   // ---------------------------------------------------------------------
00076   int* getNumbersByIndex( const int* index, int size, const int* elemNumbers=0)
00077   {
00078     int* numbers = new int[size];
00079     int* n = numbers;
00080     if ( elemNumbers ) {
00081       const int *elem = elemNumbers-1, *elemEnd = elemNumbers + size;
00082       while ( ++elem < elemEnd )
00083         *n++ = index[*elem] - index[*elem-1];
00084     }
00085     else {
00086       const int *ind = index, *indEnd = index + size + 1;
00087       while ( ++ind < indEnd )
00088         *n++ = ind[0] - ind[-1];
00089     }
00090     return numbers;
00091   }
00092   // ---------------------------------------------------------------------
00096   typedef _ValueOwner<int> TNumbers;
00097 
00098 } // namespace
00099 
00100 
00101 ENSIGHT_MESH_DRIVER::ENSIGHT_MESH_DRIVER(): _CaseFileDriver_User(), _ptrMesh((MESH *)NULL)
00102 {
00103 }
00104 
00105 ENSIGHT_MESH_DRIVER::ENSIGHT_MESH_DRIVER(const string & fileName,
00106                                          GMESH *        ptrMesh)
00107   :_CaseFileDriver_User(fileName,MED_EN::RDWR), _ptrMesh(ptrMesh)
00108 {
00109   if ( !_ptrMesh )
00110     throw MEDEXCEPTION("ENSIGHT_MESH_DRIVER(fileName, ptrMesh) : mesh is NULL");
00111   _meshName = _ptrMesh->getName();
00112 }
00113 
00114 ENSIGHT_MESH_DRIVER::ENSIGHT_MESH_DRIVER(const string & fileName,
00115                                          GMESH *        ptrMesh,
00116                                          med_mode_acces accessMode)
00117   :_CaseFileDriver_User(fileName,accessMode), _ptrMesh(ptrMesh)
00118 {
00119   if ( !_ptrMesh )
00120     throw MEDEXCEPTION("ENSIGHT_MESH_DRIVER(fileName, ptrMesh) : mesh is NULL");
00121   _meshName = _ptrMesh->getName();
00122 }
00123 
00124 ENSIGHT_MESH_DRIVER::ENSIGHT_MESH_DRIVER(const ENSIGHT_MESH_DRIVER & driver)
00125   :_CaseFileDriver_User(driver), _ptrMesh(driver._ptrMesh), _meshName(driver._meshName)
00126 {
00127 }
00128 
00129 ENSIGHT_MESH_DRIVER::~ENSIGHT_MESH_DRIVER()
00130 {
00131   MESSAGE_MED("ENSIGHT_MESH_DRIVER::~ENSIGHT_MESH_DRIVER() has been destroyed");
00132 }
00133 
00134 void    ENSIGHT_MESH_DRIVER::setMeshName(const string & meshName) { _meshName = meshName; }
00135 
00136 string  ENSIGHT_MESH_DRIVER::getMeshName() const { return _meshName; }
00137 
00138 void ENSIGHT_MESH_DRIVER::openConst(bool checkDataFile) const
00139 {
00140   const char * LOC ="ENSIGHT_MESH_DRIVER::open() : ";
00141   BEGIN_OF_MED(LOC);
00142 
00143   if ( checkDataFile )
00144   {
00145     if ( getDataFileName().empty() )
00146       throw MED_EXCEPTION
00147         ( LOCALIZED( STRING(LOC) << "Internal error, geometry file name is empty"));
00148 
00149     if (!canOpenFile( getDataFileName(), getAccessMode() ))
00150       throw MED_EXCEPTION
00151         ( LOCALIZED( STRING(LOC) << "Can not open Ensight Geometry file " << getDataFileName()
00152                      << " in access mode " << getAccessMode()));
00153   }
00154   else
00155   {
00156     if ( getCaseFileName().empty() )
00157       throw MED_EXCEPTION
00158         ( LOCALIZED( STRING(LOC) << "Case file name is empty, "
00159                      "please set a correct fileName before calling open()"));
00160 
00161     if ( !canOpenFile( getCaseFileName(), getAccessMode() ))
00162       throw MED_EXCEPTION
00163         ( LOCALIZED( STRING(LOC) << "Can not open Ensight Case file " << getCaseFileName()
00164                      << " in access mode " << getAccessMode()));
00165   }
00166 
00167   END_OF_MED(LOC);
00168 }
00169 
00170 void ENSIGHT_MESH_DRIVER::open() {
00171   openConst() ;
00172 }
00173 
00174 void ENSIGHT_MESH_DRIVER::close() {
00175 }
00176 
00177 // ================================================================================
00178 // WRONLY
00179 // ================================================================================
00180 
00181 ENSIGHT_MESH_WRONLY_DRIVER::ENSIGHT_MESH_WRONLY_DRIVER() :
00182   ENSIGHT_MESH_DRIVER(), _append(0)
00183 {
00184 }
00185 
00186 ENSIGHT_MESH_WRONLY_DRIVER::ENSIGHT_MESH_WRONLY_DRIVER(const string & fileName,
00187                                                        const GMESH *  ptrMesh,
00188                                                        bool           append)
00189   : ENSIGHT_MESH_DRIVER( fileName, (GMESH*)ptrMesh, WRONLY ), _append(append)
00190 {
00191 }
00192 
00193 ENSIGHT_MESH_WRONLY_DRIVER::ENSIGHT_MESH_WRONLY_DRIVER(const ENSIGHT_MESH_WRONLY_DRIVER & driver)
00194   : ENSIGHT_MESH_DRIVER(driver),_append(driver._append)
00195 {
00196 }
00197 
00198 ENSIGHT_MESH_WRONLY_DRIVER::~ENSIGHT_MESH_WRONLY_DRIVER()
00199 {
00200 }
00201 
00202 GENDRIVER * ENSIGHT_MESH_WRONLY_DRIVER::copy() const
00203 {
00204   return new ENSIGHT_MESH_WRONLY_DRIVER(*this) ;
00205 }
00206 
00207 void ENSIGHT_MESH_WRONLY_DRIVER::read() throw (MEDEXCEPTION) {
00208   throw MEDEXCEPTION("ENSIGHT_MESH_WRONLY_DRIVER::read : Can't read with a WRONLY driver !");
00209 }
00210 
00211 //================================================================================
00215 //================================================================================
00216 
00217 void ENSIGHT_MESH_WRONLY_DRIVER::write() const throw (MEDEXCEPTION)
00218 {
00219   const char * LOC = "ENSIGHT_MESH_WRONLY_DRIVER::write() : ";
00220   BEGIN_OF_MED(LOC);
00221 
00222   openConst(false) ; // check if can write to case file
00223 
00224   // Ensight case organization requires a main file (filename.case) which defines organization
00225 
00226   _CaseFileDriver caseFile( getCaseFileName(), this);
00227   if ( _append )
00228     caseFile.read();
00229   caseFile.addMesh( this ); 
00230   caseFile.write(); // all fields of _CaseFileDriver_User are set by this method
00231 
00232   openConst(true) ; // check if can write to data file
00233 
00234   cout << "-> creating the Ensight geometry file " << getDataFileName() << endl ;
00235 
00236   // Store mesh description and a special mark in the first two description lines each
00237   // of 79 chars length maximum, while MED mesh description is up to 200 chars
00238   const char* line1 = theDescriptionPrefix;
00239   string      line2 = _ptrMesh->getDescription();
00240   for ( unsigned i = 0; i < line2.size(); ++i ) { // protect from gabage
00241     if ( !line2[ i ] || !isascii( line2[ i ])) {
00242       line2.resize( i );
00243       break;
00244     }
00245   }
00246   if ((int) line2.size() >= MAX_LINE_LENGTH )
00247     line2.resize( MAX_LINE_LENGTH );
00248 
00249   // EnSight will assign node/element visible numbers it-self
00250   const char* line3 = "node id assign";
00251 #ifdef ELEMENT_ID_GIVEN
00252   const char* line4 = "element id given";
00253 #else
00254   const char* line4 = "element id assign";
00255 #endif
00256 
00257   MESH* mesh = const_cast<MESH*>( _ptrMesh->convertInMESH() ) ; // we write unstructured only
00258 
00259   if ( isBinaryEnSightFormatForWriting() )
00260   {
00261     // ======================================================
00262     //                          Binary
00263     // ======================================================
00264 
00265     _BinaryFileWriter ensightGeomFile( getDataFileName() ); 
00266 
00267     ensightGeomFile.addString("C Binary");
00268     ensightGeomFile.addString(line1);
00269     ensightGeomFile.addString(line2);
00270     ensightGeomFile.addString(line3);
00271     ensightGeomFile.addString(line4);
00272 
00273     // function to write a support as a part
00274     typedef void (ENSIGHT_MESH_WRONLY_DRIVER::* TWritePart) (_BinaryFileWriter&, const SUPPORT*) const;
00275     TWritePart writePart;
00276     if ( isGoldFormat() )
00277     {
00278       // GOLD
00279       writePart = & ENSIGHT_MESH_WRONLY_DRIVER::writePartGoldBinary;
00280     }
00281     else
00282     {
00283       // ENSIGHT 6. Write addionally global nodes
00284       writePart = & ENSIGHT_MESH_WRONLY_DRIVER::writePart6Binary;
00285 
00286       // All point are in 3D, so if we are in 1D or 2D, we complete by zero !
00287       int SpaceDimension = mesh->getSpaceDimension() ;
00288       int NumberOfNodes  = mesh->getNumberOfNodes() ;
00289       ensightGeomFile.addString("coordinates");
00290       ensightGeomFile.addInt( NumberOfNodes );
00291       const double *coordinate = mesh->getCoordinates(MED_FULL_INTERLACE) ;
00292       if ( SpaceDimension == 3 ) {
00293         ensightGeomFile.addReal(coordinate, NumberOfNodes * SpaceDimension );
00294       }
00295       else {
00296         typedef _ValueIterator< double > TComponentIt;
00297         vector< TComponentIt > coordCompIt( 3 );
00298         for (int j=0; j<3; j++, coordinate++)
00299           if ( j < SpaceDimension )
00300             coordCompIt[ j ] = TComponentIt( coordinate, SpaceDimension );
00301           else
00302             coordCompIt[ j ] = TComponentIt(); // to iterate on zeros
00303         ensightGeomFile.addReal( coordCompIt, NumberOfNodes, MED_FULL_INTERLACE );
00304       }
00305     }
00306 
00307     // We put connectivity
00308 
00309     if ( isToWriteEntity( MED_CELL, mesh ))
00310     {
00311       SUPPORT *allCells = const_cast<SUPPORT*>( _ptrMesh->getSupportOnAll( MED_CELL ));
00312       string oldName = allCells->getName();
00313       allCells->setName( _ptrMesh->getName() );
00314       try
00315       {
00316         (this->*writePart)( ensightGeomFile, allCells );
00317       }
00318       catch (MED_EXCEPTION& ex)
00319       {
00320         allCells->setName( oldName );
00321         throw ex;
00322       }
00323       allCells->setName( oldName );
00324     }
00325     // And meshdim-1 connectivity
00326     if ( isToWriteEntity( MED_FACE, mesh ))
00327     {
00328       const SUPPORT *allFaces = _ptrMesh->getSupportOnAll( MED_FACE );
00329       (this->*writePart)( ensightGeomFile, allFaces);
00330     }
00331     else if ( isToWriteEntity(MED_EDGE, mesh))
00332     {
00333       const SUPPORT *allEdges = _ptrMesh->getSupportOnAll( MED_EDGE );
00334       (this->*writePart)( ensightGeomFile, allEdges);
00335     }
00336 
00337     // Write all groups as parts
00338 
00339     for ( int ent = MED_CELL; ent < MED_ALL_ENTITIES; ++ent )
00340     {
00341       medEntityMesh entity = (medEntityMesh) ent;
00342       int nbGroups = mesh->getNumberOfGroups(entity);
00343       for ( int i=1; i<=nbGroups; i++)
00344       {
00345         const GROUP* group = mesh->getGroup( entity, i );
00346         (this->*writePart)( ensightGeomFile, group );
00347       }
00348     }
00349 
00350   }
00351   else
00352   {
00353     // ======================================================
00354     //                           ASCII
00355     // ======================================================
00356     ofstream ensightGeomFile( getDataFileName().c_str(), ios::out); 
00357     ensightGeomFile.setf(ios::scientific);
00358     ensightGeomFile.precision(5);
00359 
00360     ensightGeomFile << line1 << endl 
00361                     << line2 << endl
00362                     << line3 << endl
00363                     << line4 << endl;
00364 
00365     // function to write a support as a part
00366     typedef void (ENSIGHT_MESH_WRONLY_DRIVER::* TWritePart) (ofstream&, const SUPPORT*) const;
00367     TWritePart writePart;
00368     if ( isGoldFormat() )
00369     {
00370       // GOLD
00371       writePart = & ENSIGHT_MESH_WRONLY_DRIVER::writePartGoldASCII;
00372     }
00373     else
00374     {
00375       // ENSIGHT 6. Write addionally global nodes
00376       writePart = & ENSIGHT_MESH_WRONLY_DRIVER::writePart6ASCII;
00377 
00378       // Put points (all point are in 3D, so if we are in 1D or 2D, we complete by zero !
00379       int SpaceDimension = mesh->getSpaceDimension() ;
00380       int NumberOfNodes  = mesh->getNumberOfNodes() ;
00381       string zeros;
00382       if (SpaceDimension==2) zeros = " 0.00000e+00";
00383       if (SpaceDimension==1) zeros = " 0.00000e+00 0.00000e+00";
00384       ensightGeomFile << "coordinates" << endl
00385                       << setw(8) << NumberOfNodes << endl ;
00386       const double *coordinate = mesh->getCoordinates(MED_FULL_INTERLACE) ;
00387       for (int i=0; i<NumberOfNodes; i++)
00388       {
00389         //ensightGeomFile << setw(8) << i+1 ; // node id
00390         for (int j=0; j<SpaceDimension; j++, coordinate++)
00391           ensightGeomFile << setw(12) << *coordinate;
00392         ensightGeomFile << zeros << endl ;
00393       }
00394     }
00395 
00396     // We put connectivity
00397 
00398     if ( isToWriteEntity( MED_CELL, mesh ))
00399     {
00400       SUPPORT *allCells = const_cast<SUPPORT*>( _ptrMesh->getSupportOnAll( MED_CELL ));
00401       string oldName = allCells->getName();
00402       allCells->setName( _ptrMesh->getName() );
00403       try
00404       {
00405         (this->*writePart)( ensightGeomFile, allCells );
00406       }
00407       catch (MED_EXCEPTION& ex)
00408       {
00409         allCells->setName( oldName );
00410         throw ex;
00411       }
00412       allCells->setName( oldName );
00413     }
00414     // And meshdim-1 connectivity
00415     if ( isToWriteEntity( MED_FACE, mesh ))
00416     {
00417       const SUPPORT *allFaces = _ptrMesh->getSupportOnAll( MED_FACE );
00418       (this->*writePart)( ensightGeomFile, allFaces);
00419     }
00420     else if ( isToWriteEntity(MED_EDGE, mesh))
00421     {
00422       const SUPPORT *allEdges = _ptrMesh->getSupportOnAll( MED_EDGE );
00423       (this->*writePart)( ensightGeomFile, allEdges);
00424     }
00425 
00426     // Write all groups as parts
00427 
00428     for ( int ent = MED_CELL; ent < MED_ALL_ENTITIES; ++ent )
00429     {
00430       medEntityMesh entity = (medEntityMesh) ent;
00431       int nbGroups = mesh->getNumberOfGroups(entity);
00432       for ( int i=1; i<=nbGroups; i++)
00433       {
00434         const GROUP* group = mesh->getGroup( entity, i );
00435         (this->*writePart)( ensightGeomFile, group );
00436       }
00437     }
00438 
00439     ensightGeomFile.close();
00440 
00441     mesh->removeReference();
00442 
00443   } // end ASCII format
00444 
00445 } // ENSIGHT_MESH_WRONLY_DRIVER::write()
00446 
00447 //================================================================================
00451 //================================================================================
00452 
00453 void ENSIGHT_MESH_WRONLY_DRIVER::writePartGoldBinary(_BinaryFileWriter& ensightGeomFile,
00454                                                      const SUPPORT*     support) const
00455 {
00456   // part number
00457   int partNum = getPartNumber( support );
00458   if ( !partNum )
00459     throw MED_EXCEPTION ( LOCALIZED( STRING("Internal error: invalid part number")));
00460   ensightGeomFile.addString( "part" );
00461   ensightGeomFile.addInt( partNum );
00462 
00463   // group/mesh name
00464   ensightGeomFile.addString( support->getName() );
00465 
00466   // get geom types
00467   medEntityMesh entity = support->getEntity();
00468   int nbTypes = support->getNumberOfTypes();
00469   const medGeometryElement* geoType = support->getTypes();
00470 
00471   const int * connectivity = 0;
00472   const int * elemConnectivity = 0;
00473   const int * index = 0;
00474   int j;
00475 
00476   const MESH* mesh = support->getMesh()->convertInMESH() ; // we write unstructured only
00477 
00478   // COORDINATES                                                             Gold binary
00479   // ===================================================================================
00480   // In Ensight, coordinates of nodes of support elements are in MED_NO_INTERLACE mode.
00481   // We are to write only nodes belonging to elements of the support and
00482   // nodal connectivity should refer to these nodes.
00483   map<int, int> med2ensIds;
00484   map<int, int>::iterator medEnsIt;
00485   int SpaceDimension = _ptrMesh->getSpaceDimension() ;
00486   int NumberOfNodes  = _ptrMesh->getNumberOfNodes() ;
00487   // -------------------------------------------------
00488   if ( support->isOnAllElements() )
00489   {
00490     // nb of nodes
00491     ensightGeomFile.addString( "coordinates" );
00492     ensightGeomFile.addInt( NumberOfNodes );
00493 
00494     // coordinates
00495     const double *coordinate = mesh->getCoordinates(MED_FULL_INTERLACE);
00496     typedef _ValueIterator< double > TComponentIt;
00497     vector< TComponentIt > coordCompIt( 1 );
00498     for (int j=0; j<SPACE_DIM; j++, coordinate++) { // loop on dimensions
00499       if ( j < SpaceDimension )
00500         coordCompIt[ 0 ] = TComponentIt( coordinate, SpaceDimension );
00501       else
00502         coordCompIt[ 0 ] = TComponentIt(); // to iterate on zeros
00503       ensightGeomFile.addReal( coordCompIt, NumberOfNodes, MED_NO_INTERLACE );
00504     }
00505   }
00506   // -------------------------------------------------
00507   else // support is not on all elements
00508   {
00509     // nb of nodes
00510     getSupportNodes( support, med2ensIds );
00511     NumberOfNodes = med2ensIds.size();
00512     ensightGeomFile.addString( "coordinates" );
00513     ensightGeomFile.addInt( NumberOfNodes );
00514 
00515     // coordinates
00516     vector<float> floatCoords( NumberOfNodes );
00517     for ( j=0; j < SPACE_DIM; j++) { // loop on dimensions
00518       medEnsIt = med2ensIds.begin();
00519       if ( j < SpaceDimension ) {
00520         const double *coordinate = mesh->getCoordinates(MED_FULL_INTERLACE) + j;
00521         for (int i=0; i<NumberOfNodes; i++, ++medEnsIt )
00522           floatCoords[ i ] = (float) coordinate[ (medEnsIt->first-1) * SpaceDimension];
00523       }
00524       else if ( j-1 < SpaceDimension ) {
00525         for (int i=0; i<NumberOfNodes; i++)
00526           floatCoords[ i ] = 0.;
00527       }
00528       ensightGeomFile.addReal( &floatCoords[0], NumberOfNodes );
00529     }
00530     // assign local node ids
00531     for ( medEnsIt = med2ensIds.begin(), j=1; j<=NumberOfNodes; j++, ++medEnsIt )
00532       medEnsIt->second = j;
00533   }
00534 
00535   // CONNECTIVITY                                                            Gold binary
00536   // ===================================================================================
00537   // loop on types
00538   for (int i=0; i<nbTypes; i++)
00539   {
00540     const medGeometryElement    medType = geoType[i];
00541     const TEnSightElemType& ensightType = getEnSightType(medType);
00542     const int numberOfCell              = support->getNumberOfElements(medType);
00543     int nbCellNodes                     = ensightType._medIndex.size();
00544 
00545     // type name and nb cells
00546     ensightGeomFile.addString( ensightType._name );
00547     ensightGeomFile.addInt(  numberOfCell );
00548 
00549     vector<int> nodeIds;
00550 
00551     // -------------------------------------------------
00552     if ( support->isOnAllElements() )
00553     {
00554 #ifdef ELEMENT_ID_GIVEN
00555       // elem numbers
00556       nodeIds.resize( numberOfCell );
00557       for ( j = 1; j <= numberOfCell; j++)
00558         nodeIds[ j-1 ] = j;
00559       ensightGeomFile.addInt( nodeIds );
00560 
00561       if ( entity != MED_NODE ) nodeIds.clear();
00562 #endif
00563 
00564       if ( nbCellNodes > 1 ) // STANDARD ELEMENTS connectivity
00565       {
00566         connectivity = mesh->getConnectivity( MED_NODAL, entity, medType);
00567         nodeIds.reserve( numberOfCell * nbCellNodes);
00568         for (j = 0 ; j < numberOfCell; j++, connectivity += nbCellNodes)
00569           for (int k=0; k<nbCellNodes; k++)
00570             nodeIds.push_back( connectivity[ ensightType._medIndex[k] ]);
00571         ensightGeomFile.addInt( nodeIds );
00572       }
00573       else if ( entity == MED_NODE ) // NODES connectivity
00574       {
00575 #if !defined(ELEMENT_ID_GIVEN)
00576         nodeIds.resize( numberOfCell );
00577         for ( j = 1; j <= numberOfCell; j++)
00578           nodeIds[ j-1 ] = j;
00579 #endif
00580         ensightGeomFile.addInt( nodeIds );
00581       }
00582       else if ( medType == MED_POLYGON ) // POLYGONs connectivity
00583       {
00584         int nbStdCells = mesh->getGlobalNumberingIndex(entity)[i]-1;
00585         connectivity   = mesh->getConnectivity( MED_NODAL, entity, medType);
00586         int connLength = mesh->getConnectivityLength( MED_NODAL, entity, medType);
00587         index          = mesh->getConnectivityIndex(MED_NODAL, entity);
00588 
00589         // number of nodes in each element
00590         {
00591           TIntOwner nbNodesInPoly( getNumbersByIndex( index+nbStdCells, numberOfCell ));
00592           ensightGeomFile.addInt( nbNodesInPoly.myValues, numberOfCell );
00593         } // nbNodesInPoly is deleted here
00594 
00595         // connectivity
00596         ensightGeomFile.addInt( connectivity, connLength );
00597       }
00598       else // POLYHEDRA connectivity
00599       {
00600         connectivity   = mesh->getConnectivity( MED_NODAL, entity, medType);
00601         int nbStdCells = mesh->getGlobalNumberingIndex(entity)[i]-1;
00602         index          = mesh->getConnectivityIndex(MED_NODAL, entity) + nbStdCells;
00603 
00604         vector<int> nbFacesInPolyhedron, nbNodesInFace, faceConn;
00605         for ( int j = 0; j < numberOfCell; ++j )
00606         {
00607           int nbFaces = 0, nbNodes = 0;
00608           for ( int k = index[j]; k < index[j+1]; ++k )
00609             if ( connectivity[k-1] == -1 )
00610             {
00611               nbNodesInFace.push_back( nbNodes );
00612               nbNodes = 0;
00613               ++nbFaces;
00614             }
00615             else
00616             {
00617               faceConn.push_back( connectivity[k-1] );
00618               ++nbNodes;
00619             }
00620           nbNodesInFace.push_back( nbNodes );
00621           nbFacesInPolyhedron.push_back( nbFaces+1 );
00622         }
00623 
00624         // nb of faces in each polyhedron
00625         ensightGeomFile.addInt( nbFacesInPolyhedron );
00626         // number of nodes in each face
00627         ensightGeomFile.addInt( nbNodesInFace );
00628         // connectivity
00629         ensightGeomFile.addInt( faceConn );
00630       }
00631     }
00632     // -------------------------------------------------
00633     else // support is not on all elements
00634     {
00635       const int *number = support->getNumber(medType);
00636 
00637 #ifdef ELEMENT_ID_GIVEN
00638       ensightGeomFile.addInt( number, numberOfCell );
00639 #endif
00640       if ( nbCellNodes > 1 ) // STANDARD ELEMENTS connectivity
00641       {
00642         connectivity = mesh->getConnectivity( MED_NODAL, entity, MED_ALL_ELEMENTS);
00643         index = mesh->getConnectivityIndex(MED_NODAL, entity);
00644 
00645         nodeIds.reserve( numberOfCell * nbCellNodes);
00646         for (j=0; j<numberOfCell; j++) {
00647           int elem = number[j];
00648           elemConnectivity = connectivity + index[elem-1]-1;
00649           for (int k=0; k<nbCellNodes; k++)
00650           {
00651             int node = elemConnectivity[ ensightType._medIndex[k] ];
00652             nodeIds.push_back( med2ensIds[ node ]);
00653           }
00654         }
00655         ensightGeomFile.addInt( nodeIds );
00656       }
00657       else if ( entity == MED_NODE )  // NODES connectivity
00658       {
00659         nodeIds.resize( numberOfCell );
00660         for ( j = 1; j <= numberOfCell; j++)
00661           nodeIds[ j-1 ] = j;
00662         ensightGeomFile.addInt( nodeIds );
00663       }
00664       else if ( medType == MED_POLYGON ) // POLYGONs connectivity
00665       {
00666         connectivity = mesh->getConnectivity( MED_NODAL, entity, MED_ALL_ELEMENTS);
00667         index = mesh->getConnectivityIndex(MED_NODAL, entity);
00668 
00669         // number of nodes in each element
00670         {
00671           TIntOwner nbNodesInPoly( getNumbersByIndex( index, numberOfCell, number ));
00672           ensightGeomFile.addInt( nbNodesInPoly.myValues, numberOfCell );
00673         } // nbNodesInPoly is deleted here
00674 
00675         // connectivity
00676         for ( j = 0; j < numberOfCell; ++j )
00677         {
00678           int elem = number[ j ];
00679           elemConnectivity   = connectivity + index[ elem-1 ]-1;
00680           const int* connEnd = connectivity + index[ elem   ]-1;
00681           while ( elemConnectivity < connEnd )
00682             nodeIds.push_back( med2ensIds[ *elemConnectivity++ ]);
00683         }
00684         ensightGeomFile.addInt( nodeIds );
00685       }
00686       else // POLYHEDRA connectivity
00687       {
00688         connectivity = mesh->getConnectivity( MED_NODAL, entity, MED_ALL_ELEMENTS);
00689         index = mesh->getConnectivityIndex(MED_NODAL, entity);
00690         vector<int> nbFacesInPolyhedron, nbNodesInFace, faceConn;
00691         for ( int j = 0; j < numberOfCell; ++j )
00692         {
00693           int elem    = number[ j ];
00694           int nbFaces = 0, nbNodes = 0;
00695           for ( int k = index[elem]; k < index[elem+1]; ++k )
00696             if ( connectivity[k-1] == -1 )
00697             {
00698               nbNodesInFace.push_back( nbNodes );
00699               nbNodes = 0;
00700               ++nbFaces;
00701             }
00702             else
00703             {
00704               faceConn.push_back( connectivity[k-1] );
00705               ++nbNodes;
00706             }
00707           nbNodesInFace.push_back( nbNodes );
00708           nbFacesInPolyhedron.push_back( nbFaces+1 );
00709         }
00710 
00711         // nb of faces in each polyhedron
00712         ensightGeomFile.addInt( nbFacesInPolyhedron );
00713         // number of nodes in each face
00714         ensightGeomFile.addInt( nbNodesInFace );
00715         // connectivity
00716         ensightGeomFile.addInt( faceConn );
00717       }
00718     }
00719   }
00720 
00721   mesh->removeReference();
00722 
00723 } // writePartGoldBinary()
00724 
00725 //================================================================================
00729 //================================================================================
00730 
00731 void ENSIGHT_MESH_WRONLY_DRIVER::writePartGoldASCII(ofstream&      ensightGeomFile,
00732                                                     const SUPPORT* support) const
00733 {
00734   const int iw = 10;
00735 
00736   // part number
00737   int partNum = getPartNumber( support );
00738   ensightGeomFile << "part" << endl
00739                   << setw(iw) << partNum << endl;
00740   if ( !partNum )
00741     throw MED_EXCEPTION ( LOCALIZED( STRING("Internal error: invalid part number")));
00742 
00743   // group/mesh name
00744   ensightGeomFile << support->getName() << endl;
00745 
00746   // get geom types
00747   medEntityMesh entity = support->getEntity();
00748   int nbTypes = support->getNumberOfTypes();
00749   const medGeometryElement* geoType = support->getTypes();
00750 
00751   const int * connectivity = 0;
00752   const int * elemConnectivity = 0;
00753   const int * index = 0;
00754   int j;
00755 
00756   const MESH* mesh = support->getMesh()->convertInMESH() ; // we write unstructured only
00757 
00758   // COORDINATES                                                              Gold ASCII 
00759   // ===================================================================================
00760   // In Ensight, coordinates of nodes of support elements are in MED_NO_INTERLACE mode.
00761   // We are to write only nodes belonging to elements of the support and
00762   // nodal connectivity should refer to these nodes.
00763   map<int, int> med2ensIds;
00764   map<int, int>::iterator medEnsIt;
00765   int SpaceDimension = mesh->getSpaceDimension() ;
00766   int NumberOfNodes  = mesh->getNumberOfNodes() ;
00767   string zeroStr = " 0.00000e+00";
00768   // -----------------------------------
00769   if ( support->isOnAllElements() )
00770   {
00771     // nb of nodes
00772     ensightGeomFile << "coordinates" << endl
00773                     << setw(iw) << NumberOfNodes << endl ;
00774 
00775     // coordinates
00776     for (j=0; j<SPACE_DIM; j++) { // loop on dimensions
00777       if ( j < SpaceDimension ) {
00778         const double *coordinate = mesh->getCoordinates(MED_FULL_INTERLACE) + j;
00779         for (int i=0; i<NumberOfNodes; i++, coordinate += SpaceDimension)
00780           ensightGeomFile << setw(12) << (float) *coordinate << endl;
00781       }
00782       else {
00783         for (int i=0; i<NumberOfNodes; i++)
00784           ensightGeomFile << zeroStr << endl;
00785       }
00786     }
00787   }
00788   // -----------------------------------
00789   else // support is not on all elements
00790   {
00791     // nb of nodes
00792     getSupportNodes( support, med2ensIds );
00793     NumberOfNodes = med2ensIds.size();
00794     ensightGeomFile << "coordinates" << endl
00795                     << setw(iw) << NumberOfNodes << endl ;
00796 
00797     // coordinates
00798     for ( j=0; j<SPACE_DIM; j++) { // loop on dimensions
00799       medEnsIt = med2ensIds.begin();
00800       if ( j < SpaceDimension ) {
00801         const double *coordinate = mesh->getCoordinates(MED_FULL_INTERLACE) + j;
00802         for (int i=0; i<NumberOfNodes; i++, ++medEnsIt )
00803           ensightGeomFile << setw(12)
00804                           << (float) coordinate[ (medEnsIt->first-1) * SpaceDimension] << endl;
00805       }
00806       else {
00807         for (int i=0; i<NumberOfNodes; i++)
00808           ensightGeomFile << zeroStr << endl;
00809       }
00810     }
00811     // assign local node ids
00812     for ( medEnsIt = med2ensIds.begin(), j=1; j<=NumberOfNodes; j++, ++medEnsIt )
00813       medEnsIt->second = j;
00814   }
00815 
00816   // CONNECTIVITY                                                             Gold ASCII
00817   // ===================================================================================
00818   // loop on types
00819   for (int i=0; i<nbTypes; i++)
00820   {
00821     const medGeometryElement    medType = geoType[i];
00822     const TEnSightElemType& ensightType = getEnSightType(medType);
00823     const int numberOfCell              = support->getNumberOfElements(medType);
00824     int nbCellNodes                     = ensightType._medIndex.size();
00825 
00826     // type name and nb cells
00827     ensightGeomFile << ensightType._name        << endl
00828                     << setw(iw) << numberOfCell << endl;
00829 
00830     // -----------------------------------
00831     if ( support->isOnAllElements() )
00832     {
00833 #ifdef ELEMENT_ID_GIVEN
00834       for ( j = 1; j <= numberOfCell; j++)
00835         ensightGeomFile << setw(iw) << j << endl;
00836 #endif
00837 
00838       if ( nbCellNodes > 1 ) // STANDARD ELEMENTS connectivity
00839       {
00840         connectivity = mesh->getConnectivity( MED_NODAL, entity, medType);
00841         for (j = 0 ; j < numberOfCell; j++, connectivity += nbCellNodes) {
00842           for (int k=0; k<nbCellNodes; k++)
00843             ensightGeomFile << setw(iw) << connectivity[ ensightType._medIndex[k] ];
00844           ensightGeomFile << endl ;
00845         }
00846       }
00847       else if ( entity == MED_NODE ) // NODES connectivity
00848       {
00849         for ( j = 1; j <= numberOfCell; j++)
00850           ensightGeomFile << setw(iw) << j << endl;
00851       }
00852       else if ( medType == MED_POLYGON ) // POLYGONs connectivity
00853       {
00854         int nbStdCells = mesh->getGlobalNumberingIndex(entity)[i]-1;
00855         connectivity   = mesh->getConnectivity( MED_NODAL, entity, MED_ALL_ELEMENTS);
00856         index          = mesh->getConnectivityIndex(MED_NODAL, entity) + nbStdCells;
00857         // number of nodes in each element
00858         const int* ind = index;
00859         for (j = 0 ; j < numberOfCell; j++, ++ind)
00860           ensightGeomFile << setw(iw) << ( ind[1] - ind[0] ) << endl;
00861 
00862         // connectivity
00863         for (j = 0; j < numberOfCell; j++, ++index) {
00864           nbCellNodes = index[1] - index[0];
00865           for (int k=0; k<nbCellNodes; k++, ++connectivity)
00866             ensightGeomFile << setw(iw) << *connectivity;
00867           ensightGeomFile << endl;
00868         }
00869       }
00870       else // POLYHEDRA connectivity
00871       {
00872         int nbStdCells = mesh->getGlobalNumberingIndex(entity)[i]-1;
00873         connectivity   = mesh->getConnectivity( MED_NODAL, entity, medType);
00874         index          = mesh->getConnectivityIndex(MED_NODAL, entity) + nbStdCells;
00875         ostringstream nbNodesInFace, faceConn;
00876         for ( j = 0; j < numberOfCell; ++j )
00877         {
00878           int nbFaces = 0, nbNodes = 0;
00879           for ( int k = index[j]; k < index[j+1]; ++k )
00880             if ( connectivity[k-1] == -1 )
00881             {
00882               faceConn << endl;
00883               nbNodesInFace << setw(iw) << nbNodes << endl;
00884               nbNodes = 0;
00885               ++nbFaces;
00886             }
00887             else
00888             {
00889               faceConn << setw(iw) << connectivity[k-1] ;
00890               ++nbNodes;
00891             }
00892           faceConn << endl;
00893           nbNodesInFace << setw(iw) << nbNodes << endl;
00894           ensightGeomFile << setw(iw) << nbFaces+1 << endl ;// nb of faces in each polyhedron
00895         }
00896         ensightGeomFile << nbNodesInFace.str();// number of nodes in each face
00897         ensightGeomFile << faceConn.str();// connectivity of each face
00898       }
00899     }
00900     // -----------------------------------
00901     else // support is not on all elements
00902     {
00903       const int *number = support->getNumber(medType);
00904 
00905 #ifdef ELEMENT_ID_GIVEN
00906       for ( j = 0; j < numberOfCell; j++)
00907         ensightGeomFile << setw(iw) << number[j] << endl;
00908 #endif
00909       if ( nbCellNodes > 1 ) // STANDARD ELEMENTS connectivity
00910       {
00911         connectivity = mesh->getConnectivity( MED_NODAL, entity, MED_ALL_ELEMENTS);
00912         index = mesh->getConnectivityIndex(MED_NODAL, entity);
00913 
00914         for (j=0; j<numberOfCell; j++) {
00915           int elem = number[j];
00916           elemConnectivity = connectivity + index[elem-1]-1;
00917           for (int k=0; k<nbCellNodes; k++) {
00918             int node = elemConnectivity[ ensightType._medIndex[k] ];
00919             ensightGeomFile << setw(iw) << med2ensIds[ node ];
00920           }
00921           ensightGeomFile << endl;
00922         }
00923       }
00924       else if ( entity == MED_NODE )  // NODES connectivity
00925       {
00926         for (j=0; j<numberOfCell; j++) {
00927           int node = med2ensIds[ number[j] ];
00928           ensightGeomFile << setw(iw) << node << endl ;
00929         }
00930       }
00931       else if ( medType == MED_POLYGON ) // POLYGONs connectivity
00932       {
00933         connectivity   = mesh->getConnectivity( MED_NODAL, entity, MED_ALL_ELEMENTS);
00934         index          = mesh->getConnectivityIndex(MED_NODAL, entity);
00935         // number of nodes in each element
00936         for (j = 0 ; j < numberOfCell; j++) {
00937           int elem = number[j];
00938           ensightGeomFile << setw(iw) << ( index[elem] - index[elem-1] ) << endl;
00939         }
00940         // connectivity
00941         for ( j = 0; j < numberOfCell; ++j ) {
00942           int elem = number[ j ];
00943           elemConnectivity   = connectivity + index[ elem-1 ]-1;
00944           const int* connEnd = connectivity + index[ elem   ]-1;
00945           while ( elemConnectivity < connEnd )
00946             ensightGeomFile << setw(iw) << med2ensIds[ *elemConnectivity++ ];
00947           ensightGeomFile << endl;
00948         }
00949       }
00950       else // POLYHEDRA connectivity
00951       {
00952         connectivity   = mesh->getConnectivity( MED_NODAL, entity, medType);
00953         index          = mesh->getConnectivityIndex(MED_NODAL, entity);
00954         ostringstream nbNodesInFace, faceConn;
00955         for ( j = 0; j < numberOfCell; ++j )
00956         {
00957           int elem = number[j];
00958           int nbFaces = 0, nbNodes = 0;
00959           for ( int k = index[elem]; k < index[elem+1]; ++k )
00960             if ( connectivity[k-1] == -1 )
00961             {
00962               faceConn << endl;
00963               nbNodesInFace << setw(iw) << nbNodes << endl;
00964               nbNodes = 0;
00965               ++nbFaces;
00966             }
00967             else
00968             {
00969               faceConn << setw(iw) << connectivity[k-1] ;
00970               ++nbNodes;
00971             }
00972           faceConn << endl;
00973           nbNodesInFace << setw(iw) << nbNodes << endl;
00974           ensightGeomFile << setw(iw) << nbFaces+1 << endl ;// nb of faces in each polyhedron
00975         }
00976         ensightGeomFile << nbNodesInFace.str();// number of nodes in each face
00977         ensightGeomFile << faceConn.str();// connectivity of each face
00978       }
00979     }
00980   }
00981 
00982   mesh->removeReference();
00983 
00984 }  // writePartGoldASCII()
00985 
00986 //================================================================================
00990 //================================================================================
00991 
00992 void ENSIGHT_MESH_WRONLY_DRIVER::writePart6Binary(_BinaryFileWriter& ensightGeomFile,
00993                                                   const SUPPORT*     support) const
00994 {
00995   // part number
00996   int partNum = getPartNumber( support );
00997   ensightGeomFile.addString( STRING("part ") << partNum );
00998   if ( !partNum )
00999     throw MED_EXCEPTION ( LOCALIZED( STRING("Internal error: invalid part number")));
01000 
01001   // group/mesh name
01002   ensightGeomFile.addString( support->getName() );
01003 
01004   // get geom types
01005   medEntityMesh entity = support->getEntity();
01006   int nbTypes = support->getNumberOfTypes();
01007   const medGeometryElement* geoType = support->getTypes();
01008 
01009   const MESH* mesh = support->getMesh()->convertInMESH() ; // we write unstructured only
01010 
01011   int j = 1;
01012   const int * connectivity = 0;
01013   if ( entity != MED_NODE )
01014     connectivity = mesh->getConnectivity( MED_NODAL, entity, MED_ALL_ELEMENTS);
01015   const int * elemConnectivity = connectivity;
01016 
01017   // CONNECTIVITY                                                       Ensight 6 binary
01018   // ===================================================================================
01019   // loop on types
01020   for (int i=0; i<nbTypes; i++)
01021   {
01022     const medGeometryElement    medType = geoType[i];
01023     const TEnSightElemType& ensightType = getEnSightType(medType);
01024     int nbCellNodes = ensightType._medIndex.size();
01025     if ( nbCellNodes == 0 )
01026       continue; // poly?
01027 
01028     // type name and nb cells
01029     int numberOfCell = support->getNumberOfElements(medType);
01030     ensightGeomFile.addString( ensightType._name );
01031     ensightGeomFile.addInt( numberOfCell );
01032 
01033     vector<int> nodeIds;
01034     // -------------------------------------------------
01035     if ( support->isOnAllElements() )
01036     {
01037 #ifdef ELEMENT_ID_GIVEN
01038       nodeIds.resize( numberOfCell );
01039       for ( j = 1; j <= numberOfCell; j++)
01040         nodeIds[ j-1 ] = j;
01041       ensightGeomFile.addInt( nodeIds );
01042 #endif
01043       if ( entity == MED_NODE ) {
01044 #if !defined(ELEMENT_ID_GIVEN)
01045         nodeIds.resize( numberOfCell * nbCellNodes);
01046         for ( j = 1; j <= numberOfCell; j++)
01047           nodeIds[ j-1 ] = j;
01048 #endif
01049       }
01050       else {
01051         nodeIds.clear();
01052         nodeIds.reserve( numberOfCell * nbCellNodes );
01053         for (j = 0 ; j < numberOfCell; j++, elemConnectivity += nbCellNodes)
01054           for (int k=0; k<nbCellNodes; k++)
01055             nodeIds.push_back( elemConnectivity[ ensightType._medIndex[k] ]);
01056       }
01057       ensightGeomFile.addInt( nodeIds );
01058     }
01059     // -------------------------------------------------
01060     else // support is not on all elements
01061     {
01062       const int *number = support->getNumber(medType);
01063 
01064 #ifdef ELEMENT_ID_GIVEN
01065       ensightGeomFile.addInt( number, numberOfCell );
01066 #endif
01067       if ( entity == MED_NODE ) {
01068         ensightGeomFile.addInt( number, numberOfCell );
01069       }
01070       else {
01071         const int* index = mesh->getConnectivityIndex(MED_NODAL, entity);
01072 
01073         nodeIds.reserve( numberOfCell * nbCellNodes);
01074         for (j=0; j<numberOfCell; j++) {
01075           int elem = number[j];
01076           elemConnectivity = connectivity + index[elem-1]-1;
01077           for (int k=0; k<nbCellNodes; k++)
01078             nodeIds.push_back( elemConnectivity[ ensightType._medIndex[k] ]);
01079         }
01080         ensightGeomFile.addInt( nodeIds );
01081       }
01082     }
01083   } // loop on types
01084 
01085   mesh->removeReference();
01086 
01087 
01088 } // writePart6Binary()
01089 
01090 //================================================================================
01094 //================================================================================
01095 
01096 void ENSIGHT_MESH_WRONLY_DRIVER::writePart6ASCII(ofstream&      ensightGeomFile,
01097                                                  const SUPPORT* support) const
01098 {
01099   const int iw = 8;
01100 
01101   // part number
01102   int partNum = getPartNumber( support );
01103   ensightGeomFile << "part " << partNum << endl;
01104   if ( !partNum )
01105     throw MED_EXCEPTION ( LOCALIZED( STRING("Internal error: invalid part number")));
01106 
01107   // group/mesh name
01108   ensightGeomFile << support->getName() << endl;
01109 
01110   // get geom types
01111   medEntityMesh entity = support->getEntity();
01112   int nbTypes = support->getNumberOfTypes();
01113   const medGeometryElement* geoType = support->getTypes();
01114 
01115   const MESH* mesh = support->getMesh()->convertInMESH() ; // we write unstructured only
01116 
01117   int j = 1;
01118   const int * connectivity = 0;
01119   if ( entity != MED_NODE )
01120     connectivity = mesh->getConnectivity( MED_NODAL,entity, MED_ALL_ELEMENTS);
01121   const int * elemConnectivity = connectivity;
01122 
01123   // CONNECTIVITY                                                        Ensight 6 ASCII
01124   // ===================================================================================
01125   // loop on types
01126   for (int i=0; i<nbTypes; i++)
01127   {
01128     const medGeometryElement    medType = geoType[i];
01129     const TEnSightElemType& ensightType = getEnSightType(medType);
01130     int nbCellNodes = ensightType._medIndex.size();
01131     if ( nbCellNodes == 0 )
01132       continue; // poly?
01133 
01134     // type name and nb cells
01135     int numberOfCell = support->getNumberOfElements(medType);
01136     ensightGeomFile << ensightType._name       << endl
01137                     << setw(iw) << numberOfCell << endl;
01138 
01139     // -------------------------------------------------
01140     if ( support->isOnAllElements() )
01141     {
01142       if ( entity == MED_NODE ) {
01143         for ( j = 1; j <= numberOfCell; j++) {
01144 #ifdef ELEMENT_ID_GIVEN
01145           ensightGeomFile << setw(iw) << j;
01146 #endif
01147           ensightGeomFile << setw(iw) << j << endl;
01148         }
01149       }
01150       else {
01151         for (j = 1 ; j <= numberOfCell; j++, elemConnectivity += nbCellNodes) {
01152 #ifdef ELEMENT_ID_GIVEN
01153           ensightGeomFile << setw(iw) << elem++;
01154 #endif
01155           for (int k=0; k<nbCellNodes; k++)
01156           {
01157             ensightGeomFile << setw(iw) << elemConnectivity[ ensightType._medIndex[k] ];
01158           }
01159           ensightGeomFile << endl ;
01160         }
01161       }
01162     }
01163     // -------------------------------------------------
01164     else  // support is not on all elements
01165     {
01166       const int *number = support->getNumber(medType);
01167       if ( entity == MED_NODE ) {
01168         for (j=0; j<numberOfCell; j++) {
01169           int node = number[j];
01170 #ifdef ELEMENT_ID_GIVEN
01171           ensightGeomFile << setw(iw) << node;
01172 #endif
01173           ensightGeomFile << setw(iw) << node << endl ;
01174         }
01175       }
01176       else {
01177         const int* index = mesh->getConnectivityIndex(MED_NODAL, entity);
01178 
01179         for (j=0; j<numberOfCell; j++) {
01180           int elem = number[j];
01181 #ifdef ELEMENT_ID_GIVEN
01182           ensightGeomFile << setw(iw) << elem;
01183 #endif
01184           elemConnectivity = connectivity + index[elem-1]-1;
01185           for (int k=0; k<nbCellNodes; k++)
01186           {
01187             ensightGeomFile << setw(iw) << elemConnectivity[ ensightType._medIndex[k] ];
01188           }
01189           ensightGeomFile << endl ;
01190         }
01191       }
01192     }
01193   } // loop on types
01194 
01195   mesh->removeReference();
01196 
01197 
01198 } // writePart6ASCII()
01199 
01200 //================================================================================
01204 //================================================================================
01205 
01206 int ENSIGHT_MESH_WRONLY_DRIVER::nbPartsToWrite() const
01207 {
01208   int nbParts = 0;
01209   nbParts += (int) isToWriteEntity( MED_CELL, _ptrMesh );
01210   nbParts += (int) isToWriteEntity( MED_FACE, _ptrMesh );
01211   nbParts += (int) isToWriteEntity( MED_EDGE, _ptrMesh );
01212 
01213   // all groups
01214   for ( int ent = MED_CELL; ent < MED_ALL_ENTITIES; ++ent ) {
01215     int nbGroups = _ptrMesh->getNumberOfGroups(medEntityMesh(ent));
01216     nbParts += nbGroups;
01217   }
01218   return nbParts;
01219 }
01220 
01221 // ================================================================================
01222 // RDONLY
01223 // ================================================================================
01224 
01225 ENSIGHT_MESH_RDONLY_DRIVER::ENSIGHT_MESH_RDONLY_DRIVER()
01226   : ENSIGHT_MESH_DRIVER(), _indexInCaseFile(1)
01227 {
01228 }
01229 
01230 ENSIGHT_MESH_RDONLY_DRIVER::ENSIGHT_MESH_RDONLY_DRIVER(const string & fileName,
01231                                                        MESH *         ptrMesh,
01232                                                        int            index)
01233   : ENSIGHT_MESH_DRIVER(fileName,ptrMesh,RDONLY), _indexInCaseFile( index )
01234 {
01235 }
01236 
01237 ENSIGHT_MESH_RDONLY_DRIVER::ENSIGHT_MESH_RDONLY_DRIVER(const ENSIGHT_MESH_RDONLY_DRIVER & driver) : ENSIGHT_MESH_DRIVER(driver), _indexInCaseFile( driver._indexInCaseFile )
01238 {
01239 }
01240 
01241 ENSIGHT_MESH_RDONLY_DRIVER::~ENSIGHT_MESH_RDONLY_DRIVER()
01242 {
01243 }
01244 
01245 GENDRIVER * ENSIGHT_MESH_RDONLY_DRIVER::copy() const
01246 {
01247   return new ENSIGHT_MESH_RDONLY_DRIVER(*this) ;
01248 }
01249 
01250 void ENSIGHT_MESH_RDONLY_DRIVER::write() const throw (MEDEXCEPTION)
01251 {
01252   throw MEDEXCEPTION("ENSIGHT_MESH_RDONLY_DRIVER::write : Can't write with a RDONLY driver !");
01253 }
01254 
01255 void ENSIGHT_MESH_RDONLY_DRIVER::merge ( const GENDRIVER& driver )
01256 {
01257   _CaseFileDriver_User::merge( driver );
01258 
01259   const ENSIGHT_MESH_RDONLY_DRIVER* other =
01260     dynamic_cast< const ENSIGHT_MESH_RDONLY_DRIVER* >( &driver );
01261   if ( other ) {
01262     if ( _indexInCaseFile < other->_indexInCaseFile )
01263       _indexInCaseFile = other->_indexInCaseFile;
01264   }
01265 }
01266 
01267 //================================================================================
01271 //================================================================================
01272 
01273 void ENSIGHT_MESH_RDONLY_DRIVER::read() throw (MEDEXCEPTION)
01274 {
01275   const char * LOC = "ENSIGHT_MESH_RDONLY_DRIVER::read() : " ;
01276   BEGIN_OF_MED(LOC);
01277 
01278   openConst(false); // check if can read case file
01279 
01280   _CaseFileDriver caseFile( getCaseFileName(), this);
01281   caseFile.read();
01282   caseFile.setDataFileName( _indexInCaseFile, this ); // data from Case File is passed here
01283 
01284   openConst(true); // check if can read data file
01285 
01286   cout << "-> Entering into the geometry file " << getDataFileName() << endl  ;
01287 
01288   MESH* mesh = (MESH*) getMesh();
01289 
01290   _InterMed* imed = new _InterMed();
01291   imed->_medMesh = mesh;
01292   imed->_isOwnMedMesh = false;
01293   imed->_needSubParts = ( caseFile.getNbVariables() > 0 );
01294   imed->groupes.reserve(1000);
01295 
01296   // to let field drivers know eventual indices of values
01297   if ( imed->_needSubParts )
01298     setInterData( imed );
01299 
01300   if ( isBinaryDataFile( getDataFileName() )) // binary
01301   {
01302     if ( isGoldFormat() ) // Gold
01303     {
01304       readGoldBinary( *imed );
01305     }
01306     else // EnSight6
01307     {
01308       read6Binary( *imed );
01309     }
01310   }
01311   else // ASCII
01312   {
01313     if ( isGoldFormat() ) // Gold
01314     {
01315       readGoldASCII( *imed );
01316     }
01317     else // EnSight6
01318     {
01319       read6ASCII( *imed );
01320     }
01321   }
01322 
01323   if ( _isMadeByMed && !imed->groupes.empty() ) {
01324     mesh->_name = imed->groupes[0].nom;
01325     imed->groupes[0].nom = "SupportOnAll_";
01326     imed->groupes[0].nom += entNames[MED_CELL];
01327   }
01328   else {
01329     mesh->_name = theDefaultMeshName;
01330   }
01331   mesh->_spaceDimension = SPACE_DIM;
01332   mesh->_numberOfNodes  = imed->points.size() - imed->nbMerged( MED_POINT1 );
01333   mesh->_coordinate     = imed->getCoordinate();
01334 
01335   //Construction des groupes
01336   imed->getGroups(mesh->_groupCell,
01337                   mesh->_groupFace,
01338                   mesh->_groupEdge,
01339                   mesh->_groupNode, mesh);
01340 
01341   mesh->_connectivity = imed->getConnectivity();
01342 
01343   mesh->createFamilies();
01344 
01345   // add attributes to families
01346   set<string> famNames;
01347   for (medEntityMesh entity=MED_CELL; entity<MED_ALL_ENTITIES; ++entity)
01348   {
01349     int i, nb = mesh->getNumberOfFamilies(entity);
01350     for ( i = 1; i <= nb; ++i ) {
01351       FAMILY* f = const_cast<FAMILY*>( mesh->getFamily( entity, i ));
01352       f->setNumberOfAttributes( 1 );
01353       int* attIDs = new int[1];
01354       attIDs[0] = 1;
01355       f->setAttributesIdentifiers( attIDs );
01356       int* attVals = new int[1];
01357       attVals[0] = 1;
01358       f->setAttributesValues( attVals );
01359       string* attDescr = new string[1];
01360       attDescr[0] = "med_family";
01361       f->setAttributesDescriptions( attDescr );
01362       delete [] attDescr;
01363       if ( f->getName().length() > 31 ) // limit a name length
01364         f->setName( STRING("FAM_") << f->getIdentifier());
01365       // setAll() for groups
01366       nb = mesh->getNumberOfGroups(entity);
01367       for ( i = 1; i <= nb; ++i ) {
01368         GROUP * g = const_cast<GROUP*>( mesh->getGroup( entity, i ));
01369         if (mesh->getNumberOfElements( entity, MED_ALL_ELEMENTS ) ==
01370             g->getNumberOfElements( MED_ALL_ELEMENTS ))
01371         {
01372           g->setAll( true );
01373           g->update();
01374         }
01375       }
01376     }
01377   }
01378 
01379   if ( !imed->_needSubParts )
01380     delete imed;
01381 
01382   END_OF_MED(LOC);
01383 }
01384 
01385 //================================================================================
01389 //================================================================================
01390 
01391 void ENSIGHT_MESH_RDONLY_DRIVER::readGoldASCII(_InterMed & imed)
01392 {
01393   const char * LOC ="ENSIGHT_MESH_RDONLY_DRIVER::readGoldASCII() : ";
01394   BEGIN_OF_MED(LOC);
01395 
01396   _ASCIIFileReader geoFile( getDataFileName() );
01397 
01398   if ( isSingleFileMode() ) {
01399     int curTimeStep = 1;
01400     while ( curTimeStep++ < getIndexInDataFile() ) {
01401       while ( !geoFile.isTimeStepEnd())
01402         geoFile.getLine();
01403     }
01404     while ( !geoFile.isTimeStepBeginning() )
01405       geoFile.getLine();
01406   }
01407   // ----------------------
01408   // Read mesh description
01409   // ----------------------
01410   {
01411     string descriptionLine1 = geoFile.getLine();
01412     string descriptionLine2 = geoFile.getLine();
01413 
01414     // find out if the file was created by MED driver
01415     int prefixSize = strlen( theDescriptionPrefix );
01416     _isMadeByMed = ( descriptionLine1.substr(0, prefixSize ) == theDescriptionPrefix );
01417 
01418     if ( _isMadeByMed )
01419       descriptionLine1 = descriptionLine1.substr( prefixSize );
01420     _ptrMesh->setDescription( descriptionLine1 + descriptionLine2 );
01421   }
01422 
01423   // ----------------------------------------
01424   // Find out presence of node/elem numbers 
01425   // ----------------------------------------
01426 
01427   // EnSight User Manual (for v8) says:
01428   //    You do not have to assign node IDs. If you do, the element connectivities are
01429   //    based on the node numbers. If you let EnSight assign the node IDs, the nodes
01430   //    are considered to be sequential starting at node 1, and element connectivity is
01431   //    done accordingly. If node IDs are set to off, they are numbered internally;
01432   //    however, you will not be able to display or query on them. If you have node
01433   //    IDs in your data, you can have EnSight ignore them by specifying "node id
01434   //    ignore." Using this option may reduce some of the memory taken up by the
01435   //    Client and Server, but display and query on the nodes will not be available.
01436 
01437   // read "node|element id <off|given|assign|ignore>"
01438   geoFile.getWord(); geoFile.getWord();
01439   string nodeIds = geoFile.getWord();
01440   geoFile.getWord(); geoFile.getWord();
01441   string elemIds = geoFile.getWord();
01442 
01443   bool haveNodeIds = ( nodeIds == "given" || nodeIds == "ignore" );
01444   bool haveElemIds = ( elemIds == "given" || elemIds == "ignore" );
01445 
01446   // extents: xmin xmax ymin ymax zmin zmax
01447   vector<double> extents;
01448   geoFile.toNextLine();
01449   if ( strncmp( "extents", geoFile.getCurrentPtr(), 7 ) == 0 ) {
01450     geoFile.skip( /*width =*/ 7, /*nbLines =*/ 1 );
01451     extents.reserve( 6 );
01452     while ( extents.size() < extents.capacity() )
01453       extents.push_back( geoFile.getReal() );
01454   }
01455 
01456   typedef map<int,_noeud>::iterator INoeud;
01457   map<int,_noeud> & points = imed.points;
01458   INoeud firstNode;
01459 
01460   _groupe * partGroupe = 0;
01461   int partNum = 0, nbParts = 0;
01462 
01463   while ( !geoFile.isTimeStepEnd() )
01464   {
01465     string word, restLine, line = geoFile.getLine();
01466     TStrTool::split( line, word, restLine );
01467 
01468     const TEnSightElemType & partType = getEnSightType( word );
01469     if ( partType._medType != MED_ALL_ELEMENTS )
01470     {
01471       //  Unstructured element type encounters
01472       // --------------------------------------
01473       int  nbElemNodes = partType._medType % 100;
01474       int      nbElems = geoFile.getInt(); // ne
01475       bool     isGhost = isGhostType( word );
01476       int    nodeShift = points.empty() ? 0 : points.rbegin()->first;
01477 
01478       // read element ids
01479       vector<int> elemIds;
01480       if ( haveElemIds ) {
01481         elemIds.reserve( nbElems );
01482         while ((int) elemIds.size() < nbElems )
01483           elemIds.push_back( geoFile.getInt() ); // id_e
01484       }
01485       if ( isGhost ) { // do not store ghost elements (?)
01486         int nbNodes = nbElems * nbElemNodes;
01487         if ( partType._name == "nsided" ) // polygons
01488         {
01489           for ( int i = 0; i < nbElems; ++i )
01490             nbNodes += geoFile.getInt();
01491           geoFile.skip( nbNodes * INT_WIDTH_GOLD, /*nbLines = */nbElems );
01492         }
01493         else if ( partType._name == "nfaced" ) // polyhedrons
01494         {
01495           int nbFaces = 0;
01496           for ( int i = 0; i < nbElems; ++i )
01497             nbFaces += geoFile.getInt();
01498           for ( int f = 0; f < nbFaces; ++f )
01499             nbNodes += geoFile.getInt();
01500           geoFile.skip( nbNodes * INT_WIDTH_GOLD, /*nbLines = */nbFaces );
01501         }
01502         else // standard types
01503         {
01504           geoFile.skip( nbNodes, nbElemNodes, INT_WIDTH_GOLD );
01505         }
01506         continue;
01507       }
01508 
01509       // add a group corresponding to subPart (geoType)
01510       imed.groupes.push_back(_groupe());
01511       _groupe & groupe = imed.groupes.back();
01512       groupe.mailles.resize( nbElems );
01513 
01514       // find out if "coordinates" has already been encountered
01515       _SubPartDesc coordDesc( partNum , "coordinates");
01516       map< _SubPartDesc, _SubPart >::iterator descPart =
01517         imed._subPartDescribed.find( coordDesc );
01518       bool haveCoords = ( descPart != imed._subPartDescribed.end() );
01519       if ( haveCoords ) {
01520         firstNode = descPart->second.myFirstNode;
01521         nodeShift -= descPart->second.myNbNodes;
01522       }
01523 
01524       // read poly element data
01525       bool isPoly = ( !nbElemNodes );
01526       vector<int> nbElemNodesVec( 1, nbElemNodes);
01527       vector<int> nbElemFaces, nbFaceNodes;
01528       if ( partType._name == "nsided" ) // polygons
01529       {
01530         nbElemNodesVec.resize( nbElems );
01531         for ( int i = 0; i < nbElems; ++i )
01532           nbElemNodesVec[ i ] = geoFile.getInt(); // npi
01533       }
01534       else if ( partType._name == "nfaced" ) // polyhedrons
01535       {
01536         nbElemFaces.resize( nbElems );
01537         nbElemNodesVec.resize( nbElems );
01538         int totalNbFaces = 0;
01539         for ( int i = 0; i < nbElems; ++i )
01540           totalNbFaces += ( nbElemFaces[ i ] = geoFile.getInt() ); // nf_ei
01541 
01542         nbFaceNodes.resize( totalNbFaces );
01543         vector<int>::iterator nbFN = nbFaceNodes.begin();
01544         for ( int i = 0; i < nbElems; ++i ) {
01545           nbElemNodesVec[ i ] = 0;
01546           for ( int nbFaces = nbElemFaces[ i ]; nbFaces; --nbFaces, ++nbFN )
01547             nbElemNodesVec[ i ] += ( *nbFN = geoFile.getInt() ); // np(f_ei)
01548         }
01549       }
01550       // iterator returning nbElemNodes for standard elems and
01551       // next value from nbElemNodesVec for poly elements
01552       _ValueIterator<int> nbElemNodesIt( & nbElemNodesVec[0], isPoly ? 1 : 0);
01553 
01554       // iterator returning values form partType._medIndex for standard elems
01555       // and node index (n) for poly elements
01556       int n;
01557       _ValueIterator<int> medIndexIt( isPoly ? &n : &partType._medIndex[0],
01558                                       isPoly ? 0  : 1);
01559       // read connectivity
01560       _maille ma( partType._medType, nbElemNodes );
01561       ma.sommets.resize( nbElemNodes );
01562       INoeud node;
01563       for ( int i = 0; i < nbElems; ++i ) {
01564         _ValueIterator<int> medIndex = medIndexIt;
01565         nbElemNodes = nbElemNodesIt.next();
01566         if ((int) ma.sommets.size() != nbElemNodes )
01567           ma.sommets.resize( nbElemNodes );
01568         for ( n = 0; n < nbElemNodes; ++n ) {
01569           int nodeID = geoFile.getInt(); // nn_ei
01570           if ( haveCoords )
01571             node = points.find( nodeID + nodeShift );
01572           else
01573             node = points.insert( make_pair( nodeID + nodeShift, _noeud())).first;
01574           ma.sommets[ medIndex.next() ] = node;
01575         }
01576         if ( haveElemIds )
01577           ma.setOrdre( elemIds[ i ] );
01578         groupe.mailles[i] = imed.insert(ma);
01579       }
01580       // store nb nodes in polyhedron faces
01581       if ( !nbFaceNodes.empty() ) {
01582         const int* nbFaceNodesPtr = & nbFaceNodes[0];
01583         for ( int i = 0; i < nbElems; ++i ) {
01584           vector<int> & nbNodesByFace = imed.polyherdalNbFaceNodes[ &groupe.maille( i ) ];
01585           nbNodesByFace.assign( nbFaceNodesPtr, nbFaceNodesPtr + nbElemFaces[ i ] );
01586           nbFaceNodesPtr += nbElemFaces[ i ];
01587         }
01588       }
01589       // create subPart for "coordinates"
01590       if ( !haveCoords ) {
01591         _SubPart & coordSubPart = imed._subPartDescribed[ coordDesc ];
01592         coordSubPart.myFirstNode = points.insert( make_pair( 1 + nodeShift, _noeud())).first;
01593       }
01594       // add subPart group to part group
01595       int groupeIndex = imed.groupes.size();
01596       partGroupe->groupes.push_back( groupeIndex );
01597 
01598       // create subPart
01599       _SubPart subPart( partNum, partType._name );
01600       subPart.myNbCells = nbElems;
01601       subPart.myCellGroupIndex = groupeIndex;
01602       imed.addSubPart( subPart );
01603     }
01604     else if ( word == "coordinates" )
01605     {
01606       // Local node coordinates of a part
01607       // ------------------------------------
01608       int nbNodes = geoFile.getInt(); // nn
01609 
01610       // read node ids
01611       vector<int> nodeIds;
01612       if ( haveNodeIds ) {
01613         nodeIds.reserve( nbNodes );
01614         while ((int) nodeIds.size() < nbNodes )
01615           nodeIds.push_back( geoFile.getInt() ); // id_n
01616       }
01617 
01618       // find out if "coordinates" has already been add at reading connectivity
01619       _SubPartDesc coordDesc( partNum , "coordinates");
01620       map< _SubPartDesc, _SubPart >::iterator descPart =
01621         imed._subPartDescribed.find( coordDesc );
01622       bool haveCoords = ( descPart != imed._subPartDescribed.end() );
01623 
01624       if ( haveCoords ) {
01625         // check that all nodes have been added
01626         firstNode = descPart->second.myFirstNode;
01627         descPart->second.myNbNodes = nbNodes;
01628         INoeud inoeud = firstNode, inoEnd = points.end();
01629         int id = inoeud->first, idEnd = id + nbNodes;
01630         for ( ; id < idEnd; ++id ) {
01631           if ( inoeud == inoEnd || inoeud->first > id ) {
01632             INoeud in = points.insert( inoeud, make_pair( id, _noeud() ));
01633             in->second.number = id;
01634             in->second.coord.resize( SPACE_DIM );
01635           } else {
01636             ++inoeud;
01637           }
01638         }
01639       }
01640       else {
01641         // add nodes
01642         int nodeShift = points.empty() ? 0 : points.rbegin()->first;
01643         for ( int iNode = 1; iNode <= nbNodes; ++iNode ) {
01644           INoeud inoeud = points.insert( points.end(), make_pair( iNode + nodeShift, _noeud()));
01645           inoeud->second.number = inoeud->first;
01646           inoeud->second.coord.resize( SPACE_DIM );
01647         }
01648         firstNode = points.find( 1 + nodeShift );
01649         // create "coordinates" subPart
01650         _SubPart & subPart  = imed._subPartDescribed[ coordDesc ];
01651         subPart.myNbNodes   = nbNodes;
01652         subPart.myFirstNode = firstNode;
01653       }
01654 
01655       // read coordinates in no interlace mode
01656       INoeud endNode = points.end();
01657       for ( int j = 0; j < SPACE_DIM; ++j ) {
01658         for ( INoeud in = firstNode; in != endNode; ++in ) {
01659           _noeud & node = in->second;
01660           node.coord[ j ] = geoFile.getReal();
01661         }
01662       }
01663     }
01664     else if ( word == "part" )
01665     {
01666       // Another part encounters
01667       // -----------------------
01668       partNum = geoFile.getInt();
01669       nbParts++;
01670       geoFile.toNextLine();
01671 
01672       string partName = geoFile.getLine();
01673       if ( partName.empty() )
01674         partName = "Part_" + restLine;
01675 
01676       if (int( imed.groupes.capacity() - imed.groupes.size() ) < theMaxNbTypes )
01677         imed.groupes.reserve( size_t( 1.5 * imed.groupes.size() ));
01678       imed.groupes.push_back(_groupe());
01679       partGroupe = & imed.groupes.back();
01680       partGroupe->nom = partName;
01681       partGroupe->groupes.reserve( theMaxNbTypes );
01682     }
01683     else if ( word == "block" )
01684     {
01685       // Structured type
01686       // ------------------
01687       bool rectilinear = ( restLine.find( "rectilinear" ) != restLine.npos );
01688       bool uniform     = ( restLine.find( "uniform" )     != restLine.npos );
01689       bool curvilinear = ( !rectilinear && !uniform );
01690       bool iblanked    = ( restLine.find( "iblanked" )    != restLine.npos );
01691       bool with_ghost  = ( restLine.find( "with_ghost" )  != restLine.npos );
01692       bool range       = ( restLine.find( "range" )       != restLine.npos );
01693 
01694       // dimension
01695       int I = geoFile.getInt();
01696       int J = geoFile.getInt();
01697       int K = geoFile.getInt();
01698       int NumberOfNodes = I*J*K;
01699       if ( !NumberOfNodes ) continue;
01700 
01701       // range
01702       if ( range ) {
01703         vector<int> ijkRange; // imin imax jmin jmax kmin kmax
01704         ijkRange.reserve(6);
01705         while ( ijkRange.size() < 6 )
01706           ijkRange.push_back( geoFile.getInt() );
01707         I = ijkRange[1]-ijkRange[0]+1;
01708         J = ijkRange[3]-ijkRange[2]+1;
01709         K = ijkRange[5]-ijkRange[4]+1;
01710         NumberOfNodes = I*J*K;
01711       }
01712       // add nodes
01713       int nodeShift = points.empty() ? 0 : points.rbegin()->first;
01714       for ( int iNode = 1; iNode <= NumberOfNodes; ++iNode ) {
01715         INoeud inoeud = points.insert( points.end(), make_pair( iNode + nodeShift, _noeud()));
01716         _noeud & node = inoeud->second;
01717         node.number   = inoeud->first;
01718         node.coord.resize( SPACE_DIM );
01719       }
01720       INoeud firstNode = points.find( nodeShift + 1 );
01721       INoeud endNode   = points.end();
01722 
01723       GRID grid; // calculator of unstructured data
01724       grid._iArrayLength   = I;
01725       grid._jArrayLength   = J;
01726       grid._kArrayLength   = K;
01727       grid._spaceDimension = SPACE_DIM;
01728       if ( J < 2 ) { grid._spaceDimension--; grid._jArrayLength = 0; }
01729       if ( K < 2 ) { grid._spaceDimension--; grid._kArrayLength = 0; }
01730       int nbElems = grid.getNumberOfElements(MED_CELL, MED_ALL_ELEMENTS);
01731 
01732       if ( curvilinear ) // read coordinates for all nodes
01733       {
01734         for ( int j = 0; j < SPACE_DIM; ++j ) {
01735           for ( INoeud in = firstNode; in != endNode; ++in )
01736             in->second.coord[ j ] = geoFile.getReal();
01737         }
01738         grid._gridType   = MED_BODY_FITTED;
01739         grid._coordinate = new COORDINATE(3, 1, 0); // to avoid exception in convertInMESH()
01740       }
01741       else if ( rectilinear ) // read delta vectors with non-regular spacing 
01742       {
01743         grid._iArray = (double*)geoFile.convertReals<double>( I );
01744         grid._jArray = (double*)geoFile.convertReals<double>( J );
01745         grid._kArray = (double*)geoFile.convertReals<double>( K );
01746         grid._gridType = MED_CARTESIAN;
01747       }
01748       else // uniform: read grid origine and delta vectors for regular spacing grid
01749       {
01750         TDblOwner xyzOrigin( (double*)geoFile.convertReals<double>( 3 ));
01751         TDblOwner xyzDelta ( (double*)geoFile.convertReals<double>( 3 ));
01752         // compute full delta vectors
01753         grid._iArray = new double[ I ];
01754         grid._jArray = new double[ J ];
01755         grid._kArray = new double[ K ];
01756         double* coors[SPACE_DIM] = { grid._iArray, grid._jArray, grid._kArray };
01757         int     size [SPACE_DIM] = { I, J, K };
01758         for ( int j = 0; j < SPACE_DIM; ++j ) {
01759           double* coo    = coors[ j ];
01760           double* cooEnd = coo + size[ j ];
01761           coo[0]         = xyzOrigin[ j ];
01762           while ( ++coo < cooEnd )
01763             *coo = coo[-1] + xyzDelta[ j ];
01764         }
01765         grid._gridType = MED_CARTESIAN;
01766       }
01767 
01768       // iblanks
01769       if ( iblanked )
01770         geoFile.skip( NumberOfNodes, /*nbPerLine =*/ 1, INT_WIDTH_GOLD);
01771       // ghosts
01772       if ( with_ghost ) {
01773         geoFile.getWord(); // "ghost_flags"
01774         geoFile.skip( nbElems, /*nbPerLine =*/ 1, INT_WIDTH_GOLD);
01775       }
01776       // node ids
01777       if ( haveNodeIds && geoFile.lookAt( "node_ids" )) {
01778         geoFile.getWord(); // "node_ids"
01779         geoFile.skip( NumberOfNodes, /*nbPerLine =*/ 1, INT_WIDTH_GOLD);
01780       }
01781       // element ids
01782       if ( haveElemIds && geoFile.lookAt( "element_ids" ) ) {
01783         geoFile.getWord(); // "element_ids"
01784         geoFile.skip( nbElems, /*nbPerLine =*/ 1, INT_WIDTH_GOLD);
01785       }
01786 
01787       // let GRID compute all coordinates and connectivity
01788       const MESH* unstruct = grid.convertInMESH();
01789       if ( !curvilinear )
01790       {
01791         const double * coo = unstruct->getCoordinates(MED_FULL_INTERLACE);
01792         typedef _ValueIterator< double > TCoordIt;
01793         TCoordIt xCoo( coo+0, grid._spaceDimension);
01794         TCoordIt yCoo( coo+1, grid._spaceDimension);
01795         TCoordIt zCoo( coo+2, grid._spaceDimension);
01796         if ( grid._spaceDimension < 3 ) zCoo = TCoordIt( grid._kArray, 0 );
01797         if ( grid._spaceDimension < 2 ) yCoo = TCoordIt( grid._jArray, 0 );
01798         for ( INoeud in = firstNode; in != endNode; ++in ) {
01799           _noeud& node = in->second;
01800           node.coord[ 0 ] = xCoo.next();
01801           node.coord[ 1 ] = yCoo.next();
01802           node.coord[ 2 ] = zCoo.next();
01803         }
01804       }
01805 
01806       // store connectivity 
01807       const int * conn = unstruct->getConnectivity( MED_NODAL, MED_CELL, MED_ALL_ELEMENTS );
01808       medGeometryElement elemType = grid.getElementType( MED_CELL, 1 );
01809       int  nbElemNodes = elemType % 100;
01810 
01811       partGroupe->mailles.resize( nbElems );
01812       _maille ma( elemType, nbElemNodes );
01813       ma.sommets.resize( nbElemNodes );
01814 
01815       for ( int i = 0, nIndex = 0; i < nbElems; ++i ) {
01816         for ( int n = 0; n < nbElemNodes; ++n ) {
01817           int nodeID = conn[ nIndex++ ];
01818           ma.sommets[n] = points.find( nodeID + nodeShift );
01819         }
01820         //ma.ordre = ++order;
01821         partGroupe->mailles[i] = imed.insert(ma);
01822       }
01823       _SubPart subPart( partNum, "block" );
01824       subPart.myNbCells    = nbElems;
01825       subPart.myNbNodes    = NumberOfNodes;
01826       subPart.myFirstNode  = firstNode;
01827       subPart.myCellGroupIndex = imed.groupes.size();
01828       imed.addSubPart( subPart );
01829 
01830       unstruct->removeReference();
01831     }
01832     else
01833     {
01834       throw MEDEXCEPTION
01835         ( LOCALIZED( STRING(LOC) << "Unexpected word: " << word <<
01836                      " in " << getDataFileName()));
01837     }
01838   } // while ( !geoFile.eof() )
01839 
01840   if ( nbParts > 1 )
01841     imed.mergeNodesAndElements(TOLERANCE);
01842 
01843   END_OF_MED(LOC);
01844 }
01845 
01846 //================================================================================
01850 //================================================================================
01851 
01852 void ENSIGHT_MESH_RDONLY_DRIVER::readGoldBinary(_InterMed & imed)
01853 {
01854   const char * LOC ="ENSIGHT_MESH_RDONLY_DRIVER::readGoldBinary() : ";
01855   BEGIN_OF_MED(LOC);
01856 
01857   _BinaryFileReader geoFile( getDataFileName() );
01858 
01859   // check if swapping bytes needed
01860   try {
01861     countPartsBinary( geoFile, isSingleFileMode());
01862   }
01863   catch (...) {
01864     geoFile.swapBytes();
01865     geoFile.rewind();
01866   }
01867   if ( getIndexInDataFile() <= 1 )
01868     geoFile.rewind();
01869   if ( geoFile.getPosition() == 0 ) {
01870     TStrOwner format( geoFile.getLine() ); // "C|Fortran Binary"
01871     if ( !contains( "C Binary", format )) {
01872       if ( contains( "Fortran Binary", format ))
01873         throw(MEDEXCEPTION(STRING(LOC) << "Fortran Binary format not supported"));
01874       else
01875         throw(MEDEXCEPTION(STRING(LOC) << "unexpected line in " << getDataFileName()
01876                            << "\n" << format.myValues));
01877     }
01878   }
01879   if ( isSingleFileMode() ) {
01880     // one time step may be skipped by countPartsBinary
01881     int curTimeStep = geoFile.getPosition() ? 2 : 1 ;
01882     while ( curTimeStep < getIndexInDataFile() ) {
01883       countPartsBinary( geoFile, true ); // skip time step
01884       curTimeStep++;
01885     }
01886     while (1) {
01887       TStrOwner line( geoFile.getLine() );
01888       if ( isTimeStepBeginning( line.myValues ))
01889         break;
01890     }
01891   }
01892   // ----------------------
01893   // Read mesh description
01894   // ----------------------
01895   {
01896     TStrOwner descriptionLine1 ( geoFile.getLine() );
01897     TStrOwner descriptionLine2 ( geoFile.getLine() );
01898 
01899     // find out if the file was created by MED driver
01900     _isMadeByMed = contains( theDescriptionPrefix, descriptionLine1 );
01901 
01902     if ( _isMadeByMed )
01903       _ptrMesh->setDescription( descriptionLine2.myValues );
01904     else
01905       _ptrMesh->setDescription( string(descriptionLine1) + descriptionLine2.myValues );
01906   }
01907 
01908   // ----------------------------------------
01909   // Find out presence of node/elem numbers 
01910   // ----------------------------------------
01911 
01912   // EnSight User Manual (for v8) says:
01913   //    You do not have to assign node IDs. If you do, the element connectivities are
01914   //    based on the node numbers. If you let EnSight assign the node IDs, the nodes
01915   //    are considered to be sequential starting at node 1, and element connectivity is
01916   //    done accordingly. If node IDs are set to off, they are numbered internally;
01917   //    however, you will not be able to display or query on them. If you have node
01918   //    IDs in your data, you can have EnSight ignore them by specifying "node id
01919   //    ignore." Using this option may reduce some of the memory taken up by the
01920   //    Client and Server, but display and query on the nodes will not be available.
01921 
01922   // read "node|element id <off|given|assign|ignore>"
01923   bool haveNodeIds, haveElemIds;
01924   {
01925     TStrOwner nodeIds( geoFile.getLine() );
01926     TStrOwner elemIds( geoFile.getLine() );
01927 
01928     haveNodeIds = ( contains( "given", nodeIds ) || contains( "ignore", nodeIds ) );
01929     haveElemIds = ( contains( "given", elemIds ) || contains( "ignore", elemIds ) );
01930   }
01931 
01932   typedef map<int,_noeud>::iterator INoeud;
01933   map<int,_noeud> & points = imed.points;
01934   INoeud firstNode;
01935 
01936   _groupe * partGroupe = 0;
01937   int partNum = 0, nbParts = 0;
01938 
01939   TFltOwner extents(0); // extents: xmin xmax ymin ymax zmin zmax
01940 
01941   while ( !geoFile.eof() )
01942   {
01943     TStrOwner line( geoFile.getLine() );
01944     if ( isSingleFileMode() && isTimeStepEnd( line.myValues ))
01945       break;
01946     string word, restLine;
01947     TStrTool::split( line.myValues, word, restLine );
01948 
01949     const TEnSightElemType & partType = getEnSightType( word );
01950     if ( partType._medType != MED_ALL_ELEMENTS )
01951     {
01952       //  Unstructured element type encounters
01953       // --------------------------------------
01954       int      nbElems = *TIntOwner( geoFile.getInt(1) ); // ne
01955       int  nbElemNodes = partType._medType % 100;
01956       bool     isGhost = isGhostType( word );
01957       int    nodeShift = points.empty() ? 0 : points.rbegin()->first;
01958 
01959       // read element ids
01960       TIntOwner elemIds( haveElemIds ? geoFile.getInt( nbElems ): 0 ); // id_e
01961 
01962       if ( isGhost ) { // do not store ghost elements (?)
01963         int nbNodes = nbElems * nbElemNodes;
01964         if ( partType._name == "nsided" ) // polygons
01965         {
01966           TIntOwner nbNodesInFace( geoFile.getInt( nbElems ));
01967           nbNodes = std::accumulate( nbNodesInFace.myValues, nbNodesInFace.myValues + nbElems, 0 );
01968         }
01969         else if ( partType._name == "nfaced" ) // polyhedrons
01970         {
01971           TIntOwner nbElemFaces( geoFile.getInt( nbElems ));
01972           int nbFaces = accumulate( nbElemFaces.myValues, nbElemFaces.myValues + nbElems, 0 );
01973           TIntOwner nbNodesInFace( geoFile.getInt( nbFaces ));
01974           nbNodes = std::accumulate( nbNodesInFace.myValues, nbNodesInFace.myValues + nbElems, 0 );
01975         }
01976         geoFile.skip( nbNodes * sizeof(int) );
01977         continue;
01978       }
01979 
01980       // add a group corresponding to subPart (geoType)
01981       imed.groupes.push_back(_groupe());
01982       _groupe & groupe = imed.groupes.back();
01983       groupe.mailles.resize( nbElems );
01984 
01985       // find out if "coordinates" has already been encountered
01986       _SubPartDesc coordDesc( partNum , "coordinates");
01987       map< _SubPartDesc, _SubPart >::iterator descPart =
01988         imed._subPartDescribed.find( coordDesc );
01989       bool haveCoords = ( descPart != imed._subPartDescribed.end() );
01990       if ( haveCoords ) {
01991         firstNode = descPart->second.myFirstNode;
01992         nodeShift -= descPart->second.myNbNodes;
01993       }
01994 
01995       // read poly element data
01996       bool isPoly = ( !nbElemNodes );
01997       int nbNodes = 0;
01998       TIntOwner nbElemNodesVec(0), nbElemFaces(0), nbFaceNodes(0);
01999       if ( partType._name == "nsided" ) // polygons
02000       {
02001         nbElemNodesVec.myValues = geoFile.getInt( nbElems ); // npi
02002         nbNodes = accumulate( nbElemNodesVec.myValues, nbElemNodesVec.myValues + nbElems, 0 );
02003       }
02004       else if ( partType._name == "nfaced" ) // polyhedrons
02005       {
02006         nbElemFaces.myValues = geoFile.getInt( nbElems ); // nf_ei
02007         int totalNbFaces = accumulate( nbElemFaces.myValues, nbElemFaces.myValues + nbElems, 0 );
02008 
02009         nbFaceNodes.myValues = geoFile.getInt( totalNbFaces ); // np(f_ei)
02010         // calculate nb of nodes in each polyhedron
02011         int* nbEN = nbElemNodesVec.myValues = new int[ nbElems ];
02012         const int *nbFN = nbFaceNodes, *nbEF = nbElemFaces, *nbEND = nbEN + nbElems;
02013         for ( ; nbEN < nbEND; ++nbEN, ++nbEF ) {
02014           nbNodes += *nbEN = accumulate( nbFN, nbFN + *nbEF, 0 );
02015           nbFN    += *nbEF;
02016         }
02017       }
02018       else // standard types
02019       {
02020         nbElemNodesVec.myValues = new int[ 1 ];
02021         nbElemNodesVec[ 0 ] = nbElemNodes;
02022         nbNodes = nbElems * nbElemNodes;
02023       }
02024       // iterator returning nbElemNodes for standard elems and
02025       // next value from nbElemNodesVec for poly elements
02026       _ValueIterator<int> nbElemNodesIt( nbElemNodesVec, isPoly ? 1 : 0);
02027 
02028       // iterator returning values form partType._medIndex for standard elems
02029       // and node index (n) for poly elements
02030       int n;
02031       _ValueIterator<int> medIndexIt( isPoly ? &n : &partType._medIndex[0],
02032                                       isPoly ? 0  : 1);
02033       // read connectivity
02034       _maille ma( partType._medType, nbElemNodes );
02035       ma.sommets.resize( nbElemNodes );
02036       TIntOwner connectivity( geoFile.getInt( nbNodes )); // nn_ei
02037       int* nodeID = connectivity;
02038       INoeud node;
02039       for ( int i = 0; i < nbElems; ++i ) {
02040         _ValueIterator<int> medIndex = medIndexIt;
02041         nbElemNodes = nbElemNodesIt.next();
02042         if ((int) ma.sommets.size() != nbElemNodes )
02043           ma.sommets.resize( nbElemNodes );
02044         for ( n = 0; n < nbElemNodes; ++n, ++nodeID ) {
02045           if ( haveCoords )
02046             node = points.find( *nodeID + nodeShift );
02047           else
02048             node = points.insert( make_pair( *nodeID + nodeShift, _noeud())).first;
02049           ma.sommets[ medIndex.next() ] = node;
02050         }
02051         if ( haveElemIds )
02052           ma.setOrdre( elemIds[ i ] );
02053         groupe.mailles[i] = imed.insert(ma);
02054       }
02055       // store nb nodes in polyhedron faces
02056       if ( nbFaceNodes.myValues ) {
02057         const int* nbFaceNodesPtr = nbFaceNodes.myValues;
02058         for ( int i = 0; i < nbElems; ++i ) {
02059           vector<int> & nbNodesByFace = imed.polyherdalNbFaceNodes[ &groupe.maille( i ) ];
02060           nbNodesByFace.assign( nbFaceNodesPtr, nbFaceNodesPtr + nbElemFaces[ i ] );
02061           nbFaceNodesPtr += nbElemFaces[ i ];
02062         }
02063       }
02064       // create subPart for "coordinates"
02065       if ( !haveCoords ) {
02066         _SubPart & coordSubPart = imed._subPartDescribed[ coordDesc ];
02067         coordSubPart.myFirstNode = points.insert( make_pair( 1 + nodeShift, _noeud())).first;
02068       }
02069       // add subPart group to part group
02070       int groupeIndex = imed.groupes.size();
02071       partGroupe->groupes.push_back( groupeIndex );
02072 
02073       // create subPart
02074       _SubPart subPart( partNum, partType._name );
02075       subPart.myNbCells = nbElems;
02076       subPart.myCellGroupIndex = groupeIndex;
02077       imed.addSubPart( subPart );
02078     }
02079     else if ( word == "coordinates" )
02080     {
02081       // Local node coordinates of a part
02082       // ------------------------------------
02083       int nbNodes = *TIntOwner( geoFile.getInt(1) ); // nn
02084 
02085       // read node ids
02086       TIntOwner nodeIds(0);
02087       if ( haveNodeIds )
02088         nodeIds.myValues = geoFile.getInt( nbNodes ); // id_n
02089 
02090       // find out if "coordinates" has already been add at reading connectivity
02091       _SubPartDesc coordDesc( partNum , "coordinates");
02092       map< _SubPartDesc, _SubPart >::iterator descPart =
02093         imed._subPartDescribed.find( coordDesc );
02094       bool haveCoords = ( descPart != imed._subPartDescribed.end() );
02095 
02096       if ( haveCoords ) {
02097         // check that all nodes have been added
02098         firstNode = descPart->second.myFirstNode;
02099         descPart->second.myNbNodes = nbNodes;
02100         INoeud inoeud = firstNode, inoEnd = points.end();
02101         int id = inoeud->first, idEnd = id + nbNodes;
02102         for ( ; id < idEnd; ++id ) {
02103           if ( inoeud == inoEnd || inoeud->first > id ) {
02104             INoeud in = points.insert( inoeud, make_pair( id, _noeud() ));
02105             in->second.number = id;
02106           } else {
02107             ++inoeud;
02108           }
02109         }
02110       }
02111       else {
02112         // add nodes
02113         int nodeShift = points.empty() ? 0 : points.rbegin()->first;
02114         for ( int iNode = 1; iNode <= nbNodes; ++iNode ) {
02115           INoeud inoeud = points.insert( points.end(), make_pair( iNode + nodeShift, _noeud()));
02116           inoeud->second.number = inoeud->first;
02117         }
02118         firstNode = points.find( 1 + nodeShift );
02119         // create "coordinates" subPart
02120         _SubPart & subPart  = imed._subPartDescribed[ coordDesc ];
02121         subPart.myNbNodes   = nbNodes;
02122         subPart.myFirstNode = firstNode;
02123       }
02124 
02125       // read coordinates in no interlace mode
02126       TFltOwner noInterlaceCoords( geoFile.getFlt( nbNodes * SPACE_DIM ));
02127       float* x = noInterlaceCoords;
02128       float* y = x + nbNodes;
02129       float* z = y + nbNodes;
02130       INoeud endNode = points.end();
02131       for ( INoeud in = firstNode; in != endNode; ++in ) {
02132         _noeud & node = in->second;
02133         node.coord.resize( SPACE_DIM );
02134         node.coord[ 0 ] = *x++;
02135         node.coord[ 1 ] = *y++;
02136         node.coord[ 2 ] = *z++;
02137       }
02138     }
02139     else if ( word == "part" )
02140     {
02141       // Another part encounters
02142       // -----------------------
02143       partNum = *TIntOwner( geoFile.getInt(1) );
02144       nbParts++;
02145 
02146       string partName( TStrOwner( geoFile.getLine() ));
02147       if ( partName.empty() )
02148         partName = "Part_" + restLine;
02149 
02150       if (int( imed.groupes.capacity() - imed.groupes.size()) < theMaxNbTypes )
02151         imed.groupes.reserve( size_t( 1.5 * imed.groupes.size() ));
02152       imed.groupes.push_back(_groupe());
02153       partGroupe = & imed.groupes.back();
02154       partGroupe->nom = partName;
02155       partGroupe->groupes.reserve( theMaxNbTypes );
02156     }
02157     else if ( word == "block" )
02158     {
02159       // Structured type
02160       // ------------------
02161       bool rectilinear = ( restLine.find( "rectilinear" ) != restLine.npos );
02162       bool uniform     = ( restLine.find( "uniform" )     != restLine.npos );
02163       bool curvilinear = ( !rectilinear && !uniform );
02164       bool iblanked    = ( restLine.find( "iblanked" )    != restLine.npos );
02165       bool with_ghost  = ( restLine.find( "with_ghost" )  != restLine.npos );
02166       bool range       = ( restLine.find( "range" )       != restLine.npos );
02167 
02168       // dimension
02169       TIntOwner ijk( geoFile.getInt(3) );
02170       int I = ijk[0];
02171       int J = ijk[1];
02172       int K = ijk[2];
02173       int NumberOfNodes = I*J*K;
02174       if ( !NumberOfNodes ) continue;
02175 
02176       // range
02177       if ( range ) {
02178         TIntOwner ijkRange( geoFile.getInt( 6 ));// imin imax jmin jmax kmin kmax
02179         I = ijkRange[1]-ijkRange[0]+1;
02180         J = ijkRange[3]-ijkRange[2]+1;
02181         K = ijkRange[5]-ijkRange[4]+1;
02182         NumberOfNodes = I*J*K;
02183       }
02184       // add nodes
02185       int nodeShift = points.empty() ? 0 : points.rbegin()->first;
02186       for ( int iNode = 1; iNode <= NumberOfNodes; ++iNode ) {
02187         INoeud inoeud = points.insert( points.end(), make_pair( iNode + nodeShift, _noeud()));
02188         _noeud & node = inoeud->second;
02189         node.number   = inoeud->first;
02190         node.coord.resize( SPACE_DIM );
02191       }
02192       INoeud firstNode = points.find( nodeShift + 1 );
02193       INoeud endNode   = points.end();
02194 
02195       GRID grid; // calculator of unstructured data
02196       grid._iArrayLength   = I;
02197       grid._jArrayLength   = J;
02198       grid._kArrayLength   = K;
02199       grid._spaceDimension = SPACE_DIM;
02200       if ( J < 2 ) { grid._spaceDimension--; grid._jArrayLength = 0; }
02201       if ( K < 2 ) { grid._spaceDimension--; grid._kArrayLength = 0; }
02202       int nbElems = grid.getNumberOfElements(MED_CELL, MED_ALL_ELEMENTS);
02203 
02204       if ( curvilinear ) // read coordinates for all nodes
02205       {
02206         TFltOwner noInterlaceCoords( geoFile.getFlt( NumberOfNodes * SPACE_DIM ));
02207         float* x = noInterlaceCoords;
02208         float* y = x + NumberOfNodes;
02209         float* z = y + NumberOfNodes;
02210         for ( INoeud in = firstNode; in != endNode; ++in ) {
02211           _noeud & node = in->second;
02212           node.coord.resize( SPACE_DIM );
02213           node.coord[ 0 ] = *x++;
02214           node.coord[ 1 ] = *y++;
02215           node.coord[ 2 ] = *z++;
02216         }
02217         grid._gridType   = MED_BODY_FITTED;
02218         grid._coordinate = new COORDINATE(3, 1, 0); // to avoid exception in convertInMESH()
02219       }
02220       else if ( rectilinear ) // read delta vectors with non-regular spacing 
02221       {
02222         grid._iArray = (double*)geoFile.convertReals<double>( I );
02223         grid._jArray = (double*)geoFile.convertReals<double>( J );
02224         grid._kArray = (double*)geoFile.convertReals<double>( K );
02225         grid._gridType = MED_CARTESIAN;
02226       }
02227       else // uniform: read grid origine and delta vectors for regular spacing grid
02228       {
02229         TFltOwner xyzOrigin( geoFile.getFlt( 3 ));
02230         TFltOwner xyzDelta ( geoFile.getFlt( 3 ));
02231         // compute full delta vectors
02232         grid._iArray = new double[ I ];
02233         grid._jArray = new double[ J ];
02234         grid._kArray = new double[ K ];
02235         double* coors[SPACE_DIM] = { grid._iArray, grid._jArray, grid._kArray };
02236         int     size [SPACE_DIM] = { I, J, K };
02237         for ( int j = 0; j < SPACE_DIM; ++j ) {
02238           double* coo    = coors[ j ];
02239           double* cooEnd = coo + size[ j ];
02240           coo[0]         = xyzOrigin[ j ];
02241           while ( ++coo < cooEnd )
02242             *coo = coo[-1] + xyzDelta[ j ];
02243         }
02244         grid._gridType = MED_CARTESIAN;
02245       }
02246 
02247       // iblanks
02248       if ( iblanked )
02249         geoFile.skip( NumberOfNodes * sizeof(int) );
02250       // ghosts
02251       if ( with_ghost ) {
02252         TStrOwner( geoFile.getLine() ); // "ghost_flags"
02253         geoFile.skip( nbElems * sizeof(int) );
02254       }
02255       // node ids
02256       if ( haveNodeIds && !geoFile.eof() ) {
02257         TStrOwner nextLine( geoFile.getLine() ); // "node_ids"
02258         if ( contains( "node_ids", nextLine ) )
02259           geoFile.skip( NumberOfNodes * sizeof(int) );
02260         else
02261           geoFile.skip( -MAX_LINE_LENGTH );
02262       }
02263       // element ids
02264       TIntOwner elemIdOwner(0);
02265       _ValueIterator<int> elemIds;
02266       if ( haveElemIds && !geoFile.eof() ) {
02267         TStrOwner nextLine( geoFile.getLine() ); // "element_ids"
02268         if ( contains( "element_ids", nextLine ) ) {
02269           elemIdOwner.myValues = geoFile.getInt( nbElems );
02270           elemIds = _ValueIterator<int>( elemIdOwner, 1);
02271         } else {
02272           geoFile.skip( -MAX_LINE_LENGTH );
02273         }
02274       }
02275 
02276       // let GRID compute all coordinates and connectivity
02277       const MESH* unstruct = grid.convertInMESH();
02278       if ( !curvilinear )
02279       {
02280         const double * coo = unstruct->getCoordinates(MED_FULL_INTERLACE);
02281         typedef _ValueIterator< double > TCoordIt;
02282         TCoordIt xCoo( coo+0, grid._spaceDimension);
02283         TCoordIt yCoo( coo+1, grid._spaceDimension);
02284         TCoordIt zCoo( coo+2, grid._spaceDimension);
02285         if ( grid._spaceDimension < 3 ) zCoo = TCoordIt( grid._kArray, 0 );
02286         if ( grid._spaceDimension < 2 ) yCoo = TCoordIt( grid._jArray, 0 );
02287         for ( INoeud in = firstNode; in != endNode; ++in ) {
02288           _noeud& node = in->second;
02289           node.coord[ 0 ] = xCoo.next();
02290           node.coord[ 1 ] = yCoo.next();
02291           node.coord[ 2 ] = zCoo.next();
02292         }
02293       }
02294 
02295       // store connectivity 
02296       const int * conn = unstruct->getConnectivity( MED_NODAL, MED_CELL, MED_ALL_ELEMENTS );
02297       medGeometryElement elemType = grid.getElementType( MED_CELL, 1 );
02298       int  nbElemNodes = elemType % 100;
02299 
02300       partGroupe->mailles.resize( nbElems );
02301       _maille ma( elemType, nbElemNodes );
02302       ma.sommets.resize( nbElemNodes );
02303 
02304       for ( int i = 0, nIndex = 0; i < nbElems; ++i ) {
02305         for ( int n = 0; n < nbElemNodes; ++n ) {
02306           int nodeID = conn[ nIndex++ ];
02307           ma.sommets[n] = points.find( nodeID + nodeShift );
02308         }
02309         ma.setOrdre( elemIds.next() );
02310         partGroupe->mailles[i] = imed.insert(ma);
02311       }
02312 
02313       _SubPart subPart( partNum, "block" );
02314       subPart.myNbCells    = nbElems;
02315       subPart.myNbNodes    = NumberOfNodes;
02316       subPart.myFirstNode  = firstNode;
02317       subPart.myCellGroupIndex = imed.groupes.size();
02318       imed.addSubPart( subPart );
02319 
02320       unstruct->removeReference();
02321     }
02322     else if ( word == "extents" )
02323     {
02324       extents.myValues = geoFile.getFlt( 6 );
02325     }
02326     else
02327     {
02328       throw MED_EXCEPTION
02329         ( LOCALIZED( STRING(LOC) << "Unexpected word: " << word <<
02330                      " in " << getDataFileName()));
02331     }
02332   } // while ( !geoFile.eof() )
02333 
02334   if ( nbParts > 1 )
02335     imed.mergeNodesAndElements(TOLERANCE);
02336 
02337   END_OF_MED(LOC);
02338 }
02339 
02340 //================================================================================
02344 //================================================================================
02345 
02346 void ENSIGHT_MESH_RDONLY_DRIVER::read6ASCII(_InterMed & imed)
02347 {
02348   const char * LOC ="ENSIGHT_MESH_RDONLY_DRIVER::read6ASCII() : ";
02349   BEGIN_OF_MED(LOC);
02350 
02351   _ASCIIFileReader geoFile( getDataFileName() );
02352 
02353   if ( isSingleFileMode() ) {
02354     int curTimeStep = 1;
02355     while ( curTimeStep < getIndexInDataFile() ) {
02356       while ( !geoFile.isTimeStepEnd())
02357         geoFile.getLine();
02358       curTimeStep++;
02359     }
02360     while ( !geoFile.isTimeStepBeginning() )
02361       geoFile.getLine();
02362   }
02363   // ----------------------
02364   // Read mesh description
02365   // ----------------------
02366   {
02367     string descriptionLine1 = geoFile.getLine();
02368     string descriptionLine2 = geoFile.getLine();
02369 
02370     // find out if the file was created by MED driver
02371     int prefixSize = strlen( theDescriptionPrefix );
02372     _isMadeByMed = ( descriptionLine1.substr(0, prefixSize ) == theDescriptionPrefix );
02373 
02374     if ( _isMadeByMed )
02375       descriptionLine1 = descriptionLine1.substr( prefixSize );
02376     _ptrMesh->setDescription( descriptionLine1 + descriptionLine2 );
02377   }
02378 
02379   // ----------------------------------------
02380   // Find out presence of node/elem numbers 
02381   // ----------------------------------------
02382 
02383   // EnSight User Manual (for v8) says:
02384   //    You do not have to assign node IDs. If you do, the element connectivities are
02385   //    based on the node numbers. If you let EnSight assign the node IDs, the nodes
02386   //    are considered to be sequential starting at node 1, and element connectivity is
02387   //    done accordingly. If node IDs are set to off, they are numbered internally;
02388   //    however, you will not be able to display or query on them. If you have node
02389   //    IDs in your data, you can have EnSight ignore them by specifying "node id
02390   //    ignore." Using this option may reduce some of the memory taken up by the
02391   //    Client and Server, but display and query on the nodes will not be available.
02392 
02393   // read "node|element id <off|given|assign|ignore>"
02394   geoFile.getWord(); geoFile.getWord();
02395   string nodeIds = geoFile.getWord();
02396   geoFile.getWord(); geoFile.getWord();
02397   string elemIds = geoFile.getWord();
02398 
02399   bool haveNodeIds = ( nodeIds == "given" || nodeIds == "ignore" );
02400   bool haveElemIds = ( elemIds == "given" || elemIds == "ignore" );
02401 
02402   map<int,_noeud> & points = imed.points;
02403   typedef map<int,_noeud>::iterator INoeud;
02404 
02405   int haveStructuredParts = 0, haveUnstructuredParts = 0;
02406 
02407   _groupe * partGroupe = 0;
02408   int       partNum = 0;
02409 
02410   while ( !geoFile.isTimeStepEnd() )
02411   {
02412     string word, restLine, line = geoFile.getLine();
02413     TStrTool::split( line, word, restLine );
02414 
02415     const TEnSightElemType & partType = getEnSightType( word );
02416     if ( !partType._medIndex.empty() )
02417     {
02418       //  Unstructured element type encounters
02419       // --------------------------------------
02420       int  nbElemNodes = partType._medType % 100;
02421       int      nbElems = geoFile.getInt();
02422       if ( nbElems > 0 )
02423         haveUnstructuredParts++;
02424 
02425       imed.groupes.push_back(_groupe());
02426       _groupe & groupe = imed.groupes.back();
02427       groupe.mailles.resize( nbElems );
02428 
02429       // read connectivity
02430       _maille ma( partType._medType, nbElemNodes );
02431       ma.sommets.resize( nbElemNodes );
02432       INoeud node;
02433       for ( int i = 0; i < nbElems; ++i ) {
02434         if ( haveElemIds )
02435           geoFile.getInt();
02436         for ( int n = 0; n < nbElemNodes; ++n ) {
02437           int nodeID = geoFile.getInt();
02438           ma.sommets[ partType._medIndex[n] ] = points.find( nodeID );
02439         }
02440         //ma.ordre = ++order;
02441         groupe.mailles[i] = imed.insert(ma);
02442       }
02443 
02444       int groupeIndex = imed.groupes.size();
02445       partGroupe->groupes.push_back( groupeIndex );
02446 
02447       _SubPart subPart( partNum, partType._name );
02448       subPart.myNbCells = nbElems;
02449       subPart.myCellGroupIndex = groupeIndex;
02450       imed.addSubPart( subPart );
02451     }
02452     else if ( word == "part" )
02453     {
02454       // Another part encounters
02455       // -----------------------
02456       partNum = atoi( restLine.c_str() );
02457 
02458       string partName = geoFile.getLine();
02459       if ( partName.empty() )
02460         partName = "Part_" + restLine;
02461 
02462       if (int( imed.groupes.capacity() - imed.groupes.size()) < theMaxNbTypes )
02463         imed.groupes.reserve( size_t( 1.5 * imed.groupes.size() ));
02464       imed.groupes.push_back(_groupe());
02465       partGroupe = & imed.groupes.back();
02466       partGroupe->nom = partName;
02467       partGroupe->groupes.reserve( theMaxNbTypes );
02468     }
02469     else if ( word == "block" )
02470     {
02471       // Structured type
02472       // ------------------
02473       bool iblanked  = ( restLine == "iblanked" );
02474 
02475       // dimension
02476       int I = geoFile.getInt();
02477       int J = geoFile.getInt();
02478       int K = geoFile.getInt();
02479       int NumberOfNodes = I*J*K;
02480       if ( !NumberOfNodes ) continue;
02481       haveStructuredParts++;
02482 
02483       // add nodes
02484       int nodeShift = points.empty() ? 0 : points.rbegin()->first;
02485       for ( int iNode = 1; iNode <= NumberOfNodes; ++iNode ) {
02486         INoeud inoeud = points.insert( points.end(), make_pair( iNode + nodeShift, _noeud()));
02487         _noeud & node = inoeud->second;
02488         node.number   = inoeud->first;
02489         node.coord.resize( SPACE_DIM );
02490       }
02491       // read coordinates
02492       INoeud firstNode = points.find( nodeShift + 1 );
02493       INoeud endNode   = points.end();
02494       for ( int j = 0; j < SPACE_DIM; ++j ) {
02495         for ( INoeud in = firstNode; in != endNode; ++in ) {
02496           _noeud & node = in->second;
02497           node.coord[ j ] = geoFile.getReal();
02498         }
02499       }
02500       // iblanks
02501       if ( iblanked )
02502         geoFile.skip(NumberOfNodes, /*nbPerLine =*/ 10, INT_WIDTH_6);
02503 
02504       // let GRID calculate connectivity 
02505       GRID grid;
02506       grid._iArrayLength  = I;
02507       grid._jArrayLength  = J;
02508       grid._kArrayLength  = K;
02509       grid._gridType      = MED_BODY_FITTED;
02510       grid._spaceDimension= SPACE_DIM;
02511       grid._coordinate    = new COORDINATE(3, 1, 0); // to avoid exception in convertInMESH()
02512       if ( J < 2 ) { grid._spaceDimension--; grid._jArrayLength = 0; }
02513       if ( K < 2 ) { grid._spaceDimension--; grid._kArrayLength = 0; }
02514 
02515       const MESH* unstruct = grid.convertInMESH();
02516       const int * conn = unstruct->getConnectivity( MED_NODAL, MED_CELL, MED_ALL_ELEMENTS );
02517       medGeometryElement elemType = grid.getElementType( MED_CELL, 1 );
02518       int  nbElemNodes = elemType % 100;
02519       int      nbElems = grid.getNumberOfElements(MED_CELL, MED_ALL_ELEMENTS);
02520 
02521       partGroupe->mailles.resize( nbElems );
02522       _maille ma( elemType, nbElemNodes );
02523       ma.sommets.resize( nbElemNodes );
02524 
02525       for ( int i = 0, nIndex = 0; i < nbElems; ++i ) {
02526         for ( int n = 0; n < nbElemNodes; ++n ) {
02527           int nodeID = conn[ nIndex++ ];
02528           ma.sommets[n] = points.find( nodeID + nodeShift );
02529         }
02530         //ma.ordre = ++order;
02531         partGroupe->mailles[i] = imed.insert(ma);
02532       }
02533 
02534       _SubPart subPart( partNum, "block" );
02535       subPart.myNbCells    = nbElems;
02536       subPart.myNbNodes    = NumberOfNodes;
02537       subPart.myFirstNode  = firstNode;
02538       subPart.myCellGroupIndex = imed.groupes.size();
02539       imed.addSubPart( subPart );
02540 
02541       unstruct->removeReference();
02542     }
02543     else if ( word == "coordinates" )
02544     {
02545       // ------------------------------------
02546       // Unstructured global node coordinates
02547       // ------------------------------------
02548       int nbNodes = geoFile.getInt();
02549 
02550       cout << "-> loading coordinates of " << nbNodes << " nodes " << endl ;
02551 
02552       INoeud inoeud;
02553       for ( int i=0 ; i < nbNodes ; i++ )
02554       {
02555         if ( haveNodeIds ) {
02556           int nodeID = geoFile.getInt();
02557           inoeud = points.insert( make_pair( nodeID, _noeud() )).first;
02558           inoeud->second.number = nodeID;
02559         }
02560         else {
02561           int nodeID = i + 1;
02562           inoeud = points.insert( points.end(), make_pair( nodeID, _noeud()));
02563           inoeud->second.number = nodeID;
02564         }
02565         _noeud & node = inoeud->second;
02566         node.coord.resize( SPACE_DIM );
02567         node.coord[ 0 ] = geoFile.getReal();
02568         node.coord[ 1 ] = geoFile.getReal();
02569         node.coord[ 2 ] = geoFile.getReal();
02570       }
02571 
02572       _SubPartDesc cooDesc = _SubPartDesc::globalCoordDesc();
02573       _SubPart subPart( cooDesc.partNumber(), cooDesc.typeName() );
02574       subPart.myNbNodes    = nbNodes;
02575       subPart.myFirstNode  = points.begin();
02576       imed.addSubPart( subPart );
02577     }
02578     else
02579     {
02580       throw MED_EXCEPTION
02581         ( LOCALIZED( STRING(LOC) << "Unexpected word: " << word <<
02582                      " in " << getDataFileName()));
02583     }
02584   } // while ( !geoFile.eof() )
02585 
02586   if ( ( haveStructuredParts && haveUnstructuredParts ) || haveStructuredParts > 1 )
02587     imed.mergeNodesAndElements(TOLERANCE);
02588 
02589   END_OF_MED(LOC);
02590 }
02591 
02592 //================================================================================
02596 //================================================================================
02597 
02598 void ENSIGHT_MESH_RDONLY_DRIVER::read6Binary(_InterMed & imed)
02599 {
02600   const char * LOC ="ENSIGHT_MESH_RDONLY_DRIVER::read6Binary() : ";
02601   BEGIN_OF_MED(LOC);
02602 
02603   _BinaryFileReader geoFile( getDataFileName() );
02604 
02605   // check if swapping bytes needed
02606   try {
02607     countPartsBinary( geoFile, isSingleFileMode());
02608   }
02609   catch (...) {
02610     geoFile.swapBytes();
02611     geoFile.rewind();
02612   }
02613   if ( getIndexInDataFile() <= 1 )
02614     geoFile.rewind();
02615   if ( geoFile.getPosition() == 0 ) {
02616     TStrOwner format( geoFile.getLine() ); // "C|Fortran Binary"
02617     if ( !contains( "C Binary", format )) {
02618       if ( contains( "Fortran Binary", format ))
02619         throw(MEDEXCEPTION(STRING(LOC) << "Fortran Binary format not supported"));
02620       else
02621         throw(MEDEXCEPTION(STRING(LOC) << "unexpected line in " << getDataFileName()
02622                            << "\n" << format.myValues));
02623     }
02624   }
02625 
02626   if ( isSingleFileMode() ) {
02627     // one time step may be skipped by countPartsBinary
02628     int curTimeStep = geoFile.getPosition() ? 2 : 1 ;
02629     while ( curTimeStep < getIndexInDataFile() ) {
02630       countPartsBinary( geoFile, true ); // skip time step
02631       curTimeStep++;
02632     }
02633     while (1) {
02634       TStrOwner line( geoFile.getLine() );
02635       if ( isTimeStepBeginning( line.myValues ))
02636         break;
02637     }
02638   }
02639   // ----------------------
02640   // Read mesh description
02641   // ----------------------
02642   {
02643     TStrOwner descriptionLine1( geoFile.getLine() );
02644     TStrOwner descriptionLine2( geoFile.getLine() );
02645 
02646     // find out if the file was created by MED driver
02647     _isMadeByMed = contains( theDescriptionPrefix, descriptionLine1 );
02648 
02649     if ( _isMadeByMed )
02650       _ptrMesh->setDescription( descriptionLine2.myValues );
02651     else
02652       _ptrMesh->setDescription( string(descriptionLine1) + descriptionLine2.myValues );
02653   }
02654 
02655   // ----------------------------------------
02656   // Find out presence of node/elem numbers 
02657   // ----------------------------------------
02658 
02659   // EnSight User Manual (for v8) says:
02660   //    You do not have to assign node IDs. If you do, the element connectivities are
02661   //    based on the node numbers. If you let EnSight assign the node IDs, the nodes
02662   //    are considered to be sequential starting at node 1, and element connectivity is
02663   //    done accordingly. If node IDs are set to off, they are numbered internally;
02664   //    however, you will not be able to display or query on them. If you have node
02665   //    IDs in your data, you can have EnSight ignore them by specifying "node id
02666   //    ignore." Using this option may reduce some of the memory taken up by the
02667   //    Client and Server, but display and query on the nodes will not be available.
02668 
02669   // read "node|element id <off|given|assign|ignore>"
02670   bool haveNodeIds, haveElemIds;
02671   {
02672     TStrOwner nodeIds( geoFile.getLine() );
02673     TStrOwner elemIds( geoFile.getLine() );
02674 
02675     haveNodeIds = ( contains( "given", nodeIds ) || contains( "ignore", nodeIds ) );
02676     haveElemIds = ( contains( "given", elemIds ) || contains( "ignore", elemIds ) );
02677   }
02678   map<int,_noeud> & points = imed.points;
02679   typedef map<int,_noeud>::iterator INoeud;
02680 
02681   int haveStructuredParts = 0, haveUnstructuredParts = 0;
02682 
02683   _groupe * partGroupe = 0;
02684   int       partNum = 0;
02685 
02686   while ( !geoFile.eof() )
02687   {
02688     TStrOwner line( geoFile.getLine() );
02689     if ( isSingleFileMode() && isTimeStepEnd( line.myValues ))
02690       break;
02691     string word, restLine;
02692     TStrTool::split( line.myValues, word, restLine );
02693 
02694     const TEnSightElemType & partType = getEnSightType( word );
02695     if ( !partType._medIndex.empty() )
02696     {
02697       //  Unstructured element type encounters
02698       // --------------------------------------
02699       int  nbElemNodes = partType._medType % 100;
02700       int      nbElems = *TIntOwner( geoFile.getInt(1) ); // ne
02701       if ( nbElems > 0 )
02702         haveUnstructuredParts++;
02703 
02704       TIntOwner numbers(0);
02705       if ( haveElemIds )
02706         numbers.myValues = geoFile.getInt( nbElems ); // id_e
02707 
02708       imed.groupes.push_back(_groupe());
02709       _groupe & groupe = imed.groupes.back();
02710       groupe.mailles.resize( nbElems );
02711 
02712       // read connectivity
02713       _maille ma( partType._medType, nbElemNodes );
02714       ma.sommets.resize( nbElemNodes );
02715       TIntOwner connectivity( geoFile.getInt( nbElems * nbElemNodes ));
02716       int* nodeID = connectivity;
02717       INoeud node;
02718       for ( int i = 0; i < nbElems; ++i ) {
02719         for ( int n = 0; n < nbElemNodes; ++n, ++nodeID )
02720           ma.sommets[ partType._medIndex[n] ] = points.find( *nodeID );
02721         //ma.ordre = ++order;
02722         groupe.mailles[i] = imed.insert(ma);
02723       }
02724 
02725       int groupeIndex = imed.groupes.size();
02726       partGroupe->groupes.push_back( groupeIndex );
02727 
02728       _SubPart subPart( partNum, partType._name );
02729       subPart.myNbCells = nbElems;
02730       subPart.myCellGroupIndex = groupeIndex;
02731       imed.addSubPart( subPart );
02732     }
02733     else if ( word == "part" )
02734     {
02735       // Another part encounters
02736       // -----------------------
02737       partNum = atoi( restLine.c_str() );
02738 
02739       string partName( TStrOwner( geoFile.getLine() ));
02740       if ( partName.empty() )
02741         partName = "Part_" + restLine;
02742 
02743       if (int( imed.groupes.capacity() - imed.groupes.size()) < theMaxNbTypes )
02744         imed.groupes.reserve( size_t( 1.5 * imed.groupes.size() ));
02745       imed.groupes.push_back(_groupe());
02746       partGroupe = & imed.groupes.back();
02747       partGroupe->nom = partName;
02748       partGroupe->groupes.reserve( theMaxNbTypes );
02749     }
02750     else if ( word == "block" )
02751     {
02752       // Structured type
02753       // ------------------
02754       bool iblanked  = ( restLine == "iblanked" );
02755 
02756       // dimension
02757       TIntOwner ijk( geoFile.getInt(3) );
02758       int I = ijk[0];
02759       int J = ijk[1];
02760       int K = ijk[2];
02761       int NumberOfNodes = I*J*K;
02762       if ( !NumberOfNodes ) continue;
02763       haveStructuredParts++;
02764 
02765       // read coordinates
02766       int nodeShift = points.empty() ? 0 : points.rbegin()->first;
02767       {
02768         TFltOwner noInterlaceCoords( geoFile.getFlt( NumberOfNodes * SPACE_DIM ));
02769         float* x = noInterlaceCoords;
02770         float* y = x + NumberOfNodes;
02771         float* z = y + NumberOfNodes;
02772         for ( int iNode = 1; iNode <= NumberOfNodes; ++iNode ) {
02773           INoeud inoeud = points.insert( points.end(), make_pair( iNode + nodeShift, _noeud()));
02774           _noeud & node = inoeud->second;
02775           node.number   = inoeud->first;
02776           node.coord.resize( SPACE_DIM );
02777           node.coord[0] = *x++;
02778           node.coord[1] = *y++;
02779           node.coord[2] = *z++;
02780         }
02781       }
02782       // iblanks
02783       if ( iblanked )
02784         geoFile.skip(NumberOfNodes * sizeof(int));
02785 
02786       // let GRID calculate connectivity 
02787       GRID grid;
02788       grid._iArrayLength  = I;
02789       grid._jArrayLength  = J;
02790       grid._kArrayLength  = K;
02791       grid._gridType      = MED_BODY_FITTED;
02792       grid._spaceDimension= SPACE_DIM;
02793       grid._coordinate    = new COORDINATE(3, 1, 0); // to avoid exception in convertInMESH()
02794       if ( J < 2 ) { grid._spaceDimension--; grid._jArrayLength = 0; }
02795       if ( K < 2 ) { grid._spaceDimension--; grid._kArrayLength = 0; }
02796 
02797       const MESH* unstruct = grid.convertInMESH();
02798       const int * conn = unstruct->getConnectivity( MED_NODAL, MED_CELL, MED_ALL_ELEMENTS );
02799       medGeometryElement elemType = grid.getElementType( MED_CELL, 1 );
02800       int  nbElemNodes = elemType % 100;
02801       int      nbElems = grid.getNumberOfElements(MED_CELL, MED_ALL_ELEMENTS);
02802 
02803       partGroupe->mailles.resize( nbElems );
02804       _maille ma( elemType, nbElemNodes );
02805       ma.sommets.resize( nbElemNodes );
02806 
02807       for ( int i = 0, nIndex = 0; i < nbElems; ++i ) {
02808         for ( int n = 0; n < nbElemNodes; ++n ) {
02809           int nodeID = conn[ nIndex++ ];
02810           ma.sommets[n] = points.find( nodeID + nodeShift );
02811         }
02812         //ma.ordre = ++order;
02813         partGroupe->mailles[i] = imed.insert(ma);
02814       }
02815 
02816       _SubPart subPart( partNum, "block" );
02817       subPart.myNbCells    = nbElems;
02818       subPart.myNbNodes    = NumberOfNodes;
02819       subPart.myFirstNode  = points.find( nodeShift + 1 );
02820       subPart.myCellGroupIndex = imed.groupes.size();
02821       imed.addSubPart( subPart );
02822 
02823       unstruct->removeReference();
02824     }
02825     else if ( word == "coordinates" )
02826     {
02827       // ------------------------------
02828       // Unstructured node coordinates
02829       // ------------------------------
02830       int nbNodes = *TIntOwner( geoFile.getInt(1) );
02831 
02832       TIntOwner numbers(0);
02833       if ( haveNodeIds )
02834         numbers.myValues = geoFile.getInt( nbNodes );
02835 
02836       TFltOwner fullInterlaceCoords( geoFile.getFlt( nbNodes * SPACE_DIM ));
02837       float* coord = fullInterlaceCoords;
02838 
02839       cout << "-> loading coordinates of " << nbNodes << " nodes " << endl ;
02840 
02841       INoeud inoeud;
02842       for ( int i=0 ; i < nbNodes ; i++ )
02843       {
02844         if ( haveNodeIds ) {
02845           int nodeID = numbers[ i ];
02846           inoeud = points.insert( make_pair( nodeID, _noeud() )).first;
02847           inoeud->second.number = nodeID;
02848         }
02849         else {
02850           int nodeID = i + 1;
02851           inoeud = points.insert( points.end(), make_pair( nodeID, _noeud()));
02852           inoeud->second.number = nodeID;
02853         }
02854         _noeud & node = inoeud->second;
02855         node.coord.resize( SPACE_DIM );
02856         node.coord[ 0 ] = *coord++;
02857         node.coord[ 1 ] = *coord++;
02858         node.coord[ 2 ] = *coord++;
02859       }
02860 
02861       _SubPartDesc cooDesc = _SubPartDesc::globalCoordDesc();
02862       _SubPart subPart( cooDesc.partNumber(), cooDesc.typeName() );
02863       subPart.myNbNodes    = nbNodes;
02864       subPart.myFirstNode  = points.begin();
02865       imed.addSubPart( subPart );
02866     }
02867     else
02868     {
02869       throw MED_EXCEPTION
02870         ( LOCALIZED( STRING(LOC) << "Unexpected word: " << word <<
02871                      " in " << getDataFileName()));
02872     }
02873   } // while ( !geoFile.eof() )
02874 
02875   if ( ( haveStructuredParts && haveUnstructuredParts ) || haveStructuredParts > 1 )
02876     imed.mergeNodesAndElements(TOLERANCE);
02877 
02878   END_OF_MED(LOC);
02879 }
02880 
02881 //================================================================================
02885 //================================================================================
02886 
02887 int ENSIGHT_MESH_RDONLY_DRIVER::countParts(const string& geomFileName,
02888                                            const bool    isSingleFileMode)
02889 {
02890   const char * LOC ="ENSIGHT_MESH_RDONLY_DRIVER::countParts() : ";
02891 
02892   int nbParts = 0;
02893   if ( isBinaryDataFile( geomFileName ))
02894   {
02895     _BinaryFileReader geoFile(geomFileName);
02896     // check if swapping bytes needed
02897     try {
02898       return countPartsBinary( geoFile, isSingleFileMode );
02899     }
02900     catch (...) {
02901     }
02902     geoFile.swapBytes();
02903     geoFile.rewind();
02904     nbParts = countPartsBinary( geoFile, isSingleFileMode );
02905   }
02906   else
02907   {
02908     _ASCIIFileReader geoFile(geomFileName);
02909 
02910     if ( isSingleFileMode )
02911       while ( !isTimeStepBeginning( geoFile.getLine() ));
02912 
02913     geoFile.getLine(); // description line 1
02914     geoFile.getLine(); // description line 2
02915 
02916     // read "node|element id <off|given|assign|ignore>"
02917     geoFile.getWord(); geoFile.getWord();
02918     string nodeIds = geoFile.getWord();
02919     geoFile.getWord(); geoFile.getWord();
02920     string elemIds = geoFile.getWord();
02921     bool haveNodeIds = ( nodeIds == "given" || nodeIds == "ignore" );
02922     bool haveElemIds = ( elemIds == "given" || elemIds == "ignore" );
02923 
02924     bool isGold = true;
02925     while ( !geoFile.isTimeStepEnd() )
02926     {
02927       string word, restLine, line = geoFile.getLine();
02928       TStrTool::split( line, word, restLine );
02929 
02930       const TEnSightElemType & partType = getEnSightType( word );
02931       if ( partType._medType != MED_ALL_ELEMENTS )
02932       {
02933         //  Unstructured element type encounters
02934         // --------------------------------------
02935         int  nbElemNodes = partType._medType % 100;
02936         int      nbElems = geoFile.getInt(); // ne
02937 
02938         // element ids
02939         if ( haveElemIds && isGold )
02940           geoFile.skip( nbElems, /*nbPerLine =*/ 1, INT_WIDTH_GOLD ); // id_e
02941 
02942         // skip connectivity
02943         int nbNodes = nbElems * nbElemNodes;
02944         if ( partType._name == "nsided" ) // polygons
02945         {
02946           for ( int i = 0; i < nbElems; ++i )
02947             nbNodes += geoFile.getInt();
02948           geoFile.skip( nbNodes * INT_WIDTH_GOLD, /*nbLines = */nbElems );
02949         }
02950         else if ( partType._name == "nfaced" ) // polyhedrons
02951         {
02952           int nbFaces = 0;
02953           for ( int i = 0; i < nbElems; ++i )
02954             nbFaces += geoFile.getInt();
02955           for ( int f = 0; f < nbFaces; ++f )
02956             nbNodes += geoFile.getInt();
02957           geoFile.skip( nbNodes * INT_WIDTH_GOLD, /*nbLines = */nbFaces );
02958         }
02959         else // standard types
02960         {
02961           if ( isGold )
02962             geoFile.skip( nbNodes, nbElemNodes, INT_WIDTH_GOLD );
02963           else if ( haveElemIds )
02964             geoFile.skip( nbNodes + nbElems, nbElemNodes+1, INT_WIDTH_6 );
02965           else
02966             geoFile.skip( nbNodes, nbElemNodes, INT_WIDTH_6 );
02967         }
02968       }
02969       else if ( word == "coordinates" )
02970       {
02971         isGold = ( nbParts > 0 );
02972         int nbNodes = geoFile.getInt(); // nn
02973 
02974         if ( isGold )
02975         {
02976           if ( haveNodeIds )
02977             geoFile.skip( nbNodes, /*nbPerLine =*/ 1, INT_WIDTH_GOLD ); // node ids
02978           geoFile.skip( nbNodes * SPACE_DIM, /*nbPerLine =*/ 1, FLT_WIDTH ); //  coordinates
02979         }
02980         else {
02981           int coordWidth = 3 * FLT_WIDTH;
02982           if ( haveNodeIds )
02983             coordWidth += INT_WIDTH_6;
02984           geoFile.skip(nbNodes, /*nbPerLine =*/ 1, coordWidth);
02985         }
02986       }
02987       else if ( word == "part" )
02988       {
02989         nbParts++;
02990         if ( isGold )
02991           geoFile.skip( 1, /*nbPerLine =*/ 1, INT_WIDTH_GOLD ); //part number
02992         else
02993           geoFile.getWord(); // part number
02994         geoFile.toNextLine();
02995         geoFile.getLine(); // description line
02996       }
02997       else if ( word == "block" )
02998       {
02999         // Structured type
03000         // ------------------
03001         bool rectilinear = ( restLine.find( "rectilinear" ) != restLine.npos );
03002         bool uniform     = ( restLine.find( "uniform" )     != restLine.npos );
03003         bool curvilinear = ( !rectilinear && !uniform );
03004         bool iblanked    = ( restLine.find( "iblanked" )    != restLine.npos );
03005         bool with_ghost  = ( restLine.find( "with_ghost" )  != restLine.npos );
03006         bool range       = ( restLine.find( "range" )       != restLine.npos );
03007 
03008         // dimension
03009         int I = geoFile.getInt();
03010         int J = geoFile.getInt();
03011         int K = geoFile.getInt();
03012         int nbNodes = I*J*K;
03013         if ( !nbNodes ) continue;
03014 
03015         // range
03016         if ( range ) {
03017           vector<int> ijkRange; // imin imax jmin jmax kmin kmax
03018           ijkRange.reserve(6);
03019           while ( ijkRange.size() < 6 )
03020             ijkRange.push_back( geoFile.getInt() );
03021           I = ijkRange[1]-ijkRange[0]+1;
03022           J = ijkRange[3]-ijkRange[2]+1;
03023           K = ijkRange[5]-ijkRange[4]+1;
03024           nbNodes = I*J*K;
03025         }
03026         int nbElems = (I-1)*(J-1)*(K-1);
03027 
03028         if ( curvilinear ) // read coordinates for all nodes
03029         {
03030           if ( isGold )
03031             geoFile.skip( nbNodes, /*nbPerLine =*/ 1, FLT_WIDTH );
03032           else
03033             geoFile.skip( nbNodes * SPACE_DIM, /*nbPerLine =*/ 6, FLT_WIDTH );
03034         }
03035         else if ( rectilinear ) // read delta vectors with non-regular spacing 
03036         {
03037           geoFile.skip( I + J + K, /*nbPerLine =*/ 1, FLT_WIDTH );
03038         }
03039         else // uniform: read grid origine and delta vectors for regular spacing grid
03040         {
03041           geoFile.skip( 6, /*nbPerLine =*/ 1, FLT_WIDTH );
03042         }
03043 
03044         // iblanks
03045         if ( iblanked ) {
03046           if ( isGold )
03047             geoFile.skip( nbNodes, /*nbPerLine =*/ 1, INT_WIDTH_GOLD);
03048           else
03049             geoFile.skip( nbNodes, /*nbPerLine =*/ 10, INT_WIDTH_6 );
03050         }
03051         // ghosts
03052         if ( with_ghost ) {
03053           geoFile.getWord(); // "ghost_flags"
03054           geoFile.skip( nbElems, /*nbPerLine =*/ 1, INT_WIDTH_GOLD);
03055         }
03056         // node ids
03057         if ( haveNodeIds && geoFile.lookAt( "node_ids" )) {
03058           geoFile.getWord(); // "node_ids"
03059           geoFile.skip( nbNodes, /*nbPerLine =*/ 1, INT_WIDTH_GOLD);
03060         }
03061         // element ids
03062         if ( haveElemIds && geoFile.lookAt( "element_ids" ) ) {
03063           geoFile.getWord(); // "element_ids"
03064           geoFile.skip( nbElems, /*nbPerLine =*/ 1, INT_WIDTH_GOLD);
03065         }
03066       }
03067       else if ( word == "extents" ) {
03068         geoFile.getLine(); geoFile.getLine(); geoFile.getLine();// 3 x 2E12.5
03069       }
03070       else
03071       {
03072         throw MED_EXCEPTION ( LOCALIZED( STRING(LOC) << "Unexpected word: " << word));
03073       }
03074     }
03075   }
03076   return nbParts;
03077 }
03078 
03079 //================================================================================
03083 //================================================================================
03084 
03085 int ENSIGHT_MESH_RDONLY_DRIVER::countPartsBinary(_BinaryFileReader& geoFile,
03086                                                  const bool         isSingleFileMode)
03087 {
03088   const char * LOC ="ENSIGHT_MESH_RDONLY_DRIVER::countPartsBinary() : ";
03089 
03090   if ( geoFile.getPosition() == 0 ) {
03091     TStrOwner format( geoFile.getLine() ); // "C|Fortran Binary"
03092     if ( !contains( "C Binary", format )) {
03093       if ( contains( "Fortran Binary", format ))
03094         throw(MEDEXCEPTION(STRING(LOC) << "Fortran Binary format not supported"));
03095       else
03096         throw(MEDEXCEPTION(STRING(LOC) << "unexpected line: \n" << format.myValues));
03097     }
03098   }
03099 
03100   if ( isSingleFileMode ) {
03101     while (1) {
03102       TStrOwner line( geoFile.getLine() );
03103       if ( isTimeStepBeginning( line.myValues ))
03104         break;
03105     }
03106   }
03107 
03108   // 2 description lines
03109   // ----------------------
03110   geoFile.skip( 2 * MAX_LINE_LENGTH );
03111 
03112   // read "node|element id <off|given|assign|ignore>"
03113   bool haveNodeIds, haveElemIds;
03114   {
03115     TStrOwner nodeIds( geoFile.getLine() );
03116     TStrOwner elemIds( geoFile.getLine() );
03117 
03118     haveNodeIds = ( contains( "given", nodeIds ) || contains( "ignore", nodeIds ) );
03119     haveElemIds = ( contains( "given", elemIds ) || contains( "ignore", elemIds ) );
03120   }
03121 
03122   int nbParts = 0; // the result
03123   bool isGold = true;
03124 
03125   while ( !geoFile.eof() )
03126   {
03127     TStrOwner line( geoFile.getLine() );
03128     if ( isSingleFileMode && isTimeStepEnd( line.myValues ))
03129       break;
03130     string word, restLine;
03131     TStrTool::split( line.myValues, word, restLine );
03132 
03133     const TEnSightElemType & partType = getEnSightType( word );
03134     if ( partType._medType != MED_ALL_ELEMENTS )
03135     {
03136       //  Unstructured element type encounters
03137       // --------------------------------------
03138       int      nbElems = *TIntOwner( geoFile.getInt(1) ); // ne
03139       int  nbElemNodes = partType._medType % 100;
03140 
03141       // read element ids
03142       if ( haveElemIds )
03143         geoFile.skip( nbElems * sizeof(int) ); // id_e
03144 
03145       int nbNodes = nbElems * nbElemNodes;
03146       if ( partType._name == "nsided" ) // polygons
03147       {
03148         TIntOwner nbNodesInFace( geoFile.getInt( nbElems ));
03149         nbNodes = std::accumulate( nbNodesInFace.myValues, nbNodesInFace.myValues + nbElems, 0 );
03150       }
03151       else if ( partType._name == "nfaced" ) // polyhedrons
03152       {
03153         TIntOwner nbElemFaces( geoFile.getInt( nbElems ));
03154         int nbFaces = accumulate( nbElemFaces.myValues, nbElemFaces.myValues + nbElems, 0 );
03155         TIntOwner nbNodesInFace( geoFile.getInt( nbFaces ));
03156         nbNodes = std::accumulate( nbNodesInFace.myValues, nbNodesInFace.myValues + nbFaces, 0 );
03157       }
03158       geoFile.skip( nbNodes * sizeof(int) );
03159     }
03160     else if ( word == "coordinates" )
03161     {
03162       if ( nbParts == 0 )
03163         isGold = false;
03164       int nbNodes = *TIntOwner( geoFile.getInt(1) ); // nn
03165 
03166       // node ids
03167       if ( haveNodeIds )
03168         geoFile.skip( nbNodes * sizeof(int) ); // id_n
03169 
03170       // coordinates
03171       geoFile.skip( nbNodes * SPACE_DIM * sizeof(int) );
03172     }
03173     else if ( word == "part" )
03174     {
03175       ++nbParts;
03176       if ( isGold ) geoFile.skip(sizeof(int)); // part #
03177 
03178       geoFile.skip(MAX_LINE_LENGTH); // description line
03179     }
03180     else if ( word == "block" )
03181     {
03182       // Structured type
03183       // ------------------
03184       bool rectilinear = ( restLine.find( "rectilinear" ) != restLine.npos );
03185       bool uniform     = ( restLine.find( "uniform" )     != restLine.npos );
03186       bool curvilinear = ( !rectilinear && !uniform );
03187       bool iblanked    = ( restLine.find( "iblanked" )    != restLine.npos );
03188       bool with_ghost  = ( restLine.find( "with_ghost" )  != restLine.npos );
03189       bool range       = ( restLine.find( "range" )       != restLine.npos );
03190 
03191       // dimension
03192       TIntOwner ijk( geoFile.getInt(3) );
03193       int I = ijk[0];
03194       int J = ijk[1];
03195       int K = ijk[2];
03196       int NumberOfNodes = I*J*K;
03197       if ( !NumberOfNodes ) {
03198         if ( I != 0 && J != 0 && K != 0 )
03199           throw MEDEXCEPTION( "Need to swap bytes" );
03200         continue;
03201       }
03202 
03203       // range
03204       if ( range ) {
03205         TIntOwner ijkRange( geoFile.getInt( 6 ));// imin imax jmin jmax kmin kmax
03206         I = ijkRange[1]-ijkRange[0]+1;
03207         J = ijkRange[3]-ijkRange[2]+1;
03208         K = ijkRange[5]-ijkRange[4]+1;
03209         NumberOfNodes = I*J*K;
03210       }
03211       int nbElems = (I-1)*(J-1)*(K-1);
03212 
03213       if ( curvilinear ) // read coordinates for all nodes
03214       {
03215         geoFile.skip( NumberOfNodes * SPACE_DIM * sizeof(float) );
03216       }
03217       else if ( rectilinear ) // read delta vectors with non-regular spacing 
03218       {
03219         geoFile.skip( (I+J+K) * sizeof(float) );
03220       }
03221       else // uniform: read grid origine and delta vectors for regular spacing grid
03222       {
03223         geoFile.skip( 6 * sizeof(float) );
03224       }
03225 
03226       // iblanks
03227       if ( iblanked )
03228         geoFile.skip( NumberOfNodes * sizeof(int) );
03229       // ghosts
03230       if ( with_ghost ) {
03231         geoFile.skip( MAX_LINE_LENGTH ); // "ghost_flags"
03232         geoFile.skip( nbElems * sizeof(int) );
03233       }
03234       // node ids
03235       if ( haveNodeIds && isGold && !geoFile.eof()  ) {
03236         TStrOwner nextLine( geoFile.getLine() ); // "node_ids"
03237         if ( contains( "node_ids", nextLine ) )
03238           geoFile.skip( NumberOfNodes * sizeof(int) );
03239         else
03240           geoFile.skip( -MAX_LINE_LENGTH );
03241       }
03242       // element ids
03243       if ( haveElemIds && isGold && !geoFile.eof() ) {
03244         TStrOwner nextLine( geoFile.getLine() ); // "element_ids"
03245         if ( contains( "element_ids", nextLine ) )
03246           geoFile.skip( nbElems * sizeof(int) );
03247         else
03248           geoFile.skip( -MAX_LINE_LENGTH );
03249       }
03250     }
03251     else if ( word == "extents" )
03252     {
03253       geoFile.skip( 6 * sizeof(float) );
03254     }
03255     else
03256     {
03257       throw MED_EXCEPTION ( LOCALIZED( STRING(LOC) << "Unexpected word: " << word ));
03258     }
03259   } // while ( !geoFile.eof() )
03260 
03261   return nbParts;
03262 }
03263 
03264 GROUP* ENSIGHT_MESH_RDONLY_DRIVER::makeGroup( _groupe&     grp,
03265                                               _InterMed & imed)
03266 {
03267   //const char* LOC = "ENSIGHT_MESH_RDONLY_DRIVER::makeGroup(): error";
03268 
03269   // prevent creation of other groups but only this one
03270   for (size_t i=0; i < imed.groupes.size(); ++i)
03271     imed.groupes[i].nom.clear();
03272 
03273   // let _intermediateMED create a GROUP from grp
03274   grp.medGroup = 0; // the new GROUP should appear in grp.medGroup
03275   grp.nom = "TMP";
03276   vector<GROUP *> tmp;
03277   MESH* mesh = imed._isOwnMedMesh ? (MESH*) 0 : imed._medMesh;
03278   imed.getGroups( tmp, tmp, tmp, tmp, mesh );
03279   if ( !grp.medGroup )
03280     throw MEDEXCEPTION(LOCALIZED("Can't create a GROUP from _groupe"));
03281 
03282   grp.medGroup->setName(""); // to let a caller set a proper name
03283   grp.nom = "";
03284 
03285   // find pre-existing equal _groupe
03286   _groupe * equalGroupe = 0;
03287   for (unsigned int i=0; i < imed.groupes.size() && !equalGroupe; ++i) {
03288     _groupe& g = imed.groupes[i];
03289     if ( &g != &grp && g.medGroup && g.medGroup->deepCompare( *grp.medGroup ))
03290       equalGroupe = & g;
03291   }
03292   if ( equalGroupe ) {
03293     if ( grp.medGroup->getMesh() )
03294       grp.medGroup->getMesh()->addReference(); // let the mesh survive
03295     grp.medGroup->removeReference();
03296     grp.medGroup = equalGroupe->medGroup;
03297   }
03298   else { // new unique group
03299 
03300     if ( grp.medGroup->isOnAllElements() ) // on all elems
03301       grp.medGroup->setName( string("SupportOnAll_")+entNames[ grp.medGroup->getEntity() ] );
03302 
03303     // add a new group to mesh
03304     if ( !imed._isOwnMedMesh ) {
03305       vector<GROUP*> * groups = 0;
03306       switch ( grp.medGroup->getEntity() ) {
03307       case MED_CELL: groups = & imed._medMesh->_groupCell; break;
03308       case MED_FACE: groups = & imed._medMesh->_groupFace; break;
03309       case MED_EDGE: groups = & imed._medMesh->_groupEdge; break;
03310       case MED_NODE: groups = & imed._medMesh->_groupNode; break;
03311       default:;
03312       }
03313       if ( groups ) {
03314         groups->resize( groups->size() + 1 );
03315         groups->at( groups->size() - 1) = grp.medGroup;
03316       }
03317     }
03318   }
03319   return grp.medGroup;
03320 }