-// Copyright (C) 2005 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
-// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
-//
-// This library is free software; you can redistribute it and/or
-// modify it under the terms of the GNU Lesser General Public
-// License as published by the Free Software Foundation; either
-// version 2.1 of the License.
-//
-// This library is distributed in the hope that it will be useful
-// but WITHOUT ANY WARRANTY; without even the implied warranty of
-// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
-// Lesser General Public License for more details.
+// Copyright (C) 2007-2008 CEA/DEN, EDF R&D, OPEN CASCADE
//
-// You should have received a copy of the GNU Lesser General Public
-// License along with this library; if not, write to the Free Software
-// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+// Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
+// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
//
-// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License.
+//
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+// Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
//
#include "MEDMEM_EnsightMeshDriver.hxx"
#include <sstream>
-#include <strstream>
#include <iomanip>
+#include <numeric>
#include "MEDMEM_define.hxx"
#include "MEDMEM_Family.hxx"
#include "MEDMEM_CellModel.hxx"
#include "MEDMEM_Grid.hxx"
-#include "MEDMEM_Meshing.hxx"
#include "MEDMEM_MedMeshDriver.hxx"
using namespace std;
using namespace MEDMEM;
using namespace MED_EN;
+using namespace MEDMEM_ENSIGHT;
+
+#define TStrTool _ASCIIFileReader
+#define TOLERANCE 1e-15
+
+//#define ELEMENT_ID_GIVEN
+
+namespace {
+
+ // ---------------------------------------------------------------------
+ /*!
+ * \brief The beginning of mesh description used to distinguish files
+ * generated by ENSIGHT_MESH_WRONLY_DRIVER from others
+ */
+ const char* theDescriptionPrefix = "Meshing from MedMemory. ";
+
+ // ---------------------------------------------------------------------
+ /*!
+ * \brief Default name of a mesh read from EnSight
+ */
+ const char* theDefaultMeshName = "EnsightMesh";
+
+ // ---------------------------------------------------------------------
+ /*!
+ * \brief Max number of types in EnSight part
+ */
+ const int theMaxNbTypes = 20;
+
+ // ---------------------------------------------------------------------
+ /*!
+ * \brief Make array with elements == index[i+1]-index[i]
+ * \param size - result array size
+ */
+ int* getNumbersByIndex( const int* index, int size, const int* elemNumbers=0)
+ {
+ int* numbers = new int[size];
+ int* n = numbers;
+ if ( elemNumbers ) {
+ const int *elem = elemNumbers-1, *elemEnd = elemNumbers + size;
+ while ( ++elem < elemEnd )
+ *n++ = index[*elem] - index[*elem-1];
+ }
+ else {
+ const int *ind = index, *indEnd = index + size + 1;
+ while ( ++ind < indEnd )
+ *n++ = ind[0] - ind[-1];
+ }
+ return numbers;
+ }
+ // ---------------------------------------------------------------------
+ /*!
+ * \brief Type used to delete numbers returned by getNumbersByIndex()
+ */
+ typedef _ValueOwner<int> TNumbers;
-#define MED_NULL NULL
+} // namespace
-ENSIGHT_MESH_DRIVER::ENSIGHT_MESH_DRIVER(): GENDRIVER(),
- _ptrMesh((MESH *)MED_NULL)
+
+ENSIGHT_MESH_DRIVER::ENSIGHT_MESH_DRIVER(): _CaseFileDriver_User(), _ptrMesh((MESH *)NULL)
{
}
ENSIGHT_MESH_DRIVER::ENSIGHT_MESH_DRIVER(const string & fileName,
- MESH * ptrMesh) :
- GENDRIVER(fileName,MED_EN::MED_RDWR), _ptrMesh(ptrMesh)
+ MESH * ptrMesh)
+ :_CaseFileDriver_User(fileName,MED_EN::RDWR), _ptrMesh(ptrMesh),
+ _meshName(ptrMesh->getName())
{
}
ENSIGHT_MESH_DRIVER::ENSIGHT_MESH_DRIVER(const string & fileName,
- MESH * ptrMesh,
- MED_EN::med_mode_acces accessMode):
- GENDRIVER(fileName,accessMode), _ptrMesh(ptrMesh)
+ MESH * ptrMesh,
+ med_mode_acces accessMode)
+ :_CaseFileDriver_User(fileName,accessMode), _ptrMesh(ptrMesh), _meshName(ptrMesh->getName())
{
}
-ENSIGHT_MESH_DRIVER::ENSIGHT_MESH_DRIVER(const ENSIGHT_MESH_DRIVER & driver):
- GENDRIVER(driver),
- _ptrMesh(driver._ptrMesh),
- _meshName(driver._meshName)
+ENSIGHT_MESH_DRIVER::ENSIGHT_MESH_DRIVER(const ENSIGHT_MESH_DRIVER & driver)
+ :_CaseFileDriver_User(driver), _ptrMesh(driver._ptrMesh), _meshName(driver._meshName)
{
}
ENSIGHT_MESH_DRIVER::~ENSIGHT_MESH_DRIVER()
{
- MESSAGE("ENSIGHT_MESH_DRIVER::~ENSIGHT_MESH_DRIVER() has been destroyed");
+ MESSAGE_MED("ENSIGHT_MESH_DRIVER::~ENSIGHT_MESH_DRIVER() has been destroyed");
}
void ENSIGHT_MESH_DRIVER::setMeshName(const string & meshName) { _meshName = meshName; };
string ENSIGHT_MESH_DRIVER::getMeshName() const { return _meshName; };
+void ENSIGHT_MESH_DRIVER::openConst(bool checkDataFile) const
+{
+ const char * LOC ="ENSIGHT_MESH_DRIVER::open() : ";
+ BEGIN_OF_MED(LOC);
+
+ if ( checkDataFile )
+ {
+ if ( getDataFileName().empty() )
+ throw MED_EXCEPTION
+ ( LOCALIZED( STRING(LOC) << "Internal error, geometry file name is empty"));
+
+ if (!canOpenFile( getDataFileName(), getAccessMode() ))
+ throw MED_EXCEPTION
+ ( LOCALIZED( STRING(LOC) << "Can not open Ensight Geometry file " << getDataFileName()
+ << " in access mode " << getAccessMode()));
+ }
+ else
+ {
+ if ( getCaseFileName().empty() )
+ throw MED_EXCEPTION
+ ( LOCALIZED( STRING(LOC) << "Case file name is empty, "
+ "please set a correct fileName before calling open()"));
+
+ if ( !canOpenFile( getCaseFileName(), getAccessMode() ))
+ throw MED_EXCEPTION
+ ( LOCALIZED( STRING(LOC) << "Can not open Ensight Case file " << getCaseFileName()
+ << " in access mode " << getAccessMode()));
+ }
+
+ END_OF_MED(LOC);
+}
+
void ENSIGHT_MESH_DRIVER::open() {
openConst() ;
}
void ENSIGHT_MESH_DRIVER::close() {
- closeConst() ;
}
+// ================================================================================
+// WRONLY
+// ================================================================================
+
ENSIGHT_MESH_WRONLY_DRIVER::ENSIGHT_MESH_WRONLY_DRIVER() : ENSIGHT_MESH_DRIVER()
{
- _ensightFile = new ofstream();
}
-ENSIGHT_MESH_WRONLY_DRIVER::ENSIGHT_MESH_WRONLY_DRIVER(const string & fileName, MESH * ptrMesh) : ENSIGHT_MESH_DRIVER(fileName,ptrMesh)
+ENSIGHT_MESH_WRONLY_DRIVER::ENSIGHT_MESH_WRONLY_DRIVER(const string & fileName,
+ MESH * ptrMesh,
+ bool append)
+ : ENSIGHT_MESH_DRIVER( fileName, ptrMesh, WRONLY ), _append(append)
{
- _ensightFile = new ofstream();
}
-ENSIGHT_MESH_WRONLY_DRIVER::ENSIGHT_MESH_WRONLY_DRIVER(const ENSIGHT_MESH_WRONLY_DRIVER & driver) : ENSIGHT_MESH_DRIVER(driver),_support(driver._support)
+ENSIGHT_MESH_WRONLY_DRIVER::ENSIGHT_MESH_WRONLY_DRIVER(const ENSIGHT_MESH_WRONLY_DRIVER & driver)
+ : ENSIGHT_MESH_DRIVER(driver),_append(driver._append)
{
- _ensightFile = new ofstream();
}
ENSIGHT_MESH_WRONLY_DRIVER::~ENSIGHT_MESH_WRONLY_DRIVER()
{
- delete _ensightFile ;
}
GENDRIVER * ENSIGHT_MESH_WRONLY_DRIVER::copy() const
return new ENSIGHT_MESH_WRONLY_DRIVER(*this) ;
}
-void ENSIGHT_MESH_WRONLY_DRIVER::openConst() const
+void ENSIGHT_MESH_WRONLY_DRIVER::read() throw (MEDEXCEPTION) {
+ throw MEDEXCEPTION("ENSIGHT_MESH_WRONLY_DRIVER::read : Can't read with a WRONLY driver !");
+}
+
+//================================================================================
+/*!
+ * \brief writing
+ */
+//================================================================================
+
+void ENSIGHT_MESH_WRONLY_DRIVER::write() const throw (MEDEXCEPTION)
{
- const char * LOC = "ENSIGHT_MESH_WRONLY_DRIVER::open()" ;
+ const char * LOC = "ENSIGHT_MESH_WRONLY_DRIVER::write() : ";
+ BEGIN_OF_MED(LOC);
- BEGIN_OF(LOC);
+ openConst(false) ; // check if can write to case file
- if ( _fileName == "" )
- throw MED_EXCEPTION ( LOCALIZED( STRING(LOC)
- << "_fileName is |\"\"|, please set a correct fileName before calling open()"
- )
- );
+ // Ensight case organization requires a main file (filename.case) which defines organization
- if (!(*_ensightFile).is_open())
- (*_ensightFile).open(_fileName.c_str()) ;
+ _CaseFileDriver caseFile( getCaseFileName(), this);
+ if ( _append )
+ caseFile.read();
+ caseFile.addMesh( this );
+ caseFile.write(); // all fields of _CaseFileDriver_User are set by this method
- if (!(*_ensightFile))
- throw MED_EXCEPTION ( LOCALIZED( STRING(LOC) << "Could not open main Ensight file "
- << _fileName)
- );
+ openConst(true) ; // check if can write to data file
- END_OF(LOC);
-}
+ cout << "-> creating the Ensight geometry file " << getDataFileName() << endl ;
-void ENSIGHT_MESH_WRONLY_DRIVER::closeConst() const {
-
- const char * LOC = "ENSIGHT_MESH_WRONLY_DRIVER::close() : ";
- BEGIN_OF(LOC);
-
- (*_ensightFile).close();
-
- if (!(*_ensightFile))
- throw MED_EXCEPTION ( LOCALIZED( STRING(LOC) << "Could not close main Ensight file "
- << _fileName)
- );
- END_OF(LOC);
-}
+ // Store mesh description and a special mark in the first two description lines each
+ // of 79 chars length maximum, while MED mesh description is up to 200 chars
+ const char* line1 = theDescriptionPrefix;
+ string line2 = _ptrMesh->getDescription();
+ for ( int i = 0; i < line2.size(); ++i ) { // protect from gabage
+ if ( !isascii( line2[ i ])) {
+ line2.resize( i );
+ break;
+ }
+ }
+ if ( line2.size() >= MAX_LINE_LENGTH )
+ line2.resize( MAX_LINE_LENGTH );
+
+ // EnSight will assign node/element visible numbers it-self
+ const char* line3 = "node id assign";
+#ifdef ELEMENT_ID_GIVEN
+ const char* line4 = "element id given";
+#else
+ const char* line4 = "element id assign";
+#endif
-void ENSIGHT_MESH_WRONLY_DRIVER::read() throw (MEDEXCEPTION) {
- throw MEDEXCEPTION("ENSIGHT_MESH_WRONLY_DRIVER::read : Can't read with a WRONLY driver !");
-}
+ if ( isBinaryEnSightFormatForWriting() )
+ {
+ // ======================================================
+ // Binary
+ // ======================================================
-void ENSIGHT_MESH_WRONLY_DRIVER::write() const throw (MEDEXCEPTION)
-{
+ _BinaryFileWriter ensightGeomFile( getDataFileName() );
- const char * LOC = "ENSIGHT_MESH_WRONLY_DRIVER::write() : ";
- BEGIN_OF(LOC);
+ ensightGeomFile.addString("C Binary");
+ ensightGeomFile.addString(line1);
+ ensightGeomFile.addString(line2);
+ ensightGeomFile.addString(line3);
+ ensightGeomFile.addString(line4);
- // Well we must open ensight file first, because there are
- // no other driver than MESH for ENSIGHT that do it !
- openConst() ;
+ // function to write a support as a part
+ typedef void (ENSIGHT_MESH_WRONLY_DRIVER::* TWritePart) (_BinaryFileWriter&, const SUPPORT*) const;
+ TWritePart writePart;
+ if ( isGoldFormat() )
+ {
+ // GOLD
+ writePart = & ENSIGHT_MESH_WRONLY_DRIVER::writePartGoldBinary;
+ }
+ else
+ {
+ // ENSIGHT 6. Write addionally global nodes
+ writePart = & ENSIGHT_MESH_WRONLY_DRIVER::writePart6Binary;
+
+ // All point are in 3D, so if we are in 1D or 2D, we complete by zero !
+ int SpaceDimension = _ptrMesh->getSpaceDimension() ;
+ int NumberOfNodes = _ptrMesh->getNumberOfNodes() ;
+ ensightGeomFile.addString("coordinates");
+ ensightGeomFile.addInt( NumberOfNodes );
+ const double *coordinate = _ptrMesh->getCoordinates(MED_FULL_INTERLACE) ;
+ if ( SpaceDimension == 3 ) {
+ ensightGeomFile.addReal(coordinate, NumberOfNodes * SpaceDimension );
+ }
+ else {
+ typedef _ValueIterator< double > TComponentIt;
+ vector< TComponentIt > coordCompIt( 3 );
+ for (int j=0; j<3; j++, coordinate++)
+ if ( j < SpaceDimension )
+ coordCompIt[ j ] = TComponentIt( coordinate, SpaceDimension );
+ else
+ coordCompIt[ j ] = TComponentIt(); // to iterate on zeros
+ ensightGeomFile.addReal( coordCompIt, NumberOfNodes, MED_FULL_INTERLACE );
+ }
+ }
- // Ensight case organization requires a main file (filename.case) which defines organization
- // We also need a geom file (filemane.geo) to store the meshs
- // We also need a data file (filemane.data) to store the fields
-
- string MeshName = _ptrMesh->getName() ;
- // In fact, we must take care of all supports
- // We restrict Field on all nodes or cells
-
- cout << "-> creating the Ensight case file " << _fileName << endl ;
- (*_ensightFile) << "FORMAT" << endl ;
- (*_ensightFile) << "type: ensight" << endl ;
- (*_ensightFile) << endl ;
- (*_ensightFile) << "GEOMETRY" << endl ;
- int len = _fileName.size() ;
- string prefix = _fileName.substr(0,len-5); // extraction de .case
- string ensight_geomf = prefix + ".geom" ;
- string basen = getBaseName((char*)ensight_geomf.c_str());
- (*_ensightFile) << "# Mesh detected with name = " << MeshName << endl ;
- (*_ensightFile) << "model: " << basen << endl ;
- (*_ensightFile) << endl ;
-
- ofstream ensightGeomFile(ensight_geomf.c_str(),ios::out);
- cout << "-> creating the Ensight geometry file " << ensight_geomf << endl ;
-
- // ------------ create the Ensight file associated to this meshe
- ensightGeomFile << "Ensight Geometry File : " << endl
- << "Meshing from MedMemory" << endl ;
-// ensightGeomFile << "node id given " << endl ;
-// ensightGeomFile << "element id given " << endl ;
- ensightGeomFile << "node id assign " << endl ;
- ensightGeomFile << "element id assign " << endl ;
- ensightGeomFile << "coordinates" << endl ;
- // put points (all point are in 3D, so if we are in 1D or 2D, we complete by zero !
+ // We put connectivity
+
+ if ( isToWriteEntity( MED_CELL, _ptrMesh ))
+ {
+ SUPPORT allCells(_ptrMesh, getMeshName(), MED_CELL );
+ (this->*writePart)( ensightGeomFile, &allCells );
+ }
+ // And meshdim-1 connectivity
+ if ( isToWriteEntity( MED_FACE, _ptrMesh ))
+ {
+ SUPPORT allFaces(_ptrMesh, string("SupportOnAll_")+entNames[MED_FACE], MED_FACE );
+ (this->*writePart)( ensightGeomFile, & allFaces);
+ }
+ else if ( isToWriteEntity(MED_EDGE, _ptrMesh))
+ {
+ SUPPORT allEdges(_ptrMesh, string("SupportOnAll_")+entNames[MED_EDGE], MED_EDGE );
+ (this->*writePart)( ensightGeomFile, & allEdges);
+ }
+
+ // Write all groups as parts
+
+ for ( int ent = MED_CELL; ent < MED_ALL_ENTITIES; ++ent )
+ {
+ medEntityMesh entity = (medEntityMesh) ent;
+ int nbGroups = _ptrMesh->getNumberOfGroups(entity);
+ for ( int i=1; i<=nbGroups; i++)
+ {
+ const GROUP* group = _ptrMesh->getGroup( entity, i );
+ (this->*writePart)( ensightGeomFile, group );
+ }
+ }
+
+ }
+ else
+ {
+ // ======================================================
+ // ASCII
+ // ======================================================
+ ofstream ensightGeomFile( getDataFileName().c_str(), ios::out);
+ ensightGeomFile.setf(ios::scientific);
+ ensightGeomFile.precision(5);
+
+ ensightGeomFile << line1 << endl
+ << line2 << endl
+ << line3 << endl
+ << line4 << endl;
+
+ // function to write a support as a part
+ typedef void (ENSIGHT_MESH_WRONLY_DRIVER::* TWritePart) (ofstream&, const SUPPORT*) const;
+ TWritePart writePart;
+ if ( isGoldFormat() )
+ {
+ // GOLD
+ writePart = & ENSIGHT_MESH_WRONLY_DRIVER::writePartGoldASCII;
+ }
+ else
+ {
+ // ENSIGHT 6. Write addionally global nodes
+ writePart = & ENSIGHT_MESH_WRONLY_DRIVER::writePart6ASCII;
+
+ // Put points (all point are in 3D, so if we are in 1D or 2D, we complete by zero !
+ int SpaceDimension = _ptrMesh->getSpaceDimension() ;
+ int NumberOfNodes = _ptrMesh->getNumberOfNodes() ;
+ string zeros;
+ if (SpaceDimension==2) zeros = " 0.00000e+00";
+ if (SpaceDimension==1) zeros = " 0.00000e+00 0.00000e+00";
+ ensightGeomFile << "coordinates" << endl
+ << setw(8) << NumberOfNodes << endl ;
+ const double *coordinate = _ptrMesh->getCoordinates(MED_FULL_INTERLACE) ;
+ for (int i=0; i<NumberOfNodes; i++)
+ {
+ //ensightGeomFile << setw(8) << i+1 ; // node id
+ for (int j=0; j<SpaceDimension; j++, coordinate++)
+ ensightGeomFile << setw(12) << *coordinate;
+ ensightGeomFile << zeros << endl ;
+ }
+ }
+
+ // We put connectivity
+
+ if ( isToWriteEntity( MED_CELL, _ptrMesh ))
+ {
+ SUPPORT allCells(_ptrMesh, getMeshName(), MED_CELL );
+ (this->*writePart)( ensightGeomFile, &allCells );
+ }
+ // And meshdim-1 connectivity
+ if ( isToWriteEntity( MED_FACE, _ptrMesh ))
+ {
+ SUPPORT allFaces(_ptrMesh, string("SupportOnAll_")+entNames[MED_FACE], MED_FACE );
+ (this->*writePart)( ensightGeomFile, & allFaces);
+ }
+ else if ( isToWriteEntity(MED_EDGE, _ptrMesh))
+ {
+ SUPPORT allEdges(_ptrMesh, string("SupportOnAll_")+entNames[MED_EDGE], MED_EDGE );
+ (this->*writePart)( ensightGeomFile, & allEdges);
+ }
+
+ // Write all groups as parts
+
+ for ( int ent = MED_CELL; ent < MED_ALL_ENTITIES; ++ent )
+ {
+ medEntityMesh entity = (medEntityMesh) ent;
+ int nbGroups = _ptrMesh->getNumberOfGroups(entity);
+ for ( int i=1; i<=nbGroups; i++)
+ {
+ const GROUP* group = _ptrMesh->getGroup( entity, i );
+ (this->*writePart)( ensightGeomFile, group );
+ }
+ }
+
+ ensightGeomFile.close();
+
+ } // end ASCII format
+
+} // ENSIGHT_MESH_WRONLY_DRIVER::write()
+
+//================================================================================
+/*!
+ * \brief Write support as an EnSight Gold part
+ */
+//================================================================================
+
+void ENSIGHT_MESH_WRONLY_DRIVER::writePartGoldBinary(_BinaryFileWriter& ensightGeomFile,
+ const SUPPORT* support) const
+{
+ // part number
+ int partNum = getPartNumber( support );
+ if ( !partNum )
+ throw MED_EXCEPTION ( LOCALIZED( STRING("Internal error: invalid part number")));
+ ensightGeomFile.addString( "part" );
+ ensightGeomFile.addInt( partNum );
+
+ // group/mesh name
+ ensightGeomFile.addString( support->getName() );
+
+ // get geom types
+ medEntityMesh entity = support->getEntity();
+ int nbTypes = support->getNumberOfTypes();
+ const medGeometryElement* geoType = support->getTypes();
+
+ const int * connectivity = 0;
+ const int * elemConnectivity = 0;
+ const int * index = 0;
+ int j;
+
+ // COORDINATES Gold binary
+ // ===================================================================================
+ // In Ensight, coordinates of nodes of support elements are in MED_NO_INTERLACE mode.
+ // We are to write only nodes belonging to elements of the support and
+ // nodal connectivity should refer to these nodes.
+ map<int, int> med2ensIds;
+ map<int, int>::iterator medEnsIt;
int SpaceDimension = _ptrMesh->getSpaceDimension() ;
int NumberOfNodes = _ptrMesh->getNumberOfNodes() ;
- ensightGeomFile << NumberOfNodes << endl ;
- const double *coordinate = _ptrMesh->getCoordinates(MED_FULL_INTERLACE) ;
- ensightGeomFile.setf(ios::scientific);
- ensightGeomFile.precision(5);
- for (int i=0;i<NumberOfNodes;i++) {
-// ensightGeomFile << setw(8) << i+1 ;
- for (int j=0;j<SpaceDimension;j++)
- ensightGeomFile << setw(12) << coordinate[i*SpaceDimension+j] ;
- if (SpaceDimension==1)
- ensightGeomFile << "0 0" ;
- if (SpaceDimension==2)
- ensightGeomFile << "0" ;
- ensightGeomFile << endl ;
- }
-
- // we put connectivity
- // how many cells and how many value in connectivity :
- int cells_types_count = _ptrMesh->getNumberOfTypes(MED_CELL) ;
- const CELLMODEL * cells_type =_ptrMesh->getCellsTypes(MED_CELL) ;
- ensightGeomFile << "part 1 " << endl ;
- ensightGeomFile << "elements are following " << endl ;
-
- // we put connectivity
- for (int i=0;i<cells_types_count;i++) {
- int numberOfCell = _ptrMesh->getNumberOfElements(MED_CELL,cells_type[i].getType()) ;
- int *filter = (int*) NULL ;
- switch (cells_type[i].getType())
- {
- case MED_POINT1 : {
- ensightGeomFile << "point" << endl ;
- ensightGeomFile << setw(8) << numberOfCell << endl ;
- filter = new int[1] ;
- filter[0] = 0 ;
- break ;
- }
- case MED_SEG2 : {
- ensightGeomFile << "bar2" << endl ;
- ensightGeomFile << setw(8) << numberOfCell << endl ;
- filter = new int[2] ;
- filter[0] = 0 ;
- filter[1] = 1 ;
- break ;
- }
- case MED_SEG3 : {
- filter = new int[3] ;
- filter[0] = 0 ;
- filter[1] = 2 ;
- filter[2] = 1 ;
- break ;
- }
- case MED_TRIA3 : {
- ensightGeomFile << "tria3" << endl ;
- ensightGeomFile << setw(8) << numberOfCell << endl ;
- filter = new int[3] ;
- filter[0] = 0 ;
- filter[1] = 2 ;
- filter[2] = 1 ;
- break ;
- }
- case MED_QUAD4 : {
- ensightGeomFile << "quad4" << endl ;
- ensightGeomFile << setw(8) << numberOfCell << endl ;
- filter = new int[4] ;
- filter[0] = 0 ;
- filter[1] = 3 ;
- filter[2] = 2 ;
- filter[3] = 1 ;
- break ;
- }
- case MED_TRIA6 : {
- ensightGeomFile << "tria6" << endl ;
- ensightGeomFile << setw(8) << numberOfCell << endl ;
- filter = new int[6] ;
- filter[0] = 0 ;
- filter[1] = 2 ;
- filter[2] = 1 ;
- filter[3] = 5 ;
- filter[4] = 4 ;
- filter[5] = 3 ;
- break ;
- }
- case MED_QUAD8 : {
- ensightGeomFile << "quad8" << endl ;
- ensightGeomFile << setw(8) << numberOfCell << endl ;
- filter = new int[8] ;
- filter[0] = 0 ;
- filter[1] = 3 ;
- filter[2] = 2 ;
- filter[3] = 1 ;
- filter[4] = 7 ;
- filter[5] = 6 ;
- filter[6] = 5 ;
- filter[7] = 4 ;
- break ;
- }
- case MED_TETRA4 : {
- ensightGeomFile << "tetra4" << endl ;
- ensightGeomFile << setw(8) << numberOfCell << endl ;
- filter = new int[4] ;
- filter[0] = 0 ;
- filter[1] = 1 ;
- filter[2] = 3 ; // 3td element in med are 4th in vtk (array begin at 0 !)
- filter[3] = 2 ; // 4th element in med are 3rd in vtk (array begin at 0 !)
- break ;
- }
- case MED_PYRA5 : {
- ensightGeomFile << "pyramid5" << endl ;
- ensightGeomFile << setw(8) << numberOfCell << endl ;
- filter = new int[5] ;
- filter[0] = 0 ;
- filter[1] = 3 ; // 2nd element in med are 4th in vtk (array begin at 0 !)
- filter[2] = 2 ;
- filter[3] = 1 ; // 4th element in med are 2nd in vtk (array begin at 0 !)
- filter[4] = 4 ;
- break ;
- }
- case MED_PENTA6 : {
- ensightGeomFile << "penta6" << endl ;
- ensightGeomFile << setw(8) << numberOfCell << endl ;
- filter = new int[6] ;
- filter[0] = 0 ;
- filter[1] = 2 ;
- filter[2] = 1 ;
- filter[3] = 3 ;
- filter[4] = 4 ;
- filter[5] = 5 ;
- break ;
- }
- case MED_HEXA8 : {
- ensightGeomFile << "hexa8" << endl ;
- ensightGeomFile << setw(8) << numberOfCell << endl ;
- filter = new int[8] ;
- filter[0] = 0 ;
- filter[1] = 3 ;
- filter[2] = 2 ;
- filter[3] = 1 ;
- filter[4] = 4 ;
- filter[5] = 7 ;
- filter[6] = 6 ;
- filter[7] = 5 ;
- break ;
- }
- case MED_TETRA10 : {
- ensightGeomFile << "tetra10" << endl ;
- ensightGeomFile << setw(8) << numberOfCell << endl ;
- filter = new int[10] ;
- filter[0] = 0 ;
- filter[1] = 2 ;
- filter[2] = 1 ;
- filter[3] = 3 ;
- filter[4] = 6 ;
- filter[5] = 5 ;
- filter[6] = 4 ;
- filter[7] = 7 ;
- filter[8] = 9 ;
- filter[9] = 8 ;
- break ;
- }
- case MED_PYRA13 : {
- ensightGeomFile << "pyramid13" << endl ;
- ensightGeomFile << setw(8) << numberOfCell << endl ;
- break ;
- }
- case MED_PENTA15 : {
- ensightGeomFile << "penta15" << endl ;
- ensightGeomFile << setw(8) << numberOfCell << endl ;
- break ;
- }
- case MED_HEXA20 : {
- ensightGeomFile << "hexa20" << endl ;
- ensightGeomFile << setw(8) << numberOfCell << endl ;
- break ;
- }
- default : {
- break ;
- }
- }
- if (filter==NULL)
- throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<": MED element type not supported yet : " << cells_type[i].getName() ) ) ;
- int nodes_cell = cells_type[i].getNumberOfNodes();
- const int * connectivityArray = _ptrMesh->getConnectivity(MED_FULL_INTERLACE,MED_NODAL,MED_CELL,cells_type[i].getType());
- for (int j=0;j<numberOfCell;j++) {
- for (int k=0;k<nodes_cell;k++)
- ensightGeomFile << setw(8) << connectivityArray[j*nodes_cell+filter[k]] ;
- ensightGeomFile << endl ;
- }
- if (filter != NULL)
- delete[] filter ;
- }
-
- for(int i=0;i<(int)_support.size();i++){
- // we put connectivity
- // how many cells and how many value in connectivity :
- int nbTypes = _support[i]->getNumberOfTypes() ;
- const medGeometryElement * geo_type = _support[i]->getTypes() ;
-
- ensightGeomFile << "part " << i+2 << endl;
- ensightGeomFile << "connectivities description" << endl;
-
- int nodes_cell;
- // we put connectivity
- for (int i=0;i<nbTypes;i++) {
-
- int numberOfCell = _support[i]->getNumberOfElements(geo_type[i]) ;
- int *filter = (int*) NULL ; // index in ensight connectivity
- switch (geo_type[i]){
- case MED_POINT1 : {
- nodes_cell=1;
- ensightGeomFile << "point" << endl ;
- ensightGeomFile << setw(8) << numberOfCell << endl ;
- filter = new int[1] ;
- filter[0] = 0 ;
- break ;
- }
- case MED_SEG2 : {
- nodes_cell=2;
- ensightGeomFile << "bar2" << endl ;
- ensightGeomFile << setw(8) << numberOfCell << endl ;
- filter = new int[2] ;
- filter[0] = 0 ;
- filter[1] = 1 ;
- break ;
- }
- case MED_SEG3 : {
- nodes_cell=3;
- ensightGeomFile << "bar3" << endl ;
- ensightGeomFile << setw(8) << numberOfCell << endl ;
- filter = new int[3] ;
- filter[0] = 0 ;
- filter[1] = 2 ;
- filter[2] = 1 ;
- break ;
- }
- case MED_TRIA3 : {
- nodes_cell=3;
- ensightGeomFile << "tria3" << endl ;
- ensightGeomFile << setw(8) << numberOfCell << endl ;
- filter = new int[3] ;
- filter[0] = 0 ;
- filter[1] = 2 ;
- filter[2] = 1 ;
- break ;
- }
- case MED_QUAD4 : {
- nodes_cell=4;
- ensightGeomFile << "quad4" << endl ;
- ensightGeomFile << setw(8) << numberOfCell << endl ;
- filter = new int[4] ;
- filter[0] = 0 ;
- filter[1] = 3 ;
- filter[2] = 2 ;
- filter[3] = 1 ;
- break ;
- }
- case MED_TRIA6 : {
- nodes_cell=6;
- ensightGeomFile << "tria6" << endl ;
- ensightGeomFile << setw(8) << numberOfCell << endl ;
- filter = new int[6] ;
- filter[0] = 0 ;
- filter[1] = 2 ;
- filter[2] = 1 ;
- filter[3] = 5 ;
- filter[4] = 4 ;
- filter[5] = 3 ;
- break ;
- }
- case MED_QUAD8 : {
- nodes_cell=8;
- ensightGeomFile << "quad8" << endl ;
- ensightGeomFile << setw(8) << numberOfCell << endl ;
- filter = new int[8] ;
- filter[0] = 0 ;
- filter[1] = 3 ;
- filter[2] = 2 ;
- filter[3] = 1 ;
- filter[4] = 7 ;
- filter[5] = 6 ;
- filter[6] = 5 ;
- filter[7] = 4 ;
- break ;
- }
- case MED_TETRA4 : {
- nodes_cell=4;
- ensightGeomFile << "tetra4" << endl ;
- ensightGeomFile << setw(8) << numberOfCell << endl ;
- filter = new int[4] ;
- filter[0] = 0 ;
- filter[1] = 1 ;
- filter[2] = 3 ; // 3td element in med are 4th in ensight (array begin at 0 !)
- filter[3] = 2 ; // 4th element in med are 3rd in ensight (array begin at 0 !)
- break ;
- }
- case MED_PYRA5 : {
- nodes_cell=5;
- ensightGeomFile << "pyramid5" << endl ;
- ensightGeomFile << setw(8) << numberOfCell << endl ;
- filter = new int[5] ;
- filter[0] = 0 ;
- filter[1] = 3 ; // 2nd element in med are 4th in ensight (array begin at 0 !)
- filter[2] = 2 ;
- filter[3] = 1 ; // 4th element in med are 2nd in ensight (array begin at 0 !)
- filter[4] = 4 ;
- break ;
- }
- case MED_PENTA6 : {
- nodes_cell=6;
- ensightGeomFile << "penta6" << endl ;
- ensightGeomFile << setw(8) << numberOfCell << endl ;
- filter = new int[6] ;
- filter[0] = 0 ;
- filter[1] = 2 ;
- filter[2] = 1 ;
- filter[3] = 3 ;
- filter[4] = 5 ;
- filter[5] = 4 ;
- break ;
- }
- case MED_HEXA8 : {
- nodes_cell=8;
- ensightGeomFile << "hexa8" << endl ;
- ensightGeomFile << setw(8) << numberOfCell << endl ;
- filter = new int[8] ;
- filter[0] = 0 ;
- filter[1] = 3 ;
- filter[2] = 2 ;
- filter[3] = 1 ;
- filter[4] = 4 ;
- filter[5] = 7 ;
- filter[6] = 6 ;
- filter[7] = 5 ;
- break ;
- }
- case MED_TETRA10 : {
- nodes_cell=10;
- ensightGeomFile << "tetra10" << endl ;
- ensightGeomFile << setw(8) << numberOfCell << endl ;
- filter = new int[10] ;
- filter[0] = 0 ;
- filter[1] = 2 ;
- filter[2] = 1 ;
- filter[3] = 3 ;
- filter[4] = 6 ;
- filter[5] = 5 ;
- filter[6] = 4 ;
- filter[7] = 7 ;
- filter[8] = 9 ;
- filter[9] = 8 ;
- break ;
- }
- case MED_PYRA13 : {
- nodes_cell=13;
- ensightGeomFile << "pyramid13" << endl ;
- ensightGeomFile << setw(8) << numberOfCell << endl ;
- break ;
- }
- case MED_PENTA15 : {
- nodes_cell=15;
- ensightGeomFile << "penta15" << endl ;
- ensightGeomFile << setw(8) << numberOfCell << endl ;
- break ;
- }
- case MED_HEXA20 : {
- nodes_cell=20;
- ensightGeomFile << "hexa20" << endl ;
- ensightGeomFile << setw(8) << numberOfCell << endl ;
- break ;
- }
- default : {
- break ;
- }
- }
- if (filter==NULL)
- throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<": MED element type not supported yet : " << cells_type[i].getName() ) ) ;
-
- const int * connectivityArray = _ptrMesh->getConnectivity(MED_FULL_INTERLACE,MED_NODAL,_support[i]->getEntity(),geo_type[i]);
- const int *type = _support[i]->getNumber(geo_type[i]);
-
- for (int j=0;j<numberOfCell;j++) {
- for (int k=0;k<nodes_cell;k++){
- ensightGeomFile << setw(8) << connectivityArray[(type[j]-1)*nodes_cell+filter[k]];
- }
- ensightGeomFile << endl ;
- }
- if (filter != NULL)
- delete[] filter ;
-
- }
-
- }
-
- ensightGeomFile << endl ;
- return ;
-
- END_OF(LOC);
-}
+ // -------------------------------------------------
+ if ( support->isOnAllElements() )
+ {
+ // nb of nodes
+ ensightGeomFile.addString( "coordinates" );
+ ensightGeomFile.addInt( NumberOfNodes );
+
+ // coordinates
+ const double *coordinate = _ptrMesh->getCoordinates(MED_FULL_INTERLACE);
+ typedef _ValueIterator< double > TComponentIt;
+ vector< TComponentIt > coordCompIt( 1 );
+ for (int j=0; j<SPACE_DIM; j++, coordinate++) { // loop on dimensions
+ if ( j < SpaceDimension )
+ coordCompIt[ 0 ] = TComponentIt( coordinate, SpaceDimension );
+ else
+ coordCompIt[ 0 ] = TComponentIt(); // to iterate on zeros
+ ensightGeomFile.addReal( coordCompIt, NumberOfNodes, MED_NO_INTERLACE );
+ }
+ }
+ // -------------------------------------------------
+ else // support is not on all elements
+ {
+ // nb of nodes
+ getSupportNodes( support, med2ensIds );
+ NumberOfNodes = med2ensIds.size();
+ ensightGeomFile.addString( "coordinates" );
+ ensightGeomFile.addInt( NumberOfNodes );
+
+ // coordinates
+ vector<float> floatCoords( NumberOfNodes );
+ for ( j=0; j < SPACE_DIM; j++) { // loop on dimensions
+ medEnsIt = med2ensIds.begin();
+ if ( j < SpaceDimension ) {
+ const double *coordinate = _ptrMesh->getCoordinates(MED_FULL_INTERLACE) + j;
+ for (int i=0; i<NumberOfNodes; i++, ++medEnsIt )
+ floatCoords[ i ] = (float) coordinate[ (medEnsIt->first-1) * SpaceDimension];
+ }
+ else if ( j-1 < SpaceDimension ) {
+ for (int i=0; i<NumberOfNodes; i++)
+ floatCoords[ i ] = 0.;
+ }
+ ensightGeomFile.addReal( &floatCoords[0], NumberOfNodes );
+ }
+ // assign local node ids
+ for ( medEnsIt = med2ensIds.begin(), j=1; j<=NumberOfNodes; j++, ++medEnsIt )
+ medEnsIt->second = j;
+ }
+
+ // CONNECTIVITY Gold binary
+ // ===================================================================================
+ // loop on types
+ for (int i=0; i<nbTypes; i++)
+ {
+ const medGeometryElement medType = geoType[i];
+ const TEnSightElemType& ensightType = getEnSightType(medType);
+ const int numberOfCell = support->getNumberOfElements(medType);
+ int nbCellNodes = ensightType._medIndex.size();
+
+ // type name and nb cells
+ ensightGeomFile.addString( ensightType._name );
+ ensightGeomFile.addInt( numberOfCell );
+
+ vector<int> nodeIds;
+
+ // -------------------------------------------------
+ if ( support->isOnAllElements() )
+ {
+#ifdef ELEMENT_ID_GIVEN
+ // elem numbers
+ nodeIds.resize( numberOfCell );
+ for ( j = 1; j <= numberOfCell; j++)
+ nodeIds[ j-1 ] = j;
+ ensightGeomFile.addInt( nodeIds );
+
+ if ( entity != MED_NODE ) nodeIds.clear();
+#endif
+
+ if ( nbCellNodes > 1 ) // STANDARD ELEMENTS connectivity
+ {
+ connectivity = _ptrMesh->getConnectivity(MED_FULL_INTERLACE, MED_NODAL,
+ entity, medType);
+ nodeIds.reserve( numberOfCell * nbCellNodes);
+ for (j = 0 ; j < numberOfCell; j++, connectivity += nbCellNodes)
+ for (int k=0; k<nbCellNodes; k++)
+ nodeIds.push_back( connectivity[ ensightType._medIndex[k] ]);
+ ensightGeomFile.addInt( nodeIds );
+ }
+ else if ( entity == MED_NODE ) // NODES connectivity
+ {
+#if !defined(ELEMENT_ID_GIVEN)
+ nodeIds.resize( numberOfCell );
+ for ( j = 1; j <= numberOfCell; j++)
+ nodeIds[ j-1 ] = j;
+#endif
+ ensightGeomFile.addInt( nodeIds );
+ }
+ else if ( medType == MED_POLYGON ) // POLYGONs connectivity
+ {
+ connectivity = _ptrMesh->getPolygonsConnectivity(MED_NODAL, entity);
+ index = _ptrMesh->getPolygonsConnectivityIndex(MED_NODAL, entity);
+ int connLength = _ptrMesh->getPolygonsConnectivityLength(MED_NODAL, entity);
+ // number of nodes in each element
+ {
+ TIntOwner nbNodesInPoly( getNumbersByIndex( index, numberOfCell ));
+ ensightGeomFile.addInt( nbNodesInPoly.myValues, numberOfCell );
+ } // nbNodesInPoly is deleted here
+
+ // connectivity
+ ensightGeomFile.addInt( connectivity, connLength );
+ }
+ else // POLYHEDRA connectivity
+ {
+ connectivity = _ptrMesh->getPolyhedronConnectivity(MED_NODAL);
+ index = _ptrMesh->getPolyhedronIndex(MED_NODAL);
+ const int * fIndex = _ptrMesh->getPolyhedronFacesIndex();
+ int connLength = _ptrMesh->getPolyhedronConnectivityLength(MED_NODAL);
+ int nbFaces = _ptrMesh->getNumberOfPolyhedronFaces();
+ // nb of faces in each polyhedron
+ {
+ TIntOwner nbFacesInPoly( getNumbersByIndex( index, numberOfCell ));
+ ensightGeomFile.addInt( nbFacesInPoly.myValues, numberOfCell );
+ }
+ // number of nodes in each face
+ {
+ TIntOwner nbNodesInFace( getNumbersByIndex( fIndex, nbFaces ));
+ ensightGeomFile.addInt( nbNodesInFace.myValues, nbFaces );
+ }
+ // connectivity
+ ensightGeomFile.addInt( connectivity, connLength );
+ }
+ }
+ // -------------------------------------------------
+ else // support is not on all elements
+ {
+ const int *number = support->getNumber(medType);
+
+#ifdef ELEMENT_ID_GIVEN
+ ensightGeomFile.addInt( number, numberOfCell );
+#endif
+ if ( nbCellNodes > 1 ) // STANDARD ELEMENTS connectivity
+ {
+ connectivity = _ptrMesh->getConnectivity(MED_FULL_INTERLACE, MED_NODAL,
+ entity, MED_ALL_ELEMENTS);
+ index = _ptrMesh->getConnectivityIndex(MED_FULL_INTERLACE, entity);
+
+ nodeIds.reserve( numberOfCell * nbCellNodes);
+ for (j=0; j<numberOfCell; j++) {
+ int elem = number[j];
+ elemConnectivity = connectivity + index[elem-1]-1;
+ for (int k=0; k<nbCellNodes; k++)
+ {
+ int node = elemConnectivity[ ensightType._medIndex[k] ];
+ nodeIds.push_back( med2ensIds[ node ]);
+ }
+ }
+ ensightGeomFile.addInt( nodeIds );
+ }
+ else if ( entity == MED_NODE ) // NODES connectivity
+ {
+ nodeIds.resize( numberOfCell );
+ for ( j = 1; j <= numberOfCell; j++)
+ nodeIds[ j-1 ] = j;
+ ensightGeomFile.addInt( nodeIds );
+ }
+ else if ( medType == MED_POLYGON ) // POLYGONs connectivity
+ {
+ connectivity = _ptrMesh->getPolygonsConnectivity(MED_NODAL, entity);
+ index = _ptrMesh->getPolygonsConnectivityIndex(MED_NODAL, entity);
+ int connLength = _ptrMesh->getPolygonsConnectivityLength(MED_NODAL, entity);
+ int nbStdElems = _ptrMesh->getNumberOfElements(entity,MED_ALL_ELEMENTS);
+ // number of nodes in each element
+ {
+ TIntOwner nbNodesInPoly( getNumbersByIndex( index-nbStdElems, numberOfCell, number ));
+ ensightGeomFile.addInt( nbNodesInPoly.myValues, numberOfCell );
+ } // nbNodesInPoly is deleted here
+
+ // connectivity
+ nodeIds.reserve( connLength );
+ for ( j = 0; j < numberOfCell; ++j )
+ {
+ int elem = number[ j ] - nbStdElems;
+ elemConnectivity = connectivity + index[ elem-1 ]-1;
+ const int* connEnd = connectivity + index[ elem ]-1;
+ while ( elemConnectivity < connEnd )
+ nodeIds.push_back( med2ensIds[ *elemConnectivity++ ]);
+ }
+ ensightGeomFile.addInt( nodeIds );
+ }
+ else // POLYHEDRA connectivity
+ {
+ connectivity = _ptrMesh->getPolyhedronConnectivity(MED_NODAL);
+ index = _ptrMesh->getPolyhedronIndex(MED_NODAL);
+ const int * fIndex = _ptrMesh->getPolyhedronFacesIndex();
+ int connLength = _ptrMesh->getPolyhedronConnectivityLength(MED_NODAL);
+ int nbStdElems = _ptrMesh->getNumberOfElements(entity,MED_ALL_ELEMENTS);
+ // nb of faces in each polyhedron
+ {
+ TIntOwner nbFacesInPoly( getNumbersByIndex( index-nbStdElems, numberOfCell, number ));
+ ensightGeomFile.addInt( nbFacesInPoly.myValues, numberOfCell );
+ }
+ // number of nodes in each face
+ nodeIds.reserve( connLength );
+ for ( j = 0; j < numberOfCell; ++j )
+ {
+ int elem = number[ j ] - nbStdElems;
+ int f1 = index[ elem-1 ] - 1, f2 = index[ elem ] - 1;
+ int nbFaces = f2 - f1;
+ {
+ TIntOwner nbNodesInFace( getNumbersByIndex( &fIndex[ f1 ], nbFaces ));
+ ensightGeomFile.addInt( nbNodesInFace.myValues, nbFaces );
+ }
+ elemConnectivity = connectivity + fIndex[ f1 ] - 1;
+ const int* connEnd = connectivity + fIndex[ f2 ] - 1;
+ while ( elemConnectivity < connEnd )
+ nodeIds.push_back( med2ensIds[ *elemConnectivity++ ]);
+ }
+ // connectivity
+ ensightGeomFile.addInt( nodeIds );
+ }
+ }
+ }
+} // writePartGoldBinary()
+
+//================================================================================
+/*!
+ * \brief Write support as an EnSight Gold part
+ */
+//================================================================================
+
+void ENSIGHT_MESH_WRONLY_DRIVER::writePartGoldASCII(ofstream& ensightGeomFile,
+ const SUPPORT* support) const
+{
+ const int iw = 10;
+
+ // part number
+ int partNum = getPartNumber( support );
+ ensightGeomFile << "part" << endl
+ << setw(iw) << partNum << endl;
+ if ( !partNum )
+ throw MED_EXCEPTION ( LOCALIZED( STRING("Internal error: invalid part number")));
+
+ // group/mesh name
+ ensightGeomFile << support->getName() << endl;
+
+ // get geom types
+ medEntityMesh entity = support->getEntity();
+ int nbTypes = support->getNumberOfTypes();
+ const medGeometryElement* geoType = support->getTypes();
+
+ const int * connectivity = 0;
+ const int * elemConnectivity = 0;
+ const int * index = 0;
+ int j;
+
+ // COORDINATES Gold ASCII
+ // ===================================================================================
+ // In Ensight, coordinates of nodes of support elements are in MED_NO_INTERLACE mode.
+ // We are to write only nodes belonging to elements of the support and
+ // nodal connectivity should refer to these nodes.
+ map<int, int> med2ensIds;
+ map<int, int>::iterator medEnsIt;
+ int SpaceDimension = _ptrMesh->getSpaceDimension() ;
+ int NumberOfNodes = _ptrMesh->getNumberOfNodes() ;
+ string zeroStr = " 0.00000e+00";
+ // -----------------------------------
+ if ( support->isOnAllElements() )
+ {
+ // nb of nodes
+ ensightGeomFile << "coordinates" << endl
+ << setw(iw) << NumberOfNodes << endl ;
+
+ // coordinates
+ for (j=0; j<SPACE_DIM; j++) { // loop on dimensions
+ if ( j < SpaceDimension ) {
+ const double *coordinate = _ptrMesh->getCoordinates(MED_FULL_INTERLACE) + j;
+ for (int i=0; i<NumberOfNodes; i++, coordinate += SpaceDimension)
+ ensightGeomFile << setw(12) << (float) *coordinate << endl;
+ }
+ else {
+ for (int i=0; i<NumberOfNodes; i++)
+ ensightGeomFile << zeroStr << endl;
+ }
+ }
+ }
+ // -----------------------------------
+ else // support is not on all elements
+ {
+ // nb of nodes
+ getSupportNodes( support, med2ensIds );
+ NumberOfNodes = med2ensIds.size();
+ ensightGeomFile << "coordinates" << endl
+ << setw(iw) << NumberOfNodes << endl ;
+
+ // coordinates
+ for ( j=0; j<SPACE_DIM; j++) { // loop on dimensions
+ medEnsIt = med2ensIds.begin();
+ if ( j < SpaceDimension ) {
+ const double *coordinate = _ptrMesh->getCoordinates(MED_FULL_INTERLACE) + j;
+ for (int i=0; i<NumberOfNodes; i++, ++medEnsIt )
+ ensightGeomFile << setw(12)
+ << (float) coordinate[ (medEnsIt->first-1) * SpaceDimension] << endl;
+ }
+ else {
+ for (int i=0; i<NumberOfNodes; i++)
+ ensightGeomFile << zeroStr << endl;
+ }
+ }
+ // assign local node ids
+ for ( medEnsIt = med2ensIds.begin(), j=1; j<=NumberOfNodes; j++, ++medEnsIt )
+ medEnsIt->second = j;
+ }
+
+ // CONNECTIVITY Gold ASCII
+ // ===================================================================================
+ // loop on types
+ for (int i=0; i<nbTypes; i++)
+ {
+ const medGeometryElement medType = geoType[i];
+ const TEnSightElemType& ensightType = getEnSightType(medType);
+ const int numberOfCell = support->getNumberOfElements(medType);
+ int nbCellNodes = ensightType._medIndex.size();
+
+ // type name and nb cells
+ ensightGeomFile << ensightType._name << endl
+ << setw(iw) << numberOfCell << endl;
+
+ // -----------------------------------
+ if ( support->isOnAllElements() )
+ {
+#ifdef ELEMENT_ID_GIVEN
+ for ( j = 1; j <= numberOfCell; j++)
+ ensightGeomFile << setw(iw) << j << endl;
+#endif
+
+ if ( nbCellNodes > 1 ) // STANDARD ELEMENTS connectivity
+ {
+ connectivity = _ptrMesh->getConnectivity(MED_FULL_INTERLACE, MED_NODAL,
+ entity, medType);
+ for (j = 0 ; j < numberOfCell; j++, connectivity += nbCellNodes) {
+ for (int k=0; k<nbCellNodes; k++)
+ ensightGeomFile << setw(iw) << connectivity[ ensightType._medIndex[k] ];
+ ensightGeomFile << endl ;
+ }
+ }
+ else if ( entity == MED_NODE ) // NODES connectivity
+ {
+ for ( j = 1; j <= numberOfCell; j++)
+ ensightGeomFile << setw(iw) << j << endl;
+ }
+ else if ( medType == MED_POLYGON ) // POLYGONs connectivity
+ {
+ connectivity = _ptrMesh->getPolygonsConnectivity(MED_NODAL, entity);
+ index = _ptrMesh->getPolygonsConnectivityIndex(MED_NODAL, entity);
+ // number of nodes in each element
+ const int* ind = index;
+ for (j = 0 ; j < numberOfCell; j++, ++ind)
+ ensightGeomFile << setw(iw) << ( ind[1] - ind[0] ) << endl;
+
+ // connectivity
+ for (j = 0; j < numberOfCell; j++, ++index) {
+ nbCellNodes = index[1] - index[0];
+ for (int k=0; k<nbCellNodes; k++, ++connectivity)
+ ensightGeomFile << setw(iw) << *connectivity;
+ ensightGeomFile << endl;
+ }
+ }
+ else // POLYHEDRA connectivity
+ {
+ connectivity = _ptrMesh->getPolyhedronConnectivity(MED_NODAL);
+ index = _ptrMesh->getPolyhedronIndex(MED_NODAL);
+ const int * fIndex = _ptrMesh->getPolyhedronFacesIndex();
+ int nbFaces = _ptrMesh->getNumberOfPolyhedronFaces();
+ // nb of faces in each polyhedron
+ const int* ind = index;
+ for (j = 0 ; j < numberOfCell; j++, ++ind)
+ ensightGeomFile << setw(iw) << ( ind[1] - ind[0] ) << endl ;
+ // number of nodes in each face
+ ind = fIndex;
+ for (j = 0 ; j < nbFaces; j++, ++ind)
+ ensightGeomFile << setw(iw) << ( ind[1] - ind[0] ) << endl ;
+ // connectivity of each face
+ for (j = 0 ; j < nbFaces; j++, ++fIndex) {
+ int nbFaceNodes = fIndex[1] - fIndex[0];
+ for (int k=0; k<nbFaceNodes; k++, ++connectivity)
+ ensightGeomFile << setw(iw) << *connectivity;
+ ensightGeomFile << endl ;
+ }
+ }
+ }
+ // -----------------------------------
+ else // support is not on all elements
+ {
+ const int *number = support->getNumber(medType);
+
+#ifdef ELEMENT_ID_GIVEN
+ for ( j = 0; j < numberOfCell; j++)
+ ensightGeomFile << setw(iw) << number[j] << endl;
+#endif
+ if ( nbCellNodes > 1 ) // STANDARD ELEMENTS connectivity
+ {
+ connectivity = _ptrMesh->getConnectivity(MED_FULL_INTERLACE, MED_NODAL,
+ entity, MED_ALL_ELEMENTS);
+ index = _ptrMesh->getConnectivityIndex(MED_FULL_INTERLACE, entity);
+
+ for (j=0; j<numberOfCell; j++) {
+ int elem = number[j];
+ elemConnectivity = connectivity + index[elem-1]-1;
+ for (int k=0; k<nbCellNodes; k++) {
+ int node = elemConnectivity[ ensightType._medIndex[k] ];
+ ensightGeomFile << setw(iw) << med2ensIds[ node ];
+ }
+ ensightGeomFile << endl;
+ }
+ }
+ else if ( entity == MED_NODE ) // NODES connectivity
+ {
+ for (j=0; j<numberOfCell; j++) {
+ int node = med2ensIds[ number[j] ];
+ ensightGeomFile << setw(iw) << node << endl ;
+ }
+ }
+ else if ( medType == MED_POLYGON ) // POLYGONs connectivity
+ {
+ connectivity = _ptrMesh->getPolygonsConnectivity(MED_NODAL, entity);
+ index = _ptrMesh->getPolygonsConnectivityIndex(MED_NODAL, entity);
+ int nbStdElems = _ptrMesh->getNumberOfElements(entity,MED_ALL_ELEMENTS);
+ // number of nodes in each element
+ for (j = 0 ; j < numberOfCell; j++) {
+ int elem = number[j] - nbStdElems;
+ ensightGeomFile << setw(iw) << ( index[elem] - index[elem-1] ) << endl;
+ }
+ // connectivity
+ for ( j = 0; j < numberOfCell; ++j ) {
+ int elem = number[ j ] - nbStdElems;
+ elemConnectivity = connectivity + index[ elem-1 ]-1;
+ const int* connEnd = connectivity + index[ elem ]-1;
+ while ( elemConnectivity < connEnd )
+ ensightGeomFile << setw(iw) << med2ensIds[ *elemConnectivity++ ];
+ ensightGeomFile << endl;
+ }
+ }
+ else // POLYHEDRA connectivity
+ {
+ connectivity = _ptrMesh->getPolyhedronConnectivity(MED_NODAL);
+ index = _ptrMesh->getPolyhedronIndex(MED_NODAL);
+ const int * fIndex = _ptrMesh->getPolyhedronFacesIndex();
+ int nbStdElems = _ptrMesh->getNumberOfElements(entity,MED_ALL_ELEMENTS);
+ // nb of faces in each polyhedron
+ for (j = 0 ; j < numberOfCell; j++) {
+ int elem = number[j] - nbStdElems;
+ ensightGeomFile << setw(iw) << ( index[elem] - index[elem-1] ) << endl;
+ }
+ // number of nodes in each face
+ for ( j = 0; j < numberOfCell; ++j ) {
+ int elem = number[ j ] - nbStdElems;
+ int f1 = index[ elem-1 ], f2 = index[ elem ];
+ while ( f1 < f2 ) {
+ ensightGeomFile << setw(iw) << ( fIndex[f1] - fIndex[f1-1] ) << endl;
+ ++f1;
+ }
+ }
+ // connectivity of each face
+ for ( j = 0; j < numberOfCell; ++j ) {
+ int elem = number[ j ] - nbStdElems;
+ int f1 = index[ elem-1 ] - 1, f2 = index[ elem ] - 1;
+ while ( f1 < f2 ) {
+ int n1 = fIndex[f1]-1, n2 = fIndex[f1+1]-1;
+ ++f1;
+ while ( n1 < n2 )
+ ensightGeomFile << setw(iw) << connectivity[ n1++ ];
+ ensightGeomFile << endl ;
+ }
+ }
+ }
+ }
+ }
+} // writePartGoldASCII()
+
+//================================================================================
+/*!
+ * \brief Write support as an Ensight6 part
+ */
+//================================================================================
+
+void ENSIGHT_MESH_WRONLY_DRIVER::writePart6Binary(_BinaryFileWriter& ensightGeomFile,
+ const SUPPORT* support) const
+{
+ // part number
+ int partNum = getPartNumber( support );
+ ensightGeomFile.addString( STRING("part ") << partNum );
+ if ( !partNum )
+ throw MED_EXCEPTION ( LOCALIZED( STRING("Internal error: invalid part number")));
+
+ // group/mesh name
+ ensightGeomFile.addString( support->getName() );
+
+ // get geom types
+ medEntityMesh entity = support->getEntity();
+ int nbTypes = support->getNumberOfTypes();
+ const medGeometryElement* geoType = support->getTypes();
+
+ int j = 1;
+ const int * connectivity = 0;
+ if ( entity != MED_NODE )
+ connectivity = _ptrMesh->getConnectivity(MED_FULL_INTERLACE, MED_NODAL,
+ entity, MED_ALL_ELEMENTS);
+ const int * elemConnectivity = connectivity;
+
+ // CONNECTIVITY Ensight 6 binary
+ // ===================================================================================
+ // loop on types
+ for (int i=0; i<nbTypes; i++)
+ {
+ const medGeometryElement medType = geoType[i];
+ const TEnSightElemType& ensightType = getEnSightType(medType);
+ int nbCellNodes = ensightType._medIndex.size();
+ if ( nbCellNodes == 0 )
+ continue; // poly?
+
+ // type name and nb cells
+ int numberOfCell = support->getNumberOfElements(medType);
+ ensightGeomFile.addString( ensightType._name );
+ ensightGeomFile.addInt( numberOfCell );
+
+ vector<int> nodeIds;
+ // -------------------------------------------------
+ if ( support->isOnAllElements() )
+ {
+#ifdef ELEMENT_ID_GIVEN
+ nodeIds.resize( numberOfCell );
+ for ( j = 1; j <= numberOfCell; j++)
+ nodeIds[ j-1 ] = j;
+ ensightGeomFile.addInt( nodeIds );
+#endif
+ if ( entity == MED_NODE ) {
+#if !defined(ELEMENT_ID_GIVEN)
+ nodeIds.resize( numberOfCell * nbCellNodes);
+ for ( j = 1; j <= numberOfCell; j++)
+ nodeIds[ j-1 ] = j;
+#endif
+ }
+ else {
+ nodeIds.clear();
+ nodeIds.reserve( numberOfCell * nbCellNodes );
+ for (j = 0 ; j < numberOfCell; j++, elemConnectivity += nbCellNodes)
+ for (int k=0; k<nbCellNodes; k++)
+ nodeIds.push_back( elemConnectivity[ ensightType._medIndex[k] ]);
+ }
+ ensightGeomFile.addInt( nodeIds );
+ }
+ // -------------------------------------------------
+ else // support is not on all elements
+ {
+ const int *number = support->getNumber(medType);
+
+#ifdef ELEMENT_ID_GIVEN
+ ensightGeomFile.addInt( number, numberOfCell );
+#endif
+ if ( entity == MED_NODE ) {
+ ensightGeomFile.addInt( number, numberOfCell );
+ }
+ else {
+ const int* index = _ptrMesh->getConnectivityIndex(MED_FULL_INTERLACE, entity);
+
+ nodeIds.reserve( numberOfCell * nbCellNodes);
+ for (j=0; j<numberOfCell; j++) {
+ int elem = number[j];
+ elemConnectivity = connectivity + index[elem-1]-1;
+ for (int k=0; k<nbCellNodes; k++)
+ nodeIds.push_back( elemConnectivity[ ensightType._medIndex[k] ]);
+ }
+ ensightGeomFile.addInt( nodeIds );
+ }
+ }
+ } // loop on types
+
+} // writePart6Binary()
+
+//================================================================================
+/*!
+ * \brief Write support as an Ensight6 part
+ */
+//================================================================================
+
+void ENSIGHT_MESH_WRONLY_DRIVER::writePart6ASCII(ofstream& ensightGeomFile,
+ const SUPPORT* support) const
+{
+ const int iw = 8;
+
+ // part number
+ int partNum = getPartNumber( support );
+ ensightGeomFile << "part " << partNum << endl;
+ if ( !partNum )
+ throw MED_EXCEPTION ( LOCALIZED( STRING("Internal error: invalid part number")));
+
+ // group/mesh name
+ ensightGeomFile << support->getName() << endl;
+
+ // get geom types
+ medEntityMesh entity = support->getEntity();
+ int nbTypes = support->getNumberOfTypes();
+ const medGeometryElement* geoType = support->getTypes();
+
+ int j = 1;
+ const int * connectivity = 0;
+ if ( entity != MED_NODE )
+ connectivity = _ptrMesh->getConnectivity(MED_FULL_INTERLACE, MED_NODAL,
+ entity, MED_ALL_ELEMENTS);
+ const int * elemConnectivity = connectivity;
+
+ // CONNECTIVITY Ensight 6 ASCII
+ // ===================================================================================
+ // loop on types
+ for (int i=0; i<nbTypes; i++)
+ {
+ const medGeometryElement medType = geoType[i];
+ const TEnSightElemType& ensightType = getEnSightType(medType);
+ int nbCellNodes = ensightType._medIndex.size();
+ if ( nbCellNodes == 0 )
+ continue; // poly?
+
+ // type name and nb cells
+ int numberOfCell = support->getNumberOfElements(medType);
+ ensightGeomFile << ensightType._name << endl
+ << setw(iw) << numberOfCell << endl;
+
+ // -------------------------------------------------
+ if ( support->isOnAllElements() )
+ {
+ if ( entity == MED_NODE ) {
+ for ( j = 1; j <= numberOfCell; j++) {
+#ifdef ELEMENT_ID_GIVEN
+ ensightGeomFile << setw(iw) << j;
+#endif
+ ensightGeomFile << setw(iw) << j << endl;
+ }
+ }
+ else {
+ for (j = 1 ; j <= numberOfCell; j++, elemConnectivity += nbCellNodes) {
+#ifdef ELEMENT_ID_GIVEN
+ ensightGeomFile << setw(iw) << elem++;
+#endif
+ for (int k=0; k<nbCellNodes; k++)
+ {
+ ensightGeomFile << setw(iw) << elemConnectivity[ ensightType._medIndex[k] ];
+ }
+ ensightGeomFile << endl ;
+ }
+ }
+ }
+ // -------------------------------------------------
+ else // support is not on all elements
+ {
+ const int *number = support->getNumber(medType);
+ if ( entity == MED_NODE ) {
+ for (j=0; j<numberOfCell; j++) {
+ int node = number[j];
+#ifdef ELEMENT_ID_GIVEN
+ ensightGeomFile << setw(iw) << node;
+#endif
+ ensightGeomFile << setw(iw) << node << endl ;
+ }
+ }
+ else {
+ const int* index = _ptrMesh->getConnectivityIndex(MED_FULL_INTERLACE, entity);
+
+ for (j=0; j<numberOfCell; j++) {
+ int elem = number[j];
+#ifdef ELEMENT_ID_GIVEN
+ ensightGeomFile << setw(iw) << elem;
+#endif
+ elemConnectivity = connectivity + index[elem-1]-1;
+ for (int k=0; k<nbCellNodes; k++)
+ {
+ ensightGeomFile << setw(iw) << elemConnectivity[ ensightType._medIndex[k] ];
+ }
+ ensightGeomFile << endl ;
+ }
+ }
+ }
+ } // loop on types
+
+} // writePart6ASCII()
-void ENSIGHT_MESH_WRONLY_DRIVER::addSupport(SUPPORT *sup)
+//================================================================================
+/*!
+ * \brief Return nb of part to write
+ */
+//================================================================================
+
+int ENSIGHT_MESH_WRONLY_DRIVER::nbPartsToWrite() const
{
- _support.push_back(sup);
+ int nbParts = 0;
+ nbParts += (int) isToWriteEntity( MED_CELL, _ptrMesh );
+ nbParts += (int) isToWriteEntity( MED_FACE, _ptrMesh );
+ nbParts += (int) isToWriteEntity( MED_EDGE, _ptrMesh );
+
+ // all groups
+ for ( int ent = MED_CELL; ent < MED_ALL_ENTITIES; ++ent ) {
+ int nbGroups = _ptrMesh->getNumberOfGroups(medEntityMesh(ent));
+ nbParts += nbGroups;
+ }
+ return nbParts;
}
-// void ENSIGHT_MESH_WRONLY_DRIVER::writeSupport(SUPPORT * mySupport) const {
-// const char * LOC = "ENSIGHT_MESH_WRONLY_DRIVER::writeSupport(SUPPORT *) : " ;
-// BEGIN_OF(LOC) ;
-// MESSAGE(LOC << "Not yet implemented, acting on the object " << *mySupport);
-// END_OF(LOC) ;
-// }
+// ================================================================================
+// RDONLY
+// ================================================================================
-ENSIGHT_MESH_RDONLY_DRIVER::ENSIGHT_MESH_RDONLY_DRIVER() : ENSIGHT_MESH_DRIVER()
+ENSIGHT_MESH_RDONLY_DRIVER::ENSIGHT_MESH_RDONLY_DRIVER()
+ : ENSIGHT_MESH_DRIVER(), _indexInCaseFile(1)
{
- _ensightFile = new ifstream();
}
-ENSIGHT_MESH_RDONLY_DRIVER::ENSIGHT_MESH_RDONLY_DRIVER(const string & fileName, MESH * ptrMesh) : ENSIGHT_MESH_DRIVER(fileName,ptrMesh)
+ENSIGHT_MESH_RDONLY_DRIVER::ENSIGHT_MESH_RDONLY_DRIVER(const string & fileName,
+ MESH * ptrMesh,
+ int index)
+ : ENSIGHT_MESH_DRIVER(fileName,ptrMesh,RDONLY), _indexInCaseFile( index )
{
- _ensightFile = new ifstream();
}
-ENSIGHT_MESH_RDONLY_DRIVER::ENSIGHT_MESH_RDONLY_DRIVER(const ENSIGHT_MESH_RDONLY_DRIVER & driver) : ENSIGHT_MESH_DRIVER(driver)
+ENSIGHT_MESH_RDONLY_DRIVER::ENSIGHT_MESH_RDONLY_DRIVER(const ENSIGHT_MESH_RDONLY_DRIVER & driver) : ENSIGHT_MESH_DRIVER(driver), _indexInCaseFile( driver._indexInCaseFile )
{
- _ensightFile = new ifstream();
}
ENSIGHT_MESH_RDONLY_DRIVER::~ENSIGHT_MESH_RDONLY_DRIVER()
{
- delete _ensightFile ;
}
GENDRIVER * ENSIGHT_MESH_RDONLY_DRIVER::copy() const
return new ENSIGHT_MESH_RDONLY_DRIVER(*this) ;
}
-void ENSIGHT_MESH_RDONLY_DRIVER::openConst() const {
+void ENSIGHT_MESH_RDONLY_DRIVER::write() const throw (MEDEXCEPTION)
+{
+ throw MEDEXCEPTION("ENSIGHT_MESH_RDONLY_DRIVER::write : Can't write with a RDONLY driver !");
+}
- const char * LOC ="ENSIGHT_MESH_RDONLY_DRIVER::open() : ";
+void ENSIGHT_MESH_RDONLY_DRIVER::merge ( const GENDRIVER& driver )
+{
+ _CaseFileDriver_User::merge( driver );
- BEGIN_OF(LOC);
+ const ENSIGHT_MESH_RDONLY_DRIVER* other =
+ dynamic_cast< const ENSIGHT_MESH_RDONLY_DRIVER* >( &driver );
+ if ( other ) {
+ if ( _indexInCaseFile < other->_indexInCaseFile )
+ _indexInCaseFile = other->_indexInCaseFile;
+ }
+}
- if ( _fileName == "" )
- throw MED_EXCEPTION ( LOCALIZED( STRING(LOC)
- << "_fileName is |\"\"|, please set a correct fileName before calling open()"
- )
- );
+//================================================================================
+/*!
+ * \brief Read mesh in all supported formats
+ */
+//================================================================================
- if (!(*_ensightFile).is_open())
- (*_ensightFile).open(_fileName.c_str()) ;
+void ENSIGHT_MESH_RDONLY_DRIVER::read() throw (MEDEXCEPTION)
+{
+ const char * LOC = "ENSIGHT_MESH_RDONLY_DRIVER::read() : " ;
+ BEGIN_OF_MED(LOC);
- if (!(*_ensightFile))
- throw MED_EXCEPTION ( LOCALIZED( STRING(LOC) << "Could not open main Ensight file "
- << _fileName)
- );
- END_OF(LOC);
-}
+ openConst(false); // check if can read case file
-void ENSIGHT_MESH_RDONLY_DRIVER::closeConst() const {
-
- const char * LOC = "ENSIGHT_MESH_RDONLY_DRIVER::close() : ";
- BEGIN_OF(LOC);
-
- (*_ensightFile).close();
-
- if (!(*_ensightFile))
- throw MED_EXCEPTION ( LOCALIZED( STRING(LOC) << "Could not close main Ensight file "
- << _fileName)
- );
- END_OF(LOC);
-}
+ _CaseFileDriver caseFile( getCaseFileName(), this);
+ caseFile.read();
+ caseFile.setDataFileName( _indexInCaseFile, this ); // data from Case File is passed here
-void ENSIGHT_MESH_RDONLY_DRIVER::write() const throw (MEDEXCEPTION) {
- throw MEDEXCEPTION("ENSIGHT_MESH_RDONLY_DRIVER::write : Can't write with a RDONLY driver !");
-}
+ openConst(true); // check if can read data file
-void ENSIGHT_MESH_RDONLY_DRIVER::read() {
+ cout << "-> Entering into the geometry file " << getDataFileName() << endl ;
- const char * LOC = "ENSIGHT_MESH_RDONLY_DRIVER::read() : " ;
- BEGIN_OF(LOC);
+ _InterMed* imed = new _InterMed();
+ imed->_medMesh = getMesh();
+ imed->_isOwnMedMesh = false;
+ imed->_needSubParts = ( caseFile.getNbVariables() > 0 );
+ imed->groupes.reserve(1000);
- openConst() ;
+ // to let field drivers know eventual indices of values
+ if ( imed->_needSubParts )
+ setInterData( imed );
- string type_Ensight[15] = {
- "point" , "bar2", "bar3" , "tria3" , "tria6" , "quad4" , "quad8" , "tetra4" , "tetra10" , "pyramid5" ,
- "pyramid13" , "hexa8" , "hexa20" , "penta6" , "penta15" };
- int nb_edge[15] = { 1,2,3,3,6,4,8,4,10,5,13,8,20,6,15 };
-
- vector<string> type_read ;
- vector<int> nbcell_read ;
- vector< vector <int> > total_conn ;
- vector<double> var ;
-
- //int number_of_geom ;
- string geom_namefile ;
- string field_namefile ;
- string mot_lu ;
- int geom_given = 0 ;
- int num_coordinate ;
- string type_cell ;
- int number_of_cell ;
- int conn ;
- //----- ?
- int SpaceDimension = 3 ;
- //----- ?
-
- int NumberOfTypes = 0 ;
- int NumberOfNodes ;
- int iType ;
-
- string mesh_read_name = "EnsightMesh"; // defaut name for the mesh
-
- // recuperation des arguments du fichier ensight case
- // --------------------------------------------------
-
- // In this release, the following options are available :
- // For GEOMETRY -> model:
-
- cout << "****************** READ **************** starting " << endl ;
-
- ifstream ensightCaseFile(_fileName.c_str(),ios::in);
- cout << "Ensight case file name to read " << _fileName << endl ;
- string diren = getDirName((char*)_fileName.c_str());
-
- if (ensightCaseFile.is_open() )
- {
- while ( ensightCaseFile >> mot_lu )
- {
- if ( mot_lu == "GEOMETRY" ) {
- cout << "geometry detected" << endl ;
- while ( ensightCaseFile >> mot_lu ){
- if ( mot_lu == "model:" ) {
-// ensightCaseFile >> number_of_geom ;
-// cout << "number of geometries " << number_of_geom << endl ;
- ensightCaseFile >> mot_lu ;
- geom_namefile = mot_lu;
- cout << "name of geometry file : " << geom_namefile << endl ;
- break ;
- }
- }
- }
- }
+ if ( isBinaryDataFile( getDataFileName() )) // binary
+ {
+ if ( isGoldFormat() ) // Gold
+ {
+ readGoldBinary( *imed );
}
- else
+ else // EnSight6
+ {
+ read6Binary( *imed );
+ }
+ }
+ else // ASCII
+ {
+ if ( isGoldFormat() ) // Gold
+ {
+ readGoldASCII( *imed );
+ }
+ else // EnSight6
{
- cout << "Error : requested file " << ensightCaseFile << " not found " << endl;
- exit( EXIT_FAILURE );
+ read6ASCII( *imed );
}
-
- // chargement des noeuds et connectivites necessaires depuis le fichier ensight geom
- // ---------------------------------------------------------------------------------
+ }
- _ptrMesh->_name = mesh_read_name ;
+ if ( _isMadeByMed && !imed->groupes.empty() ) {
+ _ptrMesh->_name = imed->groupes[0].nom;
+ imed->groupes[0].nom = "SupportOnAll_";
+ imed->groupes[0].nom += entNames[MED_CELL];
+ }
+ else {
+ _ptrMesh->_name = theDefaultMeshName;
+ }
+ _ptrMesh->_spaceDimension = SPACE_DIM;
+ _ptrMesh->_meshDimension = _ptrMesh->_spaceDimension;
+ _ptrMesh->_numberOfNodes = imed->points.size() - imed->nbMerged( MED_POINT1 );
+ _ptrMesh->_isAGrid = 0;
+ _ptrMesh->_coordinate = imed->getCoordinate();
+
+ //Construction des groupes
+ imed->getGroups(_ptrMesh->_groupCell,
+ _ptrMesh->_groupFace,
+ _ptrMesh->_groupEdge,
+ _ptrMesh->_groupNode, _ptrMesh);
+
+ _ptrMesh->_connectivity = imed->getConnectivity();
+
+ _ptrMesh->createFamilies();
+
+ // add attributes to families
+ set<string> famNames;
+ for (medEntityMesh entity=MED_CELL; entity<MED_ALL_ENTITIES; ++entity)
+ {
+ int i, nb = _ptrMesh->getNumberOfFamilies(entity);
+ for ( i = 1; i <= nb; ++i ) {
+ FAMILY* f = const_cast<FAMILY*>( _ptrMesh->getFamily( entity, i ));
+ f->setNumberOfAttributes( 1 );
+ int* attIDs = new int[1];
+ attIDs[0] = 1;
+ f->setAttributesIdentifiers( attIDs );
+ int* attVals = new int[1];
+ attVals[0] = 1;
+ f->setAttributesValues( attVals );
+ string* attDescr = new string[1];
+ attDescr[0] = "med_family";
+ f->setAttributesDescriptions( attDescr );
+ if ( f->getName().length() > 31 ) // limit a name length
+ f->setName( STRING("FAM_") << f->getIdentifier());
+ // check if family is on the whole mesh entity
+ if (_ptrMesh->getNumberOfElements( entity, MED_ALL_ELEMENTS ) ==
+ f->getNumberOfElements( MED_ALL_ELEMENTS ))
+ {
+ f->setAll( true );
+ *(f->getnumber()) = MEDSKYLINEARRAY();
+ }
+ // setAll() for groups
+ nb = _ptrMesh->getNumberOfGroups(entity);
+ for ( i = 1; i <= nb; ++i ) {
+ GROUP * g = const_cast<GROUP*>( _ptrMesh->getGroup( entity, i ));
+ if (_ptrMesh->getNumberOfElements( entity, MED_ALL_ELEMENTS ) ==
+ g->getNumberOfElements( MED_ALL_ELEMENTS ))
+ {
+ g->setAll( true );
+ *(g->getnumber()) = MEDSKYLINEARRAY();
+ }
+ }
+ }
+ }
- string cgeom_namefile;
- if( diren.length() > 0 )
- cgeom_namefile = diren + '/' + geom_namefile;
- else
- cgeom_namefile = geom_namefile;
- cout << "-> Entering into the geometry file " << geom_namefile << endl ;
- ifstream ensightGeomFile(cgeom_namefile.c_str(),ios::in);
- if (ensightGeomFile.is_open() )
- {
- while ( ensightGeomFile >> mot_lu ){
- if ( mot_lu == "given" ) geom_given=1 ;
- if ( mot_lu == "coordinates" ) {
-//---------------- Nodes part --------------------------------------------------
- ensightGeomFile >> NumberOfNodes ;
- cout << "-> loading " << NumberOfNodes << " coordinates " << endl ;
- int NumberOfCoordinates = NumberOfNodes*SpaceDimension ;
- double* Coordinates = new double[NumberOfCoordinates];
- int iCoord = 0 ;
- //cout << "-> geom given " << geom_given << endl ;
- for ( int i=0 ; i < NumberOfNodes ; i++) {
- if( geom_given) ensightGeomFile >> setw(8) >> num_coordinate ;
- ensightGeomFile >> setw(12) >> Coordinates[iCoord] ;
- ensightGeomFile >> setw(12) >> Coordinates[iCoord+1] ;
- ensightGeomFile >> setw(12) >> Coordinates[iCoord+2] ;
-// cout << "coordinate read " << num_coordinate << " : x = " << Coordinates[iCoord] << " y = " << Coordinates[iCoord+1] << " z = " << Coordinates[iCoord+2] << endl ;
- iCoord+=3 ;
- }
- _ptrMesh->_spaceDimension = SpaceDimension;
- _ptrMesh->_numberOfNodes = NumberOfNodes;
- _ptrMesh->_coordinate = new COORDINATE(SpaceDimension,NumberOfNodes,MED_EN::MED_FULL_INTERLACE);
- _ptrMesh->_coordinate->setCoordinates(MED_EN::MED_FULL_INTERLACE,Coordinates);
- _ptrMesh->_coordinate->setCoordinatesSystem("CARTESIAN");
- delete [] Coordinates;
- }
- else if ( mot_lu == "part" ) {
-//---------------- Connectivities part --------------------------------------------
- while ( ensightGeomFile >> mot_lu ){
- for ( int j = 0 ; j < 15 ; j++){
- if( mot_lu == type_Ensight[j] ) {
- NumberOfTypes+=1;
- iType=NumberOfTypes-1 ;
- total_conn.resize(NumberOfTypes) ;
- type_read.push_back(mot_lu) ;
- ensightGeomFile >> number_of_cell ;
- nbcell_read.push_back(number_of_cell) ;
- total_conn[iType].resize(nb_edge[j]*number_of_cell);
- cout << "-> loading " << number_of_cell << " cells connectivities of type " << type_Ensight[j] << " (" << nb_edge[j]*number_of_cell << ") values " << endl ;
- for ( int k=0 ; k < nb_edge[j]*number_of_cell ; k++ ) {
- ensightGeomFile >> setw(8) >> conn ;
- total_conn[iType][k]=conn ;
- // cout << " connectivitie " << k << " read = " << total_conn[iType][k] << endl ;
- }
- }
- }
- }
- }
- }
- }
- // for compilation on WNT
-#ifndef WNT
- medGeometryElement classicalTypesCell[NumberOfTypes];
- int nbOfClassicalTypesCell[NumberOfTypes];
-#else // massive with zero size can't exist on Win32
- medGeometryElement* classicalTypesCell = new medGeometryElement(NumberOfTypes);
- int* nbOfClassicalTypesCell = new int(NumberOfTypes);
-#endif
- int ind=0 ;
- for (int k=0 ; k<NumberOfTypes ; k++){
- for (int j=0 ; j<15 ; j++)
- if(type_read[k] == type_Ensight[j] ){
- switch ( j+1 )
- {
- case 1 : {classicalTypesCell[ind] = MED_EN::MED_POINT1 ; nbOfClassicalTypesCell[ind] = nbcell_read[k] ; ind++ ; break ;}
- case 2 : {classicalTypesCell[ind] = MED_EN::MED_SEG2 ; nbOfClassicalTypesCell[ind] = nbcell_read[k] ; ind++ ; break ;}
- case 3 : {classicalTypesCell[ind] = MED_EN::MED_SEG3 ; nbOfClassicalTypesCell[ind] = nbcell_read[k] ; ind++ ; break ;}
- case 4 : {classicalTypesCell[ind] = MED_EN::MED_TRIA3 ; nbOfClassicalTypesCell[ind] = nbcell_read[k] ; ind++ ; break ;}
- case 5 : {classicalTypesCell[ind] = MED_EN::MED_TRIA6 ; nbOfClassicalTypesCell[ind] = nbcell_read[k] ; ind++ ; break ;}
- case 6 : {classicalTypesCell[ind] = MED_EN::MED_QUAD4 ; nbOfClassicalTypesCell[ind] = nbcell_read[k] ; ind++ ; break ;}
- case 7 : {classicalTypesCell[ind] = MED_EN::MED_QUAD8 ; nbOfClassicalTypesCell[ind] = nbcell_read[k] ; ind++ ; break ;}
- case 8 : {classicalTypesCell[ind] = MED_EN::MED_TETRA4 ; nbOfClassicalTypesCell[ind] = nbcell_read[k] ; ind++ ; break ;}
- case 9 : {classicalTypesCell[ind] = MED_EN::MED_TETRA10 ; nbOfClassicalTypesCell[ind] = nbcell_read[k] ; ind++ ; break ;}
- case 10 : {classicalTypesCell[ind] = MED_EN::MED_PYRA5 ; nbOfClassicalTypesCell[ind] = nbcell_read[k] ; ind++ ; break ;}
- case 11 : {classicalTypesCell[ind] = MED_EN::MED_PYRA13 ; nbOfClassicalTypesCell[ind] = nbcell_read[k] ; ind++ ; break ;}
- case 12 : {classicalTypesCell[ind] = MED_EN::MED_HEXA8 ; nbOfClassicalTypesCell[ind] = nbcell_read[k] ; ind++ ; break ;}
- case 13 : {classicalTypesCell[ind] = MED_EN::MED_HEXA20 ; nbOfClassicalTypesCell[ind] = nbcell_read[k] ; ind++ ; break ;}
- case 14 : {classicalTypesCell[ind] = MED_EN::MED_PENTA6 ; nbOfClassicalTypesCell[ind] = nbcell_read[k] ; ind++ ; break ;}
- case 15 : {classicalTypesCell[ind] = MED_EN::MED_PENTA15 ; nbOfClassicalTypesCell[ind] = nbcell_read[k] ; ind++ ; break ;}
- default : break ;
- }
- }
- }
-
- _ptrMesh->_connectivity = new CONNECTIVITY(NumberOfTypes,MED_EN::MED_CELL);
- _ptrMesh->_connectivity->setGeometricTypes(classicalTypesCell,MED_EN::MED_CELL);
- int * Count = new int[NumberOfTypes+1] ;
- Count[0]=1 ;
- for (int i=0; i<NumberOfTypes; i++)
- Count[i+1]=Count[i]+nbOfClassicalTypesCell[i] ;
- _ptrMesh->_connectivity->setCount(Count,MED_EN::MED_CELL) ;
- delete[] Count ;
- int MeshDimension ;
- MeshDimension = 2 ;
- for (int k=0 ; k<NumberOfTypes ; k++){
- for (int j=0 ; j<15 ; j++)
- if(type_read[k] == type_Ensight[j] && j>6 ) MeshDimension = 3 ;
- }
- _ptrMesh->_meshDimension = MeshDimension;
- _ptrMesh->_connectivity->setEntityDimension(MeshDimension);
-
- for (int k = 0 ; k < NumberOfTypes ; k++) {
- int nb_connectivities = total_conn[k].size();
- int *connectivities = new int[nb_connectivities];
- for (int itab=0 ; itab < nb_connectivities ; itab++) connectivities[itab]=total_conn[k][itab] ;
- for (int j=0 ; j<15 ; j++) {
- if( type_read[k] == type_Ensight[j] ) {
- switch ( j+1 )
- {
- case 1 : {_ptrMesh->_connectivity->setNodal(connectivities,MED_EN::MED_CELL,MED_EN::MED_POINT1) ; break ;}
- case 2 : {_ptrMesh->_connectivity->setNodal(connectivities,MED_EN::MED_CELL,MED_EN::MED_SEG2) ; break ;}
- case 3 : {_ptrMesh->_connectivity->setNodal(connectivities,MED_EN::MED_CELL,MED_EN::MED_SEG3) ; break ;}
- case 4 : {_ptrMesh->_connectivity->setNodal(connectivities,MED_EN::MED_CELL,MED_EN::MED_TRIA3) ; break ;}
- case 5 : {_ptrMesh->_connectivity->setNodal(connectivities,MED_EN::MED_CELL,MED_EN::MED_TRIA6) ; break ;}
- case 6 : {_ptrMesh->_connectivity->setNodal(connectivities,MED_EN::MED_CELL,MED_EN::MED_QUAD4) ; break ;}
- case 7 : {_ptrMesh->_connectivity->setNodal(connectivities,MED_EN::MED_CELL,MED_EN::MED_QUAD8) ; break ;}
- case 8 : {_ptrMesh->_connectivity->setNodal(connectivities,MED_EN::MED_CELL,MED_EN::MED_TETRA4) ; break ;}
- case 9 : {_ptrMesh->_connectivity->setNodal(connectivities,MED_EN::MED_CELL,MED_EN::MED_TETRA10) ; break ;}
- case 10 : {_ptrMesh->_connectivity->setNodal(connectivities,MED_EN::MED_CELL,MED_EN::MED_PYRA5) ; break ;}
- case 11 : {_ptrMesh->_connectivity->setNodal(connectivities,MED_EN::MED_CELL,MED_EN::MED_PYRA13) ; break ;}
- case 12 : {_ptrMesh->_connectivity->setNodal(connectivities,MED_EN::MED_CELL,MED_EN::MED_HEXA8) ; break ;}
- case 13 : {_ptrMesh->_connectivity->setNodal(connectivities,MED_EN::MED_CELL,MED_EN::MED_HEXA20) ; break ;}
- case 14 : {_ptrMesh->_connectivity->setNodal(connectivities,MED_EN::MED_CELL,MED_EN::MED_PENTA6) ; break ;}
- case 15 : {_ptrMesh->_connectivity->setNodal(connectivities,MED_EN::MED_CELL,MED_EN::MED_PENTA15) ; break ;}
- default : break ;
- }
- }
- }
- delete [] connectivities;
- }
-
-// cout << "Impression de _ptrmesh dans EnsightMeshDriver: " << endl;
-// cout << *_ptrMesh ;
-
- END_OF(LOC);
+ END_OF_MED(LOC);
}
+//================================================================================
+/*!
+ * \brief Read mesh in Gold ASCII format
+ */
+//================================================================================
+
+void ENSIGHT_MESH_RDONLY_DRIVER::readGoldASCII(_InterMed & imed)
+{
+ const char * LOC ="ENSIGHT_MESH_RDONLY_DRIVER::readGoldASCII() : ";
+ BEGIN_OF_MED(LOC);
+
+ _ASCIIFileReader geoFile( getDataFileName() );
+
+ if ( isSingleFileMode() ) {
+ int curTimeStep = 1;
+ while ( curTimeStep++ < getIndexInDataFile() ) {
+ while ( !geoFile.isTimeStepEnd())
+ geoFile.getLine();
+ }
+ while ( !geoFile.isTimeStepBeginning() )
+ geoFile.getLine();
+ }
+ // ----------------------
+ // Read mesh description
+ // ----------------------
+ {
+ string descriptionLine1 = geoFile.getLine();
+ string descriptionLine2 = geoFile.getLine();
+
+ // find out if the file was created by MED driver
+ int prefixSize = strlen( theDescriptionPrefix );
+ _isMadeByMed = ( descriptionLine1.substr(0, prefixSize ) == theDescriptionPrefix );
+
+ if ( _isMadeByMed )
+ descriptionLine1 = descriptionLine1.substr( prefixSize );
+ _ptrMesh->setDescription( descriptionLine1 + descriptionLine2 );
+ }
+
+ // ----------------------------------------
+ // Find out presence of node/elem numbers
+ // ----------------------------------------
+
+ // EnSight User Manual (for v8) says:
+// You do not have to assign node IDs. If you do, the element connectivities are
+// based on the node numbers. If you let EnSight assign the node IDs, the nodes
+// are considered to be sequential starting at node 1, and element connectivity is
+// done accordingly. If node IDs are set to off, they are numbered internally;
+// however, you will not be able to display or query on them. If you have node
+// IDs in your data, you can have EnSight ignore them by specifying "node id
+// ignore." Using this option may reduce some of the memory taken up by the
+// Client and Server, but display and query on the nodes will not be available.
+
+ // read "node|element id <off|given|assign|ignore>"
+ geoFile.getWord(); geoFile.getWord();
+ string nodeIds = geoFile.getWord();
+ geoFile.getWord(); geoFile.getWord();
+ string elemIds = geoFile.getWord();
+
+ bool haveNodeIds = ( nodeIds == "given" || nodeIds == "ignore" );
+ bool haveElemIds = ( elemIds == "given" || elemIds == "ignore" );
+
+ // extents: xmin xmax ymin ymax zmin zmax
+ vector<double> extents;
+ geoFile.toNextLine();
+ if ( strncmp( "extents", geoFile.getCurrentPtr(), 7 ) == 0 ) {
+ geoFile.skip( /*width =*/ 7, /*nbLines =*/ 1 );
+ extents.reserve( 6 );
+ while ( extents.size() < extents.capacity() )
+ extents.push_back( geoFile.getReal() );
+ }
+
+ typedef map<int,_noeud>::iterator INoeud;
+ map<int,_noeud> & points = imed.points;
+ INoeud firstNode;
+
+ _groupe * partGroupe = 0;
+ int partNum = 0, nbParts = 0;
+
+ while ( !geoFile.isTimeStepEnd() )
+ {
+ string word, restLine, line = geoFile.getLine();
+ TStrTool::split( line, word, restLine );
+
+ const TEnSightElemType & partType = getEnSightType( word );
+ if ( partType._medType != MED_ALL_ELEMENTS )
+ {
+ // Unstructured element type encounters
+ // --------------------------------------
+ int nbElemNodes = partType._medType % 100;
+ int nbElems = geoFile.getInt(); // ne
+ bool isGhost = isGhostType( word );
+ int nodeShift = points.empty() ? 0 : points.rbegin()->first;
+
+ // read element ids
+ vector<int> elemIds;
+ if ( haveElemIds ) {
+ elemIds.reserve( nbElems );
+ while ( elemIds.size() < nbElems )
+ elemIds.push_back( geoFile.getInt() ); // id_e
+ }
+ if ( isGhost ) { // do not store ghost elements (?)
+ int nbNodes = nbElems * nbElemNodes;
+ if ( partType._name == "nsided" ) // polygons
+ {
+ for ( int i = 0; i < nbElems; ++i )
+ nbNodes += geoFile.getInt();
+ geoFile.skip( nbNodes * INT_WIDTH_GOLD, /*nbLines = */nbElems );
+ }
+ else if ( partType._name == "nfaced" ) // polyhedrons
+ {
+ int nbFaces = 0;
+ for ( int i = 0; i < nbElems; ++i )
+ nbFaces += geoFile.getInt();
+ for ( int f = 0; f < nbFaces; ++f )
+ nbNodes += geoFile.getInt();
+ geoFile.skip( nbNodes * INT_WIDTH_GOLD, /*nbLines = */nbFaces );
+ }
+ else // standard types
+ {
+ geoFile.skip( nbNodes, nbElemNodes, INT_WIDTH_GOLD );
+ }
+ continue;
+ }
+
+ // add a group corresponding to subPart (geoType)
+ imed.groupes.push_back(_groupe());
+ _groupe & groupe = imed.groupes.back();
+ groupe.mailles.resize( nbElems );
+
+ // find out if "coordinates" has already been encountered
+ _SubPartDesc coordDesc( partNum , "coordinates");
+ map< _SubPartDesc, _SubPart >::iterator descPart =
+ imed._subPartDescribed.find( coordDesc );
+ bool haveCoords = ( descPart != imed._subPartDescribed.end() );
+ if ( haveCoords ) {
+ firstNode = descPart->second.myFirstNode;
+ nodeShift -= descPart->second.myNbNodes;
+ }
+
+ // read poly element data
+ bool isPoly = ( !nbElemNodes );
+ vector<int> nbElemNodesVec( 1, nbElemNodes);
+ vector<int> nbElemFaces, nbFaceNodes;
+ if ( partType._name == "nsided" ) // polygons
+ {
+ nbElemNodesVec.resize( nbElems );
+ for ( int i = 0; i < nbElems; ++i )
+ nbElemNodesVec[ i ] = geoFile.getInt(); // npi
+ }
+ else if ( partType._name == "nfaced" ) // polyhedrons
+ {
+ nbElemFaces.resize( nbElems );
+ nbElemNodesVec.resize( nbElems );
+ int totalNbFaces = 0;
+ for ( int i = 0; i < nbElems; ++i )
+ totalNbFaces += ( nbElemFaces[ i ] = geoFile.getInt() ); // nf_ei
+
+ nbFaceNodes.resize( totalNbFaces );
+ vector<int>::iterator nbFN = nbFaceNodes.begin();
+ for ( int i = 0; i < nbElems; ++i ) {
+ nbElemNodesVec[ i ] = 0;
+ for ( int nbFaces = nbElemFaces[ i ]; nbFaces; --nbFaces, ++nbFN )
+ nbElemNodesVec[ i ] += ( *nbFN = geoFile.getInt() ); // np(f_ei)
+ }
+ }
+ // iterator returning nbElemNodes for standard elems and
+ // next value from nbElemNodesVec for poly elements
+ _ValueIterator<int> nbElemNodesIt( & nbElemNodesVec[0], isPoly ? 1 : 0);
+
+ // iterator returning values form partType._medIndex for standard elems
+ // and node index (n) for poly elements
+ int n;
+ _ValueIterator<int> medIndexIt( isPoly ? &n : &partType._medIndex[0],
+ isPoly ? 0 : 1);
+ // read connectivity
+ _maille ma( partType._medType, nbElemNodes );
+ ma.sommets.resize( nbElemNodes );
+ INoeud node;
+ for ( int i = 0; i < nbElems; ++i ) {
+ _ValueIterator<int> medIndex = medIndexIt;
+ nbElemNodes = nbElemNodesIt.next();
+ if ( ma.sommets.size() != nbElemNodes )
+ ma.sommets.resize( nbElemNodes );
+ for ( n = 0; n < nbElemNodes; ++n ) {
+ int nodeID = geoFile.getInt(); // nn_ei
+ if ( haveCoords )
+ node = points.find( nodeID + nodeShift );
+ else
+ node = points.insert( make_pair( nodeID + nodeShift, _noeud())).first;
+ ma.sommets[ medIndex.next() ] = node;
+ }
+ if ( haveElemIds )
+ ma.setOrdre( elemIds[ i ] );
+ groupe.mailles[i] = imed.insert(ma);
+ }
+ // store nb nodes in polyhedron faces
+ if ( !nbFaceNodes.empty() ) {
+ const int* nbFaceNodesPtr = & nbFaceNodes[0];
+ for ( int i = 0; i < nbElems; ++i ) {
+ vector<int> & nbNodesByFace = imed.polyherdalNbFaceNodes[ &groupe.maille( i ) ];
+ nbNodesByFace.assign( nbFaceNodesPtr, nbFaceNodesPtr + nbElemFaces[ i ] );
+ nbFaceNodesPtr += nbElemFaces[ i ];
+ }
+ }
+ // create subPart for "coordinates"
+ if ( !haveCoords ) {
+ _SubPart & coordSubPart = imed._subPartDescribed[ coordDesc ];
+ coordSubPart.myFirstNode = points.insert( make_pair( 1 + nodeShift, _noeud())).first;
+ }
+ // add subPart group to part group
+ int groupeIndex = imed.groupes.size();
+ partGroupe->groupes.push_back( groupeIndex );
+
+ // create subPart
+ _SubPart subPart( partNum, partType._name );
+ subPart.myNbCells = nbElems;
+ subPart.myCellGroupIndex = groupeIndex;
+ imed.addSubPart( subPart );
+ }
+ else if ( word == "coordinates" )
+ {
+ // Local node coordinates of a part
+ // ------------------------------------
+ int nbNodes = geoFile.getInt(); // nn
+
+ // read node ids
+ vector<int> nodeIds;
+ if ( haveNodeIds ) {
+ nodeIds.reserve( nbNodes );
+ while ( nodeIds.size() < nbNodes )
+ nodeIds.push_back( geoFile.getInt() ); // id_n
+ }
+
+ // find out if "coordinates" has already been add at reading connectivity
+ _SubPartDesc coordDesc( partNum , "coordinates");
+ map< _SubPartDesc, _SubPart >::iterator descPart =
+ imed._subPartDescribed.find( coordDesc );
+ bool haveCoords = ( descPart != imed._subPartDescribed.end() );
+
+ if ( haveCoords ) {
+ // check that all nodes have been added
+ firstNode = descPart->second.myFirstNode;
+ descPart->second.myNbNodes = nbNodes;
+ INoeud inoeud = firstNode, inoEnd = points.end();
+ int id = inoeud->first, idEnd = id + nbNodes;
+ for ( ; id < idEnd; ++id ) {
+ if ( inoeud == inoEnd || inoeud->first > id ) {
+ INoeud in = points.insert( inoeud, make_pair( id, _noeud() ));
+ in->second.number = id;
+ in->second.coord.resize( SPACE_DIM );
+ } else {
+ ++inoeud;
+ }
+ }
+ }
+ else {
+ // add nodes
+ int nodeShift = points.empty() ? 0 : points.rbegin()->first;
+ for ( int iNode = 1; iNode <= nbNodes; ++iNode ) {
+ INoeud inoeud = points.insert( points.end(), make_pair( iNode + nodeShift, _noeud()));
+ inoeud->second.number = inoeud->first;
+ inoeud->second.coord.resize( SPACE_DIM );
+ }
+ firstNode = points.find( 1 + nodeShift );
+ // create "coordinates" subPart
+ _SubPart & subPart = imed._subPartDescribed[ coordDesc ];
+ subPart.myNbNodes = nbNodes;
+ subPart.myFirstNode = firstNode;
+ }
+
+ // read coordinates in no interlace mode
+ INoeud endNode = points.end();
+ for ( int j = 0; j < SPACE_DIM; ++j ) {
+ for ( INoeud in = firstNode; in != endNode; ++in ) {
+ _noeud & node = in->second;
+ node.coord[ j ] = geoFile.getReal();
+ }
+ }
+ }
+ else if ( word == "part" )
+ {
+ // Another part encounters
+ // -----------------------
+ partNum = geoFile.getInt();
+ nbParts++;
+ geoFile.toNextLine();
+
+ string partName = geoFile.getLine();
+ if ( partName.empty() )
+ partName = "Part_" + restLine;
+
+ if ( imed.groupes.capacity() - imed.groupes.size() < theMaxNbTypes )
+ imed.groupes.reserve( size_t( 1.5 * imed.groupes.size() ));
+ imed.groupes.push_back(_groupe());
+ partGroupe = & imed.groupes.back();
+ partGroupe->nom = partName;
+ partGroupe->groupes.reserve( theMaxNbTypes );
+ }
+ else if ( word == "block" )
+ {
+ // Structured type
+ // ------------------
+ bool rectilinear = ( restLine.find( "rectilinear" ) != restLine.npos );
+ bool uniform = ( restLine.find( "uniform" ) != restLine.npos );
+ bool curvilinear = ( !rectilinear && !uniform );
+ bool iblanked = ( restLine.find( "iblanked" ) != restLine.npos );
+ bool with_ghost = ( restLine.find( "with_ghost" ) != restLine.npos );
+ bool range = ( restLine.find( "range" ) != restLine.npos );
+
+ // dimension
+ int I = geoFile.getInt();
+ int J = geoFile.getInt();
+ int K = geoFile.getInt();
+ int NumberOfNodes = I*J*K;
+ if ( !NumberOfNodes ) continue;
+
+ // range
+ if ( range ) {
+ vector<int> ijkRange; // imin imax jmin jmax kmin kmax
+ ijkRange.reserve(6);
+ while ( ijkRange.size() < 6 )
+ ijkRange.push_back( geoFile.getInt() );
+ I = ijkRange[1]-ijkRange[0]+1;
+ J = ijkRange[3]-ijkRange[2]+1;
+ K = ijkRange[5]-ijkRange[4]+1;
+ NumberOfNodes = I*J*K;
+ }
+ // add nodes
+ int nodeShift = points.empty() ? 0 : points.rbegin()->first;
+ for ( int iNode = 1; iNode <= NumberOfNodes; ++iNode ) {
+ INoeud inoeud = points.insert( points.end(), make_pair( iNode + nodeShift, _noeud()));
+ _noeud & node = inoeud->second;
+ node.number = inoeud->first;
+ node.coord.resize( SPACE_DIM );
+ }
+ INoeud firstNode = points.find( nodeShift + 1 );
+ INoeud endNode = points.end();
+
+ GRID grid; // calculator of unstructured data
+ grid._iArrayLength = I;
+ grid._jArrayLength = J;
+ grid._kArrayLength = K;
+ grid._numberOfNodes = NumberOfNodes ;
+ grid._spaceDimension = SPACE_DIM;
+ if ( J < 2 ) { grid._spaceDimension--; grid._jArrayLength = 0; }
+ if ( K < 2 ) { grid._spaceDimension--; grid._kArrayLength = 0; }
+ int nbElems = grid.getNumberOfElements(MED_CELL, MED_ALL_ELEMENTS);
+
+ if ( curvilinear ) // read coordinates for all nodes
+ {
+ for ( int j = 0; j < SPACE_DIM; ++j ) {
+ for ( INoeud in = firstNode; in != endNode; ++in )
+ in->second.coord[ j ] = geoFile.getReal();
+ }
+ grid._gridType = MED_BODY_FITTED;
+ }
+ else if ( rectilinear ) // read delta vectors with non-regular spacing
+ {
+ grid._iArray = (double*)geoFile.convertReals<double>( I );
+ grid._jArray = (double*)geoFile.convertReals<double>( J );
+ grid._kArray = (double*)geoFile.convertReals<double>( K );
+ grid._gridType = MED_CARTESIAN;
+ }
+ else // uniform: read grid origine and delta vectors for regular spacing grid
+ {
+ TDblOwner xyzOrigin( (double*)geoFile.convertReals<double>( 3 ));
+ TDblOwner xyzDelta ( (double*)geoFile.convertReals<double>( 3 ));
+ // compute full delta vectors
+ grid._iArray = new double[ I ];
+ grid._jArray = new double[ J ];
+ grid._kArray = new double[ K ];
+ double* coors[SPACE_DIM] = { grid._iArray, grid._jArray, grid._kArray };
+ int size [SPACE_DIM] = { I, J, K };
+ for ( int j = 0; j < SPACE_DIM; ++j ) {
+ double* coo = coors[ j ];
+ double* cooEnd = coo + size[ j ];
+ coo[0] = xyzOrigin[ j ];
+ while ( ++coo < cooEnd )
+ *coo = coo[-1] + xyzDelta[ j ];
+ }
+ grid._gridType = MED_CARTESIAN;
+ }
+
+ // iblanks
+ if ( iblanked )
+ geoFile.skip( NumberOfNodes, /*nbPerLine =*/ 1, INT_WIDTH_GOLD);
+ // ghosts
+ if ( with_ghost ) {
+ geoFile.getWord(); // "ghost_flags"
+ geoFile.skip( nbElems, /*nbPerLine =*/ 1, INT_WIDTH_GOLD);
+ }
+ // node ids
+ if ( haveNodeIds && geoFile.lookAt( "node_ids" )) {
+ geoFile.getWord(); // "node_ids"
+ geoFile.skip( NumberOfNodes, /*nbPerLine =*/ 1, INT_WIDTH_GOLD);
+ }
+ // element ids
+ if ( haveElemIds && geoFile.lookAt( "element_ids" ) ) {
+ geoFile.getWord(); // "element_ids"
+ geoFile.skip( nbElems, /*nbPerLine =*/ 1, INT_WIDTH_GOLD);
+ }
+
+ if ( !curvilinear ) // let GRID compute all coordinates
+ {
+ grid._coordinate = new COORDINATE;
+ const double * coo = grid.getCoordinates(MED_FULL_INTERLACE);
+ typedef _ValueIterator< double > TCoordIt;
+ TCoordIt xCoo( coo+0, grid._spaceDimension);
+ TCoordIt yCoo( coo+1, grid._spaceDimension);
+ TCoordIt zCoo( coo+2, grid._spaceDimension);
+ if ( grid._spaceDimension < 3 ) zCoo = TCoordIt( grid._kArray, 0 );
+ if ( grid._spaceDimension < 2 ) yCoo = TCoordIt( grid._jArray, 0 );
+ for ( INoeud in = firstNode; in != endNode; ++in ) {
+ _noeud& node = in->second;
+ node.coord[ 0 ] = xCoo.next();
+ node.coord[ 1 ] = yCoo.next();
+ node.coord[ 2 ] = zCoo.next();
+ }
+ }
+
+ // let GRID calculate connectivity
+
+ const int * conn = grid.getConnectivity( MED_FULL_INTERLACE, MED_NODAL,
+ MED_CELL, MED_ALL_ELEMENTS );
+ medGeometryElement elemType = grid.getElementType( MED_CELL, 1 );
+ int nbElemNodes = elemType % 100;
+
+ partGroupe->mailles.resize( nbElems );
+ _maille ma( elemType, nbElemNodes );
+ ma.sommets.resize( nbElemNodes );
+
+ for ( int i = 0, nIndex = 0; i < nbElems; ++i ) {
+ for ( int n = 0; n < nbElemNodes; ++n ) {
+ int nodeID = conn[ nIndex++ ];
+ ma.sommets[n] = points.find( nodeID + nodeShift );
+ }
+ //ma.ordre = ++order;
+ partGroupe->mailles[i] = imed.insert(ma);
+ }
+
+ _SubPart subPart( partNum, "block" );
+ subPart.myNbCells = nbElems;
+ subPart.myNbNodes = NumberOfNodes;
+ subPart.myFirstNode = firstNode;
+ subPart.myCellGroupIndex = imed.groupes.size();
+ imed.addSubPart( subPart );
+ }
+ else
+ {
+ throw MEDEXCEPTION
+ ( LOCALIZED( STRING(LOC) << "Unexpected word: " << word <<
+ " in " << getDataFileName()));
+ }
+ } // while ( !geoFile.eof() )
+
+ if ( nbParts > 1 )
+ imed.mergeNodesAndElements(TOLERANCE);
+
+ END_OF_MED(LOC);
+}
+
+//================================================================================
+/*!
+ * \brief Read mesh in Gold Binary format
+ */
+//================================================================================
+
+void ENSIGHT_MESH_RDONLY_DRIVER::readGoldBinary(_InterMed & imed)
+{
+ const char * LOC ="ENSIGHT_MESH_RDONLY_DRIVER::readGoldBinary() : ";
+ BEGIN_OF_MED(LOC);
+
+ _BinaryFileReader geoFile( getDataFileName() );
+
+ // check if swapping bytes needed
+ try {
+ countPartsBinary( geoFile, isSingleFileMode());
+ }
+ catch (...) {
+ geoFile.swapBytes();
+ geoFile.rewind();
+ }
+ if ( getIndexInDataFile() <= 1 )
+ geoFile.rewind();
+ if ( geoFile.getPosition() == 0 ) {
+ TStrOwner format( geoFile.getLine() ); // "C|Fortran Binary"
+ if ( !contains( "C Binary", format )) {
+ if ( contains( "Fortran Binary", format ))
+ throw(MEDEXCEPTION(STRING(LOC) << "Fortran Binary format not supported"));
+ else
+ throw(MEDEXCEPTION(STRING(LOC) << "unexpected line in " << getDataFileName()
+ << "\n" << format.myValues));
+ }
+ }
+ if ( isSingleFileMode() ) {
+ // one time step may be skipped by countPartsBinary
+ int curTimeStep = geoFile.getPosition() ? 2 : 1 ;
+ while ( curTimeStep < getIndexInDataFile() ) {
+ countPartsBinary( geoFile, true ); // skip time step
+ curTimeStep++;
+ }
+ while (1) {
+ TStrOwner line( geoFile.getLine() );
+ if ( isTimeStepBeginning( line.myValues ))
+ break;
+ }
+ }
+ // ----------------------
+ // Read mesh description
+ // ----------------------
+ {
+ TStrOwner descriptionLine1 ( geoFile.getLine() );
+ TStrOwner descriptionLine2 ( geoFile.getLine() );
+
+ // find out if the file was created by MED driver
+ _isMadeByMed = contains( theDescriptionPrefix, descriptionLine1 );
+
+ if ( _isMadeByMed )
+ _ptrMesh->setDescription( descriptionLine2.myValues );
+ else
+ _ptrMesh->setDescription( string(descriptionLine1) + descriptionLine2.myValues );
+ }
+
+ // ----------------------------------------
+ // Find out presence of node/elem numbers
+ // ----------------------------------------
+
+ // EnSight User Manual (for v8) says:
+// You do not have to assign node IDs. If you do, the element connectivities are
+// based on the node numbers. If you let EnSight assign the node IDs, the nodes
+// are considered to be sequential starting at node 1, and element connectivity is
+// done accordingly. If node IDs are set to off, they are numbered internally;
+// however, you will not be able to display or query on them. If you have node
+// IDs in your data, you can have EnSight ignore them by specifying "node id
+// ignore." Using this option may reduce some of the memory taken up by the
+// Client and Server, but display and query on the nodes will not be available.
+
+ // read "node|element id <off|given|assign|ignore>"
+ bool haveNodeIds, haveElemIds;
+ {
+ TStrOwner nodeIds( geoFile.getLine() );
+ TStrOwner elemIds( geoFile.getLine() );
+
+ haveNodeIds = ( contains( "given", nodeIds ) || contains( "ignore", nodeIds ) );
+ haveElemIds = ( contains( "given", elemIds ) || contains( "ignore", elemIds ) );
+ }
+
+ typedef map<int,_noeud>::iterator INoeud;
+ map<int,_noeud> & points = imed.points;
+ INoeud firstNode;
+
+ _groupe * partGroupe = 0;
+ int partNum = 0, nbParts = 0;
+
+ TFltOwner extents(0); // extents: xmin xmax ymin ymax zmin zmax
+
+ while ( !geoFile.eof() )
+ {
+ TStrOwner line( geoFile.getLine() );
+ if ( isSingleFileMode() && isTimeStepEnd( line.myValues ))
+ break;
+ string word, restLine;
+ TStrTool::split( line.myValues, word, restLine );
+
+ const TEnSightElemType & partType = getEnSightType( word );
+ if ( partType._medType != MED_ALL_ELEMENTS )
+ {
+ // Unstructured element type encounters
+ // --------------------------------------
+ int nbElems = *TIntOwner( geoFile.getInt(1) ); // ne
+ int nbElemNodes = partType._medType % 100;
+ bool isGhost = isGhostType( word );
+ int nodeShift = points.empty() ? 0 : points.rbegin()->first;
+
+ // read element ids
+ TIntOwner elemIds( haveElemIds ? geoFile.getInt( nbElems ): 0 ); // id_e
+
+ if ( isGhost ) { // do not store ghost elements (?)
+ int nbNodes = nbElems * nbElemNodes;
+ if ( partType._name == "nsided" ) // polygons
+ {
+ TIntOwner nbNodesInFace( geoFile.getInt( nbElems ));
+ nbNodes = std::accumulate( nbNodesInFace.myValues, nbNodesInFace.myValues + nbElems, 0 );
+ }
+ else if ( partType._name == "nfaced" ) // polyhedrons
+ {
+ TIntOwner nbElemFaces( geoFile.getInt( nbElems ));
+ int nbFaces = accumulate( nbElemFaces.myValues, nbElemFaces.myValues + nbElems, 0 );
+ TIntOwner nbNodesInFace( geoFile.getInt( nbFaces ));
+ nbNodes = std::accumulate( nbNodesInFace.myValues, nbNodesInFace.myValues + nbElems, 0 );
+ }
+ geoFile.skip( nbNodes * sizeof(int) );
+ continue;
+ }
+
+ // add a group corresponding to subPart (geoType)
+ imed.groupes.push_back(_groupe());
+ _groupe & groupe = imed.groupes.back();
+ groupe.mailles.resize( nbElems );
+
+ // find out if "coordinates" has already been encountered
+ _SubPartDesc coordDesc( partNum , "coordinates");
+ map< _SubPartDesc, _SubPart >::iterator descPart =
+ imed._subPartDescribed.find( coordDesc );
+ bool haveCoords = ( descPart != imed._subPartDescribed.end() );
+ if ( haveCoords ) {
+ firstNode = descPart->second.myFirstNode;
+ nodeShift -= descPart->second.myNbNodes;
+ }
+
+ // read poly element data
+ bool isPoly = ( !nbElemNodes );
+ int nbNodes = 0;
+ TIntOwner nbElemNodesVec(0), nbElemFaces(0), nbFaceNodes(0);
+ if ( partType._name == "nsided" ) // polygons
+ {
+ nbElemNodesVec.myValues = geoFile.getInt( nbElems ); // npi
+ nbNodes = accumulate( nbElemNodesVec.myValues, nbElemNodesVec.myValues + nbElems, 0 );
+ }
+ else if ( partType._name == "nfaced" ) // polyhedrons
+ {
+ nbElemFaces.myValues = geoFile.getInt( nbElems ); // nf_ei
+ int totalNbFaces = accumulate( nbElemFaces.myValues, nbElemFaces.myValues + nbElems, 0 );
+
+ nbFaceNodes.myValues = geoFile.getInt( totalNbFaces ); // np(f_ei)
+ // calculate nb of nodes in each polyhedron
+ int* nbEN = nbElemNodesVec.myValues = new int[ nbElems ];
+ const int *nbFN = nbFaceNodes, *nbEF = nbElemFaces, *nbEND = nbEN + nbElems;
+ for ( ; nbEN < nbEND; ++nbEN, ++nbEF ) {
+ nbNodes += *nbEN = accumulate( nbFN, nbFN + *nbEF, 0 );
+ nbFN += *nbEF;
+ }
+ }
+ else // standard types
+ {
+ nbElemNodesVec.myValues = new int[ 1 ];
+ nbElemNodesVec[ 0 ] = nbElemNodes;
+ nbNodes = nbElems * nbElemNodes;
+ }
+ // iterator returning nbElemNodes for standard elems and
+ // next value from nbElemNodesVec for poly elements
+ _ValueIterator<int> nbElemNodesIt( nbElemNodesVec, isPoly ? 1 : 0);
+
+ // iterator returning values form partType._medIndex for standard elems
+ // and node index (n) for poly elements
+ int n;
+ _ValueIterator<int> medIndexIt( isPoly ? &n : &partType._medIndex[0],
+ isPoly ? 0 : 1);
+ // read connectivity
+ _maille ma( partType._medType, nbElemNodes );
+ ma.sommets.resize( nbElemNodes );
+ TIntOwner connectivity( geoFile.getInt( nbNodes )); // nn_ei
+ int* nodeID = connectivity;
+ INoeud node;
+ for ( int i = 0; i < nbElems; ++i ) {
+ _ValueIterator<int> medIndex = medIndexIt;
+ nbElemNodes = nbElemNodesIt.next();
+ if ( ma.sommets.size() != nbElemNodes )
+ ma.sommets.resize( nbElemNodes );
+ for ( n = 0; n < nbElemNodes; ++n, ++nodeID ) {
+ if ( haveCoords )
+ node = points.find( *nodeID + nodeShift );
+ else
+ node = points.insert( make_pair( *nodeID + nodeShift, _noeud())).first;
+ ma.sommets[ medIndex.next() ] = node;
+ }
+ if ( haveElemIds )
+ ma.setOrdre( elemIds[ i ] );
+ groupe.mailles[i] = imed.insert(ma);
+ }
+ // store nb nodes in polyhedron faces
+ if ( nbFaceNodes.myValues ) {
+ const int* nbFaceNodesPtr = nbFaceNodes.myValues;
+ for ( int i = 0; i < nbElems; ++i ) {
+ vector<int> & nbNodesByFace = imed.polyherdalNbFaceNodes[ &groupe.maille( i ) ];
+ nbNodesByFace.assign( nbFaceNodesPtr, nbFaceNodesPtr + nbElemFaces[ i ] );
+ nbFaceNodesPtr += nbElemFaces[ i ];
+ }
+ }
+ // create subPart for "coordinates"
+ if ( !haveCoords ) {
+ _SubPart & coordSubPart = imed._subPartDescribed[ coordDesc ];
+ coordSubPart.myFirstNode = points.insert( make_pair( 1 + nodeShift, _noeud())).first;
+ }
+ // add subPart group to part group
+ int groupeIndex = imed.groupes.size();
+ partGroupe->groupes.push_back( groupeIndex );
+
+ // create subPart
+ _SubPart subPart( partNum, partType._name );
+ subPart.myNbCells = nbElems;
+ subPart.myCellGroupIndex = groupeIndex;
+ imed.addSubPart( subPart );
+ }
+ else if ( word == "coordinates" )
+ {
+ // Local node coordinates of a part
+ // ------------------------------------
+ int nbNodes = *TIntOwner( geoFile.getInt(1) ); // nn
+
+ // read node ids
+ TIntOwner nodeIds(0);
+ if ( haveNodeIds )
+ nodeIds.myValues = geoFile.getInt( nbNodes ); // id_n
+
+ // find out if "coordinates" has already been add at reading connectivity
+ _SubPartDesc coordDesc( partNum , "coordinates");
+ map< _SubPartDesc, _SubPart >::iterator descPart =
+ imed._subPartDescribed.find( coordDesc );
+ bool haveCoords = ( descPart != imed._subPartDescribed.end() );
+
+ if ( haveCoords ) {
+ // check that all nodes have been added
+ firstNode = descPart->second.myFirstNode;
+ descPart->second.myNbNodes = nbNodes;
+ INoeud inoeud = firstNode, inoEnd = points.end();
+ int id = inoeud->first, idEnd = id + nbNodes;
+ for ( ; id < idEnd; ++id ) {
+ if ( inoeud == inoEnd || inoeud->first > id ) {
+ INoeud in = points.insert( inoeud, make_pair( id, _noeud() ));
+ in->second.number = id;
+ } else {
+ ++inoeud;
+ }
+ }
+ }
+ else {
+ // add nodes
+ int nodeShift = points.empty() ? 0 : points.rbegin()->first;
+ for ( int iNode = 1; iNode <= nbNodes; ++iNode ) {
+ INoeud inoeud = points.insert( points.end(), make_pair( iNode + nodeShift, _noeud()));
+ inoeud->second.number = inoeud->first;
+ }
+ firstNode = points.find( 1 + nodeShift );
+ // create "coordinates" subPart
+ _SubPart & subPart = imed._subPartDescribed[ coordDesc ];
+ subPart.myNbNodes = nbNodes;
+ subPart.myFirstNode = firstNode;
+ }
+
+ // read coordinates in no interlace mode
+ TFltOwner noInterlaceCoords( geoFile.getFlt( nbNodes * SPACE_DIM ));
+ float* x = noInterlaceCoords;
+ float* y = x + nbNodes;
+ float* z = y + nbNodes;
+ INoeud endNode = points.end();
+ for ( INoeud in = firstNode; in != endNode; ++in ) {
+ _noeud & node = in->second;
+ node.coord.resize( SPACE_DIM );
+ node.coord[ 0 ] = *x++;
+ node.coord[ 1 ] = *y++;
+ node.coord[ 2 ] = *z++;
+ }
+ }
+ else if ( word == "part" )
+ {
+ // Another part encounters
+ // -----------------------
+ partNum = *TIntOwner( geoFile.getInt(1) );
+ nbParts++;
+
+ string partName( TStrOwner( geoFile.getLine() ));
+ if ( partName.empty() )
+ partName = "Part_" + restLine;
+
+ if ( imed.groupes.capacity() - imed.groupes.size() < theMaxNbTypes )
+ imed.groupes.reserve( size_t( 1.5 * imed.groupes.size() ));
+ imed.groupes.push_back(_groupe());
+ partGroupe = & imed.groupes.back();
+ partGroupe->nom = partName;
+ partGroupe->groupes.reserve( theMaxNbTypes );
+ }
+ else if ( word == "block" )
+ {
+ // Structured type
+ // ------------------
+ bool rectilinear = ( restLine.find( "rectilinear" ) != restLine.npos );
+ bool uniform = ( restLine.find( "uniform" ) != restLine.npos );
+ bool curvilinear = ( !rectilinear && !uniform );
+ bool iblanked = ( restLine.find( "iblanked" ) != restLine.npos );
+ bool with_ghost = ( restLine.find( "with_ghost" ) != restLine.npos );
+ bool range = ( restLine.find( "range" ) != restLine.npos );
+
+ // dimension
+ TIntOwner ijk( geoFile.getInt(3) );
+ int I = ijk[0];
+ int J = ijk[1];
+ int K = ijk[2];
+ int NumberOfNodes = I*J*K;
+ if ( !NumberOfNodes ) continue;
+
+ // range
+ if ( range ) {
+ TIntOwner ijkRange( geoFile.getInt( 6 ));// imin imax jmin jmax kmin kmax
+ I = ijkRange[1]-ijkRange[0]+1;
+ J = ijkRange[3]-ijkRange[2]+1;
+ K = ijkRange[5]-ijkRange[4]+1;
+ NumberOfNodes = I*J*K;
+ }
+ // add nodes
+ int nodeShift = points.empty() ? 0 : points.rbegin()->first;
+ for ( int iNode = 1; iNode <= NumberOfNodes; ++iNode ) {
+ INoeud inoeud = points.insert( points.end(), make_pair( iNode + nodeShift, _noeud()));
+ _noeud & node = inoeud->second;
+ node.number = inoeud->first;
+ node.coord.resize( SPACE_DIM );
+ }
+ INoeud firstNode = points.find( nodeShift + 1 );
+ INoeud endNode = points.end();
+
+ GRID grid; // calculator of unstructured data
+ grid._iArrayLength = I;
+ grid._jArrayLength = J;
+ grid._kArrayLength = K;
+ grid._numberOfNodes = NumberOfNodes ;
+ grid._spaceDimension = SPACE_DIM;
+ if ( J < 2 ) { grid._spaceDimension--; grid._jArrayLength = 0; }
+ if ( K < 2 ) { grid._spaceDimension--; grid._kArrayLength = 0; }
+ int nbElems = grid.getNumberOfElements(MED_CELL, MED_ALL_ELEMENTS);
+
+ if ( curvilinear ) // read coordinates for all nodes
+ {
+ TFltOwner noInterlaceCoords( geoFile.getFlt( NumberOfNodes * SPACE_DIM ));
+ float* x = noInterlaceCoords;
+ float* y = x + NumberOfNodes;
+ float* z = y + NumberOfNodes;
+ for ( INoeud in = firstNode; in != endNode; ++in ) {
+ _noeud & node = in->second;
+ node.coord.resize( SPACE_DIM );
+ node.coord[ 0 ] = *x++;
+ node.coord[ 1 ] = *y++;
+ node.coord[ 2 ] = *z++;
+ }
+ grid._gridType = MED_BODY_FITTED;
+ }
+ else if ( rectilinear ) // read delta vectors with non-regular spacing
+ {
+ grid._iArray = (double*)geoFile.convertReals<double>( I );
+ grid._jArray = (double*)geoFile.convertReals<double>( J );
+ grid._kArray = (double*)geoFile.convertReals<double>( K );
+ grid._gridType = MED_CARTESIAN;
+ }
+ else // uniform: read grid origine and delta vectors for regular spacing grid
+ {
+ TFltOwner xyzOrigin( geoFile.getFlt( 3 ));
+ TFltOwner xyzDelta ( geoFile.getFlt( 3 ));
+ // compute full delta vectors
+ grid._iArray = new double[ I ];
+ grid._jArray = new double[ J ];
+ grid._kArray = new double[ K ];
+ double* coors[SPACE_DIM] = { grid._iArray, grid._jArray, grid._kArray };
+ int size [SPACE_DIM] = { I, J, K };
+ for ( int j = 0; j < SPACE_DIM; ++j ) {
+ double* coo = coors[ j ];
+ double* cooEnd = coo + size[ j ];
+ coo[0] = xyzOrigin[ j ];
+ while ( ++coo < cooEnd )
+ *coo = coo[-1] + xyzDelta[ j ];
+ }
+ grid._gridType = MED_CARTESIAN;
+ }
+
+ // iblanks
+ if ( iblanked )
+ geoFile.skip( NumberOfNodes * sizeof(int) );
+ // ghosts
+ if ( with_ghost ) {
+ TStrOwner( geoFile.getLine() ); // "ghost_flags"
+ geoFile.skip( nbElems * sizeof(int) );
+ }
+ // node ids
+ if ( haveNodeIds && !geoFile.eof() ) {
+ TStrOwner nextLine( geoFile.getLine() ); // "node_ids"
+ if ( contains( "node_ids", nextLine ) )
+ geoFile.skip( NumberOfNodes * sizeof(int) );
+ else
+ geoFile.skip( -MAX_LINE_LENGTH );
+ }
+ // element ids
+ TIntOwner elemIdOwner(0);
+ _ValueIterator<int> elemIds;
+ if ( haveElemIds && !geoFile.eof() ) {
+ TStrOwner nextLine( geoFile.getLine() ); // "element_ids"
+ if ( contains( "element_ids", nextLine ) ) {
+ elemIdOwner.myValues = geoFile.getInt( nbElems );
+ elemIds = _ValueIterator<int>( elemIdOwner, 1);
+ } else {
+ geoFile.skip( -MAX_LINE_LENGTH );
+ }
+ }
+
+ if ( !curvilinear ) // let GRID compute all coordinates
+ {
+ grid._coordinate = new COORDINATE;
+ const double * coo = grid.getCoordinates(MED_FULL_INTERLACE);
+ typedef _ValueIterator< double > TCoordIt;
+ TCoordIt xCoo( coo+0, grid._spaceDimension);
+ TCoordIt yCoo( coo+1, grid._spaceDimension);
+ TCoordIt zCoo( coo+2, grid._spaceDimension);
+ if ( grid._spaceDimension < 3 ) zCoo = TCoordIt( grid._kArray, 0 );
+ if ( grid._spaceDimension < 2 ) yCoo = TCoordIt( grid._jArray, 0 );
+ for ( INoeud in = firstNode; in != endNode; ++in ) {
+ _noeud& node = in->second;
+ node.coord[ 0 ] = xCoo.next();
+ node.coord[ 1 ] = yCoo.next();
+ node.coord[ 2 ] = zCoo.next();
+ }
+ }
+
+ // let GRID calculate connectivity
+
+ const int * conn = grid.getConnectivity( MED_FULL_INTERLACE, MED_NODAL,
+ MED_CELL, MED_ALL_ELEMENTS );
+ medGeometryElement elemType = grid.getElementType( MED_CELL, 1 );
+ int nbElemNodes = elemType % 100;
+
+ partGroupe->mailles.resize( nbElems );
+ _maille ma( elemType, nbElemNodes );
+ ma.sommets.resize( nbElemNodes );
+
+ for ( int i = 0, nIndex = 0; i < nbElems; ++i ) {
+ for ( int n = 0; n < nbElemNodes; ++n ) {
+ int nodeID = conn[ nIndex++ ];
+ ma.sommets[n] = points.find( nodeID + nodeShift );
+ }
+ ma.setOrdre( elemIds.next() );
+ partGroupe->mailles[i] = imed.insert(ma);
+ }
+
+ _SubPart subPart( partNum, "block" );
+ subPart.myNbCells = nbElems;
+ subPart.myNbNodes = NumberOfNodes;
+ subPart.myFirstNode = firstNode;
+ subPart.myCellGroupIndex = imed.groupes.size();
+ imed.addSubPart( subPart );
+ }
+ else if ( word == "extents" )
+ {
+ extents.myValues = geoFile.getFlt( 6 );
+ }
+ else
+ {
+ throw MED_EXCEPTION
+ ( LOCALIZED( STRING(LOC) << "Unexpected word: " << word <<
+ " in " << getDataFileName()));
+ }
+ } // while ( !geoFile.eof() )
+
+ if ( nbParts > 1 )
+ imed.mergeNodesAndElements(TOLERANCE);
+
+ END_OF_MED(LOC);
+}
+
+//================================================================================
+/*!
+ * \brief Read mesh in Ensight6 ASCII format
+ */
+//================================================================================
+
+void ENSIGHT_MESH_RDONLY_DRIVER::read6ASCII(_InterMed & imed)
+{
+ const char * LOC ="ENSIGHT_MESH_RDONLY_DRIVER::read6ASCII() : ";
+ BEGIN_OF_MED(LOC);
+
+ _ASCIIFileReader geoFile( getDataFileName() );
+
+ if ( isSingleFileMode() ) {
+ int curTimeStep = 1;
+ while ( curTimeStep < getIndexInDataFile() ) {
+ while ( !geoFile.isTimeStepEnd())
+ geoFile.getLine();
+ curTimeStep++;
+ }
+ while ( !geoFile.isTimeStepBeginning() )
+ geoFile.getLine();
+ }
+ // ----------------------
+ // Read mesh description
+ // ----------------------
+ {
+ string descriptionLine1 = geoFile.getLine();
+ string descriptionLine2 = geoFile.getLine();
+
+ // find out if the file was created by MED driver
+ int prefixSize = strlen( theDescriptionPrefix );
+ _isMadeByMed = ( descriptionLine1.substr(0, prefixSize ) == theDescriptionPrefix );
+
+ if ( _isMadeByMed )
+ descriptionLine1 = descriptionLine1.substr( prefixSize );
+ _ptrMesh->setDescription( descriptionLine1 + descriptionLine2 );
+ }
+
+ // ----------------------------------------
+ // Find out presence of node/elem numbers
+ // ----------------------------------------
+
+ // EnSight User Manual (for v8) says:
+// You do not have to assign node IDs. If you do, the element connectivities are
+// based on the node numbers. If you let EnSight assign the node IDs, the nodes
+// are considered to be sequential starting at node 1, and element connectivity is
+// done accordingly. If node IDs are set to off, they are numbered internally;
+// however, you will not be able to display or query on them. If you have node
+// IDs in your data, you can have EnSight ignore them by specifying "node id
+// ignore." Using this option may reduce some of the memory taken up by the
+// Client and Server, but display and query on the nodes will not be available.
+
+ // read "node|element id <off|given|assign|ignore>"
+ geoFile.getWord(); geoFile.getWord();
+ string nodeIds = geoFile.getWord();
+ geoFile.getWord(); geoFile.getWord();
+ string elemIds = geoFile.getWord();
+
+ bool haveNodeIds = ( nodeIds == "given" || nodeIds == "ignore" );
+ bool haveElemIds = ( elemIds == "given" || elemIds == "ignore" );
+
+ map<int,_noeud> & points = imed.points;
+ typedef map<int,_noeud>::iterator INoeud;
+
+ int haveStructuredParts = 0, haveUnstructuredParts = 0;
+
+ _groupe * partGroupe = 0;
+ int partNum = 0;
+
+ while ( !geoFile.isTimeStepEnd() )
+ {
+ string word, restLine, line = geoFile.getLine();
+ TStrTool::split( line, word, restLine );
+
+ const TEnSightElemType & partType = getEnSightType( word );
+ if ( !partType._medIndex.empty() )
+ {
+ // Unstructured element type encounters
+ // --------------------------------------
+ int nbElemNodes = partType._medType % 100;
+ int nbElems = geoFile.getInt();
+ if ( nbElems > 0 )
+ haveUnstructuredParts++;
+
+ imed.groupes.push_back(_groupe());
+ _groupe & groupe = imed.groupes.back();
+ groupe.mailles.resize( nbElems );
+
+ // read connectivity
+ _maille ma( partType._medType, nbElemNodes );
+ ma.sommets.resize( nbElemNodes );
+ INoeud node;
+ for ( int i = 0; i < nbElems; ++i ) {
+ if ( haveElemIds )
+ geoFile.getInt();
+ for ( int n = 0; n < nbElemNodes; ++n ) {
+ int nodeID = geoFile.getInt();
+ ma.sommets[ partType._medIndex[n] ] = points.find( nodeID );
+ }
+ //ma.ordre = ++order;
+ groupe.mailles[i] = imed.insert(ma);
+ }
+
+ int groupeIndex = imed.groupes.size();
+ partGroupe->groupes.push_back( groupeIndex );
+
+ _SubPart subPart( partNum, partType._name );
+ subPart.myNbCells = nbElems;
+ subPart.myCellGroupIndex = groupeIndex;
+ imed.addSubPart( subPart );
+ }
+ else if ( word == "part" )
+ {
+ // Another part encounters
+ // -----------------------
+ partNum = atoi( restLine.c_str() );
+
+ string partName = geoFile.getLine();
+ if ( partName.empty() )
+ partName = "Part_" + restLine;
+
+ if ( imed.groupes.capacity() - imed.groupes.size() < theMaxNbTypes )
+ imed.groupes.reserve( size_t( 1.5 * imed.groupes.size() ));
+ imed.groupes.push_back(_groupe());
+ partGroupe = & imed.groupes.back();
+ partGroupe->nom = partName;
+ partGroupe->groupes.reserve( theMaxNbTypes );
+ }
+ else if ( word == "block" )
+ {
+ // Structured type
+ // ------------------
+ bool iblanked = ( restLine == "iblanked" );
+
+ // dimension
+ int I = geoFile.getInt();
+ int J = geoFile.getInt();
+ int K = geoFile.getInt();
+ int NumberOfNodes = I*J*K;
+ if ( !NumberOfNodes ) continue;
+ haveStructuredParts++;
+
+ // add nodes
+ int nodeShift = points.empty() ? 0 : points.rbegin()->first;
+ for ( int iNode = 1; iNode <= NumberOfNodes; ++iNode ) {
+ INoeud inoeud = points.insert( points.end(), make_pair( iNode + nodeShift, _noeud()));
+ _noeud & node = inoeud->second;
+ node.number = inoeud->first;
+ node.coord.resize( SPACE_DIM );
+ }
+ // read coordinates
+ INoeud firstNode = points.find( nodeShift + 1 );
+ INoeud endNode = points.end();
+ for ( int j = 0; j < SPACE_DIM; ++j ) {
+ for ( INoeud in = firstNode; in != endNode; ++in ) {
+ _noeud & node = in->second;
+ node.coord[ j ] = geoFile.getReal();
+ }
+ }
+ // iblanks
+ if ( iblanked )
+ geoFile.skip(NumberOfNodes, /*nbPerLine =*/ 10, INT_WIDTH_6);
+
+ // let GRID calculate connectivity
+ GRID grid;
+ grid._numberOfNodes = NumberOfNodes ;
+ grid._iArrayLength = I;
+ grid._jArrayLength = J;
+ grid._kArrayLength = K;
+ grid._gridType = MED_BODY_FITTED;
+ grid._spaceDimension= SPACE_DIM;
+ if ( J < 2 ) { grid._spaceDimension--; grid._jArrayLength = 0; }
+ if ( K < 2 ) { grid._spaceDimension--; grid._kArrayLength = 0; }
+
+ const int * conn = grid.getConnectivity( MED_FULL_INTERLACE, MED_NODAL,
+ MED_CELL, MED_ALL_ELEMENTS );
+ medGeometryElement elemType = grid.getElementType( MED_CELL, 1 );
+ int nbElemNodes = elemType % 100;
+ int nbElems = grid.getNumberOfElements(MED_CELL, MED_ALL_ELEMENTS);
+
+ partGroupe->mailles.resize( nbElems );
+ _maille ma( elemType, nbElemNodes );
+ ma.sommets.resize( nbElemNodes );
+
+ for ( int i = 0, nIndex = 0; i < nbElems; ++i ) {
+ for ( int n = 0; n < nbElemNodes; ++n ) {
+ int nodeID = conn[ nIndex++ ];
+ ma.sommets[n] = points.find( nodeID + nodeShift );
+ }
+ //ma.ordre = ++order;
+ partGroupe->mailles[i] = imed.insert(ma);
+ }
+
+ _SubPart subPart( partNum, "block" );
+ subPart.myNbCells = nbElems;
+ subPart.myNbNodes = NumberOfNodes;
+ subPart.myFirstNode = firstNode;
+ subPart.myCellGroupIndex = imed.groupes.size();
+ imed.addSubPart( subPart );
+ }
+ else if ( word == "coordinates" )
+ {
+ // ------------------------------------
+ // Unstructured global node coordinates
+ // ------------------------------------
+ int nbNodes = geoFile.getInt();
+
+ cout << "-> loading coordinates of " << nbNodes << " nodes " << endl ;
+
+ INoeud inoeud;
+ for ( int i=0 ; i < nbNodes ; i++ )
+ {
+ if ( haveNodeIds ) {
+ int nodeID = geoFile.getInt();
+ inoeud = points.insert( make_pair( nodeID, _noeud() )).first;
+ inoeud->second.number = nodeID;
+ }
+ else {
+ int nodeID = i + 1;
+ inoeud = points.insert( points.end(), make_pair( nodeID, _noeud()));
+ inoeud->second.number = nodeID;
+ }
+ _noeud & node = inoeud->second;
+ node.coord.resize( SPACE_DIM );
+ node.coord[ 0 ] = geoFile.getReal();
+ node.coord[ 1 ] = geoFile.getReal();
+ node.coord[ 2 ] = geoFile.getReal();
+ }
+
+ _SubPartDesc cooDesc = _SubPartDesc::globalCoordDesc();
+ _SubPart subPart( cooDesc.partNumber(), cooDesc.typeName() );
+ subPart.myNbNodes = nbNodes;
+ subPart.myFirstNode = points.begin();
+ imed.addSubPart( subPart );
+ }
+ else
+ {
+ throw MED_EXCEPTION
+ ( LOCALIZED( STRING(LOC) << "Unexpected word: " << word <<
+ " in " << getDataFileName()));
+ }
+ } // while ( !geoFile.eof() )
+
+ if ( haveStructuredParts && haveUnstructuredParts || haveStructuredParts > 1 )
+ imed.mergeNodesAndElements(TOLERANCE);
+
+ END_OF_MED(LOC);
+}
+
+//================================================================================
+/*!
+ * \brief Read mesh in Ensight6 ASCII format
+ */
+//================================================================================
+
+void ENSIGHT_MESH_RDONLY_DRIVER::read6Binary(_InterMed & imed)
+{
+ const char * LOC ="ENSIGHT_MESH_RDONLY_DRIVER::read6Binary() : ";
+ BEGIN_OF_MED(LOC);
+
+ _BinaryFileReader geoFile( getDataFileName() );
+
+ // check if swapping bytes needed
+ try {
+ countPartsBinary( geoFile, isSingleFileMode());
+ }
+ catch (...) {
+ geoFile.swapBytes();
+ geoFile.rewind();
+ }
+ if ( getIndexInDataFile() <= 1 )
+ geoFile.rewind();
+ if ( geoFile.getPosition() == 0 ) {
+ TStrOwner format( geoFile.getLine() ); // "C|Fortran Binary"
+ if ( !contains( "C Binary", format )) {
+ if ( contains( "Fortran Binary", format ))
+ throw(MEDEXCEPTION(STRING(LOC) << "Fortran Binary format not supported"));
+ else
+ throw(MEDEXCEPTION(STRING(LOC) << "unexpected line in " << getDataFileName()
+ << "\n" << format.myValues));
+ }
+ }
+
+ if ( isSingleFileMode() ) {
+ // one time step may be skipped by countPartsBinary
+ int curTimeStep = geoFile.getPosition() ? 2 : 1 ;
+ while ( curTimeStep < getIndexInDataFile() ) {
+ countPartsBinary( geoFile, true ); // skip time step
+ curTimeStep++;
+ }
+ while (1) {
+ TStrOwner line( geoFile.getLine() );
+ if ( isTimeStepBeginning( line.myValues ))
+ break;
+ }
+ }
+ // ----------------------
+ // Read mesh description
+ // ----------------------
+ {
+ TStrOwner descriptionLine1( geoFile.getLine() );
+ TStrOwner descriptionLine2( geoFile.getLine() );
+
+ // find out if the file was created by MED driver
+ _isMadeByMed = contains( theDescriptionPrefix, descriptionLine1 );
+
+ if ( _isMadeByMed )
+ _ptrMesh->setDescription( descriptionLine2.myValues );
+ else
+ _ptrMesh->setDescription( string(descriptionLine1) + descriptionLine2.myValues );
+ }
+
+ // ----------------------------------------
+ // Find out presence of node/elem numbers
+ // ----------------------------------------
+
+ // EnSight User Manual (for v8) says:
+// You do not have to assign node IDs. If you do, the element connectivities are
+// based on the node numbers. If you let EnSight assign the node IDs, the nodes
+// are considered to be sequential starting at node 1, and element connectivity is
+// done accordingly. If node IDs are set to off, they are numbered internally;
+// however, you will not be able to display or query on them. If you have node
+// IDs in your data, you can have EnSight ignore them by specifying "node id
+// ignore." Using this option may reduce some of the memory taken up by the
+// Client and Server, but display and query on the nodes will not be available.
+
+ // read "node|element id <off|given|assign|ignore>"
+ bool haveNodeIds, haveElemIds;
+ {
+ TStrOwner nodeIds( geoFile.getLine() );
+ TStrOwner elemIds( geoFile.getLine() );
+
+ haveNodeIds = ( contains( "given", nodeIds ) || contains( "ignore", nodeIds ) );
+ haveElemIds = ( contains( "given", elemIds ) || contains( "ignore", elemIds ) );
+ }
+ map<int,_noeud> & points = imed.points;
+ typedef map<int,_noeud>::iterator INoeud;
+
+ int haveStructuredParts = 0, haveUnstructuredParts = 0;
+
+ _groupe * partGroupe = 0;
+ int partNum = 0;
+
+ while ( !geoFile.eof() )
+ {
+ TStrOwner line( geoFile.getLine() );
+ if ( isSingleFileMode() && isTimeStepEnd( line.myValues ))
+ break;
+ string word, restLine;
+ TStrTool::split( line.myValues, word, restLine );
+
+ const TEnSightElemType & partType = getEnSightType( word );
+ if ( !partType._medIndex.empty() )
+ {
+ // Unstructured element type encounters
+ // --------------------------------------
+ int nbElemNodes = partType._medType % 100;
+ int nbElems = *TIntOwner( geoFile.getInt(1) ); // ne
+ if ( nbElems > 0 )
+ haveUnstructuredParts++;
+
+ TIntOwner numbers(0);
+ if ( haveElemIds )
+ numbers.myValues = geoFile.getInt( nbElems ); // id_e
+
+ imed.groupes.push_back(_groupe());
+ _groupe & groupe = imed.groupes.back();
+ groupe.mailles.resize( nbElems );
+
+ // read connectivity
+ _maille ma( partType._medType, nbElemNodes );
+ ma.sommets.resize( nbElemNodes );
+ TIntOwner connectivity( geoFile.getInt( nbElems * nbElemNodes ));
+ int* nodeID = connectivity;
+ INoeud node;
+ for ( int i = 0; i < nbElems; ++i ) {
+ for ( int n = 0; n < nbElemNodes; ++n, ++nodeID )
+ ma.sommets[ partType._medIndex[n] ] = points.find( *nodeID );
+ //ma.ordre = ++order;
+ groupe.mailles[i] = imed.insert(ma);
+ }
+
+ int groupeIndex = imed.groupes.size();
+ partGroupe->groupes.push_back( groupeIndex );
+
+ _SubPart subPart( partNum, partType._name );
+ subPart.myNbCells = nbElems;
+ subPart.myCellGroupIndex = groupeIndex;
+ imed.addSubPart( subPart );
+ }
+ else if ( word == "part" )
+ {
+ // Another part encounters
+ // -----------------------
+ partNum = atoi( restLine.c_str() );
+
+ string partName( TStrOwner( geoFile.getLine() ));
+ if ( partName.empty() )
+ partName = "Part_" + restLine;
+
+ if ( imed.groupes.capacity() - imed.groupes.size() < theMaxNbTypes )
+ imed.groupes.reserve( size_t( 1.5 * imed.groupes.size() ));
+ imed.groupes.push_back(_groupe());
+ partGroupe = & imed.groupes.back();
+ partGroupe->nom = partName;
+ partGroupe->groupes.reserve( theMaxNbTypes );
+ }
+ else if ( word == "block" )
+ {
+ // Structured type
+ // ------------------
+ bool iblanked = ( restLine == "iblanked" );
+
+ // dimension
+ TIntOwner ijk( geoFile.getInt(3) );
+ int I = ijk[0];
+ int J = ijk[1];
+ int K = ijk[2];
+ int NumberOfNodes = I*J*K;
+ if ( !NumberOfNodes ) continue;
+ haveStructuredParts++;
+
+ // read coordinates
+ int nodeShift = points.empty() ? 0 : points.rbegin()->first;
+ {
+ TFltOwner noInterlaceCoords( geoFile.getFlt( NumberOfNodes * SPACE_DIM ));
+ float* x = noInterlaceCoords;
+ float* y = x + NumberOfNodes;
+ float* z = y + NumberOfNodes;
+ for ( int iNode = 1; iNode <= NumberOfNodes; ++iNode ) {
+ INoeud inoeud = points.insert( points.end(), make_pair( iNode + nodeShift, _noeud()));
+ _noeud & node = inoeud->second;
+ node.number = inoeud->first;
+ node.coord.resize( SPACE_DIM );
+ node.coord[0] = *x++;
+ node.coord[1] = *y++;
+ node.coord[2] = *z++;
+ }
+ }
+ // iblanks
+ if ( iblanked )
+ geoFile.skip(NumberOfNodes * sizeof(int));
+
+ // let GRID calculate connectivity
+ GRID grid;
+ grid._numberOfNodes = NumberOfNodes ;
+ grid._iArrayLength = I;
+ grid._jArrayLength = J;
+ grid._kArrayLength = K;
+ grid._gridType = MED_BODY_FITTED;
+ grid._spaceDimension= SPACE_DIM;
+ if ( J < 2 ) { grid._spaceDimension--; grid._jArrayLength = 0; }
+ if ( K < 2 ) { grid._spaceDimension--; grid._kArrayLength = 0; }
+
+ const int * conn = grid.getConnectivity( MED_FULL_INTERLACE, MED_NODAL,
+ MED_CELL, MED_ALL_ELEMENTS );
+ medGeometryElement elemType = grid.getElementType( MED_CELL, 1 );
+ int nbElemNodes = elemType % 100;
+ int nbElems = grid.getNumberOfElements(MED_CELL, MED_ALL_ELEMENTS);
+
+ partGroupe->mailles.resize( nbElems );
+ _maille ma( elemType, nbElemNodes );
+ ma.sommets.resize( nbElemNodes );
+
+ for ( int i = 0, nIndex = 0; i < nbElems; ++i ) {
+ for ( int n = 0; n < nbElemNodes; ++n ) {
+ int nodeID = conn[ nIndex++ ];
+ ma.sommets[n] = points.find( nodeID + nodeShift );
+ }
+ //ma.ordre = ++order;
+ partGroupe->mailles[i] = imed.insert(ma);
+ }
+
+ _SubPart subPart( partNum, "block" );
+ subPart.myNbCells = nbElems;
+ subPart.myNbNodes = NumberOfNodes;
+ subPart.myFirstNode = points.find( nodeShift + 1 );
+ subPart.myCellGroupIndex = imed.groupes.size();
+ imed.addSubPart( subPart );
+ }
+ else if ( word == "coordinates" )
+ {
+ // ------------------------------
+ // Unstructured node coordinates
+ // ------------------------------
+ int nbNodes = *TIntOwner( geoFile.getInt(1) );
+
+ TIntOwner numbers(0);
+ if ( haveNodeIds )
+ numbers.myValues = geoFile.getInt( nbNodes );
+
+ TFltOwner fullInterlaceCoords( geoFile.getFlt( nbNodes * SPACE_DIM ));
+ float* coord = fullInterlaceCoords;
+
+ cout << "-> loading coordinates of " << nbNodes << " nodes " << endl ;
+
+ INoeud inoeud;
+ for ( int i=0 ; i < nbNodes ; i++ )
+ {
+ if ( haveNodeIds ) {
+ int nodeID = numbers[ i ];
+ inoeud = points.insert( make_pair( nodeID, _noeud() )).first;
+ inoeud->second.number = nodeID;
+ }
+ else {
+ int nodeID = i + 1;
+ inoeud = points.insert( points.end(), make_pair( nodeID, _noeud()));
+ inoeud->second.number = nodeID;
+ }
+ _noeud & node = inoeud->second;
+ node.coord.resize( SPACE_DIM );
+ node.coord[ 0 ] = *coord++;
+ node.coord[ 1 ] = *coord++;
+ node.coord[ 2 ] = *coord++;
+ }
+
+ _SubPartDesc cooDesc = _SubPartDesc::globalCoordDesc();
+ _SubPart subPart( cooDesc.partNumber(), cooDesc.typeName() );
+ subPart.myNbNodes = nbNodes;
+ subPart.myFirstNode = points.begin();
+ imed.addSubPart( subPart );
+ }
+ else
+ {
+ throw MED_EXCEPTION
+ ( LOCALIZED( STRING(LOC) << "Unexpected word: " << word <<
+ " in " << getDataFileName()));
+ }
+ } // while ( !geoFile.eof() )
+
+ if ( haveStructuredParts && haveUnstructuredParts || haveStructuredParts > 1 )
+ imed.mergeNodesAndElements(TOLERANCE);
+
+ END_OF_MED(LOC);
+}
+
+//================================================================================
+/*!
+ * \brief count number of parts in EnSight geometry file
+ */
+//================================================================================
+
+int ENSIGHT_MESH_RDONLY_DRIVER::countParts(const string& geomFileName,
+ const bool isSingleFileMode)
+{
+ const char * LOC ="ENSIGHT_MESH_RDONLY_DRIVER::countParts() : ";
+
+ int nbParts = 0;
+ if ( isBinaryDataFile( geomFileName ))
+ {
+ _BinaryFileReader geoFile(geomFileName);
+ // check if swapping bytes needed
+ try {
+ return countPartsBinary( geoFile, isSingleFileMode );
+ }
+ catch (...) {
+ }
+ geoFile.swapBytes();
+ geoFile.rewind();
+ nbParts = countPartsBinary( geoFile, isSingleFileMode );
+ }
+ else
+ {
+ _ASCIIFileReader geoFile(geomFileName);
+
+ if ( isSingleFileMode )
+ while ( !isTimeStepBeginning( geoFile.getLine() ));
+
+ geoFile.getLine(); // description line 1
+ geoFile.getLine(); // description line 2
+
+ // read "node|element id <off|given|assign|ignore>"
+ geoFile.getWord(); geoFile.getWord();
+ string nodeIds = geoFile.getWord();
+ geoFile.getWord(); geoFile.getWord();
+ string elemIds = geoFile.getWord();
+ bool haveNodeIds = ( nodeIds == "given" || nodeIds == "ignore" );
+ bool haveElemIds = ( elemIds == "given" || elemIds == "ignore" );
+
+ bool isGold = true;
+ while ( !geoFile.isTimeStepEnd() )
+ {
+ string word, restLine, line = geoFile.getLine();
+ TStrTool::split( line, word, restLine );
+
+ const TEnSightElemType & partType = getEnSightType( word );
+ if ( partType._medType != MED_ALL_ELEMENTS )
+ {
+ // Unstructured element type encounters
+ // --------------------------------------
+ int nbElemNodes = partType._medType % 100;
+ int nbElems = geoFile.getInt(); // ne
+
+ // element ids
+ if ( haveElemIds && isGold )
+ geoFile.skip( nbElems, /*nbPerLine =*/ 1, INT_WIDTH_GOLD ); // id_e
+
+ // skip connectivity
+ int nbNodes = nbElems * nbElemNodes;
+ if ( partType._name == "nsided" ) // polygons
+ {
+ for ( int i = 0; i < nbElems; ++i )
+ nbNodes += geoFile.getInt();
+ geoFile.skip( nbNodes * INT_WIDTH_GOLD, /*nbLines = */nbElems );
+ }
+ else if ( partType._name == "nfaced" ) // polyhedrons
+ {
+ int nbFaces = 0;
+ for ( int i = 0; i < nbElems; ++i )
+ nbFaces += geoFile.getInt();
+ for ( int f = 0; f < nbFaces; ++f )
+ nbNodes += geoFile.getInt();
+ geoFile.skip( nbNodes * INT_WIDTH_GOLD, /*nbLines = */nbFaces );
+ }
+ else // standard types
+ {
+ if ( isGold )
+ geoFile.skip( nbNodes, nbElemNodes, INT_WIDTH_GOLD );
+ else if ( haveElemIds )
+ geoFile.skip( nbNodes + nbElems, nbElemNodes+1, INT_WIDTH_6 );
+ else
+ geoFile.skip( nbNodes, nbElemNodes, INT_WIDTH_6 );
+ }
+ }
+ else if ( word == "coordinates" )
+ {
+ isGold = ( nbParts > 0 );
+ int nbNodes = geoFile.getInt(); // nn
+
+ if ( isGold )
+ {
+ if ( haveNodeIds )
+ geoFile.skip( nbNodes, /*nbPerLine =*/ 1, INT_WIDTH_GOLD ); // node ids
+ geoFile.skip( nbNodes * SPACE_DIM, /*nbPerLine =*/ 1, FLT_WIDTH ); // coordinates
+ }
+ else {
+ int coordWidth = 3 * FLT_WIDTH;
+ if ( haveNodeIds )
+ coordWidth += INT_WIDTH_6;
+ geoFile.skip(nbNodes, /*nbPerLine =*/ 1, coordWidth);
+ }
+ }
+ else if ( word == "part" )
+ {
+ nbParts++;
+ if ( isGold )
+ geoFile.skip( 1, /*nbPerLine =*/ 1, INT_WIDTH_GOLD ); //part number
+ else
+ geoFile.getWord(); // part number
+ geoFile.toNextLine();
+ geoFile.getLine(); // description line
+ }
+ else if ( word == "block" )
+ {
+ // Structured type
+ // ------------------
+ bool rectilinear = ( restLine.find( "rectilinear" ) != restLine.npos );
+ bool uniform = ( restLine.find( "uniform" ) != restLine.npos );
+ bool curvilinear = ( !rectilinear && !uniform );
+ bool iblanked = ( restLine.find( "iblanked" ) != restLine.npos );
+ bool with_ghost = ( restLine.find( "with_ghost" ) != restLine.npos );
+ bool range = ( restLine.find( "range" ) != restLine.npos );
+
+ // dimension
+ int I = geoFile.getInt();
+ int J = geoFile.getInt();
+ int K = geoFile.getInt();
+ int nbNodes = I*J*K;
+ if ( !nbNodes ) continue;
+
+ // range
+ if ( range ) {
+ vector<int> ijkRange; // imin imax jmin jmax kmin kmax
+ ijkRange.reserve(6);
+ while ( ijkRange.size() < 6 )
+ ijkRange.push_back( geoFile.getInt() );
+ I = ijkRange[1]-ijkRange[0]+1;
+ J = ijkRange[3]-ijkRange[2]+1;
+ K = ijkRange[5]-ijkRange[4]+1;
+ nbNodes = I*J*K;
+ }
+ int nbElems = (I-1)*(J-1)*(K-1);
+
+ if ( curvilinear ) // read coordinates for all nodes
+ {
+ if ( isGold )
+ geoFile.skip( nbNodes, /*nbPerLine =*/ 1, FLT_WIDTH );
+ else
+ geoFile.skip( nbNodes * SPACE_DIM, /*nbPerLine =*/ 6, FLT_WIDTH );
+ }
+ else if ( rectilinear ) // read delta vectors with non-regular spacing
+ {
+ geoFile.skip( I + J + K, /*nbPerLine =*/ 1, FLT_WIDTH );
+ }
+ else // uniform: read grid origine and delta vectors for regular spacing grid
+ {
+ geoFile.skip( 6, /*nbPerLine =*/ 1, FLT_WIDTH );
+ }
+
+ // iblanks
+ if ( iblanked ) {
+ if ( isGold )
+ geoFile.skip( nbNodes, /*nbPerLine =*/ 1, INT_WIDTH_GOLD);
+ else
+ geoFile.skip( nbNodes, /*nbPerLine =*/ 10, INT_WIDTH_6 );
+ }
+ // ghosts
+ if ( with_ghost ) {
+ geoFile.getWord(); // "ghost_flags"
+ geoFile.skip( nbElems, /*nbPerLine =*/ 1, INT_WIDTH_GOLD);
+ }
+ // node ids
+ if ( haveNodeIds && geoFile.lookAt( "node_ids" )) {
+ geoFile.getWord(); // "node_ids"
+ geoFile.skip( nbNodes, /*nbPerLine =*/ 1, INT_WIDTH_GOLD);
+ }
+ // element ids
+ if ( haveElemIds && geoFile.lookAt( "element_ids" ) ) {
+ geoFile.getWord(); // "element_ids"
+ geoFile.skip( nbElems, /*nbPerLine =*/ 1, INT_WIDTH_GOLD);
+ }
+ }
+ else if ( word == "extents" ) {
+ geoFile.getLine(); geoFile.getLine(); geoFile.getLine();// 3 x 2E12.5
+ }
+ else
+ {
+ throw MED_EXCEPTION ( LOCALIZED( STRING(LOC) << "Unexpected word: " << word));
+ }
+ }
+ }
+ return nbParts;
+}
+
+//================================================================================
+/*!
+ * \brief count number of parts in EnSight geometry file
+ */
+//================================================================================
+
+int ENSIGHT_MESH_RDONLY_DRIVER::countPartsBinary(_BinaryFileReader& geoFile,
+ const bool isSingleFileMode)
+{
+ const char * LOC ="ENSIGHT_MESH_RDONLY_DRIVER::countPartsBinary() : ";
+
+ if ( geoFile.getPosition() == 0 ) {
+ TStrOwner format( geoFile.getLine() ); // "C|Fortran Binary"
+ if ( !contains( "C Binary", format )) {
+ if ( contains( "Fortran Binary", format ))
+ throw(MEDEXCEPTION(STRING(LOC) << "Fortran Binary format not supported"));
+ else
+ throw(MEDEXCEPTION(STRING(LOC) << "unexpected line: \n" << format.myValues));
+ }
+ }
+
+ if ( isSingleFileMode ) {
+ while (1) {
+ TStrOwner line( geoFile.getLine() );
+ if ( isTimeStepBeginning( line.myValues ))
+ break;
+ }
+ }
+
+ // 2 description lines
+ // ----------------------
+ geoFile.skip( 2 * MAX_LINE_LENGTH );
+
+ // read "node|element id <off|given|assign|ignore>"
+ bool haveNodeIds, haveElemIds;
+ {
+ TStrOwner nodeIds( geoFile.getLine() );
+ TStrOwner elemIds( geoFile.getLine() );
+
+ haveNodeIds = ( contains( "given", nodeIds ) || contains( "ignore", nodeIds ) );
+ haveElemIds = ( contains( "given", elemIds ) || contains( "ignore", elemIds ) );
+ }
+
+ int nbParts = 0; // the result
+ bool isGold = true;
+
+ while ( !geoFile.eof() )
+ {
+ TStrOwner line( geoFile.getLine() );
+ if ( isSingleFileMode && isTimeStepEnd( line.myValues ))
+ break;
+ string word, restLine;
+ TStrTool::split( line.myValues, word, restLine );
+
+ const TEnSightElemType & partType = getEnSightType( word );
+ if ( partType._medType != MED_ALL_ELEMENTS )
+ {
+ // Unstructured element type encounters
+ // --------------------------------------
+ int nbElems = *TIntOwner( geoFile.getInt(1) ); // ne
+ int nbElemNodes = partType._medType % 100;
+
+ // read element ids
+ if ( haveElemIds )
+ geoFile.skip( nbElems * sizeof(int) ); // id_e
+
+ int nbNodes = nbElems * nbElemNodes;
+ if ( partType._name == "nsided" ) // polygons
+ {
+ TIntOwner nbNodesInFace( geoFile.getInt( nbElems ));
+ nbNodes = std::accumulate( nbNodesInFace.myValues, nbNodesInFace.myValues + nbElems, 0 );
+ }
+ else if ( partType._name == "nfaced" ) // polyhedrons
+ {
+ TIntOwner nbElemFaces( geoFile.getInt( nbElems ));
+ int nbFaces = accumulate( nbElemFaces.myValues, nbElemFaces.myValues + nbElems, 0 );
+ TIntOwner nbNodesInFace( geoFile.getInt( nbFaces ));
+ nbNodes = std::accumulate( nbNodesInFace.myValues, nbNodesInFace.myValues + nbFaces, 0 );
+ }
+ geoFile.skip( nbNodes * sizeof(int) );
+ }
+ else if ( word == "coordinates" )
+ {
+ if ( nbParts == 0 )
+ isGold = false;
+ int nbNodes = *TIntOwner( geoFile.getInt(1) ); // nn
+
+ // node ids
+ if ( haveNodeIds )
+ geoFile.skip( nbNodes * sizeof(int) ); // id_n
+
+ // coordinates
+ geoFile.skip( nbNodes * SPACE_DIM * sizeof(int) );
+ }
+ else if ( word == "part" )
+ {
+ ++nbParts;
+ if ( isGold ) geoFile.skip(sizeof(int)); // part #
+
+ geoFile.skip(MAX_LINE_LENGTH); // description line
+ }
+ else if ( word == "block" )
+ {
+ // Structured type
+ // ------------------
+ bool rectilinear = ( restLine.find( "rectilinear" ) != restLine.npos );
+ bool uniform = ( restLine.find( "uniform" ) != restLine.npos );
+ bool curvilinear = ( !rectilinear && !uniform );
+ bool iblanked = ( restLine.find( "iblanked" ) != restLine.npos );
+ bool with_ghost = ( restLine.find( "with_ghost" ) != restLine.npos );
+ bool range = ( restLine.find( "range" ) != restLine.npos );
+
+ // dimension
+ TIntOwner ijk( geoFile.getInt(3) );
+ int I = ijk[0];
+ int J = ijk[1];
+ int K = ijk[2];
+ int NumberOfNodes = I*J*K;
+ if ( !NumberOfNodes ) {
+ if ( I != 0 && J != 0 && K != 0 )
+ throw MEDEXCEPTION( "Need to swap bytes" );
+ continue;
+ }
+
+ // range
+ if ( range ) {
+ TIntOwner ijkRange( geoFile.getInt( 6 ));// imin imax jmin jmax kmin kmax
+ I = ijkRange[1]-ijkRange[0]+1;
+ J = ijkRange[3]-ijkRange[2]+1;
+ K = ijkRange[5]-ijkRange[4]+1;
+ NumberOfNodes = I*J*K;
+ }
+ int nbElems = (I-1)*(J-1)*(K-1);
+
+ if ( curvilinear ) // read coordinates for all nodes
+ {
+ geoFile.skip( NumberOfNodes * SPACE_DIM * sizeof(float) );
+ }
+ else if ( rectilinear ) // read delta vectors with non-regular spacing
+ {
+ geoFile.skip( (I+J+K) * sizeof(float) );
+ }
+ else // uniform: read grid origine and delta vectors for regular spacing grid
+ {
+ geoFile.skip( 6 * sizeof(float) );
+ }
+
+ // iblanks
+ if ( iblanked )
+ geoFile.skip( NumberOfNodes * sizeof(int) );
+ // ghosts
+ if ( with_ghost ) {
+ geoFile.skip( MAX_LINE_LENGTH ); // "ghost_flags"
+ geoFile.skip( nbElems * sizeof(int) );
+ }
+ // node ids
+ if ( haveNodeIds && isGold && !geoFile.eof() ) {
+ TStrOwner nextLine( geoFile.getLine() ); // "node_ids"
+ if ( contains( "node_ids", nextLine ) )
+ geoFile.skip( NumberOfNodes * sizeof(int) );
+ else
+ geoFile.skip( -MAX_LINE_LENGTH );
+ }
+ // element ids
+ if ( haveElemIds && isGold && !geoFile.eof() ) {
+ TStrOwner nextLine( geoFile.getLine() ); // "element_ids"
+ if ( contains( "element_ids", nextLine ) )
+ geoFile.skip( nbElems * sizeof(int) );
+ else
+ geoFile.skip( -MAX_LINE_LENGTH );
+ }
+ }
+ else if ( word == "extents" )
+ {
+ geoFile.skip( 6 * sizeof(float) );
+ }
+ else
+ {
+ throw MED_EXCEPTION ( LOCALIZED( STRING(LOC) << "Unexpected word: " << word ));
+ }
+ } // while ( !geoFile.eof() )
+
+ return nbParts;
+}
+
+GROUP* ENSIGHT_MESH_RDONLY_DRIVER::makeGroup( _groupe& grp,
+ _InterMed & imed)
+{
+ //const char* LOC = "ENSIGHT_MESH_RDONLY_DRIVER::makeGroup(): error";
+
+ // prevent creation of other groups but only this one
+ for (size_t i=0; i < imed.groupes.size(); ++i)
+ imed.groupes[i].nom.clear();
+
+ // let _intermediateMED create a GROUP from grp
+ grp.medGroup = 0; // the new GROUP should appear in grp.medGroup
+ grp.nom = "TMP";
+ vector<GROUP *> tmp;
+ imed.getGroups( tmp, tmp, tmp, tmp, imed._medMesh );
+ if ( !grp.medGroup )
+ throw MEDEXCEPTION(LOCALIZED("Can't create a GROUP from _groupe"));
+
+ grp.medGroup->setName(""); // to let a caller set a proper name
+ grp.nom = "";
+
+ // find pre-existing equal _groupe
+ _groupe * equalGroupe = 0;
+ for (unsigned int i=0; i < imed.groupes.size() && !equalGroupe; ++i) {
+ _groupe& g = imed.groupes[i];
+ if ( &g != &grp && g.medGroup && g.medGroup->deepCompare( *grp.medGroup ))
+ equalGroupe = & g;
+ }
+ if ( equalGroupe ) {
+ delete grp.medGroup;
+ grp.medGroup = equalGroupe->medGroup;
+ }
+ else { // new unique group
+
+ if ( grp.medGroup->isOnAllElements() ) // on all elems
+ grp.medGroup->setName( string("SupportOnAll_")+entNames[ grp.medGroup->getEntity() ] );
+
+ // add a new group to mesh
+ if ( !imed._isOwnMedMesh ) {
+ vector<GROUP*> * groups = 0;
+ switch ( grp.medGroup->getEntity() ) {
+ case MED_CELL: groups = & imed._medMesh->_groupCell; break;
+ case MED_FACE: groups = & imed._medMesh->_groupFace; break;
+ case MED_EDGE: groups = & imed._medMesh->_groupEdge; break;
+ case MED_NODE: groups = & imed._medMesh->_groupNode; break;
+ default:;
+ }
+ if ( groups ) {
+ groups->resize( groups->size() + 1 );
+ groups->at( groups->size() - 1) = grp.medGroup;
+ }
+ }
+ }
+ return grp.medGroup;
+}
