]> SALOME platform Git repositories - modules/smesh.git/blob - src/DriverCGNS/DriverCGNS_Write.cxx
Salome HOME
bos #20256 [CEA 18523] Porting SMESH to int 64 bits
[modules/smesh.git] / src / DriverCGNS / DriverCGNS_Write.cxx
1 // Copyright (C) 2007-2021  CEA/DEN, EDF R&D, OPEN CASCADE
2 //
3 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
5 //
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, or (at your option) any later version.
10 //
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.
15 //
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
19 //
20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
21 //
22 // File      : DriverCGNS_Write.cxx
23 // Created   : Fri Aug  5 17:43:54 2011
24 // Author    : Edward AGAPOV (eap)
25
26 #include "DriverCGNS_Write.hxx"
27
28 #include "SMDS_IteratorOnIterators.hxx"
29 #include "SMDS_MeshNode.hxx"
30 #include "SMDS_VolumeTool.hxx"
31 #include "SMESHDS_GroupBase.hxx"
32 #include "SMESHDS_Mesh.hxx"
33 #include "SMESH_Comment.hxx"
34
35 #include <limits>
36 #include <cgnslib.h>
37
38 #if CGNS_VERSION < 3100
39 # define cgsize_t int
40 #endif
41
42 using namespace std;
43
44 namespace
45 {
46   //================================================================================
47   /*!
48    * \brief Return interlace and type of CGNS element for the given SMDSAbs_EntityType
49    */
50   //================================================================================
51
52   const int* getInterlaceAndType( const SMDSAbs_EntityType      smType,
53                                   CGNS_ENUMT( ElementType_t ) & cgType )
54   {
55     static vector< const int* >                 interlaces;
56     static vector< CGNS_ENUMT( ElementType_t )> cgTypes; 
57     if ( interlaces.empty() )
58     {
59       interlaces.resize( SMDSEntity_Last, 0 );
60       cgTypes.resize( SMDSEntity_Last, CGNS_ENUMV( ElementTypeNull ));
61       {
62         static int ids[] = {0};
63         interlaces[SMDSEntity_0D] = ids;
64         cgTypes   [SMDSEntity_0D] = CGNS_ENUMV( NODE );
65       }
66       {
67         static int ids[] = { 0, 1 };
68         interlaces[SMDSEntity_Edge] = ids;
69         cgTypes   [SMDSEntity_Edge] = CGNS_ENUMV( BAR_2 );
70       }
71       {
72         static int ids[] = { 0, 1, 2 };
73         interlaces[SMDSEntity_Quad_Edge] = ids;
74         cgTypes   [SMDSEntity_Quad_Edge] = CGNS_ENUMV( BAR_3 );
75       }
76       {
77         static int ids[] = { 0, 2, 1 };
78         interlaces[SMDSEntity_Triangle] = ids;
79         cgTypes   [SMDSEntity_Triangle] = CGNS_ENUMV( TRI_3 );
80       }
81       {
82         static int ids[] = { 0, 2, 1, 5, 4, 3 };
83         interlaces[SMDSEntity_Quad_Triangle] = ids;
84         cgTypes   [SMDSEntity_Quad_Triangle] = CGNS_ENUMV( TRI_6 );
85         interlaces[SMDSEntity_BiQuad_Triangle] = ids;
86         cgTypes   [SMDSEntity_BiQuad_Triangle] = CGNS_ENUMV( TRI_6 );
87       }
88       {
89         static int ids[] = { 0, 3, 2, 1 };
90         interlaces[SMDSEntity_Quadrangle] = ids;
91         cgTypes   [SMDSEntity_Quadrangle] = CGNS_ENUMV( QUAD_4 );
92       }
93       {
94         static int ids[] = { 0,3,2,1,7,6,5,4 };
95         interlaces[SMDSEntity_Quad_Quadrangle] = ids;
96         cgTypes   [SMDSEntity_Quad_Quadrangle] = CGNS_ENUMV( QUAD_8 );
97       }
98       {
99         static int ids[] = { 0,3,2,1,7,6,5,4,8 };
100         interlaces[SMDSEntity_BiQuad_Quadrangle] = ids;
101         cgTypes   [SMDSEntity_BiQuad_Quadrangle] = CGNS_ENUMV( QUAD_9 );
102       }
103       {
104         static int ids[] = { 0, 2, 1, 3 };
105         interlaces[SMDSEntity_Tetra] = ids;
106         cgTypes   [SMDSEntity_Tetra] = CGNS_ENUMV( TETRA_4 );
107       }
108       {
109         static int ids[] = { 0,2,1,3,6,5,4,7,9,8 };
110         interlaces[SMDSEntity_Quad_Tetra] = ids;
111         cgTypes   [SMDSEntity_Quad_Tetra] = CGNS_ENUMV( TETRA_10 );
112       }
113       {
114         static int ids[] = { 0,3,2,1,4 };
115         interlaces[SMDSEntity_Pyramid] = ids;
116         cgTypes   [SMDSEntity_Pyramid] = CGNS_ENUMV( PYRA_5 );
117       }
118       {
119         static int ids[] = { 0,3,2,1,4,8,7,6,5,9,12,11,10 };
120         interlaces[SMDSEntity_Quad_Pyramid] = ids;
121         cgTypes   [SMDSEntity_Quad_Pyramid] = CGNS_ENUMV( PYRA_13 );
122       }
123       {
124         static int ids[] = { 0,2,1,3,5,4 };
125         interlaces[SMDSEntity_Penta] = ids;
126         cgTypes   [SMDSEntity_Penta] = CGNS_ENUMV( PENTA_6 );
127       }
128       {
129         static int ids[] = { 0,2,1,3,5,4,8,7,6,9,11,10,14,13,12 };
130         interlaces[SMDSEntity_Quad_Penta] = ids;
131         cgTypes   [SMDSEntity_Quad_Penta] = CGNS_ENUMV( PENTA_15 );
132       }
133       {
134         static int ids[] = { 0,2,1,3,5,4,8,7,6,9,11,10,14,13,12,15,16,17 }; // TODO: check CGNS ORDER
135         interlaces[SMDSEntity_BiQuad_Penta] = ids;
136         cgTypes   [SMDSEntity_BiQuad_Penta] = CGNS_ENUMV( PENTA_18 );
137       }
138       {
139         static int ids[] = { 0,3,2,1,4,7,6,5 };
140         interlaces[SMDSEntity_Hexa] = ids;
141         cgTypes   [SMDSEntity_Hexa] = CGNS_ENUMV( HEXA_8 );
142       }
143       {
144         static int ids[] = { 0,3,2,1,4,7,6,5,11,10,9,8,12,15,14,13,19,18,17,16 };
145         interlaces[SMDSEntity_Quad_Hexa] = ids;
146         cgTypes   [SMDSEntity_Quad_Hexa] = CGNS_ENUMV( HEXA_20 );
147       }
148       {
149         static int ids[] = { 0,3,2,1,4,7,6,5,11,10,9,8,12,15,14,13,19,18,17,16,
150                              20, 24,23,22,21, 25, 26};
151         interlaces[SMDSEntity_TriQuad_Hexa] = ids;
152         cgTypes   [SMDSEntity_TriQuad_Hexa] = CGNS_ENUMV( HEXA_27 );
153       }
154       {
155         cgTypes[SMDSEntity_Polygon]         = CGNS_ENUMV( NGON_n );
156         cgTypes[SMDSEntity_Quad_Polygon]    = CGNS_ENUMV( NGON_n );
157         cgTypes[SMDSEntity_Polyhedra]       = CGNS_ENUMV( NFACE_n );
158         cgTypes[SMDSEntity_Hexagonal_Prism] = CGNS_ENUMV( NFACE_n );
159       }
160     }
161     cgType  = cgTypes[ smType ];
162     return interlaces[ smType ];
163   }
164
165   //================================================================================
166   /*!
167    * \brief Cut off type of boundary condition from the group name
168    */
169   //================================================================================
170
171   CGNS_ENUMT( BCType_t ) getBCType( string& groupName )
172   {
173     CGNS_ENUMT( BCType_t ) bcType = CGNS_ENUMV( BCGeneral ); // default type
174
175     // boundary condition type starts from "BC"
176     size_t bcBeg = groupName.find("BC");
177     if ( bcBeg != string::npos )
178     {
179       for ( int t = 0; t < NofValidBCTypes; ++t )
180       {
181         CGNS_ENUMT( BCType_t ) type = CGNS_ENUMT( BCType_t )( t );
182         string typeName = cg_BCTypeName( type );
183         if ( typeName == &groupName[0] + bcBeg )
184         {
185           bcType = type;
186           while ( bcBeg > 0 && isspace( bcBeg-1 ))
187             --bcBeg;
188           if ( bcBeg == 0 )
189             groupName = "Group";
190           else
191             groupName = groupName.substr( 0, bcBeg-1 );
192         }
193       }
194     }
195     return bcType;
196   }
197
198   //================================================================================
199   /*!
200    * \brief Sortable face of a polyhedron
201    */
202   struct TPolyhedFace
203   {
204     smIdType _id; // id of NGON_n
205     vector< smIdType > _nodes; // lowest node IDs used for sorting
206
207     TPolyhedFace( const SMDS_MeshNode** nodes, const smIdType nbNodes, smIdType ID):_id(ID)
208     {
209       set< smIdType > ids;
210       for ( smIdType i = 0; i < nbNodes; ++i )
211         ids.insert( nodes[i]->GetID() );
212
213       _nodes.resize( 3 ); // std::min( nbNodes, 4 )); hope 3 nodes is enough
214       set< smIdType >::iterator idIt = ids.begin();
215       for ( size_t j = 0; j < _nodes.size(); ++j, ++idIt )
216         _nodes[j] = *idIt;
217     }
218     bool operator< (const TPolyhedFace& o ) const
219     {
220       return _nodes < o._nodes;
221     }
222   };
223   //================================================================================
224   /*!
225    * \brief Return CGNS id of an element
226    */
227   //================================================================================
228
229   cgsize_t cgnsID( const SMDS_MeshElement*                         elem,
230                    const map< const SMDS_MeshElement*, cgsize_t >& elem2cgID )
231   {
232     map< const SMDS_MeshElement*, cgsize_t >::const_iterator e2id = elem2cgID.find( elem );
233     return ( e2id == elem2cgID.end() ? elem->GetID() : e2id->second );
234   }
235
236   //================================================================================
237   /*!
238    * \brief save nb nodes of a polygon
239    */
240   //================================================================================
241
242   void addPolySize( const cgsize_t            nbNodes,
243                     std::vector< cgsize_t >& elemData,
244                     std::vector< cgsize_t >& polyOffset )
245   {
246 #if CGNS_VERSION < 4100
247     elemData.push_back( nbNodes );
248     polyOffset.clear(); // just avoid warning: unused parameter
249 #else
250     polyOffset.push_back((cgsize_t) elemData.size() );
251     (void)nbNodes; // avoid warning: unused parameter
252 #endif
253   }
254
255 } // namespace
256
257 //================================================================================
258 /*!
259  * \brief Write the mesh into the CGNS file
260  */
261 //================================================================================
262
263 Driver_Mesh::Status DriverCGNS_Write::Perform()
264 {
265   myErrorMessages.clear();
266
267   if ( !myMesh || myMesh->GetMeshInfo().NbElements() < 1 )
268     return addMessage( !myMesh ? "NULL mesh" : "Empty mesh (no elements)", /*fatal = */true );
269
270   if ( Driver_Mesh::IsMeshTooLarge< cgsize_t >( myMesh, /*checkIDs =*/ false))
271     return DRS_TOO_LARGE_MESH;
272
273   // open the file
274   if ( cg_open(myFile.c_str(), CG_MODE_MODIFY, &_fn) != CG_OK &&
275        cg_open(myFile.c_str(), CG_MODE_WRITE,  &_fn) != CG_OK )
276     return addMessage( cg_get_error(), /*fatal = */true );
277
278   // create a Base
279   // --------------
280
281   const int spaceDim = 3;
282   int        meshDim = 1;
283   if ( myMesh->NbFaces()   > 0 ) meshDim = 2;
284   if ( myMesh->NbVolumes() > 0 ) meshDim = 3;
285
286   if ( myMeshName.empty() )
287   {
288     int nbases = 0;
289     if ( cg_nbases( _fn, &nbases) == CG_OK )
290       myMeshName = ( SMESH_Comment("Base_") << nbases+1 );
291     else
292       myMeshName = "Base_0";
293   }
294   int iBase;
295   if ( cg_base_write( _fn, myMeshName.c_str(), meshDim, spaceDim, &iBase ))
296     return addMessage( cg_get_error(), /*fatal = */true );
297
298   // create a Zone
299   // --------------
300
301   smIdType nbCells = myMesh->NbEdges();
302   if ( meshDim == 3 )
303     nbCells = myMesh->NbVolumes();
304   else if ( meshDim == 2 )
305     nbCells = myMesh->NbFaces();
306
307   cgsize_t size[9] = { FromSmIdType<cgsize_t>( myMesh->NbNodes() ),
308                        FromSmIdType<cgsize_t>( nbCells ),
309                        /*NBoundVertex=*/0, 0,0,0,0,0,0 };
310   int iZone;
311   if ( cg_zone_write( _fn, iBase, "SMESH_Mesh", size,
312                       CGNS_ENUMV( Unstructured ), &iZone) != CG_OK )
313     return addMessage( cg_get_error(), /*fatal = */true );
314
315   // Map to store only elements whose an SMDS ID differs from a CGNS one
316   typedef map< const SMDS_MeshElement*, cgsize_t > TElem2cgIDMap;
317   vector< TElem2cgIDMap > elem2cgIDByEntity( SMDSEntity_Last );
318   TElem2cgIDMap::iterator elem2cgIDIter;
319
320   TElem2cgIDMap & n2cgID = elem2cgIDByEntity[ SMDSEntity_Node ];
321
322   // Write nodes
323   // ------------
324   {
325     vector< double > coords( myMesh->NbNodes() );
326     int iC;
327     // X
328     SMDS_NodeIteratorPtr nIt = myMesh->nodesIterator();
329     for ( int i = 0; nIt->more(); ++i ) coords[i] = nIt->next()->X();
330     if ( cg_coord_write( _fn, iBase, iZone, CGNS_ENUMV(RealDouble),
331                           "CoordinateX", &coords[0], &iC) != CG_OK )
332       return addMessage( cg_get_error(), /*fatal = */true );
333     // Y
334     nIt = myMesh->nodesIterator();
335     for ( int i = 0; nIt->more(); ++i ) coords[i] = nIt->next()->Y();
336     if ( cg_coord_write( _fn, iBase, iZone, CGNS_ENUMV(RealDouble),
337                           "CoordinateY", &coords[0], &iC) != CG_OK )
338       return addMessage( cg_get_error(), /*fatal = */true );
339     // Z
340     nIt = myMesh->nodesIterator();
341     for ( int i = 0; nIt->more(); ++i ) coords[i] = nIt->next()->Z();
342     if ( cg_coord_write( _fn, iBase, iZone, CGNS_ENUMV(RealDouble),
343                           "CoordinateZ", &coords[0], &iC) != CG_OK )
344       return addMessage( cg_get_error(), /*fatal = */true );
345
346     // store CGNS ids of nodes
347     nIt = myMesh->nodesIterator();
348     for ( int i = 0; nIt->more(); ++i )
349     {
350       const SMDS_MeshElement* n = nIt->next();
351       if ( n->GetID() != i+1 )
352         n2cgID.insert( n2cgID.end(), make_pair( n, i+1 ));
353     }
354   }
355   // Write elements
356   // ---------------
357   
358   cgsize_t cgID = 1, startID;
359
360   // write into a section all successive elements of one geom type
361   int iSec;
362   vector< cgsize_t > elemData, polyOffset;
363   SMDS_ElemIteratorPtr elemIt = myMesh->elementsIterator();
364   vector< SMDS_ElemIteratorPtr > elemItVec;
365   if ( _elementsByType )
366   {
367     // create an iterator returning all elements by type
368     for ( int type = SMDSEntity_Node + 1; type < SMDSEntity_Last; ++type )
369     {
370       if ( type == SMDSEntity_Ball )
371         continue; // not supported
372       elemIt = myMesh->elementEntityIterator( SMDSAbs_EntityType( type ));
373       if ( elemIt->more() )
374         elemItVec.push_back( elemIt );
375     }
376     typedef SMDS_IteratorOnIterators< const SMDS_MeshElement*,
377                                       vector< SMDS_ElemIteratorPtr > > TVecIterator;
378     elemIt.reset( new TVecIterator( elemItVec ));
379   }
380
381   const SMDS_MeshElement* elem = elemIt->next();
382   while ( elem )
383   {
384     const SMDSAbs_EntityType elemType = elem->GetEntityType();
385     CGNS_ENUMT( ElementType_t ) cgType;
386     const int* interlace = getInterlaceAndType( elemType, cgType );
387
388     TElem2cgIDMap & elem2cgID = elem2cgIDByEntity[ elemType ];
389
390     elemData.clear();
391     startID = cgID;
392
393     if ( interlace ) // STANDARD elements
394     {
395       int cgnsNbNodes; // get nb nodes by element type, that can be less that elem->NbNodes()
396       cg_npe( cgType, &cgnsNbNodes );
397       do
398       {
399         for ( int i = 0; i < cgnsNbNodes; ++i )
400           elemData.push_back( cgnsID( elem->GetNode( interlace[i] ), n2cgID ));
401         if ( elem->GetID() != cgID )
402           elem2cgID.insert( elem2cgID.end(), make_pair( elem, cgID ));
403         ++cgID;
404         elem = elemIt->more() ? elemIt->next() : 0;
405       }
406       while ( elem && elem->GetEntityType() == elemType );
407     }
408     else if ( elemType == SMDSEntity_Polygon ) // POLYGONS
409       do
410       {
411         addPolySize( elem->NbNodes(), elemData, polyOffset );
412         for ( int i = 0, nb = elem->NbNodes(); i < nb; ++i )
413           elemData.push_back( cgnsID( elem->GetNode(i), n2cgID ));
414         if ( elem->GetID() != cgID )
415           elem2cgID.insert( elem2cgID.end(), make_pair( elem, cgID ));
416         ++cgID;
417         elem = elemIt->more() ? elemIt->next() : 0;
418       }
419       while ( elem && elem->GetEntityType() == elemType );
420
421     else if ( elemType == SMDSEntity_Quad_Polygon ) // QUADRATIC POLYGONS
422       do // write as linear NGON_n
423       {
424         addPolySize( elem->NbNodes(), elemData, polyOffset );
425         interlace = & SMDS_MeshCell::interlacedSmdsOrder( SMDSEntity_Quad_Polygon,
426                                                           elem->NbNodes() )[0];
427         for ( int i = 0, nb = elem->NbNodes(); i < nb; ++i )
428           elemData.push_back( cgnsID( elem->GetNode( interlace[i] ), n2cgID ));
429         if ( elem->GetID() != cgID )
430           elem2cgID.insert( elem2cgID.end(), make_pair( elem, cgID ));
431         ++cgID;
432         elem = elemIt->more() ? elemIt->next() : 0;
433       }
434       while ( elem && elem->GetEntityType() == elemType );
435
436     else if ( elemType == SMDSEntity_Polyhedra ||
437               elemType == SMDSEntity_Hexagonal_Prism) // POLYHEDRA
438     {
439       // to save polyhedrons after all
440       const SMDS_MeshInfo& meshInfo = myMesh->GetMeshInfo();
441       if ( meshInfo.NbPolyhedrons() == meshInfo.NbElements() - cgID + 1 )
442         break; // only polyhedrons remain
443       while ( elem && elem->GetEntityType() == elemType )
444         elem = elemIt->more() ? elemIt->next() : 0;
445       continue;
446     }
447     else // skip NOT SUPPORTED elements
448     {
449       while ( elemIt->more() )
450       {
451         elem = elemIt->next();
452         if ( elem->GetEntityType() != elemType )
453           break;
454       }
455     }
456
457     SMESH_Comment sectionName( cg_ElementTypeName( cgType ));
458     sectionName << " " << startID << " - " << cgID-1;
459
460
461 #if CGNS_VERSION >= 4000
462     if ( !polyOffset.empty() )
463     {
464       polyOffset.push_back((cgsize_t) elemData.size() );
465       if ( cg_poly_section_write(_fn, iBase, iZone, sectionName.c_str(), cgType, startID,
466                                  cgID-1, /*nbndry=*/0,
467                                  elemData.data(), polyOffset.data(), &iSec) != CG_OK )
468         return addMessage( cg_get_error(), /*fatal = */true );
469     }
470     else
471 #endif
472     {
473       if ( cg_section_write(_fn, iBase, iZone, sectionName.c_str(), cgType, startID,
474                             cgID-1, /*nbndry=*/0, elemData.data(), &iSec) != CG_OK )
475         return addMessage( cg_get_error(), /*fatal = */true );
476     }
477   } // while ( elem )
478
479   // Write polyhedral volumes
480   // -------------------------
481
482   if ( myMesh->GetMeshInfo().NbElements() > cgID-1 ) // polyhedra or hexagonal prisms remain
483   {
484     // the polyhedron (NFACE_n) is described as a set of signed face IDs,
485     // so first we are to write all polygones (NGON_n) bounding polyhedrons
486
487     vector< cgsize_t > faceData;
488     set< TPolyhedFace > faces;
489     set< TPolyhedFace >::iterator faceInSet;
490     vector<const SMDS_MeshNode *> faceNodesVec;
491     vector< cgsize_t > faceOffset;
492     int nbPolygones = 0, faceID;
493
494     SMDS_VolumeTool vol;
495
496     elemData.clear();
497     polyOffset.clear();
498
499     int nbPolyhTreated = 0;
500
501     TElem2cgIDMap * elem2cgID = 0;
502     TElem2cgIDMap & n2cgID    = elem2cgIDByEntity[ SMDSEntity_Node ];
503
504     SMDS_ElemIteratorPtr elemIt = myMesh->elementsIterator();
505     while ( elemIt->more() )
506     {
507       elem = elemIt->next();
508       SMDSAbs_EntityType type = elem->GetEntityType();
509       if ( type == SMDSEntity_Polyhedra ||
510            type == SMDSEntity_Hexagonal_Prism )
511       {
512         ++nbPolyhTreated;
513         vol.Set( elem );
514         vol.SetExternalNormal();
515         const int nbFaces = vol.NbFaces();
516         addPolySize( nbFaces, elemData, polyOffset );
517         for ( int iF = 0; iF < nbFaces; ++iF )
518         {
519           const int nbNodes = vol.NbFaceNodes( iF );
520           const SMDS_MeshNode** faceNodes = vol.GetFaceNodes( iF );
521           faceNodesVec.assign( faceNodes, faceNodes + nbNodes );
522           if (( elem = myMesh->FindElement( faceNodesVec, SMDSAbs_Face, /*noMedium=*/false)))
523           {
524             // a face of the polyhedron is present in the mesh
525             faceID = cgnsID( elem, elem2cgIDByEntity[ elem->GetEntityType() ]);
526           }
527           else if ( vol.IsFreeFace( iF ))
528           {
529             // the face is not shared by volumes
530             faceID = cgID++;
531             ++nbPolygones;
532             addPolySize( nbNodes, faceData, faceOffset );
533             for ( int i = 0; i < nbNodes; ++i )
534               faceData.push_back( cgnsID( faceNodes[i], n2cgID ));
535           }
536           else
537           {
538             TPolyhedFace face( faceNodes, nbNodes, cgID );
539             faceInSet = faces.insert( faces.end(), face );
540             if ( faceInSet->_id == cgID ) // the face encounters for the 1st time
541             {
542               faceID = cgID++;
543               ++nbPolygones;
544               addPolySize( nbNodes, faceData, faceOffset );
545               for ( int i = 0; i < nbNodes; ++i )
546                 faceData.push_back( cgnsID( faceNodes[i], n2cgID ));
547             }
548             else
549             {
550               // the face encounters for the 2nd time; we hope it won't encounter once more,
551               // for that we can erase it from the set of faces
552               faceID = -faceInSet->_id;
553               faces.erase( faceInSet );
554             }
555           }
556           elemData.push_back( faceID );
557         }
558       }
559     }
560
561 #if CGNS_VERSION >= 4000
562     if ( nbPolygones > 0 )
563     {
564       faceOffset.push_back((cgsize_t) faceData.size() );
565       if ( cg_poly_section_write(_fn, iBase, iZone, "Faces of Polyhedrons", CGNS_ENUMV( NGON_n ),
566                                  cgID - nbPolygones, cgID-1, /*nbndry=*/0,
567                                  faceData.data(), faceOffset.data(), &iSec) != CG_OK )
568         return addMessage( cg_get_error(), /*fatal = */true );
569     }
570
571     polyOffset.push_back((cgsize_t) elemData.size() );
572     if ( cg_poly_section_write(_fn, iBase, iZone, "Polyhedrons",
573                                CGNS_ENUMV( NFACE_n ), cgID, cgID+nbPolyhTreated-1,
574                                /*nbndry=*/0, &elemData[0], polyOffset.data(), &iSec) != CG_OK )
575       return addMessage( cg_get_error(), /*fatal = */true );
576 #else
577     if ( nbPolygones > 0 )
578     {
579       if ( cg_section_write(_fn, iBase, iZone, "Faces of Polyhedrons",
580                             CGNS_ENUMV( NGON_n ), cgID - nbPolygones, cgID-1,
581                             /*nbndry=*/0, &faceData[0], &iSec) != CG_OK )
582         return addMessage( cg_get_error(), /*fatal = */true );
583     }
584
585     if ( cg_section_write(_fn, iBase, iZone, "Polyhedrons",
586                           CGNS_ENUMV( NFACE_n ), cgID, cgID+nbPolyhTreated-1,
587                           /*nbndry=*/0, &elemData[0], &iSec) != CG_OK )
588       return addMessage( cg_get_error(), /*fatal = */true );
589 #endif
590
591     if ( !myMesh->GetGroups().empty() )
592     {
593       // store CGNS ids of polyhedrons
594       elem2cgID = &elem2cgIDByEntity[ SMDSEntity_Polyhedra ];
595       elemIt = myMesh->elementsIterator();
596       while ( elemIt->more() )
597       {
598         elem = elemIt->next();
599         if ( elem->GetEntityType() == SMDSEntity_Polyhedra )
600         {
601           if ( elem->GetID() != cgID )
602             elem2cgID->insert( elem2cgID->end(), make_pair( elem, cgID ));
603           ++cgID;
604         }
605       }
606     }
607   } // write polyhedral volumes
608
609
610   // Write groups as boundary conditions
611   // ------------------------------------
612
613   const set<SMESHDS_GroupBase*>& groups = myMesh->GetGroups();
614   set<SMESHDS_GroupBase*>::const_iterator grpIt = groups.begin();
615   set< string > groupNames; groupNames.insert(""); // to avoid duplicated and empty names
616   for ( ; grpIt != groups.end(); ++grpIt )
617   {
618     const SMESHDS_GroupBase* group = *grpIt;
619
620     // write BC location (default is Vertex)
621     CGNS_ENUMT( GridLocation_t ) location = CGNS_ENUMV( Vertex );
622     if ( group->GetType() != SMDSAbs_Node )
623     {
624       switch ( meshDim ) {
625       case 3:
626         switch ( group->GetType() ) {
627 #if CGNS_VERSION > 3130
628         case SMDSAbs_Volume: location = CGNS_ENUMV( CellCenter ); break;
629 #else
630         case SMDSAbs_Volume: location = CGNS_ENUMV( FaceCenter ); break;
631 #endif
632         case SMDSAbs_Face:   location = CGNS_ENUMV( FaceCenter ); break;
633         case SMDSAbs_Edge:   location = CGNS_ENUMV( EdgeCenter ); break;
634         default:;
635         }
636         break;
637       case 2:
638         switch ( group->GetType() ) {
639 #if CGNS_VERSION > 3130
640         case SMDSAbs_Face: location = CGNS_ENUMV( CellCenter ); break;
641 #else
642         case SMDSAbs_Face: location = CGNS_ENUMV( FaceCenter ); break;
643 #endif
644         case SMDSAbs_Edge: location = CGNS_ENUMV( EdgeCenter ); break;
645         default:;
646         }
647         break;
648       case 1:
649 #if CGNS_VERSION > 3130
650         location = CGNS_ENUMV( CellCenter ); break;
651 #else
652         location = CGNS_ENUMV( EdgeCenter ); break;
653 #endif
654         break;
655       }
656     }
657
658     // try to extract type of boundary condition from the group name
659     string name = group->GetStoreName();
660     CGNS_ENUMT( BCType_t ) bcType = getBCType( name );
661     if ( !groupNames.insert( name ).second ) // assure name uniqueness
662     {
663       int index = 1;
664       string newName;
665       do {
666         newName = SMESH_Comment( name ) << "_" << index++;
667       }
668       while ( !groupNames.insert( newName ).second );
669       name = newName;
670     }
671     // write IDs of elements
672     vector< cgsize_t > pnts;
673     pnts.reserve( group->Extent() );
674     SMDS_ElemIteratorPtr elemIt = group->GetElements();
675     while ( elemIt->more() )
676     {
677       const SMDS_MeshElement* elem = elemIt->next();
678       pnts.push_back( cgnsID( elem, elem2cgIDByEntity[ elem->GetEntityType() ]));
679     }
680     if ( pnts.size() == 0 )
681       continue; // can't store empty group
682     int iBC;
683     if ( cg_boco_write( _fn, iBase, iZone, name.c_str(), bcType,
684                         CGNS_ENUMV( PointList ), pnts.size(), &pnts[0], &iBC) != CG_OK )
685       return addMessage( cg_get_error(), /*fatal = */true);
686
687     // write BC location
688     if ( location != CGNS_ENUMV( Vertex ) || meshDim == 1 )
689     {
690       if ( cg_boco_gridlocation_write( _fn, iBase, iZone, iBC, location) != CG_OK )
691         return addMessage( cg_get_error(), /*fatal = */false);
692     }
693   }
694   return DRS_OK;
695 }
696
697 //================================================================================
698 /*!
699  * \brief Constructor
700  */
701 //================================================================================
702
703 DriverCGNS_Write::DriverCGNS_Write(): _fn(0), _elementsByType( false )
704 {
705 }
706
707 //================================================================================
708 /*!
709  * \brief Close the cgns file at destruction
710  */
711 //================================================================================
712
713 DriverCGNS_Write::~DriverCGNS_Write()
714 {
715   if ( _fn > 0 )
716     cg_close( _fn );
717 }