src/DriverCGNS/DriverCGNS_Read.cxx

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