1 // Copyright (C) 2007-2008 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
22 #include "MEDMEM_EnsightMeshDriver.hxx"
28 #include "MEDMEM_define.hxx"
29 #include "MEDMEM_Family.hxx"
30 #include "MEDMEM_Group.hxx"
31 #include "MEDMEM_Coordinate.hxx"
32 #include "MEDMEM_Connectivity.hxx"
33 #include "MEDMEM_Mesh.hxx"
34 #include "MEDMEM_CellModel.hxx"
35 #include "MEDMEM_Grid.hxx"
37 #include "MEDMEM_MedMeshDriver.hxx"
40 using namespace MEDMEM;
41 using namespace MED_EN;
42 using namespace MEDMEM_ENSIGHT;
44 #define TStrTool _ASCIIFileReader
45 #define TOLERANCE 1e-15
47 //#define ELEMENT_ID_GIVEN
51 // ---------------------------------------------------------------------
53 * \brief The beginning of mesh description used to distinguish files
54 * generated by ENSIGHT_MESH_WRONLY_DRIVER from others
56 const char* theDescriptionPrefix = "Meshing from MedMemory. ";
58 // ---------------------------------------------------------------------
60 * \brief Default name of a mesh read from EnSight
62 const char* theDefaultMeshName = "EnsightMesh";
64 // ---------------------------------------------------------------------
66 * \brief Max number of types in EnSight part
68 const int theMaxNbTypes = 20;
70 // ---------------------------------------------------------------------
72 * \brief Make array with elements == index[i+1]-index[i]
73 * \param size - result array size
75 int* getNumbersByIndex( const int* index, int size, const int* elemNumbers=0)
77 int* numbers = new int[size];
80 const int *elem = elemNumbers-1, *elemEnd = elemNumbers + size;
81 while ( ++elem < elemEnd )
82 *n++ = index[*elem] - index[*elem-1];
85 const int *ind = index, *indEnd = index + size + 1;
86 while ( ++ind < indEnd )
87 *n++ = ind[0] - ind[-1];
91 // ---------------------------------------------------------------------
93 * \brief Type used to delete numbers returned by getNumbersByIndex()
95 typedef _ValueOwner<int> TNumbers;
100 ENSIGHT_MESH_DRIVER::ENSIGHT_MESH_DRIVER(): _CaseFileDriver_User(), _ptrMesh((MESH *)NULL)
104 ENSIGHT_MESH_DRIVER::ENSIGHT_MESH_DRIVER(const string & fileName,
106 :_CaseFileDriver_User(fileName,MED_EN::RDWR), _ptrMesh(ptrMesh),
107 _meshName(ptrMesh->getName())
111 ENSIGHT_MESH_DRIVER::ENSIGHT_MESH_DRIVER(const string & fileName,
113 med_mode_acces accessMode)
114 :_CaseFileDriver_User(fileName,accessMode), _ptrMesh(ptrMesh), _meshName(ptrMesh->getName())
118 ENSIGHT_MESH_DRIVER::ENSIGHT_MESH_DRIVER(const ENSIGHT_MESH_DRIVER & driver)
119 :_CaseFileDriver_User(driver), _ptrMesh(driver._ptrMesh), _meshName(driver._meshName)
123 ENSIGHT_MESH_DRIVER::~ENSIGHT_MESH_DRIVER()
125 MESSAGE_MED("ENSIGHT_MESH_DRIVER::~ENSIGHT_MESH_DRIVER() has been destroyed");
128 void ENSIGHT_MESH_DRIVER::setMeshName(const string & meshName) { _meshName = meshName; };
130 string ENSIGHT_MESH_DRIVER::getMeshName() const { return _meshName; };
132 void ENSIGHT_MESH_DRIVER::openConst(bool checkDataFile) const
134 const char * LOC ="ENSIGHT_MESH_DRIVER::open() : ";
139 if ( getDataFileName().empty() )
141 ( LOCALIZED( STRING(LOC) << "Internal error, geometry file name is empty"));
143 if (!canOpenFile( getDataFileName(), getAccessMode() ))
145 ( LOCALIZED( STRING(LOC) << "Can not open Ensight Geometry file " << getDataFileName()
146 << " in access mode " << getAccessMode()));
150 if ( getCaseFileName().empty() )
152 ( LOCALIZED( STRING(LOC) << "Case file name is empty, "
153 "please set a correct fileName before calling open()"));
155 if ( !canOpenFile( getCaseFileName(), getAccessMode() ))
157 ( LOCALIZED( STRING(LOC) << "Can not open Ensight Case file " << getCaseFileName()
158 << " in access mode " << getAccessMode()));
164 void ENSIGHT_MESH_DRIVER::open() {
168 void ENSIGHT_MESH_DRIVER::close() {
171 // ================================================================================
173 // ================================================================================
175 ENSIGHT_MESH_WRONLY_DRIVER::ENSIGHT_MESH_WRONLY_DRIVER() : ENSIGHT_MESH_DRIVER()
179 ENSIGHT_MESH_WRONLY_DRIVER::ENSIGHT_MESH_WRONLY_DRIVER(const string & fileName,
182 : ENSIGHT_MESH_DRIVER( fileName, ptrMesh, WRONLY ), _append(append)
186 ENSIGHT_MESH_WRONLY_DRIVER::ENSIGHT_MESH_WRONLY_DRIVER(const ENSIGHT_MESH_WRONLY_DRIVER & driver)
187 : ENSIGHT_MESH_DRIVER(driver),_append(driver._append)
191 ENSIGHT_MESH_WRONLY_DRIVER::~ENSIGHT_MESH_WRONLY_DRIVER()
195 GENDRIVER * ENSIGHT_MESH_WRONLY_DRIVER::copy() const
197 return new ENSIGHT_MESH_WRONLY_DRIVER(*this) ;
200 void ENSIGHT_MESH_WRONLY_DRIVER::read() throw (MEDEXCEPTION) {
201 throw MEDEXCEPTION("ENSIGHT_MESH_WRONLY_DRIVER::read : Can't read with a WRONLY driver !");
204 //================================================================================
208 //================================================================================
210 void ENSIGHT_MESH_WRONLY_DRIVER::write() const throw (MEDEXCEPTION)
212 const char * LOC = "ENSIGHT_MESH_WRONLY_DRIVER::write() : ";
215 openConst(false) ; // check if can write to case file
217 // Ensight case organization requires a main file (filename.case) which defines organization
219 _CaseFileDriver caseFile( getCaseFileName(), this);
222 caseFile.addMesh( this );
223 caseFile.write(); // all fields of _CaseFileDriver_User are set by this method
225 openConst(true) ; // check if can write to data file
227 cout << "-> creating the Ensight geometry file " << getDataFileName() << endl ;
229 // Store mesh description and a special mark in the first two description lines each
230 // of 79 chars length maximum, while MED mesh description is up to 200 chars
231 const char* line1 = theDescriptionPrefix;
232 string line2 = _ptrMesh->getDescription();
233 for ( int i = 0; i < line2.size(); ++i ) { // protect from gabage
234 if ( !isascii( line2[ i ])) {
239 if ( line2.size() >= MAX_LINE_LENGTH )
240 line2.resize( MAX_LINE_LENGTH );
242 // EnSight will assign node/element visible numbers it-self
243 const char* line3 = "node id assign";
244 #ifdef ELEMENT_ID_GIVEN
245 const char* line4 = "element id given";
247 const char* line4 = "element id assign";
250 if ( isBinaryEnSightFormatForWriting() )
252 // ======================================================
254 // ======================================================
256 _BinaryFileWriter ensightGeomFile( getDataFileName() );
258 ensightGeomFile.addString("C Binary");
259 ensightGeomFile.addString(line1);
260 ensightGeomFile.addString(line2);
261 ensightGeomFile.addString(line3);
262 ensightGeomFile.addString(line4);
264 // function to write a support as a part
265 typedef void (ENSIGHT_MESH_WRONLY_DRIVER::* TWritePart) (_BinaryFileWriter&, const SUPPORT*) const;
266 TWritePart writePart;
267 if ( isGoldFormat() )
270 writePart = & ENSIGHT_MESH_WRONLY_DRIVER::writePartGoldBinary;
274 // ENSIGHT 6. Write addionally global nodes
275 writePart = & ENSIGHT_MESH_WRONLY_DRIVER::writePart6Binary;
277 // All point are in 3D, so if we are in 1D or 2D, we complete by zero !
278 int SpaceDimension = _ptrMesh->getSpaceDimension() ;
279 int NumberOfNodes = _ptrMesh->getNumberOfNodes() ;
280 ensightGeomFile.addString("coordinates");
281 ensightGeomFile.addInt( NumberOfNodes );
282 const double *coordinate = _ptrMesh->getCoordinates(MED_FULL_INTERLACE) ;
283 if ( SpaceDimension == 3 ) {
284 ensightGeomFile.addReal(coordinate, NumberOfNodes * SpaceDimension );
287 typedef _ValueIterator< double > TComponentIt;
288 vector< TComponentIt > coordCompIt( 3 );
289 for (int j=0; j<3; j++, coordinate++)
290 if ( j < SpaceDimension )
291 coordCompIt[ j ] = TComponentIt( coordinate, SpaceDimension );
293 coordCompIt[ j ] = TComponentIt(); // to iterate on zeros
294 ensightGeomFile.addReal( coordCompIt, NumberOfNodes, MED_FULL_INTERLACE );
298 // We put connectivity
300 if ( isToWriteEntity( MED_CELL, _ptrMesh ))
302 SUPPORT allCells(_ptrMesh, getMeshName(), MED_CELL );
303 (this->*writePart)( ensightGeomFile, &allCells );
305 // And meshdim-1 connectivity
306 if ( isToWriteEntity( MED_FACE, _ptrMesh ))
308 SUPPORT allFaces(_ptrMesh, string("SupportOnAll_")+entNames[MED_FACE], MED_FACE );
309 (this->*writePart)( ensightGeomFile, & allFaces);
311 else if ( isToWriteEntity(MED_EDGE, _ptrMesh))
313 SUPPORT allEdges(_ptrMesh, string("SupportOnAll_")+entNames[MED_EDGE], MED_EDGE );
314 (this->*writePart)( ensightGeomFile, & allEdges);
317 // Write all groups as parts
319 for ( int ent = MED_CELL; ent < MED_ALL_ENTITIES; ++ent )
321 medEntityMesh entity = (medEntityMesh) ent;
322 int nbGroups = _ptrMesh->getNumberOfGroups(entity);
323 for ( int i=1; i<=nbGroups; i++)
325 const GROUP* group = _ptrMesh->getGroup( entity, i );
326 (this->*writePart)( ensightGeomFile, group );
333 // ======================================================
335 // ======================================================
336 ofstream ensightGeomFile( getDataFileName().c_str(), ios::out);
337 ensightGeomFile.setf(ios::scientific);
338 ensightGeomFile.precision(5);
340 ensightGeomFile << line1 << endl
345 // function to write a support as a part
346 typedef void (ENSIGHT_MESH_WRONLY_DRIVER::* TWritePart) (ofstream&, const SUPPORT*) const;
347 TWritePart writePart;
348 if ( isGoldFormat() )
351 writePart = & ENSIGHT_MESH_WRONLY_DRIVER::writePartGoldASCII;
355 // ENSIGHT 6. Write addionally global nodes
356 writePart = & ENSIGHT_MESH_WRONLY_DRIVER::writePart6ASCII;
358 // Put points (all point are in 3D, so if we are in 1D or 2D, we complete by zero !
359 int SpaceDimension = _ptrMesh->getSpaceDimension() ;
360 int NumberOfNodes = _ptrMesh->getNumberOfNodes() ;
362 if (SpaceDimension==2) zeros = " 0.00000e+00";
363 if (SpaceDimension==1) zeros = " 0.00000e+00 0.00000e+00";
364 ensightGeomFile << "coordinates" << endl
365 << setw(8) << NumberOfNodes << endl ;
366 const double *coordinate = _ptrMesh->getCoordinates(MED_FULL_INTERLACE) ;
367 for (int i=0; i<NumberOfNodes; i++)
369 //ensightGeomFile << setw(8) << i+1 ; // node id
370 for (int j=0; j<SpaceDimension; j++, coordinate++)
371 ensightGeomFile << setw(12) << *coordinate;
372 ensightGeomFile << zeros << endl ;
376 // We put connectivity
378 if ( isToWriteEntity( MED_CELL, _ptrMesh ))
380 SUPPORT allCells(_ptrMesh, getMeshName(), MED_CELL );
381 (this->*writePart)( ensightGeomFile, &allCells );
383 // And meshdim-1 connectivity
384 if ( isToWriteEntity( MED_FACE, _ptrMesh ))
386 SUPPORT allFaces(_ptrMesh, string("SupportOnAll_")+entNames[MED_FACE], MED_FACE );
387 (this->*writePart)( ensightGeomFile, & allFaces);
389 else if ( isToWriteEntity(MED_EDGE, _ptrMesh))
391 SUPPORT allEdges(_ptrMesh, string("SupportOnAll_")+entNames[MED_EDGE], MED_EDGE );
392 (this->*writePart)( ensightGeomFile, & allEdges);
395 // Write all groups as parts
397 for ( int ent = MED_CELL; ent < MED_ALL_ENTITIES; ++ent )
399 medEntityMesh entity = (medEntityMesh) ent;
400 int nbGroups = _ptrMesh->getNumberOfGroups(entity);
401 for ( int i=1; i<=nbGroups; i++)
403 const GROUP* group = _ptrMesh->getGroup( entity, i );
404 (this->*writePart)( ensightGeomFile, group );
408 ensightGeomFile.close();
410 } // end ASCII format
412 } // ENSIGHT_MESH_WRONLY_DRIVER::write()
414 //================================================================================
416 * \brief Write support as an EnSight Gold part
418 //================================================================================
420 void ENSIGHT_MESH_WRONLY_DRIVER::writePartGoldBinary(_BinaryFileWriter& ensightGeomFile,
421 const SUPPORT* support) const
424 int partNum = getPartNumber( support );
426 throw MED_EXCEPTION ( LOCALIZED( STRING("Internal error: invalid part number")));
427 ensightGeomFile.addString( "part" );
428 ensightGeomFile.addInt( partNum );
431 ensightGeomFile.addString( support->getName() );
434 medEntityMesh entity = support->getEntity();
435 int nbTypes = support->getNumberOfTypes();
436 const medGeometryElement* geoType = support->getTypes();
438 const int * connectivity = 0;
439 const int * elemConnectivity = 0;
440 const int * index = 0;
443 // COORDINATES Gold binary
444 // ===================================================================================
445 // In Ensight, coordinates of nodes of support elements are in MED_NO_INTERLACE mode.
446 // We are to write only nodes belonging to elements of the support and
447 // nodal connectivity should refer to these nodes.
448 map<int, int> med2ensIds;
449 map<int, int>::iterator medEnsIt;
450 int SpaceDimension = _ptrMesh->getSpaceDimension() ;
451 int NumberOfNodes = _ptrMesh->getNumberOfNodes() ;
452 // -------------------------------------------------
453 if ( support->isOnAllElements() )
456 ensightGeomFile.addString( "coordinates" );
457 ensightGeomFile.addInt( NumberOfNodes );
460 const double *coordinate = _ptrMesh->getCoordinates(MED_FULL_INTERLACE);
461 typedef _ValueIterator< double > TComponentIt;
462 vector< TComponentIt > coordCompIt( 1 );
463 for (int j=0; j<SPACE_DIM; j++, coordinate++) { // loop on dimensions
464 if ( j < SpaceDimension )
465 coordCompIt[ 0 ] = TComponentIt( coordinate, SpaceDimension );
467 coordCompIt[ 0 ] = TComponentIt(); // to iterate on zeros
468 ensightGeomFile.addReal( coordCompIt, NumberOfNodes, MED_NO_INTERLACE );
471 // -------------------------------------------------
472 else // support is not on all elements
475 getSupportNodes( support, med2ensIds );
476 NumberOfNodes = med2ensIds.size();
477 ensightGeomFile.addString( "coordinates" );
478 ensightGeomFile.addInt( NumberOfNodes );
481 vector<float> floatCoords( NumberOfNodes );
482 for ( j=0; j < SPACE_DIM; j++) { // loop on dimensions
483 medEnsIt = med2ensIds.begin();
484 if ( j < SpaceDimension ) {
485 const double *coordinate = _ptrMesh->getCoordinates(MED_FULL_INTERLACE) + j;
486 for (int i=0; i<NumberOfNodes; i++, ++medEnsIt )
487 floatCoords[ i ] = (float) coordinate[ (medEnsIt->first-1) * SpaceDimension];
489 else if ( j-1 < SpaceDimension ) {
490 for (int i=0; i<NumberOfNodes; i++)
491 floatCoords[ i ] = 0.;
493 ensightGeomFile.addReal( &floatCoords[0], NumberOfNodes );
495 // assign local node ids
496 for ( medEnsIt = med2ensIds.begin(), j=1; j<=NumberOfNodes; j++, ++medEnsIt )
497 medEnsIt->second = j;
500 // CONNECTIVITY Gold binary
501 // ===================================================================================
503 for (int i=0; i<nbTypes; i++)
505 const medGeometryElement medType = geoType[i];
506 const TEnSightElemType& ensightType = getEnSightType(medType);
507 const int numberOfCell = support->getNumberOfElements(medType);
508 int nbCellNodes = ensightType._medIndex.size();
510 // type name and nb cells
511 ensightGeomFile.addString( ensightType._name );
512 ensightGeomFile.addInt( numberOfCell );
516 // -------------------------------------------------
517 if ( support->isOnAllElements() )
519 #ifdef ELEMENT_ID_GIVEN
521 nodeIds.resize( numberOfCell );
522 for ( j = 1; j <= numberOfCell; j++)
524 ensightGeomFile.addInt( nodeIds );
526 if ( entity != MED_NODE ) nodeIds.clear();
529 if ( nbCellNodes > 1 ) // STANDARD ELEMENTS connectivity
531 connectivity = _ptrMesh->getConnectivity(MED_FULL_INTERLACE, MED_NODAL,
533 nodeIds.reserve( numberOfCell * nbCellNodes);
534 for (j = 0 ; j < numberOfCell; j++, connectivity += nbCellNodes)
535 for (int k=0; k<nbCellNodes; k++)
536 nodeIds.push_back( connectivity[ ensightType._medIndex[k] ]);
537 ensightGeomFile.addInt( nodeIds );
539 else if ( entity == MED_NODE ) // NODES connectivity
541 #if !defined(ELEMENT_ID_GIVEN)
542 nodeIds.resize( numberOfCell );
543 for ( j = 1; j <= numberOfCell; j++)
546 ensightGeomFile.addInt( nodeIds );
548 else if ( medType == MED_POLYGON ) // POLYGONs connectivity
550 connectivity = _ptrMesh->getPolygonsConnectivity(MED_NODAL, entity);
551 index = _ptrMesh->getPolygonsConnectivityIndex(MED_NODAL, entity);
552 int connLength = _ptrMesh->getPolygonsConnectivityLength(MED_NODAL, entity);
553 // number of nodes in each element
555 TIntOwner nbNodesInPoly( getNumbersByIndex( index, numberOfCell ));
556 ensightGeomFile.addInt( nbNodesInPoly.myValues, numberOfCell );
557 } // nbNodesInPoly is deleted here
560 ensightGeomFile.addInt( connectivity, connLength );
562 else // POLYHEDRA connectivity
564 connectivity = _ptrMesh->getPolyhedronConnectivity(MED_NODAL);
565 index = _ptrMesh->getPolyhedronIndex(MED_NODAL);
566 const int * fIndex = _ptrMesh->getPolyhedronFacesIndex();
567 int connLength = _ptrMesh->getPolyhedronConnectivityLength(MED_NODAL);
568 int nbFaces = _ptrMesh->getNumberOfPolyhedronFaces();
569 // nb of faces in each polyhedron
571 TIntOwner nbFacesInPoly( getNumbersByIndex( index, numberOfCell ));
572 ensightGeomFile.addInt( nbFacesInPoly.myValues, numberOfCell );
574 // number of nodes in each face
576 TIntOwner nbNodesInFace( getNumbersByIndex( fIndex, nbFaces ));
577 ensightGeomFile.addInt( nbNodesInFace.myValues, nbFaces );
580 ensightGeomFile.addInt( connectivity, connLength );
583 // -------------------------------------------------
584 else // support is not on all elements
586 const int *number = support->getNumber(medType);
588 #ifdef ELEMENT_ID_GIVEN
589 ensightGeomFile.addInt( number, numberOfCell );
591 if ( nbCellNodes > 1 ) // STANDARD ELEMENTS connectivity
593 connectivity = _ptrMesh->getConnectivity(MED_FULL_INTERLACE, MED_NODAL,
594 entity, MED_ALL_ELEMENTS);
595 index = _ptrMesh->getConnectivityIndex(MED_FULL_INTERLACE, entity);
597 nodeIds.reserve( numberOfCell * nbCellNodes);
598 for (j=0; j<numberOfCell; j++) {
599 int elem = number[j];
600 elemConnectivity = connectivity + index[elem-1]-1;
601 for (int k=0; k<nbCellNodes; k++)
603 int node = elemConnectivity[ ensightType._medIndex[k] ];
604 nodeIds.push_back( med2ensIds[ node ]);
607 ensightGeomFile.addInt( nodeIds );
609 else if ( entity == MED_NODE ) // NODES connectivity
611 nodeIds.resize( numberOfCell );
612 for ( j = 1; j <= numberOfCell; j++)
614 ensightGeomFile.addInt( nodeIds );
616 else if ( medType == MED_POLYGON ) // POLYGONs connectivity
618 connectivity = _ptrMesh->getPolygonsConnectivity(MED_NODAL, entity);
619 index = _ptrMesh->getPolygonsConnectivityIndex(MED_NODAL, entity);
620 int connLength = _ptrMesh->getPolygonsConnectivityLength(MED_NODAL, entity);
621 int nbStdElems = _ptrMesh->getNumberOfElements(entity,MED_ALL_ELEMENTS);
622 // number of nodes in each element
624 TIntOwner nbNodesInPoly( getNumbersByIndex( index-nbStdElems, numberOfCell, number ));
625 ensightGeomFile.addInt( nbNodesInPoly.myValues, numberOfCell );
626 } // nbNodesInPoly is deleted here
629 nodeIds.reserve( connLength );
630 for ( j = 0; j < numberOfCell; ++j )
632 int elem = number[ j ] - nbStdElems;
633 elemConnectivity = connectivity + index[ elem-1 ]-1;
634 const int* connEnd = connectivity + index[ elem ]-1;
635 while ( elemConnectivity < connEnd )
636 nodeIds.push_back( med2ensIds[ *elemConnectivity++ ]);
638 ensightGeomFile.addInt( nodeIds );
640 else // POLYHEDRA connectivity
642 connectivity = _ptrMesh->getPolyhedronConnectivity(MED_NODAL);
643 index = _ptrMesh->getPolyhedronIndex(MED_NODAL);
644 const int * fIndex = _ptrMesh->getPolyhedronFacesIndex();
645 int connLength = _ptrMesh->getPolyhedronConnectivityLength(MED_NODAL);
646 int nbStdElems = _ptrMesh->getNumberOfElements(entity,MED_ALL_ELEMENTS);
647 // nb of faces in each polyhedron
649 TIntOwner nbFacesInPoly( getNumbersByIndex( index-nbStdElems, numberOfCell, number ));
650 ensightGeomFile.addInt( nbFacesInPoly.myValues, numberOfCell );
652 // number of nodes in each face
653 nodeIds.reserve( connLength );
654 for ( j = 0; j < numberOfCell; ++j )
656 int elem = number[ j ] - nbStdElems;
657 int f1 = index[ elem-1 ] - 1, f2 = index[ elem ] - 1;
658 int nbFaces = f2 - f1;
660 TIntOwner nbNodesInFace( getNumbersByIndex( &fIndex[ f1 ], nbFaces ));
661 ensightGeomFile.addInt( nbNodesInFace.myValues, nbFaces );
663 elemConnectivity = connectivity + fIndex[ f1 ] - 1;
664 const int* connEnd = connectivity + fIndex[ f2 ] - 1;
665 while ( elemConnectivity < connEnd )
666 nodeIds.push_back( med2ensIds[ *elemConnectivity++ ]);
669 ensightGeomFile.addInt( nodeIds );
673 } // writePartGoldBinary()
675 //================================================================================
677 * \brief Write support as an EnSight Gold part
679 //================================================================================
681 void ENSIGHT_MESH_WRONLY_DRIVER::writePartGoldASCII(ofstream& ensightGeomFile,
682 const SUPPORT* support) const
687 int partNum = getPartNumber( support );
688 ensightGeomFile << "part" << endl
689 << setw(iw) << partNum << endl;
691 throw MED_EXCEPTION ( LOCALIZED( STRING("Internal error: invalid part number")));
694 ensightGeomFile << support->getName() << endl;
697 medEntityMesh entity = support->getEntity();
698 int nbTypes = support->getNumberOfTypes();
699 const medGeometryElement* geoType = support->getTypes();
701 const int * connectivity = 0;
702 const int * elemConnectivity = 0;
703 const int * index = 0;
706 // COORDINATES Gold ASCII
707 // ===================================================================================
708 // In Ensight, coordinates of nodes of support elements are in MED_NO_INTERLACE mode.
709 // We are to write only nodes belonging to elements of the support and
710 // nodal connectivity should refer to these nodes.
711 map<int, int> med2ensIds;
712 map<int, int>::iterator medEnsIt;
713 int SpaceDimension = _ptrMesh->getSpaceDimension() ;
714 int NumberOfNodes = _ptrMesh->getNumberOfNodes() ;
715 string zeroStr = " 0.00000e+00";
716 // -----------------------------------
717 if ( support->isOnAllElements() )
720 ensightGeomFile << "coordinates" << endl
721 << setw(iw) << NumberOfNodes << endl ;
724 for (j=0; j<SPACE_DIM; j++) { // loop on dimensions
725 if ( j < SpaceDimension ) {
726 const double *coordinate = _ptrMesh->getCoordinates(MED_FULL_INTERLACE) + j;
727 for (int i=0; i<NumberOfNodes; i++, coordinate += SpaceDimension)
728 ensightGeomFile << setw(12) << (float) *coordinate << endl;
731 for (int i=0; i<NumberOfNodes; i++)
732 ensightGeomFile << zeroStr << endl;
736 // -----------------------------------
737 else // support is not on all elements
740 getSupportNodes( support, med2ensIds );
741 NumberOfNodes = med2ensIds.size();
742 ensightGeomFile << "coordinates" << endl
743 << setw(iw) << NumberOfNodes << endl ;
746 for ( j=0; j<SPACE_DIM; j++) { // loop on dimensions
747 medEnsIt = med2ensIds.begin();
748 if ( j < SpaceDimension ) {
749 const double *coordinate = _ptrMesh->getCoordinates(MED_FULL_INTERLACE) + j;
750 for (int i=0; i<NumberOfNodes; i++, ++medEnsIt )
751 ensightGeomFile << setw(12)
752 << (float) coordinate[ (medEnsIt->first-1) * SpaceDimension] << endl;
755 for (int i=0; i<NumberOfNodes; i++)
756 ensightGeomFile << zeroStr << endl;
759 // assign local node ids
760 for ( medEnsIt = med2ensIds.begin(), j=1; j<=NumberOfNodes; j++, ++medEnsIt )
761 medEnsIt->second = j;
764 // CONNECTIVITY Gold ASCII
765 // ===================================================================================
767 for (int i=0; i<nbTypes; i++)
769 const medGeometryElement medType = geoType[i];
770 const TEnSightElemType& ensightType = getEnSightType(medType);
771 const int numberOfCell = support->getNumberOfElements(medType);
772 int nbCellNodes = ensightType._medIndex.size();
774 // type name and nb cells
775 ensightGeomFile << ensightType._name << endl
776 << setw(iw) << numberOfCell << endl;
778 // -----------------------------------
779 if ( support->isOnAllElements() )
781 #ifdef ELEMENT_ID_GIVEN
782 for ( j = 1; j <= numberOfCell; j++)
783 ensightGeomFile << setw(iw) << j << endl;
786 if ( nbCellNodes > 1 ) // STANDARD ELEMENTS connectivity
788 connectivity = _ptrMesh->getConnectivity(MED_FULL_INTERLACE, MED_NODAL,
790 for (j = 0 ; j < numberOfCell; j++, connectivity += nbCellNodes) {
791 for (int k=0; k<nbCellNodes; k++)
792 ensightGeomFile << setw(iw) << connectivity[ ensightType._medIndex[k] ];
793 ensightGeomFile << endl ;
796 else if ( entity == MED_NODE ) // NODES connectivity
798 for ( j = 1; j <= numberOfCell; j++)
799 ensightGeomFile << setw(iw) << j << endl;
801 else if ( medType == MED_POLYGON ) // POLYGONs connectivity
803 connectivity = _ptrMesh->getPolygonsConnectivity(MED_NODAL, entity);
804 index = _ptrMesh->getPolygonsConnectivityIndex(MED_NODAL, entity);
805 // number of nodes in each element
806 const int* ind = index;
807 for (j = 0 ; j < numberOfCell; j++, ++ind)
808 ensightGeomFile << setw(iw) << ( ind[1] - ind[0] ) << endl;
811 for (j = 0; j < numberOfCell; j++, ++index) {
812 nbCellNodes = index[1] - index[0];
813 for (int k=0; k<nbCellNodes; k++, ++connectivity)
814 ensightGeomFile << setw(iw) << *connectivity;
815 ensightGeomFile << endl;
818 else // POLYHEDRA connectivity
820 connectivity = _ptrMesh->getPolyhedronConnectivity(MED_NODAL);
821 index = _ptrMesh->getPolyhedronIndex(MED_NODAL);
822 const int * fIndex = _ptrMesh->getPolyhedronFacesIndex();
823 int nbFaces = _ptrMesh->getNumberOfPolyhedronFaces();
824 // nb of faces in each polyhedron
825 const int* ind = index;
826 for (j = 0 ; j < numberOfCell; j++, ++ind)
827 ensightGeomFile << setw(iw) << ( ind[1] - ind[0] ) << endl ;
828 // number of nodes in each face
830 for (j = 0 ; j < nbFaces; j++, ++ind)
831 ensightGeomFile << setw(iw) << ( ind[1] - ind[0] ) << endl ;
832 // connectivity of each face
833 for (j = 0 ; j < nbFaces; j++, ++fIndex) {
834 int nbFaceNodes = fIndex[1] - fIndex[0];
835 for (int k=0; k<nbFaceNodes; k++, ++connectivity)
836 ensightGeomFile << setw(iw) << *connectivity;
837 ensightGeomFile << endl ;
841 // -----------------------------------
842 else // support is not on all elements
844 const int *number = support->getNumber(medType);
846 #ifdef ELEMENT_ID_GIVEN
847 for ( j = 0; j < numberOfCell; j++)
848 ensightGeomFile << setw(iw) << number[j] << endl;
850 if ( nbCellNodes > 1 ) // STANDARD ELEMENTS connectivity
852 connectivity = _ptrMesh->getConnectivity(MED_FULL_INTERLACE, MED_NODAL,
853 entity, MED_ALL_ELEMENTS);
854 index = _ptrMesh->getConnectivityIndex(MED_FULL_INTERLACE, entity);
856 for (j=0; j<numberOfCell; j++) {
857 int elem = number[j];
858 elemConnectivity = connectivity + index[elem-1]-1;
859 for (int k=0; k<nbCellNodes; k++) {
860 int node = elemConnectivity[ ensightType._medIndex[k] ];
861 ensightGeomFile << setw(iw) << med2ensIds[ node ];
863 ensightGeomFile << endl;
866 else if ( entity == MED_NODE ) // NODES connectivity
868 for (j=0; j<numberOfCell; j++) {
869 int node = med2ensIds[ number[j] ];
870 ensightGeomFile << setw(iw) << node << endl ;
873 else if ( medType == MED_POLYGON ) // POLYGONs connectivity
875 connectivity = _ptrMesh->getPolygonsConnectivity(MED_NODAL, entity);
876 index = _ptrMesh->getPolygonsConnectivityIndex(MED_NODAL, entity);
877 int nbStdElems = _ptrMesh->getNumberOfElements(entity,MED_ALL_ELEMENTS);
878 // number of nodes in each element
879 for (j = 0 ; j < numberOfCell; j++) {
880 int elem = number[j] - nbStdElems;
881 ensightGeomFile << setw(iw) << ( index[elem] - index[elem-1] ) << endl;
884 for ( j = 0; j < numberOfCell; ++j ) {
885 int elem = number[ j ] - nbStdElems;
886 elemConnectivity = connectivity + index[ elem-1 ]-1;
887 const int* connEnd = connectivity + index[ elem ]-1;
888 while ( elemConnectivity < connEnd )
889 ensightGeomFile << setw(iw) << med2ensIds[ *elemConnectivity++ ];
890 ensightGeomFile << endl;
893 else // POLYHEDRA connectivity
895 connectivity = _ptrMesh->getPolyhedronConnectivity(MED_NODAL);
896 index = _ptrMesh->getPolyhedronIndex(MED_NODAL);
897 const int * fIndex = _ptrMesh->getPolyhedronFacesIndex();
898 int nbStdElems = _ptrMesh->getNumberOfElements(entity,MED_ALL_ELEMENTS);
899 // nb of faces in each polyhedron
900 for (j = 0 ; j < numberOfCell; j++) {
901 int elem = number[j] - nbStdElems;
902 ensightGeomFile << setw(iw) << ( index[elem] - index[elem-1] ) << endl;
904 // number of nodes in each face
905 for ( j = 0; j < numberOfCell; ++j ) {
906 int elem = number[ j ] - nbStdElems;
907 int f1 = index[ elem-1 ], f2 = index[ elem ];
909 ensightGeomFile << setw(iw) << ( fIndex[f1] - fIndex[f1-1] ) << endl;
913 // connectivity of each face
914 for ( j = 0; j < numberOfCell; ++j ) {
915 int elem = number[ j ] - nbStdElems;
916 int f1 = index[ elem-1 ] - 1, f2 = index[ elem ] - 1;
918 int n1 = fIndex[f1]-1, n2 = fIndex[f1+1]-1;
921 ensightGeomFile << setw(iw) << connectivity[ n1++ ];
922 ensightGeomFile << endl ;
928 } // writePartGoldASCII()
930 //================================================================================
932 * \brief Write support as an Ensight6 part
934 //================================================================================
936 void ENSIGHT_MESH_WRONLY_DRIVER::writePart6Binary(_BinaryFileWriter& ensightGeomFile,
937 const SUPPORT* support) const
940 int partNum = getPartNumber( support );
941 ensightGeomFile.addString( STRING("part ") << partNum );
943 throw MED_EXCEPTION ( LOCALIZED( STRING("Internal error: invalid part number")));
946 ensightGeomFile.addString( support->getName() );
949 medEntityMesh entity = support->getEntity();
950 int nbTypes = support->getNumberOfTypes();
951 const medGeometryElement* geoType = support->getTypes();
954 const int * connectivity = 0;
955 if ( entity != MED_NODE )
956 connectivity = _ptrMesh->getConnectivity(MED_FULL_INTERLACE, MED_NODAL,
957 entity, MED_ALL_ELEMENTS);
958 const int * elemConnectivity = connectivity;
960 // CONNECTIVITY Ensight 6 binary
961 // ===================================================================================
963 for (int i=0; i<nbTypes; i++)
965 const medGeometryElement medType = geoType[i];
966 const TEnSightElemType& ensightType = getEnSightType(medType);
967 int nbCellNodes = ensightType._medIndex.size();
968 if ( nbCellNodes == 0 )
971 // type name and nb cells
972 int numberOfCell = support->getNumberOfElements(medType);
973 ensightGeomFile.addString( ensightType._name );
974 ensightGeomFile.addInt( numberOfCell );
977 // -------------------------------------------------
978 if ( support->isOnAllElements() )
980 #ifdef ELEMENT_ID_GIVEN
981 nodeIds.resize( numberOfCell );
982 for ( j = 1; j <= numberOfCell; j++)
984 ensightGeomFile.addInt( nodeIds );
986 if ( entity == MED_NODE ) {
987 #if !defined(ELEMENT_ID_GIVEN)
988 nodeIds.resize( numberOfCell * nbCellNodes);
989 for ( j = 1; j <= numberOfCell; j++)
995 nodeIds.reserve( numberOfCell * nbCellNodes );
996 for (j = 0 ; j < numberOfCell; j++, elemConnectivity += nbCellNodes)
997 for (int k=0; k<nbCellNodes; k++)
998 nodeIds.push_back( elemConnectivity[ ensightType._medIndex[k] ]);
1000 ensightGeomFile.addInt( nodeIds );
1002 // -------------------------------------------------
1003 else // support is not on all elements
1005 const int *number = support->getNumber(medType);
1007 #ifdef ELEMENT_ID_GIVEN
1008 ensightGeomFile.addInt( number, numberOfCell );
1010 if ( entity == MED_NODE ) {
1011 ensightGeomFile.addInt( number, numberOfCell );
1014 const int* index = _ptrMesh->getConnectivityIndex(MED_FULL_INTERLACE, entity);
1016 nodeIds.reserve( numberOfCell * nbCellNodes);
1017 for (j=0; j<numberOfCell; j++) {
1018 int elem = number[j];
1019 elemConnectivity = connectivity + index[elem-1]-1;
1020 for (int k=0; k<nbCellNodes; k++)
1021 nodeIds.push_back( elemConnectivity[ ensightType._medIndex[k] ]);
1023 ensightGeomFile.addInt( nodeIds );
1028 } // writePart6Binary()
1030 //================================================================================
1032 * \brief Write support as an Ensight6 part
1034 //================================================================================
1036 void ENSIGHT_MESH_WRONLY_DRIVER::writePart6ASCII(ofstream& ensightGeomFile,
1037 const SUPPORT* support) const
1042 int partNum = getPartNumber( support );
1043 ensightGeomFile << "part " << partNum << endl;
1045 throw MED_EXCEPTION ( LOCALIZED( STRING("Internal error: invalid part number")));
1048 ensightGeomFile << support->getName() << endl;
1051 medEntityMesh entity = support->getEntity();
1052 int nbTypes = support->getNumberOfTypes();
1053 const medGeometryElement* geoType = support->getTypes();
1056 const int * connectivity = 0;
1057 if ( entity != MED_NODE )
1058 connectivity = _ptrMesh->getConnectivity(MED_FULL_INTERLACE, MED_NODAL,
1059 entity, MED_ALL_ELEMENTS);
1060 const int * elemConnectivity = connectivity;
1062 // CONNECTIVITY Ensight 6 ASCII
1063 // ===================================================================================
1065 for (int i=0; i<nbTypes; i++)
1067 const medGeometryElement medType = geoType[i];
1068 const TEnSightElemType& ensightType = getEnSightType(medType);
1069 int nbCellNodes = ensightType._medIndex.size();
1070 if ( nbCellNodes == 0 )
1073 // type name and nb cells
1074 int numberOfCell = support->getNumberOfElements(medType);
1075 ensightGeomFile << ensightType._name << endl
1076 << setw(iw) << numberOfCell << endl;
1078 // -------------------------------------------------
1079 if ( support->isOnAllElements() )
1081 if ( entity == MED_NODE ) {
1082 for ( j = 1; j <= numberOfCell; j++) {
1083 #ifdef ELEMENT_ID_GIVEN
1084 ensightGeomFile << setw(iw) << j;
1086 ensightGeomFile << setw(iw) << j << endl;
1090 for (j = 1 ; j <= numberOfCell; j++, elemConnectivity += nbCellNodes) {
1091 #ifdef ELEMENT_ID_GIVEN
1092 ensightGeomFile << setw(iw) << elem++;
1094 for (int k=0; k<nbCellNodes; k++)
1096 ensightGeomFile << setw(iw) << elemConnectivity[ ensightType._medIndex[k] ];
1098 ensightGeomFile << endl ;
1102 // -------------------------------------------------
1103 else // support is not on all elements
1105 const int *number = support->getNumber(medType);
1106 if ( entity == MED_NODE ) {
1107 for (j=0; j<numberOfCell; j++) {
1108 int node = number[j];
1109 #ifdef ELEMENT_ID_GIVEN
1110 ensightGeomFile << setw(iw) << node;
1112 ensightGeomFile << setw(iw) << node << endl ;
1116 const int* index = _ptrMesh->getConnectivityIndex(MED_FULL_INTERLACE, entity);
1118 for (j=0; j<numberOfCell; j++) {
1119 int elem = number[j];
1120 #ifdef ELEMENT_ID_GIVEN
1121 ensightGeomFile << setw(iw) << elem;
1123 elemConnectivity = connectivity + index[elem-1]-1;
1124 for (int k=0; k<nbCellNodes; k++)
1126 ensightGeomFile << setw(iw) << elemConnectivity[ ensightType._medIndex[k] ];
1128 ensightGeomFile << endl ;
1134 } // writePart6ASCII()
1136 //================================================================================
1138 * \brief Return nb of part to write
1140 //================================================================================
1142 int ENSIGHT_MESH_WRONLY_DRIVER::nbPartsToWrite() const
1145 nbParts += (int) isToWriteEntity( MED_CELL, _ptrMesh );
1146 nbParts += (int) isToWriteEntity( MED_FACE, _ptrMesh );
1147 nbParts += (int) isToWriteEntity( MED_EDGE, _ptrMesh );
1150 for ( int ent = MED_CELL; ent < MED_ALL_ENTITIES; ++ent ) {
1151 int nbGroups = _ptrMesh->getNumberOfGroups(medEntityMesh(ent));
1152 nbParts += nbGroups;
1157 // ================================================================================
1159 // ================================================================================
1161 ENSIGHT_MESH_RDONLY_DRIVER::ENSIGHT_MESH_RDONLY_DRIVER()
1162 : ENSIGHT_MESH_DRIVER(), _indexInCaseFile(1)
1166 ENSIGHT_MESH_RDONLY_DRIVER::ENSIGHT_MESH_RDONLY_DRIVER(const string & fileName,
1169 : ENSIGHT_MESH_DRIVER(fileName,ptrMesh,RDONLY), _indexInCaseFile( index )
1173 ENSIGHT_MESH_RDONLY_DRIVER::ENSIGHT_MESH_RDONLY_DRIVER(const ENSIGHT_MESH_RDONLY_DRIVER & driver) : ENSIGHT_MESH_DRIVER(driver), _indexInCaseFile( driver._indexInCaseFile )
1177 ENSIGHT_MESH_RDONLY_DRIVER::~ENSIGHT_MESH_RDONLY_DRIVER()
1181 GENDRIVER * ENSIGHT_MESH_RDONLY_DRIVER::copy() const
1183 return new ENSIGHT_MESH_RDONLY_DRIVER(*this) ;
1186 void ENSIGHT_MESH_RDONLY_DRIVER::write() const throw (MEDEXCEPTION)
1188 throw MEDEXCEPTION("ENSIGHT_MESH_RDONLY_DRIVER::write : Can't write with a RDONLY driver !");
1191 void ENSIGHT_MESH_RDONLY_DRIVER::merge ( const GENDRIVER& driver )
1193 _CaseFileDriver_User::merge( driver );
1195 const ENSIGHT_MESH_RDONLY_DRIVER* other =
1196 dynamic_cast< const ENSIGHT_MESH_RDONLY_DRIVER* >( &driver );
1198 if ( _indexInCaseFile < other->_indexInCaseFile )
1199 _indexInCaseFile = other->_indexInCaseFile;
1203 //================================================================================
1205 * \brief Read mesh in all supported formats
1207 //================================================================================
1209 void ENSIGHT_MESH_RDONLY_DRIVER::read() throw (MEDEXCEPTION)
1211 const char * LOC = "ENSIGHT_MESH_RDONLY_DRIVER::read() : " ;
1214 openConst(false); // check if can read case file
1216 _CaseFileDriver caseFile( getCaseFileName(), this);
1218 caseFile.setDataFileName( _indexInCaseFile, this ); // data from Case File is passed here
1220 openConst(true); // check if can read data file
1222 cout << "-> Entering into the geometry file " << getDataFileName() << endl ;
1224 _InterMed* imed = new _InterMed();
1225 imed->_medMesh = getMesh();
1226 imed->_isOwnMedMesh = false;
1227 imed->_needSubParts = ( caseFile.getNbVariables() > 0 );
1228 imed->groupes.reserve(1000);
1230 // to let field drivers know eventual indices of values
1231 if ( imed->_needSubParts )
1232 setInterData( imed );
1234 if ( isBinaryDataFile( getDataFileName() )) // binary
1236 if ( isGoldFormat() ) // Gold
1238 readGoldBinary( *imed );
1242 read6Binary( *imed );
1247 if ( isGoldFormat() ) // Gold
1249 readGoldASCII( *imed );
1253 read6ASCII( *imed );
1257 if ( _isMadeByMed && !imed->groupes.empty() ) {
1258 _ptrMesh->_name = imed->groupes[0].nom;
1259 imed->groupes[0].nom = "SupportOnAll_";
1260 imed->groupes[0].nom += entNames[MED_CELL];
1263 _ptrMesh->_name = theDefaultMeshName;
1265 _ptrMesh->_spaceDimension = SPACE_DIM;
1266 _ptrMesh->_meshDimension = _ptrMesh->_spaceDimension;
1267 _ptrMesh->_numberOfNodes = imed->points.size() - imed->nbMerged( MED_POINT1 );
1268 _ptrMesh->_isAGrid = 0;
1269 _ptrMesh->_coordinate = imed->getCoordinate();
1271 //Construction des groupes
1272 imed->getGroups(_ptrMesh->_groupCell,
1273 _ptrMesh->_groupFace,
1274 _ptrMesh->_groupEdge,
1275 _ptrMesh->_groupNode, _ptrMesh);
1277 _ptrMesh->_connectivity = imed->getConnectivity();
1279 _ptrMesh->createFamilies();
1281 // add attributes to families
1282 set<string> famNames;
1283 for (medEntityMesh entity=MED_CELL; entity<MED_ALL_ENTITIES; ++entity)
1285 int i, nb = _ptrMesh->getNumberOfFamilies(entity);
1286 for ( i = 1; i <= nb; ++i ) {
1287 FAMILY* f = const_cast<FAMILY*>( _ptrMesh->getFamily( entity, i ));
1288 f->setNumberOfAttributes( 1 );
1289 int* attIDs = new int[1];
1291 f->setAttributesIdentifiers( attIDs );
1292 int* attVals = new int[1];
1294 f->setAttributesValues( attVals );
1295 string* attDescr = new string[1];
1296 attDescr[0] = "med_family";
1297 f->setAttributesDescriptions( attDescr );
1298 if ( f->getName().length() > 31 ) // limit a name length
1299 f->setName( STRING("FAM_") << f->getIdentifier());
1300 // check if family is on the whole mesh entity
1301 if (_ptrMesh->getNumberOfElements( entity, MED_ALL_ELEMENTS ) ==
1302 f->getNumberOfElements( MED_ALL_ELEMENTS ))
1305 *(f->getnumber()) = MEDSKYLINEARRAY();
1307 // setAll() for groups
1308 nb = _ptrMesh->getNumberOfGroups(entity);
1309 for ( i = 1; i <= nb; ++i ) {
1310 GROUP * g = const_cast<GROUP*>( _ptrMesh->getGroup( entity, i ));
1311 if (_ptrMesh->getNumberOfElements( entity, MED_ALL_ELEMENTS ) ==
1312 g->getNumberOfElements( MED_ALL_ELEMENTS ))
1315 *(g->getnumber()) = MEDSKYLINEARRAY();
1324 //================================================================================
1326 * \brief Read mesh in Gold ASCII format
1328 //================================================================================
1330 void ENSIGHT_MESH_RDONLY_DRIVER::readGoldASCII(_InterMed & imed)
1332 const char * LOC ="ENSIGHT_MESH_RDONLY_DRIVER::readGoldASCII() : ";
1335 _ASCIIFileReader geoFile( getDataFileName() );
1337 if ( isSingleFileMode() ) {
1338 int curTimeStep = 1;
1339 while ( curTimeStep++ < getIndexInDataFile() ) {
1340 while ( !geoFile.isTimeStepEnd())
1343 while ( !geoFile.isTimeStepBeginning() )
1346 // ----------------------
1347 // Read mesh description
1348 // ----------------------
1350 string descriptionLine1 = geoFile.getLine();
1351 string descriptionLine2 = geoFile.getLine();
1353 // find out if the file was created by MED driver
1354 int prefixSize = strlen( theDescriptionPrefix );
1355 _isMadeByMed = ( descriptionLine1.substr(0, prefixSize ) == theDescriptionPrefix );
1358 descriptionLine1 = descriptionLine1.substr( prefixSize );
1359 _ptrMesh->setDescription( descriptionLine1 + descriptionLine2 );
1362 // ----------------------------------------
1363 // Find out presence of node/elem numbers
1364 // ----------------------------------------
1366 // EnSight User Manual (for v8) says:
1367 // You do not have to assign node IDs. If you do, the element connectivities are
1368 // based on the node numbers. If you let EnSight assign the node IDs, the nodes
1369 // are considered to be sequential starting at node 1, and element connectivity is
1370 // done accordingly. If node IDs are set to off, they are numbered internally;
1371 // however, you will not be able to display or query on them. If you have node
1372 // IDs in your data, you can have EnSight ignore them by specifying "node id
1373 // ignore." Using this option may reduce some of the memory taken up by the
1374 // Client and Server, but display and query on the nodes will not be available.
1376 // read "node|element id <off|given|assign|ignore>"
1377 geoFile.getWord(); geoFile.getWord();
1378 string nodeIds = geoFile.getWord();
1379 geoFile.getWord(); geoFile.getWord();
1380 string elemIds = geoFile.getWord();
1382 bool haveNodeIds = ( nodeIds == "given" || nodeIds == "ignore" );
1383 bool haveElemIds = ( elemIds == "given" || elemIds == "ignore" );
1385 // extents: xmin xmax ymin ymax zmin zmax
1386 vector<double> extents;
1387 geoFile.toNextLine();
1388 if ( strncmp( "extents", geoFile.getCurrentPtr(), 7 ) == 0 ) {
1389 geoFile.skip( /*width =*/ 7, /*nbLines =*/ 1 );
1390 extents.reserve( 6 );
1391 while ( extents.size() < extents.capacity() )
1392 extents.push_back( geoFile.getReal() );
1395 typedef map<int,_noeud>::iterator INoeud;
1396 map<int,_noeud> & points = imed.points;
1399 _groupe * partGroupe = 0;
1400 int partNum = 0, nbParts = 0;
1402 while ( !geoFile.isTimeStepEnd() )
1404 string word, restLine, line = geoFile.getLine();
1405 TStrTool::split( line, word, restLine );
1407 const TEnSightElemType & partType = getEnSightType( word );
1408 if ( partType._medType != MED_ALL_ELEMENTS )
1410 // Unstructured element type encounters
1411 // --------------------------------------
1412 int nbElemNodes = partType._medType % 100;
1413 int nbElems = geoFile.getInt(); // ne
1414 bool isGhost = isGhostType( word );
1415 int nodeShift = points.empty() ? 0 : points.rbegin()->first;
1418 vector<int> elemIds;
1419 if ( haveElemIds ) {
1420 elemIds.reserve( nbElems );
1421 while ( elemIds.size() < nbElems )
1422 elemIds.push_back( geoFile.getInt() ); // id_e
1424 if ( isGhost ) { // do not store ghost elements (?)
1425 int nbNodes = nbElems * nbElemNodes;
1426 if ( partType._name == "nsided" ) // polygons
1428 for ( int i = 0; i < nbElems; ++i )
1429 nbNodes += geoFile.getInt();
1430 geoFile.skip( nbNodes * INT_WIDTH_GOLD, /*nbLines = */nbElems );
1432 else if ( partType._name == "nfaced" ) // polyhedrons
1435 for ( int i = 0; i < nbElems; ++i )
1436 nbFaces += geoFile.getInt();
1437 for ( int f = 0; f < nbFaces; ++f )
1438 nbNodes += geoFile.getInt();
1439 geoFile.skip( nbNodes * INT_WIDTH_GOLD, /*nbLines = */nbFaces );
1441 else // standard types
1443 geoFile.skip( nbNodes, nbElemNodes, INT_WIDTH_GOLD );
1448 // add a group corresponding to subPart (geoType)
1449 imed.groupes.push_back(_groupe());
1450 _groupe & groupe = imed.groupes.back();
1451 groupe.mailles.resize( nbElems );
1453 // find out if "coordinates" has already been encountered
1454 _SubPartDesc coordDesc( partNum , "coordinates");
1455 map< _SubPartDesc, _SubPart >::iterator descPart =
1456 imed._subPartDescribed.find( coordDesc );
1457 bool haveCoords = ( descPart != imed._subPartDescribed.end() );
1459 firstNode = descPart->second.myFirstNode;
1460 nodeShift -= descPart->second.myNbNodes;
1463 // read poly element data
1464 bool isPoly = ( !nbElemNodes );
1465 vector<int> nbElemNodesVec( 1, nbElemNodes);
1466 vector<int> nbElemFaces, nbFaceNodes;
1467 if ( partType._name == "nsided" ) // polygons
1469 nbElemNodesVec.resize( nbElems );
1470 for ( int i = 0; i < nbElems; ++i )
1471 nbElemNodesVec[ i ] = geoFile.getInt(); // npi
1473 else if ( partType._name == "nfaced" ) // polyhedrons
1475 nbElemFaces.resize( nbElems );
1476 nbElemNodesVec.resize( nbElems );
1477 int totalNbFaces = 0;
1478 for ( int i = 0; i < nbElems; ++i )
1479 totalNbFaces += ( nbElemFaces[ i ] = geoFile.getInt() ); // nf_ei
1481 nbFaceNodes.resize( totalNbFaces );
1482 vector<int>::iterator nbFN = nbFaceNodes.begin();
1483 for ( int i = 0; i < nbElems; ++i ) {
1484 nbElemNodesVec[ i ] = 0;
1485 for ( int nbFaces = nbElemFaces[ i ]; nbFaces; --nbFaces, ++nbFN )
1486 nbElemNodesVec[ i ] += ( *nbFN = geoFile.getInt() ); // np(f_ei)
1489 // iterator returning nbElemNodes for standard elems and
1490 // next value from nbElemNodesVec for poly elements
1491 _ValueIterator<int> nbElemNodesIt( & nbElemNodesVec[0], isPoly ? 1 : 0);
1493 // iterator returning values form partType._medIndex for standard elems
1494 // and node index (n) for poly elements
1496 _ValueIterator<int> medIndexIt( isPoly ? &n : &partType._medIndex[0],
1498 // read connectivity
1499 _maille ma( partType._medType, nbElemNodes );
1500 ma.sommets.resize( nbElemNodes );
1502 for ( int i = 0; i < nbElems; ++i ) {
1503 _ValueIterator<int> medIndex = medIndexIt;
1504 nbElemNodes = nbElemNodesIt.next();
1505 if ( ma.sommets.size() != nbElemNodes )
1506 ma.sommets.resize( nbElemNodes );
1507 for ( n = 0; n < nbElemNodes; ++n ) {
1508 int nodeID = geoFile.getInt(); // nn_ei
1510 node = points.find( nodeID + nodeShift );
1512 node = points.insert( make_pair( nodeID + nodeShift, _noeud())).first;
1513 ma.sommets[ medIndex.next() ] = node;
1516 ma.setOrdre( elemIds[ i ] );
1517 groupe.mailles[i] = imed.insert(ma);
1519 // store nb nodes in polyhedron faces
1520 if ( !nbFaceNodes.empty() ) {
1521 const int* nbFaceNodesPtr = & nbFaceNodes[0];
1522 for ( int i = 0; i < nbElems; ++i ) {
1523 vector<int> & nbNodesByFace = imed.polyherdalNbFaceNodes[ &groupe.maille( i ) ];
1524 nbNodesByFace.assign( nbFaceNodesPtr, nbFaceNodesPtr + nbElemFaces[ i ] );
1525 nbFaceNodesPtr += nbElemFaces[ i ];
1528 // create subPart for "coordinates"
1529 if ( !haveCoords ) {
1530 _SubPart & coordSubPart = imed._subPartDescribed[ coordDesc ];
1531 coordSubPart.myFirstNode = points.insert( make_pair( 1 + nodeShift, _noeud())).first;
1533 // add subPart group to part group
1534 int groupeIndex = imed.groupes.size();
1535 partGroupe->groupes.push_back( groupeIndex );
1538 _SubPart subPart( partNum, partType._name );
1539 subPart.myNbCells = nbElems;
1540 subPart.myCellGroupIndex = groupeIndex;
1541 imed.addSubPart( subPart );
1543 else if ( word == "coordinates" )
1545 // Local node coordinates of a part
1546 // ------------------------------------
1547 int nbNodes = geoFile.getInt(); // nn
1550 vector<int> nodeIds;
1551 if ( haveNodeIds ) {
1552 nodeIds.reserve( nbNodes );
1553 while ( nodeIds.size() < nbNodes )
1554 nodeIds.push_back( geoFile.getInt() ); // id_n
1557 // find out if "coordinates" has already been add at reading connectivity
1558 _SubPartDesc coordDesc( partNum , "coordinates");
1559 map< _SubPartDesc, _SubPart >::iterator descPart =
1560 imed._subPartDescribed.find( coordDesc );
1561 bool haveCoords = ( descPart != imed._subPartDescribed.end() );
1564 // check that all nodes have been added
1565 firstNode = descPart->second.myFirstNode;
1566 descPart->second.myNbNodes = nbNodes;
1567 INoeud inoeud = firstNode, inoEnd = points.end();
1568 int id = inoeud->first, idEnd = id + nbNodes;
1569 for ( ; id < idEnd; ++id ) {
1570 if ( inoeud == inoEnd || inoeud->first > id ) {
1571 INoeud in = points.insert( inoeud, make_pair( id, _noeud() ));
1572 in->second.number = id;
1573 in->second.coord.resize( SPACE_DIM );
1581 int nodeShift = points.empty() ? 0 : points.rbegin()->first;
1582 for ( int iNode = 1; iNode <= nbNodes; ++iNode ) {
1583 INoeud inoeud = points.insert( points.end(), make_pair( iNode + nodeShift, _noeud()));
1584 inoeud->second.number = inoeud->first;
1585 inoeud->second.coord.resize( SPACE_DIM );
1587 firstNode = points.find( 1 + nodeShift );
1588 // create "coordinates" subPart
1589 _SubPart & subPart = imed._subPartDescribed[ coordDesc ];
1590 subPart.myNbNodes = nbNodes;
1591 subPart.myFirstNode = firstNode;
1594 // read coordinates in no interlace mode
1595 INoeud endNode = points.end();
1596 for ( int j = 0; j < SPACE_DIM; ++j ) {
1597 for ( INoeud in = firstNode; in != endNode; ++in ) {
1598 _noeud & node = in->second;
1599 node.coord[ j ] = geoFile.getReal();
1603 else if ( word == "part" )
1605 // Another part encounters
1606 // -----------------------
1607 partNum = geoFile.getInt();
1609 geoFile.toNextLine();
1611 string partName = geoFile.getLine();
1612 if ( partName.empty() )
1613 partName = "Part_" + restLine;
1615 if ( imed.groupes.capacity() - imed.groupes.size() < theMaxNbTypes )
1616 imed.groupes.reserve( size_t( 1.5 * imed.groupes.size() ));
1617 imed.groupes.push_back(_groupe());
1618 partGroupe = & imed.groupes.back();
1619 partGroupe->nom = partName;
1620 partGroupe->groupes.reserve( theMaxNbTypes );
1622 else if ( word == "block" )
1625 // ------------------
1626 bool rectilinear = ( restLine.find( "rectilinear" ) != restLine.npos );
1627 bool uniform = ( restLine.find( "uniform" ) != restLine.npos );
1628 bool curvilinear = ( !rectilinear && !uniform );
1629 bool iblanked = ( restLine.find( "iblanked" ) != restLine.npos );
1630 bool with_ghost = ( restLine.find( "with_ghost" ) != restLine.npos );
1631 bool range = ( restLine.find( "range" ) != restLine.npos );
1634 int I = geoFile.getInt();
1635 int J = geoFile.getInt();
1636 int K = geoFile.getInt();
1637 int NumberOfNodes = I*J*K;
1638 if ( !NumberOfNodes ) continue;
1642 vector<int> ijkRange; // imin imax jmin jmax kmin kmax
1643 ijkRange.reserve(6);
1644 while ( ijkRange.size() < 6 )
1645 ijkRange.push_back( geoFile.getInt() );
1646 I = ijkRange[1]-ijkRange[0]+1;
1647 J = ijkRange[3]-ijkRange[2]+1;
1648 K = ijkRange[5]-ijkRange[4]+1;
1649 NumberOfNodes = I*J*K;
1652 int nodeShift = points.empty() ? 0 : points.rbegin()->first;
1653 for ( int iNode = 1; iNode <= NumberOfNodes; ++iNode ) {
1654 INoeud inoeud = points.insert( points.end(), make_pair( iNode + nodeShift, _noeud()));
1655 _noeud & node = inoeud->second;
1656 node.number = inoeud->first;
1657 node.coord.resize( SPACE_DIM );
1659 INoeud firstNode = points.find( nodeShift + 1 );
1660 INoeud endNode = points.end();
1662 GRID grid; // calculator of unstructured data
1663 grid._iArrayLength = I;
1664 grid._jArrayLength = J;
1665 grid._kArrayLength = K;
1666 grid._numberOfNodes = NumberOfNodes ;
1667 grid._spaceDimension = SPACE_DIM;
1668 if ( J < 2 ) { grid._spaceDimension--; grid._jArrayLength = 0; }
1669 if ( K < 2 ) { grid._spaceDimension--; grid._kArrayLength = 0; }
1670 int nbElems = grid.getNumberOfElements(MED_CELL, MED_ALL_ELEMENTS);
1672 if ( curvilinear ) // read coordinates for all nodes
1674 for ( int j = 0; j < SPACE_DIM; ++j ) {
1675 for ( INoeud in = firstNode; in != endNode; ++in )
1676 in->second.coord[ j ] = geoFile.getReal();
1678 grid._gridType = MED_BODY_FITTED;
1680 else if ( rectilinear ) // read delta vectors with non-regular spacing
1682 grid._iArray = (double*)geoFile.convertReals<double>( I );
1683 grid._jArray = (double*)geoFile.convertReals<double>( J );
1684 grid._kArray = (double*)geoFile.convertReals<double>( K );
1685 grid._gridType = MED_CARTESIAN;
1687 else // uniform: read grid origine and delta vectors for regular spacing grid
1689 TDblOwner xyzOrigin( (double*)geoFile.convertReals<double>( 3 ));
1690 TDblOwner xyzDelta ( (double*)geoFile.convertReals<double>( 3 ));
1691 // compute full delta vectors
1692 grid._iArray = new double[ I ];
1693 grid._jArray = new double[ J ];
1694 grid._kArray = new double[ K ];
1695 double* coors[SPACE_DIM] = { grid._iArray, grid._jArray, grid._kArray };
1696 int size [SPACE_DIM] = { I, J, K };
1697 for ( int j = 0; j < SPACE_DIM; ++j ) {
1698 double* coo = coors[ j ];
1699 double* cooEnd = coo + size[ j ];
1700 coo[0] = xyzOrigin[ j ];
1701 while ( ++coo < cooEnd )
1702 *coo = coo[-1] + xyzDelta[ j ];
1704 grid._gridType = MED_CARTESIAN;
1709 geoFile.skip( NumberOfNodes, /*nbPerLine =*/ 1, INT_WIDTH_GOLD);
1712 geoFile.getWord(); // "ghost_flags"
1713 geoFile.skip( nbElems, /*nbPerLine =*/ 1, INT_WIDTH_GOLD);
1716 if ( haveNodeIds && geoFile.lookAt( "node_ids" )) {
1717 geoFile.getWord(); // "node_ids"
1718 geoFile.skip( NumberOfNodes, /*nbPerLine =*/ 1, INT_WIDTH_GOLD);
1721 if ( haveElemIds && geoFile.lookAt( "element_ids" ) ) {
1722 geoFile.getWord(); // "element_ids"
1723 geoFile.skip( nbElems, /*nbPerLine =*/ 1, INT_WIDTH_GOLD);
1726 if ( !curvilinear ) // let GRID compute all coordinates
1728 grid._coordinate = new COORDINATE;
1729 const double * coo = grid.getCoordinates(MED_FULL_INTERLACE);
1730 typedef _ValueIterator< double > TCoordIt;
1731 TCoordIt xCoo( coo+0, grid._spaceDimension);
1732 TCoordIt yCoo( coo+1, grid._spaceDimension);
1733 TCoordIt zCoo( coo+2, grid._spaceDimension);
1734 if ( grid._spaceDimension < 3 ) zCoo = TCoordIt( grid._kArray, 0 );
1735 if ( grid._spaceDimension < 2 ) yCoo = TCoordIt( grid._jArray, 0 );
1736 for ( INoeud in = firstNode; in != endNode; ++in ) {
1737 _noeud& node = in->second;
1738 node.coord[ 0 ] = xCoo.next();
1739 node.coord[ 1 ] = yCoo.next();
1740 node.coord[ 2 ] = zCoo.next();
1744 // let GRID calculate connectivity
1746 const int * conn = grid.getConnectivity( MED_FULL_INTERLACE, MED_NODAL,
1747 MED_CELL, MED_ALL_ELEMENTS );
1748 medGeometryElement elemType = grid.getElementType( MED_CELL, 1 );
1749 int nbElemNodes = elemType % 100;
1751 partGroupe->mailles.resize( nbElems );
1752 _maille ma( elemType, nbElemNodes );
1753 ma.sommets.resize( nbElemNodes );
1755 for ( int i = 0, nIndex = 0; i < nbElems; ++i ) {
1756 for ( int n = 0; n < nbElemNodes; ++n ) {
1757 int nodeID = conn[ nIndex++ ];
1758 ma.sommets[n] = points.find( nodeID + nodeShift );
1760 //ma.ordre = ++order;
1761 partGroupe->mailles[i] = imed.insert(ma);
1764 _SubPart subPart( partNum, "block" );
1765 subPart.myNbCells = nbElems;
1766 subPart.myNbNodes = NumberOfNodes;
1767 subPart.myFirstNode = firstNode;
1768 subPart.myCellGroupIndex = imed.groupes.size();
1769 imed.addSubPart( subPart );
1774 ( LOCALIZED( STRING(LOC) << "Unexpected word: " << word <<
1775 " in " << getDataFileName()));
1777 } // while ( !geoFile.eof() )
1780 imed.mergeNodesAndElements(TOLERANCE);
1785 //================================================================================
1787 * \brief Read mesh in Gold Binary format
1789 //================================================================================
1791 void ENSIGHT_MESH_RDONLY_DRIVER::readGoldBinary(_InterMed & imed)
1793 const char * LOC ="ENSIGHT_MESH_RDONLY_DRIVER::readGoldBinary() : ";
1796 _BinaryFileReader geoFile( getDataFileName() );
1798 // check if swapping bytes needed
1800 countPartsBinary( geoFile, isSingleFileMode());
1803 geoFile.swapBytes();
1806 if ( getIndexInDataFile() <= 1 )
1808 if ( geoFile.getPosition() == 0 ) {
1809 TStrOwner format( geoFile.getLine() ); // "C|Fortran Binary"
1810 if ( !contains( "C Binary", format )) {
1811 if ( contains( "Fortran Binary", format ))
1812 throw(MEDEXCEPTION(STRING(LOC) << "Fortran Binary format not supported"));
1814 throw(MEDEXCEPTION(STRING(LOC) << "unexpected line in " << getDataFileName()
1815 << "\n" << format.myValues));
1818 if ( isSingleFileMode() ) {
1819 // one time step may be skipped by countPartsBinary
1820 int curTimeStep = geoFile.getPosition() ? 2 : 1 ;
1821 while ( curTimeStep < getIndexInDataFile() ) {
1822 countPartsBinary( geoFile, true ); // skip time step
1826 TStrOwner line( geoFile.getLine() );
1827 if ( isTimeStepBeginning( line.myValues ))
1831 // ----------------------
1832 // Read mesh description
1833 // ----------------------
1835 TStrOwner descriptionLine1 ( geoFile.getLine() );
1836 TStrOwner descriptionLine2 ( geoFile.getLine() );
1838 // find out if the file was created by MED driver
1839 _isMadeByMed = contains( theDescriptionPrefix, descriptionLine1 );
1842 _ptrMesh->setDescription( descriptionLine2.myValues );
1844 _ptrMesh->setDescription( string(descriptionLine1) + descriptionLine2.myValues );
1847 // ----------------------------------------
1848 // Find out presence of node/elem numbers
1849 // ----------------------------------------
1851 // EnSight User Manual (for v8) says:
1852 // You do not have to assign node IDs. If you do, the element connectivities are
1853 // based on the node numbers. If you let EnSight assign the node IDs, the nodes
1854 // are considered to be sequential starting at node 1, and element connectivity is
1855 // done accordingly. If node IDs are set to off, they are numbered internally;
1856 // however, you will not be able to display or query on them. If you have node
1857 // IDs in your data, you can have EnSight ignore them by specifying "node id
1858 // ignore." Using this option may reduce some of the memory taken up by the
1859 // Client and Server, but display and query on the nodes will not be available.
1861 // read "node|element id <off|given|assign|ignore>"
1862 bool haveNodeIds, haveElemIds;
1864 TStrOwner nodeIds( geoFile.getLine() );
1865 TStrOwner elemIds( geoFile.getLine() );
1867 haveNodeIds = ( contains( "given", nodeIds ) || contains( "ignore", nodeIds ) );
1868 haveElemIds = ( contains( "given", elemIds ) || contains( "ignore", elemIds ) );
1871 typedef map<int,_noeud>::iterator INoeud;
1872 map<int,_noeud> & points = imed.points;
1875 _groupe * partGroupe = 0;
1876 int partNum = 0, nbParts = 0;
1878 TFltOwner extents(0); // extents: xmin xmax ymin ymax zmin zmax
1880 while ( !geoFile.eof() )
1882 TStrOwner line( geoFile.getLine() );
1883 if ( isSingleFileMode() && isTimeStepEnd( line.myValues ))
1885 string word, restLine;
1886 TStrTool::split( line.myValues, word, restLine );
1888 const TEnSightElemType & partType = getEnSightType( word );
1889 if ( partType._medType != MED_ALL_ELEMENTS )
1891 // Unstructured element type encounters
1892 // --------------------------------------
1893 int nbElems = *TIntOwner( geoFile.getInt(1) ); // ne
1894 int nbElemNodes = partType._medType % 100;
1895 bool isGhost = isGhostType( word );
1896 int nodeShift = points.empty() ? 0 : points.rbegin()->first;
1899 TIntOwner elemIds( haveElemIds ? geoFile.getInt( nbElems ): 0 ); // id_e
1901 if ( isGhost ) { // do not store ghost elements (?)
1902 int nbNodes = nbElems * nbElemNodes;
1903 if ( partType._name == "nsided" ) // polygons
1905 TIntOwner nbNodesInFace( geoFile.getInt( nbElems ));
1906 nbNodes = std::accumulate( nbNodesInFace.myValues, nbNodesInFace.myValues + nbElems, 0 );
1908 else if ( partType._name == "nfaced" ) // polyhedrons
1910 TIntOwner nbElemFaces( geoFile.getInt( nbElems ));
1911 int nbFaces = accumulate( nbElemFaces.myValues, nbElemFaces.myValues + nbElems, 0 );
1912 TIntOwner nbNodesInFace( geoFile.getInt( nbFaces ));
1913 nbNodes = std::accumulate( nbNodesInFace.myValues, nbNodesInFace.myValues + nbElems, 0 );
1915 geoFile.skip( nbNodes * sizeof(int) );
1919 // add a group corresponding to subPart (geoType)
1920 imed.groupes.push_back(_groupe());
1921 _groupe & groupe = imed.groupes.back();
1922 groupe.mailles.resize( nbElems );
1924 // find out if "coordinates" has already been encountered
1925 _SubPartDesc coordDesc( partNum , "coordinates");
1926 map< _SubPartDesc, _SubPart >::iterator descPart =
1927 imed._subPartDescribed.find( coordDesc );
1928 bool haveCoords = ( descPart != imed._subPartDescribed.end() );
1930 firstNode = descPart->second.myFirstNode;
1931 nodeShift -= descPart->second.myNbNodes;
1934 // read poly element data
1935 bool isPoly = ( !nbElemNodes );
1937 TIntOwner nbElemNodesVec(0), nbElemFaces(0), nbFaceNodes(0);
1938 if ( partType._name == "nsided" ) // polygons
1940 nbElemNodesVec.myValues = geoFile.getInt( nbElems ); // npi
1941 nbNodes = accumulate( nbElemNodesVec.myValues, nbElemNodesVec.myValues + nbElems, 0 );
1943 else if ( partType._name == "nfaced" ) // polyhedrons
1945 nbElemFaces.myValues = geoFile.getInt( nbElems ); // nf_ei
1946 int totalNbFaces = accumulate( nbElemFaces.myValues, nbElemFaces.myValues + nbElems, 0 );
1948 nbFaceNodes.myValues = geoFile.getInt( totalNbFaces ); // np(f_ei)
1949 // calculate nb of nodes in each polyhedron
1950 int* nbEN = nbElemNodesVec.myValues = new int[ nbElems ];
1951 const int *nbFN = nbFaceNodes, *nbEF = nbElemFaces, *nbEND = nbEN + nbElems;
1952 for ( ; nbEN < nbEND; ++nbEN, ++nbEF ) {
1953 nbNodes += *nbEN = accumulate( nbFN, nbFN + *nbEF, 0 );
1957 else // standard types
1959 nbElemNodesVec.myValues = new int[ 1 ];
1960 nbElemNodesVec[ 0 ] = nbElemNodes;
1961 nbNodes = nbElems * nbElemNodes;
1963 // iterator returning nbElemNodes for standard elems and
1964 // next value from nbElemNodesVec for poly elements
1965 _ValueIterator<int> nbElemNodesIt( nbElemNodesVec, isPoly ? 1 : 0);
1967 // iterator returning values form partType._medIndex for standard elems
1968 // and node index (n) for poly elements
1970 _ValueIterator<int> medIndexIt( isPoly ? &n : &partType._medIndex[0],
1972 // read connectivity
1973 _maille ma( partType._medType, nbElemNodes );
1974 ma.sommets.resize( nbElemNodes );
1975 TIntOwner connectivity( geoFile.getInt( nbNodes )); // nn_ei
1976 int* nodeID = connectivity;
1978 for ( int i = 0; i < nbElems; ++i ) {
1979 _ValueIterator<int> medIndex = medIndexIt;
1980 nbElemNodes = nbElemNodesIt.next();
1981 if ( ma.sommets.size() != nbElemNodes )
1982 ma.sommets.resize( nbElemNodes );
1983 for ( n = 0; n < nbElemNodes; ++n, ++nodeID ) {
1985 node = points.find( *nodeID + nodeShift );
1987 node = points.insert( make_pair( *nodeID + nodeShift, _noeud())).first;
1988 ma.sommets[ medIndex.next() ] = node;
1991 ma.setOrdre( elemIds[ i ] );
1992 groupe.mailles[i] = imed.insert(ma);
1994 // store nb nodes in polyhedron faces
1995 if ( nbFaceNodes.myValues ) {
1996 const int* nbFaceNodesPtr = nbFaceNodes.myValues;
1997 for ( int i = 0; i < nbElems; ++i ) {
1998 vector<int> & nbNodesByFace = imed.polyherdalNbFaceNodes[ &groupe.maille( i ) ];
1999 nbNodesByFace.assign( nbFaceNodesPtr, nbFaceNodesPtr + nbElemFaces[ i ] );
2000 nbFaceNodesPtr += nbElemFaces[ i ];
2003 // create subPart for "coordinates"
2004 if ( !haveCoords ) {
2005 _SubPart & coordSubPart = imed._subPartDescribed[ coordDesc ];
2006 coordSubPart.myFirstNode = points.insert( make_pair( 1 + nodeShift, _noeud())).first;
2008 // add subPart group to part group
2009 int groupeIndex = imed.groupes.size();
2010 partGroupe->groupes.push_back( groupeIndex );
2013 _SubPart subPart( partNum, partType._name );
2014 subPart.myNbCells = nbElems;
2015 subPart.myCellGroupIndex = groupeIndex;
2016 imed.addSubPart( subPart );
2018 else if ( word == "coordinates" )
2020 // Local node coordinates of a part
2021 // ------------------------------------
2022 int nbNodes = *TIntOwner( geoFile.getInt(1) ); // nn
2025 TIntOwner nodeIds(0);
2027 nodeIds.myValues = geoFile.getInt( nbNodes ); // id_n
2029 // find out if "coordinates" has already been add at reading connectivity
2030 _SubPartDesc coordDesc( partNum , "coordinates");
2031 map< _SubPartDesc, _SubPart >::iterator descPart =
2032 imed._subPartDescribed.find( coordDesc );
2033 bool haveCoords = ( descPart != imed._subPartDescribed.end() );
2036 // check that all nodes have been added
2037 firstNode = descPart->second.myFirstNode;
2038 descPart->second.myNbNodes = nbNodes;
2039 INoeud inoeud = firstNode, inoEnd = points.end();
2040 int id = inoeud->first, idEnd = id + nbNodes;
2041 for ( ; id < idEnd; ++id ) {
2042 if ( inoeud == inoEnd || inoeud->first > id ) {
2043 INoeud in = points.insert( inoeud, make_pair( id, _noeud() ));
2044 in->second.number = id;
2052 int nodeShift = points.empty() ? 0 : points.rbegin()->first;
2053 for ( int iNode = 1; iNode <= nbNodes; ++iNode ) {
2054 INoeud inoeud = points.insert( points.end(), make_pair( iNode + nodeShift, _noeud()));
2055 inoeud->second.number = inoeud->first;
2057 firstNode = points.find( 1 + nodeShift );
2058 // create "coordinates" subPart
2059 _SubPart & subPart = imed._subPartDescribed[ coordDesc ];
2060 subPart.myNbNodes = nbNodes;
2061 subPart.myFirstNode = firstNode;
2064 // read coordinates in no interlace mode
2065 TFltOwner noInterlaceCoords( geoFile.getFlt( nbNodes * SPACE_DIM ));
2066 float* x = noInterlaceCoords;
2067 float* y = x + nbNodes;
2068 float* z = y + nbNodes;
2069 INoeud endNode = points.end();
2070 for ( INoeud in = firstNode; in != endNode; ++in ) {
2071 _noeud & node = in->second;
2072 node.coord.resize( SPACE_DIM );
2073 node.coord[ 0 ] = *x++;
2074 node.coord[ 1 ] = *y++;
2075 node.coord[ 2 ] = *z++;
2078 else if ( word == "part" )
2080 // Another part encounters
2081 // -----------------------
2082 partNum = *TIntOwner( geoFile.getInt(1) );
2085 string partName( TStrOwner( geoFile.getLine() ));
2086 if ( partName.empty() )
2087 partName = "Part_" + restLine;
2089 if ( imed.groupes.capacity() - imed.groupes.size() < theMaxNbTypes )
2090 imed.groupes.reserve( size_t( 1.5 * imed.groupes.size() ));
2091 imed.groupes.push_back(_groupe());
2092 partGroupe = & imed.groupes.back();
2093 partGroupe->nom = partName;
2094 partGroupe->groupes.reserve( theMaxNbTypes );
2096 else if ( word == "block" )
2099 // ------------------
2100 bool rectilinear = ( restLine.find( "rectilinear" ) != restLine.npos );
2101 bool uniform = ( restLine.find( "uniform" ) != restLine.npos );
2102 bool curvilinear = ( !rectilinear && !uniform );
2103 bool iblanked = ( restLine.find( "iblanked" ) != restLine.npos );
2104 bool with_ghost = ( restLine.find( "with_ghost" ) != restLine.npos );
2105 bool range = ( restLine.find( "range" ) != restLine.npos );
2108 TIntOwner ijk( geoFile.getInt(3) );
2112 int NumberOfNodes = I*J*K;
2113 if ( !NumberOfNodes ) continue;
2117 TIntOwner ijkRange( geoFile.getInt( 6 ));// imin imax jmin jmax kmin kmax
2118 I = ijkRange[1]-ijkRange[0]+1;
2119 J = ijkRange[3]-ijkRange[2]+1;
2120 K = ijkRange[5]-ijkRange[4]+1;
2121 NumberOfNodes = I*J*K;
2124 int nodeShift = points.empty() ? 0 : points.rbegin()->first;
2125 for ( int iNode = 1; iNode <= NumberOfNodes; ++iNode ) {
2126 INoeud inoeud = points.insert( points.end(), make_pair( iNode + nodeShift, _noeud()));
2127 _noeud & node = inoeud->second;
2128 node.number = inoeud->first;
2129 node.coord.resize( SPACE_DIM );
2131 INoeud firstNode = points.find( nodeShift + 1 );
2132 INoeud endNode = points.end();
2134 GRID grid; // calculator of unstructured data
2135 grid._iArrayLength = I;
2136 grid._jArrayLength = J;
2137 grid._kArrayLength = K;
2138 grid._numberOfNodes = NumberOfNodes ;
2139 grid._spaceDimension = SPACE_DIM;
2140 if ( J < 2 ) { grid._spaceDimension--; grid._jArrayLength = 0; }
2141 if ( K < 2 ) { grid._spaceDimension--; grid._kArrayLength = 0; }
2142 int nbElems = grid.getNumberOfElements(MED_CELL, MED_ALL_ELEMENTS);
2144 if ( curvilinear ) // read coordinates for all nodes
2146 TFltOwner noInterlaceCoords( geoFile.getFlt( NumberOfNodes * SPACE_DIM ));
2147 float* x = noInterlaceCoords;
2148 float* y = x + NumberOfNodes;
2149 float* z = y + NumberOfNodes;
2150 for ( INoeud in = firstNode; in != endNode; ++in ) {
2151 _noeud & node = in->second;
2152 node.coord.resize( SPACE_DIM );
2153 node.coord[ 0 ] = *x++;
2154 node.coord[ 1 ] = *y++;
2155 node.coord[ 2 ] = *z++;
2157 grid._gridType = MED_BODY_FITTED;
2159 else if ( rectilinear ) // read delta vectors with non-regular spacing
2161 grid._iArray = (double*)geoFile.convertReals<double>( I );
2162 grid._jArray = (double*)geoFile.convertReals<double>( J );
2163 grid._kArray = (double*)geoFile.convertReals<double>( K );
2164 grid._gridType = MED_CARTESIAN;
2166 else // uniform: read grid origine and delta vectors for regular spacing grid
2168 TFltOwner xyzOrigin( geoFile.getFlt( 3 ));
2169 TFltOwner xyzDelta ( geoFile.getFlt( 3 ));
2170 // compute full delta vectors
2171 grid._iArray = new double[ I ];
2172 grid._jArray = new double[ J ];
2173 grid._kArray = new double[ K ];
2174 double* coors[SPACE_DIM] = { grid._iArray, grid._jArray, grid._kArray };
2175 int size [SPACE_DIM] = { I, J, K };
2176 for ( int j = 0; j < SPACE_DIM; ++j ) {
2177 double* coo = coors[ j ];
2178 double* cooEnd = coo + size[ j ];
2179 coo[0] = xyzOrigin[ j ];
2180 while ( ++coo < cooEnd )
2181 *coo = coo[-1] + xyzDelta[ j ];
2183 grid._gridType = MED_CARTESIAN;
2188 geoFile.skip( NumberOfNodes * sizeof(int) );
2191 TStrOwner( geoFile.getLine() ); // "ghost_flags"
2192 geoFile.skip( nbElems * sizeof(int) );
2195 if ( haveNodeIds && !geoFile.eof() ) {
2196 TStrOwner nextLine( geoFile.getLine() ); // "node_ids"
2197 if ( contains( "node_ids", nextLine ) )
2198 geoFile.skip( NumberOfNodes * sizeof(int) );
2200 geoFile.skip( -MAX_LINE_LENGTH );
2203 TIntOwner elemIdOwner(0);
2204 _ValueIterator<int> elemIds;
2205 if ( haveElemIds && !geoFile.eof() ) {
2206 TStrOwner nextLine( geoFile.getLine() ); // "element_ids"
2207 if ( contains( "element_ids", nextLine ) ) {
2208 elemIdOwner.myValues = geoFile.getInt( nbElems );
2209 elemIds = _ValueIterator<int>( elemIdOwner, 1);
2211 geoFile.skip( -MAX_LINE_LENGTH );
2215 if ( !curvilinear ) // let GRID compute all coordinates
2217 grid._coordinate = new COORDINATE;
2218 const double * coo = grid.getCoordinates(MED_FULL_INTERLACE);
2219 typedef _ValueIterator< double > TCoordIt;
2220 TCoordIt xCoo( coo+0, grid._spaceDimension);
2221 TCoordIt yCoo( coo+1, grid._spaceDimension);
2222 TCoordIt zCoo( coo+2, grid._spaceDimension);
2223 if ( grid._spaceDimension < 3 ) zCoo = TCoordIt( grid._kArray, 0 );
2224 if ( grid._spaceDimension < 2 ) yCoo = TCoordIt( grid._jArray, 0 );
2225 for ( INoeud in = firstNode; in != endNode; ++in ) {
2226 _noeud& node = in->second;
2227 node.coord[ 0 ] = xCoo.next();
2228 node.coord[ 1 ] = yCoo.next();
2229 node.coord[ 2 ] = zCoo.next();
2233 // let GRID calculate connectivity
2235 const int * conn = grid.getConnectivity( MED_FULL_INTERLACE, MED_NODAL,
2236 MED_CELL, MED_ALL_ELEMENTS );
2237 medGeometryElement elemType = grid.getElementType( MED_CELL, 1 );
2238 int nbElemNodes = elemType % 100;
2240 partGroupe->mailles.resize( nbElems );
2241 _maille ma( elemType, nbElemNodes );
2242 ma.sommets.resize( nbElemNodes );
2244 for ( int i = 0, nIndex = 0; i < nbElems; ++i ) {
2245 for ( int n = 0; n < nbElemNodes; ++n ) {
2246 int nodeID = conn[ nIndex++ ];
2247 ma.sommets[n] = points.find( nodeID + nodeShift );
2249 ma.setOrdre( elemIds.next() );
2250 partGroupe->mailles[i] = imed.insert(ma);
2253 _SubPart subPart( partNum, "block" );
2254 subPart.myNbCells = nbElems;
2255 subPart.myNbNodes = NumberOfNodes;
2256 subPart.myFirstNode = firstNode;
2257 subPart.myCellGroupIndex = imed.groupes.size();
2258 imed.addSubPart( subPart );
2260 else if ( word == "extents" )
2262 extents.myValues = geoFile.getFlt( 6 );
2267 ( LOCALIZED( STRING(LOC) << "Unexpected word: " << word <<
2268 " in " << getDataFileName()));
2270 } // while ( !geoFile.eof() )
2273 imed.mergeNodesAndElements(TOLERANCE);
2278 //================================================================================
2280 * \brief Read mesh in Ensight6 ASCII format
2282 //================================================================================
2284 void ENSIGHT_MESH_RDONLY_DRIVER::read6ASCII(_InterMed & imed)
2286 const char * LOC ="ENSIGHT_MESH_RDONLY_DRIVER::read6ASCII() : ";
2289 _ASCIIFileReader geoFile( getDataFileName() );
2291 if ( isSingleFileMode() ) {
2292 int curTimeStep = 1;
2293 while ( curTimeStep < getIndexInDataFile() ) {
2294 while ( !geoFile.isTimeStepEnd())
2298 while ( !geoFile.isTimeStepBeginning() )
2301 // ----------------------
2302 // Read mesh description
2303 // ----------------------
2305 string descriptionLine1 = geoFile.getLine();
2306 string descriptionLine2 = geoFile.getLine();
2308 // find out if the file was created by MED driver
2309 int prefixSize = strlen( theDescriptionPrefix );
2310 _isMadeByMed = ( descriptionLine1.substr(0, prefixSize ) == theDescriptionPrefix );
2313 descriptionLine1 = descriptionLine1.substr( prefixSize );
2314 _ptrMesh->setDescription( descriptionLine1 + descriptionLine2 );
2317 // ----------------------------------------
2318 // Find out presence of node/elem numbers
2319 // ----------------------------------------
2321 // EnSight User Manual (for v8) says:
2322 // You do not have to assign node IDs. If you do, the element connectivities are
2323 // based on the node numbers. If you let EnSight assign the node IDs, the nodes
2324 // are considered to be sequential starting at node 1, and element connectivity is
2325 // done accordingly. If node IDs are set to off, they are numbered internally;
2326 // however, you will not be able to display or query on them. If you have node
2327 // IDs in your data, you can have EnSight ignore them by specifying "node id
2328 // ignore." Using this option may reduce some of the memory taken up by the
2329 // Client and Server, but display and query on the nodes will not be available.
2331 // read "node|element id <off|given|assign|ignore>"
2332 geoFile.getWord(); geoFile.getWord();
2333 string nodeIds = geoFile.getWord();
2334 geoFile.getWord(); geoFile.getWord();
2335 string elemIds = geoFile.getWord();
2337 bool haveNodeIds = ( nodeIds == "given" || nodeIds == "ignore" );
2338 bool haveElemIds = ( elemIds == "given" || elemIds == "ignore" );
2340 map<int,_noeud> & points = imed.points;
2341 typedef map<int,_noeud>::iterator INoeud;
2343 int haveStructuredParts = 0, haveUnstructuredParts = 0;
2345 _groupe * partGroupe = 0;
2348 while ( !geoFile.isTimeStepEnd() )
2350 string word, restLine, line = geoFile.getLine();
2351 TStrTool::split( line, word, restLine );
2353 const TEnSightElemType & partType = getEnSightType( word );
2354 if ( !partType._medIndex.empty() )
2356 // Unstructured element type encounters
2357 // --------------------------------------
2358 int nbElemNodes = partType._medType % 100;
2359 int nbElems = geoFile.getInt();
2361 haveUnstructuredParts++;
2363 imed.groupes.push_back(_groupe());
2364 _groupe & groupe = imed.groupes.back();
2365 groupe.mailles.resize( nbElems );
2367 // read connectivity
2368 _maille ma( partType._medType, nbElemNodes );
2369 ma.sommets.resize( nbElemNodes );
2371 for ( int i = 0; i < nbElems; ++i ) {
2374 for ( int n = 0; n < nbElemNodes; ++n ) {
2375 int nodeID = geoFile.getInt();
2376 ma.sommets[ partType._medIndex[n] ] = points.find( nodeID );
2378 //ma.ordre = ++order;
2379 groupe.mailles[i] = imed.insert(ma);
2382 int groupeIndex = imed.groupes.size();
2383 partGroupe->groupes.push_back( groupeIndex );
2385 _SubPart subPart( partNum, partType._name );
2386 subPart.myNbCells = nbElems;
2387 subPart.myCellGroupIndex = groupeIndex;
2388 imed.addSubPart( subPart );
2390 else if ( word == "part" )
2392 // Another part encounters
2393 // -----------------------
2394 partNum = atoi( restLine.c_str() );
2396 string partName = geoFile.getLine();
2397 if ( partName.empty() )
2398 partName = "Part_" + restLine;
2400 if ( imed.groupes.capacity() - imed.groupes.size() < theMaxNbTypes )
2401 imed.groupes.reserve( size_t( 1.5 * imed.groupes.size() ));
2402 imed.groupes.push_back(_groupe());
2403 partGroupe = & imed.groupes.back();
2404 partGroupe->nom = partName;
2405 partGroupe->groupes.reserve( theMaxNbTypes );
2407 else if ( word == "block" )
2410 // ------------------
2411 bool iblanked = ( restLine == "iblanked" );
2414 int I = geoFile.getInt();
2415 int J = geoFile.getInt();
2416 int K = geoFile.getInt();
2417 int NumberOfNodes = I*J*K;
2418 if ( !NumberOfNodes ) continue;
2419 haveStructuredParts++;
2422 int nodeShift = points.empty() ? 0 : points.rbegin()->first;
2423 for ( int iNode = 1; iNode <= NumberOfNodes; ++iNode ) {
2424 INoeud inoeud = points.insert( points.end(), make_pair( iNode + nodeShift, _noeud()));
2425 _noeud & node = inoeud->second;
2426 node.number = inoeud->first;
2427 node.coord.resize( SPACE_DIM );
2430 INoeud firstNode = points.find( nodeShift + 1 );
2431 INoeud endNode = points.end();
2432 for ( int j = 0; j < SPACE_DIM; ++j ) {
2433 for ( INoeud in = firstNode; in != endNode; ++in ) {
2434 _noeud & node = in->second;
2435 node.coord[ j ] = geoFile.getReal();
2440 geoFile.skip(NumberOfNodes, /*nbPerLine =*/ 10, INT_WIDTH_6);
2442 // let GRID calculate connectivity
2444 grid._numberOfNodes = NumberOfNodes ;
2445 grid._iArrayLength = I;
2446 grid._jArrayLength = J;
2447 grid._kArrayLength = K;
2448 grid._gridType = MED_BODY_FITTED;
2449 grid._spaceDimension= SPACE_DIM;
2450 if ( J < 2 ) { grid._spaceDimension--; grid._jArrayLength = 0; }
2451 if ( K < 2 ) { grid._spaceDimension--; grid._kArrayLength = 0; }
2453 const int * conn = grid.getConnectivity( MED_FULL_INTERLACE, MED_NODAL,
2454 MED_CELL, MED_ALL_ELEMENTS );
2455 medGeometryElement elemType = grid.getElementType( MED_CELL, 1 );
2456 int nbElemNodes = elemType % 100;
2457 int nbElems = grid.getNumberOfElements(MED_CELL, MED_ALL_ELEMENTS);
2459 partGroupe->mailles.resize( nbElems );
2460 _maille ma( elemType, nbElemNodes );
2461 ma.sommets.resize( nbElemNodes );
2463 for ( int i = 0, nIndex = 0; i < nbElems; ++i ) {
2464 for ( int n = 0; n < nbElemNodes; ++n ) {
2465 int nodeID = conn[ nIndex++ ];
2466 ma.sommets[n] = points.find( nodeID + nodeShift );
2468 //ma.ordre = ++order;
2469 partGroupe->mailles[i] = imed.insert(ma);
2472 _SubPart subPart( partNum, "block" );
2473 subPart.myNbCells = nbElems;
2474 subPart.myNbNodes = NumberOfNodes;
2475 subPart.myFirstNode = firstNode;
2476 subPart.myCellGroupIndex = imed.groupes.size();
2477 imed.addSubPart( subPart );
2479 else if ( word == "coordinates" )
2481 // ------------------------------------
2482 // Unstructured global node coordinates
2483 // ------------------------------------
2484 int nbNodes = geoFile.getInt();
2486 cout << "-> loading coordinates of " << nbNodes << " nodes " << endl ;
2489 for ( int i=0 ; i < nbNodes ; i++ )
2491 if ( haveNodeIds ) {
2492 int nodeID = geoFile.getInt();
2493 inoeud = points.insert( make_pair( nodeID, _noeud() )).first;
2494 inoeud->second.number = nodeID;
2498 inoeud = points.insert( points.end(), make_pair( nodeID, _noeud()));
2499 inoeud->second.number = nodeID;
2501 _noeud & node = inoeud->second;
2502 node.coord.resize( SPACE_DIM );
2503 node.coord[ 0 ] = geoFile.getReal();
2504 node.coord[ 1 ] = geoFile.getReal();
2505 node.coord[ 2 ] = geoFile.getReal();
2508 _SubPartDesc cooDesc = _SubPartDesc::globalCoordDesc();
2509 _SubPart subPart( cooDesc.partNumber(), cooDesc.typeName() );
2510 subPart.myNbNodes = nbNodes;
2511 subPart.myFirstNode = points.begin();
2512 imed.addSubPart( subPart );
2517 ( LOCALIZED( STRING(LOC) << "Unexpected word: " << word <<
2518 " in " << getDataFileName()));
2520 } // while ( !geoFile.eof() )
2522 if ( haveStructuredParts && haveUnstructuredParts || haveStructuredParts > 1 )
2523 imed.mergeNodesAndElements(TOLERANCE);
2528 //================================================================================
2530 * \brief Read mesh in Ensight6 ASCII format
2532 //================================================================================
2534 void ENSIGHT_MESH_RDONLY_DRIVER::read6Binary(_InterMed & imed)
2536 const char * LOC ="ENSIGHT_MESH_RDONLY_DRIVER::read6Binary() : ";
2539 _BinaryFileReader geoFile( getDataFileName() );
2541 // check if swapping bytes needed
2543 countPartsBinary( geoFile, isSingleFileMode());
2546 geoFile.swapBytes();
2549 if ( getIndexInDataFile() <= 1 )
2551 if ( geoFile.getPosition() == 0 ) {
2552 TStrOwner format( geoFile.getLine() ); // "C|Fortran Binary"
2553 if ( !contains( "C Binary", format )) {
2554 if ( contains( "Fortran Binary", format ))
2555 throw(MEDEXCEPTION(STRING(LOC) << "Fortran Binary format not supported"));
2557 throw(MEDEXCEPTION(STRING(LOC) << "unexpected line in " << getDataFileName()
2558 << "\n" << format.myValues));
2562 if ( isSingleFileMode() ) {
2563 // one time step may be skipped by countPartsBinary
2564 int curTimeStep = geoFile.getPosition() ? 2 : 1 ;
2565 while ( curTimeStep < getIndexInDataFile() ) {
2566 countPartsBinary( geoFile, true ); // skip time step
2570 TStrOwner line( geoFile.getLine() );
2571 if ( isTimeStepBeginning( line.myValues ))
2575 // ----------------------
2576 // Read mesh description
2577 // ----------------------
2579 TStrOwner descriptionLine1( geoFile.getLine() );
2580 TStrOwner descriptionLine2( geoFile.getLine() );
2582 // find out if the file was created by MED driver
2583 _isMadeByMed = contains( theDescriptionPrefix, descriptionLine1 );
2586 _ptrMesh->setDescription( descriptionLine2.myValues );
2588 _ptrMesh->setDescription( string(descriptionLine1) + descriptionLine2.myValues );
2591 // ----------------------------------------
2592 // Find out presence of node/elem numbers
2593 // ----------------------------------------
2595 // EnSight User Manual (for v8) says:
2596 // You do not have to assign node IDs. If you do, the element connectivities are
2597 // based on the node numbers. If you let EnSight assign the node IDs, the nodes
2598 // are considered to be sequential starting at node 1, and element connectivity is
2599 // done accordingly. If node IDs are set to off, they are numbered internally;
2600 // however, you will not be able to display or query on them. If you have node
2601 // IDs in your data, you can have EnSight ignore them by specifying "node id
2602 // ignore." Using this option may reduce some of the memory taken up by the
2603 // Client and Server, but display and query on the nodes will not be available.
2605 // read "node|element id <off|given|assign|ignore>"
2606 bool haveNodeIds, haveElemIds;
2608 TStrOwner nodeIds( geoFile.getLine() );
2609 TStrOwner elemIds( geoFile.getLine() );
2611 haveNodeIds = ( contains( "given", nodeIds ) || contains( "ignore", nodeIds ) );
2612 haveElemIds = ( contains( "given", elemIds ) || contains( "ignore", elemIds ) );
2614 map<int,_noeud> & points = imed.points;
2615 typedef map<int,_noeud>::iterator INoeud;
2617 int haveStructuredParts = 0, haveUnstructuredParts = 0;
2619 _groupe * partGroupe = 0;
2622 while ( !geoFile.eof() )
2624 TStrOwner line( geoFile.getLine() );
2625 if ( isSingleFileMode() && isTimeStepEnd( line.myValues ))
2627 string word, restLine;
2628 TStrTool::split( line.myValues, word, restLine );
2630 const TEnSightElemType & partType = getEnSightType( word );
2631 if ( !partType._medIndex.empty() )
2633 // Unstructured element type encounters
2634 // --------------------------------------
2635 int nbElemNodes = partType._medType % 100;
2636 int nbElems = *TIntOwner( geoFile.getInt(1) ); // ne
2638 haveUnstructuredParts++;
2640 TIntOwner numbers(0);
2642 numbers.myValues = geoFile.getInt( nbElems ); // id_e
2644 imed.groupes.push_back(_groupe());
2645 _groupe & groupe = imed.groupes.back();
2646 groupe.mailles.resize( nbElems );
2648 // read connectivity
2649 _maille ma( partType._medType, nbElemNodes );
2650 ma.sommets.resize( nbElemNodes );
2651 TIntOwner connectivity( geoFile.getInt( nbElems * nbElemNodes ));
2652 int* nodeID = connectivity;
2654 for ( int i = 0; i < nbElems; ++i ) {
2655 for ( int n = 0; n < nbElemNodes; ++n, ++nodeID )
2656 ma.sommets[ partType._medIndex[n] ] = points.find( *nodeID );
2657 //ma.ordre = ++order;
2658 groupe.mailles[i] = imed.insert(ma);
2661 int groupeIndex = imed.groupes.size();
2662 partGroupe->groupes.push_back( groupeIndex );
2664 _SubPart subPart( partNum, partType._name );
2665 subPart.myNbCells = nbElems;
2666 subPart.myCellGroupIndex = groupeIndex;
2667 imed.addSubPart( subPart );
2669 else if ( word == "part" )
2671 // Another part encounters
2672 // -----------------------
2673 partNum = atoi( restLine.c_str() );
2675 string partName( TStrOwner( geoFile.getLine() ));
2676 if ( partName.empty() )
2677 partName = "Part_" + restLine;
2679 if ( imed.groupes.capacity() - imed.groupes.size() < theMaxNbTypes )
2680 imed.groupes.reserve( size_t( 1.5 * imed.groupes.size() ));
2681 imed.groupes.push_back(_groupe());
2682 partGroupe = & imed.groupes.back();
2683 partGroupe->nom = partName;
2684 partGroupe->groupes.reserve( theMaxNbTypes );
2686 else if ( word == "block" )
2689 // ------------------
2690 bool iblanked = ( restLine == "iblanked" );
2693 TIntOwner ijk( geoFile.getInt(3) );
2697 int NumberOfNodes = I*J*K;
2698 if ( !NumberOfNodes ) continue;
2699 haveStructuredParts++;
2702 int nodeShift = points.empty() ? 0 : points.rbegin()->first;
2704 TFltOwner noInterlaceCoords( geoFile.getFlt( NumberOfNodes * SPACE_DIM ));
2705 float* x = noInterlaceCoords;
2706 float* y = x + NumberOfNodes;
2707 float* z = y + NumberOfNodes;
2708 for ( int iNode = 1; iNode <= NumberOfNodes; ++iNode ) {
2709 INoeud inoeud = points.insert( points.end(), make_pair( iNode + nodeShift, _noeud()));
2710 _noeud & node = inoeud->second;
2711 node.number = inoeud->first;
2712 node.coord.resize( SPACE_DIM );
2713 node.coord[0] = *x++;
2714 node.coord[1] = *y++;
2715 node.coord[2] = *z++;
2720 geoFile.skip(NumberOfNodes * sizeof(int));
2722 // let GRID calculate connectivity
2724 grid._numberOfNodes = NumberOfNodes ;
2725 grid._iArrayLength = I;
2726 grid._jArrayLength = J;
2727 grid._kArrayLength = K;
2728 grid._gridType = MED_BODY_FITTED;
2729 grid._spaceDimension= SPACE_DIM;
2730 if ( J < 2 ) { grid._spaceDimension--; grid._jArrayLength = 0; }
2731 if ( K < 2 ) { grid._spaceDimension--; grid._kArrayLength = 0; }
2733 const int * conn = grid.getConnectivity( MED_FULL_INTERLACE, MED_NODAL,
2734 MED_CELL, MED_ALL_ELEMENTS );
2735 medGeometryElement elemType = grid.getElementType( MED_CELL, 1 );
2736 int nbElemNodes = elemType % 100;
2737 int nbElems = grid.getNumberOfElements(MED_CELL, MED_ALL_ELEMENTS);
2739 partGroupe->mailles.resize( nbElems );
2740 _maille ma( elemType, nbElemNodes );
2741 ma.sommets.resize( nbElemNodes );
2743 for ( int i = 0, nIndex = 0; i < nbElems; ++i ) {
2744 for ( int n = 0; n < nbElemNodes; ++n ) {
2745 int nodeID = conn[ nIndex++ ];
2746 ma.sommets[n] = points.find( nodeID + nodeShift );
2748 //ma.ordre = ++order;
2749 partGroupe->mailles[i] = imed.insert(ma);
2752 _SubPart subPart( partNum, "block" );
2753 subPart.myNbCells = nbElems;
2754 subPart.myNbNodes = NumberOfNodes;
2755 subPart.myFirstNode = points.find( nodeShift + 1 );
2756 subPart.myCellGroupIndex = imed.groupes.size();
2757 imed.addSubPart( subPart );
2759 else if ( word == "coordinates" )
2761 // ------------------------------
2762 // Unstructured node coordinates
2763 // ------------------------------
2764 int nbNodes = *TIntOwner( geoFile.getInt(1) );
2766 TIntOwner numbers(0);
2768 numbers.myValues = geoFile.getInt( nbNodes );
2770 TFltOwner fullInterlaceCoords( geoFile.getFlt( nbNodes * SPACE_DIM ));
2771 float* coord = fullInterlaceCoords;
2773 cout << "-> loading coordinates of " << nbNodes << " nodes " << endl ;
2776 for ( int i=0 ; i < nbNodes ; i++ )
2778 if ( haveNodeIds ) {
2779 int nodeID = numbers[ i ];
2780 inoeud = points.insert( make_pair( nodeID, _noeud() )).first;
2781 inoeud->second.number = nodeID;
2785 inoeud = points.insert( points.end(), make_pair( nodeID, _noeud()));
2786 inoeud->second.number = nodeID;
2788 _noeud & node = inoeud->second;
2789 node.coord.resize( SPACE_DIM );
2790 node.coord[ 0 ] = *coord++;
2791 node.coord[ 1 ] = *coord++;
2792 node.coord[ 2 ] = *coord++;
2795 _SubPartDesc cooDesc = _SubPartDesc::globalCoordDesc();
2796 _SubPart subPart( cooDesc.partNumber(), cooDesc.typeName() );
2797 subPart.myNbNodes = nbNodes;
2798 subPart.myFirstNode = points.begin();
2799 imed.addSubPart( subPart );
2804 ( LOCALIZED( STRING(LOC) << "Unexpected word: " << word <<
2805 " in " << getDataFileName()));
2807 } // while ( !geoFile.eof() )
2809 if ( haveStructuredParts && haveUnstructuredParts || haveStructuredParts > 1 )
2810 imed.mergeNodesAndElements(TOLERANCE);
2815 //================================================================================
2817 * \brief count number of parts in EnSight geometry file
2819 //================================================================================
2821 int ENSIGHT_MESH_RDONLY_DRIVER::countParts(const string& geomFileName,
2822 const bool isSingleFileMode)
2824 const char * LOC ="ENSIGHT_MESH_RDONLY_DRIVER::countParts() : ";
2827 if ( isBinaryDataFile( geomFileName ))
2829 _BinaryFileReader geoFile(geomFileName);
2830 // check if swapping bytes needed
2832 return countPartsBinary( geoFile, isSingleFileMode );
2836 geoFile.swapBytes();
2838 nbParts = countPartsBinary( geoFile, isSingleFileMode );
2842 _ASCIIFileReader geoFile(geomFileName);
2844 if ( isSingleFileMode )
2845 while ( !isTimeStepBeginning( geoFile.getLine() ));
2847 geoFile.getLine(); // description line 1
2848 geoFile.getLine(); // description line 2
2850 // read "node|element id <off|given|assign|ignore>"
2851 geoFile.getWord(); geoFile.getWord();
2852 string nodeIds = geoFile.getWord();
2853 geoFile.getWord(); geoFile.getWord();
2854 string elemIds = geoFile.getWord();
2855 bool haveNodeIds = ( nodeIds == "given" || nodeIds == "ignore" );
2856 bool haveElemIds = ( elemIds == "given" || elemIds == "ignore" );
2859 while ( !geoFile.isTimeStepEnd() )
2861 string word, restLine, line = geoFile.getLine();
2862 TStrTool::split( line, word, restLine );
2864 const TEnSightElemType & partType = getEnSightType( word );
2865 if ( partType._medType != MED_ALL_ELEMENTS )
2867 // Unstructured element type encounters
2868 // --------------------------------------
2869 int nbElemNodes = partType._medType % 100;
2870 int nbElems = geoFile.getInt(); // ne
2873 if ( haveElemIds && isGold )
2874 geoFile.skip( nbElems, /*nbPerLine =*/ 1, INT_WIDTH_GOLD ); // id_e
2876 // skip connectivity
2877 int nbNodes = nbElems * nbElemNodes;
2878 if ( partType._name == "nsided" ) // polygons
2880 for ( int i = 0; i < nbElems; ++i )
2881 nbNodes += geoFile.getInt();
2882 geoFile.skip( nbNodes * INT_WIDTH_GOLD, /*nbLines = */nbElems );
2884 else if ( partType._name == "nfaced" ) // polyhedrons
2887 for ( int i = 0; i < nbElems; ++i )
2888 nbFaces += geoFile.getInt();
2889 for ( int f = 0; f < nbFaces; ++f )
2890 nbNodes += geoFile.getInt();
2891 geoFile.skip( nbNodes * INT_WIDTH_GOLD, /*nbLines = */nbFaces );
2893 else // standard types
2896 geoFile.skip( nbNodes, nbElemNodes, INT_WIDTH_GOLD );
2897 else if ( haveElemIds )
2898 geoFile.skip( nbNodes + nbElems, nbElemNodes+1, INT_WIDTH_6 );
2900 geoFile.skip( nbNodes, nbElemNodes, INT_WIDTH_6 );
2903 else if ( word == "coordinates" )
2905 isGold = ( nbParts > 0 );
2906 int nbNodes = geoFile.getInt(); // nn
2911 geoFile.skip( nbNodes, /*nbPerLine =*/ 1, INT_WIDTH_GOLD ); // node ids
2912 geoFile.skip( nbNodes * SPACE_DIM, /*nbPerLine =*/ 1, FLT_WIDTH ); // coordinates
2915 int coordWidth = 3 * FLT_WIDTH;
2917 coordWidth += INT_WIDTH_6;
2918 geoFile.skip(nbNodes, /*nbPerLine =*/ 1, coordWidth);
2921 else if ( word == "part" )
2925 geoFile.skip( 1, /*nbPerLine =*/ 1, INT_WIDTH_GOLD ); //part number
2927 geoFile.getWord(); // part number
2928 geoFile.toNextLine();
2929 geoFile.getLine(); // description line
2931 else if ( word == "block" )
2934 // ------------------
2935 bool rectilinear = ( restLine.find( "rectilinear" ) != restLine.npos );
2936 bool uniform = ( restLine.find( "uniform" ) != restLine.npos );
2937 bool curvilinear = ( !rectilinear && !uniform );
2938 bool iblanked = ( restLine.find( "iblanked" ) != restLine.npos );
2939 bool with_ghost = ( restLine.find( "with_ghost" ) != restLine.npos );
2940 bool range = ( restLine.find( "range" ) != restLine.npos );
2943 int I = geoFile.getInt();
2944 int J = geoFile.getInt();
2945 int K = geoFile.getInt();
2946 int nbNodes = I*J*K;
2947 if ( !nbNodes ) continue;
2951 vector<int> ijkRange; // imin imax jmin jmax kmin kmax
2952 ijkRange.reserve(6);
2953 while ( ijkRange.size() < 6 )
2954 ijkRange.push_back( geoFile.getInt() );
2955 I = ijkRange[1]-ijkRange[0]+1;
2956 J = ijkRange[3]-ijkRange[2]+1;
2957 K = ijkRange[5]-ijkRange[4]+1;
2960 int nbElems = (I-1)*(J-1)*(K-1);
2962 if ( curvilinear ) // read coordinates for all nodes
2965 geoFile.skip( nbNodes, /*nbPerLine =*/ 1, FLT_WIDTH );
2967 geoFile.skip( nbNodes * SPACE_DIM, /*nbPerLine =*/ 6, FLT_WIDTH );
2969 else if ( rectilinear ) // read delta vectors with non-regular spacing
2971 geoFile.skip( I + J + K, /*nbPerLine =*/ 1, FLT_WIDTH );
2973 else // uniform: read grid origine and delta vectors for regular spacing grid
2975 geoFile.skip( 6, /*nbPerLine =*/ 1, FLT_WIDTH );
2981 geoFile.skip( nbNodes, /*nbPerLine =*/ 1, INT_WIDTH_GOLD);
2983 geoFile.skip( nbNodes, /*nbPerLine =*/ 10, INT_WIDTH_6 );
2987 geoFile.getWord(); // "ghost_flags"
2988 geoFile.skip( nbElems, /*nbPerLine =*/ 1, INT_WIDTH_GOLD);
2991 if ( haveNodeIds && geoFile.lookAt( "node_ids" )) {
2992 geoFile.getWord(); // "node_ids"
2993 geoFile.skip( nbNodes, /*nbPerLine =*/ 1, INT_WIDTH_GOLD);
2996 if ( haveElemIds && geoFile.lookAt( "element_ids" ) ) {
2997 geoFile.getWord(); // "element_ids"
2998 geoFile.skip( nbElems, /*nbPerLine =*/ 1, INT_WIDTH_GOLD);
3001 else if ( word == "extents" ) {
3002 geoFile.getLine(); geoFile.getLine(); geoFile.getLine();// 3 x 2E12.5
3006 throw MED_EXCEPTION ( LOCALIZED( STRING(LOC) << "Unexpected word: " << word));
3013 //================================================================================
3015 * \brief count number of parts in EnSight geometry file
3017 //================================================================================
3019 int ENSIGHT_MESH_RDONLY_DRIVER::countPartsBinary(_BinaryFileReader& geoFile,
3020 const bool isSingleFileMode)
3022 const char * LOC ="ENSIGHT_MESH_RDONLY_DRIVER::countPartsBinary() : ";
3024 if ( geoFile.getPosition() == 0 ) {
3025 TStrOwner format( geoFile.getLine() ); // "C|Fortran Binary"
3026 if ( !contains( "C Binary", format )) {
3027 if ( contains( "Fortran Binary", format ))
3028 throw(MEDEXCEPTION(STRING(LOC) << "Fortran Binary format not supported"));
3030 throw(MEDEXCEPTION(STRING(LOC) << "unexpected line: \n" << format.myValues));
3034 if ( isSingleFileMode ) {
3036 TStrOwner line( geoFile.getLine() );
3037 if ( isTimeStepBeginning( line.myValues ))
3042 // 2 description lines
3043 // ----------------------
3044 geoFile.skip( 2 * MAX_LINE_LENGTH );
3046 // read "node|element id <off|given|assign|ignore>"
3047 bool haveNodeIds, haveElemIds;
3049 TStrOwner nodeIds( geoFile.getLine() );
3050 TStrOwner elemIds( geoFile.getLine() );
3052 haveNodeIds = ( contains( "given", nodeIds ) || contains( "ignore", nodeIds ) );
3053 haveElemIds = ( contains( "given", elemIds ) || contains( "ignore", elemIds ) );
3056 int nbParts = 0; // the result
3059 while ( !geoFile.eof() )
3061 TStrOwner line( geoFile.getLine() );
3062 if ( isSingleFileMode && isTimeStepEnd( line.myValues ))
3064 string word, restLine;
3065 TStrTool::split( line.myValues, word, restLine );
3067 const TEnSightElemType & partType = getEnSightType( word );
3068 if ( partType._medType != MED_ALL_ELEMENTS )
3070 // Unstructured element type encounters
3071 // --------------------------------------
3072 int nbElems = *TIntOwner( geoFile.getInt(1) ); // ne
3073 int nbElemNodes = partType._medType % 100;
3077 geoFile.skip( nbElems * sizeof(int) ); // id_e
3079 int nbNodes = nbElems * nbElemNodes;
3080 if ( partType._name == "nsided" ) // polygons
3082 TIntOwner nbNodesInFace( geoFile.getInt( nbElems ));
3083 nbNodes = std::accumulate( nbNodesInFace.myValues, nbNodesInFace.myValues + nbElems, 0 );
3085 else if ( partType._name == "nfaced" ) // polyhedrons
3087 TIntOwner nbElemFaces( geoFile.getInt( nbElems ));
3088 int nbFaces = accumulate( nbElemFaces.myValues, nbElemFaces.myValues + nbElems, 0 );
3089 TIntOwner nbNodesInFace( geoFile.getInt( nbFaces ));
3090 nbNodes = std::accumulate( nbNodesInFace.myValues, nbNodesInFace.myValues + nbFaces, 0 );
3092 geoFile.skip( nbNodes * sizeof(int) );
3094 else if ( word == "coordinates" )
3098 int nbNodes = *TIntOwner( geoFile.getInt(1) ); // nn
3102 geoFile.skip( nbNodes * sizeof(int) ); // id_n
3105 geoFile.skip( nbNodes * SPACE_DIM * sizeof(int) );
3107 else if ( word == "part" )
3110 if ( isGold ) geoFile.skip(sizeof(int)); // part #
3112 geoFile.skip(MAX_LINE_LENGTH); // description line
3114 else if ( word == "block" )
3117 // ------------------
3118 bool rectilinear = ( restLine.find( "rectilinear" ) != restLine.npos );
3119 bool uniform = ( restLine.find( "uniform" ) != restLine.npos );
3120 bool curvilinear = ( !rectilinear && !uniform );
3121 bool iblanked = ( restLine.find( "iblanked" ) != restLine.npos );
3122 bool with_ghost = ( restLine.find( "with_ghost" ) != restLine.npos );
3123 bool range = ( restLine.find( "range" ) != restLine.npos );
3126 TIntOwner ijk( geoFile.getInt(3) );
3130 int NumberOfNodes = I*J*K;
3131 if ( !NumberOfNodes ) {
3132 if ( I != 0 && J != 0 && K != 0 )
3133 throw MEDEXCEPTION( "Need to swap bytes" );
3139 TIntOwner ijkRange( geoFile.getInt( 6 ));// imin imax jmin jmax kmin kmax
3140 I = ijkRange[1]-ijkRange[0]+1;
3141 J = ijkRange[3]-ijkRange[2]+1;
3142 K = ijkRange[5]-ijkRange[4]+1;
3143 NumberOfNodes = I*J*K;
3145 int nbElems = (I-1)*(J-1)*(K-1);
3147 if ( curvilinear ) // read coordinates for all nodes
3149 geoFile.skip( NumberOfNodes * SPACE_DIM * sizeof(float) );
3151 else if ( rectilinear ) // read delta vectors with non-regular spacing
3153 geoFile.skip( (I+J+K) * sizeof(float) );
3155 else // uniform: read grid origine and delta vectors for regular spacing grid
3157 geoFile.skip( 6 * sizeof(float) );
3162 geoFile.skip( NumberOfNodes * sizeof(int) );
3165 geoFile.skip( MAX_LINE_LENGTH ); // "ghost_flags"
3166 geoFile.skip( nbElems * sizeof(int) );
3169 if ( haveNodeIds && isGold && !geoFile.eof() ) {
3170 TStrOwner nextLine( geoFile.getLine() ); // "node_ids"
3171 if ( contains( "node_ids", nextLine ) )
3172 geoFile.skip( NumberOfNodes * sizeof(int) );
3174 geoFile.skip( -MAX_LINE_LENGTH );
3177 if ( haveElemIds && isGold && !geoFile.eof() ) {
3178 TStrOwner nextLine( geoFile.getLine() ); // "element_ids"
3179 if ( contains( "element_ids", nextLine ) )
3180 geoFile.skip( nbElems * sizeof(int) );
3182 geoFile.skip( -MAX_LINE_LENGTH );
3185 else if ( word == "extents" )
3187 geoFile.skip( 6 * sizeof(float) );
3191 throw MED_EXCEPTION ( LOCALIZED( STRING(LOC) << "Unexpected word: " << word ));
3193 } // while ( !geoFile.eof() )
3198 GROUP* ENSIGHT_MESH_RDONLY_DRIVER::makeGroup( _groupe& grp,
3201 //const char* LOC = "ENSIGHT_MESH_RDONLY_DRIVER::makeGroup(): error";
3203 // prevent creation of other groups but only this one
3204 for (size_t i=0; i < imed.groupes.size(); ++i)
3205 imed.groupes[i].nom.clear();
3207 // let _intermediateMED create a GROUP from grp
3208 grp.medGroup = 0; // the new GROUP should appear in grp.medGroup
3210 vector<GROUP *> tmp;
3211 imed.getGroups( tmp, tmp, tmp, tmp, imed._medMesh );
3212 if ( !grp.medGroup )
3213 throw MEDEXCEPTION(LOCALIZED("Can't create a GROUP from _groupe"));
3215 grp.medGroup->setName(""); // to let a caller set a proper name
3218 // find pre-existing equal _groupe
3219 _groupe * equalGroupe = 0;
3220 for (unsigned int i=0; i < imed.groupes.size() && !equalGroupe; ++i) {
3221 _groupe& g = imed.groupes[i];
3222 if ( &g != &grp && g.medGroup && g.medGroup->deepCompare( *grp.medGroup ))
3225 if ( equalGroupe ) {
3226 delete grp.medGroup;
3227 grp.medGroup = equalGroupe->medGroup;
3229 else { // new unique group
3231 if ( grp.medGroup->isOnAllElements() ) // on all elems
3232 grp.medGroup->setName( string("SupportOnAll_")+entNames[ grp.medGroup->getEntity() ] );
3234 // add a new group to mesh
3235 if ( !imed._isOwnMedMesh ) {
3236 vector<GROUP*> * groups = 0;
3237 switch ( grp.medGroup->getEntity() ) {
3238 case MED_CELL: groups = & imed._medMesh->_groupCell; break;
3239 case MED_FACE: groups = & imed._medMesh->_groupFace; break;
3240 case MED_EDGE: groups = & imed._medMesh->_groupEdge; break;
3241 case MED_NODE: groups = & imed._medMesh->_groupNode; break;
3245 groups->resize( groups->size() + 1 );
3246 groups->at( groups->size() - 1) = grp.medGroup;
3250 return grp.medGroup;