src/DriverCGNS/DriverCGNS_Read.cxx

   1 // Copyright (C) 2007-2011  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.
  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_Read.cxx
  23 // Created   : Thu Jun 30 10:33:31 2011
  24 // Author    : Edward AGAPOV (eap)
  25
  26 #include "DriverCGNS_Read.hxx"
  27
  28 #include "SMDS_MeshNode.hxx"
  29 #include "SMESHDS_Group.hxx"
  30 #include "SMESHDS_Mesh.hxx"
  31 #include "SMESH_Comment.hxx"
  32
  33 #include <gp_XYZ.hxx>
  34
  35 #include <cgnslib.h>
  36
  37 #include <map>
  38
  39 #if CGNS_VERSION < 3100
  40 # define cgsize_t int
  41 #endif
  42
  43 #define NB_ZONE_SIZE_VAL 9
  44 #define CGNS_NAME_SIZE 33
  45 #define CGNS_STRUCT_RANGE_SZ 6
  46
  47 using namespace std;
  48
  49 namespace
  50 {
  51   //================================================================================
  52   /*!
  53    * \brief Data of a zone
  54    */
  55   struct TZoneData
  56   {
  57     int                    _id;
  58     int                    _nodeIdShift; // nb nodes in previously read zones
  59     int                    _elemIdShift; // nb faces in previously read zones
  60     int                    _nbNodes, _nbElems;
  61     int                    _meshDim;
  62     int                    _sizeX, _sizeY, _sizeZ, _nbCells; // structured
  63     cgsize_t               _sizes[NB_ZONE_SIZE_VAL];
  64     CGNS_ENUMT(ZoneType_t) _type;
  65     map< int, int >        _nodeReplacementMap;/* key:   id of node to replace (in this zone),
  66                                                   value: id of node to replace by (in another zone)
  67                                                   id values include _nodeIdShift of the zones */
  68     void SetSizeAndDim( cgsize_t* sizes, int meshDim )
  69     {
  70       _meshDim = meshDim;
  71       memcpy( _sizes, sizes, NB_ZONE_SIZE_VAL*sizeof(cgsize_t));
  72       _sizeX = _sizes[0];
  73       _sizeY = _meshDim > 1 ? _sizes[1] : 0;
  74       _sizeZ = _meshDim > 2 ? _sizes[2] : 0;
  75       _nbCells = (_sizeX - 1) * ( _meshDim > 1 ? _sizeY : 1 ) * ( _meshDim > 2 ? _sizeZ : 1 );
  76     }
  77     bool IsStructured() const { return ( _type == CGNS_ENUMV( Structured )); }
  78     int IndexSize() const { return IsStructured() ? _meshDim : 1; }
  79     string ReadZonesConnection(int file, int base, const map< string, TZoneData >& zonesByName);
  80     void ReplaceNodes( cgsize_t* ids, int nbIds, int idShift = 0 ) const;
  81
  82     // Methods for a structured zone
  83
  84     int NodeID( int i, int j, int k = 1 ) const
  85     {
  86       return _nodeIdShift + (k-1)*_sizeX*_sizeY + (j-1)*_sizeX + i;
  87     }
  88     int NodeID( const gp_XYZ& ijk ) const
  89     {
  90       return NodeID( int(ijk.X()), int(ijk.Y()), int(ijk.Z()));
  91     }
  92     void CellNodes( int i, int j, int k, cgsize_t* ids ) const
  93     {
  94       ids[0] = NodeID( i  , j  , k  );
  95       ids[1] = NodeID( i  , j+1, k  );
  96       ids[2] = NodeID( i+1, j+1, k  );
  97       ids[3] = NodeID( i+1, j  , k  );
  98       ids[4] = NodeID( i  , j  , k+1);
  99       ids[5] = NodeID( i  , j+1, k+1);
 100       ids[6] = NodeID( i+1, j+1, k+1);
 101       ids[7] = NodeID( i+1, j  , k+1);
 102     }
 103     void CellNodes( int i, int j, cgsize_t* ids ) const
 104     {
 105       ids[0] = NodeID( i  , j   );
 106       ids[1] = NodeID( i  , j+1 );
 107       ids[2] = NodeID( i+1, j+1 );
 108       ids[3] = NodeID( i+1, j   );
 109     }
 110     void IFaceNodes( int i, int j, int k, cgsize_t* ids ) const // face perpendiculaire to X (3D)
 111     {
 112       ids[0] = NodeID( i, j, k );
 113       ids[1] = ids[0] + _sizeX*( i==_sizeX ? 1 : _sizeY );
 114       ids[2] = ids[0] + _sizeX*( _sizeY + 1 );
 115       ids[3] = ids[0] + _sizeX*( i==_sizeX ? _sizeY : 1 );
 116     }
 117     void JFaceNodes( int i, int j, int k, cgsize_t* ids ) const
 118     {
 119       ids[0] = NodeID( i, j, k );
 120       ids[1] = ids[0] + ( j==_sizeY ? _sizeX*_sizeY : 1);
 121       ids[2] = ids[0] + _sizeX*_sizeY + 1;
 122       ids[3] = ids[0] + ( j==_sizeY ? 1 : _sizeX*_sizeY);
 123     }
 124     void KFaceNodes( int i, int j, int k, cgsize_t* ids ) const
 125     {
 126       ids[0] = NodeID( i, j, k );
 127       ids[1] = ids[0] + ( k==_sizeZ ? 1 : _sizeX);
 128       ids[2] = ids[0] + _sizeX + 1;
 129       ids[3] = ids[0] + ( k==_sizeZ ? _sizeX : 1);
 130     }
 131     void IEdgeNodes( int i, int j, int k, cgsize_t* ids ) const // edge perpendiculaire to X (2D)
 132     {
 133       ids[0] = NodeID( i, j, 0 );
 134       ids[1] = ids[0] + _sizeX;
 135     }
 136     void JEdgeNodes( int i, int j, int k, cgsize_t* ids ) const
 137     {
 138       ids[0] = NodeID( i, j, 0 );
 139       ids[1] = ids[0] + 1;
 140     }
 141 #define gpXYZ2IJK(METHOD)                                     \
 142     void METHOD( const gp_XYZ& ijk, cgsize_t* ids ) const {        \
 143       METHOD( int(ijk.X()), int(ijk.Y()), int(ijk.Z()), ids); \
 144     }
 145     gpXYZ2IJK( IFaceNodes )
 146     gpXYZ2IJK( JFaceNodes )
 147     gpXYZ2IJK( KFaceNodes )
 148     gpXYZ2IJK( IEdgeNodes )
 149     gpXYZ2IJK( JEdgeNodes )
 150   };
 151
 152   //================================================================================
 153   /*!
 154    * \brief Iterator over nodes of the structired grid using FORTRAN multidimensional
 155    * array ordering.
 156    */
 157   class TPointRangeIterator
 158   {
 159     int _beg[3], _end[3], _cur[3], _dir[3], _dim;
 160     bool _more;
 161   public:
 162     TPointRangeIterator( const cgsize_t* range, int dim ):_dim(dim)
 163     {
 164       _more = false;
 165       for ( int i = 0; i < dim; ++i )
 166       {
 167         _beg[i] = range[i];
 168         _end[i] = range[i+dim];
 169         _dir[i] = _end[i] < _beg[i] ? -1 : 1;
 170         _end[i] += _dir[i];
 171         _cur[i] = _beg[i];
 172         if ( _end[i] - _beg[i] )
 173           _more = true;
 174       }
 175 //       for ( int i = dim; i < 3; ++i )
 176 //         _cur[i] = _beg[i] = _end[i] = _dir[i] = 0;
 177     }
 178     bool More() const
 179     {
 180       return _more;
 181     }
 182     gp_XYZ Next()
 183     {
 184       gp_XYZ res( _cur[0], _cur[1], _cur[2] );
 185       for ( int i = 0; i < _dim; ++i )
 186       {
 187         _cur[i] += _dir[i];
 188         if ( _cur[i]*_dir[i] < _end[i]*_dir[i] )
 189           break;
 190         if ( i+1 < _dim )
 191           _cur[i] = _beg[i];
 192         else
 193           _more = false;
 194       }
 195       return res;
 196     }
 197     size_t Size()  const
 198     {
 199       size_t size = 1;
 200       for ( int i = 0; i < _dim; ++i )
 201         size *= _dir[i]*(_end[i]-_beg[i]);
 202       return size;
 203     }
 204     gp_XYZ Begin() const { return gp_XYZ( _beg[0], _beg[1], _beg[2] ); }
 205     //gp_XYZ End() const { return gp_XYZ( _end[0]-1, _end[1]-1, _end[2]-1 ); }
 206   };
 207
 208   //================================================================================
 209   /*!
 210    * \brief Reads zone interface connectivity
 211    *  \param file - file to read
 212    *  \param base - base to read
 213    *  \param zone - zone to replace nodes in
 214    *  \param zonesByName - TZoneData by name
 215    *  \retval string - warning message
 216    *
 217    * see // http://www.grc.nasa.gov/WWW/cgns/CGNS_docs_current/sids/cnct.html
 218    */
 219   //================================================================================
 220
 221   string TZoneData::ReadZonesConnection( int                             file,
 222                                          int                             base,
 223                                          const map< string, TZoneData >& zonesByName)
 224   {
 225     string error;
 226
 227     char connectName[ CGNS_NAME_SIZE ], donorName [ CGNS_NAME_SIZE ];
 228
 229     // ----------------------------
 230     // read zone 1 to 1 interfaces
 231     // ----------------------------
 232     if ( IsStructured() )
 233     {
 234       int nb1to1 = 0;
 235       if ( cg_n1to1 ( file, base, _id, &nb1to1) == CG_OK )
 236       {
 237         cgsize_t range[CGNS_STRUCT_RANGE_SZ], donorRange[CGNS_STRUCT_RANGE_SZ];
 238         int transform[3] = {0,0,0};
 239
 240         for ( int I = 1; I <= nb1to1; ++I )
 241         {
 242           if ( cg_1to1_read(file, base, _id, I, connectName,
 243                             donorName, range, donorRange, transform) == CG_OK )
 244           {
 245             map< string, TZoneData>::const_iterator n_z = zonesByName.find( donorName );
 246             if ( n_z == zonesByName.end() )
 247               continue; // donor zone not yet read
 248             const TZoneData& zone2 = n_z->second;
 249
 250             // set up matrix to transform ijk of the zone to ijk of the zone2
 251             gp_Mat T;
 252             for ( int i = 0; i < _meshDim; ++i )
 253               if ( transform[i] )
 254               {
 255                 int row = Abs(transform[i]);
 256                 int col = i+1;
 257                 int val = transform[i] > 0 ? +1 : -1;
 258                 T( row, col ) = val;
 259               }
 260
 261             // fill nodeReplacementMap
 262             TPointRangeIterator rangeIt1( range, _meshDim );
 263             TPointRangeIterator rangeIt2( donorRange, _meshDim );
 264             gp_XYZ begin1 = rangeIt1.Begin(), begin2 = rangeIt2.Begin(), index1, index2;
 265             if ( &zone2 == this )
 266             {
 267               // not to read twice the same interface with self
 268               TPointRangeIterator rangeIt1bis( range, _meshDim );
 269               if ( rangeIt1bis.More() )
 270               {
 271                 index1 = rangeIt1bis.Next();
 272                 index2 = T * ( index1 - begin1 ) + begin2;
 273                 int node1 = NodeID( index1 );
 274                 int node2 = zone2.NodeID( index2 );
 275                 if ( _nodeReplacementMap.count( node2 ) &&
 276                      _nodeReplacementMap[ node2 ] == node1 )
 277                   continue; // this interface already read
 278               }
 279             }
 280             while ( rangeIt1.More() )
 281             {
 282               index1 = rangeIt1.Next();
 283               index2 = T * ( index1 - begin1 ) + begin2;
 284               int node1 = NodeID( index1 );
 285               int node2 = zone2.NodeID( index2 );
 286               _nodeReplacementMap.insert( make_pair( node1, node2 ));
 287             }
 288           }
 289           else
 290           {
 291             error = cg_get_error();
 292           }
 293         }
 294       }
 295       else
 296       {
 297         error = cg_get_error();
 298       }
 299     }
 300
 301     // ---------------------------------
 302     // read general zone connectivities
 303     // ---------------------------------
 304     int nbConn = 0;
 305     if ( cg_nconns( file, base, _id, &nbConn) == CG_OK )
 306     {
 307       cgsize_t nb, donorNb;
 308       CGNS_ENUMT(GridLocation_t) location;
 309       CGNS_ENUMT(GridConnectivityType_t) connectType;
 310       CGNS_ENUMT(PointSetType_t) ptype, donorPtype;
 311       CGNS_ENUMT(ZoneType_t) donorZonetype;
 312       CGNS_ENUMT(DataType_t) donorDatatype;
 313
 314       for ( int I = 1; I <= nbConn; ++I )
 315       {
 316         if ( cg_conn_info(file, base, _id, I, connectName, &location, &connectType,
 317                           &ptype, &nb, donorName, &donorZonetype, &donorPtype,
 318                           &donorDatatype, &donorNb ) == CG_OK )
 319         {
 320           if ( location != CGNS_ENUMV( Vertex ))
 321             continue; // we do not support cell-to-cell connectivity
 322           if ( ptype != CGNS_ENUMV( PointList ) &&
 323                ptype != CGNS_ENUMV( PointRange ))
 324             continue;
 325           if ( donorPtype != CGNS_ENUMV( PointList ) &&
 326                donorPtype != CGNS_ENUMV( PointRange ))
 327             continue;
 328
 329           map< string, TZoneData>::const_iterator n_z = zonesByName.find( donorName );
 330           if ( n_z == zonesByName.end() )
 331             continue; // donor zone not yet read
 332           const TZoneData& zone2 = n_z->second;
 333
 334           vector< cgsize_t > ids( nb * IndexSize() );
 335           vector< cgsize_t > donorIds( donorNb * zone2.IndexSize() );
 336           if (cg_conn_read ( file, base, _id, I,
 337                              &ids[0], CGNS_ENUMV(Integer), &donorIds[0]) == CG_OK )
 338           {
 339             for ( int isThisZone = 0; isThisZone < 2; ++isThisZone )
 340             {
 341               const TZoneData&            zone = isThisZone ? *this : zone2;
 342               CGNS_ENUMT(PointSetType_t) type = isThisZone ? ptype : donorPtype;
 343               vector< cgsize_t >&      points = isThisZone ? ids : donorIds;
 344               if ( type == CGNS_ENUMV( PointRange ))
 345               {
 346                 TPointRangeIterator rangeIt( &points[0], zone._meshDim );
 347                 points.clear();
 348                 while ( rangeIt.More() )
 349                   points.push_back ( NodeID( rangeIt.Next() ));
 350               }
 351               else if ( zone.IsStructured() )
 352               {
 353                 vector< cgsize_t > resIDs; resIDs.reserve( points.size() / IndexSize() );
 354                 for ( size_t i = 0; i < points.size(); i += IndexSize() )
 355                   resIDs.push_back( zone.NodeID( points[i+0], points[i+1], points[i+2] ));
 356                 resIDs.swap( points );
 357               }
 358               else if ( zone._nodeIdShift > 0 )
 359               {
 360                 for ( size_t i = 0; i < points.size(); ++i )
 361                   points[i] += zone._nodeIdShift;
 362               }
 363             }
 364             for ( size_t i = 0; i < ids.size() && i < donorIds.size(); ++i )
 365               _nodeReplacementMap.insert( make_pair( ids[i], donorIds[i] ));
 366           }
 367           else
 368           {
 369             error = cg_get_error();
 370           }
 371         }
 372         else
 373         {
 374           error = cg_get_error();
 375         }
 376       }
 377     }
 378     else
 379     {
 380       error = cg_get_error();
 381     }
 382     return error;
 383   }
 384
 385   //================================================================================
 386   /*!
 387    * \brief Replaces node ids according to nodeReplacementMap to take into account
 388    *        connection of zones
 389    */
 390   //================================================================================
 391
 392   void TZoneData::ReplaceNodes( cgsize_t* ids, int nbIds, int idShift/* = 0*/ ) const
 393   {
 394     if ( !_nodeReplacementMap.empty() )
 395     {
 396       map< int, int >::const_iterator it, end = _nodeReplacementMap.end();
 397       for ( size_t i = 0; i < nbIds; ++i )
 398         if (( it = _nodeReplacementMap.find( ids[i] + idShift)) != end )
 399           ids[i] = it->second;
 400         else
 401           ids[i] += idShift;
 402     }
 403     else if ( idShift )
 404     {
 405       for ( size_t i = 0; i < nbIds; ++i )
 406         ids[i] += idShift;
 407     }
 408   }
 409   //================================================================================
 410   /*!
 411    * \brief functions adding an element of a particular type
 412    */
 413   SMDS_MeshElement* add_0D(cgsize_t* ids, SMESHDS_Mesh* mesh, int ID)
 414   {
 415     return mesh->Add0DElementWithID( ids[0], ID );
 416   }
 417   SMDS_MeshElement* add_BAR_2(cgsize_t* ids, SMESHDS_Mesh* mesh, int ID)
 418   {
 419     return mesh->AddEdgeWithID( ids[0], ids[1], ID );
 420   }
 421   SMDS_MeshElement* add_BAR_3(cgsize_t* ids, SMESHDS_Mesh* mesh, int ID)
 422   {
 423     return mesh->AddEdgeWithID( ids[0], ids[1], ids[2], ID );
 424   }
 425   SMDS_MeshElement* add_TRI_3(cgsize_t* ids, SMESHDS_Mesh* mesh, int ID)
 426   {
 427     return mesh->AddFaceWithID( ids[0], ids[2], ids[1], ID );
 428   }
 429   SMDS_MeshElement* add_TRI_6(cgsize_t* ids, SMESHDS_Mesh* mesh, int ID)
 430   {
 431     return mesh->AddFaceWithID( ids[0], ids[2], ids[1], ids[5], ids[4], ids[3], ID );
 432   }
 433   SMDS_MeshElement* add_QUAD_4(cgsize_t* ids, SMESHDS_Mesh* mesh, int ID)
 434   {
 435     return mesh->AddFaceWithID( ids[0], ids[3], ids[2], ids[1], ID );
 436   }
 437   SMDS_MeshElement* add_QUAD_8(cgsize_t* ids, SMESHDS_Mesh* mesh, int ID)
 438   {
 439     return mesh->AddFaceWithID( ids[0],ids[3],ids[2],ids[1],ids[7],ids[6],ids[5],ids[4], ID );
 440   }
 441   SMDS_MeshElement* add_TETRA_4(cgsize_t* ids, SMESHDS_Mesh* mesh, int ID)
 442   {
 443     return mesh->AddVolumeWithID( ids[0], ids[2], ids[1], ids[3], ID );
 444   }
 445   SMDS_MeshElement* add_TETRA_10(cgsize_t* ids, SMESHDS_Mesh* mesh, int ID)
 446   {
 447     return mesh->AddVolumeWithID( ids[0],ids[2],ids[1],ids[3],ids[6],
 448                                   ids[5],ids[4],ids[7],ids[9],ids[8], ID );
 449   }
 450   SMDS_MeshElement* add_PYRA_5(cgsize_t* ids, SMESHDS_Mesh* mesh, int ID)
 451   {
 452     return mesh->AddVolumeWithID( ids[0],ids[3],ids[2],ids[1],ids[4],ID );
 453   }
 454   SMDS_MeshElement* add_PYRA_13(cgsize_t* ids, SMESHDS_Mesh* mesh, int ID)
 455   {
 456     return mesh->AddVolumeWithID( ids[0],ids[3],ids[2],ids[1],ids[4],ids[8],ids[7],
 457                                   ids[6],ids[5],ids[9],ids[12],ids[11],ids[10], ID );
 458   }
 459   SMDS_MeshElement* add_PENTA_6(cgsize_t* ids, SMESHDS_Mesh* mesh, int ID)
 460   {
 461     return mesh->AddVolumeWithID( ids[0],ids[2],ids[1],ids[3],ids[5],ids[4], ID );
 462   }
 463   SMDS_MeshElement* add_PENTA_15(cgsize_t* ids, SMESHDS_Mesh* mesh, int ID)
 464   {
 465     return mesh->AddVolumeWithID( ids[0],ids[2],ids[1],ids[3],ids[5],ids[4],ids[8],ids[7],
 466                                   ids[6],ids[9],ids[11],ids[10],ids[14],ids[13],ids[12], ID );
 467   }
 468   SMDS_MeshElement* add_HEXA_8(cgsize_t* ids, SMESHDS_Mesh* mesh, int ID)
 469   {
 470     return mesh->AddVolumeWithID( ids[0],ids[3],ids[2],ids[1],ids[4],ids[7],ids[6],ids[5], ID );
 471   }
 472   SMDS_MeshElement* add_HEXA_20(cgsize_t* ids, SMESHDS_Mesh* mesh, int ID)
 473   {
 474     return mesh->AddVolumeWithID( ids[0],ids[3],ids[2],ids[1],ids[4],ids[7],ids[6],
 475                                   ids[5],ids[11],ids[10],ids[9],ids[8],ids[12],ids[15],
 476                                   ids[14],ids[13],ids[19],ids[18],ids[17],ids[16], ID );
 477   }
 478   SMDS_MeshElement* add_NGON(cgsize_t* ids, SMESHDS_Mesh* mesh, int ID)
 479   {
 480     vector<int> idVec( ids[0] );
 481     for ( int i = 0; i < ids[0]; ++i )
 482       idVec[ i ] = (int) ids[ i + 1];
 483     return mesh->AddPolygonalFaceWithID( idVec, ID );
 484   }
 485
 486   typedef SMDS_MeshElement* (* PAddElemFun) (cgsize_t* ids, SMESHDS_Mesh* mesh, int ID);
 487
 488   //================================================================================
 489   /*!
 490    * \brief Return an array of functions each adding an element of a particular type
 491    */
 492   //================================================================================
 493
 494   PAddElemFun* getAddElemFunTable()
 495   {
 496     static vector< PAddElemFun > funVec;
 497     if ( funVec.empty() )
 498     {
 499       funVec.resize( NofValidElementTypes, (PAddElemFun)0 );
 500       funVec[ CGNS_ENUMV( NODE     )] = add_0D      ;
 501       funVec[ CGNS_ENUMV( BAR_2    )] = add_BAR_2   ;
 502       funVec[ CGNS_ENUMV( BAR_3    )] = add_BAR_3   ;
 503       funVec[ CGNS_ENUMV( TRI_3    )] = add_TRI_3   ;
 504       funVec[ CGNS_ENUMV( TRI_6    )] = add_TRI_6   ;
 505       funVec[ CGNS_ENUMV( QUAD_4   )] = add_QUAD_4  ;
 506       funVec[ CGNS_ENUMV( QUAD_8   )] = add_QUAD_8  ;
 507       funVec[ CGNS_ENUMV( QUAD_9   )] = add_QUAD_8  ;
 508       funVec[ CGNS_ENUMV( TETRA_4  )] = add_TETRA_4 ;
 509       funVec[ CGNS_ENUMV( TETRA_10 )] = add_TETRA_10;
 510       funVec[ CGNS_ENUMV( PYRA_5   )] = add_PYRA_5  ;
 511       funVec[ CGNS_ENUMV( PYRA_13  )] = add_PYRA_13 ;
 512       funVec[ CGNS_ENUMV( PYRA_14  )] = add_PYRA_13 ;
 513       funVec[ CGNS_ENUMV( PENTA_6  )] = add_PENTA_6 ;
 514       funVec[ CGNS_ENUMV( PENTA_15 )] = add_PENTA_15;
 515       funVec[ CGNS_ENUMV( PENTA_18 )] = add_PENTA_15;
 516       funVec[ CGNS_ENUMV( HEXA_8   )] = add_HEXA_8  ;
 517       funVec[ CGNS_ENUMV( HEXA_20  )] = add_HEXA_20 ;
 518       funVec[ CGNS_ENUMV( HEXA_27  )] = add_HEXA_20 ;
 519       funVec[ CGNS_ENUMV( NGON_n   )] = add_NGON    ;
 520     }
 521     return &funVec[0];
 522   }
 523
 524   //================================================================================
 525   /*!
 526    * \brief Finds an existing boundary element
 527    */
 528   //================================================================================
 529
 530   const SMDS_MeshElement* findElement(const cgsize_t*     nodeIDs,
 531                                       const int           nbNodes,
 532                                       const SMESHDS_Mesh* mesh)
 533   {
 534     const SMDS_MeshNode* nn[4]; // look for quad4 or seg2
 535     if (( nn[0] = mesh->FindNode( nodeIDs[0] )))
 536     {
 537       SMDSAbs_ElementType eType = nbNodes==4 ? SMDSAbs_Face : SMDSAbs_Edge;
 538       SMDS_ElemIteratorPtr eIt = nn[0]->GetInverseElementIterator( eType );
 539       if ( eIt->more() )
 540         for ( int i = 1; i < nbNodes; ++i )
 541           nn[i] = mesh->FindNode( nodeIDs[i] );
 542       while ( eIt->more() )
 543       {
 544         const SMDS_MeshElement* e = eIt->next();
 545         if ( e->NbNodes() == nbNodes )
 546         {
 547           bool elemOK = true;
 548           for ( int i = 1; i < nbNodes && elemOK; ++i )
 549             elemOK = ( e->GetNodeIndex( nn[i] ) >= 0 );
 550           if ( elemOK )
 551             return e;
 552         }
 553       }
 554     }
 555     return 0;
 556   }
 557
 558 } // namespace
 559
 560 //================================================================================
 561 /*!
 562  * \brief Perform reading a myMeshId-th mesh
 563  */
 564 //================================================================================
 565
 566 Driver_Mesh::Status DriverCGNS_Read::Perform()
 567 {
 568   myErrorMessages.clear();
 569
 570   Status aResult;
 571   if (( aResult = open() ) != DRS_OK )
 572     return aResult;
 573
 574   // read nb of meshes (CGNSBase_t)
 575   if ( myMeshId < 0 || myMeshId >= GetNbMeshes(aResult))
 576     return addMessage( SMESH_Comment("Invalid mesh index :") << myMeshId );
 577
 578   // read a name and a dimension of the mesh
 579   const int cgnsBase = myMeshId + 1;
 580   char meshName[CGNS_NAME_SIZE];
 581   int meshDim, spaceDim;
 582   if ( cg_base_read( _fn, cgnsBase, meshName, &meshDim, &spaceDim) != CG_OK )
 583     return addMessage( cg_get_error() );
 584
 585   if ( spaceDim < 1 || spaceDim > 3 )
 586     return addMessage( SMESH_Comment("Invalid space dimension: ") << spaceDim
 587                        << " in mesh '" << meshName << "'");
 588
 589   myMeshName = meshName;
 590
 591   // read nb of domains (Zone_t) in the mesh
 592   int nbZones = 0;
 593   if ( cg_nzones (_fn, cgnsBase, &nbZones) != CG_OK )
 594     return addMessage( cg_get_error() );
 595
 596   if ( nbZones < 1 )
 597     return addMessage( SMESH_Comment("Empty mesh: '") << meshName << "'");
 598
 599   // read the domains (zones)
 600   // ------------------------
 601   map< string, TZoneData > zonesByName;
 602   char name[CGNS_NAME_SIZE];
 603   cgsize_t sizes[NB_ZONE_SIZE_VAL];
 604   memset(sizes, 0, NB_ZONE_SIZE_VAL * sizeof(cgsize_t));
 605
 606   const SMDS_MeshInfo& meshInfo = myMesh->GetMeshInfo();
 607   int groupID = myMesh->GetGroups().size();
 608
 609   for ( int iZone = 1; iZone <= nbZones; ++iZone )
 610   {
 611     // size and name of a zone
 612     if ( cg_zone_read( _fn, cgnsBase, iZone, name, sizes) != CG_OK) {
 613       addMessage( cg_get_error() );
 614       continue;
 615     }
 616     TZoneData& zone = zonesByName[ name ];
 617     zone._id          = iZone;
 618     zone._nodeIdShift = meshInfo.NbNodes();
 619     zone._elemIdShift = meshInfo.NbElements();
 620     zone.SetSizeAndDim( sizes, meshDim );
 621
 622     // mesh type of the zone
 623     if ( cg_zone_type ( _fn, cgnsBase, iZone, &zone._type) != CG_OK) {
 624       addMessage( cg_get_error() );
 625       continue;
 626     }
 627
 628     switch ( zone._type )
 629     {
 630     case CGNS_ENUMV( Unstructured ):
 631     case CGNS_ENUMV( Structured ):
 632       break;
 633     case CGNS_ENUMV( ZoneTypeNull ):
 634       addMessage( "Meshes with ZoneTypeNull are not supported");
 635       continue;
 636     case CGNS_ENUMV( ZoneTypeUserDefined ):
 637       addMessage( "Meshes with ZoneTypeUserDefined are not supported");
 638       continue;
 639     default:
 640       addMessage( "Unknown ZoneType_t");
 641       continue;
 642     }
 643
 644     // -----------
 645     // Read nodes
 646     // -----------
 647
 648     if ( cg_ncoords( _fn, cgnsBase, iZone, &spaceDim) != CG_OK ) {
 649       addMessage( cg_get_error() );
 650       continue;
 651     }
 652     if ( spaceDim < 1 ) {
 653       addMessage( SMESH_Comment("No coordinates defined in zone ")
 654                   << iZone << " of Mesh " << myMeshId );
 655       continue;
 656     }
 657     // read coordinates
 658
 659     cgsize_t rmin[3] = {1,1,1}; // range of nodes to read
 660     cgsize_t rmax[3] = {1,1,1};
 661     int nbNodes = rmax[0] = zone._sizes[0];
 662     if ( zone.IsStructured())
 663       for ( int i = 1; i < meshDim; ++i )
 664         nbNodes *= rmax[i] = zone._sizes[i];
 665
 666     vector<double> coords[3];
 667     for ( int c = 1; c <= spaceDim; ++c)
 668     {
 669       coords[c-1].resize( nbNodes );
 670
 671       CGNS_ENUMV( DataType_t ) type;
 672       if ( cg_coord_info( _fn, cgnsBase, iZone, c, &type, name) != CG_OK ||
 673            cg_coord_read( _fn, cgnsBase, iZone, name, CGNS_ENUMV(RealDouble),
 674                           rmin, rmax, (void*)&(coords[c-1][0])) != CG_OK)
 675       {
 676         addMessage( cg_get_error() );
 677         coords[c-1].clear();
 678         break;
 679       }
 680     }
 681     if ( coords[ spaceDim-1 ].empty() )
 682       continue; // there was an error while reading coordinates
 683
 684     // fill coords with zero if spaceDim < 3
 685     for ( int c = 2; c <= 3; ++c)
 686       if ( coords[ c-1 ].empty() )
 687         coords[ c-1 ].resize( nbNodes, 0.0 );
 688
 689     // create nodes
 690     try {
 691       for ( int i = 0; i < nbNodes; ++i )
 692         myMesh->AddNodeWithID( coords[0][i], coords[1][i], coords[2][i], i+1+zone._nodeIdShift );
 693     }
 694     catch ( std::exception& exc ) // expect std::bad_alloc
 695     {
 696       addMessage( exc.what() );
 697       break;
 698     }
 699
 700     // Read connectivity between zones. Nodes of the zone interface will be
 701     // replaced withing the zones read later
 702     string err = zone.ReadZonesConnection( _fn, cgnsBase, zonesByName );
 703     if ( !err.empty() )
 704       addMessage( err );
 705
 706     // --------------
 707     // Read elements
 708     // --------------
 709     if ( zone.IsStructured())
 710     {
 711       int nbI = zone._sizeX - 1, nbJ = zone._sizeY - 1, nbK = zone._sizeZ - 1;
 712       cgsize_t nID[8];
 713       if ( meshDim > 2 && nbK > 0 )
 714       {
 715         for ( int k = 1; k <= nbK; ++k )
 716           for ( int j = 1; j <= nbJ; ++j )
 717             for ( int i = 1; i <= nbI; ++i )
 718             {
 719               zone.CellNodes( i, j, k, nID );
 720               zone.ReplaceNodes( nID, 8 );
 721               myMesh->AddVolumeWithID(nID[0],nID[1],nID[2],nID[3],nID[4],nID[5],nID[6],nID[7],
 722                                       meshInfo.NbElements()+1);
 723             }
 724       }
 725       else if ( meshDim > 1 && nbJ > 0 )
 726       {
 727         for ( int j = 1; j <= nbJ; ++j )
 728           for ( int i = 1; i <= nbI; ++i )
 729           {
 730             zone.CellNodes( i, j, nID );
 731             zone.ReplaceNodes( nID, 4 );
 732             myMesh->AddFaceWithID(nID[0],nID[1],nID[2],nID[3], meshInfo.NbElements()+1);
 733           }
 734       }
 735       else if ( meshDim > 0 && nbI > 0 )
 736       {
 737         nID[0] = zone.NodeID( 1, 0, 0 );
 738         for ( int i = 1; i <= nbI; ++i, ++nID[0] )
 739         {
 740           nID[1] = nID[0]+1;
 741           zone.ReplaceNodes( nID, 2 );
 742           myMesh->AddEdgeWithID(nID[0],nID[1], meshInfo.NbElements()+1);
 743         }
 744       }
 745     }
 746     else
 747     {
 748       // elements can be stored in different sections each dedicated to one element type
 749       int nbSections = 0;
 750       if ( cg_nsections( _fn, cgnsBase, iZone, &nbSections) != CG_OK)
 751       {
 752         addMessage( cg_get_error() );
 753         continue;
 754       }
 755       PAddElemFun* addElemFuns = getAddElemFunTable(), curAddElemFun = 0;
 756       int nbNotSuppElem = 0; // nb elements of not supported types
 757       bool polyhedError = false; // error at polyhedron creation
 758
 759       // read element data
 760
 761       CGNS_ENUMT( ElementType_t ) elemType;
 762       cgsize_t start, end; // range of ids of elements of a zone
 763       cgsize_t eDataSize = 0;
 764       int nbBnd, parent_flag;
 765       for ( int iSec = 1; iSec <= nbSections; ++iSec )
 766       {
 767         if ( cg_section_read( _fn, cgnsBase, iZone, iSec, name, &elemType,
 768                               &start, &end, &nbBnd, &parent_flag) != CG_OK ||
 769              cg_ElementDataSize( _fn, cgnsBase, iZone, iSec, &eDataSize ) != CG_OK )
 770         {
 771           addMessage( cg_get_error() );
 772           continue;
 773         }
 774         vector< cgsize_t > elemData( eDataSize );
 775         if ( cg_elements_read( _fn, cgnsBase, iZone, iSec, &elemData[0], NULL ) != CG_OK )
 776         {
 777           addMessage( cg_get_error() );
 778           continue;
 779         }
 780         // store elements
 781
 782         int pos = 0, cgnsNbNodes = 0, elemID = start + zone._elemIdShift;
 783         cg_npe( elemType, &cgnsNbNodes ); // get nb nodes by element type
 784         curAddElemFun = addElemFuns[ elemType ];
 785         SMDS_MeshElement* newElem = 0;
 786         const SMDS_MeshElement* face;
 787
 788         while ( pos < eDataSize )
 789         {
 790           CGNS_ENUMT( ElementType_t ) currentType = elemType;
 791           if ( currentType == CGNS_ENUMV( MIXED )) {
 792             //ElementConnectivity = Etype1, Node11, Node21, ... NodeN1,
 793             //                      Etype2, Node12, Node22, ... NodeN2,
 794             //                      ...
 795             //                      EtypeM, Node1M, Node2M, ... NodeNM
 796             currentType = (CGNS_ENUMT(ElementType_t)) elemData[ pos++ ];
 797             cg_npe( currentType, &cgnsNbNodes );
 798             curAddElemFun = addElemFuns[ currentType ];
 799           }
 800           if ( cgnsNbNodes < 1 ) // poly elements
 801           {
 802             if ( currentType == CGNS_ENUMV( NFACE_n )) // polyhedron
 803             {
 804               //ElementConnectivity = Nfaces1, Face11, Face21, ... FaceN1,
 805               //                      Nfaces2, Face12, Face22, ... FaceN2,
 806               //                      ...
 807               //                      NfacesM, Face1M, Face2M, ... FaceNM
 808               const int nbFaces = elemData[ pos++ ];
 809               vector<int> quantities( nbFaces );
 810               vector<const SMDS_MeshNode*> nodes, faceNodes;
 811               nodes.reserve( nbFaces * 4 );
 812               for ( int iF = 0; iF < nbFaces; ++iF )
 813               {
 814                 const int faceID = std::abs( elemData[ pos++ ]) + zone._elemIdShift;
 815                 if (( face = myMesh->FindElement( faceID )) && face->GetType() == SMDSAbs_Face )
 816                 {
 817                   const bool reverse = ( elemData[ pos-1 ] < 0 );
 818                   const int    iQuad = face->IsQuadratic() ? 1 : 0;
 819                   SMDS_ElemIteratorPtr nIter = face->interlacedNodesElemIterator();
 820                   faceNodes.assign( SMDS_MeshElement::iterator( nIter ),
 821                                     SMDS_MeshElement::iterator());
 822                   if ( iQuad && reverse )
 823                     nodes.push_back( faceNodes[0] );
 824                   if ( reverse )
 825                     nodes.insert( nodes.end(), faceNodes.rbegin(), faceNodes.rend() - iQuad );
 826                   else
 827                     nodes.insert( nodes.end(), faceNodes.begin(), faceNodes.end() );
 828
 829                   quantities[ iF ] = face->NbNodes();
 830                 }
 831                 else {
 832                   polyhedError = true;
 833                   break;
 834                 }
 835               }
 836               if ( quantities.back() )
 837               {
 838                 myMesh->AddPolyhedralVolumeWithID( nodes, quantities, elemID );
 839               }
 840             }
 841             else if ( currentType == CGNS_ENUMV( NGON_n )) // polygon
 842             {
 843               // ElementConnectivity = Nnodes1, Node11, Node21, ... NodeN1,
 844               //                       Nnodes2, Node12, Node22, ... NodeN2,
 845               //                       ...
 846               //                       NnodesM, Node1M, Node2M, ... NodeNM
 847               const int nbNodes = elemData[ pos ];
 848               zone.ReplaceNodes( &elemData[pos+1], nbNodes, zone._nodeIdShift );
 849               newElem = add_NGON( &elemData[pos  ], myMesh, elemID );
 850               pos += nbNodes + 1;
 851             }
 852           }
 853           else // standard elements
 854           {
 855             zone.ReplaceNodes( &elemData[pos], cgnsNbNodes, zone._nodeIdShift );
 856             newElem = curAddElemFun( &elemData[pos], myMesh, elemID );
 857             pos += cgnsNbNodes;
 858             nbNotSuppElem += int( newElem && newElem->NbNodes() != cgnsNbNodes );
 859           }
 860           elemID++;
 861
 862         } // loop on elemData
 863       } // loop on cgns sections
 864
 865       if ( nbNotSuppElem > 0 )
 866         addMessage( SMESH_Comment(nbNotSuppElem) << " elements of not supported types"
 867                     << " have beem converted to close types");
 868       if ( polyhedError )
 869         addMessage( "Some polyhedral elements have been skipped due to internal(?) errors" );
 870
 871     } // reading unstructured elements
 872
 873     zone._nbNodes = meshInfo.NbNodes() - zone._nodeIdShift;
 874     zone._nbElems = meshInfo.NbElements() - zone._elemIdShift;
 875
 876     // -------------------------------------------
 877     // Read Boundary Conditions into SMESH groups
 878     // -------------------------------------------
 879     int nbBC = 0;
 880     if ( cg_nbocos( _fn, cgnsBase, iZone, &nbBC) == CG_OK )
 881     {
 882       CGNS_ENUMT( BCType_t ) bcType;
 883       CGNS_ENUMT( PointSetType_t ) psType;
 884       CGNS_ENUMT( DataType_t ) normDataType;
 885       cgsize_t nbPnt, normFlag;
 886       int normIndex[3], nbDS;
 887       for ( int iBC = 1; iBC <= nbBC; ++iBC )
 888       {
 889         if ( cg_boco_info( _fn, cgnsBase, iZone, iBC, name, &bcType, &psType,
 890                            &nbPnt, normIndex, &normFlag, &normDataType, &nbDS ) != CG_OK )
 891         {
 892           addMessage( cg_get_error() );
 893           continue;
 894         }
 895         vector< cgsize_t > ids( nbPnt * zone.IndexSize() );
 896         CGNS_ENUMT( GridLocation_t ) location;
 897         if ( cg_boco_read( _fn, cgnsBase, iZone, iBC, &ids[0], NULL ) != CG_OK ||
 898              cg_boco_gridlocation_read( _fn, cgnsBase, iZone, iBC, &location) != CG_OK )
 899         {
 900           addMessage( cg_get_error() );
 901           continue;
 902         }
 903         SMDSAbs_ElementType elemType = SMDSAbs_All;
 904         switch ( location ) {
 905         case CGNS_ENUMV( Vertex      ): elemType = SMDSAbs_Node; break;
 906         case CGNS_ENUMV( FaceCenter  ): elemType = SMDSAbs_Face; break;
 907         case CGNS_ENUMV( IFaceCenter ): elemType = SMDSAbs_Face; break;
 908         case CGNS_ENUMV( JFaceCenter ): elemType = SMDSAbs_Face; break;
 909         case CGNS_ENUMV( KFaceCenter ): elemType = SMDSAbs_Face; break;
 910         case CGNS_ENUMV( EdgeCenter  ): elemType = SMDSAbs_Edge; break;
 911         default:;
 912         }
 913         SMESHDS_Group* group = new SMESHDS_Group ( groupID++, myMesh, elemType );
 914         myMesh->AddGroup( group );
 915         SMESH_Comment groupName( name ); groupName << " " << cg_BCTypeName( bcType );
 916         group->SetStoreName( groupName.c_str() );
 917         SMDS_MeshGroup& groupDS = group->SMDSGroup();
 918
 919         if ( elemType == SMDSAbs_Node )
 920         {
 921           if ( zone.IsStructured() )
 922           {
 923             vector< cgsize_t > nodeIds;
 924             if ( psType == CGNS_ENUMV( PointRange ))
 925             {
 926               // nodes are given as (ijkMin, ijkMax)
 927               TPointRangeIterator idIt( & ids[0], meshDim );
 928               nodeIds.reserve( idIt.Size() );
 929               while ( idIt.More() )
 930                 nodeIds.push_back( zone.NodeID( idIt.Next() ));
 931             }
 932             else
 933             {
 934               // nodes are given as (ijk1, ijk2, ..., ijkN)
 935               nodeIds.reserve( ids.size() / meshDim );
 936               for ( size_t i = 0; i < ids.size(); i += meshDim )
 937                 nodeIds.push_back( zone.NodeID( ids[i], ids[i+1], ids[i+2] ));
 938             }
 939             ids.swap( nodeIds );
 940           }
 941           else if ( zone._nodeIdShift )
 942           {
 943             for ( size_t i = 0; i < ids.size(); ++i )
 944               ids[i] += zone._nodeIdShift;
 945           }
 946           zone.ReplaceNodes( &ids[0], ids.size() );
 947
 948           for ( size_t i = 0; i < ids.size(); ++i )
 949             if ( const SMDS_MeshNode* n = myMesh->FindNode( ids[i] ))
 950               groupDS.Add( n );
 951         }
 952         else // BC applied to elements
 953         {
 954           if ( zone.IsStructured() )
 955           {
 956             int axis = 0; // axis perpendiculaire to which boundary elements are oriented
 957             if ( ids.size() >= meshDim * 2 )
 958             {
 959               for ( ; axis < meshDim; ++axis )
 960                 if ( ids[axis] - ids[axis+meshDim] == 0 )
 961                   break;
 962             }
 963             else
 964             {
 965               for ( ; axis < meshDim; ++axis )
 966                 if ( normIndex[axis] != 0 )
 967                   break;
 968             }
 969             if ( axis == meshDim )
 970             {
 971               addMessage( SMESH_Comment("Invalid NormalIndex in BC ") << name );
 972               continue;
 973             }
 974             const int nbElemNodesByDim[] = { 1, 2, 4, 8 };
 975             const int nbElemNodes = nbElemNodesByDim[ meshDim ];
 976
 977             if ( psType == CGNS_ENUMV( PointRange ) ||
 978                  psType == CGNS_ENUMV( ElementRange ))
 979             {
 980               // elements are given as (ijkMin, ijkMax)
 981               typedef void (TZoneData::*PGetNodesFun)( const gp_XYZ& ijk, cgsize_t* ids ) const;
 982               PGetNodesFun getNodesFun = 0;
 983               if ( elemType == SMDSAbs_Face  && meshDim == 3 )
 984                 switch ( axis ) {
 985                 case 0: getNodesFun = & TZoneData::IFaceNodes;
 986                 case 1: getNodesFun = & TZoneData::JFaceNodes;
 987                 case 2: getNodesFun = & TZoneData::KFaceNodes;
 988                 }
 989               else if ( elemType == SMDSAbs_Edge && meshDim == 2 )
 990                 switch ( axis ) {
 991                 case 0: getNodesFun = & TZoneData::IEdgeNodes;
 992                 case 1: getNodesFun = & TZoneData::JEdgeNodes;
 993                 }
 994               if ( !getNodesFun )
 995               {
 996                 addMessage( SMESH_Comment("Unsupported BC location in BC ") << name
 997                             << " " << cg_GridLocationName( location )
 998                             << " in " << meshDim << " mesh");
 999                 continue;
1000               }
1001               TPointRangeIterator rangeIt( & ids[0], meshDim );
1002               vector< cgsize_t > elemNodeIds( rangeIt.Size() * nbElemNodes );
1003               for ( int i = 0; rangeIt.More(); i+= nbElemNodes )
1004                 (zone.*getNodesFun)( rangeIt.Next(), &elemNodeIds[i] );
1005
1006               ids.swap( elemNodeIds );
1007             }
1008             else
1009             {
1010               // elements are given as (ijk1, ijk2, ..., ijkN)
1011               typedef void (TZoneData::*PGetNodesFun)( int i, int j, int k, cgsize_t* ids ) const;
1012               PGetNodesFun getNodesFun = 0;
1013               if ( elemType == SMDSAbs_Face )
1014                 switch ( axis ) {
1015                 case 0: getNodesFun = & TZoneData::IFaceNodes;
1016                 case 1: getNodesFun = & TZoneData::JFaceNodes;
1017                 case 2: getNodesFun = & TZoneData::KFaceNodes;
1018                 }
1019               else if ( elemType == SMDSAbs_Edge && meshDim == 2 )
1020                 switch ( axis ) {
1021                 case 0: getNodesFun = & TZoneData::IEdgeNodes;
1022                 case 1: getNodesFun = & TZoneData::JEdgeNodes;
1023                 }
1024               if ( !getNodesFun )
1025               {
1026                 addMessage( SMESH_Comment("Unsupported BC location in BC ") << name
1027                             << " " << cg_GridLocationName( location )
1028                             << " in " << meshDim << " mesh");
1029                 continue;
1030               }
1031               vector< cgsize_t > elemNodeIds( ids.size()/meshDim * nbElemNodes );
1032               for ( size_t i = 0, j = 0; i < ids.size(); i += meshDim, j += nbElemNodes )
1033                 (zone.*getNodesFun)( ids[i], ids[i+1], ids[i+2], &elemNodeIds[j] );
1034
1035               ids.swap( elemNodeIds );
1036             }
1037             zone.ReplaceNodes( &ids[0], ids.size() );
1038
1039             PAddElemFun addElemFun = 0;
1040             switch ( meshDim ) {
1041             case 1: addElemFun = & add_BAR_2;
1042             case 2: addElemFun = & add_QUAD_4;
1043             case 3: addElemFun = & add_HEXA_8;
1044             }
1045             int elemID = meshInfo.NbElements();
1046             const SMDS_MeshElement* elem = 0;
1047             for ( size_t i = 0; i < ids.size(); i += nbElemNodes )
1048             {
1049               if ( iZone == 1 || !( elem = findElement( &ids[i], nbElemNodes, myMesh )))
1050                 elem = addElemFun( &ids[i], myMesh, ++elemID );
1051               groupDS.Add( elem );
1052             }
1053           }
1054           else // unstructured zone
1055           {
1056             if ( zone._elemIdShift )
1057               for ( size_t i = 0; i < ids.size(); ++i )
1058                 ids[i] += zone._elemIdShift;
1059
1060             if ( psType == CGNS_ENUMV( PointRange ) && ids.size() == 2 )
1061             {
1062               for ( size_t i = ids[0]; i <= ids[1]; ++i )
1063                 if ( const SMDS_MeshElement* e = myMesh->FindElement( i ))
1064                   groupDS.Add( e );
1065             }
1066             else
1067             {
1068               for ( size_t i = 0; i < ids.size(); ++i )
1069                 if ( const SMDS_MeshElement* e = myMesh->FindElement( ids[i] ))
1070                   groupDS.Add( e );
1071             }
1072           }
1073         } // end "BC applied to elements"
1074
1075         // to have group type according to a real elem type
1076         group->SetType( groupDS.GetType() );
1077
1078       } // loop on BCs of the zone
1079     }
1080     else
1081     {
1082       addMessage( cg_get_error() );
1083     }
1084   } // loop on the zones of a mesh
1085
1086
1087   // ------------------------------------------------------------------------
1088   // Make groups for multiple zones and remove free nodes at zone interfaces
1089   // ------------------------------------------------------------------------
1090   map< string, TZoneData >::iterator nameZoneIt = zonesByName.begin();
1091   for ( ; nameZoneIt != zonesByName.end(); ++nameZoneIt )
1092   {
1093     TZoneData& zone = nameZoneIt->second;
1094     if ( zone._nbElems == 0 ) continue;
1095     if ( zone._nbElems == meshInfo.NbElements() ) break; // there is only one non-empty zone
1096
1097     // make a group
1098     SMDSAbs_ElementType elemType = myMesh->GetElementType( zone._elemIdShift + 1,
1099                                                            /*iselem=*/true );
1100     SMESHDS_Group* group = new SMESHDS_Group ( groupID++, myMesh, elemType );
1101     myMesh->AddGroup( group );
1102     group->SetStoreName( nameZoneIt->first.c_str() );
1103     SMDS_MeshGroup& groupDS = group->SMDSGroup();
1104
1105     for ( int i = 1; i <= zone._nbElems; ++i )
1106       if ( const SMDS_MeshElement* e = myMesh->FindElement( i + zone._elemIdShift ))
1107         groupDS.Add( e );
1108
1109     // remove free nodes
1110     map< int, int >::iterator nnRmKeepIt = zone._nodeReplacementMap.begin();
1111     for ( ; nnRmKeepIt != zone._nodeReplacementMap.end(); ++nnRmKeepIt )
1112       if ( const SMDS_MeshNode* n = myMesh->FindNode( nnRmKeepIt->first ))
1113         if ( n->NbInverseElements() == 0 )
1114           myMesh->RemoveFreeNode( n, (SMESHDS_SubMesh *)0, /*fromGroups=*/false );
1115   }
1116
1117   aResult = myErrorMessages.empty() ? DRS_OK : DRS_WARN_SKIP_ELEM;
1118
1119   return aResult;
1120 }
1121
1122 //================================================================================
1123 /*!
1124  * \brief Constructor
1125  */
1126 //================================================================================
1127
1128 DriverCGNS_Read::DriverCGNS_Read()
1129 {
1130   _fn = -1;
1131 }
1132 //================================================================================
1133 /*!
1134  * \brief Close the cgns file at destruction
1135  */
1136 //================================================================================
1137
1138 DriverCGNS_Read::~DriverCGNS_Read()
1139 {
1140   if ( _fn > 0 )
1141     cg_close( _fn );
1142 }
1143
1144 //================================================================================
1145 /*!
1146  * \brief Opens myFile
1147  */
1148 //================================================================================
1149
1150 Driver_Mesh::Status DriverCGNS_Read::open()
1151 {
1152   if ( _fn < 0 )
1153   {
1154
1155 #ifdef CG_MODE_READ
1156     int res = cg_open(myFile.c_str(), CG_MODE_READ, &_fn);
1157 #else
1158     int res = cg_open(myFile.c_str(), MODE_READ, &_fn);
1159 #endif
1160     if ( res != CG_OK)
1161     {
1162       addMessage( cg_get_error(), /*fatal = */true );
1163     }
1164   }
1165   return _fn >= 0 ? DRS_OK : DRS_FAIL;
1166 }
1167
1168 //================================================================================
1169 /*!
1170  * \brief Reads nb of meshes in myFile
1171  */
1172 //================================================================================
1173
1174 int DriverCGNS_Read::GetNbMeshes(Status& theStatus)
1175 {
1176   if (( theStatus = open()) != DRS_OK )
1177     return 0;
1178
1179   int nbases = 0;
1180   if(cg_nbases( _fn, &nbases) != CG_OK)
1181     theStatus = addMessage( cg_get_error(), /*fatal = */true );
1182
1183   return nbases;
1184 }