1 // Copyright (C) 2007-2012 CEA/DEN, EDF R&D, OPEN CASCADE
3 // Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
6 // This library is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 2.1 of the License.
11 // This library is distributed in the hope that it will be useful,
12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 // Lesser General Public License for more details.
16 // You should have received a copy of the GNU Lesser General Public
17 // License along with this library; if not, write to the Free Software
18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
23 #include "MEDMEM_EnsightMeshDriver.hxx"
29 #include "MEDMEM_define.hxx"
30 #include "MEDMEM_Family.hxx"
31 #include "MEDMEM_Group.hxx"
32 #include "MEDMEM_Coordinate.hxx"
33 #include "MEDMEM_Connectivity.hxx"
34 #include "MEDMEM_Mesh.hxx"
35 #include "MEDMEM_CellModel.hxx"
36 #include "MEDMEM_Grid.hxx"
38 #include "MEDMEM_MedMeshDriver.hxx"
41 using namespace MEDMEM;
42 using namespace MED_EN;
43 using namespace MEDMEM_ENSIGHT;
45 #define TStrTool _ASCIIFileReader
46 #define TOLERANCE 1e-15
48 //#define ELEMENT_ID_GIVEN
52 // ---------------------------------------------------------------------
54 * \brief The beginning of mesh description used to distinguish files
55 * generated by ENSIGHT_MESH_WRONLY_DRIVER from others
57 const char* theDescriptionPrefix = "Meshing from MedMemory. ";
59 // ---------------------------------------------------------------------
61 * \brief Default name of a mesh read from EnSight
63 const char* theDefaultMeshName = "EnsightMesh";
65 // ---------------------------------------------------------------------
67 * \brief Max number of types in EnSight part
69 const int theMaxNbTypes = 20;
71 // ---------------------------------------------------------------------
73 * \brief Make array with elements == index[i+1]-index[i]
74 * \param size - result array size
76 int* getNumbersByIndex( const int* index, int size, const int* elemNumbers=0)
78 int* numbers = new int[size];
81 const int *elem = elemNumbers-1, *elemEnd = elemNumbers + size;
82 while ( ++elem < elemEnd )
83 *n++ = index[*elem] - index[*elem-1];
86 const int *ind = index, *indEnd = index + size + 1;
87 while ( ++ind < indEnd )
88 *n++ = ind[0] - ind[-1];
92 // ---------------------------------------------------------------------
94 * \brief Type used to delete numbers returned by getNumbersByIndex()
96 typedef _ValueOwner<int> TNumbers;
101 ENSIGHT_MESH_DRIVER::ENSIGHT_MESH_DRIVER(): _CaseFileDriver_User(), _ptrMesh((MESH *)NULL)
105 ENSIGHT_MESH_DRIVER::ENSIGHT_MESH_DRIVER(const string & fileName,
107 :_CaseFileDriver_User(fileName,MED_EN::RDWR), _ptrMesh(ptrMesh)
110 throw MEDEXCEPTION("ENSIGHT_MESH_DRIVER(fileName, ptrMesh) : mesh is NULL");
111 _meshName = _ptrMesh->getName();
114 ENSIGHT_MESH_DRIVER::ENSIGHT_MESH_DRIVER(const string & fileName,
116 med_mode_acces accessMode)
117 :_CaseFileDriver_User(fileName,accessMode), _ptrMesh(ptrMesh)
120 throw MEDEXCEPTION("ENSIGHT_MESH_DRIVER(fileName, ptrMesh) : mesh is NULL");
121 _meshName = _ptrMesh->getName();
124 ENSIGHT_MESH_DRIVER::ENSIGHT_MESH_DRIVER(const ENSIGHT_MESH_DRIVER & driver)
125 :_CaseFileDriver_User(driver), _ptrMesh(driver._ptrMesh), _meshName(driver._meshName)
129 ENSIGHT_MESH_DRIVER::~ENSIGHT_MESH_DRIVER()
131 MESSAGE_MED("ENSIGHT_MESH_DRIVER::~ENSIGHT_MESH_DRIVER() has been destroyed");
134 void ENSIGHT_MESH_DRIVER::setMeshName(const string & meshName) { _meshName = meshName; }
136 string ENSIGHT_MESH_DRIVER::getMeshName() const { return _meshName; }
138 void ENSIGHT_MESH_DRIVER::openConst(bool checkDataFile) const
140 const char * LOC ="ENSIGHT_MESH_DRIVER::open() : ";
145 if ( getDataFileName().empty() )
147 ( LOCALIZED( STRING(LOC) << "Internal error, geometry file name is empty"));
149 if (!canOpenFile( getDataFileName(), getAccessMode() ))
151 ( LOCALIZED( STRING(LOC) << "Can not open Ensight Geometry file " << getDataFileName()
152 << " in access mode " << getAccessMode()));
156 if ( getCaseFileName().empty() )
158 ( LOCALIZED( STRING(LOC) << "Case file name is empty, "
159 "please set a correct fileName before calling open()"));
161 if ( !canOpenFile( getCaseFileName(), getAccessMode() ))
163 ( LOCALIZED( STRING(LOC) << "Can not open Ensight Case file " << getCaseFileName()
164 << " in access mode " << getAccessMode()));
170 void ENSIGHT_MESH_DRIVER::open() {
174 void ENSIGHT_MESH_DRIVER::close() {
177 // ================================================================================
179 // ================================================================================
181 ENSIGHT_MESH_WRONLY_DRIVER::ENSIGHT_MESH_WRONLY_DRIVER() :
182 ENSIGHT_MESH_DRIVER(), _append(0)
186 ENSIGHT_MESH_WRONLY_DRIVER::ENSIGHT_MESH_WRONLY_DRIVER(const string & fileName,
187 const GMESH * ptrMesh,
189 : ENSIGHT_MESH_DRIVER( fileName, (GMESH*)ptrMesh, WRONLY ), _append(append)
193 ENSIGHT_MESH_WRONLY_DRIVER::ENSIGHT_MESH_WRONLY_DRIVER(const ENSIGHT_MESH_WRONLY_DRIVER & driver)
194 : ENSIGHT_MESH_DRIVER(driver),_append(driver._append)
198 ENSIGHT_MESH_WRONLY_DRIVER::~ENSIGHT_MESH_WRONLY_DRIVER()
202 GENDRIVER * ENSIGHT_MESH_WRONLY_DRIVER::copy() const
204 return new ENSIGHT_MESH_WRONLY_DRIVER(*this) ;
207 void ENSIGHT_MESH_WRONLY_DRIVER::read() throw (MEDEXCEPTION) {
208 throw MEDEXCEPTION("ENSIGHT_MESH_WRONLY_DRIVER::read : Can't read with a WRONLY driver !");
211 //================================================================================
215 //================================================================================
217 void ENSIGHT_MESH_WRONLY_DRIVER::write() const throw (MEDEXCEPTION)
219 const char * LOC = "ENSIGHT_MESH_WRONLY_DRIVER::write() : ";
222 openConst(false) ; // check if can write to case file
224 // Ensight case organization requires a main file (filename.case) which defines organization
226 _CaseFileDriver caseFile( getCaseFileName(), this);
229 caseFile.addMesh( this );
230 caseFile.write(); // all fields of _CaseFileDriver_User are set by this method
232 openConst(true) ; // check if can write to data file
234 cout << "-> creating the Ensight geometry file " << getDataFileName() << endl ;
236 // Store mesh description and a special mark in the first two description lines each
237 // of 79 chars length maximum, while MED mesh description is up to 200 chars
238 const char* line1 = theDescriptionPrefix;
239 string line2 = _ptrMesh->getDescription();
240 for ( unsigned i = 0; i < line2.size(); ++i ) { // protect from gabage
241 if ( !line2[ i ] || !isascii( line2[ i ])) {
246 if ((int) line2.size() >= MAX_LINE_LENGTH )
247 line2.resize( MAX_LINE_LENGTH );
249 // EnSight will assign node/element visible numbers it-self
250 const char* line3 = "node id assign";
251 #ifdef ELEMENT_ID_GIVEN
252 const char* line4 = "element id given";
254 const char* line4 = "element id assign";
257 MESH* mesh = const_cast<MESH*>( _ptrMesh->convertInMESH() ) ; // we write unstructured only
259 if ( isBinaryEnSightFormatForWriting() )
261 // ======================================================
263 // ======================================================
265 _BinaryFileWriter ensightGeomFile( getDataFileName() );
267 ensightGeomFile.addString("C Binary");
268 ensightGeomFile.addString(line1);
269 ensightGeomFile.addString(line2);
270 ensightGeomFile.addString(line3);
271 ensightGeomFile.addString(line4);
273 // function to write a support as a part
274 typedef void (ENSIGHT_MESH_WRONLY_DRIVER::* TWritePart) (_BinaryFileWriter&, const SUPPORT*) const;
275 TWritePart writePart;
276 if ( isGoldFormat() )
279 writePart = & ENSIGHT_MESH_WRONLY_DRIVER::writePartGoldBinary;
283 // ENSIGHT 6. Write addionally global nodes
284 writePart = & ENSIGHT_MESH_WRONLY_DRIVER::writePart6Binary;
286 // All point are in 3D, so if we are in 1D or 2D, we complete by zero !
287 int SpaceDimension = mesh->getSpaceDimension() ;
288 int NumberOfNodes = mesh->getNumberOfNodes() ;
289 ensightGeomFile.addString("coordinates");
290 ensightGeomFile.addInt( NumberOfNodes );
291 const double *coordinate = mesh->getCoordinates(MED_FULL_INTERLACE) ;
292 if ( SpaceDimension == 3 ) {
293 ensightGeomFile.addReal(coordinate, NumberOfNodes * SpaceDimension );
296 typedef _ValueIterator< double > TComponentIt;
297 vector< TComponentIt > coordCompIt( 3 );
298 for (int j=0; j<3; j++, coordinate++)
299 if ( j < SpaceDimension )
300 coordCompIt[ j ] = TComponentIt( coordinate, SpaceDimension );
302 coordCompIt[ j ] = TComponentIt(); // to iterate on zeros
303 ensightGeomFile.addReal( coordCompIt, NumberOfNodes, MED_FULL_INTERLACE );
307 // We put connectivity
309 if ( isToWriteEntity( MED_CELL, mesh ))
311 SUPPORT *allCells = const_cast<SUPPORT*>( _ptrMesh->getSupportOnAll( MED_CELL ));
312 string oldName = allCells->getName();
313 allCells->setName( _ptrMesh->getName() );
316 (this->*writePart)( ensightGeomFile, allCells );
318 catch (MED_EXCEPTION& ex)
320 allCells->setName( oldName );
323 allCells->setName( oldName );
325 // And meshdim-1 connectivity
326 if ( isToWriteEntity( MED_FACE, mesh ))
328 const SUPPORT *allFaces = _ptrMesh->getSupportOnAll( MED_FACE );
329 (this->*writePart)( ensightGeomFile, allFaces);
331 else if ( isToWriteEntity(MED_EDGE, mesh))
333 const SUPPORT *allEdges = _ptrMesh->getSupportOnAll( MED_EDGE );
334 (this->*writePart)( ensightGeomFile, allEdges);
337 // Write all groups as parts
339 for ( int ent = MED_CELL; ent < MED_ALL_ENTITIES; ++ent )
341 medEntityMesh entity = (medEntityMesh) ent;
342 int nbGroups = mesh->getNumberOfGroups(entity);
343 for ( int i=1; i<=nbGroups; i++)
345 const GROUP* group = mesh->getGroup( entity, i );
346 (this->*writePart)( ensightGeomFile, group );
353 // ======================================================
355 // ======================================================
356 ofstream ensightGeomFile( getDataFileName().c_str(), ios::out);
357 ensightGeomFile.setf(ios::scientific);
358 ensightGeomFile.precision(5);
360 ensightGeomFile << line1 << endl
365 // function to write a support as a part
366 typedef void (ENSIGHT_MESH_WRONLY_DRIVER::* TWritePart) (ofstream&, const SUPPORT*) const;
367 TWritePart writePart;
368 if ( isGoldFormat() )
371 writePart = & ENSIGHT_MESH_WRONLY_DRIVER::writePartGoldASCII;
375 // ENSIGHT 6. Write addionally global nodes
376 writePart = & ENSIGHT_MESH_WRONLY_DRIVER::writePart6ASCII;
378 // Put points (all point are in 3D, so if we are in 1D or 2D, we complete by zero !
379 int SpaceDimension = mesh->getSpaceDimension() ;
380 int NumberOfNodes = mesh->getNumberOfNodes() ;
382 if (SpaceDimension==2) zeros = " 0.00000e+00";
383 if (SpaceDimension==1) zeros = " 0.00000e+00 0.00000e+00";
384 ensightGeomFile << "coordinates" << endl
385 << setw(8) << NumberOfNodes << endl ;
386 const double *coordinate = mesh->getCoordinates(MED_FULL_INTERLACE) ;
387 for (int i=0; i<NumberOfNodes; i++)
389 //ensightGeomFile << setw(8) << i+1 ; // node id
390 for (int j=0; j<SpaceDimension; j++, coordinate++)
391 ensightGeomFile << setw(12) << *coordinate;
392 ensightGeomFile << zeros << endl ;
396 // We put connectivity
398 if ( isToWriteEntity( MED_CELL, mesh ))
400 SUPPORT *allCells = const_cast<SUPPORT*>( _ptrMesh->getSupportOnAll( MED_CELL ));
401 string oldName = allCells->getName();
402 allCells->setName( _ptrMesh->getName() );
405 (this->*writePart)( ensightGeomFile, allCells );
407 catch (MED_EXCEPTION& ex)
409 allCells->setName( oldName );
412 allCells->setName( oldName );
414 // And meshdim-1 connectivity
415 if ( isToWriteEntity( MED_FACE, mesh ))
417 const SUPPORT *allFaces = _ptrMesh->getSupportOnAll( MED_FACE );
418 (this->*writePart)( ensightGeomFile, allFaces);
420 else if ( isToWriteEntity(MED_EDGE, mesh))
422 const SUPPORT *allEdges = _ptrMesh->getSupportOnAll( MED_EDGE );
423 (this->*writePart)( ensightGeomFile, allEdges);
426 // Write all groups as parts
428 for ( int ent = MED_CELL; ent < MED_ALL_ENTITIES; ++ent )
430 medEntityMesh entity = (medEntityMesh) ent;
431 int nbGroups = mesh->getNumberOfGroups(entity);
432 for ( int i=1; i<=nbGroups; i++)
434 const GROUP* group = mesh->getGroup( entity, i );
435 (this->*writePart)( ensightGeomFile, group );
439 ensightGeomFile.close();
441 mesh->removeReference();
443 } // end ASCII format
445 } // ENSIGHT_MESH_WRONLY_DRIVER::write()
447 //================================================================================
449 * \brief Write support as an EnSight Gold part
451 //================================================================================
453 void ENSIGHT_MESH_WRONLY_DRIVER::writePartGoldBinary(_BinaryFileWriter& ensightGeomFile,
454 const SUPPORT* support) const
457 int partNum = getPartNumber( support );
459 throw MED_EXCEPTION ( LOCALIZED( STRING("Internal error: invalid part number")));
460 ensightGeomFile.addString( "part" );
461 ensightGeomFile.addInt( partNum );
464 ensightGeomFile.addString( support->getName() );
467 medEntityMesh entity = support->getEntity();
468 int nbTypes = support->getNumberOfTypes();
469 const medGeometryElement* geoType = support->getTypes();
471 const int * connectivity = 0;
472 const int * elemConnectivity = 0;
473 const int * index = 0;
476 const MESH* mesh = support->getMesh()->convertInMESH() ; // we write unstructured only
478 // COORDINATES Gold binary
479 // ===================================================================================
480 // In Ensight, coordinates of nodes of support elements are in MED_NO_INTERLACE mode.
481 // We are to write only nodes belonging to elements of the support and
482 // nodal connectivity should refer to these nodes.
483 map<int, int> med2ensIds;
484 map<int, int>::iterator medEnsIt;
485 int SpaceDimension = _ptrMesh->getSpaceDimension() ;
486 int NumberOfNodes = _ptrMesh->getNumberOfNodes() ;
487 // -------------------------------------------------
488 if ( support->isOnAllElements() )
491 ensightGeomFile.addString( "coordinates" );
492 ensightGeomFile.addInt( NumberOfNodes );
495 const double *coordinate = mesh->getCoordinates(MED_FULL_INTERLACE);
496 typedef _ValueIterator< double > TComponentIt;
497 vector< TComponentIt > coordCompIt( 1 );
498 for (int j=0; j<SPACE_DIM; j++, coordinate++) { // loop on dimensions
499 if ( j < SpaceDimension )
500 coordCompIt[ 0 ] = TComponentIt( coordinate, SpaceDimension );
502 coordCompIt[ 0 ] = TComponentIt(); // to iterate on zeros
503 ensightGeomFile.addReal( coordCompIt, NumberOfNodes, MED_NO_INTERLACE );
506 // -------------------------------------------------
507 else // support is not on all elements
510 getSupportNodes( support, med2ensIds );
511 NumberOfNodes = med2ensIds.size();
512 ensightGeomFile.addString( "coordinates" );
513 ensightGeomFile.addInt( NumberOfNodes );
516 vector<float> floatCoords( NumberOfNodes );
517 for ( j=0; j < SPACE_DIM; j++) { // loop on dimensions
518 medEnsIt = med2ensIds.begin();
519 if ( j < SpaceDimension ) {
520 const double *coordinate = mesh->getCoordinates(MED_FULL_INTERLACE) + j;
521 for (int i=0; i<NumberOfNodes; i++, ++medEnsIt )
522 floatCoords[ i ] = (float) coordinate[ (medEnsIt->first-1) * SpaceDimension];
524 else if ( j-1 < SpaceDimension ) {
525 for (int i=0; i<NumberOfNodes; i++)
526 floatCoords[ i ] = 0.;
528 ensightGeomFile.addReal( &floatCoords[0], NumberOfNodes );
530 // assign local node ids
531 for ( medEnsIt = med2ensIds.begin(), j=1; j<=NumberOfNodes; j++, ++medEnsIt )
532 medEnsIt->second = j;
535 // CONNECTIVITY Gold binary
536 // ===================================================================================
538 for (int i=0; i<nbTypes; i++)
540 const medGeometryElement medType = geoType[i];
541 const TEnSightElemType& ensightType = getEnSightType(medType);
542 const int numberOfCell = support->getNumberOfElements(medType);
543 int nbCellNodes = ensightType._medIndex.size();
545 // type name and nb cells
546 ensightGeomFile.addString( ensightType._name );
547 ensightGeomFile.addInt( numberOfCell );
551 // -------------------------------------------------
552 if ( support->isOnAllElements() )
554 #ifdef ELEMENT_ID_GIVEN
556 nodeIds.resize( numberOfCell );
557 for ( j = 1; j <= numberOfCell; j++)
559 ensightGeomFile.addInt( nodeIds );
561 if ( entity != MED_NODE ) nodeIds.clear();
564 if ( nbCellNodes > 1 ) // STANDARD ELEMENTS connectivity
566 connectivity = mesh->getConnectivity( MED_NODAL, entity, medType);
567 nodeIds.reserve( numberOfCell * nbCellNodes);
568 for (j = 0 ; j < numberOfCell; j++, connectivity += nbCellNodes)
569 for (int k=0; k<nbCellNodes; k++)
570 nodeIds.push_back( connectivity[ ensightType._medIndex[k] ]);
571 ensightGeomFile.addInt( nodeIds );
573 else if ( entity == MED_NODE ) // NODES connectivity
575 #if !defined(ELEMENT_ID_GIVEN)
576 nodeIds.resize( numberOfCell );
577 for ( j = 1; j <= numberOfCell; j++)
580 ensightGeomFile.addInt( nodeIds );
582 else if ( medType == MED_POLYGON ) // POLYGONs connectivity
584 int nbStdCells = mesh->getGlobalNumberingIndex(entity)[i]-1;
585 connectivity = mesh->getConnectivity( MED_NODAL, entity, medType);
586 int connLength = mesh->getConnectivityLength( MED_NODAL, entity, medType);
587 index = mesh->getConnectivityIndex(MED_NODAL, entity);
589 // number of nodes in each element
591 TIntOwner nbNodesInPoly( getNumbersByIndex( index+nbStdCells, numberOfCell ));
592 ensightGeomFile.addInt( nbNodesInPoly.myValues, numberOfCell );
593 } // nbNodesInPoly is deleted here
596 ensightGeomFile.addInt( connectivity, connLength );
598 else // POLYHEDRA connectivity
600 connectivity = mesh->getConnectivity( MED_NODAL, entity, medType);
601 int nbStdCells = mesh->getGlobalNumberingIndex(entity)[i]-1;
602 index = mesh->getConnectivityIndex(MED_NODAL, entity) + nbStdCells;
604 vector<int> nbFacesInPolyhedron, nbNodesInFace, faceConn;
605 for ( int j = 0; j < numberOfCell; ++j )
607 int nbFaces = 0, nbNodes = 0;
608 for ( int k = index[j]; k < index[j+1]; ++k )
609 if ( connectivity[k-1] == -1 )
611 nbNodesInFace.push_back( nbNodes );
617 faceConn.push_back( connectivity[k-1] );
620 nbNodesInFace.push_back( nbNodes );
621 nbFacesInPolyhedron.push_back( nbFaces+1 );
624 // nb of faces in each polyhedron
625 ensightGeomFile.addInt( nbFacesInPolyhedron );
626 // number of nodes in each face
627 ensightGeomFile.addInt( nbNodesInFace );
629 ensightGeomFile.addInt( faceConn );
632 // -------------------------------------------------
633 else // support is not on all elements
635 const int *number = support->getNumber(medType);
637 #ifdef ELEMENT_ID_GIVEN
638 ensightGeomFile.addInt( number, numberOfCell );
640 if ( nbCellNodes > 1 ) // STANDARD ELEMENTS connectivity
642 connectivity = mesh->getConnectivity( MED_NODAL, entity, MED_ALL_ELEMENTS);
643 index = mesh->getConnectivityIndex(MED_NODAL, entity);
645 nodeIds.reserve( numberOfCell * nbCellNodes);
646 for (j=0; j<numberOfCell; j++) {
647 int elem = number[j];
648 elemConnectivity = connectivity + index[elem-1]-1;
649 for (int k=0; k<nbCellNodes; k++)
651 int node = elemConnectivity[ ensightType._medIndex[k] ];
652 nodeIds.push_back( med2ensIds[ node ]);
655 ensightGeomFile.addInt( nodeIds );
657 else if ( entity == MED_NODE ) // NODES connectivity
659 nodeIds.resize( numberOfCell );
660 for ( j = 1; j <= numberOfCell; j++)
662 ensightGeomFile.addInt( nodeIds );
664 else if ( medType == MED_POLYGON ) // POLYGONs connectivity
666 connectivity = mesh->getConnectivity( MED_NODAL, entity, MED_ALL_ELEMENTS);
667 index = mesh->getConnectivityIndex(MED_NODAL, entity);
669 // number of nodes in each element
671 TIntOwner nbNodesInPoly( getNumbersByIndex( index, numberOfCell, number ));
672 ensightGeomFile.addInt( nbNodesInPoly.myValues, numberOfCell );
673 } // nbNodesInPoly is deleted here
676 for ( j = 0; j < numberOfCell; ++j )
678 int elem = number[ j ];
679 elemConnectivity = connectivity + index[ elem-1 ]-1;
680 const int* connEnd = connectivity + index[ elem ]-1;
681 while ( elemConnectivity < connEnd )
682 nodeIds.push_back( med2ensIds[ *elemConnectivity++ ]);
684 ensightGeomFile.addInt( nodeIds );
686 else // POLYHEDRA connectivity
688 connectivity = mesh->getConnectivity( MED_NODAL, entity, MED_ALL_ELEMENTS);
689 index = mesh->getConnectivityIndex(MED_NODAL, entity);
690 vector<int> nbFacesInPolyhedron, nbNodesInFace, faceConn;
691 for ( int j = 0; j < numberOfCell; ++j )
693 int elem = number[ j ];
694 int nbFaces = 0, nbNodes = 0;
695 for ( int k = index[elem]; k < index[elem+1]; ++k )
696 if ( connectivity[k-1] == -1 )
698 nbNodesInFace.push_back( nbNodes );
704 faceConn.push_back( connectivity[k-1] );
707 nbNodesInFace.push_back( nbNodes );
708 nbFacesInPolyhedron.push_back( nbFaces+1 );
711 // nb of faces in each polyhedron
712 ensightGeomFile.addInt( nbFacesInPolyhedron );
713 // number of nodes in each face
714 ensightGeomFile.addInt( nbNodesInFace );
716 ensightGeomFile.addInt( faceConn );
721 mesh->removeReference();
723 } // writePartGoldBinary()
725 //================================================================================
727 * \brief Write support as an EnSight Gold part
729 //================================================================================
731 void ENSIGHT_MESH_WRONLY_DRIVER::writePartGoldASCII(ofstream& ensightGeomFile,
732 const SUPPORT* support) const
737 int partNum = getPartNumber( support );
738 ensightGeomFile << "part" << endl
739 << setw(iw) << partNum << endl;
741 throw MED_EXCEPTION ( LOCALIZED( STRING("Internal error: invalid part number")));
744 ensightGeomFile << support->getName() << endl;
747 medEntityMesh entity = support->getEntity();
748 int nbTypes = support->getNumberOfTypes();
749 const medGeometryElement* geoType = support->getTypes();
751 const int * connectivity = 0;
752 const int * elemConnectivity = 0;
753 const int * index = 0;
756 const MESH* mesh = support->getMesh()->convertInMESH() ; // we write unstructured only
758 // COORDINATES Gold ASCII
759 // ===================================================================================
760 // In Ensight, coordinates of nodes of support elements are in MED_NO_INTERLACE mode.
761 // We are to write only nodes belonging to elements of the support and
762 // nodal connectivity should refer to these nodes.
763 map<int, int> med2ensIds;
764 map<int, int>::iterator medEnsIt;
765 int SpaceDimension = mesh->getSpaceDimension() ;
766 int NumberOfNodes = mesh->getNumberOfNodes() ;
767 string zeroStr = " 0.00000e+00";
768 // -----------------------------------
769 if ( support->isOnAllElements() )
772 ensightGeomFile << "coordinates" << endl
773 << setw(iw) << NumberOfNodes << endl ;
776 for (j=0; j<SPACE_DIM; j++) { // loop on dimensions
777 if ( j < SpaceDimension ) {
778 const double *coordinate = mesh->getCoordinates(MED_FULL_INTERLACE) + j;
779 for (int i=0; i<NumberOfNodes; i++, coordinate += SpaceDimension)
780 ensightGeomFile << setw(12) << (float) *coordinate << endl;
783 for (int i=0; i<NumberOfNodes; i++)
784 ensightGeomFile << zeroStr << endl;
788 // -----------------------------------
789 else // support is not on all elements
792 getSupportNodes( support, med2ensIds );
793 NumberOfNodes = med2ensIds.size();
794 ensightGeomFile << "coordinates" << endl
795 << setw(iw) << NumberOfNodes << endl ;
798 for ( j=0; j<SPACE_DIM; j++) { // loop on dimensions
799 medEnsIt = med2ensIds.begin();
800 if ( j < SpaceDimension ) {
801 const double *coordinate = mesh->getCoordinates(MED_FULL_INTERLACE) + j;
802 for (int i=0; i<NumberOfNodes; i++, ++medEnsIt )
803 ensightGeomFile << setw(12)
804 << (float) coordinate[ (medEnsIt->first-1) * SpaceDimension] << endl;
807 for (int i=0; i<NumberOfNodes; i++)
808 ensightGeomFile << zeroStr << endl;
811 // assign local node ids
812 for ( medEnsIt = med2ensIds.begin(), j=1; j<=NumberOfNodes; j++, ++medEnsIt )
813 medEnsIt->second = j;
816 // CONNECTIVITY Gold ASCII
817 // ===================================================================================
819 for (int i=0; i<nbTypes; i++)
821 const medGeometryElement medType = geoType[i];
822 const TEnSightElemType& ensightType = getEnSightType(medType);
823 const int numberOfCell = support->getNumberOfElements(medType);
824 int nbCellNodes = ensightType._medIndex.size();
826 // type name and nb cells
827 ensightGeomFile << ensightType._name << endl
828 << setw(iw) << numberOfCell << endl;
830 // -----------------------------------
831 if ( support->isOnAllElements() )
833 #ifdef ELEMENT_ID_GIVEN
834 for ( j = 1; j <= numberOfCell; j++)
835 ensightGeomFile << setw(iw) << j << endl;
838 if ( nbCellNodes > 1 ) // STANDARD ELEMENTS connectivity
840 connectivity = mesh->getConnectivity( MED_NODAL, entity, medType);
841 for (j = 0 ; j < numberOfCell; j++, connectivity += nbCellNodes) {
842 for (int k=0; k<nbCellNodes; k++)
843 ensightGeomFile << setw(iw) << connectivity[ ensightType._medIndex[k] ];
844 ensightGeomFile << endl ;
847 else if ( entity == MED_NODE ) // NODES connectivity
849 for ( j = 1; j <= numberOfCell; j++)
850 ensightGeomFile << setw(iw) << j << endl;
852 else if ( medType == MED_POLYGON ) // POLYGONs connectivity
854 int nbStdCells = mesh->getGlobalNumberingIndex(entity)[i]-1;
855 connectivity = mesh->getConnectivity( MED_NODAL, entity, MED_ALL_ELEMENTS);
856 index = mesh->getConnectivityIndex(MED_NODAL, entity) + nbStdCells;
857 // number of nodes in each element
858 const int* ind = index;
859 for (j = 0 ; j < numberOfCell; j++, ++ind)
860 ensightGeomFile << setw(iw) << ( ind[1] - ind[0] ) << endl;
863 for (j = 0; j < numberOfCell; j++, ++index) {
864 nbCellNodes = index[1] - index[0];
865 for (int k=0; k<nbCellNodes; k++, ++connectivity)
866 ensightGeomFile << setw(iw) << *connectivity;
867 ensightGeomFile << endl;
870 else // POLYHEDRA connectivity
872 int nbStdCells = mesh->getGlobalNumberingIndex(entity)[i]-1;
873 connectivity = mesh->getConnectivity( MED_NODAL, entity, medType);
874 index = mesh->getConnectivityIndex(MED_NODAL, entity) + nbStdCells;
875 ostringstream nbNodesInFace, faceConn;
876 for ( j = 0; j < numberOfCell; ++j )
878 int nbFaces = 0, nbNodes = 0;
879 for ( int k = index[j]; k < index[j+1]; ++k )
880 if ( connectivity[k-1] == -1 )
883 nbNodesInFace << setw(iw) << nbNodes << endl;
889 faceConn << setw(iw) << connectivity[k-1] ;
893 nbNodesInFace << setw(iw) << nbNodes << endl;
894 ensightGeomFile << setw(iw) << nbFaces+1 << endl ;// nb of faces in each polyhedron
896 ensightGeomFile << nbNodesInFace.str();// number of nodes in each face
897 ensightGeomFile << faceConn.str();// connectivity of each face
900 // -----------------------------------
901 else // support is not on all elements
903 const int *number = support->getNumber(medType);
905 #ifdef ELEMENT_ID_GIVEN
906 for ( j = 0; j < numberOfCell; j++)
907 ensightGeomFile << setw(iw) << number[j] << endl;
909 if ( nbCellNodes > 1 ) // STANDARD ELEMENTS connectivity
911 connectivity = mesh->getConnectivity( MED_NODAL, entity, MED_ALL_ELEMENTS);
912 index = mesh->getConnectivityIndex(MED_NODAL, entity);
914 for (j=0; j<numberOfCell; j++) {
915 int elem = number[j];
916 elemConnectivity = connectivity + index[elem-1]-1;
917 for (int k=0; k<nbCellNodes; k++) {
918 int node = elemConnectivity[ ensightType._medIndex[k] ];
919 ensightGeomFile << setw(iw) << med2ensIds[ node ];
921 ensightGeomFile << endl;
924 else if ( entity == MED_NODE ) // NODES connectivity
926 for (j=0; j<numberOfCell; j++) {
927 int node = med2ensIds[ number[j] ];
928 ensightGeomFile << setw(iw) << node << endl ;
931 else if ( medType == MED_POLYGON ) // POLYGONs connectivity
933 connectivity = mesh->getConnectivity( MED_NODAL, entity, MED_ALL_ELEMENTS);
934 index = mesh->getConnectivityIndex(MED_NODAL, entity);
935 // number of nodes in each element
936 for (j = 0 ; j < numberOfCell; j++) {
937 int elem = number[j];
938 ensightGeomFile << setw(iw) << ( index[elem] - index[elem-1] ) << endl;
941 for ( j = 0; j < numberOfCell; ++j ) {
942 int elem = number[ j ];
943 elemConnectivity = connectivity + index[ elem-1 ]-1;
944 const int* connEnd = connectivity + index[ elem ]-1;
945 while ( elemConnectivity < connEnd )
946 ensightGeomFile << setw(iw) << med2ensIds[ *elemConnectivity++ ];
947 ensightGeomFile << endl;
950 else // POLYHEDRA connectivity
952 connectivity = mesh->getConnectivity( MED_NODAL, entity, medType);
953 index = mesh->getConnectivityIndex(MED_NODAL, entity);
954 ostringstream nbNodesInFace, faceConn;
955 for ( j = 0; j < numberOfCell; ++j )
957 int elem = number[j];
958 int nbFaces = 0, nbNodes = 0;
959 for ( int k = index[elem]; k < index[elem+1]; ++k )
960 if ( connectivity[k-1] == -1 )
963 nbNodesInFace << setw(iw) << nbNodes << endl;
969 faceConn << setw(iw) << connectivity[k-1] ;
973 nbNodesInFace << setw(iw) << nbNodes << endl;
974 ensightGeomFile << setw(iw) << nbFaces+1 << endl ;// nb of faces in each polyhedron
976 ensightGeomFile << nbNodesInFace.str();// number of nodes in each face
977 ensightGeomFile << faceConn.str();// connectivity of each face
982 mesh->removeReference();
984 } // writePartGoldASCII()
986 //================================================================================
988 * \brief Write support as an Ensight6 part
990 //================================================================================
992 void ENSIGHT_MESH_WRONLY_DRIVER::writePart6Binary(_BinaryFileWriter& ensightGeomFile,
993 const SUPPORT* support) const
996 int partNum = getPartNumber( support );
997 ensightGeomFile.addString( STRING("part ") << partNum );
999 throw MED_EXCEPTION ( LOCALIZED( STRING("Internal error: invalid part number")));
1002 ensightGeomFile.addString( support->getName() );
1005 medEntityMesh entity = support->getEntity();
1006 int nbTypes = support->getNumberOfTypes();
1007 const medGeometryElement* geoType = support->getTypes();
1009 const MESH* mesh = support->getMesh()->convertInMESH() ; // we write unstructured only
1012 const int * connectivity = 0;
1013 if ( entity != MED_NODE )
1014 connectivity = mesh->getConnectivity( MED_NODAL, entity, MED_ALL_ELEMENTS);
1015 const int * elemConnectivity = connectivity;
1017 // CONNECTIVITY Ensight 6 binary
1018 // ===================================================================================
1020 for (int i=0; i<nbTypes; i++)
1022 const medGeometryElement medType = geoType[i];
1023 const TEnSightElemType& ensightType = getEnSightType(medType);
1024 int nbCellNodes = ensightType._medIndex.size();
1025 if ( nbCellNodes == 0 )
1028 // type name and nb cells
1029 int numberOfCell = support->getNumberOfElements(medType);
1030 ensightGeomFile.addString( ensightType._name );
1031 ensightGeomFile.addInt( numberOfCell );
1033 vector<int> nodeIds;
1034 // -------------------------------------------------
1035 if ( support->isOnAllElements() )
1037 #ifdef ELEMENT_ID_GIVEN
1038 nodeIds.resize( numberOfCell );
1039 for ( j = 1; j <= numberOfCell; j++)
1041 ensightGeomFile.addInt( nodeIds );
1043 if ( entity == MED_NODE ) {
1044 #if !defined(ELEMENT_ID_GIVEN)
1045 nodeIds.resize( numberOfCell * nbCellNodes);
1046 for ( j = 1; j <= numberOfCell; j++)
1052 nodeIds.reserve( numberOfCell * nbCellNodes );
1053 for (j = 0 ; j < numberOfCell; j++, elemConnectivity += nbCellNodes)
1054 for (int k=0; k<nbCellNodes; k++)
1055 nodeIds.push_back( elemConnectivity[ ensightType._medIndex[k] ]);
1057 ensightGeomFile.addInt( nodeIds );
1059 // -------------------------------------------------
1060 else // support is not on all elements
1062 const int *number = support->getNumber(medType);
1064 #ifdef ELEMENT_ID_GIVEN
1065 ensightGeomFile.addInt( number, numberOfCell );
1067 if ( entity == MED_NODE ) {
1068 ensightGeomFile.addInt( number, numberOfCell );
1071 const int* index = mesh->getConnectivityIndex(MED_NODAL, entity);
1073 nodeIds.reserve( numberOfCell * nbCellNodes);
1074 for (j=0; j<numberOfCell; j++) {
1075 int elem = number[j];
1076 elemConnectivity = connectivity + index[elem-1]-1;
1077 for (int k=0; k<nbCellNodes; k++)
1078 nodeIds.push_back( elemConnectivity[ ensightType._medIndex[k] ]);
1080 ensightGeomFile.addInt( nodeIds );
1085 mesh->removeReference();
1088 } // writePart6Binary()
1090 //================================================================================
1092 * \brief Write support as an Ensight6 part
1094 //================================================================================
1096 void ENSIGHT_MESH_WRONLY_DRIVER::writePart6ASCII(ofstream& ensightGeomFile,
1097 const SUPPORT* support) const
1102 int partNum = getPartNumber( support );
1103 ensightGeomFile << "part " << partNum << endl;
1105 throw MED_EXCEPTION ( LOCALIZED( STRING("Internal error: invalid part number")));
1108 ensightGeomFile << support->getName() << endl;
1111 medEntityMesh entity = support->getEntity();
1112 int nbTypes = support->getNumberOfTypes();
1113 const medGeometryElement* geoType = support->getTypes();
1115 const MESH* mesh = support->getMesh()->convertInMESH() ; // we write unstructured only
1118 const int * connectivity = 0;
1119 if ( entity != MED_NODE )
1120 connectivity = mesh->getConnectivity( MED_NODAL,entity, MED_ALL_ELEMENTS);
1121 const int * elemConnectivity = connectivity;
1123 // CONNECTIVITY Ensight 6 ASCII
1124 // ===================================================================================
1126 for (int i=0; i<nbTypes; i++)
1128 const medGeometryElement medType = geoType[i];
1129 const TEnSightElemType& ensightType = getEnSightType(medType);
1130 int nbCellNodes = ensightType._medIndex.size();
1131 if ( nbCellNodes == 0 )
1134 // type name and nb cells
1135 int numberOfCell = support->getNumberOfElements(medType);
1136 ensightGeomFile << ensightType._name << endl
1137 << setw(iw) << numberOfCell << endl;
1139 // -------------------------------------------------
1140 if ( support->isOnAllElements() )
1142 if ( entity == MED_NODE ) {
1143 for ( j = 1; j <= numberOfCell; j++) {
1144 #ifdef ELEMENT_ID_GIVEN
1145 ensightGeomFile << setw(iw) << j;
1147 ensightGeomFile << setw(iw) << j << endl;
1151 for (j = 1 ; j <= numberOfCell; j++, elemConnectivity += nbCellNodes) {
1152 #ifdef ELEMENT_ID_GIVEN
1153 ensightGeomFile << setw(iw) << elem++;
1155 for (int k=0; k<nbCellNodes; k++)
1157 ensightGeomFile << setw(iw) << elemConnectivity[ ensightType._medIndex[k] ];
1159 ensightGeomFile << endl ;
1163 // -------------------------------------------------
1164 else // support is not on all elements
1166 const int *number = support->getNumber(medType);
1167 if ( entity == MED_NODE ) {
1168 for (j=0; j<numberOfCell; j++) {
1169 int node = number[j];
1170 #ifdef ELEMENT_ID_GIVEN
1171 ensightGeomFile << setw(iw) << node;
1173 ensightGeomFile << setw(iw) << node << endl ;
1177 const int* index = mesh->getConnectivityIndex(MED_NODAL, entity);
1179 for (j=0; j<numberOfCell; j++) {
1180 int elem = number[j];
1181 #ifdef ELEMENT_ID_GIVEN
1182 ensightGeomFile << setw(iw) << elem;
1184 elemConnectivity = connectivity + index[elem-1]-1;
1185 for (int k=0; k<nbCellNodes; k++)
1187 ensightGeomFile << setw(iw) << elemConnectivity[ ensightType._medIndex[k] ];
1189 ensightGeomFile << endl ;
1195 mesh->removeReference();
1198 } // writePart6ASCII()
1200 //================================================================================
1202 * \brief Return nb of part to write
1204 //================================================================================
1206 int ENSIGHT_MESH_WRONLY_DRIVER::nbPartsToWrite() const
1209 nbParts += (int) isToWriteEntity( MED_CELL, _ptrMesh );
1210 nbParts += (int) isToWriteEntity( MED_FACE, _ptrMesh );
1211 nbParts += (int) isToWriteEntity( MED_EDGE, _ptrMesh );
1214 for ( int ent = MED_CELL; ent < MED_ALL_ENTITIES; ++ent ) {
1215 int nbGroups = _ptrMesh->getNumberOfGroups(medEntityMesh(ent));
1216 nbParts += nbGroups;
1221 // ================================================================================
1223 // ================================================================================
1225 ENSIGHT_MESH_RDONLY_DRIVER::ENSIGHT_MESH_RDONLY_DRIVER()
1226 : ENSIGHT_MESH_DRIVER(), _indexInCaseFile(1)
1230 ENSIGHT_MESH_RDONLY_DRIVER::ENSIGHT_MESH_RDONLY_DRIVER(const string & fileName,
1233 : ENSIGHT_MESH_DRIVER(fileName,ptrMesh,RDONLY), _indexInCaseFile( index )
1237 ENSIGHT_MESH_RDONLY_DRIVER::ENSIGHT_MESH_RDONLY_DRIVER(const ENSIGHT_MESH_RDONLY_DRIVER & driver) : ENSIGHT_MESH_DRIVER(driver), _indexInCaseFile( driver._indexInCaseFile )
1241 ENSIGHT_MESH_RDONLY_DRIVER::~ENSIGHT_MESH_RDONLY_DRIVER()
1245 GENDRIVER * ENSIGHT_MESH_RDONLY_DRIVER::copy() const
1247 return new ENSIGHT_MESH_RDONLY_DRIVER(*this) ;
1250 void ENSIGHT_MESH_RDONLY_DRIVER::write() const throw (MEDEXCEPTION)
1252 throw MEDEXCEPTION("ENSIGHT_MESH_RDONLY_DRIVER::write : Can't write with a RDONLY driver !");
1255 void ENSIGHT_MESH_RDONLY_DRIVER::merge ( const GENDRIVER& driver )
1257 _CaseFileDriver_User::merge( driver );
1259 const ENSIGHT_MESH_RDONLY_DRIVER* other =
1260 dynamic_cast< const ENSIGHT_MESH_RDONLY_DRIVER* >( &driver );
1262 if ( _indexInCaseFile < other->_indexInCaseFile )
1263 _indexInCaseFile = other->_indexInCaseFile;
1267 //================================================================================
1269 * \brief Read mesh in all supported formats
1271 //================================================================================
1273 void ENSIGHT_MESH_RDONLY_DRIVER::read() throw (MEDEXCEPTION)
1275 const char * LOC = "ENSIGHT_MESH_RDONLY_DRIVER::read() : " ;
1278 openConst(false); // check if can read case file
1280 _CaseFileDriver caseFile( getCaseFileName(), this);
1282 caseFile.setDataFileName( _indexInCaseFile, this ); // data from Case File is passed here
1284 openConst(true); // check if can read data file
1286 cout << "-> Entering into the geometry file " << getDataFileName() << endl ;
1288 MESH* mesh = (MESH*) getMesh();
1290 _InterMed* imed = new _InterMed();
1291 imed->_medMesh = mesh;
1292 imed->_isOwnMedMesh = false;
1293 imed->_needSubParts = ( caseFile.getNbVariables() > 0 );
1294 imed->groupes.reserve(1000);
1296 // to let field drivers know eventual indices of values
1297 if ( imed->_needSubParts )
1298 setInterData( imed );
1300 if ( isBinaryDataFile( getDataFileName() )) // binary
1302 if ( isGoldFormat() ) // Gold
1304 readGoldBinary( *imed );
1308 read6Binary( *imed );
1313 if ( isGoldFormat() ) // Gold
1315 readGoldASCII( *imed );
1319 read6ASCII( *imed );
1323 if ( _isMadeByMed && !imed->groupes.empty() ) {
1324 mesh->_name = imed->groupes[0].nom;
1325 imed->groupes[0].nom = "SupportOnAll_";
1326 imed->groupes[0].nom += entNames[MED_CELL];
1329 mesh->_name = theDefaultMeshName;
1331 mesh->_spaceDimension = SPACE_DIM;
1332 mesh->_numberOfNodes = imed->points.size() - imed->nbMerged( MED_POINT1 );
1333 mesh->_coordinate = imed->getCoordinate();
1335 //Construction des groupes
1336 imed->getGroups(mesh->_groupCell,
1339 mesh->_groupNode, mesh);
1341 mesh->_connectivity = imed->getConnectivity();
1343 mesh->createFamilies();
1345 // add attributes to families
1346 set<string> famNames;
1347 for (medEntityMesh entity=MED_CELL; entity<MED_ALL_ENTITIES; ++entity)
1349 int i, nb = mesh->getNumberOfFamilies(entity);
1350 for ( i = 1; i <= nb; ++i ) {
1351 FAMILY* f = const_cast<FAMILY*>( mesh->getFamily( entity, i ));
1352 f->setNumberOfAttributes( 1 );
1353 int* attIDs = new int[1];
1355 f->setAttributesIdentifiers( attIDs );
1356 int* attVals = new int[1];
1358 f->setAttributesValues( attVals );
1359 string* attDescr = new string[1];
1360 attDescr[0] = "med_family";
1361 f->setAttributesDescriptions( attDescr );
1363 if ( f->getName().length() > 31 ) // limit a name length
1364 f->setName( STRING("FAM_") << f->getIdentifier());
1365 // setAll() for groups
1366 nb = mesh->getNumberOfGroups(entity);
1367 for ( i = 1; i <= nb; ++i ) {
1368 GROUP * g = const_cast<GROUP*>( mesh->getGroup( entity, i ));
1369 if (mesh->getNumberOfElements( entity, MED_ALL_ELEMENTS ) ==
1370 g->getNumberOfElements( MED_ALL_ELEMENTS ))
1379 if ( !imed->_needSubParts )
1385 //================================================================================
1387 * \brief Read mesh in Gold ASCII format
1389 //================================================================================
1391 void ENSIGHT_MESH_RDONLY_DRIVER::readGoldASCII(_InterMed & imed)
1393 const char * LOC ="ENSIGHT_MESH_RDONLY_DRIVER::readGoldASCII() : ";
1396 _ASCIIFileReader geoFile( getDataFileName() );
1398 if ( isSingleFileMode() ) {
1399 int curTimeStep = 1;
1400 while ( curTimeStep++ < getIndexInDataFile() ) {
1401 while ( !geoFile.isTimeStepEnd())
1404 while ( !geoFile.isTimeStepBeginning() )
1407 // ----------------------
1408 // Read mesh description
1409 // ----------------------
1411 string descriptionLine1 = geoFile.getLine();
1412 string descriptionLine2 = geoFile.getLine();
1414 // find out if the file was created by MED driver
1415 int prefixSize = strlen( theDescriptionPrefix );
1416 _isMadeByMed = ( descriptionLine1.substr(0, prefixSize ) == theDescriptionPrefix );
1419 descriptionLine1 = descriptionLine1.substr( prefixSize );
1420 _ptrMesh->setDescription( descriptionLine1 + descriptionLine2 );
1423 // ----------------------------------------
1424 // Find out presence of node/elem numbers
1425 // ----------------------------------------
1427 // EnSight User Manual (for v8) says:
1428 // You do not have to assign node IDs. If you do, the element connectivities are
1429 // based on the node numbers. If you let EnSight assign the node IDs, the nodes
1430 // are considered to be sequential starting at node 1, and element connectivity is
1431 // done accordingly. If node IDs are set to off, they are numbered internally;
1432 // however, you will not be able to display or query on them. If you have node
1433 // IDs in your data, you can have EnSight ignore them by specifying "node id
1434 // ignore." Using this option may reduce some of the memory taken up by the
1435 // Client and Server, but display and query on the nodes will not be available.
1437 // read "node|element id <off|given|assign|ignore>"
1438 geoFile.getWord(); geoFile.getWord();
1439 string nodeIds = geoFile.getWord();
1440 geoFile.getWord(); geoFile.getWord();
1441 string elemIds = geoFile.getWord();
1443 bool haveNodeIds = ( nodeIds == "given" || nodeIds == "ignore" );
1444 bool haveElemIds = ( elemIds == "given" || elemIds == "ignore" );
1446 // extents: xmin xmax ymin ymax zmin zmax
1447 vector<double> extents;
1448 geoFile.toNextLine();
1449 if ( strncmp( "extents", geoFile.getCurrentPtr(), 7 ) == 0 ) {
1450 geoFile.skip( /*width =*/ 7, /*nbLines =*/ 1 );
1451 extents.reserve( 6 );
1452 while ( extents.size() < extents.capacity() )
1453 extents.push_back( geoFile.getReal() );
1456 typedef map<int,_noeud>::iterator INoeud;
1457 map<int,_noeud> & points = imed.points;
1460 _groupe * partGroupe = 0;
1461 int partNum = 0, nbParts = 0;
1463 while ( !geoFile.isTimeStepEnd() )
1465 string word, restLine, line = geoFile.getLine();
1466 TStrTool::split( line, word, restLine );
1468 const TEnSightElemType & partType = getEnSightType( word );
1469 if ( partType._medType != MED_ALL_ELEMENTS )
1471 // Unstructured element type encounters
1472 // --------------------------------------
1473 int nbElemNodes = partType._medType % 100;
1474 int nbElems = geoFile.getInt(); // ne
1475 bool isGhost = isGhostType( word );
1476 int nodeShift = points.empty() ? 0 : points.rbegin()->first;
1479 vector<int> elemIds;
1480 if ( haveElemIds ) {
1481 elemIds.reserve( nbElems );
1482 while ((int) elemIds.size() < nbElems )
1483 elemIds.push_back( geoFile.getInt() ); // id_e
1485 if ( isGhost ) { // do not store ghost elements (?)
1486 int nbNodes = nbElems * nbElemNodes;
1487 if ( partType._name == "nsided" ) // polygons
1489 for ( int i = 0; i < nbElems; ++i )
1490 nbNodes += geoFile.getInt();
1491 geoFile.skip( nbNodes * INT_WIDTH_GOLD, /*nbLines = */nbElems );
1493 else if ( partType._name == "nfaced" ) // polyhedrons
1496 for ( int i = 0; i < nbElems; ++i )
1497 nbFaces += geoFile.getInt();
1498 for ( int f = 0; f < nbFaces; ++f )
1499 nbNodes += geoFile.getInt();
1500 geoFile.skip( nbNodes * INT_WIDTH_GOLD, /*nbLines = */nbFaces );
1502 else // standard types
1504 geoFile.skip( nbNodes, nbElemNodes, INT_WIDTH_GOLD );
1509 // add a group corresponding to subPart (geoType)
1510 imed.groupes.push_back(_groupe());
1511 _groupe & groupe = imed.groupes.back();
1512 groupe.mailles.resize( nbElems );
1514 // find out if "coordinates" has already been encountered
1515 _SubPartDesc coordDesc( partNum , "coordinates");
1516 map< _SubPartDesc, _SubPart >::iterator descPart =
1517 imed._subPartDescribed.find( coordDesc );
1518 bool haveCoords = ( descPart != imed._subPartDescribed.end() );
1520 firstNode = descPart->second.myFirstNode;
1521 nodeShift -= descPart->second.myNbNodes;
1524 // read poly element data
1525 bool isPoly = ( !nbElemNodes );
1526 vector<int> nbElemNodesVec( 1, nbElemNodes);
1527 vector<int> nbElemFaces, nbFaceNodes;
1528 if ( partType._name == "nsided" ) // polygons
1530 nbElemNodesVec.resize( nbElems );
1531 for ( int i = 0; i < nbElems; ++i )
1532 nbElemNodesVec[ i ] = geoFile.getInt(); // npi
1534 else if ( partType._name == "nfaced" ) // polyhedrons
1536 nbElemFaces.resize( nbElems );
1537 nbElemNodesVec.resize( nbElems );
1538 int totalNbFaces = 0;
1539 for ( int i = 0; i < nbElems; ++i )
1540 totalNbFaces += ( nbElemFaces[ i ] = geoFile.getInt() ); // nf_ei
1542 nbFaceNodes.resize( totalNbFaces );
1543 vector<int>::iterator nbFN = nbFaceNodes.begin();
1544 for ( int i = 0; i < nbElems; ++i ) {
1545 nbElemNodesVec[ i ] = 0;
1546 for ( int nbFaces = nbElemFaces[ i ]; nbFaces; --nbFaces, ++nbFN )
1547 nbElemNodesVec[ i ] += ( *nbFN = geoFile.getInt() ); // np(f_ei)
1550 // iterator returning nbElemNodes for standard elems and
1551 // next value from nbElemNodesVec for poly elements
1552 _ValueIterator<int> nbElemNodesIt( & nbElemNodesVec[0], isPoly ? 1 : 0);
1554 // iterator returning values form partType._medIndex for standard elems
1555 // and node index (n) for poly elements
1557 _ValueIterator<int> medIndexIt( isPoly ? &n : &partType._medIndex[0],
1559 // read connectivity
1560 _maille ma( partType._medType, nbElemNodes );
1561 ma.sommets.resize( nbElemNodes );
1563 for ( int i = 0; i < nbElems; ++i ) {
1564 _ValueIterator<int> medIndex = medIndexIt;
1565 nbElemNodes = nbElemNodesIt.next();
1566 if ((int) ma.sommets.size() != nbElemNodes )
1567 ma.sommets.resize( nbElemNodes );
1568 for ( n = 0; n < nbElemNodes; ++n ) {
1569 int nodeID = geoFile.getInt(); // nn_ei
1571 node = points.find( nodeID + nodeShift );
1573 node = points.insert( make_pair( nodeID + nodeShift, _noeud())).first;
1574 ma.sommets[ medIndex.next() ] = node;
1577 ma.setOrdre( elemIds[ i ] );
1578 groupe.mailles[i] = imed.insert(ma);
1580 // store nb nodes in polyhedron faces
1581 if ( !nbFaceNodes.empty() ) {
1582 const int* nbFaceNodesPtr = & nbFaceNodes[0];
1583 for ( int i = 0; i < nbElems; ++i ) {
1584 vector<int> & nbNodesByFace = imed.polyherdalNbFaceNodes[ &groupe.maille( i ) ];
1585 nbNodesByFace.assign( nbFaceNodesPtr, nbFaceNodesPtr + nbElemFaces[ i ] );
1586 nbFaceNodesPtr += nbElemFaces[ i ];
1589 // create subPart for "coordinates"
1590 if ( !haveCoords ) {
1591 _SubPart & coordSubPart = imed._subPartDescribed[ coordDesc ];
1592 coordSubPart.myFirstNode = points.insert( make_pair( 1 + nodeShift, _noeud())).first;
1594 // add subPart group to part group
1595 int groupeIndex = imed.groupes.size();
1596 partGroupe->groupes.push_back( groupeIndex );
1599 _SubPart subPart( partNum, partType._name );
1600 subPart.myNbCells = nbElems;
1601 subPart.myCellGroupIndex = groupeIndex;
1602 imed.addSubPart( subPart );
1604 else if ( word == "coordinates" )
1606 // Local node coordinates of a part
1607 // ------------------------------------
1608 int nbNodes = geoFile.getInt(); // nn
1611 vector<int> nodeIds;
1612 if ( haveNodeIds ) {
1613 nodeIds.reserve( nbNodes );
1614 while ((int) nodeIds.size() < nbNodes )
1615 nodeIds.push_back( geoFile.getInt() ); // id_n
1618 // find out if "coordinates" has already been add at reading connectivity
1619 _SubPartDesc coordDesc( partNum , "coordinates");
1620 map< _SubPartDesc, _SubPart >::iterator descPart =
1621 imed._subPartDescribed.find( coordDesc );
1622 bool haveCoords = ( descPart != imed._subPartDescribed.end() );
1625 // check that all nodes have been added
1626 firstNode = descPart->second.myFirstNode;
1627 descPart->second.myNbNodes = nbNodes;
1628 INoeud inoeud = firstNode, inoEnd = points.end();
1629 int id = inoeud->first, idEnd = id + nbNodes;
1630 for ( ; id < idEnd; ++id ) {
1631 if ( inoeud == inoEnd || inoeud->first > id ) {
1632 INoeud in = points.insert( inoeud, make_pair( id, _noeud() ));
1633 in->second.number = id;
1634 in->second.coord.resize( SPACE_DIM );
1642 int nodeShift = points.empty() ? 0 : points.rbegin()->first;
1643 for ( int iNode = 1; iNode <= nbNodes; ++iNode ) {
1644 INoeud inoeud = points.insert( points.end(), make_pair( iNode + nodeShift, _noeud()));
1645 inoeud->second.number = inoeud->first;
1646 inoeud->second.coord.resize( SPACE_DIM );
1648 firstNode = points.find( 1 + nodeShift );
1649 // create "coordinates" subPart
1650 _SubPart & subPart = imed._subPartDescribed[ coordDesc ];
1651 subPart.myNbNodes = nbNodes;
1652 subPart.myFirstNode = firstNode;
1655 // read coordinates in no interlace mode
1656 INoeud endNode = points.end();
1657 for ( int j = 0; j < SPACE_DIM; ++j ) {
1658 for ( INoeud in = firstNode; in != endNode; ++in ) {
1659 _noeud & node = in->second;
1660 node.coord[ j ] = geoFile.getReal();
1664 else if ( word == "part" )
1666 // Another part encounters
1667 // -----------------------
1668 partNum = geoFile.getInt();
1670 geoFile.toNextLine();
1672 string partName = geoFile.getLine();
1673 if ( partName.empty() )
1674 partName = "Part_" + restLine;
1676 if (int( imed.groupes.capacity() - imed.groupes.size() ) < theMaxNbTypes )
1677 imed.groupes.reserve( size_t( 1.5 * imed.groupes.size() ));
1678 imed.groupes.push_back(_groupe());
1679 partGroupe = & imed.groupes.back();
1680 partGroupe->nom = partName;
1681 partGroupe->groupes.reserve( theMaxNbTypes );
1683 else if ( word == "block" )
1686 // ------------------
1687 bool rectilinear = ( restLine.find( "rectilinear" ) != restLine.npos );
1688 bool uniform = ( restLine.find( "uniform" ) != restLine.npos );
1689 bool curvilinear = ( !rectilinear && !uniform );
1690 bool iblanked = ( restLine.find( "iblanked" ) != restLine.npos );
1691 bool with_ghost = ( restLine.find( "with_ghost" ) != restLine.npos );
1692 bool range = ( restLine.find( "range" ) != restLine.npos );
1695 int I = geoFile.getInt();
1696 int J = geoFile.getInt();
1697 int K = geoFile.getInt();
1698 int NumberOfNodes = I*J*K;
1699 if ( !NumberOfNodes ) continue;
1703 vector<int> ijkRange; // imin imax jmin jmax kmin kmax
1704 ijkRange.reserve(6);
1705 while ( ijkRange.size() < 6 )
1706 ijkRange.push_back( geoFile.getInt() );
1707 I = ijkRange[1]-ijkRange[0]+1;
1708 J = ijkRange[3]-ijkRange[2]+1;
1709 K = ijkRange[5]-ijkRange[4]+1;
1710 NumberOfNodes = I*J*K;
1713 int nodeShift = points.empty() ? 0 : points.rbegin()->first;
1714 for ( int iNode = 1; iNode <= NumberOfNodes; ++iNode ) {
1715 INoeud inoeud = points.insert( points.end(), make_pair( iNode + nodeShift, _noeud()));
1716 _noeud & node = inoeud->second;
1717 node.number = inoeud->first;
1718 node.coord.resize( SPACE_DIM );
1720 INoeud firstNode = points.find( nodeShift + 1 );
1721 INoeud endNode = points.end();
1723 GRID grid; // calculator of unstructured data
1724 grid._iArrayLength = I;
1725 grid._jArrayLength = J;
1726 grid._kArrayLength = K;
1727 grid._spaceDimension = SPACE_DIM;
1728 if ( J < 2 ) { grid._spaceDimension--; grid._jArrayLength = 0; }
1729 if ( K < 2 ) { grid._spaceDimension--; grid._kArrayLength = 0; }
1730 int nbElems = grid.getNumberOfElements(MED_CELL, MED_ALL_ELEMENTS);
1732 if ( curvilinear ) // read coordinates for all nodes
1734 for ( int j = 0; j < SPACE_DIM; ++j ) {
1735 for ( INoeud in = firstNode; in != endNode; ++in )
1736 in->second.coord[ j ] = geoFile.getReal();
1738 grid._gridType = MED_BODY_FITTED;
1739 grid._coordinate = new COORDINATE(3, 1, 0); // to avoid exception in convertInMESH()
1741 else if ( rectilinear ) // read delta vectors with non-regular spacing
1743 grid._iArray = (double*)geoFile.convertReals<double>( I );
1744 grid._jArray = (double*)geoFile.convertReals<double>( J );
1745 grid._kArray = (double*)geoFile.convertReals<double>( K );
1746 grid._gridType = MED_CARTESIAN;
1748 else // uniform: read grid origine and delta vectors for regular spacing grid
1750 TDblOwner xyzOrigin( (double*)geoFile.convertReals<double>( 3 ));
1751 TDblOwner xyzDelta ( (double*)geoFile.convertReals<double>( 3 ));
1752 // compute full delta vectors
1753 grid._iArray = new double[ I ];
1754 grid._jArray = new double[ J ];
1755 grid._kArray = new double[ K ];
1756 double* coors[SPACE_DIM] = { grid._iArray, grid._jArray, grid._kArray };
1757 int size [SPACE_DIM] = { I, J, K };
1758 for ( int j = 0; j < SPACE_DIM; ++j ) {
1759 double* coo = coors[ j ];
1760 double* cooEnd = coo + size[ j ];
1761 coo[0] = xyzOrigin[ j ];
1762 while ( ++coo < cooEnd )
1763 *coo = coo[-1] + xyzDelta[ j ];
1765 grid._gridType = MED_CARTESIAN;
1770 geoFile.skip( NumberOfNodes, /*nbPerLine =*/ 1, INT_WIDTH_GOLD);
1773 geoFile.getWord(); // "ghost_flags"
1774 geoFile.skip( nbElems, /*nbPerLine =*/ 1, INT_WIDTH_GOLD);
1777 if ( haveNodeIds && geoFile.lookAt( "node_ids" )) {
1778 geoFile.getWord(); // "node_ids"
1779 geoFile.skip( NumberOfNodes, /*nbPerLine =*/ 1, INT_WIDTH_GOLD);
1782 if ( haveElemIds && geoFile.lookAt( "element_ids" ) ) {
1783 geoFile.getWord(); // "element_ids"
1784 geoFile.skip( nbElems, /*nbPerLine =*/ 1, INT_WIDTH_GOLD);
1787 // let GRID compute all coordinates and connectivity
1788 const MESH* unstruct = grid.convertInMESH();
1791 const double * coo = unstruct->getCoordinates(MED_FULL_INTERLACE);
1792 typedef _ValueIterator< double > TCoordIt;
1793 TCoordIt xCoo( coo+0, grid._spaceDimension);
1794 TCoordIt yCoo( coo+1, grid._spaceDimension);
1795 TCoordIt zCoo( coo+2, grid._spaceDimension);
1796 if ( grid._spaceDimension < 3 ) zCoo = TCoordIt( grid._kArray, 0 );
1797 if ( grid._spaceDimension < 2 ) yCoo = TCoordIt( grid._jArray, 0 );
1798 for ( INoeud in = firstNode; in != endNode; ++in ) {
1799 _noeud& node = in->second;
1800 node.coord[ 0 ] = xCoo.next();
1801 node.coord[ 1 ] = yCoo.next();
1802 node.coord[ 2 ] = zCoo.next();
1806 // store connectivity
1807 const int * conn = unstruct->getConnectivity( MED_NODAL, MED_CELL, MED_ALL_ELEMENTS );
1808 medGeometryElement elemType = grid.getElementType( MED_CELL, 1 );
1809 int nbElemNodes = elemType % 100;
1811 partGroupe->mailles.resize( nbElems );
1812 _maille ma( elemType, nbElemNodes );
1813 ma.sommets.resize( nbElemNodes );
1815 for ( int i = 0, nIndex = 0; i < nbElems; ++i ) {
1816 for ( int n = 0; n < nbElemNodes; ++n ) {
1817 int nodeID = conn[ nIndex++ ];
1818 ma.sommets[n] = points.find( nodeID + nodeShift );
1820 //ma.ordre = ++order;
1821 partGroupe->mailles[i] = imed.insert(ma);
1823 _SubPart subPart( partNum, "block" );
1824 subPart.myNbCells = nbElems;
1825 subPart.myNbNodes = NumberOfNodes;
1826 subPart.myFirstNode = firstNode;
1827 subPart.myCellGroupIndex = imed.groupes.size();
1828 imed.addSubPart( subPart );
1830 unstruct->removeReference();
1835 ( LOCALIZED( STRING(LOC) << "Unexpected word: " << word <<
1836 " in " << getDataFileName()));
1838 } // while ( !geoFile.eof() )
1841 imed.mergeNodesAndElements(TOLERANCE);
1846 //================================================================================
1848 * \brief Read mesh in Gold Binary format
1850 //================================================================================
1852 void ENSIGHT_MESH_RDONLY_DRIVER::readGoldBinary(_InterMed & imed)
1854 const char * LOC ="ENSIGHT_MESH_RDONLY_DRIVER::readGoldBinary() : ";
1857 _BinaryFileReader geoFile( getDataFileName() );
1859 // check if swapping bytes needed
1861 countPartsBinary( geoFile, isSingleFileMode());
1864 geoFile.swapBytes();
1867 if ( getIndexInDataFile() <= 1 )
1869 if ( geoFile.getPosition() == 0 ) {
1870 TStrOwner format( geoFile.getLine() ); // "C|Fortran Binary"
1871 if ( !contains( "C Binary", format )) {
1872 if ( contains( "Fortran Binary", format ))
1873 throw(MEDEXCEPTION(STRING(LOC) << "Fortran Binary format not supported"));
1875 throw(MEDEXCEPTION(STRING(LOC) << "unexpected line in " << getDataFileName()
1876 << "\n" << format.myValues));
1879 if ( isSingleFileMode() ) {
1880 // one time step may be skipped by countPartsBinary
1881 int curTimeStep = geoFile.getPosition() ? 2 : 1 ;
1882 while ( curTimeStep < getIndexInDataFile() ) {
1883 countPartsBinary( geoFile, true ); // skip time step
1887 TStrOwner line( geoFile.getLine() );
1888 if ( isTimeStepBeginning( line.myValues ))
1892 // ----------------------
1893 // Read mesh description
1894 // ----------------------
1896 TStrOwner descriptionLine1 ( geoFile.getLine() );
1897 TStrOwner descriptionLine2 ( geoFile.getLine() );
1899 // find out if the file was created by MED driver
1900 _isMadeByMed = contains( theDescriptionPrefix, descriptionLine1 );
1903 _ptrMesh->setDescription( descriptionLine2.myValues );
1905 _ptrMesh->setDescription( string(descriptionLine1) + descriptionLine2.myValues );
1908 // ----------------------------------------
1909 // Find out presence of node/elem numbers
1910 // ----------------------------------------
1912 // EnSight User Manual (for v8) says:
1913 // You do not have to assign node IDs. If you do, the element connectivities are
1914 // based on the node numbers. If you let EnSight assign the node IDs, the nodes
1915 // are considered to be sequential starting at node 1, and element connectivity is
1916 // done accordingly. If node IDs are set to off, they are numbered internally;
1917 // however, you will not be able to display or query on them. If you have node
1918 // IDs in your data, you can have EnSight ignore them by specifying "node id
1919 // ignore." Using this option may reduce some of the memory taken up by the
1920 // Client and Server, but display and query on the nodes will not be available.
1922 // read "node|element id <off|given|assign|ignore>"
1923 bool haveNodeIds, haveElemIds;
1925 TStrOwner nodeIds( geoFile.getLine() );
1926 TStrOwner elemIds( geoFile.getLine() );
1928 haveNodeIds = ( contains( "given", nodeIds ) || contains( "ignore", nodeIds ) );
1929 haveElemIds = ( contains( "given", elemIds ) || contains( "ignore", elemIds ) );
1932 typedef map<int,_noeud>::iterator INoeud;
1933 map<int,_noeud> & points = imed.points;
1936 _groupe * partGroupe = 0;
1937 int partNum = 0, nbParts = 0;
1939 TFltOwner extents(0); // extents: xmin xmax ymin ymax zmin zmax
1941 while ( !geoFile.eof() )
1943 TStrOwner line( geoFile.getLine() );
1944 if ( isSingleFileMode() && isTimeStepEnd( line.myValues ))
1946 string word, restLine;
1947 TStrTool::split( line.myValues, word, restLine );
1949 const TEnSightElemType & partType = getEnSightType( word );
1950 if ( partType._medType != MED_ALL_ELEMENTS )
1952 // Unstructured element type encounters
1953 // --------------------------------------
1954 int nbElems = *TIntOwner( geoFile.getInt(1) ); // ne
1955 int nbElemNodes = partType._medType % 100;
1956 bool isGhost = isGhostType( word );
1957 int nodeShift = points.empty() ? 0 : points.rbegin()->first;
1960 TIntOwner elemIds( haveElemIds ? geoFile.getInt( nbElems ): 0 ); // id_e
1962 if ( isGhost ) { // do not store ghost elements (?)
1963 int nbNodes = nbElems * nbElemNodes;
1964 if ( partType._name == "nsided" ) // polygons
1966 TIntOwner nbNodesInFace( geoFile.getInt( nbElems ));
1967 nbNodes = std::accumulate( nbNodesInFace.myValues, nbNodesInFace.myValues + nbElems, 0 );
1969 else if ( partType._name == "nfaced" ) // polyhedrons
1971 TIntOwner nbElemFaces( geoFile.getInt( nbElems ));
1972 int nbFaces = accumulate( nbElemFaces.myValues, nbElemFaces.myValues + nbElems, 0 );
1973 TIntOwner nbNodesInFace( geoFile.getInt( nbFaces ));
1974 nbNodes = std::accumulate( nbNodesInFace.myValues, nbNodesInFace.myValues + nbElems, 0 );
1976 geoFile.skip( nbNodes * sizeof(int) );
1980 // add a group corresponding to subPart (geoType)
1981 imed.groupes.push_back(_groupe());
1982 _groupe & groupe = imed.groupes.back();
1983 groupe.mailles.resize( nbElems );
1985 // find out if "coordinates" has already been encountered
1986 _SubPartDesc coordDesc( partNum , "coordinates");
1987 map< _SubPartDesc, _SubPart >::iterator descPart =
1988 imed._subPartDescribed.find( coordDesc );
1989 bool haveCoords = ( descPart != imed._subPartDescribed.end() );
1991 firstNode = descPart->second.myFirstNode;
1992 nodeShift -= descPart->second.myNbNodes;
1995 // read poly element data
1996 bool isPoly = ( !nbElemNodes );
1998 TIntOwner nbElemNodesVec(0), nbElemFaces(0), nbFaceNodes(0);
1999 if ( partType._name == "nsided" ) // polygons
2001 nbElemNodesVec.myValues = geoFile.getInt( nbElems ); // npi
2002 nbNodes = accumulate( nbElemNodesVec.myValues, nbElemNodesVec.myValues + nbElems, 0 );
2004 else if ( partType._name == "nfaced" ) // polyhedrons
2006 nbElemFaces.myValues = geoFile.getInt( nbElems ); // nf_ei
2007 int totalNbFaces = accumulate( nbElemFaces.myValues, nbElemFaces.myValues + nbElems, 0 );
2009 nbFaceNodes.myValues = geoFile.getInt( totalNbFaces ); // np(f_ei)
2010 // calculate nb of nodes in each polyhedron
2011 int* nbEN = nbElemNodesVec.myValues = new int[ nbElems ];
2012 const int *nbFN = nbFaceNodes, *nbEF = nbElemFaces, *nbEND = nbEN + nbElems;
2013 for ( ; nbEN < nbEND; ++nbEN, ++nbEF ) {
2014 nbNodes += *nbEN = accumulate( nbFN, nbFN + *nbEF, 0 );
2018 else // standard types
2020 nbElemNodesVec.myValues = new int[ 1 ];
2021 nbElemNodesVec[ 0 ] = nbElemNodes;
2022 nbNodes = nbElems * nbElemNodes;
2024 // iterator returning nbElemNodes for standard elems and
2025 // next value from nbElemNodesVec for poly elements
2026 _ValueIterator<int> nbElemNodesIt( nbElemNodesVec, isPoly ? 1 : 0);
2028 // iterator returning values form partType._medIndex for standard elems
2029 // and node index (n) for poly elements
2031 _ValueIterator<int> medIndexIt( isPoly ? &n : &partType._medIndex[0],
2033 // read connectivity
2034 _maille ma( partType._medType, nbElemNodes );
2035 ma.sommets.resize( nbElemNodes );
2036 TIntOwner connectivity( geoFile.getInt( nbNodes )); // nn_ei
2037 int* nodeID = connectivity;
2039 for ( int i = 0; i < nbElems; ++i ) {
2040 _ValueIterator<int> medIndex = medIndexIt;
2041 nbElemNodes = nbElemNodesIt.next();
2042 if ((int) ma.sommets.size() != nbElemNodes )
2043 ma.sommets.resize( nbElemNodes );
2044 for ( n = 0; n < nbElemNodes; ++n, ++nodeID ) {
2046 node = points.find( *nodeID + nodeShift );
2048 node = points.insert( make_pair( *nodeID + nodeShift, _noeud())).first;
2049 ma.sommets[ medIndex.next() ] = node;
2052 ma.setOrdre( elemIds[ i ] );
2053 groupe.mailles[i] = imed.insert(ma);
2055 // store nb nodes in polyhedron faces
2056 if ( nbFaceNodes.myValues ) {
2057 const int* nbFaceNodesPtr = nbFaceNodes.myValues;
2058 for ( int i = 0; i < nbElems; ++i ) {
2059 vector<int> & nbNodesByFace = imed.polyherdalNbFaceNodes[ &groupe.maille( i ) ];
2060 nbNodesByFace.assign( nbFaceNodesPtr, nbFaceNodesPtr + nbElemFaces[ i ] );
2061 nbFaceNodesPtr += nbElemFaces[ i ];
2064 // create subPart for "coordinates"
2065 if ( !haveCoords ) {
2066 _SubPart & coordSubPart = imed._subPartDescribed[ coordDesc ];
2067 coordSubPart.myFirstNode = points.insert( make_pair( 1 + nodeShift, _noeud())).first;
2069 // add subPart group to part group
2070 int groupeIndex = imed.groupes.size();
2071 partGroupe->groupes.push_back( groupeIndex );
2074 _SubPart subPart( partNum, partType._name );
2075 subPart.myNbCells = nbElems;
2076 subPart.myCellGroupIndex = groupeIndex;
2077 imed.addSubPart( subPart );
2079 else if ( word == "coordinates" )
2081 // Local node coordinates of a part
2082 // ------------------------------------
2083 int nbNodes = *TIntOwner( geoFile.getInt(1) ); // nn
2086 TIntOwner nodeIds(0);
2088 nodeIds.myValues = geoFile.getInt( nbNodes ); // id_n
2090 // find out if "coordinates" has already been add at reading connectivity
2091 _SubPartDesc coordDesc( partNum , "coordinates");
2092 map< _SubPartDesc, _SubPart >::iterator descPart =
2093 imed._subPartDescribed.find( coordDesc );
2094 bool haveCoords = ( descPart != imed._subPartDescribed.end() );
2097 // check that all nodes have been added
2098 firstNode = descPart->second.myFirstNode;
2099 descPart->second.myNbNodes = nbNodes;
2100 INoeud inoeud = firstNode, inoEnd = points.end();
2101 int id = inoeud->first, idEnd = id + nbNodes;
2102 for ( ; id < idEnd; ++id ) {
2103 if ( inoeud == inoEnd || inoeud->first > id ) {
2104 INoeud in = points.insert( inoeud, make_pair( id, _noeud() ));
2105 in->second.number = id;
2113 int nodeShift = points.empty() ? 0 : points.rbegin()->first;
2114 for ( int iNode = 1; iNode <= nbNodes; ++iNode ) {
2115 INoeud inoeud = points.insert( points.end(), make_pair( iNode + nodeShift, _noeud()));
2116 inoeud->second.number = inoeud->first;
2118 firstNode = points.find( 1 + nodeShift );
2119 // create "coordinates" subPart
2120 _SubPart & subPart = imed._subPartDescribed[ coordDesc ];
2121 subPart.myNbNodes = nbNodes;
2122 subPart.myFirstNode = firstNode;
2125 // read coordinates in no interlace mode
2126 TFltOwner noInterlaceCoords( geoFile.getFlt( nbNodes * SPACE_DIM ));
2127 float* x = noInterlaceCoords;
2128 float* y = x + nbNodes;
2129 float* z = y + nbNodes;
2130 INoeud endNode = points.end();
2131 for ( INoeud in = firstNode; in != endNode; ++in ) {
2132 _noeud & node = in->second;
2133 node.coord.resize( SPACE_DIM );
2134 node.coord[ 0 ] = *x++;
2135 node.coord[ 1 ] = *y++;
2136 node.coord[ 2 ] = *z++;
2139 else if ( word == "part" )
2141 // Another part encounters
2142 // -----------------------
2143 partNum = *TIntOwner( geoFile.getInt(1) );
2146 string partName( TStrOwner( geoFile.getLine() ));
2147 if ( partName.empty() )
2148 partName = "Part_" + restLine;
2150 if (int( imed.groupes.capacity() - imed.groupes.size()) < theMaxNbTypes )
2151 imed.groupes.reserve( size_t( 1.5 * imed.groupes.size() ));
2152 imed.groupes.push_back(_groupe());
2153 partGroupe = & imed.groupes.back();
2154 partGroupe->nom = partName;
2155 partGroupe->groupes.reserve( theMaxNbTypes );
2157 else if ( word == "block" )
2160 // ------------------
2161 bool rectilinear = ( restLine.find( "rectilinear" ) != restLine.npos );
2162 bool uniform = ( restLine.find( "uniform" ) != restLine.npos );
2163 bool curvilinear = ( !rectilinear && !uniform );
2164 bool iblanked = ( restLine.find( "iblanked" ) != restLine.npos );
2165 bool with_ghost = ( restLine.find( "with_ghost" ) != restLine.npos );
2166 bool range = ( restLine.find( "range" ) != restLine.npos );
2169 TIntOwner ijk( geoFile.getInt(3) );
2173 int NumberOfNodes = I*J*K;
2174 if ( !NumberOfNodes ) continue;
2178 TIntOwner ijkRange( geoFile.getInt( 6 ));// imin imax jmin jmax kmin kmax
2179 I = ijkRange[1]-ijkRange[0]+1;
2180 J = ijkRange[3]-ijkRange[2]+1;
2181 K = ijkRange[5]-ijkRange[4]+1;
2182 NumberOfNodes = I*J*K;
2185 int nodeShift = points.empty() ? 0 : points.rbegin()->first;
2186 for ( int iNode = 1; iNode <= NumberOfNodes; ++iNode ) {
2187 INoeud inoeud = points.insert( points.end(), make_pair( iNode + nodeShift, _noeud()));
2188 _noeud & node = inoeud->second;
2189 node.number = inoeud->first;
2190 node.coord.resize( SPACE_DIM );
2192 INoeud firstNode = points.find( nodeShift + 1 );
2193 INoeud endNode = points.end();
2195 GRID grid; // calculator of unstructured data
2196 grid._iArrayLength = I;
2197 grid._jArrayLength = J;
2198 grid._kArrayLength = K;
2199 grid._spaceDimension = SPACE_DIM;
2200 if ( J < 2 ) { grid._spaceDimension--; grid._jArrayLength = 0; }
2201 if ( K < 2 ) { grid._spaceDimension--; grid._kArrayLength = 0; }
2202 int nbElems = grid.getNumberOfElements(MED_CELL, MED_ALL_ELEMENTS);
2204 if ( curvilinear ) // read coordinates for all nodes
2206 TFltOwner noInterlaceCoords( geoFile.getFlt( NumberOfNodes * SPACE_DIM ));
2207 float* x = noInterlaceCoords;
2208 float* y = x + NumberOfNodes;
2209 float* z = y + NumberOfNodes;
2210 for ( INoeud in = firstNode; in != endNode; ++in ) {
2211 _noeud & node = in->second;
2212 node.coord.resize( SPACE_DIM );
2213 node.coord[ 0 ] = *x++;
2214 node.coord[ 1 ] = *y++;
2215 node.coord[ 2 ] = *z++;
2217 grid._gridType = MED_BODY_FITTED;
2218 grid._coordinate = new COORDINATE(3, 1, 0); // to avoid exception in convertInMESH()
2220 else if ( rectilinear ) // read delta vectors with non-regular spacing
2222 grid._iArray = (double*)geoFile.convertReals<double>( I );
2223 grid._jArray = (double*)geoFile.convertReals<double>( J );
2224 grid._kArray = (double*)geoFile.convertReals<double>( K );
2225 grid._gridType = MED_CARTESIAN;
2227 else // uniform: read grid origine and delta vectors for regular spacing grid
2229 TFltOwner xyzOrigin( geoFile.getFlt( 3 ));
2230 TFltOwner xyzDelta ( geoFile.getFlt( 3 ));
2231 // compute full delta vectors
2232 grid._iArray = new double[ I ];
2233 grid._jArray = new double[ J ];
2234 grid._kArray = new double[ K ];
2235 double* coors[SPACE_DIM] = { grid._iArray, grid._jArray, grid._kArray };
2236 int size [SPACE_DIM] = { I, J, K };
2237 for ( int j = 0; j < SPACE_DIM; ++j ) {
2238 double* coo = coors[ j ];
2239 double* cooEnd = coo + size[ j ];
2240 coo[0] = xyzOrigin[ j ];
2241 while ( ++coo < cooEnd )
2242 *coo = coo[-1] + xyzDelta[ j ];
2244 grid._gridType = MED_CARTESIAN;
2249 geoFile.skip( NumberOfNodes * sizeof(int) );
2252 TStrOwner( geoFile.getLine() ); // "ghost_flags"
2253 geoFile.skip( nbElems * sizeof(int) );
2256 if ( haveNodeIds && !geoFile.eof() ) {
2257 TStrOwner nextLine( geoFile.getLine() ); // "node_ids"
2258 if ( contains( "node_ids", nextLine ) )
2259 geoFile.skip( NumberOfNodes * sizeof(int) );
2261 geoFile.skip( -MAX_LINE_LENGTH );
2264 TIntOwner elemIdOwner(0);
2265 _ValueIterator<int> elemIds;
2266 if ( haveElemIds && !geoFile.eof() ) {
2267 TStrOwner nextLine( geoFile.getLine() ); // "element_ids"
2268 if ( contains( "element_ids", nextLine ) ) {
2269 elemIdOwner.myValues = geoFile.getInt( nbElems );
2270 elemIds = _ValueIterator<int>( elemIdOwner, 1);
2272 geoFile.skip( -MAX_LINE_LENGTH );
2276 // let GRID compute all coordinates and connectivity
2277 const MESH* unstruct = grid.convertInMESH();
2280 const double * coo = unstruct->getCoordinates(MED_FULL_INTERLACE);
2281 typedef _ValueIterator< double > TCoordIt;
2282 TCoordIt xCoo( coo+0, grid._spaceDimension);
2283 TCoordIt yCoo( coo+1, grid._spaceDimension);
2284 TCoordIt zCoo( coo+2, grid._spaceDimension);
2285 if ( grid._spaceDimension < 3 ) zCoo = TCoordIt( grid._kArray, 0 );
2286 if ( grid._spaceDimension < 2 ) yCoo = TCoordIt( grid._jArray, 0 );
2287 for ( INoeud in = firstNode; in != endNode; ++in ) {
2288 _noeud& node = in->second;
2289 node.coord[ 0 ] = xCoo.next();
2290 node.coord[ 1 ] = yCoo.next();
2291 node.coord[ 2 ] = zCoo.next();
2295 // store connectivity
2296 const int * conn = unstruct->getConnectivity( MED_NODAL, MED_CELL, MED_ALL_ELEMENTS );
2297 medGeometryElement elemType = grid.getElementType( MED_CELL, 1 );
2298 int nbElemNodes = elemType % 100;
2300 partGroupe->mailles.resize( nbElems );
2301 _maille ma( elemType, nbElemNodes );
2302 ma.sommets.resize( nbElemNodes );
2304 for ( int i = 0, nIndex = 0; i < nbElems; ++i ) {
2305 for ( int n = 0; n < nbElemNodes; ++n ) {
2306 int nodeID = conn[ nIndex++ ];
2307 ma.sommets[n] = points.find( nodeID + nodeShift );
2309 ma.setOrdre( elemIds.next() );
2310 partGroupe->mailles[i] = imed.insert(ma);
2313 _SubPart subPart( partNum, "block" );
2314 subPart.myNbCells = nbElems;
2315 subPart.myNbNodes = NumberOfNodes;
2316 subPart.myFirstNode = firstNode;
2317 subPart.myCellGroupIndex = imed.groupes.size();
2318 imed.addSubPart( subPart );
2320 unstruct->removeReference();
2322 else if ( word == "extents" )
2324 extents.myValues = geoFile.getFlt( 6 );
2329 ( LOCALIZED( STRING(LOC) << "Unexpected word: " << word <<
2330 " in " << getDataFileName()));
2332 } // while ( !geoFile.eof() )
2335 imed.mergeNodesAndElements(TOLERANCE);
2340 //================================================================================
2342 * \brief Read mesh in Ensight6 ASCII format
2344 //================================================================================
2346 void ENSIGHT_MESH_RDONLY_DRIVER::read6ASCII(_InterMed & imed)
2348 const char * LOC ="ENSIGHT_MESH_RDONLY_DRIVER::read6ASCII() : ";
2351 _ASCIIFileReader geoFile( getDataFileName() );
2353 if ( isSingleFileMode() ) {
2354 int curTimeStep = 1;
2355 while ( curTimeStep < getIndexInDataFile() ) {
2356 while ( !geoFile.isTimeStepEnd())
2360 while ( !geoFile.isTimeStepBeginning() )
2363 // ----------------------
2364 // Read mesh description
2365 // ----------------------
2367 string descriptionLine1 = geoFile.getLine();
2368 string descriptionLine2 = geoFile.getLine();
2370 // find out if the file was created by MED driver
2371 int prefixSize = strlen( theDescriptionPrefix );
2372 _isMadeByMed = ( descriptionLine1.substr(0, prefixSize ) == theDescriptionPrefix );
2375 descriptionLine1 = descriptionLine1.substr( prefixSize );
2376 _ptrMesh->setDescription( descriptionLine1 + descriptionLine2 );
2379 // ----------------------------------------
2380 // Find out presence of node/elem numbers
2381 // ----------------------------------------
2383 // EnSight User Manual (for v8) says:
2384 // You do not have to assign node IDs. If you do, the element connectivities are
2385 // based on the node numbers. If you let EnSight assign the node IDs, the nodes
2386 // are considered to be sequential starting at node 1, and element connectivity is
2387 // done accordingly. If node IDs are set to off, they are numbered internally;
2388 // however, you will not be able to display or query on them. If you have node
2389 // IDs in your data, you can have EnSight ignore them by specifying "node id
2390 // ignore." Using this option may reduce some of the memory taken up by the
2391 // Client and Server, but display and query on the nodes will not be available.
2393 // read "node|element id <off|given|assign|ignore>"
2394 geoFile.getWord(); geoFile.getWord();
2395 string nodeIds = geoFile.getWord();
2396 geoFile.getWord(); geoFile.getWord();
2397 string elemIds = geoFile.getWord();
2399 bool haveNodeIds = ( nodeIds == "given" || nodeIds == "ignore" );
2400 bool haveElemIds = ( elemIds == "given" || elemIds == "ignore" );
2402 map<int,_noeud> & points = imed.points;
2403 typedef map<int,_noeud>::iterator INoeud;
2405 int haveStructuredParts = 0, haveUnstructuredParts = 0;
2407 _groupe * partGroupe = 0;
2410 while ( !geoFile.isTimeStepEnd() )
2412 string word, restLine, line = geoFile.getLine();
2413 TStrTool::split( line, word, restLine );
2415 const TEnSightElemType & partType = getEnSightType( word );
2416 if ( !partType._medIndex.empty() )
2418 // Unstructured element type encounters
2419 // --------------------------------------
2420 int nbElemNodes = partType._medType % 100;
2421 int nbElems = geoFile.getInt();
2423 haveUnstructuredParts++;
2425 imed.groupes.push_back(_groupe());
2426 _groupe & groupe = imed.groupes.back();
2427 groupe.mailles.resize( nbElems );
2429 // read connectivity
2430 _maille ma( partType._medType, nbElemNodes );
2431 ma.sommets.resize( nbElemNodes );
2433 for ( int i = 0; i < nbElems; ++i ) {
2436 for ( int n = 0; n < nbElemNodes; ++n ) {
2437 int nodeID = geoFile.getInt();
2438 ma.sommets[ partType._medIndex[n] ] = points.find( nodeID );
2440 //ma.ordre = ++order;
2441 groupe.mailles[i] = imed.insert(ma);
2444 int groupeIndex = imed.groupes.size();
2445 partGroupe->groupes.push_back( groupeIndex );
2447 _SubPart subPart( partNum, partType._name );
2448 subPart.myNbCells = nbElems;
2449 subPart.myCellGroupIndex = groupeIndex;
2450 imed.addSubPart( subPart );
2452 else if ( word == "part" )
2454 // Another part encounters
2455 // -----------------------
2456 partNum = atoi( restLine.c_str() );
2458 string partName = geoFile.getLine();
2459 if ( partName.empty() )
2460 partName = "Part_" + restLine;
2462 if (int( imed.groupes.capacity() - imed.groupes.size()) < theMaxNbTypes )
2463 imed.groupes.reserve( size_t( 1.5 * imed.groupes.size() ));
2464 imed.groupes.push_back(_groupe());
2465 partGroupe = & imed.groupes.back();
2466 partGroupe->nom = partName;
2467 partGroupe->groupes.reserve( theMaxNbTypes );
2469 else if ( word == "block" )
2472 // ------------------
2473 bool iblanked = ( restLine == "iblanked" );
2476 int I = geoFile.getInt();
2477 int J = geoFile.getInt();
2478 int K = geoFile.getInt();
2479 int NumberOfNodes = I*J*K;
2480 if ( !NumberOfNodes ) continue;
2481 haveStructuredParts++;
2484 int nodeShift = points.empty() ? 0 : points.rbegin()->first;
2485 for ( int iNode = 1; iNode <= NumberOfNodes; ++iNode ) {
2486 INoeud inoeud = points.insert( points.end(), make_pair( iNode + nodeShift, _noeud()));
2487 _noeud & node = inoeud->second;
2488 node.number = inoeud->first;
2489 node.coord.resize( SPACE_DIM );
2492 INoeud firstNode = points.find( nodeShift + 1 );
2493 INoeud endNode = points.end();
2494 for ( int j = 0; j < SPACE_DIM; ++j ) {
2495 for ( INoeud in = firstNode; in != endNode; ++in ) {
2496 _noeud & node = in->second;
2497 node.coord[ j ] = geoFile.getReal();
2502 geoFile.skip(NumberOfNodes, /*nbPerLine =*/ 10, INT_WIDTH_6);
2504 // let GRID calculate connectivity
2506 grid._iArrayLength = I;
2507 grid._jArrayLength = J;
2508 grid._kArrayLength = K;
2509 grid._gridType = MED_BODY_FITTED;
2510 grid._spaceDimension= SPACE_DIM;
2511 grid._coordinate = new COORDINATE(3, 1, 0); // to avoid exception in convertInMESH()
2512 if ( J < 2 ) { grid._spaceDimension--; grid._jArrayLength = 0; }
2513 if ( K < 2 ) { grid._spaceDimension--; grid._kArrayLength = 0; }
2515 const MESH* unstruct = grid.convertInMESH();
2516 const int * conn = unstruct->getConnectivity( MED_NODAL, MED_CELL, MED_ALL_ELEMENTS );
2517 medGeometryElement elemType = grid.getElementType( MED_CELL, 1 );
2518 int nbElemNodes = elemType % 100;
2519 int nbElems = grid.getNumberOfElements(MED_CELL, MED_ALL_ELEMENTS);
2521 partGroupe->mailles.resize( nbElems );
2522 _maille ma( elemType, nbElemNodes );
2523 ma.sommets.resize( nbElemNodes );
2525 for ( int i = 0, nIndex = 0; i < nbElems; ++i ) {
2526 for ( int n = 0; n < nbElemNodes; ++n ) {
2527 int nodeID = conn[ nIndex++ ];
2528 ma.sommets[n] = points.find( nodeID + nodeShift );
2530 //ma.ordre = ++order;
2531 partGroupe->mailles[i] = imed.insert(ma);
2534 _SubPart subPart( partNum, "block" );
2535 subPart.myNbCells = nbElems;
2536 subPart.myNbNodes = NumberOfNodes;
2537 subPart.myFirstNode = firstNode;
2538 subPart.myCellGroupIndex = imed.groupes.size();
2539 imed.addSubPart( subPart );
2541 unstruct->removeReference();
2543 else if ( word == "coordinates" )
2545 // ------------------------------------
2546 // Unstructured global node coordinates
2547 // ------------------------------------
2548 int nbNodes = geoFile.getInt();
2550 cout << "-> loading coordinates of " << nbNodes << " nodes " << endl ;
2553 for ( int i=0 ; i < nbNodes ; i++ )
2555 if ( haveNodeIds ) {
2556 int nodeID = geoFile.getInt();
2557 inoeud = points.insert( make_pair( nodeID, _noeud() )).first;
2558 inoeud->second.number = nodeID;
2562 inoeud = points.insert( points.end(), make_pair( nodeID, _noeud()));
2563 inoeud->second.number = nodeID;
2565 _noeud & node = inoeud->second;
2566 node.coord.resize( SPACE_DIM );
2567 node.coord[ 0 ] = geoFile.getReal();
2568 node.coord[ 1 ] = geoFile.getReal();
2569 node.coord[ 2 ] = geoFile.getReal();
2572 _SubPartDesc cooDesc = _SubPartDesc::globalCoordDesc();
2573 _SubPart subPart( cooDesc.partNumber(), cooDesc.typeName() );
2574 subPart.myNbNodes = nbNodes;
2575 subPart.myFirstNode = points.begin();
2576 imed.addSubPart( subPart );
2581 ( LOCALIZED( STRING(LOC) << "Unexpected word: " << word <<
2582 " in " << getDataFileName()));
2584 } // while ( !geoFile.eof() )
2586 if ( ( haveStructuredParts && haveUnstructuredParts ) || haveStructuredParts > 1 )
2587 imed.mergeNodesAndElements(TOLERANCE);
2592 //================================================================================
2594 * \brief Read mesh in Ensight6 ASCII format
2596 //================================================================================
2598 void ENSIGHT_MESH_RDONLY_DRIVER::read6Binary(_InterMed & imed)
2600 const char * LOC ="ENSIGHT_MESH_RDONLY_DRIVER::read6Binary() : ";
2603 _BinaryFileReader geoFile( getDataFileName() );
2605 // check if swapping bytes needed
2607 countPartsBinary( geoFile, isSingleFileMode());
2610 geoFile.swapBytes();
2613 if ( getIndexInDataFile() <= 1 )
2615 if ( geoFile.getPosition() == 0 ) {
2616 TStrOwner format( geoFile.getLine() ); // "C|Fortran Binary"
2617 if ( !contains( "C Binary", format )) {
2618 if ( contains( "Fortran Binary", format ))
2619 throw(MEDEXCEPTION(STRING(LOC) << "Fortran Binary format not supported"));
2621 throw(MEDEXCEPTION(STRING(LOC) << "unexpected line in " << getDataFileName()
2622 << "\n" << format.myValues));
2626 if ( isSingleFileMode() ) {
2627 // one time step may be skipped by countPartsBinary
2628 int curTimeStep = geoFile.getPosition() ? 2 : 1 ;
2629 while ( curTimeStep < getIndexInDataFile() ) {
2630 countPartsBinary( geoFile, true ); // skip time step
2634 TStrOwner line( geoFile.getLine() );
2635 if ( isTimeStepBeginning( line.myValues ))
2639 // ----------------------
2640 // Read mesh description
2641 // ----------------------
2643 TStrOwner descriptionLine1( geoFile.getLine() );
2644 TStrOwner descriptionLine2( geoFile.getLine() );
2646 // find out if the file was created by MED driver
2647 _isMadeByMed = contains( theDescriptionPrefix, descriptionLine1 );
2650 _ptrMesh->setDescription( descriptionLine2.myValues );
2652 _ptrMesh->setDescription( string(descriptionLine1) + descriptionLine2.myValues );
2655 // ----------------------------------------
2656 // Find out presence of node/elem numbers
2657 // ----------------------------------------
2659 // EnSight User Manual (for v8) says:
2660 // You do not have to assign node IDs. If you do, the element connectivities are
2661 // based on the node numbers. If you let EnSight assign the node IDs, the nodes
2662 // are considered to be sequential starting at node 1, and element connectivity is
2663 // done accordingly. If node IDs are set to off, they are numbered internally;
2664 // however, you will not be able to display or query on them. If you have node
2665 // IDs in your data, you can have EnSight ignore them by specifying "node id
2666 // ignore." Using this option may reduce some of the memory taken up by the
2667 // Client and Server, but display and query on the nodes will not be available.
2669 // read "node|element id <off|given|assign|ignore>"
2670 bool haveNodeIds, haveElemIds;
2672 TStrOwner nodeIds( geoFile.getLine() );
2673 TStrOwner elemIds( geoFile.getLine() );
2675 haveNodeIds = ( contains( "given", nodeIds ) || contains( "ignore", nodeIds ) );
2676 haveElemIds = ( contains( "given", elemIds ) || contains( "ignore", elemIds ) );
2678 map<int,_noeud> & points = imed.points;
2679 typedef map<int,_noeud>::iterator INoeud;
2681 int haveStructuredParts = 0, haveUnstructuredParts = 0;
2683 _groupe * partGroupe = 0;
2686 while ( !geoFile.eof() )
2688 TStrOwner line( geoFile.getLine() );
2689 if ( isSingleFileMode() && isTimeStepEnd( line.myValues ))
2691 string word, restLine;
2692 TStrTool::split( line.myValues, word, restLine );
2694 const TEnSightElemType & partType = getEnSightType( word );
2695 if ( !partType._medIndex.empty() )
2697 // Unstructured element type encounters
2698 // --------------------------------------
2699 int nbElemNodes = partType._medType % 100;
2700 int nbElems = *TIntOwner( geoFile.getInt(1) ); // ne
2702 haveUnstructuredParts++;
2704 TIntOwner numbers(0);
2706 numbers.myValues = geoFile.getInt( nbElems ); // id_e
2708 imed.groupes.push_back(_groupe());
2709 _groupe & groupe = imed.groupes.back();
2710 groupe.mailles.resize( nbElems );
2712 // read connectivity
2713 _maille ma( partType._medType, nbElemNodes );
2714 ma.sommets.resize( nbElemNodes );
2715 TIntOwner connectivity( geoFile.getInt( nbElems * nbElemNodes ));
2716 int* nodeID = connectivity;
2718 for ( int i = 0; i < nbElems; ++i ) {
2719 for ( int n = 0; n < nbElemNodes; ++n, ++nodeID )
2720 ma.sommets[ partType._medIndex[n] ] = points.find( *nodeID );
2721 //ma.ordre = ++order;
2722 groupe.mailles[i] = imed.insert(ma);
2725 int groupeIndex = imed.groupes.size();
2726 partGroupe->groupes.push_back( groupeIndex );
2728 _SubPart subPart( partNum, partType._name );
2729 subPart.myNbCells = nbElems;
2730 subPart.myCellGroupIndex = groupeIndex;
2731 imed.addSubPart( subPart );
2733 else if ( word == "part" )
2735 // Another part encounters
2736 // -----------------------
2737 partNum = atoi( restLine.c_str() );
2739 string partName( TStrOwner( geoFile.getLine() ));
2740 if ( partName.empty() )
2741 partName = "Part_" + restLine;
2743 if (int( imed.groupes.capacity() - imed.groupes.size()) < theMaxNbTypes )
2744 imed.groupes.reserve( size_t( 1.5 * imed.groupes.size() ));
2745 imed.groupes.push_back(_groupe());
2746 partGroupe = & imed.groupes.back();
2747 partGroupe->nom = partName;
2748 partGroupe->groupes.reserve( theMaxNbTypes );
2750 else if ( word == "block" )
2753 // ------------------
2754 bool iblanked = ( restLine == "iblanked" );
2757 TIntOwner ijk( geoFile.getInt(3) );
2761 int NumberOfNodes = I*J*K;
2762 if ( !NumberOfNodes ) continue;
2763 haveStructuredParts++;
2766 int nodeShift = points.empty() ? 0 : points.rbegin()->first;
2768 TFltOwner noInterlaceCoords( geoFile.getFlt( NumberOfNodes * SPACE_DIM ));
2769 float* x = noInterlaceCoords;
2770 float* y = x + NumberOfNodes;
2771 float* z = y + NumberOfNodes;
2772 for ( int iNode = 1; iNode <= NumberOfNodes; ++iNode ) {
2773 INoeud inoeud = points.insert( points.end(), make_pair( iNode + nodeShift, _noeud()));
2774 _noeud & node = inoeud->second;
2775 node.number = inoeud->first;
2776 node.coord.resize( SPACE_DIM );
2777 node.coord[0] = *x++;
2778 node.coord[1] = *y++;
2779 node.coord[2] = *z++;
2784 geoFile.skip(NumberOfNodes * sizeof(int));
2786 // let GRID calculate connectivity
2788 grid._iArrayLength = I;
2789 grid._jArrayLength = J;
2790 grid._kArrayLength = K;
2791 grid._gridType = MED_BODY_FITTED;
2792 grid._spaceDimension= SPACE_DIM;
2793 grid._coordinate = new COORDINATE(3, 1, 0); // to avoid exception in convertInMESH()
2794 if ( J < 2 ) { grid._spaceDimension--; grid._jArrayLength = 0; }
2795 if ( K < 2 ) { grid._spaceDimension--; grid._kArrayLength = 0; }
2797 const MESH* unstruct = grid.convertInMESH();
2798 const int * conn = unstruct->getConnectivity( MED_NODAL, MED_CELL, MED_ALL_ELEMENTS );
2799 medGeometryElement elemType = grid.getElementType( MED_CELL, 1 );
2800 int nbElemNodes = elemType % 100;
2801 int nbElems = grid.getNumberOfElements(MED_CELL, MED_ALL_ELEMENTS);
2803 partGroupe->mailles.resize( nbElems );
2804 _maille ma( elemType, nbElemNodes );
2805 ma.sommets.resize( nbElemNodes );
2807 for ( int i = 0, nIndex = 0; i < nbElems; ++i ) {
2808 for ( int n = 0; n < nbElemNodes; ++n ) {
2809 int nodeID = conn[ nIndex++ ];
2810 ma.sommets[n] = points.find( nodeID + nodeShift );
2812 //ma.ordre = ++order;
2813 partGroupe->mailles[i] = imed.insert(ma);
2816 _SubPart subPart( partNum, "block" );
2817 subPart.myNbCells = nbElems;
2818 subPart.myNbNodes = NumberOfNodes;
2819 subPart.myFirstNode = points.find( nodeShift + 1 );
2820 subPart.myCellGroupIndex = imed.groupes.size();
2821 imed.addSubPart( subPart );
2823 unstruct->removeReference();
2825 else if ( word == "coordinates" )
2827 // ------------------------------
2828 // Unstructured node coordinates
2829 // ------------------------------
2830 int nbNodes = *TIntOwner( geoFile.getInt(1) );
2832 TIntOwner numbers(0);
2834 numbers.myValues = geoFile.getInt( nbNodes );
2836 TFltOwner fullInterlaceCoords( geoFile.getFlt( nbNodes * SPACE_DIM ));
2837 float* coord = fullInterlaceCoords;
2839 cout << "-> loading coordinates of " << nbNodes << " nodes " << endl ;
2842 for ( int i=0 ; i < nbNodes ; i++ )
2844 if ( haveNodeIds ) {
2845 int nodeID = numbers[ i ];
2846 inoeud = points.insert( make_pair( nodeID, _noeud() )).first;
2847 inoeud->second.number = nodeID;
2851 inoeud = points.insert( points.end(), make_pair( nodeID, _noeud()));
2852 inoeud->second.number = nodeID;
2854 _noeud & node = inoeud->second;
2855 node.coord.resize( SPACE_DIM );
2856 node.coord[ 0 ] = *coord++;
2857 node.coord[ 1 ] = *coord++;
2858 node.coord[ 2 ] = *coord++;
2861 _SubPartDesc cooDesc = _SubPartDesc::globalCoordDesc();
2862 _SubPart subPart( cooDesc.partNumber(), cooDesc.typeName() );
2863 subPart.myNbNodes = nbNodes;
2864 subPart.myFirstNode = points.begin();
2865 imed.addSubPart( subPart );
2870 ( LOCALIZED( STRING(LOC) << "Unexpected word: " << word <<
2871 " in " << getDataFileName()));
2873 } // while ( !geoFile.eof() )
2875 if ( ( haveStructuredParts && haveUnstructuredParts ) || haveStructuredParts > 1 )
2876 imed.mergeNodesAndElements(TOLERANCE);
2881 //================================================================================
2883 * \brief count number of parts in EnSight geometry file
2885 //================================================================================
2887 int ENSIGHT_MESH_RDONLY_DRIVER::countParts(const string& geomFileName,
2888 const bool isSingleFileMode)
2890 const char * LOC ="ENSIGHT_MESH_RDONLY_DRIVER::countParts() : ";
2893 if ( isBinaryDataFile( geomFileName ))
2895 _BinaryFileReader geoFile(geomFileName);
2896 // check if swapping bytes needed
2898 return countPartsBinary( geoFile, isSingleFileMode );
2902 geoFile.swapBytes();
2904 nbParts = countPartsBinary( geoFile, isSingleFileMode );
2908 _ASCIIFileReader geoFile(geomFileName);
2910 if ( isSingleFileMode )
2911 while ( !isTimeStepBeginning( geoFile.getLine() ));
2913 geoFile.getLine(); // description line 1
2914 geoFile.getLine(); // description line 2
2916 // read "node|element id <off|given|assign|ignore>"
2917 geoFile.getWord(); geoFile.getWord();
2918 string nodeIds = geoFile.getWord();
2919 geoFile.getWord(); geoFile.getWord();
2920 string elemIds = geoFile.getWord();
2921 bool haveNodeIds = ( nodeIds == "given" || nodeIds == "ignore" );
2922 bool haveElemIds = ( elemIds == "given" || elemIds == "ignore" );
2925 while ( !geoFile.isTimeStepEnd() )
2927 string word, restLine, line = geoFile.getLine();
2928 TStrTool::split( line, word, restLine );
2930 const TEnSightElemType & partType = getEnSightType( word );
2931 if ( partType._medType != MED_ALL_ELEMENTS )
2933 // Unstructured element type encounters
2934 // --------------------------------------
2935 int nbElemNodes = partType._medType % 100;
2936 int nbElems = geoFile.getInt(); // ne
2939 if ( haveElemIds && isGold )
2940 geoFile.skip( nbElems, /*nbPerLine =*/ 1, INT_WIDTH_GOLD ); // id_e
2942 // skip connectivity
2943 int nbNodes = nbElems * nbElemNodes;
2944 if ( partType._name == "nsided" ) // polygons
2946 for ( int i = 0; i < nbElems; ++i )
2947 nbNodes += geoFile.getInt();
2948 geoFile.skip( nbNodes * INT_WIDTH_GOLD, /*nbLines = */nbElems );
2950 else if ( partType._name == "nfaced" ) // polyhedrons
2953 for ( int i = 0; i < nbElems; ++i )
2954 nbFaces += geoFile.getInt();
2955 for ( int f = 0; f < nbFaces; ++f )
2956 nbNodes += geoFile.getInt();
2957 geoFile.skip( nbNodes * INT_WIDTH_GOLD, /*nbLines = */nbFaces );
2959 else // standard types
2962 geoFile.skip( nbNodes, nbElemNodes, INT_WIDTH_GOLD );
2963 else if ( haveElemIds )
2964 geoFile.skip( nbNodes + nbElems, nbElemNodes+1, INT_WIDTH_6 );
2966 geoFile.skip( nbNodes, nbElemNodes, INT_WIDTH_6 );
2969 else if ( word == "coordinates" )
2971 isGold = ( nbParts > 0 );
2972 int nbNodes = geoFile.getInt(); // nn
2977 geoFile.skip( nbNodes, /*nbPerLine =*/ 1, INT_WIDTH_GOLD ); // node ids
2978 geoFile.skip( nbNodes * SPACE_DIM, /*nbPerLine =*/ 1, FLT_WIDTH ); // coordinates
2981 int coordWidth = 3 * FLT_WIDTH;
2983 coordWidth += INT_WIDTH_6;
2984 geoFile.skip(nbNodes, /*nbPerLine =*/ 1, coordWidth);
2987 else if ( word == "part" )
2991 geoFile.skip( 1, /*nbPerLine =*/ 1, INT_WIDTH_GOLD ); //part number
2993 geoFile.getWord(); // part number
2994 geoFile.toNextLine();
2995 geoFile.getLine(); // description line
2997 else if ( word == "block" )
3000 // ------------------
3001 bool rectilinear = ( restLine.find( "rectilinear" ) != restLine.npos );
3002 bool uniform = ( restLine.find( "uniform" ) != restLine.npos );
3003 bool curvilinear = ( !rectilinear && !uniform );
3004 bool iblanked = ( restLine.find( "iblanked" ) != restLine.npos );
3005 bool with_ghost = ( restLine.find( "with_ghost" ) != restLine.npos );
3006 bool range = ( restLine.find( "range" ) != restLine.npos );
3009 int I = geoFile.getInt();
3010 int J = geoFile.getInt();
3011 int K = geoFile.getInt();
3012 int nbNodes = I*J*K;
3013 if ( !nbNodes ) continue;
3017 vector<int> ijkRange; // imin imax jmin jmax kmin kmax
3018 ijkRange.reserve(6);
3019 while ( ijkRange.size() < 6 )
3020 ijkRange.push_back( geoFile.getInt() );
3021 I = ijkRange[1]-ijkRange[0]+1;
3022 J = ijkRange[3]-ijkRange[2]+1;
3023 K = ijkRange[5]-ijkRange[4]+1;
3026 int nbElems = (I-1)*(J-1)*(K-1);
3028 if ( curvilinear ) // read coordinates for all nodes
3031 geoFile.skip( nbNodes, /*nbPerLine =*/ 1, FLT_WIDTH );
3033 geoFile.skip( nbNodes * SPACE_DIM, /*nbPerLine =*/ 6, FLT_WIDTH );
3035 else if ( rectilinear ) // read delta vectors with non-regular spacing
3037 geoFile.skip( I + J + K, /*nbPerLine =*/ 1, FLT_WIDTH );
3039 else // uniform: read grid origine and delta vectors for regular spacing grid
3041 geoFile.skip( 6, /*nbPerLine =*/ 1, FLT_WIDTH );
3047 geoFile.skip( nbNodes, /*nbPerLine =*/ 1, INT_WIDTH_GOLD);
3049 geoFile.skip( nbNodes, /*nbPerLine =*/ 10, INT_WIDTH_6 );
3053 geoFile.getWord(); // "ghost_flags"
3054 geoFile.skip( nbElems, /*nbPerLine =*/ 1, INT_WIDTH_GOLD);
3057 if ( haveNodeIds && geoFile.lookAt( "node_ids" )) {
3058 geoFile.getWord(); // "node_ids"
3059 geoFile.skip( nbNodes, /*nbPerLine =*/ 1, INT_WIDTH_GOLD);
3062 if ( haveElemIds && geoFile.lookAt( "element_ids" ) ) {
3063 geoFile.getWord(); // "element_ids"
3064 geoFile.skip( nbElems, /*nbPerLine =*/ 1, INT_WIDTH_GOLD);
3067 else if ( word == "extents" ) {
3068 geoFile.getLine(); geoFile.getLine(); geoFile.getLine();// 3 x 2E12.5
3072 throw MED_EXCEPTION ( LOCALIZED( STRING(LOC) << "Unexpected word: " << word));
3079 //================================================================================
3081 * \brief count number of parts in EnSight geometry file
3083 //================================================================================
3085 int ENSIGHT_MESH_RDONLY_DRIVER::countPartsBinary(_BinaryFileReader& geoFile,
3086 const bool isSingleFileMode)
3088 const char * LOC ="ENSIGHT_MESH_RDONLY_DRIVER::countPartsBinary() : ";
3090 if ( geoFile.getPosition() == 0 ) {
3091 TStrOwner format( geoFile.getLine() ); // "C|Fortran Binary"
3092 if ( !contains( "C Binary", format )) {
3093 if ( contains( "Fortran Binary", format ))
3094 throw(MEDEXCEPTION(STRING(LOC) << "Fortran Binary format not supported"));
3096 throw(MEDEXCEPTION(STRING(LOC) << "unexpected line: \n" << format.myValues));
3100 if ( isSingleFileMode ) {
3102 TStrOwner line( geoFile.getLine() );
3103 if ( isTimeStepBeginning( line.myValues ))
3108 // 2 description lines
3109 // ----------------------
3110 geoFile.skip( 2 * MAX_LINE_LENGTH );
3112 // read "node|element id <off|given|assign|ignore>"
3113 bool haveNodeIds, haveElemIds;
3115 TStrOwner nodeIds( geoFile.getLine() );
3116 TStrOwner elemIds( geoFile.getLine() );
3118 haveNodeIds = ( contains( "given", nodeIds ) || contains( "ignore", nodeIds ) );
3119 haveElemIds = ( contains( "given", elemIds ) || contains( "ignore", elemIds ) );
3122 int nbParts = 0; // the result
3125 while ( !geoFile.eof() )
3127 TStrOwner line( geoFile.getLine() );
3128 if ( isSingleFileMode && isTimeStepEnd( line.myValues ))
3130 string word, restLine;
3131 TStrTool::split( line.myValues, word, restLine );
3133 const TEnSightElemType & partType = getEnSightType( word );
3134 if ( partType._medType != MED_ALL_ELEMENTS )
3136 // Unstructured element type encounters
3137 // --------------------------------------
3138 int nbElems = *TIntOwner( geoFile.getInt(1) ); // ne
3139 int nbElemNodes = partType._medType % 100;
3143 geoFile.skip( nbElems * sizeof(int) ); // id_e
3145 int nbNodes = nbElems * nbElemNodes;
3146 if ( partType._name == "nsided" ) // polygons
3148 TIntOwner nbNodesInFace( geoFile.getInt( nbElems ));
3149 nbNodes = std::accumulate( nbNodesInFace.myValues, nbNodesInFace.myValues + nbElems, 0 );
3151 else if ( partType._name == "nfaced" ) // polyhedrons
3153 TIntOwner nbElemFaces( geoFile.getInt( nbElems ));
3154 int nbFaces = accumulate( nbElemFaces.myValues, nbElemFaces.myValues + nbElems, 0 );
3155 TIntOwner nbNodesInFace( geoFile.getInt( nbFaces ));
3156 nbNodes = std::accumulate( nbNodesInFace.myValues, nbNodesInFace.myValues + nbFaces, 0 );
3158 geoFile.skip( nbNodes * sizeof(int) );
3160 else if ( word == "coordinates" )
3164 int nbNodes = *TIntOwner( geoFile.getInt(1) ); // nn
3168 geoFile.skip( nbNodes * sizeof(int) ); // id_n
3171 geoFile.skip( nbNodes * SPACE_DIM * sizeof(int) );
3173 else if ( word == "part" )
3176 if ( isGold ) geoFile.skip(sizeof(int)); // part #
3178 geoFile.skip(MAX_LINE_LENGTH); // description line
3180 else if ( word == "block" )
3183 // ------------------
3184 bool rectilinear = ( restLine.find( "rectilinear" ) != restLine.npos );
3185 bool uniform = ( restLine.find( "uniform" ) != restLine.npos );
3186 bool curvilinear = ( !rectilinear && !uniform );
3187 bool iblanked = ( restLine.find( "iblanked" ) != restLine.npos );
3188 bool with_ghost = ( restLine.find( "with_ghost" ) != restLine.npos );
3189 bool range = ( restLine.find( "range" ) != restLine.npos );
3192 TIntOwner ijk( geoFile.getInt(3) );
3196 int NumberOfNodes = I*J*K;
3197 if ( !NumberOfNodes ) {
3198 if ( I != 0 && J != 0 && K != 0 )
3199 throw MEDEXCEPTION( "Need to swap bytes" );
3205 TIntOwner ijkRange( geoFile.getInt( 6 ));// imin imax jmin jmax kmin kmax
3206 I = ijkRange[1]-ijkRange[0]+1;
3207 J = ijkRange[3]-ijkRange[2]+1;
3208 K = ijkRange[5]-ijkRange[4]+1;
3209 NumberOfNodes = I*J*K;
3211 int nbElems = (I-1)*(J-1)*(K-1);
3213 if ( curvilinear ) // read coordinates for all nodes
3215 geoFile.skip( NumberOfNodes * SPACE_DIM * sizeof(float) );
3217 else if ( rectilinear ) // read delta vectors with non-regular spacing
3219 geoFile.skip( (I+J+K) * sizeof(float) );
3221 else // uniform: read grid origine and delta vectors for regular spacing grid
3223 geoFile.skip( 6 * sizeof(float) );
3228 geoFile.skip( NumberOfNodes * sizeof(int) );
3231 geoFile.skip( MAX_LINE_LENGTH ); // "ghost_flags"
3232 geoFile.skip( nbElems * sizeof(int) );
3235 if ( haveNodeIds && isGold && !geoFile.eof() ) {
3236 TStrOwner nextLine( geoFile.getLine() ); // "node_ids"
3237 if ( contains( "node_ids", nextLine ) )
3238 geoFile.skip( NumberOfNodes * sizeof(int) );
3240 geoFile.skip( -MAX_LINE_LENGTH );
3243 if ( haveElemIds && isGold && !geoFile.eof() ) {
3244 TStrOwner nextLine( geoFile.getLine() ); // "element_ids"
3245 if ( contains( "element_ids", nextLine ) )
3246 geoFile.skip( nbElems * sizeof(int) );
3248 geoFile.skip( -MAX_LINE_LENGTH );
3251 else if ( word == "extents" )
3253 geoFile.skip( 6 * sizeof(float) );
3257 throw MED_EXCEPTION ( LOCALIZED( STRING(LOC) << "Unexpected word: " << word ));
3259 } // while ( !geoFile.eof() )
3264 GROUP* ENSIGHT_MESH_RDONLY_DRIVER::makeGroup( _groupe& grp,
3267 //const char* LOC = "ENSIGHT_MESH_RDONLY_DRIVER::makeGroup(): error";
3269 // prevent creation of other groups but only this one
3270 for (size_t i=0; i < imed.groupes.size(); ++i)
3271 imed.groupes[i].nom.clear();
3273 // let _intermediateMED create a GROUP from grp
3274 grp.medGroup = 0; // the new GROUP should appear in grp.medGroup
3276 vector<GROUP *> tmp;
3277 MESH* mesh = imed._isOwnMedMesh ? (MESH*) 0 : imed._medMesh;
3278 imed.getGroups( tmp, tmp, tmp, tmp, mesh );
3279 if ( !grp.medGroup )
3280 throw MEDEXCEPTION(LOCALIZED("Can't create a GROUP from _groupe"));
3282 grp.medGroup->setName(""); // to let a caller set a proper name
3285 // find pre-existing equal _groupe
3286 _groupe * equalGroupe = 0;
3287 for (unsigned int i=0; i < imed.groupes.size() && !equalGroupe; ++i) {
3288 _groupe& g = imed.groupes[i];
3289 if ( &g != &grp && g.medGroup && g.medGroup->deepCompare( *grp.medGroup ))
3292 if ( equalGroupe ) {
3293 if ( grp.medGroup->getMesh() )
3294 grp.medGroup->getMesh()->addReference(); // let the mesh survive
3295 grp.medGroup->removeReference();
3296 grp.medGroup = equalGroupe->medGroup;
3298 else { // new unique group
3300 if ( grp.medGroup->isOnAllElements() ) // on all elems
3301 grp.medGroup->setName( string("SupportOnAll_")+entNames[ grp.medGroup->getEntity() ] );
3303 // add a new group to mesh
3304 if ( !imed._isOwnMedMesh ) {
3305 vector<GROUP*> * groups = 0;
3306 switch ( grp.medGroup->getEntity() ) {
3307 case MED_CELL: groups = & imed._medMesh->_groupCell; break;
3308 case MED_FACE: groups = & imed._medMesh->_groupFace; break;
3309 case MED_EDGE: groups = & imed._medMesh->_groupEdge; break;
3310 case MED_NODE: groups = & imed._medMesh->_groupNode; break;
3314 groups->resize( groups->size() + 1 );
3315 groups->at( groups->size() - 1) = grp.medGroup;
3319 return grp.medGroup;