0021543: EDF 1978 SMESH: Viscous layer for 2D meshes

Add StdMeshers_ViscousLayers2D hypothesis

resources/NETGENPlugin.xml

@@ -79,7 +79,7 @@
                label-id="Netgen 2D"
                icon-id="mesh_algo_netgen_2d.png"
                hypos="LengthFromEdges,MaxElementArea,NETGEN_Parameters_2D_ONLY"
-               opt-hypos="QuadranglePreference"
+               opt-hypos="QuadranglePreference,ViscousLayers2D"
                input="EDGE"
                output="TRIA,QUAD"
                dim="2">
@@ -96,6 +96,7 @@
                label-id="Netgen 1D-2D"
                icon-id="mesh_algo_netgen_2d.png"
                hypos="NETGEN_Parameters_2D, NETGEN_SimpleParameters_2D"
+               opt-hypos="ViscousLayers2D"
                output="TRIA,QUAD"
                dim="2"
                support-submeshes="true">

src/NETGENPlugin/NETGENPlugin_Mesher.cxx

@@ -44,6 +44,7 @@
 #include <SMESH_MesherHelper.hxx>
 #include <SMESH_subMesh.hxx>
 #include <StdMeshers_QuadToTriaAdaptor.hxx>
+#include <StdMeshers_ViscousLayers2D.hxx>
 
 #include <utilities.h>
 
@@ -96,7 +97,7 @@ using namespace std;
 // dump elements added to ng mesh
 //#define DUMP_SEGMENTS
 //#define DUMP_TRIANGLES
-//#define DUMP_TRIANGLES_SCRIPT "/tmp/trias.py" //!< debug addIntVerticesInSolids()
+//#define DUMP_TRIANGLES_SCRIPT "/tmp/trias.py" //!< debug AddIntVerticesInSolids()
 
 TopTools_IndexedMapOfShape ShapesWithLocalSize;
 std::map<int,double> VertexId2LocalSize;
@@ -109,17 +110,18 @@ std::map<int,double> FaceId2LocalSize;
  */
 //=============================================================================
 
-NETGENPlugin_Mesher::NETGENPlugin_Mesher (SMESH_Mesh* mesh,
+NETGENPlugin_Mesher::NETGENPlugin_Mesher (SMESH_Mesh*         mesh,
                                           const TopoDS_Shape& aShape,
-                                          const bool isVolume)
+                                          const bool          isVolume)
   : _mesh    (mesh),
     _shape   (aShape),
     _isVolume(isVolume),
     _optimize(true),
     _fineness(NETGENPlugin_Hypothesis::GetDefaultFineness()),
+    _isViscousLayers2D(false),
     _simpleHyp(NULL)
 {
-  defaultParameters();
+  SetDefaultParameters();
   ShapesWithLocalSize.Clear();
   VertexId2LocalSize.clear();
   EdgeId2LocalSize.clear();
@@ -132,7 +134,7 @@ NETGENPlugin_Mesher::NETGENPlugin_Mesher (SMESH_Mesh* mesh,
  */
 //================================================================================
 
-void NETGENPlugin_Mesher::defaultParameters()
+void NETGENPlugin_Mesher::SetDefaultParameters()
 {
   netgen::MeshingParameters& mparams = netgen::mparam;
   // maximal mesh edge size
@@ -253,7 +255,7 @@ void NETGENPlugin_Mesher::SetParameters(const NETGENPlugin_SimpleHypothesis_2D*
 {
   _simpleHyp = hyp;
   if ( _simpleHyp )
-    defaultParameters();
+    SetDefaultParameters();
 }
 
 //=============================================================================
@@ -607,7 +609,7 @@ void NETGENPlugin_Mesher::RestrictLocalSize(netgen::Mesh& ngMesh, const gp_XYZ&
  */
 //================================================================================
 
-bool NETGENPlugin_Mesher::fillNgMesh(netgen::OCCGeometry&           occgeom,
+bool NETGENPlugin_Mesher::FillNgMesh(netgen::OCCGeometry&           occgeom,
                                      netgen::Mesh&                  ngMesh,
                                      vector<const SMDS_MeshNode*>&  nodeVec,
                                      const list< SMESH_subMesh* > & meshedSM,
@@ -941,7 +943,7 @@ bool NETGENPlugin_Mesher::fillNgMesh(netgen::OCCGeometry&           occgeom,
  */
 //================================================================================
 
-void NETGENPlugin_Mesher::fixIntFaces(const netgen::OCCGeometry& occgeom,
+void NETGENPlugin_Mesher::FixIntFaces(const netgen::OCCGeometry& occgeom,
                                       netgen::Mesh&              ngMesh,
                                       NETGENPlugin_Internals&    internalShapes)
 {
@@ -1048,7 +1050,7 @@ namespace
  */
 //================================================================================
 
-void NETGENPlugin_Mesher::addIntVerticesInFaces(const netgen::OCCGeometry&     occgeom,
+void NETGENPlugin_Mesher::AddIntVerticesInFaces(const netgen::OCCGeometry&     occgeom,
                                                 netgen::Mesh&                  ngMesh,
                                                 vector<const SMDS_MeshNode*>&  nodeVec,
                                                 NETGENPlugin_Internals&        internalShapes)
@@ -1234,7 +1236,7 @@ void NETGENPlugin_Mesher::addIntVerticesInFaces(const netgen::OCCGeometry&     o
  */
 //================================================================================
 
-void NETGENPlugin_Mesher::addIntVerticesInSolids(const netgen::OCCGeometry&     occgeom,
+void NETGENPlugin_Mesher::AddIntVerticesInSolids(const netgen::OCCGeometry&     occgeom,
                                                  netgen::Mesh&                  ngMesh,
                                                  vector<const SMDS_MeshNode*>&  nodeVec,
                                                  NETGENPlugin_Internals&        internalShapes)
@@ -1451,6 +1453,237 @@ void NETGENPlugin_Mesher::addIntVerticesInSolids(const netgen::OCCGeometry&
   } // loop on solids with internal vertices
 }
 
+//================================================================================
+/*!
+ * \brief Fill netgen mesh with segments of a FACE
+ *  \param ngMesh - netgen mesh
+ *  \param geom - container of OCCT geometry to mesh
+ *  \param wires - data of nodes on FACE boundary
+ *  \param helper - mesher helper holding the FACE
+ *  \param nodeVec - vector of nodes in which node index == netgen ID
+ *  \retval SMESH_ComputeErrorPtr - error description 
+ */
+//================================================================================
+
+SMESH_ComputeErrorPtr
+NETGENPlugin_Mesher::AddSegmentsToMesh(netgen::Mesh&                    ngMesh,
+                                       netgen::OCCGeometry&             geom,
+                                       const TSideVector&               wires,
+                                       SMESH_MesherHelper&              helper,
+                                       vector< const SMDS_MeshNode* > & nodeVec)
+{
+  // ----------------------------
+  // Check wires and count nodes
+  // ----------------------------
+  int nbNodes = 0;
+  for ( int iW = 0; iW < wires.size(); ++iW )
+  {
+    StdMeshers_FaceSidePtr wire = wires[ iW ];
+    if ( wire->MissVertexNode() )
+    {
+      // Commented for issue 0020960. It worked for the case, let's wait for case where it doesn't.
+      // It seems that there is no reason for this limitation
+//       return TError
+//         (new SMESH_ComputeError(COMPERR_BAD_INPUT_MESH, "Missing nodes on vertices"));
+    }
+    const vector<UVPtStruct>& uvPtVec = wire->GetUVPtStruct();
+    if ( uvPtVec.size() != wire->NbPoints() )
+      return SMESH_ComputeError::New(COMPERR_BAD_INPUT_MESH,
+                                     SMESH_Comment("Unexpected nb of points on wire ") << iW
+                                     << ": " << uvPtVec.size()<<" != "<<wire->NbPoints());
+    nbNodes += wire->NbPoints();
+  }
+  nodeVec.reserve( nodeVec.size() + nbNodes + 1 );
+  if ( nodeVec.empty() )
+    nodeVec.push_back( 0 );
+
+  // -----------------
+  // Fill netgen mesh
+  // -----------------
+
+  const bool wasNgMeshEmpty = ( ngMesh.GetNP() < 1 ); /* true => this method is called by
+                                                         NETGENPlugin_NETGEN_2D_ONLY */
+
+  // map for nodes on vertices since they can be shared between wires
+  // ( issue 0020676, face_int_box.brep) and nodes built by NETGEN
+  map<const SMDS_MeshNode*, int > node2ngID;
+  if ( !wasNgMeshEmpty ) // fill node2ngID with nodes built by NETGEN
+  {
+    set< int > subIDs; // ids of sub-shapes of the FACE
+    for ( int iW = 0; iW < wires.size(); ++iW )
+    {
+      StdMeshers_FaceSidePtr wire = wires[ iW ];
+      for ( int iE = 0, nbE = wire->NbEdges(); iE < nbE; ++iE )
+      {
+        subIDs.insert( wire->EdgeID( iE ));
+        subIDs.insert( helper.GetMeshDS()->ShapeToIndex( wire->FirstVertex( iE )));
+      }
+    }
+    for ( size_t ngID = 1; ngID < nodeVec.size(); ++ngID )
+      if ( subIDs.count( nodeVec[ngID]->getshapeId() ))
+        node2ngID.insert( make_pair( nodeVec[ngID], ngID ));
+  }
+
+  const int solidID = 0, faceID = geom.fmap.FindIndex( helper.GetSubShape() );
+  if ( ngMesh.GetNFD() < 1 )
+    ngMesh.AddFaceDescriptor (netgen::FaceDescriptor(faceID, solidID, solidID, 0));
+
+  for ( int iW = 0; iW < wires.size(); ++iW )
+  {
+    StdMeshers_FaceSidePtr wire = wires[ iW ];
+    const vector<UVPtStruct>& uvPtVec = wire->GetUVPtStruct();
+    const int nbSegments = wire->NbPoints() - 1;
+
+    // assure the 1st node to be in node2ngID, which is needed to correctly
+    // "close chain of segments" (see below) in case if the 1st node is not
+    // onVertex because it is on a Viscous layer
+    node2ngID.insert( make_pair( uvPtVec[ 0 ].node, ngMesh.GetNP() + 1 ));
+
+    // compute length of every segment
+    vector<double> segLen( nbSegments );
+    for ( int i = 0; i < nbSegments; ++i )
+      segLen[i] = SMESH_TNodeXYZ( uvPtVec[ i ].node ).Distance( uvPtVec[ i+1 ].node );
+
+    int edgeID = 1, posID = -2;
+    bool isInternalWire = false;
+    double vertexNormPar = 0;
+    for ( int i = 0; i < nbSegments; ++i ) // loop on segments
+    {
+      // Add the first point of a segment
+
+      const SMDS_MeshNode * n = uvPtVec[ i ].node;
+      const int posShapeID = n->getshapeId();
+      bool onVertex = ( n->GetPosition()->GetTypeOfPosition() == SMDS_TOP_VERTEX );
+      bool onEdge   = ( n->GetPosition()->GetTypeOfPosition() == SMDS_TOP_EDGE   );
+
+      // skip nodes on degenerated edges
+      if ( helper.IsDegenShape( posShapeID ) &&
+           helper.IsDegenShape( uvPtVec[ i+1 ].node->getshapeId() ))
+        continue;
+
+      int ngID1 = ngMesh.GetNP() + 1, ngID2 = ngID1+1;
+      if ( onVertex || ( !wasNgMeshEmpty && onEdge ))
+        ngID1 = node2ngID.insert( make_pair( n, ngID1 )).first->second;
+      if ( ngID1 > ngMesh.GetNP() )
+      {
+        netgen::MeshPoint mp( netgen::Point<3> (n->X(), n->Y(), n->Z()) );
+        ngMesh.AddPoint ( mp, 1, netgen::EDGEPOINT );
+        nodeVec.push_back( n );
+      }
+      else // n is in ngMesh already, and ngID2 in prev segment is wrong
+      {
+        ngID2 = ngMesh.GetNP() + 1;
+        if ( i > 0 ) // prev segment belongs to same wire
+        {
+          netgen::Segment& prevSeg = ngMesh.LineSegment( ngMesh.GetNSeg() );
+          prevSeg[1] = ngID1;
+        }
+      }
+
+      // Add the segment
+
+      netgen::Segment seg;
+
+      seg[0]     = ngID1;                // ng node id
+      seg[1]     = ngID2;                // ng node id
+      seg.edgenr = ngMesh.GetNSeg() + 1; // ng segment id
+      seg.si     = faceID;               // = geom.fmap.FindIndex (face);
+
+      for ( int iEnd = 0; iEnd < 2; ++iEnd)
+      {
+        const UVPtStruct& pnt = uvPtVec[ i + iEnd ];
+
+        seg.epgeominfo[ iEnd ].dist = pnt.param; // param on curve
+        seg.epgeominfo[ iEnd ].u    = pnt.u;
+        seg.epgeominfo[ iEnd ].v    = pnt.v;
+
+        // find out edge id and node parameter on edge
+        onVertex = ( pnt.normParam + 1e-10 > vertexNormPar );
+        if ( onVertex || posShapeID != posID )
+        {
+          // get edge id
+          double normParam = pnt.normParam;
+          if ( onVertex )
+            normParam = 0.5 * ( uvPtVec[ i ].normParam + uvPtVec[ i+1 ].normParam );
+          int edgeIndexInWire = wire->EdgeIndex( normParam );
+          vertexNormPar = wire->LastParameter( edgeIndexInWire );
+          const TopoDS_Edge& edge = wire->Edge( edgeIndexInWire );
+          edgeID = geom.emap.FindIndex( edge );
+          posID  = posShapeID;
+          isInternalWire = ( edge.Orientation() == TopAbs_INTERNAL );
+          // if ( onVertex ) // param on curve is different on each of two edges
+          //   seg.epgeominfo[ iEnd ].dist = helper.GetNodeU( edge, pnt.node );
+        }
+        seg.epgeominfo[ iEnd ].edgenr = edgeID; //  = geom.emap.FindIndex(edge);
+      }
+
+      ngMesh.AddSegment (seg);
+      {
+        // restrict size of elements near the segment
+        SMESH_TNodeXYZ np1( n ), np2( uvPtVec[ i+1 ].node );
+        // get an average size of adjacent segments to avoid sharp change of
+        // element size (regression on issue 0020452, note 0010898)
+        int iPrev = SMESH_MesherHelper::WrapIndex( i-1, nbSegments );
+        int iNext = SMESH_MesherHelper::WrapIndex( i+1, nbSegments );
+        double avgH = ( segLen[ iPrev ] + segLen[ i ] + segLen[ iNext ]) / 3;
+
+        RestrictLocalSize( ngMesh, 0.5*(np1+np2), avgH );
+      }
+#ifdef DUMP_SEGMENTS
+        cout << "Segment: " << seg.edgenr << endl
+           << "\tp1: " << seg[0] << endl
+           << "\tp2: " << seg[1] << endl
+           << "\tp0 param: " << seg.epgeominfo[ 0 ].dist << endl
+           << "\tp0 uv: " << seg.epgeominfo[ 0 ].u <<", "<< seg.epgeominfo[ 0 ].v << endl
+           << "\tp0 edge: " << seg.epgeominfo[ 0 ].edgenr << endl
+           << "\tp1 param: " << seg.epgeominfo[ 1 ].dist << endl
+           << "\tp1 uv: " << seg.epgeominfo[ 1 ].u <<", "<< seg.epgeominfo[ 1 ].v << endl
+           << "\tp1 edge: " << seg.epgeominfo[ 1 ].edgenr << endl;
+#endif
+      if ( isInternalWire )
+      {
+        swap (seg[0], seg[1]);
+        swap( seg.epgeominfo[0], seg.epgeominfo[1] );
+        seg.edgenr = ngMesh.GetNSeg() + 1; // segment id
+        ngMesh.AddSegment (seg);
+#ifdef DUMP_SEGMENTS
+        cout << "Segment: " << seg.edgenr << endl << "\tis REVRESE of the previous one" << endl;
+#endif
+      }
+    } // loop on segments on a wire
+
+    // close chain of segments
+    if ( nbSegments > 0 )
+    {
+      netgen::Segment& lastSeg = ngMesh.LineSegment( ngMesh.GetNSeg() - int( isInternalWire));
+      const SMDS_MeshNode * lastNode = uvPtVec.back().node;
+      lastSeg[1] = node2ngID.insert( make_pair( lastNode, lastSeg[1] )).first->second;
+      if ( lastSeg[1] > ngMesh.GetNP() )
+      {
+        netgen::MeshPoint mp( netgen::Point<3> (lastNode->X(), lastNode->Y(), lastNode->Z()) );
+        ngMesh.AddPoint ( mp, 1, netgen::EDGEPOINT );
+        nodeVec.push_back( lastNode );
+      }
+      if ( isInternalWire )
+      {
+        netgen::Segment& realLastSeg = ngMesh.LineSegment( ngMesh.GetNSeg() );
+        realLastSeg[0] = lastSeg[1];
+      }
+    }
+
+  } // loop on WIREs of a FACE
+
+  // add a segment instead of an internal vertex
+  if ( wasNgMeshEmpty )
+  {
+    NETGENPlugin_Internals intShapes( *helper.GetMesh(), helper.GetSubShape(), /*is3D=*/false );
+    AddIntVerticesInFaces( geom, ngMesh, nodeVec, intShapes );
+  }
+  ngMesh.CalcSurfacesOfNode();
+
+  return TError();
+}
+
 //================================================================================
 /*!
  * \brief Fill SMESH mesh according to contents of netgen mesh
@@ -1822,11 +2055,12 @@ bool NETGENPlugin_Mesher::Compute()
           " elements per radius = " << mparams.curvaturesafety << "\n"
           " second order = " << mparams.secondorder << "\n"
           " quad allowed = " << mparams.quad);
-  cout << " quad allowed = " << mparams.quad<<endl;
+  //cout << " quad allowed = " << mparams.quad<<endl;
 
   SMESH_ComputeErrorPtr error = SMESH_ComputeError::New();
 
-
+  static string debugFile = "/tmp/ngMesh.py"; /* to call toPython( ngMesh, debugFile )
+                                                 while debugging netgen */
   // -------------------------
   // Prepare OCC geometry
   // -------------------------
@@ -2002,9 +2236,8 @@ bool NETGENPlugin_Mesher::Compute()
     // Fill ngMesh with nodes and segments of computed submeshes
     if ( !err )
     {
-      _faceDescriptors.clear();
+      err = ! ( FillNgMesh(occgeo, *ngMesh, nodeVec, meshedSM[ MeshDim_0D ]) &&
-      err = ! ( fillNgMesh(occgeo, *ngMesh, nodeVec, meshedSM[ MeshDim_0D ]) &&
+                FillNgMesh(occgeo, *ngMesh, nodeVec, meshedSM[ MeshDim_1D ]));
-                fillNgMesh(occgeo, *ngMesh, nodeVec, meshedSM[ MeshDim_1D ]));
     }
     initState = NETGENPlugin_ngMeshInfo(ngMesh);
 
@@ -2072,10 +2305,44 @@ bool NETGENPlugin_Mesher::Compute()
       if ( internals.hasInternalVertexInFace() )
       {
         // store computed segments in SMESH in order not to create SMESH
-        // edges for ng segments added by addIntVerticesInFaces()
+        // edges for ng segments added by AddIntVerticesInFaces()
         FillSMesh( occgeo, *ngMesh, initState, *_mesh, nodeVec, comment );
         // add segments to faces with internal vertices
-        addIntVerticesInFaces( occgeo, *ngMesh, nodeVec, internals );
+        AddIntVerticesInFaces( occgeo, *ngMesh, nodeVec, internals );
+        initState = NETGENPlugin_ngMeshInfo(ngMesh);
+      }
+
+      // Build viscous layers
+      if ( _isViscousLayers2D )
+      {
+        if ( !internals.hasInternalVertexInFace() ) {
+          FillSMesh( occgeo, *ngMesh, initState, *_mesh, nodeVec, comment );
+          initState = NETGENPlugin_ngMeshInfo(ngMesh);
+        }
+        SMESH_ProxyMesh::Ptr viscousMesh;
+        SMESH_MesherHelper   helper( *_mesh );
+        for ( int faceID = 1; faceID <= occgeo.fmap.Extent(); ++faceID )
+        {
+          const TopoDS_Face& F = TopoDS::Face( occgeo.fmap( faceID ));
+          viscousMesh = StdMeshers_ViscousLayers2D::Compute( *_mesh, F );
+          if ( !viscousMesh )
+            return false;
+          // exclude from computation ng segments built on EDGEs of F
+          for (int i = 1; i <= ngMesh->GetNSeg(); i++)
+          {
+            netgen::Segment & seg = ngMesh->LineSegment(i);
+            if (seg.si == faceID)
+              seg.si = 0;
+          }
+          // add new segments to ngMesh instead of excluded ones
+          helper.SetSubShape( F );
+          TSideVector wires =
+            StdMeshers_FaceSide::GetFaceWires( F, *_mesh, /*skipMediumNodes=*/true,
+                                               error, viscousMesh );
+          error = AddSegmentsToMesh( *ngMesh, occgeo, wires, helper, nodeVec );
+
+          if ( !error ) error = SMESH_ComputeError::New();
+        }
         initState = NETGENPlugin_ngMeshInfo(ngMesh);
       }
 
@@ -2131,11 +2398,11 @@ bool NETGENPlugin_Mesher::Compute()
                 quadFaceSM.push_back( _mesh->GetSubMesh( face.Current() ));
                 meshedSM[ MeshDim_2D ].remove( quadFaceSM.back() );
               }
-            fillNgMesh(occgeo, *ngMesh, nodeVec, quadFaceSM, proxyMesh);
+            FillNgMesh(occgeo, *ngMesh, nodeVec, quadFaceSM, proxyMesh);
           }
         }
       // fill ngMesh with faces of sub-meshes
-      err = ! ( fillNgMesh(occgeo, *ngMesh, nodeVec, meshedSM[ MeshDim_2D ]));
+      err = ! ( FillNgMesh(occgeo, *ngMesh, nodeVec, meshedSM[ MeshDim_2D ]));
       initState = NETGENPlugin_ngMeshInfo(ngMesh);
       //toPython( ngMesh, "/tmp/ngPython.py");
     }
@@ -2165,9 +2432,9 @@ bool NETGENPlugin_Mesher::Compute()
         // faces for ng faces added here
         FillSMesh( occgeo, *ngMesh, initState, *_mesh, nodeVec, comment );
         // add ng faces to solids with internal vertices
-        addIntVerticesInSolids( occgeo, *ngMesh, nodeVec, internals );
+        AddIntVerticesInSolids( occgeo, *ngMesh, nodeVec, internals );
         // duplicate mesh faces on internal faces
-        fixIntFaces( occgeo, *ngMesh, internals );
+        FixIntFaces( occgeo, *ngMesh, internals );
         initState = NETGENPlugin_ngMeshInfo(ngMesh);
       }
       // Let netgen compute 3D mesh
@@ -2258,7 +2525,7 @@ bool NETGENPlugin_Mesher::Compute()
   if ( true /*isOK*/ ) // get whatever built
     FillSMesh( occgeo, *ngMesh, initState, *_mesh, nodeVec, comment ); //!< 
 
-  SMESH_ComputeErrorPtr readErr = readErrors(nodeVec);
+  SMESH_ComputeErrorPtr readErr = ReadErrors(nodeVec);
   if ( readErr && !readErr->myBadElements.empty() )
     error = readErr;
 
@@ -2605,7 +2872,7 @@ void NETGENPlugin_Mesher::RemoveTmpFiles()
 //================================================================================
 
 SMESH_ComputeErrorPtr
-NETGENPlugin_Mesher::readErrors(const vector<const SMDS_MeshNode* >& nodeVec)
+NETGENPlugin_Mesher::ReadErrors(const vector<const SMDS_MeshNode* >& nodeVec)
 {
   SMESH_ComputeErrorPtr err = SMESH_ComputeError::New
     (COMPERR_BAD_INPUT_MESH, "Some edges multiple times in surface mesh");
@@ -2688,7 +2955,7 @@ void NETGENPlugin_Mesher::toPython( const netgen::Mesh* ngMesh,
     outfile << "mesh.AddNode( ";
     outfile << (*ngMesh)[pi](0) << ", ";
     outfile << (*ngMesh)[pi](1) << ", ";
-    outfile << (*ngMesh)[pi](2) << ")" << endl;
+    outfile << (*ngMesh)[pi](2) << ") ## "<< pi << endl;
   }
 
   int nbDom = ngMesh->GetNDomains();
@@ -2702,6 +2969,7 @@ void NETGENPlugin_Mesher::toPython( const netgen::Mesh* ngMesh,
     Element2d sel = (*ngMesh)[sei];
     for (int j = 0; j < sel.GetNP(); j++)
       outfile << sel[j] << ( j+1 < sel.GetNP() ? ", " : " ])");
+    if ( sel.IsDeleted() ) outfile << " ## IsDeleted ";
     outfile << endl;
 
     if ((*ngMesh)[sei].GetIndex())

src/NETGENPlugin/NETGENPlugin_Mesher.hxx

@@ -88,6 +88,7 @@ class NETGENPLUGIN_EXPORT NETGENPlugin_Mesher
 
   void SetParameters(const NETGENPlugin_Hypothesis*          hyp);
   void SetParameters(const NETGENPlugin_SimpleHypothesis_2D* hyp);
+  void SetViscousLayers2DAssigned(bool isAssigned) { _isViscousLayers2D = isAssigned; }
 
   bool Compute();
 
@@ -111,31 +112,38 @@ class NETGENPLUGIN_EXPORT NETGENPlugin_Mesher
                        std::vector<const SMDS_MeshNode*>&  nodeVec,
                        SMESH_Comment&                      comment);
 
-  bool fillNgMesh(netgen::OCCGeometry&                occgeom,
+  bool FillNgMesh(netgen::OCCGeometry&                occgeom,
                   netgen::Mesh&                       ngMesh,
                   std::vector<const SMDS_MeshNode*>&  nodeVec,
                   const std::list< SMESH_subMesh* > & meshedSM,
                   SMESH_ProxyMesh::Ptr                proxyMesh=SMESH_ProxyMesh::Ptr());
 
-  static void fixIntFaces(const netgen::OCCGeometry& occgeom,
+  static void FixIntFaces(const netgen::OCCGeometry& occgeom,
                           netgen::Mesh&              ngMesh,
                           NETGENPlugin_Internals&    internalShapes);
 
-  static void addIntVerticesInFaces(const netgen::OCCGeometry&          occgeom,
+  static void AddIntVerticesInFaces(const netgen::OCCGeometry&          occgeom,
                                     netgen::Mesh&                       ngMesh,
                                     std::vector<const SMDS_MeshNode*>&  nodeVec,
                                     NETGENPlugin_Internals&             internalShapes);
 
-  static void addIntVerticesInSolids(const netgen::OCCGeometry&         occgeom,
+  static void AddIntVerticesInSolids(const netgen::OCCGeometry&         occgeom,
                                     netgen::Mesh&                       ngMesh,
                                     std::vector<const SMDS_MeshNode*>&  nodeVec,
                                     NETGENPlugin_Internals&             internalShapes);
 
-  void defaultParameters();
+  static SMESH_ComputeErrorPtr
+    AddSegmentsToMesh(netgen::Mesh&                         ngMesh,
+                      netgen::OCCGeometry&                  geom,
+                      const TSideVector&                    wires,
+                      SMESH_MesherHelper&                   helper,
+                      std::vector< const SMDS_MeshNode* > & nodeVec);
+
+  void SetDefaultParameters();
 
   static void RemoveTmpFiles();
 
-  static SMESH_ComputeErrorPtr readErrors(const std::vector< const SMDS_MeshNode* >& nodeVec);
+  static SMESH_ComputeErrorPtr ReadErrors(const std::vector< const SMDS_MeshNode* >& nodeVec);
 
 
   static void toPython( const netgen::Mesh* ngMesh,
@@ -147,9 +155,9 @@ class NETGENPLUGIN_EXPORT NETGENPlugin_Mesher
   bool                 _isVolume;
   bool                 _optimize;
   int                  _fineness;
+  bool                 _isViscousLayers2D;
 
   const NETGENPlugin_SimpleHypothesis_2D * _simpleHyp;
-  std::map< int, std::pair<int,int> >      _faceDescriptors;
 };
 
 //=============================================================================

src/NETGENPlugin/NETGENPlugin_NETGEN_2D.cxx

@@ -32,10 +32,11 @@
 #include "NETGENPlugin_SimpleHypothesis_2D.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_ViscousLayers2D.hxx>
-#include <SMESHDS_Mesh.hxx>
 #include <utilities.h>
 
 #include <list>
@@ -64,10 +65,11 @@ NETGENPlugin_NETGEN_2D::NETGENPlugin_NETGEN_2D(int hypId, int studyId,
   _shapeType = (1 << TopAbs_FACE); // 1 bit /shape type
   _compatibleHypothesis.push_back("NETGEN_Parameters_2D");
   _compatibleHypothesis.push_back("NETGEN_SimpleParameters_2D");
+  _compatibleHypothesis.push_back( StdMeshers_ViscousLayers2D::GetHypType() );
   _requireDiscreteBoundary = false;
-  _onlyUnaryInput = false;
+  _onlyUnaryInput          = false;
-  _hypothesis = NULL;
+  _hypothesis              = NULL;
-  _supportSubmeshes = true;
+  _supportSubmeshes        = true;
 }
 
 //=============================================================================
@@ -87,33 +89,28 @@ NETGENPlugin_NETGEN_2D::~NETGENPlugin_NETGEN_2D()
  */
 //=============================================================================
 
-bool NETGENPlugin_NETGEN_2D::CheckHypothesis
+bool NETGENPlugin_NETGEN_2D::CheckHypothesis (SMESH_Mesh&         aMesh,
-                         (SMESH_Mesh& aMesh,
+                                              const TopoDS_Shape& aShape,
-                          const TopoDS_Shape& aShape,
+                                              Hypothesis_Status&  aStatus)
-                          SMESH_Hypothesis::Hypothesis_Status& aStatus)
 {
-  _hypothesis = NULL;
+  _hypothesis        = NULL;
+  _isViscousLayers2D = false;
 
-  const list<const SMESHDS_Hypothesis*>& hyps = GetUsedHypothesis(aMesh, aShape);
+  // can work with no hypothesis
-  int nbHyp = hyps.size();
+  aStatus = SMESH_Hypothesis::HYP_OK;
-  if (!nbHyp)
-  {
-    aStatus = SMESH_Hypothesis::HYP_OK;
-    return true;  // can work with no hypothesis
-  }
-  // use only the first hypothesis
-  const SMESHDS_Hypothesis* theHyp = hyps.front();
 
-  string hypName = theHyp->GetName();
+  const list<const SMESHDS_Hypothesis*>& hyps = GetUsedHypothesis(aMesh, aShape, /*skipAux=*/false);
-  if ( find( _compatibleHypothesis.begin(), _compatibleHypothesis.end(),
+  list<const SMESHDS_Hypothesis*>::const_iterator h = hyps.begin();
-             hypName ) != _compatibleHypothesis.end() )
+  for ( ; h != hyps.end(); ++h )
   {
-    _hypothesis = theHyp;
+    const SMESHDS_Hypothesis* theHyp = *h;
-    aStatus = SMESH_Hypothesis::HYP_OK;
+    string hypName = theHyp->GetName();
-  }
+    if ( hypName == StdMeshers_ViscousLayers2D::GetHypType() )
-  else
+      _isViscousLayers2D = true;
-  {
+    else if ( _hypothesis )
-    aStatus = SMESH_Hypothesis::HYP_INCOMPATIBLE;
+      aStatus = SMESH_Hypothesis::HYP_INCOMPATIBLE;
+    else
+      _hypothesis = theHyp;
   }
 
   return aStatus == SMESH_Hypothesis::HYP_OK;
@@ -132,9 +129,10 @@ bool NETGENPlugin_NETGEN_2D::Compute(SMESH_Mesh&         aMesh,
   netgen::multithread.terminate = 0;
 #endif
 
-  NETGENPlugin_Mesher mesher(&aMesh, aShape, false);
+  NETGENPlugin_Mesher mesher(&aMesh, aShape, /*is3D = */false);
   mesher.SetParameters(dynamic_cast<const NETGENPlugin_Hypothesis*>(_hypothesis));
   mesher.SetParameters(dynamic_cast<const NETGENPlugin_SimpleHypothesis_2D*>(_hypothesis));
+  mesher.SetViscousLayers2DAssigned( _isViscousLayers2D );
   return mesher.Compute();
 }
 

src/NETGENPlugin/NETGENPlugin_NETGEN_2D.hxx

@@ -59,6 +59,7 @@ public:
 
 protected:
   const SMESHDS_Hypothesis* _hypothesis;
+  bool                      _isViscousLayers2D;
 };
 
 #endif

src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.cxx

@@ -38,6 +38,7 @@
 #include <StdMeshers_LengthFromEdges.hxx>
 #include <StdMeshers_MaxElementArea.hxx>
 #include <StdMeshers_QuadranglePreference.hxx>
+#include <StdMeshers_ViscousLayers2D.hxx>
 
 #include <Precision.hxx>
 #include <Standard_ErrorHandler.hxx>
@@ -91,6 +92,7 @@ NETGENPlugin_NETGEN_2D_ONLY::NETGENPlugin_NETGEN_2D_ONLY(int hypId, int studyId,
   _compatibleHypothesis.push_back("LengthFromEdges");
   _compatibleHypothesis.push_back("QuadranglePreference");
   _compatibleHypothesis.push_back("NETGEN_Parameters_2D");
+  _compatibleHypothesis.push_back("ViscousLayers2D");
 
   _hypMaxElementArea = 0;
   _hypLengthFromEdges = 0;
@@ -148,6 +150,8 @@ bool NETGENPlugin_NETGEN_2D_ONLY::CheckHypothesis (SMESH_Mesh&         aMesh,
       _hypQuadranglePreference = static_cast<const StdMeshers_QuadranglePreference*>(hyp);
     else if ( hypName == "NETGEN_Parameters_2D" )
       _hypParameters = static_cast<const NETGENPlugin_Hypothesis_2D*>(hyp);
+    else if ( hypName == StdMeshers_ViscousLayers2D::GetHypType() )
+      continue;
     else {
       aStatus = HYP_INCOMPATIBLE;
       return false;
@@ -157,213 +161,12 @@ bool NETGENPlugin_NETGEN_2D_ONLY::CheckHypothesis (SMESH_Mesh&         aMesh,
   int nbHyps = bool(_hypMaxElementArea) + bool(_hypLengthFromEdges) + bool(_hypParameters );
   if ( nbHyps > 1 )
     aStatus = HYP_CONCURENT;
-  else if ( nbHyps == 1)
+  else
     aStatus = HYP_OK;
 
   return ( aStatus == HYP_OK );
 }
 
-//================================================================================
-/*!
- * \brief Fill netgen mesh with segments
-  * \retval SMESH_ComputeErrorPtr - error description
- */
-//================================================================================
-
-static TError addSegmentsToMesh(netgen::Mesh&                    ngMesh,
-                                OCCGeometry&                     geom,
-                                const TSideVector&               wires,
-                                SMESH_MesherHelper&              helper,
-                                vector< const SMDS_MeshNode* > & nodeVec)
-{
-  // ----------------------------
-  // Check wires and count nodes
-  // ----------------------------
-  int nbNodes = 0;
-  double totalLength = 0;
-  for ( int iW = 0; iW < wires.size(); ++iW )
-  {
-    StdMeshers_FaceSidePtr wire = wires[ iW ];
-    if ( wire->MissVertexNode() )
-    {
-      // Commented for issue 0020960. It worked for the case, let's wait for case where it doesn't.
-      // It seems that there is no reason for this limitation
-//       return TError
-//         (new SMESH_ComputeError(COMPERR_BAD_INPUT_MESH, "Missing nodes on vertices"));
-      if (getenv("USER") && string("eap")==getenv("USER"))
-        cout << "Warning: NETGENPlugin_NETGEN_2D_ONLY : try to work with missing nodes on vertices"<<endl;
-    }
-    const vector<UVPtStruct>& uvPtVec = wire->GetUVPtStruct();
-    if ( uvPtVec.size() != wire->NbPoints() )
-      return TError
-        (new SMESH_ComputeError(COMPERR_BAD_INPUT_MESH,
-                                SMESH_Comment("Unexpected nb of points on wire ") << iW
-                                << ": " << uvPtVec.size()<<" != "<<wire->NbPoints()));
-    nbNodes += wire->NbPoints();
-    totalLength += wire->Length();
-  }
-  nodeVec.reserve( nbNodes );
-
-  // -----------------
-  // Fill netgen mesh
-  // -----------------
-
-//   netgen::Box<3> bb = geom.GetBoundingBox();
-//   bb.Increase (bb.Diam()/10);
-//   ngMesh.SetLocalH (bb.PMin(), bb.PMax(), 0.5); // set grading
-
-  // map for nodes on vertices since they can be shared between wires
-  // ( issue 0020676, face_int_box.brep)
-  map<const SMDS_MeshNode*, int > node2ngID;
-
-  const int faceID = 1, solidID = 0;
-  if ( ngMesh.GetNFD() < 1 )
-    ngMesh.AddFaceDescriptor (FaceDescriptor(faceID, solidID, solidID, 0));
-
-  for ( int iW = 0; iW < wires.size(); ++iW )
-  {
-    StdMeshers_FaceSidePtr wire = wires[ iW ];
-    const vector<UVPtStruct>& uvPtVec = wire->GetUVPtStruct();
-    const int nbSegments = wire->NbPoints() - 1;
-
-    // compute length of every segment
-    vector<double> segLen( nbSegments );
-    for ( int i = 0; i < nbSegments; ++i )
-      segLen[i] = SMESH_TNodeXYZ( uvPtVec[ i ].node ).Distance( uvPtVec[ i+1 ].node );
-
-    int edgeID = 1, posID = -2;
-    bool isInternalWire = false;
-    for ( int i = 0; i < nbSegments; ++i ) // loop on segments
-    {
-      // Add the first point of a segment
-      const SMDS_MeshNode * n = uvPtVec[ i ].node;
-      const int posShapeID = n->getshapeId();
-      bool onVertex = ( n->GetPosition()->GetTypeOfPosition() == SMDS_TOP_VERTEX );
-
-      // skip nodes on degenerated edges
-      if ( helper.IsDegenShape( posShapeID ) &&
-           helper.IsDegenShape( uvPtVec[ i+1 ].node->getshapeId() ))
-        continue;
-
-      int ngID1 = ngMesh.GetNP() + 1, ngID2 = ngID1+1;
-      if ( onVertex )
-        ngID1 = node2ngID.insert( make_pair( n, ngID1 )).first->second;
-      if ( ngID1 > ngMesh.GetNP() )
-      {
-        MeshPoint mp( Point<3> (n->X(), n->Y(), n->Z()) );
-        ngMesh.AddPoint ( mp, 1, EDGEPOINT );
-        nodeVec.push_back( n );
-      }
-      else
-      {
-        ngID2 = ngMesh.GetNP() + 1;
-        if ( i > 0 ) // prev segment belongs to same wire
-        {
-          Segment& prevSeg = ngMesh.LineSegment( ngMesh.GetNSeg() );
-          prevSeg[1] = ngID1;
-        }
-      }
-
-      // Add the segment
-      Segment seg;
-
-      seg[0] = ngID1;  // ng node id
-      seg[1] = ngID2;  // ng node id
-
-      seg.edgenr = ngMesh.GetNSeg() + 1;// segment id
-      seg.si = faceID;                  // = geom.fmap.FindIndex (face);
-
-      for ( int iEnd = 0; iEnd < 2; ++iEnd)
-      {
-        const UVPtStruct& pnt = uvPtVec[ i + iEnd ];
-
-        seg.epgeominfo[ iEnd ].dist = pnt.param; // param on curve
-        seg.epgeominfo[ iEnd ].u    = pnt.u;
-        seg.epgeominfo[ iEnd ].v    = pnt.v;
-
-        // find out edge id and node parameter on edge
-        onVertex = ( pnt.node->GetPosition()->GetTypeOfPosition() == SMDS_TOP_VERTEX );
-        if ( onVertex || posShapeID != posID )
-        {
-          // get edge id
-          double normParam = pnt.normParam;
-          if ( onVertex )
-            normParam = 0.5 * ( uvPtVec[ i ].normParam + uvPtVec[ i+1 ].normParam );
-          const TopoDS_Edge& edge = wire->Edge( wire->EdgeIndex( normParam ));
-          edgeID = geom.emap.FindIndex( edge );
-          posID  = posShapeID;
-          isInternalWire = ( edge.Orientation() == TopAbs_INTERNAL );
-          // if ( onVertex ) // param on curve is different on each of two edges
-          //   seg.epgeominfo[ iEnd ].dist = helper.GetNodeU( edge, pnt.node );
-        }
-        seg.epgeominfo[ iEnd ].edgenr = edgeID; //  = geom.emap.FindIndex(edge);
-      }
-
-      ngMesh.AddSegment (seg);
-      {
-        // restrict size of elements near the segment
-        SMESH_TNodeXYZ np1( n ), np2( uvPtVec[ i+1 ].node );
-        // get an average size of adjacent segments to avoid sharp change of
-        // element size (regression on issue 0020452, note 0010898)
-        int iPrev = SMESH_MesherHelper::WrapIndex( i-1, nbSegments );
-        int iNext = SMESH_MesherHelper::WrapIndex( i+1, nbSegments );
-        double avgH = ( segLen[ iPrev ] + segLen[ i ] + segLen[ iNext ]) / 3;
-
-        NETGENPlugin_Mesher::RestrictLocalSize( ngMesh, 0.5*(np1+np2), avgH );
-      }
-#ifdef DUMP_SEGMENTS
-        cout << "Segment: " << seg.edgenr << endl
-           << "\tp1: " << seg[0] << endl
-           << "\tp2: " << seg[1] << endl
-           << "\tp0 param: " << seg.epgeominfo[ 0 ].dist << endl
-           << "\tp0 uv: " << seg.epgeominfo[ 0 ].u <<", "<< seg.epgeominfo[ 0 ].v << endl
-           << "\tp0 edge: " << seg.epgeominfo[ 0 ].edgenr << endl
-           << "\tp1 param: " << seg.epgeominfo[ 1 ].dist << endl
-           << "\tp1 uv: " << seg.epgeominfo[ 1 ].u <<", "<< seg.epgeominfo[ 1 ].v << endl
-           << "\tp1 edge: " << seg.epgeominfo[ 1 ].edgenr << endl;
-#endif
-      if ( isInternalWire )
-      {
-        swap (seg[0], seg[1]);
-        swap( seg.epgeominfo[0], seg.epgeominfo[1] );
-        seg.edgenr = ngMesh.GetNSeg() + 1; // segment id
-        ngMesh.AddSegment (seg);
-#ifdef DUMP_SEGMENTS
-        cout << "Segment: " << seg.edgenr << endl << "\tis REVRESE of the previous one" << endl;
-#endif
-      }
-    } // loop on segments on a wire
-
-    // close chain of segments
-    if ( nbSegments > 0 )
-    {
-      Segment& lastSeg = ngMesh.LineSegment( ngMesh.GetNSeg() - int( isInternalWire));
-      const SMDS_MeshNode * lastNode = uvPtVec.back().node;
-      lastSeg[1] = node2ngID.insert( make_pair( lastNode, lastSeg[1] )).first->second;
-      if ( lastSeg[1] > ngMesh.GetNP() )
-      {
-        MeshPoint mp( Point<3> (lastNode->X(), lastNode->Y(), lastNode->Z()) );
-        ngMesh.AddPoint ( mp, 1, EDGEPOINT );
-        nodeVec.push_back( lastNode );
-      }
-      if ( isInternalWire )
-      {
-        Segment& realLastSeg = ngMesh.LineSegment( ngMesh.GetNSeg() );
-        realLastSeg[0] = lastSeg[1];
-      }
-    }
-
-  } // loop on wires of a face
-
-  // add a segment instead of internal vertex
-  NETGENPlugin_Internals intShapes( *helper.GetMesh(), helper.GetSubShape(), /*is3D=*/false );
-  NETGENPlugin_Mesher::addIntVerticesInFaces( geom, ngMesh, nodeVec, intShapes );
-
-  ngMesh.CalcSurfacesOfNode();
-
-  return TError();
-}
-
 //=============================================================================
 /*!
  *Here we are going to use the NETGEN mesher
@@ -386,12 +189,18 @@ bool NETGENPlugin_NETGEN_2D_ONLY::Compute(SMESH_Mesh&         aMesh,
   helper.SetElementsOnShape( true );
   const bool ignoreMediumNodes = _quadraticMesh;
   
+  // build viscous layers if required
+  const TopoDS_Face F = TopoDS::Face( aShape.Oriented( TopAbs_FORWARD ));
+  SMESH_ProxyMesh::Ptr proxyMesh = StdMeshers_ViscousLayers2D::Compute( aMesh, F );
+  if ( !proxyMesh )
+    return false;
+
   // ------------------------
   // get all edges of a face
   // ------------------------
-  const TopoDS_Face F = TopoDS::Face( aShape.Oriented( TopAbs_FORWARD ));
   TError problem;
-  TSideVector wires = StdMeshers_FaceSide::GetFaceWires( F, aMesh, ignoreMediumNodes, problem );
+  TSideVector wires =
+    StdMeshers_FaceSide::GetFaceWires( F, aMesh, ignoreMediumNodes, problem, proxyMesh );
   if ( problem && !problem->IsOK() )
     return error( problem );
   int nbWires = wires.size();
@@ -459,7 +268,7 @@ bool NETGENPlugin_NETGEN_2D_ONLY::Compute(SMESH_Mesh&         aMesh,
   ngMesh->SetGlobalH (netgen::mparam.maxh);
 
   vector< const SMDS_MeshNode* > nodeVec;
-  problem = addSegmentsToMesh( *ngMesh, occgeo, wires, helper, nodeVec );
+  problem = aMesher.AddSegmentsToMesh( *ngMesh, occgeo, wires, helper, nodeVec );
   if ( problem && !problem->IsOK() )
     return error( problem );
 
@@ -506,19 +315,19 @@ bool NETGENPlugin_NETGEN_2D_ONLY::Compute(SMESH_Mesh&         aMesh,
   int nbNodes = ngMesh->GetNP();
   int nbFaces = ngMesh->GetNSE();
 
-  int nbInputNodes = nodeVec.size();
+  int nbInputNodes = nodeVec.size()-1;
-  nodeVec.resize( nbNodes, 0 );
+  nodeVec.resize( nbNodes+1, 0 );
 
   // add nodes
-  for ( int i = nbInputNodes + 1; i <= nbNodes; ++i )
+  for ( int ngID = nbInputNodes + 1; ngID <= nbNodes; ++ngID )
   {
-    const MeshPoint& ngPoint = ngMesh->Point(i);
+    const MeshPoint& ngPoint = ngMesh->Point( ngID );
 #ifdef NETGEN_NEW
     SMDS_MeshNode * node = meshDS->AddNode(ngPoint(0), ngPoint(1), ngPoint(2));
 #else
     SMDS_MeshNode * node = meshDS->AddNode(ngPoint.X(), ngPoint.Y(), ngPoint.Z());
 #endif
-    nodeVec[ i-1 ] = node;
+    nodeVec[ ngID ] = node;
   }
 
   // create faces
@@ -531,9 +340,9 @@ bool NETGENPlugin_NETGEN_2D_ONLY::Compute(SMESH_Mesh&         aMesh,
     for (j=1; j <= elem.GetNP(); ++j)
     {
       int pind = elem.PNum(j);
-      if ( pind-1 < 0 )
+      if ( pind < 1 )
         break;
-      const SMDS_MeshNode* node = nodeVec.at(pind-1);
+      const SMDS_MeshNode* node = nodeVec[ pind ];
       if ( reverse )
         nodes[ nodes.size()-j ] = node;
       else

src/NETGENPlugin/NETGENPlugin_NETGEN_3D.cxx

@@ -314,7 +314,7 @@ bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh&         aMesh,
     if ( internals.hasInternalVertexInSolid() )
     {
       netgen::OCCGeometry occgeo;
-      NETGENPlugin_Mesher::addIntVerticesInSolids( occgeo,
+      NETGENPlugin_Mesher::AddIntVerticesInSolids( occgeo,
                                                    (netgen::Mesh&) *Netgen_mesh,
                                                    nodeVec,
                                                    internals);
@@ -423,7 +423,7 @@ bool NETGENPlugin_NETGEN_3D::compute(SMESH_Mesh&                     aMesh,
 
   if ( err )
   {
-    SMESH_ComputeErrorPtr ce = NETGENPlugin_Mesher::readErrors(nodeVec);
+    SMESH_ComputeErrorPtr ce = NETGENPlugin_Mesher::ReadErrors(nodeVec);
     if ( ce && !ce->myBadElements.empty() )
       error( ce );
   }