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   }
 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 }