52426, 22582: EDF 8036 SMESH: ConvertToQuadratic fails with theForce3d off

For ConvertToQuadratic: fix a case with a seam edge located not at zero parameter
 For Viscous layers: fix getting normal at a corner classified as 'Reversed' (+-PI)

src/SMDS/SMDS_MeshElement.hxx

@@ -154,7 +154,7 @@ public:
   struct Filter
   {
     virtual bool operator()(const SMDS_MeshElement* e) const = 0;
-    ~Filter() {}
+    virtual ~Filter() {}
   };
   struct NonNullFilter: public Filter
   {

src/SMESH/SMESH_MeshEditor.cxx

@@ -8751,6 +8751,7 @@ void SMESH_MeshEditor::ConvertToQuadratic(const bool theForce3d, const bool theT
   aHelper.SetIsQuadratic( true );
   aHelper.SetIsBiQuadratic( theToBiQuad );
   aHelper.SetElementsOnShape(true);
+  aHelper.ToFixNodeParameters( true );
 
   // convert elements assigned to sub-meshes
   int nbCheckedElems = 0;

src/SMESH/SMESH_MesherHelper.cxx

@@ -242,34 +242,41 @@ void SMESH_MesherHelper::SetSubShape(const TopoDS_Shape& aSh)
   for ( TopExp_Explorer eF( aSh, TopAbs_FACE ); eF.More(); eF.Next() )
   {
     const TopoDS_Face& face = TopoDS::Face( eF.Current() );
-    TopLoc_Location loc;
+    // TopLoc_Location loc;
-    Handle(Geom_Surface) surface = BRep_Tool::Surface( face, loc );
+    // Handle(Geom_Surface) surface = BRep_Tool::Surface( face, loc );
 
     // if ( surface->IsUPeriodic() || surface->IsVPeriodic() ||
     //      surface->IsUClosed()   || surface->IsVClosed() )
     {
       //while ( surface->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface )))
       //surface = Handle(Geom_RectangularTrimmedSurface)::DownCast( surface )->BasisSurface();
-      GeomAdaptor_Surface surf( surface );
 
       for (TopExp_Explorer exp( face, TopAbs_EDGE ); exp.More(); exp.Next())
       {
         // look for a seam edge
-        const TopoDS_Edge& edge = TopoDS::Edge( exp.Current() );
+        TopoDS_Edge edge = TopoDS::Edge( exp.Current() );
         if ( BRep_Tool::IsClosed( edge, face )) {
           // initialize myPar1, myPar2 and myParIndex
           gp_Pnt2d uv1, uv2;
           BRep_Tool::UVPoints( edge, face, uv1, uv2 );
           if ( Abs( uv1.Coord(1) - uv2.Coord(1) ) < Abs( uv1.Coord(2) - uv2.Coord(2) ))
           {
+            double u1 = uv1.Coord(1);
+            edge.Reverse();
+            BRep_Tool::UVPoints( edge, face, uv1, uv2 );
+            double u2 = uv1.Coord(1);
             myParIndex |= U_periodic;
-            myPar1[0] = surf.FirstUParameter();
+            myPar1[0] = Min( u1, u2 );
-            myPar2[0] = surf.LastUParameter();
+            myPar2[0] = Max( u1, u2 );
           }
           else {
+            double v1 = uv1.Coord(2);
+            edge.Reverse();
+            BRep_Tool::UVPoints( edge, face, uv1, uv2 );
+            double v2 = uv1.Coord(2);
             myParIndex |= V_periodic;
-            myPar1[1] = surf.FirstVParameter();
+            myPar1[1] = Min( v1, v2 );
-            myPar2[1] = surf.LastVParameter();
+            myPar2[1] = Max( v1, v2 );
           }
           // store seam shape indices, negative if shape encounters twice
           int edgeID = meshDS->ShapeToIndex( edge );
@@ -287,13 +294,15 @@ void SMESH_MesherHelper::SetSubShape(const TopoDS_Shape& aSh)
             myDegenShapeIds.insert( meshDS->ShapeToIndex( v.Current() ));
         }
       }
-      if ( !myDegenShapeIds.empty() && !myParIndex ) {
+      if ( !myDegenShapeIds.empty() && !myParIndex )
-        if ( surface->IsUPeriodic() || surface->IsUClosed() ) {
+      {
+        BRepAdaptor_Surface surf( face, false );
+        if ( surf.IsUPeriodic() || surf.IsUClosed() ) {
           myParIndex |= U_periodic;
           myPar1[0] = surf.FirstUParameter();
           myPar2[0] = surf.LastUParameter();
         }
-        else if ( surface->IsVPeriodic() || surface->IsVClosed() ) {
+        else if ( surf.IsVPeriodic() || surf.IsVClosed() ) {
           myParIndex |= V_periodic;
           myPar1[1] = surf.FirstVParameter();
           myPar2[1] = surf.LastVParameter();
@@ -513,8 +522,8 @@ gp_Pnt2d SMESH_MesherHelper::GetUVOnSeam( const gp_Pnt2d& uv1, const gp_Pnt2d& u
       double p1 = uv1.Coord( i );
       double dp1 = Abs( p1-myPar1[i-1]), dp2 = Abs( p1-myPar2[i-1]);
       if ( myParIndex == i ||
-           dp1 < ( myPar2[i-1] - myPar2[i-1] ) / 100. ||
+           dp1 < ( myPar2[i-1] - myPar1[i-1] ) / 100. ||
-           dp2 < ( myPar2[i-1] - myPar2[i-1] ) / 100. )
+           dp2 < ( myPar2[i-1] - myPar1[i-1] ) / 100. )
       {
         double p2 = uv2.Coord( i );
         double p1Alt = ( dp1 < dp2 ) ? myPar2[i-1] : myPar1[i-1];
@@ -544,7 +553,7 @@ gp_XY SMESH_MesherHelper::GetNodeUV(const TopoDS_Face&   F,
   {
     // node has position on face
     const SMDS_FacePosition* fpos =
-      static_cast<const SMDS_FacePosition*>(n->GetPosition());
+      static_cast<const SMDS_FacePosition*>( Pos );
     uv.SetCoord(fpos->GetUParameter(),fpos->GetVParameter());
     if ( check )
       uvOK = CheckNodeUV( F, n, uv.ChangeCoord(), 10*MaxTolerance( F ));
@@ -555,7 +564,7 @@ gp_XY SMESH_MesherHelper::GetNodeUV(const TopoDS_Face&   F,
     // corresponding edge from face, get pcurve for this
     // edge and retrieve value from this pcurve
     const SMDS_EdgePosition* epos =
-      static_cast<const SMDS_EdgePosition*>(n->GetPosition());
+      static_cast<const SMDS_EdgePosition*>( Pos );
     int edgeID = n->getshapeId();
     TopoDS_Edge E = TopoDS::Edge(GetMeshDS()->IndexToShape(edgeID));
     double f, l, u = epos->GetUParameter();
@@ -1457,7 +1466,6 @@ const SMDS_MeshNode* SMESH_MesherHelper::GetMediumNode(const SMDS_MeshNode* n1,
           if ( myParIndex & U_periodic ) uv[1].SetCoord( 1, uv[0].Coord( 1 ));
           else                           uv[1].SetCoord( 2, uv[0].Coord( 2 ));
         }
-
         TopLoc_Location loc;
         Handle(Geom_Surface) S = BRep_Tool::Surface(F,loc);
         gp_XY UV = GetMiddleUV( S, uv[0], uv[1] );
@@ -2728,7 +2736,8 @@ double SMESH_MesherHelper::MaxTolerance( const TopoDS_Shape& shape )
 
 double SMESH_MesherHelper::GetAngle( const TopoDS_Edge & theE1,
                                      const TopoDS_Edge & theE2,
-                                     const TopoDS_Face & theFace)
+                                     const TopoDS_Face & theFace,
+                                     gp_Vec*             theFaceNormal)
 {
   double angle = 1e100;
   try
@@ -2768,16 +2777,33 @@ double SMESH_MesherHelper::GetAngle( const TopoDS_Edge & theE1,
     }
     if ( theFace.Orientation() == TopAbs_REVERSED )
       vecRef.Reverse();
+    if ( theFaceNormal ) *theFaceNormal = vecRef;
+
     c1->D1( p1, p, vec1 );
     c2->D1( p2, p, vec2 );
-    TopoDS_Face F = theFace;
+    // TopoDS_Face F = theFace;
-    if ( F.Orientation() == TopAbs_INTERNAL )
+    // if ( F.Orientation() == TopAbs_INTERNAL )
-      F.Orientation( TopAbs_FORWARD );
+    //   F.Orientation( TopAbs_FORWARD );
     if ( theE1.Orientation() /*GetSubShapeOri( F, theE1 )*/ == TopAbs_REVERSED )
       vec1.Reverse();
     if ( theE2.Orientation() /*GetSubShapeOri( F, theE2 )*/ == TopAbs_REVERSED )
       vec2.Reverse();
     angle = vec1.AngleWithRef( vec2, vecRef );
+
+    if ( Abs ( angle ) >= 0.99 * M_PI )
+    {
+      BRep_Tool::Range( theE1, f, l );
+      p1 += 1e-7 * ( p1-f < l-p1 ? +1. : -1. );
+      c1->D1( p1, p, vec1 );
+      if ( theE1.Orientation() == TopAbs_REVERSED )
+        vec1.Reverse();
+      BRep_Tool::Range( theE2, f, l );
+      p2 += 1e-7 * ( p2-f < l-p2 ? +1. : -1. );
+      c2->D1( p2, p, vec2 );
+      if ( theE2.Orientation() == TopAbs_REVERSED )
+        vec2.Reverse();
+      angle = vec1.AngleWithRef( vec2, vecRef );
+    }
   }
   catch (...)
   {
@@ -4136,9 +4162,11 @@ namespace { // Structures used by FixQuadraticElements()
       //BRepClass3d_SolidClassifier solidClassifier( shape );
 
       TIDSortedElemSet checkedVols, movedNodes;
-      for ( faceIt.ReInit(); faceIt.More(); faceIt.Next() ) // loop on FACEs of a SOLID
+      //for ( faceIt.ReInit(); faceIt.More(); faceIt.Next() ) // loop on FACEs of a SOLID
+      for ( size_t iF = 0; iF < concaveFaces.size(); ++iF ) // loop on concave FACEs
       {
-        const TopoDS_Shape& face = faceIt.Current();
+        //const TopoDS_Shape& face = faceIt.Current();
+        const TopoDS_Shape& face = concaveFaces[ iF ];
         SMESHDS_SubMesh*  faceSM = meshDS->MeshElements( face );
         if ( !faceSM ) continue;
 
@@ -4153,11 +4181,16 @@ namespace { // Structures used by FixQuadraticElements()
           if ( !vertexSM ) continue;
           nodeIt = vertexSM->GetNodes();
         }
+        // get ids of sub-shapes of the FACE
+        set< int > subIDs;
+        SMESH_subMeshIteratorPtr smIt =
+          theHelper.GetMesh()->GetSubMesh( face )->getDependsOnIterator(/*includeSelf=*/true);
+        while ( smIt->more() )
+          subIDs.insert( smIt->next()->GetId() );
 
         // find suspicious volumes adjacent to the FACE
         vector< const SMDS_MeshNode* >    nOnFace( 4 );
         const SMDS_MeshNode*              nInSolid;
-        //vector< const SMDS_MeshElement* > intersectedFaces;
         while ( nodeIt->more() )
         {
           const SMDS_MeshNode* n     = nodeIt->next();
@@ -4179,7 +4212,7 @@ namespace { // Structures used by FixQuadraticElements()
               n = *volNode;
               if ( n->GetPosition()->GetDim() == 3 )
                 nInSolid = n;
-              else
+              else if ( subIDs.count( n->getshapeId() ))
                 nOnFace.push_back( n );
             }
             if ( !nInSolid || nOnFace.size() != nbN - 1 )

src/SMESH/SMESH_MesherHelper.hxx

@@ -217,7 +217,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);
+  static double GetAngle( const TopoDS_Edge & E1, const TopoDS_Edge & E2,
+                          const TopoDS_Face & F, gp_Vec* faceNormal=0);
 
   static bool IsClosedEdge( const TopoDS_Edge& anEdge );
 

src/StdMeshers/StdMeshers_ViscousLayers.cxx

@@ -83,7 +83,7 @@
 #include <cmath>
 #include <limits>
 
-//#define __myDEBUG
+#define __myDEBUG
 
 using namespace std;
 
@@ -486,6 +486,11 @@ namespace VISCOUS_3D
     bool makeLayer(_SolidData& data);
     bool setEdgeData(_LayerEdge& edge, const set<TGeomID>& subIds,
                      SMESH_MesherHelper& helper, _SolidData& data);
+    gp_XYZ getFaceNormal(const SMDS_MeshNode* n,
+                         const TopoDS_Face&   face,
+                         SMESH_MesherHelper&  helper,
+                         bool&                isOK,
+                         bool                 shiftInside=false);
     gp_XYZ getWeigthedNormal( const SMDS_MeshNode*         n,
                               std::pair< TGeomID, gp_XYZ > fId2Normal[],
                               const int                    nbFaces );
@@ -501,6 +506,8 @@ namespace VISCOUS_3D
                                vector< vector<_LayerEdge*> >& edgesByGeom);
     bool sortEdges( _SolidData&                    data,
                     vector< vector<_LayerEdge*> >& edgesByGeom);
+    void limitStepSizeByCurvature( _SolidData&                    data,
+                                   vector< vector<_LayerEdge*> >& edgesByGeom);
     void limitStepSize( _SolidData&             data,
                         const SMDS_MeshElement* face,
                         const double            cosin);
@@ -778,8 +785,8 @@ namespace
     if ( node->GetPosition()->GetTypeOfPosition() == SMDS_TOP_VERTEX &&
          helper.IsClosedEdge( fromE ))
     {
-      if ( fabs(u-f) < fabs(u-l )) c->D1( l, p, dv );
+      if ( fabs(u-f) < fabs(u-l)) c->D1( l, p, dv );
-      else                         c->D1( f, p, dv );
+      else                        c->D1( f, p, dv );
       if ( o == TopAbs_REVERSED )
         dv.Reverse();
       gp_Vec dir2 = norm ^ dv;
@@ -792,51 +799,71 @@ namespace
                      const SMDS_MeshNode* node, SMESH_MesherHelper& helper,
                      bool& ok, double* cosin=0)
   {
+    TopoDS_Face faceFrw = F;
+    faceFrw.Orientation( TopAbs_FORWARD );
     double f,l; TopLoc_Location loc;
-    vector< TopoDS_Edge > edges; // sharing a vertex
+    TopoDS_Edge edges[2]; // sharing a vertex
-    PShapeIteratorPtr eIt = helper.GetAncestors( fromV, *helper.GetMesh(), TopAbs_EDGE);
+    int nbEdges = 0;
-    while ( eIt->more())
     {
-      const TopoDS_Edge* e = static_cast<const TopoDS_Edge*>( eIt->next() );
+      TopoDS_Vertex VV[2];
-      if ( helper.IsSubShape( *e, F ) && !BRep_Tool::Curve( *e, loc,f,l).IsNull() )
+      TopExp_Explorer exp( faceFrw, TopAbs_EDGE );
-        edges.push_back( *e );
+      for ( ; exp.More() && nbEdges < 2; exp.Next() )
+      {
+        const TopoDS_Edge& e = TopoDS::Edge( exp.Current() );
+        if ( SMESH_Algo::isDegenerated( e )) continue;
+        TopExp::Vertices( e, VV[0], VV[1], /*CumOri=*/true );
+        if ( VV[1].IsSame( fromV )) {
+          edges[ 0 ] = e;
+          nbEdges++;
+        }
+        else if ( VV[0].IsSame( fromV )) {
+          edges[ 1 ] = e;
+          nbEdges++;
+        }
+      }
     }
-    gp_XYZ dir(0,0,0);
+    gp_XYZ dir(0,0,0), edgeDir[2];
-    if ( !( ok = ( edges.size() > 0 ))) return dir;
+    if ( nbEdges == 2 )
-    // get average dir of edges going fromV
-    gp_XYZ edgeDir;
-    //if ( edges.size() > 1 )
-    for ( size_t i = 0; i < edges.size(); ++i )
     {
-      edgeDir = getEdgeDir( edges[i], fromV );
+      // get dirs of edges going fromV
-      double size2 = edgeDir.SquareModulus();
+      ok = true;
-      if ( size2 > numeric_limits<double>::min() )
+      for ( size_t i = 0; i < nbEdges && ok; ++i )
-        edgeDir /= sqrt( size2 );
+      {
+        edgeDir[i] = getEdgeDir( edges[i], fromV );
+        double size2 = edgeDir[i].SquareModulus();
+        if (( ok = size2 > numeric_limits<double>::min() ))
+          edgeDir[i] /= sqrt( size2 );
+      }
+      if ( !ok ) return dir;
+
+      // get angle between the 2 edges
+      gp_Vec faceNormal;
+      double angle = helper.GetAngle( edges[0], edges[1], faceFrw, &faceNormal );
+      if ( Abs( angle ) < 5 * M_PI/180 )
+      {
+        dir = ( faceNormal.XYZ() ^ edgeDir[0].Reversed()) + ( faceNormal.XYZ() ^ edgeDir[1] );
+      }
       else
-        ok = false;
+      {
-      dir += edgeDir;
+        dir = edgeDir[0] + edgeDir[1];
+        if ( angle < 0 )
+          dir.Reverse();
+      }
+      if ( cosin ) {
+        double angle = gp_Vec( edgeDir[0] ).Angle( dir );
+        *cosin = Cos( angle );
+      }
     }
-    gp_XYZ fromEdgeDir = getFaceDir( F, edges[0], node, helper, ok );
+    else if ( nbEdges == 1 )
-    try {
+    {
-      if ( edges.size() == 1 )
+      dir = getFaceDir( faceFrw, edges[0], node, helper, ok );
-        dir = fromEdgeDir;
+      if ( cosin ) *cosin = 1.;
-      else if ( dir.SquareModulus() < 0.1 ) // ~< 20 degrees
-        dir = fromEdgeDir.Normalized() + getFaceDir( F, edges[1], node, helper, ok ).Normalized();
-      else if ( dir * fromEdgeDir < 0 )
-        dir *= -1;
     }
-    catch ( Standard_Failure )
+    else
     {
       ok = false;
     }
-    if ( ok )
+
-    {
-      //dir /= edges.size();
-      if ( cosin ) {
-        double angle = gp_Vec( edgeDir ).Angle( dir );
-        *cosin = cos( angle );
-      }
-    }
     return dir;
   }
   //================================================================================
@@ -1557,6 +1584,9 @@ bool _ViscousBuilder::makeLayer(_SolidData& data)
   // fill data._simplexTestEdges
   findSimplexTestEdges( data, edgesByGeom );
 
+  // limit data._stepSize depending on surface curvature
+  limitStepSizeByCurvature( data, edgesByGeom );
+
   // Put _LayerEdge's into the vector data._edges
   if ( !sortEdges( data, edgesByGeom ))
     return false;
@@ -1644,6 +1674,62 @@ void _ViscousBuilder::limitStepSize( _SolidData& data, const double minSize )
   }
 }
 
+//================================================================================
+/*!
+ * \brief Limit data._stepSize by evaluating curvature of shapes
+ */
+//================================================================================
+
+void _ViscousBuilder::limitStepSizeByCurvature( _SolidData&                    data,
+                                                vector< vector<_LayerEdge*> >& edgesByGeom)
+{
+  const int nbTestPnt = 5;
+  const double minCurvature = 0.9 / data._hyp->GetTotalThickness();
+
+  BRepLProp_SLProps surfProp( 2, 1e-6 );
+  SMESH_MesherHelper helper( *_mesh );
+
+  TopExp_Explorer face( data._solid, TopAbs_FACE );
+  for ( ; face.More(); face.Next() )
+  {
+    const TopoDS_Face& F = TopoDS::Face( face.Current() );
+    BRepAdaptor_Surface surface( F, false );
+    surfProp.SetSurface( surface );
+
+    SMESH_subMesh *   sm = _mesh->GetSubMesh( F );
+    SMESH_subMeshIteratorPtr smIt = sm->getDependsOnIterator(/*includeSelf=*/true);
+    while ( smIt->more() )
+    {
+      sm = smIt->next();
+      const vector<_LayerEdge*>& ledges = edgesByGeom[ sm->GetId() ];
+      int step = Max( 1, int( ledges.size()) / nbTestPnt );
+      for ( size_t i = 0; i < ledges.size(); i += step )
+      {
+        gp_XY uv = helper.GetNodeUV( F, ledges[i]->_nodes[0] );
+        surfProp.SetParameters( uv.X(), uv.Y() );
+        if ( !surfProp.IsCurvatureDefined() )
+          continue;
+        double surfCurvature = Max( Abs( surfProp.MaxCurvature() ),
+                                    Abs( surfProp.MinCurvature() ));
+        if ( surfCurvature < minCurvature )
+          continue;
+
+        gp_Dir minDir, maxDir;
+        surfProp.CurvatureDirections( maxDir, minDir );
+        if ( F.Orientation() == TopAbs_REVERSED ) {
+          maxDir.Reverse(); minDir.Reverse();
+        }
+        const gp_XYZ& inDir = ledges[i]->_normal;
+        if ( inDir * maxDir.XYZ() < 0 &&
+             inDir * minDir.XYZ() < 0 )
+          continue;
+
+        limitStepSize( data, 0.9 / surfCurvature );
+      }
+    }
+  }
+}
+
 //================================================================================
 /*!
  * Fill data._simplexTestEdges. These _LayerEdge's are used to stop inflation
@@ -1666,7 +1752,7 @@ void _ViscousBuilder::findSimplexTestEdges( _SolidData&                    data,
 
   const double minCurvature = 1. / data._hyp->GetTotalThickness();
 
-  for ( size_t iS = 0; iS < edgesByGeom.size(); ++iS )
+  for ( size_t iS = 1; iS < edgesByGeom.size(); ++iS )
   {
     // look for a FACE with layers and w/o _LayerEdge's
     const vector<_LayerEdge*>& eS = edgesByGeom[iS];
@@ -1935,37 +2021,9 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge&         edge,
       if ( s.IsNull() || s.ShapeType() != TopAbs_FACE || !subIds.count( *id ))
         continue;
       F = TopoDS::Face( s );
+      geomNorm = getFaceNormal( node, F, helper, normOK );
+      if ( !normOK ) continue;
 
-      gp_XY uv = helper.GetNodeUV( F, node, 0, &normOK );
-      Handle(Geom_Surface) surface = BRep_Tool::Surface( F );
-      {
-        gp_Dir normal;
-        if ( GeomLib::NormEstim( surface, uv, 1e-10, normal ) < 3 )
-        {
-          geomNorm = normal;
-          normOK = true;
-        }
-        else // hard singularity
-        {
-          SMDS_ElemIteratorPtr fIt = node->GetInverseElementIterator(SMDSAbs_Face);
-          while ( fIt->more() )
-          {
-            const SMDS_MeshElement* f = fIt->next();
-            if ( editor.FindShape( f ) == *id )
-            {
-              SMESH_MeshAlgos::FaceNormal( f, (gp_XYZ&) geomNorm.XYZ(), /*normalized=*/false );
-              if ( helper.IsReversedSubMesh( F ))
-                geomNorm.Reverse();
-              break;
-            }
-          }
-          double size2 = geomNorm.SquareMagnitude();
-          if ( size2 > numeric_limits<double>::min() )
-            geomNorm /= sqrt( size2 );
-          else
-            normOK = false;
-        }
-      }
       if ( helper.GetSubShapeOri( data._solid, F ) != TopAbs_REVERSED )
         geomNorm.Reverse();
       id2Norm[ totalNbFaces ].first  = *id;
@@ -1976,6 +2034,22 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge&         edge,
     if ( totalNbFaces == 0 )
       return error(SMESH_Comment("Can't get normal to node ") << node->GetID(), data._index);
 
+    if ( normOK && edge._normal.Modulus() < 1e-3 && totalNbFaces > 1 )
+    {
+      // opposite normals, re-get normals at shifted positions (IPAL 52426)
+      edge._normal.SetCoord( 0,0,0 );
+      for ( int i = 0; i < totalNbFaces; ++i )
+      {
+        const TopoDS_Face& F = TopoDS::Face( getMeshDS()->IndexToShape( id2Norm[i].first ));
+        geomNorm = getFaceNormal( node, F, helper, normOK, /*shiftInside=*/true );
+        if ( helper.GetSubShapeOri( data._solid, F ) != TopAbs_REVERSED )
+          geomNorm.Reverse();
+        if ( normOK )
+          id2Norm[ i ].second = geomNorm.XYZ();
+        edge._normal += id2Norm[ i ].second;
+      }
+    }
+
     if ( totalNbFaces < 3 )
     {
       //edge._normal /= totalNbFaces;
@@ -1992,7 +2066,7 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge&         edge,
 
     case SMDS_TOP_EDGE: {
       TopoDS_Edge E = TopoDS::Edge( helper.GetSubShapeByNode( node, getMeshDS()));
-      gp_Vec inFaceDir = getFaceDir( F, E, node, helper, normOK);
+      gp_Vec inFaceDir = getFaceDir( F, E, node, helper, normOK );
       double angle = inFaceDir.Angle( edge._normal ); // [0,PI]
       edge._cosin = cos( angle );
       //cout << "Cosin on EDGE " << edge._cosin << " node " << node->GetID() << endl;
@@ -2000,8 +2074,8 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge&         edge,
     }
     case SMDS_TOP_VERTEX: {
       TopoDS_Vertex V = TopoDS::Vertex( helper.GetSubShapeByNode( node, getMeshDS()));
-      gp_XYZ inFaceDir = getFaceDir( F, V, node, helper, normOK);
+      gp_Vec inFaceDir = getFaceDir( F, V, node, helper, normOK );
-      double angle = gp_Vec( inFaceDir).Angle( edge._normal ); // [0,PI]
+      double angle = inFaceDir.Angle( edge._normal ); // [0,PI]
       edge._cosin = cos( angle );
       //cout << "Cosin on VERTEX " << edge._cosin << " node " << node->GetID() << endl;
       break;
@@ -2009,7 +2083,7 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge&         edge,
     default:
       return error(SMESH_Comment("Invalid shape position of node ")<<node, data._index);
     }
-  }
+  } // case _sWOL.IsNull()
 
   double normSize = edge._normal.SquareModulus();
   if ( normSize < numeric_limits<double>::min() )
@@ -2087,7 +2161,102 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge&         edge,
 
 //================================================================================
 /*!
- * \brief Return normal at a node weighted with angles taken by FACEs
+ * \brief Return normal to a FACE at a node
+ *  \param [in] n - node
+ *  \param [in] face - FACE
+ *  \param [in] helper - helper
+ *  \param [out] isOK - true or false
+ *  \param [in] shiftInside - to find normal at a position shifted inside the face
+ *  \return gp_XYZ - normal
+ */
+//================================================================================
+
+gp_XYZ _ViscousBuilder::getFaceNormal(const SMDS_MeshNode* node,
+                                      const TopoDS_Face&   face,
+                                      SMESH_MesherHelper&  helper,
+                                      bool&                isOK,
+                                      bool                 shiftInside)
+{
+  gp_XY uv;
+  if ( shiftInside )
+  {
+    // get a shifted position
+    gp_Pnt p = SMESH_TNodeXYZ( node );
+    gp_XYZ shift( 0,0,0 );
+    TopoDS_Shape S = helper.GetSubShapeByNode( node, helper.GetMeshDS() );
+    switch ( S.ShapeType() ) {
+    case TopAbs_VERTEX:
+    {
+      shift = getFaceDir( face, TopoDS::Vertex( S ), node, helper, isOK );
+      break;
+    }
+    case TopAbs_EDGE:
+    {
+      shift = getFaceDir( face, TopoDS::Edge( S ), node, helper, isOK );
+      break;
+    }
+    default:
+      isOK = false;
+    }
+    if ( isOK )
+      shift.Normalize();
+    p.Translate( shift * 1e-5 );
+
+    TopLoc_Location loc;
+    GeomAPI_ProjectPointOnSurf& projector = helper.GetProjector( face, loc, 1e-7 );
+
+    if ( !loc.IsIdentity() ) p.Transform( loc.Transformation().Inverted() );
+    
+    projector.Perform( p );
+    if ( !projector.IsDone() || projector.NbPoints() < 1 )
+    {
+      isOK = false;
+      return p.XYZ();
+    }
+    Quantity_Parameter U,V;
+    projector.LowerDistanceParameters(U,V);
+    uv.SetCoord( U,V );
+  }
+  else
+  {
+    uv = helper.GetNodeUV( face, node, 0, &isOK );
+  }
+
+  gp_Dir normal;
+  isOK = false;
+
+  Handle(Geom_Surface) surface = BRep_Tool::Surface( face );
+  if ( GeomLib::NormEstim( surface, uv, 1e-10, normal ) < 3 )
+  {
+    normal;
+    isOK = true;
+  }
+  else // hard singularity
+  {
+    const TGeomID faceID = helper.GetMeshDS()->ShapeToIndex( face );
+
+    SMDS_ElemIteratorPtr fIt = node->GetInverseElementIterator(SMDSAbs_Face);
+    while ( fIt->more() )
+    {
+      const SMDS_MeshElement* f = fIt->next();
+      if ( f->getshapeId() == faceID )
+      {
+        isOK = SMESH_MeshAlgos::FaceNormal( f, (gp_XYZ&) normal.XYZ(), /*normalized=*/true );
+        if ( isOK )
+        {
+          if ( helper.IsReversedSubMesh( face ))
+            normal.Reverse();
+          break;
+        }
+      }
+    }
+  }
+  return normal.XYZ();
+}
+
+//================================================================================
+/*!
+ * \brief Return a normal at a node weighted with angles taken by FACEs
  *  \param [in] n - the node
  *  \param [in] fId2Normal - FACE ids and normals
  *  \param [in] nbFaces - nb of FACEs meeting at the node
@@ -2311,7 +2480,8 @@ gp_XYZ _LayerEdge::Copy( _LayerEdge& other, SMESH_MesherHelper& helper )
 void _LayerEdge::SetCosin( double cosin )
 {
   _cosin = cosin;
-  _lenFactor = ( Abs( _cosin ) > 0.1 ) ?  1./sqrt(1-_cosin*_cosin) : 1.0;
+  cosin = Abs( _cosin );
+  _lenFactor = ( /*0.1 < cosin &&*/ cosin < 1-1e-12 ) ?  1./sqrt(1-cosin*cosin) : 1.0;
 }
 
 //================================================================================
@@ -2939,7 +3109,7 @@ bool _ViscousBuilder::smoothAnalyticEdge( _SolidData&           data,
       gp_Vec2d vec1( center, uv1 );
       double uLast = vec0.Angle( vec1 ); // -PI - +PI
       double uMidl = vec0.Angle( vecM );
-      if ( uLast * uMidl < 0. )
+      if ( uLast * uMidl <= 0. )
         uLast += ( uMidl > 0 ? +2. : -2. ) * M_PI;
       const double radius = 0.5 * ( vec0.Magnitude() + vec1.Magnitude() );
 
@@ -3145,7 +3315,7 @@ bool _ViscousBuilder::updateNormals( _SolidData&         data,
       if ( edge1->_cosin < 0 )
         dirInFace = dir1;
       else
-        getFaceDir( FF1[0], E1, edge1->_nodes[0], helper, ok );
+        dirInFace = getFaceDir( FF1[0], E1, edge1->_nodes[0], helper, ok );
       double angle = dir1.Angle( edge1->_normal ); // [0,PI]
       edge1->SetCosin( cos( angle ));