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