IPAL0052448: Too thin viscous layers are constructed

diff --git a/src/SMESH/SMESH_MesherHelper.cxx b/src/SMESH/SMESH_MesherHelper.cxx
index d41ff05787741281e1ecb822fb3b89566b6aed2d..819b13c4cab20857dbef1a9df3c93fb400c110eb 100644 (file)
--- a/src/SMESH/SMESH_MesherHelper.cxx
+++ b/src/SMESH/SMESH_MesherHelper.cxx
@@ -2732,24 +2732,22 @@ double SMESH_MesherHelper::MaxTolerance( const TopoDS_Shape& shape )
  */
 //================================================================================
 
-double SMESH_MesherHelper::GetAngle( const TopoDS_Edge & theE1,
-                                     const TopoDS_Edge & theE2,
-                                     const TopoDS_Face & theFace,
-                                     gp_Vec*             theFaceNormal)
+double SMESH_MesherHelper::GetAngle( const TopoDS_Edge &   theE1,
+                                     const TopoDS_Edge &   theE2,
+                                     const TopoDS_Face &   theFace,
+                                     const TopoDS_Vertex & theCommonV,
+                                     gp_Vec*               theFaceNormal)
 {
   double angle = 1e100;
   try
   {
-    TopoDS_Vertex vCommon;
-    if ( !TopExp::CommonVertex( theE1, theE2, vCommon ))
-      return angle;
     double f,l;
     Handle(Geom_Curve)     c1 = BRep_Tool::Curve( theE1, f,l );
     Handle(Geom_Curve)     c2 = BRep_Tool::Curve( theE2, f,l );
     Handle(Geom2d_Curve) c2d1 = BRep_Tool::CurveOnSurface( theE1, theFace, f,l );
     Handle(Geom_Surface) surf = BRep_Tool::Surface( theFace );
-    double                 p1 = BRep_Tool::Parameter( vCommon, theE1 );
-    double                 p2 = BRep_Tool::Parameter( vCommon, theE2 );
+    double                 p1 = BRep_Tool::Parameter( theCommonV, theE1 );
+    double                 p2 = BRep_Tool::Parameter( theCommonV, theE2 );
     if ( c1.IsNull() || c2.IsNull() )
       return angle;
     gp_XY uv = c2d1->Value( p1 ).XY();
@@ -2758,10 +2756,10 @@ double SMESH_MesherHelper::GetAngle( const TopoDS_Edge & theE1,
     gp_Vec vec1, vec2, vecRef = du ^ dv;
     int  nbLoops = 0;
     double p1tmp = p1;
-    while ( vecRef.SquareMagnitude() < std::numeric_limits<double>::min() )
+    while ( vecRef.SquareMagnitude() < 1e-25 )
     {
       double dp = ( l - f ) / 1000.;
-      p1tmp += dp * (( Abs( p1 - f ) > Abs( p1 - l )) ? +1. : -1.);
+      p1tmp += dp * (( Abs( p1 - f ) > Abs( p1 - l )) ? -1. : +1.);
       uv = c2d1->Value( p1tmp ).XY();
       surf->D1( uv.X(), uv.Y(), p, du, dv );
       vecRef = du ^ dv;

diff --git a/src/SMESH/SMESH_MesherHelper.hxx b/src/SMESH/SMESH_MesherHelper.hxx
index 91d7d10411ac35a193e61ac081b50fdd59fc7846..7c8ddbdcf3ddd7dd4a2137b19f6c43846f75d062 100644 (file)
--- a/src/SMESH/SMESH_MesherHelper.hxx
+++ b/src/SMESH/SMESH_MesherHelper.hxx
@@ -218,7 +218,8 @@ class SMESH_EXPORT SMESH_MesherHelper
   static double MaxTolerance( const TopoDS_Shape& shape );
 
   static double GetAngle( const TopoDS_Edge & E1, const TopoDS_Edge & E2,
-                          const TopoDS_Face & F, gp_Vec* faceNormal=0);
+                          const TopoDS_Face & F,  const TopoDS_Vertex & V,
+                          gp_Vec* faceNormal=0);
 
   static bool IsClosedEdge( const TopoDS_Edge& anEdge );
 

diff --git a/src/StdMeshers/StdMeshers_Prism_3D.cxx b/src/StdMeshers/StdMeshers_Prism_3D.cxx
index f06261e78323bdb848af4d6fc03c0508461d1315..99fdbe90fc90f2844f79c88dd166a62dde373ff6 100644 (file)
--- a/src/StdMeshers/StdMeshers_Prism_3D.cxx
+++ b/src/StdMeshers/StdMeshers_Prism_3D.cxx
@@ -427,10 +427,10 @@ namespace {
    */
   //================================================================================
 
-  double normAngle(const TopoDS_Edge & E1, const TopoDS_Edge & E2, const TopoDS_Face & F)
-  {
-    return SMESH_MesherHelper::GetAngle( E1, E2, F ) / ( 0.5 * M_PI );
-  }
+  // double normAngle(const TopoDS_Edge & E1, const TopoDS_Edge & E2, const TopoDS_Face & F)
+  // {
+  //   return SMESH_MesherHelper::GetAngle( E1, E2, F ) / ( 0.5 * M_PI );
+  // }
 
   //================================================================================
   /*!

diff --git a/src/StdMeshers/StdMeshers_Quadrangle_2D.cxx b/src/StdMeshers/StdMeshers_Quadrangle_2D.cxx
index 6ff2210495b413f2bdfe3fc496809c0c47696e6b..f720468f7bd3bfce81d71f3dee855186610ac29c 100644 (file)
--- a/src/StdMeshers/StdMeshers_Quadrangle_2D.cxx
+++ b/src/StdMeshers/StdMeshers_Quadrangle_2D.cxx
@@ -4070,7 +4070,7 @@ bool StdMeshers_Quadrangle_2D::check()
       int iPrev = myHelper->WrapIndex( i-1, wire->NbEdges() );
       const TopoDS_Edge& e1 = wire->Edge( iPrev );
       const TopoDS_Edge& e2 = wire->Edge( i );
-      double angle = myHelper->GetAngle( e1, e2, geomFace );
+      double angle = myHelper->GetAngle( e1, e2, geomFace, wire->FirstVertex( i ));
       if (( maxAngle < angle ) &&
           ( 5.* M_PI/180 < angle && angle < 175.* M_PI/180  ))
       {
@@ -4224,7 +4224,7 @@ int StdMeshers_Quadrangle_2D::getCorners(const TopoDS_Face&          theFace,
     TopoDS_Vertex v = helper.IthVertex( 0, *edge );
     if ( !theConsiderMesh || SMESH_Algo::VertexNode( v, helper.GetMeshDS() ))
     {
-      double angle = SMESH_MesherHelper::GetAngle( prevE, *edge, theFace );
+      double angle = SMESH_MesherHelper::GetAngle( prevE, *edge, theFace, v );
       vertexByAngle.insert( make_pair( angle, v ));
       angleByVertex.Bind( v, angle );
     }

diff --git a/src/StdMeshers/StdMeshers_ViscousLayers.cxx b/src/StdMeshers/StdMeshers_ViscousLayers.cxx
index 10006f4810f3025fad58816c168ac630b7116adb..241a21636aa89f76dc89d610accf9120ac302e8c 100644 (file)
--- a/src/StdMeshers/StdMeshers_ViscousLayers.cxx
+++ b/src/StdMeshers/StdMeshers_ViscousLayers.cxx
@@ -922,7 +922,7 @@ namespace
 
       // get angle between the 2 edges
       gp_Vec faceNormal;
-      double angle = helper.GetAngle( edges[0], edges[1], faceFrw, &faceNormal );
+      double angle = helper.GetAngle( edges[0], edges[1], faceFrw, fromV, &faceNormal );
       if ( Abs( angle ) < 5 * M_PI/180 )
       {
         dir = ( faceNormal.XYZ() ^ edgeDir[0].Reversed()) + ( faceNormal.XYZ() ^ edgeDir[1] );
@@ -1004,7 +1004,8 @@ namespace
         while ( SMESH_Algo::isDegenerated( wires[iW]->Edge( iE2 )))
           iE2 = ( iE2 + 1 ) % nbEdges;
         double angle = helper.GetAngle( wires[iW]->Edge( iE1 ),
-                                        wires[iW]->Edge( iE2 ), F );
+                                        wires[iW]->Edge( iE2 ), F,
+                                        wires[iW]->FirstVertex( iE2 ));
         if ( angle < -5. * M_PI / 180. )
           return true;
       }
@@ -2283,12 +2284,19 @@ gp_XYZ _ViscousBuilder::getFaceNormal(const SMDS_MeshNode* node,
   isOK = false;
 
   Handle(Geom_Surface) surface = BRep_Tool::Surface( face );
-  if ( GeomLib::NormEstim( surface, uv, 1e-10, normal ) < 3 )
+  int pointKind = GeomLib::NormEstim( surface, uv, 1e-5, normal );
+  enum { REGULAR = 0, QUASYSINGULAR, CONICAL, IMPOSSIBLE };
+  if ( pointKind < IMPOSSIBLE )
   {
-    normal;
+    if ( pointKind != REGULAR && !shiftInside )
+    {
+      gp_XYZ normShift = getFaceNormal( node, face, helper, isOK, /*shiftInside=*/true );
+      if ( normShift * normal.XYZ() < 0. )
+        normal = normShift;
+    }
     isOK = true;
   }
-  else // hard singularity
+  else // hard singularity, to call with shiftInside=true ?
   {
     const TGeomID faceID = helper.GetMeshDS()->ShapeToIndex( face );
 
@@ -2371,12 +2379,10 @@ gp_XYZ _ViscousBuilder::getWeigthedNormal( const SMDS_MeshNode*         n,
     }
     else
     {
-      TopoDS_Vertex v10 = SMESH_MesherHelper::IthVertex( 1, ee[ 0 ]);
-      TopoDS_Vertex v01 = SMESH_MesherHelper::IthVertex( 0, ee[ 1 ]);
-      if ( !v10.IsSame( v01 ))
+      if ( !V.IsSame( SMESH_MesherHelper::IthVertex( 0, ee[ 1 ] )))
         std::swap( ee[0], ee[1] );
     }
-    angles[i] = SMESH_MesherHelper::GetAngle( ee[0], ee[1], F );
+    angles[i] = SMESH_MesherHelper::GetAngle( ee[0], ee[1], F, TopoDS::Vertex( V ));
   }
 
   // compute a weighted normal