src/DriverCGNS/DriverCGNS_Read.cxx

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