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