23427: [CEA 2073] No hypothesis "Viscous Layers" with Netgen 1D-2D-3D

resources/NETGENPlugin.xml

@@ -152,6 +152,7 @@
                group-id         ="1"
                priority         ="10"
                hypos            ="NETGEN_Parameters, NETGEN_SimpleParameters_3D"
+               opt-hypos        ="ViscousLayers"
                output           ="TETRA,PYRAMID"
                dim              ="3"
                support-submeshes="true">
@@ -159,6 +160,7 @@
         <algo>NETGEN_2D3D=Tetrahedron(algo=smeshBuilder.NETGEN_1D2D3D)</algo>
         <hypo>NETGEN_Parameters=Parameters()</hypo>
         <hypo>NETGEN_SimpleParameters_3D=Parameters(smeshBuilder.SIMPLE)</hypo>
+        <hypo>ViscousLayers=ViscousLayers(SetTotalThickness(),SetNumberLayers(),SetStretchFactor(),SetFaces(1),SetFaces(2),SetMethod())</hypo>
       </python-wrap>
     </algorithm>
 

src/NETGENPlugin/NETGENPlugin_Mesher.cxx

@@ -157,6 +157,7 @@ NETGENPlugin_Mesher::NETGENPlugin_Mesher (SMESH_Mesh*         mesh,
     _progressTic(1),
     _totalTime(1.0),
     _simpleHyp(NULL),
+    _viscousLayersHyp(NULL),
     _ptrToMe(NULL)
 {
   SetDefaultParameters();
@@ -336,6 +337,17 @@ void NETGENPlugin_Mesher::SetParameters(const NETGENPlugin_SimpleHypothesis_2D*
     SetDefaultParameters();
 }
 
+//================================================================================
+/*!
+ * \brief Store a Viscous Layers hypothesis
+ */
+//================================================================================
+
+void NETGENPlugin_Mesher::SetParameters(const StdMeshers_ViscousLayers* hyp )
+{
+  _viscousLayersHyp = hyp;
+}
+
 //=============================================================================
 /*!
  *  Link - a pair of integer numbers
@@ -1054,13 +1066,6 @@ bool NETGENPlugin_Mesher::FillNgMesh(netgen::OCCGeometry&           occgeom,
 
       // Find solids the geomFace bounds
       int solidID1 = 0, solidID2 = 0;
-      StdMeshers_QuadToTriaAdaptor* quadAdaptor =
-        dynamic_cast<StdMeshers_QuadToTriaAdaptor*>( proxyMesh.get() );
-      if ( quadAdaptor )
-      {
-        solidID1 = occgeom.somap.FindIndex( quadAdaptor->GetShape() );
-      }
-      else
       {
         PShapeIteratorPtr solidIt = helper.GetAncestors( geomFace, *sm->GetFather(), TopAbs_SOLID);
         while ( const TopoDS_Shape * solid = solidIt->next() )
@@ -1070,6 +1075,81 @@ bool NETGENPlugin_Mesher::FillNgMesh(netgen::OCCGeometry&           occgeom,
           else                              solidID1 = id;
         }
       }
+      if ( proxyMesh && proxyMesh->GetProxySubMesh( geomFace ))
+      {
+        // if a proxy sub-mesh contains temporary faces, then these faces
+        // should be used to mesh only one SOLID
+        bool hasTmp = false;
+        smDS = proxyMesh->GetSubMesh( geomFace );
+        SMDS_ElemIteratorPtr faces = smDS->GetElements();
+        while ( faces->more() )
+        {
+          const SMDS_MeshElement* f = faces->next();
+          if ( proxyMesh->IsTemporary( f ))
+          {
+            hasTmp = true;
+            std::vector<const SMDS_MeshNode*> fNodes( f->begin_nodes(), f->end_nodes() );
+            std::vector<const SMDS_MeshElement*> vols;
+            if ( _mesh->GetMeshDS()->GetElementsByNodes( fNodes, vols, SMDSAbs_Volume ) == 1 )
+            {
+              int geomID = vols[0]->getshapeId();
+              const TopoDS_Shape& solid =  helper.GetMeshDS()->IndexToShape( geomID );
+              if ( !solid.IsNull() )
+                solidID1 = occgeom.somap.FindIndex ( solid );
+              solidID2 = 0;
+              break;
+            }
+          }
+        }
+        // exclude faces generated by NETGEN from computation of 3D mesh
+        const int fID = occgeom.fmap.FindIndex( geomFace );
+        if ( !hasTmp ) // shrunk mesh
+        {
+          // move netgen points according to moved nodes
+          SMESH_subMeshIteratorPtr smIt = sm->getDependsOnIterator(/*includeSelf=*/true);
+          while ( smIt->more() )
+          {
+            SMESH_subMesh* sub = smIt->next();
+            if ( !sub->GetSubMeshDS() ) continue;
+            SMDS_NodeIteratorPtr nodeIt = sub->GetSubMeshDS()->GetNodes();
+            while ( nodeIt->more() )
+            {
+              const SMDS_MeshNode* n = nodeIt->next();
+              int ngID = ngNodeId( n, ngMesh, nodeNgIdMap );
+              netgen::MeshPoint& ngPoint = ngMesh.Point( ngID );
+              ngPoint(0) = n->X();
+              ngPoint(1) = n->Y();
+              ngPoint(2) = n->Z();
+            }
+          }
+          // remove faces near boundary to avoid their overlapping
+          // with shrunk faces
+          for ( int i = 1; i <= ngMesh.GetNSE(); ++i )
+          {
+            const netgen::Element2d& elem = ngMesh.SurfaceElement(i);
+            if ( elem.GetIndex() == fID )
+            {
+              for ( int iN = 0; iN < elem.GetNP(); ++iN )
+                if ( ngMesh[ elem[ iN ]].Type() != netgen::SURFACEPOINT )
+                {
+                  ngMesh.DeleteSurfaceElement( i );
+                  break;
+                }
+            }
+          }
+        }
+        //if ( hasTmp )
+        {
+          faceNgID++;
+          ngMesh.AddFaceDescriptor( netgen::FaceDescriptor( faceNgID,/*solid1=*/0,/*solid2=*/0,0 ));
+          for (int i = 1; i <= ngMesh.GetNSE(); ++i )
+          {
+            const netgen::Element2d& elem = ngMesh.SurfaceElement(i);
+            if ( elem.GetIndex() == fID )
+              const_cast< netgen::Element2d& >( elem ).SetIndex( faceNgID );
+          }
+        }
+      }
       // Add ng face descriptors of meshed faces
       faceNgID++;
       ngMesh.AddFaceDescriptor( netgen::FaceDescriptor( faceNgID, solidID1, solidID2, 0 ));
@@ -1112,8 +1192,6 @@ bool NETGENPlugin_Mesher::FillNgMesh(netgen::OCCGeometry&           occgeom,
       cout << "SMESH face " << helper.GetMeshDS()->ShapeToIndex( geomFace )
            << " internal="<<isInternalFace << endl;
 #endif
-      if ( proxyMesh )
-        smDS = proxyMesh->GetSubMesh( geomFace );
 
       SMDS_ElemIteratorPtr faces = smDS->GetElements();
       while ( faces->more() )
@@ -1430,10 +1508,15 @@ namespace
     int   ngIdCloseN; //!< ng id of closest node of the closest 2d mesh element
   };
 
-  inline double dist2(const netgen::MeshPoint& p1, const netgen::MeshPoint& p2)
+  inline double dist2( const netgen::MeshPoint& p1, const netgen::MeshPoint& p2 )
   {
     return gp_Pnt( NGPOINT_COORDS(p1)).SquareDistance( gp_Pnt( NGPOINT_COORDS(p2)));
   }
+
+  // inline double dist2(const netgen::MeshPoint& p, const SMDS_MeshNode* n )
+  // {
+  //   return gp_Pnt( NGPOINT_COORDS(p)).SquareDistance( SMESH_NodeXYZ(n));
+  // }
 }
 
 //================================================================================
@@ -2136,13 +2219,35 @@ int NETGENPlugin_Mesher::FillSMesh(const netgen::OCCGeometry&          occgeo,
   if ( quadHelper && !quadHelper->GetIsQuadratic() && quadHelper->GetTLinkNodeMap().empty() )
     quadHelper = 0;
 
+  int i, nbInitNod = initState._nbNodes;
+  if ( initState._elementsRemoved )
+  {
+    // PAL23427. Update nodeVec to track removal of netgen free points as a result
+    // of removal of faces in FillNgMesh() in the case of a shrunk sub-mesh
+    int ngID, nodeVecSize = nodeVec.size();
+    const double eps = std::numeric_limits<double>::min();
+    for ( ngID = i = 1; i < nodeVecSize; ++ngID, ++i )
+    {
+      gp_Pnt ngPnt( NGPOINT_COORDS( ngMesh.Point( ngID )));
+      gp_Pnt node ( SMESH_NodeXYZ ( nodeVec[ i ]));
+      if ( ngPnt.SquareDistance( node ) < eps )
+      {
+        nodeVec[ ngID ] = nodeVec[ i ];
+      }
+      else
+      {
+        --ngID;
+      }
+    }
+    nodeVec.resize( ngID );
+    nbInitNod = ngID - 1;
+  }
   // -------------------------------------
   // Create and insert nodes into nodeVec
   // -------------------------------------
 
   nodeVec.resize( nbNod + 1 );
-  int i, nbInitNod = initState._nbNodes;
+  for ( i = nbInitNod+1; i <= nbNod; ++i )
-  for (i = nbInitNod+1; i <= nbNod; ++i )
   {
     const netgen::MeshPoint& ngPoint = ngMesh.Point(i);
     SMDS_MeshNode* node = NULL;
@@ -2833,37 +2938,58 @@ bool NETGENPlugin_Mesher::Compute()
     // generate volume mesh
     // ---------------------
     // Fill _ngMesh with nodes and faces of computed 2D submeshes
-    if ( !err && _isVolume && ( !meshedSM[ MeshDim_2D ].empty() || mparams.quad ))
+    if ( !err && _isVolume &&
+         ( !meshedSM[ MeshDim_2D ].empty() || mparams.quad || _viscousLayersHyp ))
     {
       // load SMESH with computed segments and faces
       FillSMesh( occgeo, *_ngMesh, initState, *_mesh, nodeVec, comment, &quadHelper );
 
+      // compute prismatic boundary volumes
+      int nbQuad = _mesh->NbQuadrangles();
+      SMESH_ProxyMesh::Ptr viscousMesh;
+      if ( _viscousLayersHyp )
+      {
+        viscousMesh = _viscousLayersHyp->Compute( *_mesh, _shape );
+        if ( !viscousMesh )
+          return false;
+      }
       // compute pyramids on quadrangles
-      SMESH_ProxyMesh::Ptr proxyMesh;
+      vector<SMESH_ProxyMesh::Ptr> pyramidMeshes( occgeo.somap.Extent() );
-      if ( _mesh->NbQuadrangles() > 0 )
+      if ( nbQuad > 0 )
         for ( int iS = 1; iS <= occgeo.somap.Extent(); ++iS )
         {
-          StdMeshers_QuadToTriaAdaptor* Adaptor = new StdMeshers_QuadToTriaAdaptor;
+          StdMeshers_QuadToTriaAdaptor* adaptor = new StdMeshers_QuadToTriaAdaptor;
-          proxyMesh.reset( Adaptor );
+          pyramidMeshes[ iS-1 ].reset( adaptor );
-
+          bool ok = adaptor->Compute( *_mesh, occgeo.somap(iS), viscousMesh.get() );
-          int nbPyrams = _mesh->NbPyramids();
+          if ( !ok )
-          Adaptor->Compute( *_mesh, occgeo.somap(iS) );
+            return false;
-          if ( nbPyrams != _mesh->NbPyramids() )
+        }
+      // add proxy faces to NG mesh
+      list< SMESH_subMesh* > viscousSM;
+      for ( int iS = 1; iS <= occgeo.somap.Extent(); ++iS )
       {
         list< SMESH_subMesh* > quadFaceSM;
         for (TopExp_Explorer face(occgeo.somap(iS), TopAbs_FACE); face.More(); face.Next())
-              if ( Adaptor->GetProxySubMesh( face.Current() ))
+          if ( pyramidMeshes[iS-1] && pyramidMeshes[iS-1]->GetProxySubMesh( face.Current() ))
           {
             quadFaceSM.push_back( _mesh->GetSubMesh( face.Current() ));
             meshedSM[ MeshDim_2D ].remove( quadFaceSM.back() );
           }
-            FillNgMesh(occgeo, *_ngMesh, nodeVec, quadFaceSM, &quadHelper, proxyMesh);
+          else if ( viscousMesh && viscousMesh->GetProxySubMesh( face.Current() ))
+          {
+            viscousSM.push_back( _mesh->GetSubMesh( face.Current() ));
+            meshedSM[ MeshDim_2D ].remove( viscousSM.back() );
           }
+        if ( !quadFaceSM.empty() )
+          FillNgMesh(occgeo, *_ngMesh, nodeVec, quadFaceSM, &quadHelper, pyramidMeshes[iS-1]);
       }
+      if ( !viscousSM.empty() )
+        FillNgMesh(occgeo, *_ngMesh, nodeVec, viscousSM, &quadHelper, viscousMesh );
+
       // fill _ngMesh with faces of sub-meshes
       err = ! ( FillNgMesh(occgeo, *_ngMesh, nodeVec, meshedSM[ MeshDim_2D ], &quadHelper));
-      initState = NETGENPlugin_ngMeshInfo(_ngMesh);
+      initState = NETGENPlugin_ngMeshInfo(_ngMesh, /*checkRemovedElems=*/true);
-      //toPython( _ngMesh, "/tmp/ngPython.py");
+      // toPython( _ngMesh );
     }
     if (!err && _isVolume)
     {
@@ -3554,8 +3680,9 @@ void NETGENPlugin_Mesher::toPython( const netgen::Mesh* ngMesh )
  */
 //================================================================================
 
-NETGENPlugin_ngMeshInfo::NETGENPlugin_ngMeshInfo( netgen::Mesh* ngMesh):
+NETGENPlugin_ngMeshInfo::NETGENPlugin_ngMeshInfo( netgen::Mesh* ngMesh,
-  _copyOfLocalH(0)
+                                                  bool          checkRemovedElems):
+  _elementsRemoved( false ), _copyOfLocalH(0)
 {
   if ( ngMesh )
   {
@@ -3563,6 +3690,10 @@ NETGENPlugin_ngMeshInfo::NETGENPlugin_ngMeshInfo( netgen::Mesh* ngMesh):
     _nbSegments = ngMesh->GetNSeg();
     _nbFaces    = ngMesh->GetNSE();
     _nbVolumes  = ngMesh->GetNE();
+
+    if ( checkRemovedElems )
+      for ( int i = 1; i <= ngMesh->GetNSE() &&  !_elementsRemoved; ++i )
+        _elementsRemoved = ngMesh->SurfaceElement(i).IsDeleted();
   }
   else
   {

src/NETGENPlugin/NETGENPlugin_Mesher.hxx

@@ -46,14 +46,15 @@ namespace nglib {
 #include <vector>
 #include <set>
 
+class NETGENPlugin_Hypothesis;
+class NETGENPlugin_Internals;
+class NETGENPlugin_SimpleHypothesis_2D;
 class SMESHDS_Mesh;
 class SMESH_Comment;
 class SMESH_Mesh;
 class SMESH_MesherHelper;
+class StdMeshers_ViscousLayers;
 class TopoDS_Shape;
-class NETGENPlugin_Hypothesis;
-class NETGENPlugin_SimpleHypothesis_2D;
-class NETGENPlugin_Internals;
 namespace netgen {
   class OCCGeometry;
   class Mesh;
@@ -67,8 +68,9 @@ namespace netgen {
 struct NETGENPlugin_ngMeshInfo
 {
   int   _nbNodes, _nbSegments, _nbFaces, _nbVolumes;
+  bool  _elementsRemoved; // case where netgen can remove free nodes
   char* _copyOfLocalH;
-  NETGENPlugin_ngMeshInfo( netgen::Mesh* ngMesh=0);
+  NETGENPlugin_ngMeshInfo( netgen::Mesh* ngMesh=0, bool checkRemovedElems=false );
   void transferLocalH( netgen::Mesh* fromMesh, netgen::Mesh* toMesh );
   void restoreLocalH ( netgen::Mesh* ngMesh);
 };
@@ -120,6 +122,7 @@ class NETGENPLUGIN_EXPORT NETGENPlugin_Mesher
 
   void SetParameters(const NETGENPlugin_Hypothesis*          hyp);
   void SetParameters(const NETGENPlugin_SimpleHypothesis_2D* hyp);
+  void SetParameters(const StdMeshers_ViscousLayers*         hyp );
   void SetViscousLayers2DAssigned(bool isAssigned) { _isViscousLayers2D = isAssigned; }
   static void SetLocalSize( netgen::OCCGeometry& occgeo, netgen::Mesh& ngMesh );
 
@@ -210,6 +213,7 @@ class NETGENPLUGIN_EXPORT NETGENPlugin_Mesher
   volatile double      _totalTime;
 
   const NETGENPlugin_SimpleHypothesis_2D * _simpleHyp;
+  const StdMeshers_ViscousLayers*   _viscousLayersHyp;
 
   // a pointer to NETGENPlugin_Mesher* field of the holder, that will be
   // nullified at destruction of this

src/NETGENPlugin/NETGENPlugin_NETGEN_2D3D.cxx

@@ -32,10 +32,12 @@
 #include "NETGENPlugin_SimpleHypothesis_3D.hxx"
 #include "NETGENPlugin_Mesher.hxx"
 
+#include <SMESHDS_Mesh.hxx>
+#include <SMESH_ControlsDef.hxx>
 #include <SMESH_Gen.hxx>
 #include <SMESH_Mesh.hxx>
-#include <SMESH_ControlsDef.hxx>
+#include <StdMeshers_ViscousLayers.hxx>
-#include <SMESHDS_Mesh.hxx>
+
 #include <utilities.h>
 
 #include <list>
@@ -62,6 +64,7 @@ NETGENPlugin_NETGEN_2D3D::NETGENPlugin_NETGEN_2D3D(int hypId, int studyId,
   _shapeType = (1 << TopAbs_SHELL) | (1 << TopAbs_SOLID);// 1 bit /shape type
   _compatibleHypothesis.push_back("NETGEN_Parameters");
   _compatibleHypothesis.push_back("NETGEN_SimpleParameters_3D");
+  _compatibleHypothesis.push_back( StdMeshers_ViscousLayers::GetHypType() );
   _requireDiscreteBoundary = false;
   _onlyUnaryInput = false;
   _hypothesis = NULL;
@@ -85,35 +88,41 @@ NETGENPlugin_NETGEN_2D3D::~NETGENPlugin_NETGEN_2D3D()
  */
 //=============================================================================
 
-bool NETGENPlugin_NETGEN_2D3D::CheckHypothesis
+bool NETGENPlugin_NETGEN_2D3D::CheckHypothesis (SMESH_Mesh&         aMesh,
-                         (SMESH_Mesh& aMesh,
                                                 const TopoDS_Shape& aShape,
-                          SMESH_Hypothesis::Hypothesis_Status& aStatus)
+                                                Hypothesis_Status&  aStatus)
 {
   _hypothesis       = NULL;
+  _viscousLayersHyp = NULL;
   _mesher           = NULL;
 
-  const list<const SMESHDS_Hypothesis*>& hyps = GetUsedHypothesis(aMesh, aShape);
+  const list<const SMESHDS_Hypothesis*>& hyps = GetUsedHypothesis(aMesh, aShape, /*noAux=*/false);
-  int nbHyp = hyps.size();
+  if ( hyps.empty() )
-  if (!nbHyp)
   {
     aStatus = SMESH_Hypothesis::HYP_OK;
     return true;  // can work with no hypothesis
   }
 
-  const SMESHDS_Hypothesis* theHyp = hyps.front(); // use only the first hypothesis
+  list<const SMESHDS_Hypothesis*>::const_iterator h = hyps.begin();
+  for ( ; h != hyps.end(); ++h )
+  {
+    const SMESHDS_Hypothesis* aHyp = *h;
+    std::string            hypName = aHyp->GetName();
 
-  string hypName = theHyp->GetName();
+    if ( std::find( _compatibleHypothesis.begin(), _compatibleHypothesis.end(),
-
-  if ( find( _compatibleHypothesis.begin(), _compatibleHypothesis.end(),
                     hypName ) != _compatibleHypothesis.end() )
     {
-    _hypothesis = theHyp;
+      if ( hypName == StdMeshers_ViscousLayers::GetHypType() )
+        _viscousLayersHyp = dynamic_cast<const StdMeshers_ViscousLayers*>( aHyp );
+      else
+        _hypothesis = aHyp;
       aStatus = SMESH_Hypothesis::HYP_OK;
     }
     else
     {
       aStatus = SMESH_Hypothesis::HYP_INCOMPATIBLE;
+      break;
+    }
   }
 
   return aStatus == SMESH_Hypothesis::HYP_OK;
@@ -133,6 +142,7 @@ bool NETGENPlugin_NETGEN_2D3D::Compute(SMESH_Mesh&         aMesh,
   NETGENPlugin_Mesher mesher(&aMesh, aShape, true);
   mesher.SetParameters(dynamic_cast<const NETGENPlugin_Hypothesis*>(_hypothesis));
   mesher.SetParameters(dynamic_cast<const NETGENPlugin_SimpleHypothesis_3D*>(_hypothesis));
+  mesher.SetParameters(_viscousLayersHyp);
   mesher.SetSelfPointer( &_mesher );
   return mesher.Compute();
 }

src/NETGENPlugin/NETGENPlugin_NETGEN_2D3D.hxx

@@ -35,6 +35,7 @@
 #include <SMESH_Algo.hxx>
 
 class NETGENPlugin_Mesher;
+class StdMeshers_ViscousLayers;
 
 class NETGENPLUGIN_EXPORT NETGENPlugin_NETGEN_2D3D: public SMESH_3D_Algo
 {
@@ -44,7 +45,7 @@ public:
 
   virtual bool CheckHypothesis(SMESH_Mesh&         aMesh,
                                const TopoDS_Shape& aShape,
-                               SMESH_Hypothesis::Hypothesis_Status& aStatus);
+                               Hypothesis_Status&  aStatus);
 
   virtual bool Compute(SMESH_Mesh&         aMesh,
                        const TopoDS_Shape& aShape);
@@ -60,6 +61,7 @@ public:
 
 protected:
   const SMESHDS_Hypothesis*       _hypothesis;
+  const StdMeshers_ViscousLayers* _viscousLayersHyp;
   NETGENPlugin_Mesher *           _mesher;
 };