Salome HOME
IPAL0052448: Too thin viscous layers are constructed
[modules/smesh.git] / src / StdMeshers / StdMeshers_ViscousLayers.cxx
index db924388118797749ed66a9752507a5a40ce8f65..241a21636aa89f76dc89d610accf9120ac302e8c 100644 (file)
@@ -286,6 +286,7 @@ namespace VISCOUS_3D
         c = new _Curvature;
         c->_r = avgDist * avgDist / avgNormProj;
         c->_k = avgDist * avgDist / c->_r / c->_r;
+        //c->_k = avgNormProj / c->_r;
         c->_k *= ( c->_r < 0 ? 1/1.1 : 1.1 ); // not to be too restrictive
         c->_h2lenRatio = avgNormProj / ( avgDist + avgDist );
       }
@@ -419,7 +420,7 @@ namespace VISCOUS_3D
     TopoDS_Shape                    _hypShape;
     _MeshOfSolid*                   _proxyMesh;
     set<TGeomID>                    _reversedFaceIds;
-    set<TGeomID>                    _ignoreFaceIds;
+    set<TGeomID>                    _ignoreFaceIds; // WOL FACEs and FACEs of other SOLIDS
 
     double                          _stepSize, _stepSizeCoeff;
     const SMDS_MeshNode*            _stepSizeNodes[2];
@@ -486,7 +487,7 @@ namespace VISCOUS_3D
 
     bool GetShapeEdges(const TGeomID shapeID, size_t& edgeEnd, int* iBeg=0, int* iEnd=0 ) const;
 
-    void AddFacesToSmooth( const set< TGeomID >& faceIDs );
+    void AddShapesToSmooth( const set< TGeomID >& shapeIDs );
   };
   //--------------------------------------------------------------------------------
   /*!
@@ -588,7 +589,7 @@ namespace VISCOUS_3D
     void limitStepSizeByCurvature( _SolidData&  data );
     void limitStepSize( _SolidData&             data,
                         const SMDS_MeshElement* face,
-                        const double            cosin);
+                        const _LayerEdge*       maxCosinEdge );
     void limitStepSize( _SolidData& data, const double minSize);
     bool inflate(_SolidData& data);
     bool smoothAndCheck(_SolidData& data, const int nbSteps, double & distToIntersection);
@@ -896,12 +897,12 @@ namespace
         if ( SMESH_Algo::isDegenerated( e )) continue;
         TopExp::Vertices( e, VV[0], VV[1], /*CumOri=*/true );
         if ( VV[1].IsSame( fromV )) {
+          nbEdges += edges[ 0 ].IsNull();
           edges[ 0 ] = e;
-          nbEdges++;
         }
         else if ( VV[0].IsSame( fromV )) {
+          nbEdges += edges[ 1 ].IsNull();
           edges[ 1 ] = e;
-          nbEdges++;
         }
       }
     }
@@ -921,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] );
@@ -939,7 +940,7 @@ namespace
     }
     else if ( nbEdges == 1 )
     {
-      dir = getFaceDir( faceFrw, edges[0], node, helper, ok );
+      dir = getFaceDir( faceFrw, edges[ edges[0].IsNull() ], node, helper, ok );
       if ( cosin ) *cosin = 1.;
     }
     else
@@ -1003,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;
       }
@@ -1030,8 +1032,12 @@ namespace
       theNbFunc = 0;
     }
     void Finish() {
-      if (py)
-        *py << "mesh.MakeGroup('Viscous Prisms',SMESH.VOLUME,SMESH.FT_ElemGeomType,'=',SMESH.Geom_PENTA)"<<endl;
+      if (py) {
+        *py << "mesh.GroupOnFilter(SMESH.VOLUME,'Viscous Prisms',"
+          "smesh.GetFilter(SMESH.VOLUME,SMESH.FT_ElemGeomType,'=',SMESH.Geom_PENTA))"<<endl;
+        *py << "mesh.GroupOnFilter(SMESH.VOLUME,'Neg Volumes',"
+          "smesh.GetFilter(SMESH.VOLUME,SMESH.FT_Volume3D,'<',0))"<<endl;
+      }
       delete py; py=0;
     }
     ~PyDump() { Finish(); cout << "NB FUNCTIONS: " << theNbFunc << endl; }
@@ -1592,6 +1598,7 @@ bool _ViscousBuilder::makeLayer(_SolidData& data)
       const SMDS_MeshElement* face = eIt->next();
       newNodes.resize( face->NbCornerNodes() );
       double faceMaxCosin = -1;
+      _LayerEdge* maxCosinEdge = 0;
       for ( int i = 0 ; i < face->NbCornerNodes(); ++i )
       {
         const SMDS_MeshNode* n = face->GetNode(i);
@@ -1626,10 +1633,11 @@ bool _ViscousBuilder::makeLayer(_SolidData& data)
               return false;
           }
           dumpMove(edge->_nodes.back());
-          if ( edge->_cosin > 0.01 )
+
+          if ( edge->_cosin > faceMaxCosin )
           {
-            if ( edge->_cosin > faceMaxCosin )
-              faceMaxCosin = edge->_cosin;
+            faceMaxCosin = edge->_cosin;
+            maxCosinEdge = edge;
           }
         }
         newNodes[ i ] = n2e->second->_nodes.back();
@@ -1641,7 +1649,7 @@ bool _ViscousBuilder::makeLayer(_SolidData& data)
 
       // compute inflation step size by min size of element on a convex surface
       if ( faceMaxCosin > theMinSmoothCosin )
-        limitStepSize( data, face, faceMaxCosin );
+        limitStepSize( data, face, maxCosinEdge );
     } // loop on 2D elements on a FACE
   } // loop on FACEs of a SOLID
 
@@ -1692,7 +1700,7 @@ bool _ViscousBuilder::makeLayer(_SolidData& data)
 
 void _ViscousBuilder::limitStepSize( _SolidData&             data,
                                      const SMDS_MeshElement* face,
-                                     const double            cosin)
+                                     const _LayerEdge*       maxCosinEdge )
 {
   int iN = 0;
   double minSize = 10 * data._stepSize;
@@ -1700,20 +1708,20 @@ void _ViscousBuilder::limitStepSize( _SolidData&             data,
   for ( int i = 0; i < nbNodes; ++i )
   {
     const SMDS_MeshNode* nextN = face->GetNode( SMESH_MesherHelper::WrapIndex( i+1, nbNodes ));
-    const SMDS_MeshNode* curN = face->GetNode( i );
+    const SMDS_MeshNode*  curN = face->GetNode( i );
     if ( nextN->GetPosition()->GetTypeOfPosition() == SMDS_TOP_FACE ||
-         curN->GetPosition()->GetTypeOfPosition() == SMDS_TOP_FACE )
+         curN-> GetPosition()->GetTypeOfPosition() == SMDS_TOP_FACE )
     {
-      double dist = SMESH_TNodeXYZ( face->GetNode(i)).Distance( nextN );
+      double dist = SMESH_TNodeXYZ( curN ).Distance( nextN );
       if ( dist < minSize )
         minSize = dist, iN = i;
     }
   }
-  double newStep = 0.8 * minSize / cosin;
+  double newStep = 0.8 * minSize / maxCosinEdge->_lenFactor;
   if ( newStep < data._stepSize )
   {
     data._stepSize = newStep;
-    data._stepSizeCoeff = 0.8 / cosin;
+    data._stepSizeCoeff = 0.8 / maxCosinEdge->_lenFactor;
     data._stepSizeNodes[0] = face->GetNode( iN );
     data._stepSizeNodes[1] = face->GetNode( SMESH_MesherHelper::WrapIndex( iN+1, nbNodes ));
   }
@@ -1898,21 +1906,24 @@ bool _ViscousBuilder::sortEdges( _SolidData&                    data,
         TGeomID iV = getMeshDS()->ShapeToIndex( vIt.Value() );
         vector<_LayerEdge*>& eV = edgesByGeom[ iV ];
         if ( eV.empty() ) continue;
-        double cosin = eV[0]->_cosin;
-        bool badCosin =
-          ( !eV[0]->_sWOL.IsNull() && ( eV[0]->_sWOL.ShapeType() == TopAbs_EDGE || !isShrinkEdge));
-        if ( badCosin )
-        {
-          gp_Vec dir1, dir2;
-          if ( eV[0]->_sWOL.ShapeType() == TopAbs_EDGE )
-            dir1 = getEdgeDir( TopoDS::Edge( eV[0]->_sWOL ), TopoDS::Vertex( vIt.Value() ));
-          else
-            dir1 = getFaceDir( TopoDS::Face( eV[0]->_sWOL ), TopoDS::Vertex( vIt.Value() ),
-                               eV[0]->_nodes[0], helper, ok);
-          dir2 = getEdgeDir( TopoDS::Edge( S ), TopoDS::Vertex( vIt.Value() ));
-          double angle = dir1.Angle( dir2 );
-          cosin = cos( angle );
-        }
+        // double cosin = eV[0]->_cosin;
+        // bool badCosin =
+        //   ( !eV[0]->_sWOL.IsNull() && ( eV[0]->_sWOL.ShapeType() == TopAbs_EDGE || !isShrinkEdge));
+        // if ( badCosin )
+        // {
+        //   gp_Vec dir1, dir2;
+        //   if ( eV[0]->_sWOL.ShapeType() == TopAbs_EDGE )
+        //     dir1 = getEdgeDir( TopoDS::Edge( eV[0]->_sWOL ), TopoDS::Vertex( vIt.Value() ));
+        //   else
+        //     dir1 = getFaceDir( TopoDS::Face( eV[0]->_sWOL ), TopoDS::Vertex( vIt.Value() ),
+        //                        eV[0]->_nodes[0], helper, ok);
+        //   dir2 = getEdgeDir( TopoDS::Edge( S ), TopoDS::Vertex( vIt.Value() ));
+        //   double angle = dir1.Angle( dir2 );
+        //   cosin = cos( angle );
+        // }
+        gp_Vec eDir = getEdgeDir( TopoDS::Edge( S ), TopoDS::Vertex( vIt.Value() ));
+        double angle = eDir.Angle( eV[0]->_normal );
+        double cosin = Cos( angle );
         needSmooth = ( cosin > theMinSmoothCosin );
       }
       break;
@@ -2004,8 +2015,8 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge&         edge,
   const SMDS_MeshNode* node = edge._nodes[0]; // source node
   SMDS_TypeOfPosition posType = node->GetPosition()->GetTypeOfPosition();
 
-  edge._len = 0;
-  edge._2neibors = 0;
+  edge._len       = 0;
+  edge._2neibors  = 0;
   edge._curvature = 0;
 
   // --------------------------
@@ -2016,18 +2027,16 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge&         edge,
   edge._normal.SetCoord(0,0,0);
 
   int totalNbFaces = 0;
-  gp_Pnt p;
-  gp_Vec du, dv, geomNorm;
+  gp_Vec geomNorm;
   bool normOK = true;
 
-  TGeomID shapeInd = node->getshapeId();
+  const TGeomID shapeInd = node->getshapeId();
   map< TGeomID, TopoDS_Shape >::const_iterator s2s = data._shrinkShape2Shape.find( shapeInd );
-  bool onShrinkShape ( s2s != data._shrinkShape2Shape.end() );
-  TopoDS_Shape vertEdge;
+  const bool onShrinkShape ( s2s != data._shrinkShape2Shape.end() );
 
   if ( onShrinkShape ) // one of faces the node is on has no layers
   {
-    vertEdge = getMeshDS()->IndexToShape( s2s->first ); // vertex or edge
+    TopoDS_Shape vertEdge = getMeshDS()->IndexToShape( s2s->first ); // vertex or edge
     if ( s2s->second.ShapeType() == TopAbs_EDGE )
     {
       // inflate from VERTEX along EDGE
@@ -2275,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 );
 
@@ -2363,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
@@ -2822,11 +2836,15 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data,
         int oldBadNb = badNb;
         badNb = 0;
         moved = false;
-        for ( int i = iBeg; i < iEnd; ++i )
-          moved |= data._edges[i]->Smooth(badNb);
+        if ( step % 2 )
+          for ( int i = iBeg; i < iEnd; ++i )
+            moved |= data._edges[i]->Smooth(badNb);
+        else
+          for ( int i = iEnd-1; i >= iBeg; --i )
+            moved |= data._edges[i]->Smooth(badNb);
         improved = ( badNb < oldBadNb );
 
-        // issue 22576. no bad faces but still there are intersections to fix
+        // issue 22576 -- no bad faces but still there are intersections to fix
         if ( improved && badNb == 0 )
           stepLimit = step + 3;
 
@@ -2881,6 +2899,8 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data,
   int iLE = 0;
   for ( size_t i = 0; i < data._edges.size(); ++i )
   {
+    if ( !data._edges[i]->_sWOL.IsNull() )
+      continue;
     if ( data._edges[i]->FindIntersection( *searcher, dist, data._epsilon, &intFace ))
       return false;
     if ( distToIntersection > dist )
@@ -3084,7 +3104,7 @@ bool _SolidData::GetShapeEdges(const TGeomID shapeID,
  */
 //================================================================================
 
-void _SolidData::AddFacesToSmooth( const set< TGeomID >& faceIDs )
+void _SolidData::AddShapesToSmooth( const set< TGeomID >& faceIDs )
 {
   // convert faceIDs to indices in _endEdgeOnShape
   set< size_t > iEnds;
@@ -3097,7 +3117,7 @@ void _SolidData::AddFacesToSmooth( const set< TGeomID >& faceIDs )
   set< size_t >::iterator endsIt = iEnds.begin();
 
   // "add" by move of _nbShapesToSmooth
-  int nbFacesToAdd = faceIDs.size();
+  int nbFacesToAdd = iEnds.size();
   while ( endsIt != iEnds.end() && *endsIt == _nbShapesToSmooth )
   {
     ++endsIt;
@@ -3344,7 +3364,9 @@ bool _ViscousBuilder::updateNormals( _SolidData&         data,
   for ( size_t i = 0; i < data._edges.size(); ++i )
   {
     _LayerEdge* edge = data._edges[i];
-    if ( !edge->IsOnEdge() || !edge->_sWOL.IsNull() ) continue;
+    if (( !edge->IsOnEdge() ) &&
+        ( edge->_sWOL.IsNull() || edge->_sWOL.ShapeType() != TopAbs_FACE ))
+      continue;
     if ( edge->FindIntersection( *searcher, dist, eps, &face ))
     {
       const _TmpMeshFaceOnEdge* f = (const _TmpMeshFaceOnEdge*) face;
@@ -3364,107 +3386,163 @@ bool _ViscousBuilder::updateNormals( _SolidData&         data,
   {
     dumpFunction(SMESH_Comment("updateNormals")<<data._index);
 
+    set< TGeomID > shapesToSmooth;
+
+    // vector to store new _normal and _cosin for each edge in edge2CloseEdge
+    vector< pair< _LayerEdge*, _LayerEdge > > edge2newEdge( edge2CloseEdge.size() );
+
     TLEdge2LEdgeSet::iterator e2ee = edge2CloseEdge.begin();
-    for ( ; e2ee != edge2CloseEdge.end(); ++e2ee )
+    for ( size_t iE = 0; e2ee != edge2CloseEdge.end(); ++e2ee, ++iE )
     {
-      _LayerEdge* edge1       = e2ee->first;
-      _LayerEdge* edge2       = 0;
-      set< _LayerEdge*, _LayerEdgeCmp >& ee  = e2ee->second;
+      _LayerEdge* edge1 = e2ee->first;
+      _LayerEdge* edge2 = 0;
+      set< _LayerEdge*, _LayerEdgeCmp >& ee = e2ee->second;
+
+      edge2newEdge[ iE ].first = NULL;
 
       // find EDGEs the edges reside
-      TopoDS_Edge E1, E2;
-      TopoDS_Shape S = helper.GetSubShapeByNode( edge1->_nodes[0], getMeshDS() );
-      if ( S.ShapeType() != TopAbs_EDGE )
-        continue; // TODO: find EDGE by VERTEX
-      E1 = TopoDS::Edge( S );
+      // TopoDS_Edge E1, E2;
+      // TopoDS_Shape S = helper.GetSubShapeByNode( edge1->_nodes[0], getMeshDS() );
+      // if ( S.ShapeType() != TopAbs_EDGE )
+      //   continue; // TODO: find EDGE by VERTEX
+      // E1 = TopoDS::Edge( S );
       set< _LayerEdge*, _LayerEdgeCmp >::iterator eIt = ee.begin();
-      while ( E2.IsNull() && eIt != ee.end())
+      for ( ; !edge2 && eIt != ee.end(); ++eIt )
       {
-        _LayerEdge* e2 = *eIt++;
-        TopoDS_Shape S = helper.GetSubShapeByNode( e2->_nodes[0], getMeshDS() );
-        if ( S.ShapeType() == TopAbs_EDGE )
-          E2 = TopoDS::Edge( S ), edge2 = e2;
+        if ( edge1->_sWOL == (*eIt)->_sWOL )
+          edge2 = *eIt;
       }
-      if ( E2.IsNull() ) continue; // TODO: find EDGE by VERTEX
+      if ( !edge2 ) continue;
+
+      edge2newEdge[ iE ].first = edge1;
+      _LayerEdge& newEdge = edge2newEdge[ iE ].second;
+      // while ( E2.IsNull() && eIt != ee.end())
+      // {
+      //   _LayerEdge* e2 = *eIt++;
+      //   TopoDS_Shape S = helper.GetSubShapeByNode( e2->_nodes[0], getMeshDS() );
+      //   if ( S.ShapeType() == TopAbs_EDGE )
+      //     E2 = TopoDS::Edge( S ), edge2 = e2;
+      // }
+      // if ( E2.IsNull() ) continue; // TODO: find EDGE by VERTEX
 
       // find 3 FACEs sharing 2 EDGEs
 
-      TopoDS_Face FF1[2], FF2[2];
-      PShapeIteratorPtr fIt = helper.GetAncestors(E1, *_mesh, TopAbs_FACE);
-      while ( fIt->more() && FF1[1].IsNull())
+      // TopoDS_Face FF1[2], FF2[2];
+      // PShapeIteratorPtr fIt = helper.GetAncestors(E1, *_mesh, TopAbs_FACE);
+      // while ( fIt->more() && FF1[1].IsNull() )
+      // {
+      //   const TopoDS_Face *F = (const TopoDS_Face*) fIt->next();
+      //   if ( helper.IsSubShape( *F, data._solid))
+      //     FF1[ FF1[0].IsNull() ? 0 : 1 ] = *F;
+      // }
+      // fIt = helper.GetAncestors(E2, *_mesh, TopAbs_FACE);
+      // while ( fIt->more() && FF2[1].IsNull())
+      // {
+      //   const TopoDS_Face *F = (const TopoDS_Face*) fIt->next();
+      //   if ( helper.IsSubShape( *F, data._solid))
+      //     FF2[ FF2[0].IsNull() ? 0 : 1 ] = *F;
+      // }
+      // // exclude a FACE common to E1 and E2 (put it to FFn[1] )
+      // if ( FF1[0].IsSame( FF2[0]) || FF1[0].IsSame( FF2[1]))
+      //   std::swap( FF1[0], FF1[1] );
+      // if ( FF2[0].IsSame( FF1[0]) )
+      //   std::swap( FF2[0], FF2[1] );
+      // if ( FF1[0].IsNull() || FF2[0].IsNull() )
+      //   continue;
+
+      // get a new normal for edge1
+      //bool ok;
+      gp_Vec dir1 = edge1->_normal, dir2 = edge2->_normal;
+      // if ( edge1->_cosin < 0 )
+      //   dir1 = getFaceDir( FF1[0], E1, edge1->_nodes[0], helper, ok ).Normalized();
+      // if ( edge2->_cosin < 0 )
+      //   dir2 = getFaceDir( FF2[0], E2, edge2->_nodes[0], helper, ok ).Normalized();
+
+      double cos1 = Abs( edge1->_cosin ), cos2 = Abs( edge2->_cosin );
+      double wgt1 = ( cos1 + 0.001 ) / ( cos1 + cos2 + 0.002 );
+      double wgt2 = ( cos2 + 0.001 ) / ( cos1 + cos2 + 0.002 );
+      newEdge._normal = ( wgt1 * dir1 + wgt2 * dir2 ).XYZ();
+      newEdge._normal.Normalize();
+
+      // cout << edge1->_nodes[0]->GetID() << " "
+      //      << edge2->_nodes[0]->GetID() << " NORM: "
+      //      << newEdge._normal.X() << ", " << newEdge._normal.Y() << ", " << newEdge._normal.Z() << endl;
+
+      // get new cosin
+      if ( cos1 < theMinSmoothCosin )
       {
-        const TopoDS_Face *F = (const TopoDS_Face*) fIt->next();
-        if ( helper.IsSubShape( *F, data._solid))
-          FF1[ FF1[0].IsNull() ? 0 : 1 ] = *F;
+        newEdge._cosin = edge2->_cosin;
       }
-      fIt = helper.GetAncestors(E2, *_mesh, TopAbs_FACE);
-      while ( fIt->more() && FF2[1].IsNull())
+      else if ( cos2 > theMinSmoothCosin ) // both cos1 and cos2 > theMinSmoothCosin
       {
-        const TopoDS_Face *F = (const TopoDS_Face*) fIt->next();
-        if ( helper.IsSubShape( *F, data._solid))
-          FF2[ FF2[0].IsNull() ? 0 : 1 ] = *F;
+        // gp_Vec dirInFace;
+        // if ( edge1->_cosin < 0 )
+        //   dirInFace = dir1;
+        // else
+        //   dirInFace = getFaceDir( FF1[0], E1, edge1->_nodes[0], helper, ok );
+        // double angle = dirInFace.Angle( edge1->_normal ); // [0,PI]
+        // edge1->SetCosin( Cos( angle ));
+        //newEdge._cosin = 0; // ???????????
+        newEdge._cosin = ( wgt1 * cos1 + wgt2 * cos2 ) * edge1->_cosin / cos1;
       }
-      // exclude a FACE common to E1 and E2 (put it at [1] in FF* )
-      if ( FF1[0].IsSame( FF2[0]) || FF1[0].IsSame( FF2[1]))
-        std::swap( FF1[0], FF1[1] );
-      if ( FF2[0].IsSame( FF1[0]) )
-        std::swap( FF2[0], FF2[1] );
-      if ( FF1[0].IsNull() || FF2[0].IsNull() )
-        continue;
-
-      // get a new normal for edge1
-      bool ok;
-      gp_Vec dir1 = edge1->_normal, dir2 = edge2->_normal;
-      if ( edge1->_cosin < 0 )
-        dir1 = getFaceDir( FF1[0], E1, edge1->_nodes[0], helper, ok ).Normalized();
-      if ( edge2->_cosin < 0 )
-        dir2 = getFaceDir( FF2[0], E2, edge2->_nodes[0], helper, ok ).Normalized();
-      // gp_Vec dir1 = getFaceDir( FF1[0], E1, edge1->_nodes[0], helper, ok );
-      // gp_Vec dir2 = getFaceDir( FF2[0], E2, edge2->_nodes[0], helper, ok2 );
-      // double wgt1 = ( edge1->_cosin + 1 ) / ( edge1->_cosin + edge2->_cosin + 2 );
-      // double wgt2 = ( edge2->_cosin + 1 ) / ( edge1->_cosin + edge2->_cosin + 2 );
-      // gp_Vec newNorm = wgt1 * dir1 + wgt2 * dir2;
-      // newNorm.Normalize();
-
-      double wgt1 = ( edge1->_cosin + 1 ) / ( edge1->_cosin + edge2->_cosin + 2 );
-      double wgt2 = ( edge2->_cosin + 1 ) / ( edge1->_cosin + edge2->_cosin + 2 );
-      gp_Vec newNorm = wgt1 * dir1 + wgt2 * dir2;
-      newNorm.Normalize();
-
-      edge1->_normal = newNorm.XYZ();
-
-      // update data of edge1 depending on _normal
-      const SMDS_MeshNode *n1, *n2;
-      n1 = edge1->_2neibors->_edges[0]->_nodes[0];
-      n2 = edge1->_2neibors->_edges[1]->_nodes[0];
-      edge1->SetDataByNeighbors( n1, n2, helper );
-      gp_Vec dirInFace;
-      if ( edge1->_cosin < 0 )
-        dirInFace = dir1;
       else
-        dirInFace = getFaceDir( FF1[0], E1, edge1->_nodes[0], helper, ok );
-      double angle = dir1.Angle( edge1->_normal ); // [0,PI]
-      edge1->SetCosin( cos( angle ));
+      {
+        newEdge._cosin = edge1->_cosin;
+      }
 
-      // limit data._stepSize
-      if ( edge1->_cosin > theMinSmoothCosin )
+      // find shapes that need smoothing due to change of _normal
+      if ( edge1->_cosin  < theMinSmoothCosin &&
+           newEdge._cosin > theMinSmoothCosin )
       {
-        SMDS_ElemIteratorPtr fIt = edge1->_nodes[0]->GetInverseElementIterator(SMDSAbs_Face);
-        while ( fIt->more() )
-          limitStepSize( data, fIt->next(), edge1->_cosin );
+        if ( edge1->_sWOL.IsNull() )
+        {
+          SMDS_ElemIteratorPtr fIt = edge1->_nodes[0]->GetInverseElementIterator(SMDSAbs_Face);
+          while ( fIt->more() )
+            shapesToSmooth.insert( fIt->next()->getshapeId() );
+          //limitStepSize( data, fIt->next(), edge1->_cosin ); // too late
+        }
+        else // edge1 inflates along a FACE
+        {
+          TopoDS_Shape V = helper.GetSubShapeByNode( edge1->_nodes[0], getMeshDS() );
+          PShapeIteratorPtr eIt = helper.GetAncestors( V, *_mesh, TopAbs_EDGE );
+          while ( const TopoDS_Shape* E = eIt->next() )
+          {
+            if ( !helper.IsSubShape( *E, /*FACE=*/edge1->_sWOL ))
+              continue;
+            gp_Vec edgeDir = getEdgeDir( TopoDS::Edge( *E ), TopoDS::Vertex( V ));
+            double   angle = edgeDir.Angle( newEdge._normal ); // [0,PI]
+            if ( angle < M_PI / 2 )
+              shapesToSmooth.insert( getMeshDS()->ShapeToIndex( *E ));
+          }
+        }
       }
-      // set new XYZ of target node
+    }
+
+    data.AddShapesToSmooth( shapesToSmooth );
+
+    // Update data of edges depending on a new _normal
+
+    for ( size_t iE = 0; iE < edge2newEdge.size(); ++iE )
+    {
+      _LayerEdge*   edge1 = edge2newEdge[ iE ].first;
+      _LayerEdge& newEdge = edge2newEdge[ iE ].second;
+      if ( !edge1 ) continue;
+
+      edge1->_normal = newEdge._normal;
+      edge1->SetCosin( newEdge._cosin );
       edge1->InvalidateStep( 1 );
       edge1->_len = 0;
       edge1->SetNewLength( data._stepSize, helper );
-    }
+      if ( edge1->IsOnEdge() )
+      {
+        const SMDS_MeshNode * n1 = edge1->_2neibors->_edges[0]->_nodes[0];
+        const SMDS_MeshNode * n2 = edge1->_2neibors->_edges[1]->_nodes[0];
+        edge1->SetDataByNeighbors( n1, n2, helper );
+      }
 
-    // Update normals and other dependent data of not intersecting _LayerEdge's
-    // neighboring the intersecting ones
+      // Update normals and other dependent data of not intersecting _LayerEdge's
+      // neighboring the intersecting ones
 
-    for ( e2ee = edge2CloseEdge.begin(); e2ee != edge2CloseEdge.end(); ++e2ee )
-    {
-      _LayerEdge* edge1 = e2ee->first;
       if ( !edge1->_2neibors )
         continue;
       for ( int j = 0; j < 2; ++j ) // loop on 2 neighbors
@@ -3473,7 +3551,7 @@ bool _ViscousBuilder::updateNormals( _SolidData&         data,
         if ( edge2CloseEdge.count ( neighbor ))
           continue; // j-th neighbor is also intersected
         _LayerEdge* prevEdge = edge1;
-        const int nbSteps = 6;
+        const int nbSteps = 10;
         for ( int step = nbSteps; step; --step ) // step from edge1 in j-th direction
         {
           if ( !neighbor->_2neibors )
@@ -3867,7 +3945,7 @@ bool _ViscousBuilder::updateNormalsOfConvexFaces( _SolidData&         data,
         }
       }
     }
-    data.AddFacesToSmooth( adjFacesToSmooth );
+    data.AddShapesToSmooth( adjFacesToSmooth );
 
     dumpFunctionEnd();
 
@@ -4023,7 +4101,7 @@ bool _LayerEdge::FindIntersection( SMESH_ElementSearcher&   searcher,
   bool segmentIntersected = false;
   distance = Precision::Infinite();
   int iFace = -1; // intersected face
-  for ( size_t j = 0 ; j < suspectFaces.size() && !segmentIntersected; ++j )
+  for ( size_t j = 0 ; j < suspectFaces.size() /*&& !segmentIntersected*/; ++j )
   {
     const SMDS_MeshElement* face = suspectFaces[j];
     if ( face->GetNodeIndex( _nodes.back() ) >= 0 ||
@@ -4041,7 +4119,7 @@ bool _LayerEdge::FindIntersection( SMESH_ElementSearcher&   searcher,
     {
       const SMDS_MeshNode* tria[3];
       tria[0] = *nIt++;
-      tria[1] = *nIt++;;
+      tria[1] = *nIt++;
       for ( int n2 = 2; n2 < nbNodes && !intFound; ++n2 )
       {
         tria[2] = *nIt++;
@@ -4051,7 +4129,7 @@ bool _LayerEdge::FindIntersection( SMESH_ElementSearcher&   searcher,
     }
     if ( intFound )
     {
-      if ( dist < segLen*(1.01) && dist > -(_len-segLen) )
+      if ( dist < segLen*(1.01) && dist > -(_len*_lenFactor-segLen) )
         segmentIntersected = true;
       if ( distance > dist )
         distance = dist, iFace = j;
@@ -4160,9 +4238,9 @@ bool _LayerEdge::SegTriaInter( const gp_Ax1&        lastSegment,
   /* calculate distance from vert0 to ray origin */
   gp_XYZ tvec = orig - vert0;
 
-  if ( tvec * dir > EPSILON )
+  //if ( tvec * dir > EPSILON )
     // intersected face is at back side of the temporary face this _LayerEdge belongs to
-    return false;
+    //return false;
 
   gp_XYZ edge1 = vert1 - vert0;
   gp_XYZ edge2 = vert2 - vert0;
@@ -4174,26 +4252,29 @@ bool _LayerEdge::SegTriaInter( const gp_Ax1&        lastSegment,
   double det = edge1 * pvec;
 
   if (det > -EPSILON && det < EPSILON)
-    return 0;
+    return false;
   double inv_det = 1.0 / det;
 
   /* calculate U parameter and test bounds */
   double u = ( tvec * pvec ) * inv_det;
-  if (u < 0.0 || u > 1.0)
-    return 0;
+  //if (u < 0.0 || u > 1.0)
+  if (u < -EPSILON || u > 1.0 + EPSILON)
+    return false;
 
   /* prepare to test V parameter */
   gp_XYZ qvec = tvec ^ edge1;
 
   /* calculate V parameter and test bounds */
   double v = (dir * qvec) * inv_det;
-  if ( v < 0.0 || u + v > 1.0 )
-    return 0;
+  //if ( v < 0.0 || u + v > 1.0 )
+  if ( v < -EPSILON || u + v > 1.0 + EPSILON)
+    return false;
 
   /* calculate t, ray intersects triangle */
   t = (edge2 * qvec) * inv_det;
 
-  return true;
+  //return true;
+  return t > 0.;
 }
 
 //================================================================================
@@ -4282,16 +4363,27 @@ bool _LayerEdge::Smooth(int& badNb)
     newPos += SMESH_TNodeXYZ( _simplices[i]._nPrev );
   newPos /= _simplices.size();
 
+  const gp_XYZ& curPos ( _pos.back() );
+  const gp_Pnt  prevPos( _pos[ _pos.size()-2 ]);
   if ( _curvature )
-    newPos += _normal * _curvature->lenDelta( _len );
-
-  gp_Pnt prevPos( _pos[ _pos.size()-2 ]);
+  {
+    double delta  = _curvature->lenDelta( _len );
+    if ( delta > 0 )
+      newPos += _normal * delta;
+    else
+    {
+      double segLen = _normal * ( newPos - prevPos.XYZ() );
+      if ( segLen + delta > 0 )
+        newPos += _normal * delta;
+    }
+    // double segLenChange = _normal * ( curPos - newPos );
+    // newPos += 0.5 * _normal * segLenChange;
+  }
 
   // count quality metrics (orientation) of tetras around _tgtNode
   int nbOkBefore = 0;
-  SMESH_TNodeXYZ tgtXYZ( _nodes.back() );
   for ( size_t i = 0; i < _simplices.size(); ++i )
-    nbOkBefore += _simplices[i].IsForward( _nodes[0], &tgtXYZ );
+    nbOkBefore += _simplices[i].IsForward( _nodes[0], &curPos );
 
   int nbOkAfter = 0;
   for ( size_t i = 0; i < _simplices.size(); ++i )