Salome HOME
Copyright update 2022
[modules/smesh.git] / src / DriverCGNS / DriverCGNS_Read.cxx
index 31c7ef6db08ff729999723bc4299cd8bd6145f80..2b28c9f6ba2afc6e97f191e57131263f65d5da3f 100644 (file)
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2019  CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2022  CEA/DEN, EDF R&D, OPEN CASCADE
 //
 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
@@ -22,6 +22,8 @@
 // File      : DriverCGNS_Read.cxx
 // Created   : Thu Jun 30 10:33:31 2011
 // Author    : Edward AGAPOV (eap)
+//#define _DEBUG_
+#include <utilities.h>
 
 #include "DriverCGNS_Read.hxx"
 
 #include "SMESH_Comment.hxx"
 #include "SMESH_TypeDefs.hxx"
 
+#include <smIdType.hxx>
+
 #include <gp_XYZ.hxx>
 
 #include <cgnslib.h>
 
 #include <map>
 
+
 #if CGNS_VERSION < 3100
 # define cgsize_t int
 #endif
@@ -56,9 +61,9 @@ namespace
   struct TZoneData
   {
     int                    _id;
-    int                    _nodeIdShift; // nb nodes in previously read zones
-    int                    _elemIdShift; // nb faces in previously read zones
-    int                    _nbNodes, _nbElems;
+    smIdType               _nodeIdShift; // nb nodes in previously read zones
+    smIdType               _elemIdShift; // nb faces in previously read zones
+    smIdType               _nbNodes, _nbElems;
     int                    _meshDim;
     int                    _sizeX, _sizeY, _sizeZ, _nbCells; // structured
     cgsize_t               _sizes[NB_ZONE_SIZE_VAL];
@@ -131,12 +136,12 @@ namespace
       ids[2] = ids[0] + _sizeX + 1;
       ids[3] = ids[0] + ( k==_sizeZ ? _sizeX : 1);
     }
-    void IEdgeNodes( int i, int j, int k, cgsize_t* ids ) const // edge perpendiculaire to X (2D)
+    void IEdgeNodes( int i, int j, int /*k*/, cgsize_t* ids ) const // edge perpendiculaire to X (2D)
     {
       ids[0] = NodeID( i, j, 0 );
       ids[1] = ids[0] + _sizeX;
     }
-    void JEdgeNodes( int i, int j, int k, cgsize_t* ids ) const
+    void JEdgeNodes( int i, int j, int /*k*/, cgsize_t* ids ) const
     {
       ids[0] = NodeID( i, j, 0 );
       ids[1] = ids[0] + 1;
@@ -162,9 +167,9 @@ namespace
     int _beg[3], _end[3], _cur[3], _dir[3], _dim;
     bool _more;
   public:
-    TPointRangeIterator( const cgsize_t* range, int dim ):_dim(dim)
+    TPointRangeIterator( const cgsize_t* range, int dim ):
+      _beg{0,0,0}, _end{0,0,0}, _cur{0,0,0}, _dir{0,0,0}, _dim(dim), _more(false)
     {
-      _more = false;
       for ( int i = 0; i < dim; ++i )
       {
         _beg[i] = range[i];
@@ -175,8 +180,6 @@ namespace
         if ( _end[i] - _beg[i] )
           _more = true;
       }
-//       for ( int i = dim; i < 3; ++i )
-//         _cur[i] = _beg[i] = _end[i] = _dir[i] = 0;
     }
     bool More() const
     {
@@ -214,7 +217,8 @@ namespace
    */
   //================================================================================
 
-  bool isEqualNodes( const int* nIds1, const int* nIds2, int nbNodes, SMESHDS_Mesh* mesh )
+  bool isEqualNodes( const cgsize_t* nIds1, const cgsize_t* nIds2, size_t nbNodes,
+                     SMESHDS_Mesh* mesh )
   {
     if ( nbNodes > 0 )
     {
@@ -314,7 +318,7 @@ namespace
             }
             // check if range and donorRange describe the same nodes
             {
-              int ids1[2], ids2[2], nbN = 0;
+              cgsize_t ids1[2], ids2[2], nbN = 0;
               TPointRangeIterator rangeIt1bis( range, _meshDim );
               index1 = rangeIt1bis.Next();
               index2 = T * ( index1 - begin1 ) + begin2;
@@ -543,12 +547,16 @@ namespace
                                   ids[14],ids[13],ids[19],ids[18],ids[17],ids[16],
                                   ids[20],ids[24],ids[23],ids[22],ids[21],ids[25],ids[26], ID );
   }
-  SMDS_MeshElement* add_NGON(cgsize_t* ids, SMESHDS_Mesh* mesh, int ID)
+  SMDS_MeshElement* add_NGON(cgsize_t* ids, int nbNodes, SMESHDS_Mesh* mesh, int ID)
   {
-    vector<int> idVec( ids[0] );
-    for ( int i = 0; i < ids[0]; ++i )
-      idVec[ i ] = (int) ids[ i + 1];
-    return mesh->AddPolygonalFaceWithID( idVec, ID );
+#if CGNS_VERSION < 4000
+    nbNodes = ids[0];
+    ++ids;
+#endif
+    vector<smIdType> idVec( nbNodes );
+    for ( int i = 0; i < nbNodes; ++i )
+      idVec[ i ] = ToSmIdType( ids[ i ]);
+    return mesh->AddPolygonalFaceWithID( idVec, ToSmIdType(ID) );
   }
 
   typedef SMDS_MeshElement* (* PAddElemFun) (cgsize_t* ids, SMESHDS_Mesh* mesh, int ID);
@@ -584,7 +592,7 @@ namespace
       funVec[ CGNS_ENUMV( HEXA_8   )] = add_HEXA_8  ;
       funVec[ CGNS_ENUMV( HEXA_20  )] = add_HEXA_20 ;
       funVec[ CGNS_ENUMV( HEXA_27  )] = add_HEXA_27 ;
-      funVec[ CGNS_ENUMV( NGON_n   )] = add_NGON    ;
+      //funVec[ CGNS_ENUMV( NGON_n   )] = add_NGON    ;
     }
     return &funVec[0];
   }
@@ -633,6 +641,7 @@ namespace
 
 Driver_Mesh::Status DriverCGNS_Read::Perform()
 {
+  MESSAGE("DriverCGNS_Read::Perform");
   myErrorMessages.clear();
 
   Status aResult;
@@ -642,6 +651,7 @@ Driver_Mesh::Status DriverCGNS_Read::Perform()
   // read nb of meshes (CGNSBase_t)
   if ( myMeshId < 0 || myMeshId >= GetNbMeshes(aResult))
     return addMessage( SMESH_Comment("Invalid mesh index :") << myMeshId );
+  MESSAGE("NbMeshes: " << GetNbMeshes(aResult));
 
   // read a name and a dimension of the mesh
   const int cgnsBase = myMeshId + 1;
@@ -655,6 +665,8 @@ Driver_Mesh::Status DriverCGNS_Read::Perform()
                        << " in mesh '" << meshName << "'");
 
   myMeshName = meshName;
+  MESSAGE("myMeshName: " << myMeshName);
+
 
   // read nb of domains (Zone_t) in the mesh
   int nbZones = 0;
@@ -663,6 +675,7 @@ Driver_Mesh::Status DriverCGNS_Read::Perform()
 
   if ( nbZones < 1 )
     return addMessage( SMESH_Comment("Empty mesh: '") << meshName << "'");
+  MESSAGE("nbZones: " << nbZones);
 
   // read the domains (zones)
   // ------------------------
@@ -686,6 +699,7 @@ Driver_Mesh::Status DriverCGNS_Read::Perform()
     zone._nodeIdShift = meshInfo.NbNodes();
     zone._elemIdShift = meshInfo.NbElements();
     zone.SetSizeAndDim( sizes, meshDim );
+    MESSAGE("  zone name: " << name);
 
     // mesh type of the zone
     if ( cg_zone_type ( _fn, cgnsBase, iZone, &zone._type) != CG_OK) {
@@ -696,7 +710,10 @@ Driver_Mesh::Status DriverCGNS_Read::Perform()
     switch ( zone._type )
     {
     case CGNS_ENUMV( Unstructured ):
+      MESSAGE("  zone type: unstructured");
+      break;
     case CGNS_ENUMV( Structured ):
+      MESSAGE("  zone type: structured");
       break;
     case CGNS_ENUMV( ZoneTypeNull ):
       addMessage( "Meshes with ZoneTypeNull are not supported");
@@ -712,7 +729,7 @@ Driver_Mesh::Status DriverCGNS_Read::Perform()
     // -----------
     // Read nodes
     // -----------
-
+    MESSAGE("  Read nodes");
     if ( cg_ncoords( _fn, cgnsBase, iZone, &spaceDim) != CG_OK ) {
       addMessage( cg_get_error() );
       continue;
@@ -724,6 +741,7 @@ Driver_Mesh::Status DriverCGNS_Read::Perform()
     }
     // read coordinates
 
+    MESSAGE("  Read coordinates");
     cgsize_t rmin[3] = {1,1,1}; // range of nodes to read
     cgsize_t rmax[3] = {1,1,1};
     int nbNodes = rmax[0] = zone._sizes[0];
@@ -755,6 +773,7 @@ Driver_Mesh::Status DriverCGNS_Read::Perform()
         coords[ c-1 ].resize( nbNodes, 0.0 );
 
     // create nodes
+    MESSAGE("  create nodes");
     try {
       for ( int i = 0; i < nbNodes; ++i )
         myMesh->AddNodeWithID( coords[0][i], coords[1][i], coords[2][i], i+1+zone._nodeIdShift );
@@ -774,6 +793,7 @@ Driver_Mesh::Status DriverCGNS_Read::Perform()
     // --------------
     // Read elements
     // --------------
+    MESSAGE("  read elements");
     if ( zone.IsStructured())
     {
       int nbI = zone._sizeX - 1, nbJ = zone._sizeY - 1, nbK = zone._sizeZ - 1;
@@ -826,12 +846,14 @@ Driver_Mesh::Status DriverCGNS_Read::Perform()
 
       // read element data
 
+      MESSAGE("  read element data");
       CGNS_ENUMT( ElementType_t ) elemType;
       cgsize_t start, end; // range of ids of elements of a zone
       cgsize_t eDataSize = 0;
       int nbBnd, parent_flag;
       for ( int iSec = 1; iSec <= nbSections; ++iSec )
       {
+        MESSAGE("  section " << iSec << " of " << nbSections);
         if ( cg_section_read( _fn, cgnsBase, iZone, iSec, name, &elemType,
                               &start, &end, &nbBnd, &parent_flag) != CG_OK ||
              cg_ElementDataSize( _fn, cgnsBase, iZone, iSec, &eDataSize ) != CG_OK )
@@ -839,19 +861,40 @@ Driver_Mesh::Status DriverCGNS_Read::Perform()
           addMessage( cg_get_error() );
           continue;
         }
-        vector< cgsize_t > elemData( eDataSize );
-        if ( cg_elements_read( _fn, cgnsBase, iZone, iSec, &elemData[0], NULL ) != CG_OK )
+        vector< cgsize_t > elemData( eDataSize ), polyOffset;
+#if CGNS_VERSION >= 4000
+        if ( elemType == CGNS_ENUMV( MIXED ) ||
+             elemType == CGNS_ENUMV( NGON_n ) ||
+             elemType == CGNS_ENUMV( NFACE_n ))
         {
-          addMessage( cg_get_error() );
-          continue;
+          polyOffset.resize( end - start + 2 );
+          if ( cg_poly_elements_read( _fn, cgnsBase, iZone, iSec,
+                                      elemData.data(), polyOffset.data(), NULL ) != CG_OK )
+          {
+            addMessage( cg_get_error() );
+            continue;
+          }
+        }
+        else
+#endif
+        {
+          if ( cg_elements_read( _fn, cgnsBase, iZone, iSec, elemData.data(), NULL ) != CG_OK )
+          {
+            addMessage( cg_get_error() );
+            continue;
+          }
         }
         // store elements
 
+        MESSAGE("   store elements");
         int pos = 0, cgnsNbNodes = 0, elemID = start + zone._elemIdShift;
+        size_t iElem = 0;
         cg_npe( elemType, &cgnsNbNodes ); // get nb nodes by element type
         curAddElemFun = addElemFuns[ elemType ];
         SMDS_MeshElement* newElem = 0;
         const SMDS_MeshElement* face;
+        vector<int> quantities;
+        vector<const SMDS_MeshNode*> nodes, faceNodes;
 
         while ( pos < eDataSize )
         {
@@ -869,17 +912,25 @@ Driver_Mesh::Status DriverCGNS_Read::Perform()
           {
             if ( currentType == CGNS_ENUMV( NFACE_n )) // polyhedron
             {
-              //ElementConnectivity = Nfaces1, Face11, Face21, ... FaceN1,
-              //                      Nfaces2, Face12, Face22, ... FaceN2,
-              //                      ...
-              //                      NfacesM, Face1M, Face2M, ... FaceNM
-              const int nbFaces = elemData[ pos++ ];
-              vector<int> quantities( nbFaces );
-              vector<const SMDS_MeshNode*> nodes, faceNodes;
-              nodes.reserve( nbFaces * 4 );
+              int nbFaces = 0;
+              if ( polyOffset.empty() )
+                //ElementConnectivity = Nfaces1, Face11, Face21, ... FaceN1,
+                //                      Nfaces2, Face12, Face22, ... FaceN2,
+                //                      ...
+                //                      NfacesM, Face1M, Face2M, ... FaceNM
+                nbFaces = elemData[ pos++ ];
+              else // CGNS_VERSION >= 4000
+                // ElementConnectivity = Face11, Face21, ... FaceN1,
+                //                       Face12, Face22, ... FaceN2,
+                //                       ...
+                //                       Face1M, Face2M, ... FaceNM
+                nbFaces = polyOffset[ iElem + 1 ] - polyOffset[ iElem ];
+
+              quantities.resize( nbFaces ); quantities.back() = 0;
+              nodes.clear();                nodes.reserve( nbFaces * 4 );
               for ( int iF = 0; iF < nbFaces; ++iF )
               {
-                const int faceID = std::abs( elemData[ pos++ ]) + zone._elemIdShift; 
+                const int faceID = std::abs( elemData[ pos++ ]) + zone._elemIdShift;
                 if (( face = myMesh->FindElement( faceID )) && face->GetType() == SMDSAbs_Face )
                 {
                   const bool reverse = ( elemData[ pos-1 ] < 0 );
@@ -898,6 +949,7 @@ Driver_Mesh::Status DriverCGNS_Read::Perform()
                 }
                 else {
                   polyhedError = true;
+                  pos += nbFaces - iF - 1;
                   break;
                 }
               }
@@ -908,14 +960,23 @@ Driver_Mesh::Status DriverCGNS_Read::Perform()
             }
             else if ( currentType == CGNS_ENUMV( NGON_n )) // polygon
             {
-              // ElementConnectivity = Nnodes1, Node11, Node21, ... NodeN1,
-              //                       Nnodes2, Node12, Node22, ... NodeN2,
-              //                       ...
-              //                       NnodesM, Node1M, Node2M, ... NodeNM
-              const int nbNodes = elemData[ pos ];
-              zone.ReplaceNodes( &elemData[pos+1], nbNodes, zone._nodeIdShift );
-              newElem = add_NGON( &elemData[pos  ], myMesh, elemID );
-              pos += nbNodes + 1;
+              int nbNodes;
+              if ( polyOffset.empty() )
+                // ElementConnectivity = Nnodes1, Node11, Node21, ... NodeN1,
+                //                       Nnodes2, Node12, Node22, ... NodeN2,
+                //                       ...
+                //                       NnodesM, Node1M, Node2M, ... NodeNM
+                nbNodes = elemData[ pos ];
+              else // CGNS_VERSION >= 4000
+                // ElementConnectivity = Node11, Node21, ... NodeN1,
+                //                       Node12, Node22, ... NodeN2,
+                //                       ...
+                //                       Node1M, Node2M, ... NodeNM
+                nbNodes = polyOffset[ iElem + 1 ] - polyOffset[ iElem ];
+
+              zone.ReplaceNodes( &elemData[ pos + polyOffset.empty()], nbNodes, zone._nodeIdShift );
+              newElem = add_NGON( &elemData[ pos ], nbNodes, myMesh, elemID );
+              pos += nbNodes + polyOffset.empty();
             }
           }
           else // standard elements
@@ -926,6 +987,7 @@ Driver_Mesh::Status DriverCGNS_Read::Perform()
             nbNotSuppElem += int( newElem && newElem->NbNodes() != cgnsNbNodes );
           }
           elemID++;
+          iElem++;
 
         } // loop on elemData
       } // loop on cgns sections
@@ -944,6 +1006,8 @@ Driver_Mesh::Status DriverCGNS_Read::Perform()
     // -------------------------------------------
     // Read Boundary Conditions into SMESH groups
     // -------------------------------------------
+
+    MESSAGE("  read Boundary Conditions");
     int nbBC = 0;
     if ( cg_nbocos( _fn, cgnsBase, iZone, &nbBC) == CG_OK )
     {
@@ -952,14 +1016,17 @@ Driver_Mesh::Status DriverCGNS_Read::Perform()
       CGNS_ENUMT( DataType_t ) normDataType;
       cgsize_t nbPnt, normFlag;
       int normIndex[3], nbDS;
+      MESSAGE("  nbBC: " << nbBC);
       for ( int iBC = 1; iBC <= nbBC; ++iBC )
       {
+        MESSAGE("  iBC: " << iBC);
         if ( cg_boco_info( _fn, cgnsBase, iZone, iBC, name, &bcType, &psType,
                            &nbPnt, normIndex, &normFlag, &normDataType, &nbDS ) != CG_OK )
         {
           addMessage( cg_get_error() );
           continue;
         }
+        MESSAGE("  iBC info OK: " << iBC);
         vector< cgsize_t > ids( nbPnt * zone.IndexSize() );
         CGNS_ENUMT( GridLocation_t ) location;
         if ( cg_boco_read( _fn, cgnsBase, iZone, iBC, &ids[0], NULL ) != CG_OK ||
@@ -1050,14 +1117,14 @@ Driver_Mesh::Status DriverCGNS_Read::Perform()
               PGetNodesFun getNodesFun = 0;
               if ( elemType == SMDSAbs_Face  && meshDim == 3 )
                 switch ( axis ) {
-                case 0: getNodesFun = & TZoneData::IFaceNodes;
-                case 1: getNodesFun = & TZoneData::JFaceNodes;
-                case 2: getNodesFun = & TZoneData::KFaceNodes;
+                case 0: getNodesFun = & TZoneData::IFaceNodes; break;
+                case 1: getNodesFun = & TZoneData::JFaceNodes; break;
+                case 2: getNodesFun = & TZoneData::KFaceNodes; break;
                 }
               else if ( elemType == SMDSAbs_Edge && meshDim == 2 )
                 switch ( axis ) {
-                case 0: getNodesFun = & TZoneData::IEdgeNodes;
-                case 1: getNodesFun = & TZoneData::JEdgeNodes;
+                case 0: getNodesFun = & TZoneData::IEdgeNodes; break;
+                case 1: getNodesFun = & TZoneData::JEdgeNodes; break;
                 }
               if ( !getNodesFun )
               {
@@ -1080,14 +1147,14 @@ Driver_Mesh::Status DriverCGNS_Read::Perform()
               PGetNodesFun getNodesFun = 0;
               if ( elemType == SMDSAbs_Face )
                 switch ( axis ) {
-                case 0: getNodesFun = & TZoneData::IFaceNodes;
-                case 1: getNodesFun = & TZoneData::JFaceNodes;
-                case 2: getNodesFun = & TZoneData::KFaceNodes;
+                case 0: getNodesFun = & TZoneData::IFaceNodes; break;
+                case 1: getNodesFun = & TZoneData::JFaceNodes; break;
+                case 2: getNodesFun = & TZoneData::KFaceNodes; break;
                 }
               else if ( elemType == SMDSAbs_Edge && meshDim == 2 )
                 switch ( axis ) {
-                case 0: getNodesFun = & TZoneData::IEdgeNodes;
-                case 1: getNodesFun = & TZoneData::JEdgeNodes;
+                case 0: getNodesFun = & TZoneData::IEdgeNodes; break;
+                case 1: getNodesFun = & TZoneData::JEdgeNodes; break;
                 }
               if ( !getNodesFun )
               {
@@ -1106,11 +1173,11 @@ Driver_Mesh::Status DriverCGNS_Read::Perform()
 
             PAddElemFun addElemFun = 0;
             switch ( meshDim ) {
-            case 1: addElemFun = & add_BAR_2;
-            case 2: addElemFun = & add_QUAD_4;
-            case 3: addElemFun = & add_HEXA_8;
+            case 1: addElemFun = & add_BAR_2;  break;
+            case 2: addElemFun = & add_QUAD_4; break;
+            case 3: addElemFun = & add_HEXA_8; break;
             }
-            int elemID = meshInfo.NbElements();
+            smIdType elemID = meshInfo.NbElements();
             const SMDS_MeshElement* elem = 0;
             for ( size_t i = 0; i < ids.size(); i += nbElemNodes )
             {
@@ -1145,12 +1212,62 @@ Driver_Mesh::Status DriverCGNS_Read::Perform()
 
       } // loop on BCs of the zone
     }
-    else
+    else addMessage( cg_get_error() );
+
+    
+    MESSAGE("  read flow solutions");
+    int nsols = 0;
+    if ( cg_nsols( _fn, cgnsBase, iZone, &nsols) == CG_OK )
     {
-      addMessage( cg_get_error() );
+       MESSAGE("  nb flow solutions: " << nsols);
+    }
+    else addMessage( cg_get_error() );
+    
+    MESSAGE("  read discrete data");
+    int nbdiscrete = 0;
+    if ( cg_ndiscrete( _fn, cgnsBase, iZone, &nbdiscrete) == CG_OK )
+    {
+        MESSAGE("  nb discrete data: " << nbdiscrete);
+        char nameDiscrete[CGNS_NAME_SIZE];
+        for (int idisc = 1; idisc <= nbdiscrete; idisc++)
+        {
+            if ( cg_discrete_read( _fn, cgnsBase, iZone, idisc, nameDiscrete) == CG_OK )
+            {
+                MESSAGE("  discrete data #"<< idisc << " name: " << nameDiscrete);
+                PointSetType_t ptset_type;
+                cgsize_t npnts;
+                if ( cg_discrete_ptset_info( _fn, cgnsBase, iZone, idisc, &ptset_type, &npnts) == CG_OK )
+                {
+                    MESSAGE("  discrete data #"<< idisc << " npnts: " << npnts);
+                }
+                else addMessage( cg_get_error() );
+            }
+            else addMessage( cg_get_error() );
+        }
     }
+    else addMessage( cg_get_error() );
+
+
+    MESSAGE("  read subregions");
+    int nbSubrg = 0;
+    if ( cg_nsubregs( _fn, cgnsBase, iZone, &nbSubrg) == CG_OK )
+    {
+       MESSAGE("  nb subregions: " << nbSubrg);
+    }
+    else addMessage( cg_get_error() );
+
+    MESSAGE("  end zone");
   } // loop on the zones of a mesh
 
+    MESSAGE("read families");
+    int nbFam = 0;
+    if ( cg_nfamilies( _fn, cgnsBase, &nbFam) == CG_OK )
+    {
+        MESSAGE("nb families: " << nbFam);
+    }
+    else addMessage( cg_get_error() );
+
+    
 
   // ------------------------------------------------------------------------
   // Make groups for multiple zones and remove free nodes at zone interfaces
@@ -1158,6 +1275,7 @@ Driver_Mesh::Status DriverCGNS_Read::Perform()
   map< string, TZoneData >::iterator nameZoneIt = zonesByName.begin();
   for ( ; nameZoneIt != zonesByName.end(); ++nameZoneIt )
   {
+    MESSAGE("nameZone: " << nameZoneIt->first);
     TZoneData& zone = nameZoneIt->second;
     if ( zone._nbElems == 0 ) continue;
     if ( zone._nbElems == meshInfo.NbElements() ) break; // there is only one non-empty zone
@@ -1186,7 +1304,7 @@ Driver_Mesh::Status DriverCGNS_Read::Perform()
 
   myMesh->Modified();
   myMesh->CompactMesh();
-
+  MESSAGE("end perform");
   return aResult;
 }