src/DriverCGNS/DriverCGNS_Read.cxx

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