Salome HOME
Merge branch 'occ/shaper2smesh'
[modules/smesh.git] / src / DriverCGNS / DriverCGNS_Read.cxx
index 492c1eb62a75ca2162673a512d8f5e9c22ed8cf2..48b1a3a286bab86aeec86801e68f9b8da42e3f3f 100644 (file)
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2013  CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2019  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
@@ -6,7 +6,7 @@
 // This library is free software; you can redistribute it and/or
 // modify it under the terms of the GNU Lesser General Public
 // License as published by the Free Software Foundation; either
-// version 2.1 of the License.
+// version 2.1 of the License, or (at your option) any later version.
 //
 // This library is distributed in the hope that it will be useful,
 // but WITHOUT ANY WARRANTY; without even the implied warranty of
@@ -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"
 
@@ -29,6 +31,7 @@
 #include "SMESHDS_Group.hxx"
 #include "SMESHDS_Mesh.hxx"
 #include "SMESH_Comment.hxx"
+#include "SMESH_TypeDefs.hxx"
 
 #include <gp_XYZ.hxx>
 
@@ -36,6 +39,7 @@
 
 #include <map>
 
+
 #if CGNS_VERSION < 3100
 # define cgsize_t int
 #endif
@@ -76,7 +80,9 @@ namespace
     }
     bool IsStructured() const { return ( _type == CGNS_ENUMV( Structured )); }
     int IndexSize() const { return IsStructured() ? _meshDim : 1; }
-    string ReadZonesConnection(int file, int base, const map< string, TZoneData >& zonesByName);
+    string ReadZonesConnection(int file, int base,
+                               const map< string, TZoneData >& zonesByName,
+                               SMESHDS_Mesh*                   mesh);
     void ReplaceNodes( cgsize_t* ids, int nbIds, int idShift = 0 ) const;
 
     // Methods for a structured zone
@@ -205,6 +211,37 @@ namespace
     //gp_XYZ End() const { return gp_XYZ( _end[0]-1, _end[1]-1, _end[2]-1 ); }
   };
 
+  //================================================================================
+  /*!
+   * \brief Checks if the two arrays of node IDs describe nodes with equal coordinates
+   */
+  //================================================================================
+
+  bool isEqualNodes( const int* nIds1, const int* nIds2, int nbNodes, SMESHDS_Mesh* mesh )
+  {
+    if ( nbNodes > 0 )
+    {
+      SMESH_TNodeXYZ nn1[2], nn2[2];
+      nn1[0] = mesh->FindNode( nIds1[0] );
+      nn2[0] = mesh->FindNode( nIds2[0] );
+      if ( !nn1[0]._node || !nn2[0]._node )
+        return false;
+      double dist1 = ( nn1[0] - nn2[0] ).Modulus();
+      double dist2 = 0, tol = 1e-7;
+      if ( nbNodes > 1 )
+      {
+        nn1[1] = mesh->FindNode( nIds1[1] );
+        nn2[1] = mesh->FindNode( nIds2[1] );
+        if ( !nn1[1]._node || !nn2[1]._node )
+          return false;
+        dist2 = ( nn1[1] - nn2[1] ).Modulus();
+        tol   = 1e-5 * ( nn1[0] - nn1[1] ).Modulus();
+      }
+      return ( dist1 < tol && dist2 < tol );
+    }
+    return false;
+  }
+
   //================================================================================
   /*!
    * \brief Reads zone interface connectivity
@@ -220,7 +257,8 @@ namespace
 
   string TZoneData::ReadZonesConnection( int                             file,
                                          int                             base,
-                                         const map< string, TZoneData >& zonesByName)
+                                         const map< string, TZoneData >& zonesByName,
+                                         SMESHDS_Mesh*                   mesh)
   {
     string error;
 
@@ -277,6 +315,26 @@ namespace
                   continue; // this interface already read
               }
             }
+            // check if range and donorRange describe the same nodes
+            {
+              int ids1[2], ids2[2], nbN = 0;
+              TPointRangeIterator rangeIt1bis( range, _meshDim );
+              index1 = rangeIt1bis.Next();
+              index2 = T * ( index1 - begin1 ) + begin2;
+              ids1[0] = NodeID( index1 );
+              ids2[0] = zone2.NodeID( index2 );
+              ++nbN;
+              if ( rangeIt1bis.More() )
+              {
+                index1 = rangeIt1bis.Next();
+                index2 = T * ( index1 - begin1 ) + begin2;
+                ids1[1] = NodeID( index1 );
+                ids2[1] = zone2.NodeID( index2 );
+                ++nbN;
+              }
+              if ( !isEqualNodes( &ids1[0], &ids2[0], nbN, mesh ))
+                continue;
+            }
             while ( rangeIt1.More() )
             {
               index1 = rangeIt1.Next();
@@ -338,7 +396,7 @@ namespace
           {
             for ( int isThisZone = 0; isThisZone < 2; ++isThisZone )
             {
-              const TZoneData&            zone = isThisZone ? *this : zone2;
+              const TZoneData&           zone = isThisZone ? *this : zone2;
               CGNS_ENUMT(PointSetType_t) type = isThisZone ? ptype : donorPtype;
               vector< cgsize_t >&      points = isThisZone ? ids : donorIds;
               if ( type == CGNS_ENUMV( PointRange ))
@@ -361,8 +419,10 @@ namespace
                   points[i] += zone._nodeIdShift;
               }
             }
-            for ( size_t i = 0; i < ids.size() && i < donorIds.size(); ++i )
-              _nodeReplacementMap.insert( make_pair( ids[i], donorIds[i] ));
+            size_t nbN = std::min( ids.size(), donorIds.size());
+            if ( isEqualNodes( &ids[0], &donorIds[0], nbN, mesh ))
+              for ( size_t i = 0; i < nbN; ++i )
+                _nodeReplacementMap.insert( make_pair( ids[i], donorIds[i] ));
           }
           else
           {
@@ -394,7 +454,7 @@ namespace
     if ( !_nodeReplacementMap.empty() )
     {
       map< int, int >::const_iterator it, end = _nodeReplacementMap.end();
-      for ( size_t i = 0; i < nbIds; ++i )
+      for ( int i = 0; i < nbIds; ++i )
         if (( it = _nodeReplacementMap.find( ids[i] + idShift)) != end )
           ids[i] = it->second;
         else
@@ -402,7 +462,7 @@ namespace
     }
     else if ( idShift )
     {
-      for ( size_t i = 0; i < nbIds; ++i )
+      for ( int i = 0; i < nbIds; ++i )
         ids[i] += idShift;
     }
   }
@@ -576,6 +636,7 @@ namespace
 
 Driver_Mesh::Status DriverCGNS_Read::Perform()
 {
+  MESSAGE("DriverCGNS_Read::Perform");
   myErrorMessages.clear();
 
   Status aResult;
@@ -585,6 +646,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;
@@ -598,6 +660,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;
@@ -606,6 +670,7 @@ Driver_Mesh::Status DriverCGNS_Read::Perform()
 
   if ( nbZones < 1 )
     return addMessage( SMESH_Comment("Empty mesh: '") << meshName << "'");
+  MESSAGE("nbZones: " << nbZones);
 
   // read the domains (zones)
   // ------------------------
@@ -629,6 +694,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) {
@@ -639,7 +705,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");
@@ -655,7 +724,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;
@@ -667,6 +736,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];
@@ -698,6 +768,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 );
@@ -709,14 +780,15 @@ Driver_Mesh::Status DriverCGNS_Read::Perform()
     }
 
     // Read connectivity between zones. Nodes of the zone interface will be
-    // replaced withing the zones read later
-    string err = zone.ReadZonesConnection( _fn, cgnsBase, zonesByName );
+    // replaced within the zones read later
+    string err = zone.ReadZonesConnection( _fn, cgnsBase, zonesByName, myMesh );
     if ( !err.empty() )
       addMessage( err );
 
     // --------------
     // Read elements
     // --------------
+    MESSAGE("  read elements");
     if ( zone.IsStructured())
     {
       int nbI = zone._sizeX - 1, nbJ = zone._sizeY - 1, nbK = zone._sizeZ - 1;
@@ -769,12 +841,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 )
@@ -790,6 +864,7 @@ Driver_Mesh::Status DriverCGNS_Read::Perform()
         }
         // store elements
 
+        MESSAGE("   store elements");
         int pos = 0, cgnsNbNodes = 0, elemID = start + zone._elemIdShift;
         cg_npe( elemType, &cgnsNbNodes ); // get nb nodes by element type
         curAddElemFun = addElemFuns[ elemType ];
@@ -827,7 +902,7 @@ Driver_Mesh::Status DriverCGNS_Read::Perform()
                 {
                   const bool reverse = ( elemData[ pos-1 ] < 0 );
                   const int    iQuad = face->IsQuadratic() ? 1 : 0;
-                  SMDS_ElemIteratorPtr nIter = face->interlacedNodesElemIterator();
+                  SMDS_NodeIteratorPtr nIter = face->interlacedNodesIterator();
                   faceNodes.assign( SMDS_MeshElement::iterator( nIter ),
                                     SMDS_MeshElement::iterator());
                   if ( iQuad && reverse )
@@ -887,6 +962,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 )
     {
@@ -895,14 +972,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 ||
@@ -965,7 +1045,7 @@ Driver_Mesh::Status DriverCGNS_Read::Perform()
           if ( zone.IsStructured() )
           {
             int axis = 0; // axis perpendiculaire to which boundary elements are oriented
-            if ( ids.size() >= meshDim * 2 )
+            if ( (int) ids.size() >= meshDim * 2 )
             {
               for ( ; axis < meshDim; ++axis )
                 if ( ids[axis] - ids[axis+meshDim] == 0 )
@@ -1070,7 +1150,7 @@ Driver_Mesh::Status DriverCGNS_Read::Perform()
 
             if ( psType == CGNS_ENUMV( PointRange ) && ids.size() == 2 )
             {
-              for ( size_t i = ids[0]; i <= ids[1]; ++i )
+              for ( cgsize_t i = ids[0]; i <= ids[1]; ++i )
                 if ( const SMDS_MeshElement* e = myMesh->FindElement( i ))
                   groupDS.Add( e );
             }
@@ -1088,12 +1168,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
@@ -1101,6 +1231,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
@@ -1127,6 +1258,9 @@ Driver_Mesh::Status DriverCGNS_Read::Perform()
 
   aResult = myErrorMessages.empty() ? DRS_OK : DRS_WARN_SKIP_ELEM;
 
+  myMesh->Modified();
+  myMesh->CompactMesh();
+  MESSAGE("end perform");
   return aResult;
 }