1 // Copyright (C) 2007-2011 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 // File : DriverCGNS_Write.cxx
23 // Created : Fri Aug 5 17:43:54 2011
24 // Author : Edward AGAPOV (eap)
26 #include "DriverCGNS_Write.hxx"
28 #include "SMDS_MeshNode.hxx"
29 #include "SMDS_VolumeTool.hxx"
30 #include "SMESHDS_GroupBase.hxx"
31 #include "SMESHDS_Mesh.hxx"
32 #include "SMESH_Comment.hxx"
37 #if CGNS_VERSION < 3100
45 //================================================================================
47 * \brief Return interlace and type of CGNS element for the given SMDSAbs_EntityType
49 //================================================================================
51 const int* getInterlaceAndType( const SMDSAbs_EntityType smType,
52 CGNS_ENUMT( ElementType_t ) & cgType )
54 static vector< const int* > interlaces;
55 static vector< CGNS_ENUMT( ElementType_t )> cgTypes;
56 if ( interlaces.empty() )
58 interlaces.resize( SMDSEntity_Last, 0 );
59 cgTypes.resize( SMDSEntity_Last, CGNS_ENUMV( ElementTypeNull ));
61 static int ids[] = {0};
62 interlaces[SMDSEntity_0D] = ids;
63 cgTypes [SMDSEntity_0D] = CGNS_ENUMV( NODE );
66 static int ids[] = { 0, 1 };
67 interlaces[SMDSEntity_Edge] = ids;
68 cgTypes [SMDSEntity_Edge] = CGNS_ENUMV( BAR_2 );
71 static int ids[] = { 0, 1, 2 };
72 interlaces[SMDSEntity_Quad_Edge] = ids;
73 cgTypes [SMDSEntity_Quad_Edge] = CGNS_ENUMV( BAR_3 );
76 static int ids[] = { 0, 2, 1 };
77 interlaces[SMDSEntity_Triangle] = ids;
78 cgTypes [SMDSEntity_Triangle] = CGNS_ENUMV( TRI_3 );
81 static int ids[] = { 0, 2, 1, 5, 4, 3 };
82 interlaces[SMDSEntity_Quad_Triangle] = ids;
83 cgTypes [SMDSEntity_Quad_Triangle] = CGNS_ENUMV( TRI_6 );
86 static int ids[] = { 0, 3, 2, 1 };
87 interlaces[SMDSEntity_Quadrangle] = ids;
88 cgTypes [SMDSEntity_Quadrangle] = CGNS_ENUMV( QUAD_4 );
91 static int ids[] = { 0,3,2,1,7,6,5,4 };
92 interlaces[SMDSEntity_Quad_Quadrangle] = ids;
93 cgTypes [SMDSEntity_Quad_Quadrangle] = CGNS_ENUMV( QUAD_8 );
96 static int ids[] = { 0,3,2,1,7,6,5,4,8 };
97 interlaces[SMDSEntity_BiQuad_Quadrangle] = ids;
98 cgTypes [SMDSEntity_BiQuad_Quadrangle] = CGNS_ENUMV( QUAD_9 );
101 static int ids[] = { 0, 2, 1, 3 };
102 interlaces[SMDSEntity_Tetra] = ids;
103 cgTypes [SMDSEntity_Tetra] = CGNS_ENUMV( TETRA_4 );
106 static int ids[] = { 0,2,1,3,6,5,4,7,9,8 };
107 interlaces[SMDSEntity_Quad_Tetra] = ids;
108 cgTypes [SMDSEntity_Quad_Tetra] = CGNS_ENUMV( TETRA_10 );
111 static int ids[] = { 0,3,2,1,4 };
112 interlaces[SMDSEntity_Pyramid] = ids;
113 cgTypes [SMDSEntity_Pyramid] = CGNS_ENUMV( PYRA_5 );
116 static int ids[] = { 0,3,2,1,4,8,7,6,5,9,12,11,10 };
117 interlaces[SMDSEntity_Quad_Pyramid] = ids;
118 cgTypes [SMDSEntity_Quad_Pyramid] = CGNS_ENUMV( PYRA_13 );
121 static int ids[] = { 0,2,1,3,5,4 };
122 interlaces[SMDSEntity_Penta] = ids;
123 cgTypes [SMDSEntity_Penta] = CGNS_ENUMV( PENTA_6 );
126 static int ids[] = { 0,2,1,3,5,4,8,7,6,9,11,10,14,13,12 };
127 interlaces[SMDSEntity_Quad_Penta] = ids;
128 cgTypes [SMDSEntity_Quad_Penta] = CGNS_ENUMV( PENTA_15 );
131 static int ids[] = { 0,3,2,1,4,7,6,5 };
132 interlaces[SMDSEntity_Hexa] = ids;
133 cgTypes [SMDSEntity_Hexa] = CGNS_ENUMV( HEXA_8 );
136 static int ids[] = { 0,3,2,1,4,7,6,5,11,10,9,8,12,15,14,13,19,18,17,16 };
137 interlaces[SMDSEntity_Quad_Hexa] = ids;
138 cgTypes [SMDSEntity_Quad_Hexa] = CGNS_ENUMV( HEXA_20 );
141 static int ids[] = { 0,3,2,1,4,7,6,5,11,10,9,8,12,15,14,13,19,18,17,16,
142 20, 24,23,22,21, 25};
143 interlaces[SMDSEntity_TriQuad_Hexa] = ids;
144 cgTypes [SMDSEntity_TriQuad_Hexa] = CGNS_ENUMV( HEXA_27 );
147 cgTypes[SMDSEntity_Polygon] = CGNS_ENUMV( NGON_n );
148 cgTypes[SMDSEntity_Polyhedra] = CGNS_ENUMV( NFACE_n );
149 cgTypes[SMDSEntity_Hexagonal_Prism] = CGNS_ENUMV( NFACE_n );
152 cgType = cgTypes[ smType ];
153 return interlaces[ smType ];
156 //================================================================================
158 * \brief Cut off type of boundary condition from the group name
160 //================================================================================
162 CGNS_ENUMT( BCType_t ) getBCType( string& groupName )
164 CGNS_ENUMT( BCType_t ) bcType = CGNS_ENUMV( BCGeneral ); // default type
166 // boundary condition type starts from "BC"
167 size_t bcBeg = groupName.find("BC");
168 if ( bcBeg != string::npos )
170 for ( int t = 0; t < NofValidBCTypes; ++t )
172 CGNS_ENUMT( BCType_t ) type = CGNS_ENUMT( BCType_t)( t );
173 string typeName = cg_BCTypeName( type );
174 if ( typeName == &groupName[0] + bcBeg )
177 while ( bcBeg > 0 && isspace( bcBeg-1 ))
182 groupName = groupName.substr( 0, bcBeg-1 );
189 //================================================================================
191 * \brief Sortable face of a polyhedron
195 int _id; // id of NGON_n
196 vector< int > _nodes; // lowest node IDs used for sorting
198 TPolyhedFace( const SMDS_MeshNode** nodes, const int nbNodes, int ID):_id(ID)
201 for ( int i = 0; i < nbNodes; ++i )
202 ids.insert( nodes[i]->GetID() );
204 _nodes.resize( 3 ); // std::min( nbNodes, 4 )); hope 3 nodes is enough
205 set< int >::iterator idIt = ids.begin();
206 for ( size_t j = 0; j < _nodes.size(); ++j, ++idIt )
209 bool operator< (const TPolyhedFace& o ) const
211 return _nodes < o._nodes;
214 //================================================================================
216 * \brief Return CGNS id of an element
218 //================================================================================
220 cgsize_t cgnsID( const SMDS_MeshElement* elem,
221 const map< const SMDS_MeshElement*, cgsize_t >& elem2cgID )
223 map< const SMDS_MeshElement*, cgsize_t >::const_iterator e2id = elem2cgID.find( elem );
224 return ( e2id == elem2cgID.end() ? elem->GetID() : e2id->second );
229 //================================================================================
231 * \brief Write the mesh into the CGNS file
233 //================================================================================
235 Driver_Mesh::Status DriverCGNS_Write::Perform()
237 myErrorMessages.clear();
239 if ( !myMesh || myMesh->GetMeshInfo().NbElements() < 1 )
240 return addMessage( !myMesh ? "NULL mesh" : "Empty mesh (no elements)", /*fatal = */true );
243 if ( cg_open(myFile.c_str(), CG_MODE_MODIFY, &_fn) != CG_OK &&
244 cg_open(myFile.c_str(), CG_MODE_WRITE, &_fn) != CG_OK )
245 return addMessage( cg_get_error(), /*fatal = */true );
250 const int spaceDim = 3;
252 if ( myMesh->NbFaces() > 0 ) meshDim = 2;
253 if ( myMesh->NbVolumes() > 0 ) meshDim = 3;
255 if ( myMeshName.empty() )
258 if ( cg_nbases( _fn, &nbases) == CG_OK)
259 myMeshName = ( SMESH_Comment("Base_") << nbases+1 );
261 myMeshName = "Base_0";
264 if ( cg_base_write( _fn, myMeshName.c_str(), meshDim, spaceDim, &iBase))
265 return addMessage( cg_get_error(), /*fatal = */true );
270 int nbCells = myMesh->NbEdges();
272 nbCells = myMesh->NbVolumes();
273 else if ( meshDim == 2 )
274 nbCells = myMesh->NbFaces();
276 cgsize_t size[9] = { myMesh->NbNodes(), nbCells, /*NBoundVertex=*/0, 0,0,0,0,0,0 };
278 if ( cg_zone_write( _fn, iBase, "SMESH_Mesh", size,
279 CGNS_ENUMV( Unstructured ), &iZone) != CG_OK )
280 return addMessage( cg_get_error(), /*fatal = */true );
282 // Map to store only elements whose an SMDS ID differs from a CGNS one
283 typedef map< const SMDS_MeshElement*, cgsize_t > TElem2cgIDMap;
284 vector< TElem2cgIDMap > elem2cgIDByEntity( SMDSEntity_Last );
285 TElem2cgIDMap::iterator elem2cgIDIter;
287 TElem2cgIDMap & n2cgID = elem2cgIDByEntity[ SMDSEntity_Node ];
292 vector< double > coords( myMesh->NbNodes() );
295 SMDS_NodeIteratorPtr nIt = myMesh->nodesIterator( /*idInceasingOrder=*/true );
296 for ( int i = 0; nIt->more(); ++i ) coords[i] = nIt->next()->X();
297 if ( cg_coord_write( _fn, iBase, iZone, CGNS_ENUMV(RealDouble),
298 "CoordinateX", &coords[0], &iC) != CG_OK )
299 return addMessage( cg_get_error(), /*fatal = */true );
301 nIt = myMesh->nodesIterator( /*idInceasingOrder=*/true );
302 for ( int i = 0; nIt->more(); ++i ) coords[i] = nIt->next()->Y();
303 if ( cg_coord_write( _fn, iBase, iZone, CGNS_ENUMV(RealDouble),
304 "CoordinateY", &coords[0], &iC) != CG_OK )
305 return addMessage( cg_get_error(), /*fatal = */true );
307 nIt = myMesh->nodesIterator( /*idInceasingOrder=*/true );
308 for ( int i = 0; nIt->more(); ++i ) coords[i] = nIt->next()->Z();
309 if ( cg_coord_write( _fn, iBase, iZone, CGNS_ENUMV(RealDouble),
310 "CoordinateZ", &coords[0], &iC) != CG_OK )
311 return addMessage( cg_get_error(), /*fatal = */true );
313 // store CGNS ids of nodes
314 nIt = myMesh->nodesIterator( /*idInceasingOrder=*/true );
315 for ( int i = 0; nIt->more(); ++i )
317 const SMDS_MeshElement* n = nIt->next();
318 if ( n->GetID() != i+1 )
319 n2cgID.insert( n2cgID.end(), make_pair( n, i+1 ));
325 cgsize_t cgID = 1, startID;
327 // write into a section all successive elements of one geom type
329 vector< cgsize_t > elemData;
330 SMDS_ElemIteratorPtr elemIt = myMesh->elementsIterator();
331 const SMDS_MeshElement* elem = elemIt->next();
334 const SMDSAbs_EntityType elemType = elem->GetEntityType();
335 CGNS_ENUMT( ElementType_t ) cgType;
336 const int* interlace = getInterlaceAndType( elemType, cgType );
338 TElem2cgIDMap & elem2cgID = elem2cgIDByEntity[ elemType ];
343 if ( interlace ) // STANDARD elements
346 for ( int i = 0, nb = elem->NbNodes(); i < nb; ++i )
347 elemData.push_back( cgnsID( elem->GetNode( interlace[i] ), n2cgID ));
348 if ( elem->GetID() != cgID )
349 elem2cgID.insert( elem2cgID.end(), make_pair( elem, cgID ));
351 elem = elemIt->more() ? elemIt->next() : 0;
353 while ( elem && elem->GetEntityType() == elemType );
355 else if ( elemType == SMDSEntity_Polygon ) // POLYGONS
358 elemData.push_back( elem->NbNodes() );
359 for ( int i = 0, nb = elem->NbNodes(); i < nb; ++i )
360 elemData.push_back( cgnsID( elem->GetNode(i), n2cgID ));
361 if ( elem->GetID() != cgID )
362 elem2cgID.insert( elem2cgID.end(), make_pair( elem, cgID ));
364 elem = elemIt->more() ? elemIt->next() : 0;
366 while ( elem && elem->GetEntityType() == elemType );
368 else if ( elemType == SMDSEntity_Polyhedra ||
369 elemType == SMDSEntity_Hexagonal_Prism) // POLYHEDRA
371 // to save polyhedrons after all
372 const SMDS_MeshInfo& meshInfo = myMesh->GetMeshInfo();
373 if ( meshInfo.NbPolyhedrons() == meshInfo.NbElements() - cgID + 1 )
374 break; // only polyhedrons remain
375 while ( elem && elem->GetEntityType() == elemType )
376 elem = elemIt->more() ? elemIt->next() : 0;
380 SMESH_Comment sectionName( cg_ElementTypeName( cgType ));
381 sectionName << " " << startID << " - " << cgID-1;
383 if ( cg_section_write(_fn, iBase, iZone, sectionName.c_str(), cgType, startID,
384 cgID-1, /*nbndry=*/0, &elemData[0], &iSec) != CG_OK )
385 return addMessage( cg_get_error(), /*fatal = */true );
387 // Write polyhedral volumes
388 // -------------------------
390 if ( myMesh->GetMeshInfo().NbElements() != cgID ) // polyhedra or hexagonal prisms remain
392 // the polyhedron (NFACE_n) is described as a set of signed face IDs,
393 // so first we are to write all polygones (NGON_n) bounding polyhedrons
395 vector< cgsize_t > faceData;
396 set< TPolyhedFace > faces;
397 set< TPolyhedFace >::iterator faceInSet;
398 vector<const SMDS_MeshNode *> faceNodesVec;
399 int nbPolygones = 0, faceID;
405 int nbPolyhTreated = 0;
407 TElem2cgIDMap * elem2cgID = 0;
408 TElem2cgIDMap & n2cgID = elem2cgIDByEntity[ SMDSEntity_Node ];
410 SMDS_ElemIteratorPtr elemIt = myMesh->elementsIterator();
411 while ( elemIt->more() )
413 elem = elemIt->next();
414 SMDSAbs_EntityType type = elem->GetEntityType();
415 if ( type == SMDSEntity_Polyhedra ||
416 type == SMDSEntity_Hexagonal_Prism )
420 vol.SetExternalNormal();
421 const int nbFaces = vol.NbFaces();
422 elemData.push_back( nbFaces );
423 for ( int iF = 0; iF < nbFaces; ++iF )
425 const int nbNodes = vol.NbFaceNodes( iF );
426 const SMDS_MeshNode** faceNodes = vol.GetFaceNodes( iF );
427 faceNodesVec.assign( faceNodes, faceNodes + nbNodes );
428 if (( elem = myMesh->FindElement( faceNodesVec, SMDSAbs_Face, /*noMedium=*/false)))
430 // a face of the polyhedron is present in the mesh
431 faceID = cgnsID( elem, elem2cgIDByEntity[ elem->GetEntityType() ]);
433 else if ( vol.IsFreeFace( iF ))
435 // the face is not shared by volumes
438 faceData.push_back( nbNodes );
439 for ( int i = 0; i < nbNodes; ++i )
440 faceData.push_back( cgnsID( faceNodes[i], n2cgID ));
444 TPolyhedFace face( faceNodes, nbNodes, cgID );
445 faceInSet = faces.insert( faces.end(), face );
446 if ( faceInSet->_id == cgID ) // the face encounters for the 1st time
450 faceData.push_back( nbNodes );
451 for ( int i = 0; i < nbNodes; ++i )
452 faceData.push_back( cgnsID( faceNodes[i], n2cgID ));
456 // the face encounters for the 2nd time; we hope it won't encounter once more,
457 // for that we can erase it from the set of faces
458 faceID = -faceInSet->_id;
459 faces.erase( faceInSet );
462 elemData.push_back( faceID );
467 if ( nbPolygones > 0 )
469 if ( cg_section_write(_fn, iBase, iZone, "Faces of Polyhedrons",
470 CGNS_ENUMV( NGON_n ), cgID - nbPolygones, cgID-1,
471 /*nbndry=*/0, &faceData[0], &iSec) != CG_OK )
472 return addMessage( cg_get_error(), /*fatal = */true );
475 if ( cg_section_write(_fn, iBase, iZone, "Polyhedrons",
476 CGNS_ENUMV( NFACE_n ), cgID, cgID+nbPolyhTreated-1,
477 /*nbndry=*/0, &elemData[0], &iSec) != CG_OK )
478 return addMessage( cg_get_error(), /*fatal = */true );
480 if ( !myMesh->GetGroups().empty() )
482 // store CGNS ids of polyhedrons
483 elem2cgID = &elem2cgIDByEntity[ SMDSEntity_Polyhedra ];
484 elemIt = myMesh->elementsIterator();
485 while ( elemIt->more() )
487 elem = elemIt->next();
488 if ( elem->GetEntityType() == SMDSEntity_Polyhedra )
490 if ( elem->GetID() != cgID )
491 elem2cgID->insert( elem2cgID->end(), make_pair( elem, cgID ));
496 } // write polyhedral volumes
499 // Write groups as boundary conditions
500 // ------------------------------------
502 const set<SMESHDS_GroupBase*>& groups = myMesh->GetGroups();
503 set<SMESHDS_GroupBase*>::const_iterator grpIt = groups.begin();
504 set< string > groupNames; groupNames.insert(""); // to avoid duplicated and empty names
505 for ( ; grpIt != groups.end(); ++grpIt )
507 const SMESHDS_GroupBase* group = *grpIt;
509 // write BC location (default is Vertex)
510 CGNS_ENUMT( GridLocation_t ) location = CGNS_ENUMV( Vertex );
511 if ( group->GetType() != SMDSAbs_Node )
515 switch ( group->GetType() ) {
516 case SMDSAbs_Volume: location = CGNS_ENUMV( FaceCenter ); break; // !!!
517 case SMDSAbs_Face: location = CGNS_ENUMV( FaceCenter ); break; // OK
518 case SMDSAbs_Edge: location = CGNS_ENUMV( EdgeCenter ); break; // OK
523 switch ( group->GetType() ) {
524 case SMDSAbs_Face: location = CGNS_ENUMV( FaceCenter ); break; // ???
525 case SMDSAbs_Edge: location = CGNS_ENUMV( EdgeCenter ); break; // OK
530 location = CGNS_ENUMV( EdgeCenter ); break; // ???
535 // try to extract type of boundary condition from the group name
536 string name = group->GetStoreName();
537 CGNS_ENUMT( BCType_t ) bcType = getBCType( name );
538 while ( !groupNames.insert( name ).second )
539 name = (SMESH_Comment( "Group_") << groupNames.size());
541 // write IDs of elements
542 vector< cgsize_t > pnts;
543 pnts.reserve( group->Extent() );
544 SMDS_ElemIteratorPtr elemIt = group->GetElements();
545 while ( elemIt->more() )
547 const SMDS_MeshElement* elem = elemIt->next();
548 pnts.push_back( cgnsID( elem, elem2cgIDByEntity[ elem->GetEntityType() ]));
551 if ( cg_boco_write( _fn, iBase, iZone, name.c_str(), bcType,
552 CGNS_ENUMV( PointList ), pnts.size(), &pnts[0], &iBC) != CG_OK )
553 return addMessage( cg_get_error(), /*fatal = */true);
556 if ( location != CGNS_ENUMV( Vertex ))
558 if ( cg_boco_gridlocation_write( _fn, iBase, iZone, iBC, location) != CG_OK )
559 return addMessage( cg_get_error(), /*fatal = */false);
565 //================================================================================
569 //================================================================================
571 DriverCGNS_Write::DriverCGNS_Write(): _fn(0)
575 //================================================================================
577 * \brief Close the cgns file at destruction
579 //================================================================================
581 DriverCGNS_Write::~DriverCGNS_Write()