Salome HOME
http://www.salome-platform.org/forum/forum_10/450300019
[modules/smesh.git] / src / StdMeshers / StdMeshers_ViscousLayers.cxx
index c875ca104a4f622be907664d77f9c5f2329f6597..b2f6855b01c065ec0245b23b6a550c94e1552f25 100644 (file)
@@ -1,9 +1,9 @@
-// Copyright (C) 2007-2012  CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2014  CEA/DEN, EDF R&D, OPEN CASCADE
 //
 // This library is free software; you can redistribute it and/or
 // modify it under the terms of the GNU Lesser General Public
 // License as published by the Free Software Foundation; either
-// version 2.1 of the License.
+// version 2.1 of the License, or (at your option) any later version.
 //
 // This library is distributed in the hope that it will be useful,
 // but WITHOUT ANY WARRANTY; without even the implied warranty of
 #include "SMESH_ControlsDef.hxx"
 #include "SMESH_Gen.hxx"
 #include "SMESH_Group.hxx"
+#include "SMESH_HypoFilter.hxx"
 #include "SMESH_Mesh.hxx"
+#include "SMESH_MeshAlgos.hxx"
 #include "SMESH_MesherHelper.hxx"
 #include "SMESH_ProxyMesh.hxx"
 #include "SMESH_subMesh.hxx"
 #include "SMESH_subMeshEventListener.hxx"
-
-#include "utilities.h"
+#include "StdMeshers_FaceSide.hxx"
 
 #include <BRepAdaptor_Curve2d.hxx>
+#include <BRepAdaptor_Surface.hxx>
+#include <BRepLProp_SLProps.hxx>
 #include <BRep_Tool.hxx>
 #include <Bnd_B2d.hxx>
 #include <Bnd_B3d.hxx>
 #include <Geom2d_Line.hxx>
 #include <Geom2d_TrimmedCurve.hxx>
 #include <GeomAdaptor_Curve.hxx>
+#include <GeomLib.hxx>
 #include <Geom_Circle.hxx>
 #include <Geom_Curve.hxx>
 #include <Geom_Line.hxx>
 #include <Geom_TrimmedCurve.hxx>
 #include <Precision.hxx>
 #include <Standard_ErrorHandler.hxx>
+#include <Standard_Failure.hxx>
 #include <TColStd_Array1OfReal.hxx>
 #include <TopExp.hxx>
 #include <TopExp_Explorer.hxx>
 #include <TopTools_IndexedMapOfShape.hxx>
+#include <TopTools_ListOfShape.hxx>
 #include <TopTools_MapOfShape.hxx>
 #include <TopoDS.hxx>
 #include <TopoDS_Edge.hxx>
@@ -82,7 +88,7 @@
 using namespace std;
 
 //================================================================================
-namespace VISCOUS
+namespace VISCOUS_3D
 {
   typedef int TGeomID;
 
@@ -121,13 +127,13 @@ namespace VISCOUS
    * \brief Listener of events of 3D sub-meshes computed with viscous layers.
    * It is used to clear an inferior dim sub-meshes modified by viscous layers
    */
-  class _SrinkShapeListener : SMESH_subMeshEventListener
+  class _ShrinkShapeListener : SMESH_subMeshEventListener
   {
-    _SrinkShapeListener()
+    _ShrinkShapeListener()
       : SMESH_subMeshEventListener(/*isDeletable=*/false,
-                                   "StdMeshers_ViscousLayers::_SrinkShapeListener") {}
-    static SMESH_subMeshEventListener* Get() { static _SrinkShapeListener l; return &l; }
+                                   "StdMeshers_ViscousLayers::_ShrinkShapeListener") {}
   public:
+    static SMESH_subMeshEventListener* Get() { static _ShrinkShapeListener l; return &l; }
     virtual void ProcessEvent(const int                       event,
                               const int                       eventType,
                               SMESH_subMesh*                  solidSM,
@@ -139,23 +145,6 @@ namespace VISCOUS
         SMESH_subMeshEventListener::ProcessEvent(event,eventType,solidSM,data,hyp);
       }
     }
-    static void ToClearSubMeshWithSolid( SMESH_subMesh*      sm,
-                                         const TopoDS_Shape& solid)
-    {
-      SMESH_subMesh* solidSM = sm->GetFather()->GetSubMesh( solid );
-      SMESH_subMeshEventListenerData* data = solidSM->GetEventListenerData( Get());
-      if ( data )
-      {
-        if ( find( data->mySubMeshes.begin(), data->mySubMeshes.end(), sm ) ==
-             data->mySubMeshes.end())
-          data->mySubMeshes.push_back( sm );
-      }
-      else
-      {
-        data = SMESH_subMeshEventListenerData::MakeData( /*dependent=*/sm );
-        sm->SetEventListener( Get(), data, /*whereToListenTo=*/solidSM );
-      }
-    }
   };
   //--------------------------------------------------------------------------------
   /*!
@@ -205,6 +194,32 @@ namespace VISCOUS
     }
   };
   
+  //================================================================================
+  /*!
+   * \brief sets a sub-mesh event listener to clear sub-meshes of sub-shapes of
+   * the main shape when sub-mesh of the main shape is cleared,
+   * for example to clear sub-meshes of FACEs when sub-mesh of a SOLID
+   * is cleared
+   */
+  //================================================================================
+
+  void ToClearSubWithMain( SMESH_subMesh* sub, const TopoDS_Shape& main)
+  {
+    SMESH_subMesh* mainSM = sub->GetFather()->GetSubMesh( main );
+    SMESH_subMeshEventListenerData* data =
+      mainSM->GetEventListenerData( _ShrinkShapeListener::Get());
+    if ( data )
+    {
+      if ( find( data->mySubMeshes.begin(), data->mySubMeshes.end(), sub ) ==
+           data->mySubMeshes.end())
+        data->mySubMeshes.push_back( sub );
+    }
+    else
+    {
+      data = SMESH_subMeshEventListenerData::MakeData( /*dependent=*/sub );
+      sub->SetEventListener( _ShrinkShapeListener::Get(), data, /*whereToListenTo=*/mainSM );
+    }
+  }
   //--------------------------------------------------------------------------------
   /*!
    * \brief Simplex (triangle or tetrahedron) based on 1 (tria) or 2 (tet) nodes of
@@ -215,8 +230,11 @@ namespace VISCOUS
   struct _Simplex
   {
     const SMDS_MeshNode *_nPrev, *_nNext; // nodes on a smoothed mesh surface
-    _Simplex(const SMDS_MeshNode* nPrev=0, const SMDS_MeshNode* nNext=0)
-      : _nPrev(nPrev), _nNext(nNext) {}
+    const SMDS_MeshNode *_nOpp; // in 2D case, a node opposite to a smoothed node in QUAD
+    _Simplex(const SMDS_MeshNode* nPrev=0,
+             const SMDS_MeshNode* nNext=0,
+             const SMDS_MeshNode* nOpp=0)
+      : _nPrev(nPrev), _nNext(nNext), _nOpp(nOpp) {}
     bool IsForward(const SMDS_MeshNode* nSrc, const gp_XYZ* pntTgt) const
     {
       const double M[3][3] =
@@ -256,6 +274,7 @@ namespace VISCOUS
   {
     double _r; // radius
     double _k; // factor to correct node smoothed position
+    double _h2lenRatio; // avgNormProj / (2*avgDist)
   public:
     static _Curvature* New( double avgNormProj, double avgDist )
     {
@@ -266,10 +285,12 @@ namespace VISCOUS
         c->_r = avgDist * avgDist / avgNormProj;
         c->_k = avgDist * avgDist / c->_r / c->_r;
         c->_k *= ( c->_r < 0 ? 1/1.1 : 1.1 ); // not to be too restrictive
+        c->_h2lenRatio = avgNormProj / ( avgDist + avgDist );
       }
       return c;
     }
     double lenDelta(double len) const { return _k * ( _r + len ); }
+    double lenDeltaByDist(double dist) const { return dist * _h2lenRatio; }
   };
   struct _LayerEdge;
   //--------------------------------------------------------------------------------
@@ -305,7 +326,7 @@ namespace VISCOUS
 
     gp_XYZ              _normal; // to solid surface
     vector<gp_XYZ>      _pos; // points computed during inflation
-    double              _len; // length achived with the last step
+    double              _len; // length achived with the last inflation step
     double              _cosin; // of angle (_normal ^ surface)
     double              _lenFactor; // to compute _len taking _cosin into account
 
@@ -344,7 +365,7 @@ namespace VISCOUS
                        const double&        epsilon) const;
     gp_Ax1 LastSegment(double& segLen) const;
     bool IsOnEdge() const { return _2neibors; }
-    void Copy( _LayerEdge& other, SMESH_MesherHelper& helper );
+    gp_XYZ Copy( _LayerEdge& other, SMESH_MesherHelper& helper );
     void SetCosin( double cosin );
   };
   struct _LayerEdgeCmp
@@ -367,21 +388,27 @@ namespace VISCOUS
   {
     TopoDS_Shape                    _solid;
     const StdMeshers_ViscousLayers* _hyp;
+    TopoDS_Shape                    _hypShape;
     _MeshOfSolid*                   _proxyMesh;
     set<TGeomID>                    _reversedFaceIds;
+    set<TGeomID>                    _ignoreFaceIds;
 
     double                          _stepSize, _stepSizeCoeff;
     const SMDS_MeshNode*            _stepSizeNodes[2];
 
     TNode2Edge                      _n2eMap;
+    // map to find _n2eMap of another _SolidData by a shrink shape shared by two _SolidData's
+    map< TGeomID, TNode2Edge* >     _s2neMap;
     // edges of _n2eMap. We keep same data in two containers because
     // iteration over the map is 5 time longer than over the vector
     vector< _LayerEdge* >           _edges;
+    // edges on EDGE's with null _sWOL, whose _simplices are used to stop inflation
+    vector< _LayerEdge* >           _simplexTestEdges;
 
-    // key: an id of shape (EDGE or VERTEX) shared by a FACE with
-    // layers and a FACE w/o layers
+    // key:   an id of shape (EDGE or VERTEX) shared by a FACE with
+    //        layers and a FACE w/o layers
     // value: the shape (FACE or EDGE) to shrink mesh on.
-    // _LayerEdge's basing on nodes on key shape are inflated along the value shape
+    //       _LayerEdge's basing on nodes on key shape are inflated along the value shape
     map< TGeomID, TopoDS_Shape >     _shrinkShape2Shape;
 
     // FACE's WOL, srink on which is forbiden due to algo on the adjacent SOLID
@@ -399,7 +426,9 @@ namespace VISCOUS
 
     _SolidData(const TopoDS_Shape&             s=TopoDS_Shape(),
                const StdMeshers_ViscousLayers* h=0,
-               _MeshOfSolid*                   m=0) :_solid(s), _hyp(h), _proxyMesh(m) {}
+               const TopoDS_Shape&             hs=TopoDS_Shape(),
+               _MeshOfSolid*                   m=0)
+      :_solid(s), _hyp(h), _hypShape(hs), _proxyMesh(m) {}
     ~_SolidData();
 
     Handle(Geom_Curve) CurveForSmooth( const TopoDS_Edge&    E,
@@ -419,12 +448,18 @@ namespace VISCOUS
     //vector<const SMDS_MeshNode*> _nodesAround;
     vector<_Simplex>             _simplices; // for quality check
 
+    enum SmoothType { LAPLACIAN, CENTROIDAL, ANGULAR, TFI };
+
     bool Smooth(int&                  badNb,
                 Handle(Geom_Surface)& surface,
                 SMESH_MesherHelper&   helper,
                 const double          refSign,
-                bool                  isCentroidal,
+                SmoothType            how,
                 bool                  set3D);
+
+    gp_XY computeAngularPos(vector<gp_XY>& uv,
+                            const gp_XY&   uvToFix,
+                            const double   refSign );
   };
   //--------------------------------------------------------------------------------
   /*!
@@ -451,6 +486,9 @@ namespace VISCOUS
     bool makeLayer(_SolidData& data);
     bool setEdgeData(_LayerEdge& edge, const set<TGeomID>& subIds,
                      SMESH_MesherHelper& helper, _SolidData& data);
+    gp_XYZ getWeigthedNormal( const SMDS_MeshNode*         n,
+                              std::pair< TGeomID, gp_XYZ > fId2Normal[],
+                              const int                    nbFaces );
     bool findNeiborsOnEdge(const _LayerEdge*     edge,
                            const SMDS_MeshNode*& n1,
                            const SMDS_MeshNode*& n2,
@@ -459,6 +497,8 @@ namespace VISCOUS
                        const set<TGeomID>& ingnoreShapes,
                        const _SolidData*   dataToCheckOri = 0,
                        const bool          toSort = false);
+    void findSimplexTestEdges( _SolidData&                    data,
+                               vector< vector<_LayerEdge*> >& edgesByGeom);
     bool sortEdges( _SolidData&                    data,
                     vector< vector<_LayerEdge*> >& edgesByGeom);
     void limitStepSize( _SolidData&             data,
@@ -479,7 +519,11 @@ namespace VISCOUS
     bool prepareEdgeToShrink( _LayerEdge& edge, const TopoDS_Face& F,
                               SMESH_MesherHelper& helper,
                               const SMESHDS_SubMesh* faceSubMesh );
-    void fixBadFaces(const TopoDS_Face& F, SMESH_MesherHelper& helper);
+    void fixBadFaces(const TopoDS_Face&          F,
+                     SMESH_MesherHelper&         helper,
+                     const bool                  is2D,
+                     const int                   step,
+                     set<const SMDS_MeshNode*> * involvedNodes=NULL);
     bool addBoundaryElements();
 
     bool error( const string& text, int solidID=-1 );
@@ -492,7 +536,6 @@ namespace VISCOUS
     SMESH_ComputeErrorPtr _error;
 
     vector< _SolidData >  _sdVec;
-    set<TGeomID>          _ignoreShapeIds;
     int                   _tmpFaceID;
   };
   //--------------------------------------------------------------------------------
@@ -528,7 +571,7 @@ namespace VISCOUS
     virtual vtkIdType GetVtkType() const                      { return -1; }
     virtual SMDSAbs_EntityType   GetEntityType() const        { return SMDSEntity_Last; }
     virtual SMDSAbs_GeometryType GetGeomType() const          { return SMDSGeom_TRIANGLE; }
-virtual SMDS_ElemIteratorPtr elementsIterator(SMDSAbs_ElementType type) const
+    virtual SMDS_ElemIteratorPtr elementsIterator(SMDSAbs_ElementType) const
     { return SMDS_ElemIteratorPtr( new SMDS_NodeVectorElemIterator( _nn.begin(), _nn.end()));}
   };
   //--------------------------------------------------------------------------------
@@ -547,22 +590,62 @@ virtual SMDS_ElemIteratorPtr elementsIterator(SMDSAbs_ElementType type) const
       _nn[3]=_le2->_nodes[0];
     }
   };
-} // namespace VISCOUS
+  //--------------------------------------------------------------------------------
+  /*!
+   * \brief Retriever of node coordinates either directly of from a surface by node UV.
+   * \warning Location of a surface is ignored
+   */
+  struct NodeCoordHelper
+  {
+    SMESH_MesherHelper&        _helper;
+    const TopoDS_Face&         _face;
+    Handle(Geom_Surface)       _surface;
+    gp_XYZ (NodeCoordHelper::* _fun)(const SMDS_MeshNode* n) const;
+
+    NodeCoordHelper(const TopoDS_Face& F, SMESH_MesherHelper& helper, bool is2D)
+      : _helper( helper ), _face( F )
+    {
+      if ( is2D )
+      {
+        TopLoc_Location loc;
+        _surface = BRep_Tool::Surface( _face, loc );
+      }
+      if ( _surface.IsNull() )
+        _fun = & NodeCoordHelper::direct;
+      else
+        _fun = & NodeCoordHelper::byUV;
+    }
+    gp_XYZ operator()(const SMDS_MeshNode* n) const { return (this->*_fun)( n ); }
+
+  private:
+    gp_XYZ direct(const SMDS_MeshNode* n) const
+    {
+      return SMESH_TNodeXYZ( n );
+    }
+    gp_XYZ byUV  (const SMDS_MeshNode* n) const
+    {
+      gp_XY uv = _helper.GetNodeUV( _face, n );
+      return _surface->Value( uv.X(), uv.Y() ).XYZ();
+    }
+  };
+} // namespace VISCOUS_3D
 
 //================================================================================
 // StdMeshers_ViscousLayers hypothesis
 //
 StdMeshers_ViscousLayers::StdMeshers_ViscousLayers(int hypId, int studyId, SMESH_Gen* gen)
   :SMESH_Hypothesis(hypId, studyId, gen),
-   _nbLayers(1), _thickness(1), _stretchFactor(1)
+   _isToIgnoreShapes(1), _nbLayers(1), _thickness(1), _stretchFactor(1)
 {
   _name = StdMeshers_ViscousLayers::GetHypType();
   _param_algo_dim = -3; // auxiliary hyp used by 3D algos
 } // --------------------------------------------------------------------------------
-void StdMeshers_ViscousLayers::SetIgnoreFaces(const std::vector<int>& faceIds)
+void StdMeshers_ViscousLayers::SetBndShapes(const std::vector<int>& faceIds, bool toIgnore)
 {
-  if ( faceIds != _ignoreFaceIds )
-    _ignoreFaceIds = faceIds, NotifySubMeshesHypothesisModification();
+  if ( faceIds != _shapeIds )
+    _shapeIds = faceIds, NotifySubMeshesHypothesisModification();
+  if ( _isToIgnoreShapes != toIgnore )
+    _isToIgnoreShapes = toIgnore, NotifySubMeshesHypothesisModification();
 } // --------------------------------------------------------------------------------
 void StdMeshers_ViscousLayers::SetTotalThickness(double thickness)
 {
@@ -584,7 +667,7 @@ StdMeshers_ViscousLayers::Compute(SMESH_Mesh&         theMesh,
                                   const TopoDS_Shape& theShape,
                                   const bool          toMakeN2NMap) const
 {
-  using namespace VISCOUS;
+  using namespace VISCOUS_3D;
   _ViscousBuilder bulder;
   SMESH_ComputeErrorPtr err = bulder.Compute( theMesh, theShape );
   if ( err && !err->IsOK() )
@@ -620,17 +703,22 @@ std::ostream & StdMeshers_ViscousLayers::SaveTo(std::ostream & save)
   save << " " << _nbLayers
        << " " << _thickness
        << " " << _stretchFactor
-       << " " << _ignoreFaceIds.size();
-  for ( unsigned i = 0; i < _ignoreFaceIds.size(); ++i )
-    save << " " << _ignoreFaceIds[i];
+       << " " << _shapeIds.size();
+  for ( size_t i = 0; i < _shapeIds.size(); ++i )
+    save << " " << _shapeIds[i];
+  save << " " << !_isToIgnoreShapes; // negate to keep the behavior in old studies.
   return save;
 } // --------------------------------------------------------------------------------
 std::istream & StdMeshers_ViscousLayers::LoadFrom(std::istream & load)
 {
-  int nbFaces, faceID;
+  int nbFaces, faceID, shapeToTreat;
   load >> _nbLayers >> _thickness >> _stretchFactor >> nbFaces;
-  while ( _ignoreFaceIds.size() < nbFaces && load >> faceID )
-    _ignoreFaceIds.push_back( faceID );
+  while ( _shapeIds.size() < nbFaces && load >> faceID )
+    _shapeIds.push_back( faceID );
+  if ( load >> shapeToTreat )
+    _isToIgnoreShapes = !shapeToTreat;
+  else
+    _isToIgnoreShapes = true; // old behavior
   return load;
 } // --------------------------------------------------------------------------------
 bool StdMeshers_ViscousLayers::SetParametersByMesh(const SMESH_Mesh*   theMesh,
@@ -716,27 +804,36 @@ namespace
     gp_XYZ dir(0,0,0);
     if ( !( ok = ( edges.size() > 0 ))) return dir;
     // get average dir of edges going fromV
-    gp_Vec edgeDir;
-    for ( unsigned i = 0; i < edges.size(); ++i )
+    gp_XYZ edgeDir;
+    //if ( edges.size() > 1 )
+    for ( size_t i = 0; i < edges.size(); ++i )
     {
       edgeDir = getEdgeDir( edges[i], fromV );
-      double size2 = edgeDir.SquareMagnitude();
+      double size2 = edgeDir.SquareModulus();
       if ( size2 > numeric_limits<double>::min() )
         edgeDir /= sqrt( size2 );
       else
         ok = false;
-      dir += edgeDir.XYZ();
+      dir += edgeDir;
     }
     gp_XYZ fromEdgeDir = getFaceDir( F, edges[0], node, helper, ok );
-    if ( edges.size() == 1 || dir.SquareModulus() < 1e-10)
-      dir = fromEdgeDir;
-    else if ( dir * fromEdgeDir < 0 )
-      dir *= -1;
+    try {
+      if ( edges.size() == 1 )
+        dir = fromEdgeDir;
+      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 )
+    {
+      ok = false;
+    }
     if ( ok )
     {
       //dir /= edges.size();
       if ( cosin ) {
-        double angle = edgeDir.Angle( dir );
+        double angle = gp_Vec( edgeDir ).Angle( dir );
         *cosin = cos( angle );
       }
     }
@@ -750,13 +847,15 @@ namespace
 
   bool isConcave( const TopoDS_Face& F, SMESH_MesherHelper& helper )
   {
+    // if ( helper.Count( F, TopAbs_WIRE, /*useMap=*/false) > 1 )
+    //   return true;
     gp_Vec2d drv1, drv2;
     gp_Pnt2d p;
     TopExp_Explorer eExp( F.Oriented( TopAbs_FORWARD ), TopAbs_EDGE );
     for ( ; eExp.More(); eExp.Next() )
     {
       const TopoDS_Edge& E = TopoDS::Edge( eExp.Current() );
-      if ( BRep_Tool::Degenerated( E )) continue;
+      if ( SMESH_Algo::isDegenerated( E )) continue;
       // check if 2D curve is concave
       BRepAdaptor_Curve2d curve( E, F );
       const int nbIntervals = curve.NbIntervals( GeomAbs_C2 );
@@ -771,7 +870,7 @@ namespace
         double cross = drv2 ^ drv1;
         if ( E.Orientation() == TopAbs_REVERSED )
           cross = -cross;
-        isConvex = ( cross 1e-9 );
+        isConvex = ( cross > 0.1 ); //-1e-9 );
       }
       // check if concavity is strong enough to care about it
       //const double maxAngle = 5 * Standard_PI180;
@@ -801,6 +900,26 @@ namespace
         // }
       }
     }
+    // check angles at VERTEXes
+    TError error;
+    TSideVector wires = StdMeshers_FaceSide::GetFaceWires( F, *helper.GetMesh(), 0, error );
+    for ( size_t iW = 0; iW < wires.size(); ++iW )
+    {
+      const int nbEdges = wires[iW]->NbEdges();
+      if ( nbEdges < 2 && SMESH_Algo::isDegenerated( wires[iW]->Edge(0)))
+        continue;
+      for ( int iE1 = 0; iE1 < nbEdges; ++iE1 )
+      {
+        if ( SMESH_Algo::isDegenerated( wires[iW]->Edge( iE1 ))) continue;
+        int iE2 = ( iE1 + 1 ) % nbEdges;
+        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 );
+        if ( angle < -5. * M_PI / 180. )
+          return true;
+      }
+    }
     return false;
   }
   //--------------------------------------------------------------------------------
@@ -812,13 +931,15 @@ namespace
       const char* fname = "/tmp/viscous.py";
       cout << "execfile('"<<fname<<"')"<<endl;
       py = new ofstream(fname);
-      *py << "from smesh import *" << endl
-          << "meshSO = GetCurrentStudy().FindObjectID('0:1:2:3')" << endl
-          << "mesh = Mesh( meshSO.GetObject() )"<<endl;
+      *py << "import SMESH" << endl
+          << "from salome.smesh import smeshBuilder" << endl
+          << "smesh  = smeshBuilder.New(salome.myStudy)" << endl
+          << "meshSO = smesh.GetCurrentStudy().FindObjectID('0:1:2:3')" << endl
+          << "mesh   = smesh.Mesh( meshSO.GetObject() )"<<endl;
     }
     void Finish() {
       if (py)
-        *py << "mesh.MakeGroup('Viscous Prisms',VOLUME,FT_ElemGeomType,'=',Geom_PENTA)"<<endl;
+        *py << "mesh.MakeGroup('Viscous Prisms',SMESH.VOLUME,SMESH.FT_ElemGeomType,'=',SMESH.Geom_PENTA)"<<endl;
       delete py; py=0;
     }
     ~PyDump() { Finish(); }
@@ -849,7 +970,7 @@ namespace
 #endif
 }
 
-using namespace VISCOUS;
+using namespace VISCOUS_3D;
 
 //================================================================================
 /*!
@@ -972,7 +1093,7 @@ SMESH_ComputeErrorPtr _ViscousBuilder::Compute(SMESH_Mesh&         theMesh,
   if ( !findFacesWithLayers() )
     return _error;
 
-  for ( unsigned i = 0; i < _sdVec.size(); ++i )
+  for ( size_t i = 0; i < _sdVec.size(); ++i )
   {
     if ( ! makeLayer(_sdVec[i]) )
       return _error;
@@ -1011,6 +1132,7 @@ bool _ViscousBuilder::findSolidsWithLayers()
   _sdVec.reserve( allSolids.Extent());
 
   SMESH_Gen* gen = _mesh->GetGen();
+  SMESH_HypoFilter filter;
   for ( int i = 1; i <= allSolids.Extent(); ++i )
   {
     // find StdMeshers_ViscousLayers hyp assigned to the i-th solid
@@ -1025,10 +1147,14 @@ bool _ViscousBuilder::findSolidsWithLayers()
       viscHyp = dynamic_cast<const StdMeshers_ViscousLayers*>( *hyp );
     if ( viscHyp )
     {
+      TopoDS_Shape hypShape;
+      filter.Init( filter.Is( viscHyp ));
+      _mesh->GetHypothesis( allSolids(i), filter, true, &hypShape );
+
       _MeshOfSolid* proxyMesh = _ViscousListener::GetSolidMesh( _mesh,
                                                                 allSolids(i),
                                                                 /*toCreate=*/true);
-      _sdVec.push_back( _SolidData( allSolids(i), viscHyp, proxyMesh ));
+      _sdVec.push_back( _SolidData( allSolids(i), viscHyp, hypShape, proxyMesh ));
       _sdVec.back()._index = getMeshDS()->ShapeToIndex( allSolids(i));
     }
   }
@@ -1047,44 +1173,69 @@ bool _ViscousBuilder::findSolidsWithLayers()
 
 bool _ViscousBuilder::findFacesWithLayers()
 {
+  SMESH_MesherHelper helper( *_mesh );
+  TopExp_Explorer exp;
+  TopTools_IndexedMapOfShape solids;
+
   // collect all faces to ignore defined by hyp
-  vector<TopoDS_Shape> ignoreFaces;
-  for ( unsigned i = 0; i < _sdVec.size(); ++i )
+  for ( size_t i = 0; i < _sdVec.size(); ++i )
   {
-    vector<TGeomID> ids = _sdVec[i]._hyp->GetIgnoreFaces();
-    for ( unsigned i = 0; i < ids.size(); ++i )
+    solids.Add( _sdVec[i]._solid );
+
+    vector<TGeomID> ids = _sdVec[i]._hyp->GetBndShapes();
+    if ( _sdVec[i]._hyp->IsToIgnoreShapes() ) // FACEs to ignore are given
     {
-      const TopoDS_Shape& s = getMeshDS()->IndexToShape( ids[i] );
-      if ( !s.IsNull() && s.ShapeType() == TopAbs_FACE )
+      for ( size_t ii = 0; ii < ids.size(); ++ii )
       {
-        _ignoreShapeIds.insert( ids[i] );
-        ignoreFaces.push_back( s );
+        const TopoDS_Shape& s = getMeshDS()->IndexToShape( ids[ii] );
+        if ( !s.IsNull() && s.ShapeType() == TopAbs_FACE )
+          _sdVec[i]._ignoreFaceIds.insert( ids[ii] );
+      }
+    }
+    else // FACEs with layers are given
+    {
+      exp.Init( _sdVec[i]._solid, TopAbs_FACE );
+      for ( ; exp.More(); exp.Next() )
+      {
+        TGeomID faceInd = getMeshDS()->ShapeToIndex( exp.Current() );
+        if ( find( ids.begin(), ids.end(), faceInd ) == ids.end() )
+          _sdVec[i]._ignoreFaceIds.insert( faceInd );
       }
     }
-  }
 
-  // ignore internal faces
-  SMESH_MesherHelper helper( *_mesh );
-  TopExp_Explorer exp;
-  for ( unsigned i = 0; i < _sdVec.size(); ++i )
-  {
-    exp.Init( _sdVec[i]._solid.Oriented( TopAbs_FORWARD ), TopAbs_FACE );
-    for ( ; exp.More(); exp.Next() )
-    {
-      TGeomID faceInd = getMeshDS()->ShapeToIndex( exp.Current() );
-      if ( helper.NbAncestors( exp.Current(), *_mesh, TopAbs_SOLID ) > 1 )
-      {     
-        _ignoreShapeIds.insert( faceInd );
-        ignoreFaces.push_back( exp.Current() );
-        if ( SMESH_Algo::IsReversedSubMesh( TopoDS::Face( exp.Current() ), getMeshDS()))
+    // ignore internal FACEs if inlets and outlets are specified
+    {
+      TopTools_IndexedDataMapOfShapeListOfShape solidsOfFace;
+      if ( _sdVec[i]._hyp->IsToIgnoreShapes() )
+        TopExp::MapShapesAndAncestors( _sdVec[i]._hypShape,
+                                       TopAbs_FACE, TopAbs_SOLID, solidsOfFace);
+
+      exp.Init( _sdVec[i]._solid.Oriented( TopAbs_FORWARD ), TopAbs_FACE );
+      for ( ; exp.More(); exp.Next() )
+      {
+        const TopoDS_Face& face = TopoDS::Face( exp.Current() );
+        if ( helper.NbAncestors( face, *_mesh, TopAbs_SOLID ) < 2 )
+          continue;
+
+        const TGeomID faceInd = getMeshDS()->ShapeToIndex( face );
+        if ( _sdVec[i]._hyp->IsToIgnoreShapes() )
+        {
+          int nbSolids = solidsOfFace.FindFromKey( face ).Extent();
+          if ( nbSolids > 1 )
+            _sdVec[i]._ignoreFaceIds.insert( faceInd );
+        }
+
+        if ( helper.IsReversedSubMesh( face ))
+        {
           _sdVec[i]._reversedFaceIds.insert( faceInd );
+        }
       }
     }
   }
 
   // Find faces to shrink mesh on (solution 2 in issue 0020832);
   TopTools_IndexedMapOfShape shapes;
-  for ( unsigned i = 0; i < _sdVec.size(); ++i )
+  for ( size_t i = 0; i < _sdVec.size(); ++i )
   {
     shapes.Clear();
     TopExp::MapShapes(_sdVec[i]._solid, TopAbs_EDGE, shapes);
@@ -1104,18 +1255,35 @@ bool _ViscousBuilder::findFacesWithLayers()
       // check presence of layers on them
       int ignore[2];
       for ( int j = 0; j < 2; ++j )
-        ignore[j] = _ignoreShapeIds.count ( getMeshDS()->ShapeToIndex( FF[j] ));
-      if ( ignore[0] == ignore[1] ) continue; // nothing interesting
+        ignore[j] = _sdVec[i]._ignoreFaceIds.count ( getMeshDS()->ShapeToIndex( FF[j] ));
+      if ( ignore[0] == ignore[1] )
+        continue; // nothing interesting
       TopoDS_Shape fWOL = FF[ ignore[0] ? 0 : 1 ];
+      // check presence of layers on fWOL within an adjacent SOLID
+      PShapeIteratorPtr sIt = helper.GetAncestors( fWOL, *_mesh, TopAbs_SOLID );
+      while ( const TopoDS_Shape* solid = sIt->next() )
+        if ( !solid->IsSame( _sdVec[i]._solid ))
+        {
+          int iSolid = solids.FindIndex( *solid );
+          int  iFace = getMeshDS()->ShapeToIndex( fWOL );
+          if ( iSolid > 0 && !_sdVec[ iSolid-1 ]._ignoreFaceIds.count( iFace ))
+          {
+            _sdVec[i]._noShrinkFaces.insert( iFace );
+            fWOL.Nullify();
+          }
+        }
       // add edge to maps
-      TGeomID edgeInd = getMeshDS()->ShapeToIndex( edge );
-      _sdVec[i]._shrinkShape2Shape.insert( make_pair( edgeInd, fWOL ));
+      if ( !fWOL.IsNull())
+      {
+        TGeomID edgeInd = getMeshDS()->ShapeToIndex( edge );
+        _sdVec[i]._shrinkShape2Shape.insert( make_pair( edgeInd, fWOL ));
+      }
     }
   }
   // Exclude from _shrinkShape2Shape FACE's that can't be shrinked since
   // the algo of the SOLID sharing the FACE does not support it
   set< string > notSupportAlgos; notSupportAlgos.insert("Hexa_3D");
-  for ( unsigned i = 0; i < _sdVec.size(); ++i )
+  for ( size_t i = 0; i < _sdVec.size(); ++i )
   {
     TopTools_MapOfShape noShrinkVertices;
     map< TGeomID, TopoDS_Shape >::iterator e2f = _sdVec[i]._shrinkShape2Shape.begin();
@@ -1132,7 +1300,7 @@ bool _ViscousBuilder::findFacesWithLayers()
         SMESH_Algo* algo = _mesh->GetGen()->GetAlgo( *_mesh, *solid );
         if ( !algo || !notSupportAlgos.count( algo->GetName() )) continue;
         notShrinkFace = true;
-        for ( unsigned j = 0; j < _sdVec.size(); ++j )
+        for ( size_t j = 0; j < _sdVec.size(); ++j )
         {
           if ( _sdVec[j]._solid.IsSame( *solid ) )
             if ( _sdVec[j]._shrinkShape2Shape.count( edgeID ))
@@ -1165,10 +1333,10 @@ bool _ViscousBuilder::findFacesWithLayers()
       }
     }
   }
-      
+
   // Find the SHAPE along which to inflate _LayerEdge based on VERTEX
 
-  for ( unsigned i = 0; i < _sdVec.size(); ++i )
+  for ( size_t i = 0; i < _sdVec.size(); ++i )
   {
     shapes.Clear();
     TopExp::MapShapes(_sdVec[i]._solid, TopAbs_VERTEX, shapes);
@@ -1182,11 +1350,12 @@ bool _ViscousBuilder::findFacesWithLayers()
       while ( fIt->more())
       {
         const TopoDS_Shape* f = fIt->next();
-        const int         fID = getMeshDS()->ShapeToIndex( *f );
         if ( helper.IsSubShape( *f, _sdVec[i]._solid ) )
         {
           totalNbFaces++;
-          if ( _ignoreShapeIds.count ( fID ) && ! _sdVec[i]._noShrinkFaces.count( fID ))
+          const int fID = getMeshDS()->ShapeToIndex( *f );
+          if ( _sdVec[i]._ignoreFaceIds.count ( fID ) &&
+               !_sdVec[i]._noShrinkFaces.count( fID ))
             facesWOL.push_back( *f );
         }
       }
@@ -1196,48 +1365,61 @@ bool _ViscousBuilder::findFacesWithLayers()
       switch ( facesWOL.size() )
       {
       case 1:
+      {
+        helper.SetSubShape( facesWOL[0] );
+        if ( helper.IsRealSeam( vInd )) // inflate along a seam edge?
         {
-          helper.SetSubShape( facesWOL[0] );
-          if ( helper.IsRealSeam( vInd )) // inflate along a seam edge?
+          TopoDS_Shape seamEdge;
+          PShapeIteratorPtr eIt = helper.GetAncestors(vertex, *_mesh, TopAbs_EDGE);
+          while ( eIt->more() && seamEdge.IsNull() )
           {
-            TopoDS_Shape seamEdge;
-            PShapeIteratorPtr eIt = helper.GetAncestors(vertex, *_mesh, TopAbs_EDGE);
-            while ( eIt->more() && seamEdge.IsNull() )
-            {
-              const TopoDS_Shape* e = eIt->next();
-              if ( helper.IsRealSeam( *e ) )
-                seamEdge = *e;
-            }
-            if ( !seamEdge.IsNull() )
-            {
-              _sdVec[i]._shrinkShape2Shape.insert( make_pair( vInd, seamEdge ));
-              break;
-            }
+            const TopoDS_Shape* e = eIt->next();
+            if ( helper.IsRealSeam( *e ) )
+              seamEdge = *e;
+          }
+          if ( !seamEdge.IsNull() )
+          {
+            _sdVec[i]._shrinkShape2Shape.insert( make_pair( vInd, seamEdge ));
+            break;
           }
-          _sdVec[i]._shrinkShape2Shape.insert( make_pair( vInd, facesWOL[0] ));
-          break;
         }
+        _sdVec[i]._shrinkShape2Shape.insert( make_pair( vInd, facesWOL[0] ));
+        break;
+      }
       case 2:
+      {
+        // find an edge shared by 2 faces
+        PShapeIteratorPtr eIt = helper.GetAncestors(vertex, *_mesh, TopAbs_EDGE);
+        while ( eIt->more())
         {
-          // find an edge shared by 2 faces
-          PShapeIteratorPtr eIt = helper.GetAncestors(vertex, *_mesh, TopAbs_EDGE);
-          while ( eIt->more())
+          const TopoDS_Shape* e = eIt->next();
+          if ( helper.IsSubShape( *e, facesWOL[0]) &&
+               helper.IsSubShape( *e, facesWOL[1]))
           {
-            const TopoDS_Shape* e = eIt->next();
-            if ( helper.IsSubShape( *e, facesWOL[0]) &&
-                 helper.IsSubShape( *e, facesWOL[1]))
-            {
-              _sdVec[i]._shrinkShape2Shape.insert( make_pair( vInd, *e )); break;
-            }
+            _sdVec[i]._shrinkShape2Shape.insert( make_pair( vInd, *e )); break;
           }
-          break;
         }
+        break;
+      }
       default:
         return error("Not yet supported case", _sdVec[i]._index);
       }
     }
   }
 
+  // add FACEs of other SOLIDs to _ignoreFaceIds
+  for ( size_t i = 0; i < _sdVec.size(); ++i )
+  {
+    shapes.Clear();
+    TopExp::MapShapes(_sdVec[i]._solid, TopAbs_FACE, shapes);
+
+    for ( exp.Init( _mesh->GetShapeToMesh(), TopAbs_FACE ); exp.More(); exp.Next() )
+    {
+      if ( !shapes.Contains( exp.Current() ))
+        _sdVec[i]._ignoreFaceIds.insert( getMeshDS()->ShapeToIndex( exp.Current() ));
+    }
+  }
+
   return true;
 }
 
@@ -1254,10 +1436,10 @@ bool _ViscousBuilder::makeLayer(_SolidData& data)
   subIds = data._noShrinkFaces;
   TopExp_Explorer exp( data._solid, TopAbs_FACE );
   for ( ; exp.More(); exp.Next() )
-    if ( ! _ignoreShapeIds.count( getMeshDS()->ShapeToIndex( exp.Current() )))
     {
       SMESH_subMesh* fSubM = _mesh->GetSubMesh( exp.Current() );
-      faceIds.insert( fSubM->GetId() );
+      if ( ! data._ignoreFaceIds.count( getMeshDS()->ShapeToIndex( exp.Current() )))
+        faceIds.insert( fSubM->GetId() );
       SMESH_subMeshIteratorPtr subIt =
         fSubM->getDependsOnIterator(/*includeSelf=*/true, /*complexShapeFirst=*/false);
       while ( subIt->more() )
@@ -1265,20 +1447,19 @@ bool _ViscousBuilder::makeLayer(_SolidData& data)
     }
 
   // make a map to find new nodes on sub-shapes shared with other SOLID
-  map< TGeomID, TNode2Edge* > s2neMap;
   map< TGeomID, TNode2Edge* >::iterator s2ne;
   map< TGeomID, TopoDS_Shape >::iterator s2s = data._shrinkShape2Shape.begin();
   for (; s2s != data._shrinkShape2Shape.end(); ++s2s )
   {
     TGeomID shapeInd = s2s->first;
-    for ( unsigned i = 0; i < _sdVec.size(); ++i )
+    for ( size_t i = 0; i < _sdVec.size(); ++i )
     {
       if ( _sdVec[i]._index == data._index ) continue;
       map< TGeomID, TopoDS_Shape >::iterator s2s2 = _sdVec[i]._shrinkShape2Shape.find( shapeInd );
       if ( s2s2 != _sdVec[i]._shrinkShape2Shape.end() &&
            *s2s == *s2s2 && !_sdVec[i]._n2eMap.empty() )
       {
-        s2neMap.insert( make_pair( shapeInd, &_sdVec[i]._n2eMap ));
+        data._s2neMap.insert( make_pair( shapeInd, &_sdVec[i]._n2eMap ));
         break;
       }
     }
@@ -1330,21 +1511,23 @@ bool _ViscousBuilder::makeLayer(_SolidData& data)
           const int shapeID = n->getshapeId();
           edgesByGeom[ shapeID ].push_back( edge );
 
+          SMESH_TNodeXYZ xyz( n );
+
           // set edge data or find already refined _LayerEdge and get data from it
           if ( n->GetPosition()->GetTypeOfPosition() != SMDS_TOP_FACE &&
-               ( s2ne = s2neMap.find( shapeID )) != s2neMap.end() &&
+               ( s2ne = data._s2neMap.find( shapeID )) != data._s2neMap.end() &&
                ( n2e2 = (*s2ne).second->find( n )) != s2ne->second->end())
           {
             _LayerEdge* foundEdge = (*n2e2).second;
-            edge->Copy( *foundEdge, helper );
-            // location of the last node is modified but we can restore
-            // it by node position on _sWOL stored by the node
+            gp_XYZ        lastPos = edge->Copy( *foundEdge, helper );
+            foundEdge->_pos.push_back( lastPos );
+            // location of the last node is modified and we restore it by foundEdge->_pos.back()
             const_cast< SMDS_MeshNode* >
-              ( edge->_nodes.back() )->setXYZ( n->X(), n->Y(), n->Z() );
+              ( edge->_nodes.back() )->setXYZ( xyz.X(), xyz.Y(), xyz.Z() );
           }
           else
           {
-            edge->_nodes.push_back( helper.AddNode( n->X(), n->Y(), n->Z() ));
+            edge->_nodes.push_back( helper.AddNode( xyz.X(), xyz.Y(), xyz.Z() ));
             if ( !setEdgeData( *edge, subIds, helper, data ))
               return false;
           }
@@ -1371,14 +1554,16 @@ bool _ViscousBuilder::makeLayer(_SolidData& data)
   if ( data._stepSize < 1. )
     data._epsilon *= data._stepSize;
 
-  // Put _LayerEdge's into a vector
+  // fill data._simplexTestEdges
+  findSimplexTestEdges( data, edgesByGeom );
 
+  // Put _LayerEdge's into the vector data._edges
   if ( !sortEdges( data, edgesByGeom ))
     return false;
 
-  // Set target nodes into _Simplex and _2NearEdges
+  // Set target nodes into _Simplex and _2NearEdges of _LayerEdge's
   TNode2Edge::iterator n2e;
-  for ( unsigned i = 0; i < data._edges.size(); ++i )
+  for ( size_t i = 0; i < data._edges.size(); ++i )
   {
     if ( data._edges[i]->IsOnEdge())
       for ( int j = 0; j < 2; ++j )
@@ -1391,13 +1576,13 @@ bool _ViscousBuilder::makeLayer(_SolidData& data)
         n = (*n2e).second->_nodes.back();
         data._edges[i]->_2neibors->_edges[j] = n2e->second;
       }
-    else
-      for ( unsigned j = 0; j < data._edges[i]->_simplices.size(); ++j )
-      {
-        _Simplex& s = data._edges[i]->_simplices[j];
-        s._nNext = data._n2eMap[ s._nNext ]->_nodes.back();
-        s._nPrev = data._n2eMap[ s._nPrev ]->_nodes.back();
-      }
+    //else
+    for ( size_t j = 0; j < data._edges[i]->_simplices.size(); ++j )
+    {
+      _Simplex& s = data._edges[i]->_simplices[j];
+      s._nNext = data._n2eMap[ s._nNext ]->_nodes.back();
+      s._nPrev = data._n2eMap[ s._nPrev ]->_nodes.back();
+    }
   }
 
   dumpFunctionEnd();
@@ -1445,7 +1630,7 @@ void _ViscousBuilder::limitStepSize( _SolidData&             data,
  */
 //================================================================================
 
-void _ViscousBuilder::limitStepSize( _SolidData& data, const double minSize)
+void _ViscousBuilder::limitStepSize( _SolidData& data, const double minSize )
 {
   if ( minSize < data._stepSize )
   {
@@ -1459,6 +1644,96 @@ void _ViscousBuilder::limitStepSize( _SolidData& data, const double minSize)
   }
 }
 
+//================================================================================
+/*!
+ * Fill data._simplexTestEdges. These _LayerEdge's are used to stop inflation
+ * in the case where there are no _LayerEdge's on a curved convex FACE,
+ * as e.g. on a fillet surface with no internal nodes - issue 22580,
+ * so that collision of viscous internal faces is not detected by check of
+ * intersection of _LayerEdge's with the viscous internal faces.
+ */
+//================================================================================
+
+void _ViscousBuilder::findSimplexTestEdges( _SolidData&                    data,
+                                            vector< vector<_LayerEdge*> >& edgesByGeom)
+{
+  data._simplexTestEdges.clear();
+
+  SMESH_MesherHelper helper( *_mesh );
+
+  vector< vector<_LayerEdge*> * > ledgesOnEdges;
+  set< const SMDS_MeshNode* > usedNodes;
+
+  const double minCurvature = 1. / data._hyp->GetTotalThickness();
+
+  for ( size_t iS = 0; iS < edgesByGeom.size(); ++iS )
+  {
+    // look for a FACE with layers and w/o _LayerEdge's
+    const vector<_LayerEdge*>& eS = edgesByGeom[iS];
+    if ( !eS.empty() ) continue;
+    const TopoDS_Shape& S = getMeshDS()->IndexToShape( iS );
+    if ( S.IsNull() || S.ShapeType() != TopAbs_FACE ) continue;
+    if ( data._ignoreFaceIds.count( iS )) continue;
+
+    const TopoDS_Face& F = TopoDS::Face( S );
+
+    // look for _LayerEdge's on EDGEs with null _sWOL
+    ledgesOnEdges.clear();
+    TopExp_Explorer eExp( F, TopAbs_EDGE );
+    for ( ; eExp.More(); eExp.Next() )
+    {
+      TGeomID iE = getMeshDS()->ShapeToIndex( eExp.Current() );
+      vector<_LayerEdge*>& eE = edgesByGeom[iE];
+      if ( !eE.empty() && eE[0]->_sWOL.IsNull() )
+        ledgesOnEdges.push_back( & eE );
+    }
+    if ( ledgesOnEdges.empty() ) continue;
+
+    // check FACE convexity
+    const _LayerEdge* le = ledgesOnEdges[0]->back();
+    gp_XY uv = helper.GetNodeUV( F, le->_nodes[0] );
+    BRepAdaptor_Surface surf( F );
+    BRepLProp_SLProps surfProp( surf, uv.X(), uv.Y(), 2, 1e-6 );
+    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 = le->_normal;
+    if ( inDir * maxDir.XYZ() < 0 &&
+         inDir * minDir.XYZ() < 0 )
+      continue;
+
+    limitStepSize( data, 0.9 / surfCurvature );
+
+    // add _simplices to the _LayerEdge's
+    for ( size_t iE = 0; iE < ledgesOnEdges.size(); ++iE )
+    {
+      const vector<_LayerEdge*>& ledges = *ledgesOnEdges[iE];
+      for ( size_t iLE = 0; iLE < ledges.size(); ++iLE )
+      {
+        _LayerEdge* ledge = ledges[iLE];
+        const SMDS_MeshNode* srcNode = ledge->_nodes[0];
+        if ( !usedNodes.insert( srcNode ).second ) continue;
+
+        getSimplices( srcNode, ledge->_simplices, data._ignoreFaceIds, &data );
+        for ( size_t i = 0; i < ledge->_simplices.size(); ++i )
+        {
+          usedNodes.insert( ledge->_simplices[i]._nPrev );
+          usedNodes.insert( ledge->_simplices[i]._nNext );
+        }
+        data._simplexTestEdges.push_back( ledge );
+      }
+    }
+  }
+}
+
 //================================================================================
 /*!
  * \brief Separate shapes (and _LayerEdge's on them) to smooth from the rest ones
@@ -1476,11 +1751,11 @@ bool _ViscousBuilder::sortEdges( _SolidData&                    data,
   SMESH_MesherHelper helper( *_mesh );
   bool ok = true;
 
-  for ( unsigned iS = 0; iS < edgesByGeom.size(); ++iS )
+  for ( size_t iS = 0; iS < edgesByGeom.size(); ++iS )
   {
     vector<_LayerEdge*>& eS = edgesByGeom[iS];
     if ( eS.empty() ) continue;
-    TopoDS_Shape S = getMeshDS()->IndexToShape( iS );
+    const TopoDS_Shape& S = getMeshDS()->IndexToShape( iS );
     bool needSmooth = false;
     switch ( S.ShapeType() )
     {
@@ -1520,7 +1795,7 @@ bool _ViscousBuilder::sortEdges( _SolidData&                    data,
         if ( eE.empty() ) continue;
         if ( eE[0]->_sWOL.IsNull() )
         {
-          for ( unsigned i = 0; i < eE.size() && !needSmooth; ++i )
+          for ( size_t i = 0; i < eE.size() && !needSmooth; ++i )
             needSmooth = ( eE[i]->_cosin > 0.1 );
         }
         else
@@ -1528,7 +1803,7 @@ bool _ViscousBuilder::sortEdges( _SolidData&                    data,
           const TopoDS_Face& F1 = TopoDS::Face( S );
           const TopoDS_Face& F2 = TopoDS::Face( eE[0]->_sWOL );
           const TopoDS_Edge& E  = TopoDS::Edge( eExp.Current() );
-          for ( unsigned i = 0; i < eE.size() && !needSmooth; ++i )
+          for ( size_t i = 0; i < eE.size() && !needSmooth; ++i )
           {
             gp_Vec dir1 = getFaceDir( F1, E, eE[i]->_nodes[0], helper, ok );
             gp_Vec dir2 = getFaceDir( F2, E, eE[i]->_nodes[0], helper, ok );
@@ -1567,7 +1842,7 @@ bool _ViscousBuilder::sortEdges( _SolidData&                    data,
   }
 
   // then the rest _LayerEdge's
-  for ( unsigned iS = 0; iS < edgesByGeom.size(); ++iS )
+  for ( size_t iS = 0; iS < edgesByGeom.size(); ++iS )
   {
     vector<_LayerEdge*>& eVec = edgesByGeom[iS];
     data._edges.insert( data._edges.end(), eVec.begin(), eVec.end() );
@@ -1653,32 +1928,62 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge&         edge,
 
     set<TGeomID>::iterator id = faceIds.begin();
     TopoDS_Face F;
+    std::pair< TGeomID, gp_XYZ > id2Norm[20];
     for ( ; id != faceIds.end(); ++id )
     {
       const TopoDS_Shape& s = getMeshDS()->IndexToShape( *id );
       if ( s.IsNull() || s.ShapeType() != TopAbs_FACE || !subIds.count( *id ))
         continue;
-      totalNbFaces++;
-      //nbLayerFaces += subIds.count( *id );
       F = TopoDS::Face( s );
 
       gp_XY uv = helper.GetNodeUV( F, node, 0, &normOK );
       Handle(Geom_Surface) surface = BRep_Tool::Surface( F );
-      surface->D1( uv.X(),uv.Y(), p, du,dv );
-      geomNorm = du ^ dv;
-      double size2 = geomNorm.SquareMagnitude();
-      if ( size2 > numeric_limits<double>::min() )
-        geomNorm /= sqrt( size2 );
-      else
-        normOK = false;
+      {
+        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;
+      id2Norm[ totalNbFaces ].second = geomNorm.XYZ();
+      totalNbFaces++;
       edge._normal += geomNorm.XYZ();
     }
     if ( totalNbFaces == 0 )
       return error(SMESH_Comment("Can't get normal to node ") << node->GetID(), data._index);
 
-    edge._normal /= totalNbFaces;
+    if ( totalNbFaces < 3 )
+    {
+      //edge._normal /= totalNbFaces;
+    }
+    else
+    {
+      edge._normal = getWeigthedNormal( node, id2Norm, totalNbFaces );
+    }
 
     switch ( posType )
     {
@@ -1695,8 +2000,8 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge&         edge,
     }
     case SMDS_TOP_VERTEX: {
       TopoDS_Vertex V = TopoDS::Vertex( helper.GetSubShapeByNode( node, getMeshDS()));
-      gp_Vec inFaceDir = getFaceDir( F, V, node, helper, normOK);
-      double angle = inFaceDir.Angle( edge._normal ); // [0,PI]
+      gp_XYZ inFaceDir = getFaceDir( F, V, node, helper, normOK);
+      double angle = gp_Vec( inFaceDir).Angle( edge._normal ); // [0,PI]
       edge._cosin = cos( angle );
       //cout << "Cosin on VERTEX " << edge._cosin << " node " << node->GetID() << endl;
       break;
@@ -1744,9 +2049,9 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge&         edge,
 
     if ( posType == SMDS_TOP_FACE )
     {
-      getSimplices( node, edge._simplices, _ignoreShapeIds, &data );
+      getSimplices( node, edge._simplices, data._ignoreFaceIds, &data );
       double avgNormProj = 0, avgLen = 0;
-      for ( unsigned i = 0; i < edge._simplices.size(); ++i )
+      for ( size_t i = 0; i < edge._simplices.size(); ++i )
       {
         gp_XYZ vec = edge._pos.back() - SMESH_TNodeXYZ( edge._simplices[i]._nPrev );
         avgNormProj += edge._normal * vec;
@@ -1780,6 +2085,87 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge&         edge,
   return true;
 }
 
+//================================================================================
+/*!
+ * \brief Return 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
+ *  \return gp_XYZ - computed normal
+ */
+//================================================================================
+
+gp_XYZ _ViscousBuilder::getWeigthedNormal( const SMDS_MeshNode*         n,
+                                           std::pair< TGeomID, gp_XYZ > fId2Normal[],
+                                           const int                    nbFaces )
+{
+  gp_XYZ resNorm(0,0,0);
+  TopoDS_Shape V = SMESH_MesherHelper::GetSubShapeByNode( n, getMeshDS() );
+  if ( V.ShapeType() != TopAbs_VERTEX )
+  {
+    for ( int i = 0; i < nbFaces; ++i )
+      resNorm += fId2Normal[i].second / nbFaces ;
+    return resNorm;
+  }
+
+  double angles[30];
+  for ( int i = 0; i < nbFaces; ++i )
+  {
+    const TopoDS_Face& F = TopoDS::Face( getMeshDS()->IndexToShape( fId2Normal[i].first ));
+
+    // look for two EDGEs shared by F and other FACEs within fId2Normal
+    TopoDS_Edge ee[2];
+    int nbE = 0;
+    PShapeIteratorPtr eIt = SMESH_MesherHelper::GetAncestors( V, *_mesh, TopAbs_EDGE );
+    while ( const TopoDS_Shape* E = eIt->next() )
+    {
+      if ( !SMESH_MesherHelper::IsSubShape( *E, F ))
+        continue;
+      bool isSharedEdge = false;
+      for ( int j = 0; j < nbFaces && !isSharedEdge; ++j )
+      {
+        if ( i == j ) continue;
+        const TopoDS_Shape& otherF = getMeshDS()->IndexToShape( fId2Normal[j].first );
+        isSharedEdge = SMESH_MesherHelper::IsSubShape( *E, otherF );
+      }
+      if ( !isSharedEdge )
+        continue;
+      ee[ nbE ] = TopoDS::Edge( *E );
+      ee[ nbE ].Orientation( SMESH_MesherHelper::GetSubShapeOri( F, *E ));
+      if ( ++nbE == 2 )
+        break;
+    }
+
+    // get an angle between the two EDGEs
+    angles[i] = 0;
+    if ( nbE < 1 ) continue;
+    if ( nbE == 1 )
+    {
+      ee[ 1 ] == ee[ 0 ];
+    }
+    else
+    {
+      TopoDS_Vertex v10 = SMESH_MesherHelper::IthVertex( 1, ee[ 0 ]);
+      TopoDS_Vertex v01 = SMESH_MesherHelper::IthVertex( 0, ee[ 1 ]);
+      if ( !v10.IsSame( v01 ))
+        std::swap( ee[0], ee[1] );
+    }
+    angles[i] = SMESH_MesherHelper::GetAngle( ee[0], ee[1], F );
+  }
+
+  // compute a weighted normal
+  double sumAngle = 0;
+  for ( int i = 0; i < nbFaces; ++i )
+  {
+    angles[i] = ( angles[i] > 2*M_PI )  ?  0  :  M_PI - angles[i];
+    sumAngle += angles[i];
+  }
+  for ( int i = 0; i < nbFaces; ++i )
+    resNorm += angles[i] / sumAngle * fId2Normal[i].second;
+
+  return resNorm;
+}
+
 //================================================================================
 /*!
  * \brief Find 2 neigbor nodes of a node on EDGE
@@ -1883,7 +2269,7 @@ void _LayerEdge::SetDataByNeighbors( const SMDS_MeshNode* n1,
  */
 //================================================================================
 
-void _LayerEdge::Copy( _LayerEdge& other, SMESH_MesherHelper& helper )
+gp_XYZ _LayerEdge::Copy( _LayerEdge& other, SMESH_MesherHelper& helper )
 {
   _nodes     = other._nodes;
   _normal    = other._normal;
@@ -1895,16 +2281,25 @@ void _LayerEdge::Copy( _LayerEdge& other, SMESH_MesherHelper& helper )
   _curvature = 0; std::swap( _curvature, other._curvature );
   _2neibors  = 0; std::swap( _2neibors,  other._2neibors );
 
+  gp_XYZ lastPos( 0,0,0 );
   if ( _sWOL.ShapeType() == TopAbs_EDGE )
   {
     double u = helper.GetNodeU( TopoDS::Edge( _sWOL ), _nodes[0] );
     _pos.push_back( gp_XYZ( u, 0, 0));
+
+    u = helper.GetNodeU( TopoDS::Edge( _sWOL ), _nodes.back() );
+    lastPos.SetX( u );
   }
   else // TopAbs_FACE
   {
     gp_XY uv = helper.GetNodeUV( TopoDS::Face( _sWOL ), _nodes[0]);
     _pos.push_back( gp_XYZ( uv.X(), uv.Y(), 0));
+
+    uv = helper.GetNodeUV( TopoDS::Face( _sWOL ), _nodes.back() );
+    lastPos.SetX( uv.X() );
+    lastPos.SetY( uv.Y() );
   }
+  return lastPos;
 }
 
 //================================================================================
@@ -1916,7 +2311,7 @@ void _LayerEdge::Copy( _LayerEdge& other, SMESH_MesherHelper& helper )
 void _LayerEdge::SetCosin( double cosin )
 {
   _cosin = cosin;
-  _lenFactor = ( _cosin > 0.1 ) ?  1./sqrt(1-_cosin*_cosin) : 1.0;
+  _lenFactor = ( Abs( _cosin ) > 0.1 ) ?  1./sqrt(1-_cosin*_cosin) : 1.0;
 }
 
 //================================================================================
@@ -1931,19 +2326,21 @@ void _ViscousBuilder::getSimplices( const SMDS_MeshNode* node,
                                     const _SolidData*    dataToCheckOri,
                                     const bool           toSort)
 {
+  simplices.clear();
   SMDS_ElemIteratorPtr fIt = node->GetInverseElementIterator(SMDSAbs_Face);
   while ( fIt->more() )
   {
     const SMDS_MeshElement* f = fIt->next();
-    const TGeomID shapeInd = f->getshapeId();
+    const TGeomID    shapeInd = f->getshapeId();
     if ( ingnoreShapes.count( shapeInd )) continue;
     const int nbNodes = f->NbCornerNodes();
-    int srcInd = f->GetNodeIndex( node );
+    const int  srcInd = f->GetNodeIndex( node );
     const SMDS_MeshNode* nPrev = f->GetNode( SMESH_MesherHelper::WrapIndex( srcInd-1, nbNodes ));
     const SMDS_MeshNode* nNext = f->GetNode( SMESH_MesherHelper::WrapIndex( srcInd+1, nbNodes ));
+    const SMDS_MeshNode* nOpp  = f->GetNode( SMESH_MesherHelper::WrapIndex( srcInd+2, nbNodes ));
     if ( dataToCheckOri && dataToCheckOri->_reversedFaceIds.count( shapeInd ))
       std::swap( nPrev, nNext );
-    simplices.push_back( _Simplex( nPrev, nNext ));
+    simplices.push_back( _Simplex( nPrev, nNext, nOpp ));
   }
 
   if ( toSort )
@@ -1975,7 +2372,7 @@ void _ViscousBuilder::getSimplices( const SMDS_MeshNode* node,
 void _ViscousBuilder::makeGroupOfLE()
 {
 #ifdef _DEBUG_
-  for ( unsigned i = 0 ; i < _sdVec.size(); ++i )
+  for ( size_t i = 0 ; i < _sdVec.size(); ++i )
   {
     if ( _sdVec[i]._edges.empty() ) continue;
 //     string name = SMESH_Comment("_LayerEdge's_") << i;
@@ -1985,10 +2382,10 @@ void _ViscousBuilder::makeGroupOfLE()
 //     SMESHDS_Mesh* mDS = _mesh->GetMeshDS();
 
     dumpFunction( SMESH_Comment("make_LayerEdge_") << i );
-    for ( unsigned j = 0 ; j < _sdVec[i]._edges.size(); ++j )
+    for ( size_t j = 0 ; j < _sdVec[i]._edges.size(); ++j )
     {
       _LayerEdge* le = _sdVec[i]._edges[j];
-      for ( unsigned iN = 1; iN < le->_nodes.size(); ++iN )
+      for ( size_t iN = 1; iN < le->_nodes.size(); ++iN )
         dumpCmd(SMESH_Comment("mesh.AddEdge([ ") <<le->_nodes[iN-1]->GetID()
                 << ", " << le->_nodes[iN]->GetID() <<"])");
       //gDS->SMDSGroup().Add( mDS->AddEdge( le->_nodes[iN-1], le->_nodes[iN]));
@@ -1996,7 +2393,7 @@ void _ViscousBuilder::makeGroupOfLE()
     dumpFunctionEnd();
 
     dumpFunction( SMESH_Comment("makeNormals") << i );
-    for ( unsigned j = 0 ; j < _sdVec[i]._edges.size(); ++j )
+    for ( size_t j = 0 ; j < _sdVec[i]._edges.size(); ++j )
     {
       _LayerEdge& edge = *_sdVec[i]._edges[j];
       SMESH_TNodeXYZ nXYZ( edge._nodes[0] );
@@ -2047,14 +2444,14 @@ bool _ViscousBuilder::inflate(_SolidData& data)
   // Limit inflation step size by geometry size found by itersecting
   // normals of _LayerEdge's with mesh faces
   double geomSize = Precision::Infinite(), intersecDist;
-  SMESH_MeshEditor editor( _mesh );
   auto_ptr<SMESH_ElementSearcher> searcher
-    ( editor.GetElementSearcher( data._proxyMesh->GetFaces( data._solid )) );
-  for ( unsigned i = 0; i < data._edges.size(); ++i )
+    ( SMESH_MeshAlgos::GetElementSearcher( *getMeshDS(),
+                                           data._proxyMesh->GetFaces( data._solid )) );
+  for ( size_t i = 0; i < data._edges.size(); ++i )
   {
     if ( data._edges[i]->IsOnEdge() ) continue;
     data._edges[i]->FindIntersection( *searcher, intersecDist, data._epsilon );
-    if ( geomSize > intersecDist )
+    if ( geomSize > intersecDist && intersecDist > 0 )
       geomSize = intersecDist;
   }
   if ( data._stepSize > 0.3 * geomSize )
@@ -2085,7 +2482,7 @@ bool _ViscousBuilder::inflate(_SolidData& data)
 
     // Elongate _LayerEdge's
     dumpFunction(SMESH_Comment("inflate")<<data._index<<"_step"<<nbSteps); // debug
-    for ( unsigned i = 0; i < data._edges.size(); ++i )
+    for ( size_t i = 0; i < data._edges.size(); ++i )
     {
       data._edges[i]->SetNewLength( curThick, helper );
     }
@@ -2101,7 +2498,7 @@ bool _ViscousBuilder::inflate(_SolidData& data)
       if ( nbSteps > 0 )
       {
         dumpFunction(SMESH_Comment("invalidate")<<data._index<<"_step"<<nbSteps); // debug
-        for ( unsigned i = 0; i < data._edges.size(); ++i )
+        for ( size_t i = 0; i < data._edges.size(); ++i )
         {
           data._edges[i]->InvalidateStep( nbSteps+1 );
         }
@@ -2113,7 +2510,7 @@ bool _ViscousBuilder::inflate(_SolidData& data)
 
     // Evaluate achieved thickness
     avgThick = 0;
-    for ( unsigned i = 0; i < data._edges.size(); ++i )
+    for ( size_t i = 0; i < data._edges.size(); ++i )
       avgThick += data._edges[i]->_len;
     avgThick /= data._edges.size();
 #ifdef __myDEBUG
@@ -2161,7 +2558,7 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data,
   TopoDS_Face F;
 
   int iBeg, iEnd = 0;
-  for ( unsigned iS = 0; iS < data._endEdgeToSmooth.size(); ++iS )
+  for ( size_t iS = 0; iS < data._endEdgeToSmooth.size(); ++iS )
   {
     iBeg = iEnd;
     iEnd = data._endEdgeToSmooth[ iS ];
@@ -2206,8 +2603,8 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data,
     else
     {
       // smooth on FACE's
-      int step = 0, badNb = 0; moved = true;
-      while (( ++step <= 5 && moved ) || improved )
+      int step = 0, stepLimit = 5, badNb = 0; moved = true;
+      while (( ++step <= stepLimit && moved ) || improved )
       {
         dumpFunction(SMESH_Comment("smooth")<<data._index<<"_Fa"<<sInd
                      <<"_InfStep"<<nbSteps<<"_"<<step); // debug
@@ -2218,6 +2615,10 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data,
           moved |= data._edges[i]->Smooth(badNb);
         improved = ( badNb < oldBadNb );
 
+        // issue 22576. no bad faces but still there are intersections to fix
+        if ( improved && badNb == 0 )
+          stepLimit = step + 3;
+
         dumpFunctionEnd();
       }
       if ( badNb > 0 )
@@ -2227,7 +2628,7 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data,
         {
           _LayerEdge* edge = data._edges[i];
           SMESH_TNodeXYZ tgtXYZ( edge->_nodes.back() );
-          for ( unsigned j = 0; j < edge->_simplices.size(); ++j )
+          for ( size_t j = 0; j < edge->_simplices.size(); ++j )
             if ( !edge->_simplices[j].IsForward( edge->_nodes[0], &tgtXYZ ))
             {
               cout << "Bad simplex ( " << edge->_nodes[0]->GetID()<< " "<< tgtXYZ._node->GetID()
@@ -2242,12 +2643,30 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data,
     }
   } // loop on shapes to smooth
 
+  // Check orientation of simplices of _simplexTestEdges
+  for ( size_t i = 0; i < data._simplexTestEdges.size(); ++i )
+  {
+    const _LayerEdge* edge = data._simplexTestEdges[i];
+    SMESH_TNodeXYZ tgtXYZ( edge->_nodes.back() );
+    for ( size_t j = 0; j < edge->_simplices.size(); ++j )
+      if ( !edge->_simplices[j].IsForward( edge->_nodes[0], &tgtXYZ ))
+      {
+#ifdef __myDEBUG
+        cout << "Bad simplex of _simplexTestEdges ("
+             << " "<< edge->_nodes[0]->GetID()<< " "<< tgtXYZ._node->GetID()
+             << " "<< edge->_simplices[j]._nPrev->GetID()
+             << " "<< edge->_simplices[j]._nNext->GetID() << " )" << endl;
+#endif
+        return false;
+      }
+  }
+
   // Check if the last segments of _LayerEdge intersects 2D elements;
   // checked elements are either temporary faces or faces on surfaces w/o the layers
 
-  SMESH_MeshEditor editor( _mesh );
   auto_ptr<SMESH_ElementSearcher> searcher
-    ( editor.GetElementSearcher( data._proxyMesh->GetFaces( data._solid )) );
+    ( SMESH_MeshAlgos::GetElementSearcher( *getMeshDS(),
+                                           data._proxyMesh->GetFaces( data._solid )) );
 
   distToIntersection = Precision::Infinite();
   double dist;
@@ -2256,7 +2675,7 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data,
   const SMDS_MeshElement* closestFace = 0;
   int iLE = 0;
 #endif
-  for ( unsigned i = 0; i < data._edges.size(); ++i )
+  for ( size_t i = 0; i < data._edges.size(); ++i )
   {
     if ( data._edges[i]->FindIntersection( *searcher, dist, data._epsilon, &intFace ))
       return false;
@@ -2502,6 +2921,10 @@ bool _ViscousBuilder::smoothAnalyticEdge( _SolidData&           data,
 
     if ( F.IsNull() ) // 3D
     {
+      if ( data._edges[iFrom]->_2neibors->_nodes[0] ==
+           data._edges[iTo-1]->_2neibors->_nodes[1] )
+        return true; // closed EDGE - nothing to do
+
       return false; // TODO ???
     }
     else // 2D
@@ -2563,7 +2986,7 @@ bool _ViscousBuilder::updateNormals( _SolidData&         data,
     vector< const SMDS_MeshNode*> nodes(4); // of a tmp mesh face
 
     dumpFunction(SMESH_Comment("makeTmpFacesOnEdges")<<data._index);
-    for ( unsigned i = 0; i < data._edges.size(); ++i )
+    for ( size_t i = 0; i < data._edges.size(); ++i )
     {
       _LayerEdge* edge = data._edges[i];
       if ( !edge->IsOnEdge() || !edge->_sWOL.IsNull() ) continue;
@@ -2580,7 +3003,7 @@ bool _ViscousBuilder::updateNormals( _SolidData&         data,
         }
         // look for a _LayerEdge containg tgt2
 //         _LayerEdge* neiborEdge = 0;
-//         unsigned di = 0; // check _edges[i+di] and _edges[i-di]
+//         size_t di = 0; // check _edges[i+di] and _edges[i-di]
 //         while ( !neiborEdge && ++di <= data._edges.size() )
 //         {
 //           if ( i+di < data._edges.size() && data._edges[i+di]->_nodes.back() == tgt2 )
@@ -2607,10 +3030,10 @@ bool _ViscousBuilder::updateNormals( _SolidData&         data,
   // 1) to find and fix intersection
   // 2) to check that no new intersection appears as result of 1)
 
-  SMESH_MeshEditor editor( _mesh );
   SMDS_ElemIteratorPtr fIt( new SMDS_ElementVectorIterator( tmpFaces.begin(),
                                                             tmpFaces.end()));
-  auto_ptr<SMESH_ElementSearcher> searcher ( editor.GetElementSearcher( fIt ));
+  auto_ptr<SMESH_ElementSearcher> searcher
+    ( SMESH_MeshAlgos::GetElementSearcher( *getMeshDS(), fIt ));
 
   // 1) Find intersections
   double dist;
@@ -2619,7 +3042,7 @@ bool _ViscousBuilder::updateNormals( _SolidData&         data,
   TLEdge2LEdgeSet edge2CloseEdge;
 
   const double eps = data._epsilon * data._epsilon;
-  for ( unsigned i = 0; i < data._edges.size(); ++i )
+  for ( size_t i = 0; i < data._edges.size(); ++i )
   {
     _LayerEdge* edge = data._edges[i];
     if ( !edge->IsOnEdge() || !edge->_sWOL.IsNull() ) continue;
@@ -2690,19 +3113,19 @@ bool _ViscousBuilder::updateNormals( _SolidData&         data,
       if ( FF1[0].IsNull() || FF2[0].IsNull() )
         continue;
 
-//       // get a new normal for edge1
+      // 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();
+      // 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 );
@@ -2716,7 +3139,7 @@ bool _ViscousBuilder::updateNormals( _SolidData&         data,
       n1 = edge1->_2neibors->_edges[0]->_nodes[0];
       n2 = edge1->_2neibors->_edges[1]->_nodes[0];
       //if ( !findNeiborsOnEdge( edge1, n1, n2, data ))
-      //continue;
+      //  continue;
       edge1->SetDataByNeighbors( n1, n2, helper );
       gp_Vec dirInFace;
       if ( edge1->_cosin < 0 )
@@ -2790,7 +3213,7 @@ bool _ViscousBuilder::updateNormals( _SolidData&         data,
   // 2) Check absence of intersections
   // TODO?
 
-  for ( unsigned i = 0 ; i < tmpFaces.size(); ++i )
+  for ( size_t i = 0 ; i < tmpFaces.size(); ++i )
     delete tmpFaces[i];
 
   return true;
@@ -2816,7 +3239,7 @@ bool _LayerEdge::FindIntersection( SMESH_ElementSearcher&   searcher,
   bool segmentIntersected = false;
   distance = Precision::Infinite();
   int iFace = -1; // intersected face
-  for ( unsigned 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 ||
@@ -2844,7 +3267,7 @@ bool _LayerEdge::FindIntersection( SMESH_ElementSearcher&   searcher,
     }
     if ( intFound )
     {
-      if ( dist < segLen*(1.01))
+      if ( dist < segLen*(1.01) && dist > -(_len-segLen) )
         segmentIntersected = true;
       if ( distance > dist )
         distance = dist, iFace = j;
@@ -3047,7 +3470,8 @@ bool _LayerEdge::SmoothOnEdge(Handle(Geom_Surface)& surface,
   double lenDelta = 0;
   if ( _curvature )
   {
-    lenDelta = _curvature->lenDelta( _len );
+    //lenDelta = _curvature->lenDelta( _len );
+    lenDelta = _curvature->lenDeltaByDist( dist01 );
     newPos.ChangeCoord() += _normal * lenDelta;
   }
 
@@ -3103,7 +3527,7 @@ bool _LayerEdge::Smooth(int& badNb)
 
   // compute new position for the last _pos
   gp_XYZ newPos (0,0,0);
-  for ( unsigned i = 0; i < _simplices.size(); ++i )
+  for ( size_t i = 0; i < _simplices.size(); ++i )
     newPos += SMESH_TNodeXYZ( _simplices[i]._nPrev );
   newPos /= _simplices.size();
 
@@ -3126,11 +3550,11 @@ bool _LayerEdge::Smooth(int& badNb)
   // count quality metrics (orientation) of tetras around _tgtNode
   int nbOkBefore = 0;
   SMESH_TNodeXYZ tgtXYZ( _nodes.back() );
-  for ( unsigned i = 0; i < _simplices.size(); ++i )
+  for ( size_t i = 0; i < _simplices.size(); ++i )
     nbOkBefore += _simplices[i].IsForward( _nodes[0], &tgtXYZ );
 
   int nbOkAfter = 0;
-  for ( unsigned i = 0; i < _simplices.size(); ++i )
+  for ( size_t i = 0; i < _simplices.size(); ++i )
     nbOkAfter += _simplices[i].IsForward( _nodes[0], &newPos );
 
   if ( nbOkAfter < nbOkBefore )
@@ -3251,19 +3675,23 @@ bool _ViscousBuilder::refine(_SolidData& data)
   Handle(Geom_Surface) surface;
   TopoDS_Edge geomEdge;
   TopoDS_Face geomFace;
+  TopoDS_Shape prevSWOL;
   TopLoc_Location loc;
-  double f,l, u/*, distXYZ[4]*/;
+  double f,l, u;
   gp_XY uv;
   bool isOnEdge;
+  TGeomID prevBaseId = -1;
+  TNode2Edge* n2eMap = 0;
+  TNode2Edge::iterator n2e;
 
-  for ( unsigned i = 0; i < data._edges.size(); ++i )
+  for ( size_t i = 0; i < data._edges.size(); ++i )
   {
     _LayerEdge& edge = *data._edges[i];
 
     // get accumulated length of segments
     vector< double > segLen( edge._pos.size() );
     segLen[0] = 0.0;
-    for ( unsigned j = 1; j < edge._pos.size(); ++j )
+    for ( size_t j = 1; j < edge._pos.size(); ++j )
       segLen[j] = segLen[j-1] + (edge._pos[j-1] - edge._pos[j] ).Modulus();
 
     // allocate memory for new nodes if it is not yet refined
@@ -3274,27 +3702,46 @@ bool _ViscousBuilder::refine(_SolidData& data)
       edge._nodes[1] = 0;
       edge._nodes.back() = tgtNode;
     }
-    if ( !edge._sWOL.IsNull() )
+    // get data of a shrink shape
+    if ( !edge._sWOL.IsNull() && edge._sWOL != prevSWOL )
     {
       isOnEdge = ( edge._sWOL.ShapeType() == TopAbs_EDGE );
-      // restore position of the last node
-//       gp_Pnt p;
       if ( isOnEdge )
       {
         geomEdge = TopoDS::Edge( edge._sWOL );
-        curve = BRep_Tool::Curve( geomEdge, loc, f,l);
-//         double u = helper.GetNodeU( tgtNode );
-//         p = curve->Value( u );
+        curve    = BRep_Tool::Curve( geomEdge, loc, f,l);
       }
       else
       {
         geomFace = TopoDS::Face( edge._sWOL );
-        surface = BRep_Tool::Surface( geomFace, loc );
-//         gp_XY uv = helper.GetNodeUV( tgtNode );
-//         p = surface->Value( uv.X(), uv.Y() );
+        surface  = BRep_Tool::Surface( geomFace, loc );
+      }
+      prevSWOL = edge._sWOL;
+    }
+    // restore shapePos of the last node by already treated _LayerEdge of another _SolidData
+    const TGeomID baseShapeId = edge._nodes[0]->getshapeId();
+    if ( baseShapeId != prevBaseId )
+    {
+      map< TGeomID, TNode2Edge* >::iterator s2ne = data._s2neMap.find( baseShapeId );
+      n2eMap = ( s2ne == data._s2neMap.end() ) ? 0 : n2eMap = s2ne->second;
+      prevBaseId = baseShapeId;
+    }
+    if ( n2eMap && (( n2e = n2eMap->find( edge._nodes[0] )) != n2eMap->end() ))
+    {
+      _LayerEdge*  foundEdge = n2e->second;
+      const gp_XYZ& foundPos = foundEdge->_pos.back();
+      SMDS_PositionPtr lastPos = tgtNode->GetPosition();
+      if ( isOnEdge )
+      {
+        SMDS_EdgePosition* epos = static_cast<SMDS_EdgePosition*>( lastPos );
+        epos->SetUParameter( foundPos.X() );
+      }
+      else
+      {
+        SMDS_FacePosition* fpos = static_cast<SMDS_FacePosition*>( lastPos );
+        fpos->SetUParameter( foundPos.X() );
+        fpos->SetVParameter( foundPos.Y() );
       }
-//       p.Transform( loc );
-//       const_cast< SMDS_MeshNode* >( tgtNode )->setXYZ( p.X(), p.Y(), p.Z() );
     }
     // calculate height of the first layer
     double h0;
@@ -3311,8 +3758,8 @@ bool _ViscousBuilder::refine(_SolidData& data)
 
     // create intermediate nodes
     double hSum = 0, hi = h0/f;
-    unsigned iSeg = 1;
-    for ( unsigned iStep = 1; iStep < edge._nodes.size(); ++iStep )
+    size_t iSeg = 1;
+    for ( size_t iStep = 1; iStep < edge._nodes.size(); ++iStep )
     {
       // compute an intermediate position
       hi *= f;
@@ -3332,12 +3779,14 @@ bool _ViscousBuilder::refine(_SolidData& data)
         if ( isOnEdge )
         {
           u = pos.X();
-          pos = curve->Value( u ).Transformed(loc);
+          if ( !node )
+            pos = curve->Value( u ).Transformed(loc);
         }
         else
         {
           uv.SetCoord( pos.X(), pos.Y() );
-          pos = surface->Value( pos.X(), pos.Y() ).Transformed(loc);
+          if ( !node )
+            pos = surface->Value( pos.X(), pos.Y() ).Transformed(loc);
         }
       }
       // create or update the node
@@ -3365,11 +3814,18 @@ bool _ViscousBuilder::refine(_SolidData& data)
           {
             u = 0.5 * ( u + helper.GetNodeU( geomEdge, node ));
             pos = curve->Value( u ).Transformed(loc);
+
+            SMDS_EdgePosition* epos = static_cast<SMDS_EdgePosition*>( node->GetPosition() );
+            epos->SetUParameter( u );
           }
           else
           {
             uv = 0.5 * ( uv + helper.GetNodeUV( geomFace, node ));
             pos = surface->Value( uv.X(), uv.Y()).Transformed(loc);
+
+            SMDS_FacePosition* fpos = static_cast<SMDS_FacePosition*>( node->GetPosition() );
+            fpos->SetUParameter( uv.X() );
+            fpos->SetVParameter( uv.Y() );
           }
         }
         node->setXYZ( pos.X(), pos.Y(), pos.Z() );
@@ -3379,7 +3835,7 @@ bool _ViscousBuilder::refine(_SolidData& data)
 
   if ( !getMeshDS()->IsEmbeddedMode() )
     // Log node movement
-    for ( unsigned i = 0; i < data._edges.size(); ++i )
+    for ( size_t i = 0; i < data._edges.size(); ++i )
     {
       _LayerEdge& edge = *data._edges[i];
       SMESH_TNodeXYZ p ( edge._nodes.back() );
@@ -3393,7 +3849,7 @@ bool _ViscousBuilder::refine(_SolidData& data)
   TopExp_Explorer exp( data._solid, TopAbs_FACE );
   for ( ; exp.More(); exp.Next() )
   {
-    if ( _ignoreShapeIds.count( getMeshDS()->ShapeToIndex( exp.Current() )))
+    if ( data._ignoreFaceIds.count( getMeshDS()->ShapeToIndex( exp.Current() )))
       continue;
     SMESHDS_SubMesh* fSubM = getMeshDS()->MeshElements( exp.Current() );
     SMDS_ElemIteratorPtr fIt = fSubM->GetElements();
@@ -3444,7 +3900,7 @@ bool _ViscousBuilder::shrink()
   // make map of (ids of FACEs to shrink mesh on) to (_SolidData containing _LayerEdge's
   // inflated along FACE or EDGE)
   map< TGeomID, _SolidData* > f2sdMap;
-  for ( unsigned i = 0 ; i < _sdVec.size(); ++i )
+  for ( size_t i = 0 ; i < _sdVec.size(); ++i )
   {
     _SolidData& data = _sdVec[i];
     TopTools_MapOfShape FFMap;
@@ -3529,7 +3985,7 @@ bool _ViscousBuilder::shrink()
         sm->getDependsOnIterator(/*includeSelf=*/false, /*complexShapeFirst=*/false);
       while ( subIt->more() )
       {
-        SMESH_subMesh* sub = subIt->next();
+        SMESH_subMesh*     sub = subIt->next();
         SMESHDS_SubMesh* subDS = sub->GetSubMeshDS();
         if ( subDS->NbNodes() == 0 || !n2eMap.count( subDS->GetNodes()->next() ))
           continue;
@@ -3543,9 +3999,15 @@ bool _ViscousBuilder::shrink()
       }
     }
 
+    dumpFunction(SMESH_Comment("beforeShrinkFace")<<f2sd->first); // debug
+    SMDS_ElemIteratorPtr fIt = smDS->GetElements();
+    while ( fIt->more() )
+      if ( const SMDS_MeshElement* f = fIt->next() )
+        dumpChangeNodes( f );
+
     // Replace source nodes by target nodes in mesh faces to shrink
     const SMDS_MeshNode* nodes[20];
-    for ( unsigned i = 0; i < lEdges.size(); ++i )
+    for ( size_t i = 0; i < lEdges.size(); ++i )
     {
       _LayerEdge& edge = *lEdges[i];
       const SMDS_MeshNode* srcNode = edge._nodes[0];
@@ -3556,10 +4018,10 @@ bool _ViscousBuilder::shrink()
         const SMDS_MeshElement* f = fIt->next();
         if ( !smDS->Contains( f ))
           continue;
-        SMDS_ElemIteratorPtr nIt = f->nodesIterator();
-        for ( int iN = 0; iN < f->NbNodes(); ++iN )
+        SMDS_NodeIteratorPtr nIt = f->nodeIterator();
+        for ( int iN = 0; nIt->more(); ++iN )
         {
-          const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nIt->next() );
+          const SMDS_MeshNode* n = nIt->next();
           nodes[iN] = ( n == srcNode ? tgtNode : n );
         }
         helper.GetMeshDS()->ChangeElementNodes( f, nodes, f->NbNodes() );
@@ -3572,25 +4034,24 @@ bool _ViscousBuilder::shrink()
     // Create _SmoothNode's on face F
     vector< _SmoothNode > nodesToSmooth( smoothNodes.size() );
     {
-      dumpFunction(SMESH_Comment("beforeShrinkFace")<<f2sd->first); // debug
-      for ( unsigned i = 0; i < smoothNodes.size(); ++i )
+      const bool sortSimplices = isConcaveFace;
+      for ( size_t i = 0; i < smoothNodes.size(); ++i )
       {
         const SMDS_MeshNode* n = smoothNodes[i];
         nodesToSmooth[ i ]._node = n;
         // src nodes must be replaced by tgt nodes to have tgt nodes in _simplices
-        getSimplices( n, nodesToSmooth[ i ]._simplices, ignoreShapes, NULL, isConcaveFace );
+        getSimplices( n, nodesToSmooth[ i ]._simplices, ignoreShapes, NULL, sortSimplices );
         // fix up incorrect uv of nodes on the FACE
         helper.GetNodeUV( F, n, 0, &isOkUV);
         dumpMove( n );
       }
-      dumpFunctionEnd();
     }
     //if ( nodesToSmooth.empty() ) continue;
 
-    // Find EDGE's to shrink
+    // Find EDGE's to shrink and set simpices to LayerEdge's
     set< _Shrinker1D* > eShri1D;
     {
-      for ( unsigned i = 0; i < lEdges.size(); ++i )
+      for ( size_t i = 0; i < lEdges.size(); ++i )
       {
         _LayerEdge* edge = lEdges[i];
         if ( edge->_sWOL.ShapeType() == TopAbs_EDGE )
@@ -3599,10 +4060,28 @@ bool _ViscousBuilder::shrink()
           _Shrinker1D& srinker = e2shrMap[ edgeIndex ];
           eShri1D.insert( & srinker );
           srinker.AddEdge( edge, helper );
+          VISCOUS_3D::ToClearSubWithMain( _mesh->GetSubMesh( edge->_sWOL ), data._solid );
           // restore params of nodes on EGDE if the EDGE has been already
           // srinked while srinking another FACE
           srinker.RestoreParams();
         }
+        getSimplices( /*tgtNode=*/edge->_nodes.back(), edge->_simplices, ignoreShapes );
+      }
+    }
+
+    bool toFixTria = false; // to improve quality of trias by diagonal swap
+    if ( isConcaveFace )
+    {
+      const bool hasTria = _mesh->NbTriangles(), hasQuad = _mesh->NbQuadrangles();
+      if ( hasTria != hasQuad ) {
+        toFixTria = hasTria;
+      }
+      else {
+        set<int> nbNodesSet;
+        SMDS_ElemIteratorPtr fIt = smDS->GetElements();
+        while ( fIt->more() && nbNodesSet.size() < 2 )
+          nbNodesSet.insert( fIt->next()->NbCornerNodes() );
+        toFixTria = ( *nbNodesSet.begin() == 3 );
       }
     }
 
@@ -3612,19 +4091,23 @@ bool _ViscousBuilder::shrink()
 
     bool shrinked = true;
     int badNb, shriStep=0, smooStep=0;
+    _SmoothNode::SmoothType smoothType
+      = isConcaveFace ? _SmoothNode::ANGULAR : _SmoothNode::LAPLACIAN;
     while ( shrinked )
     {
+      shriStep++;
       // Move boundary nodes (actually just set new UV)
       // -----------------------------------------------
-      dumpFunction(SMESH_Comment("moveBoundaryOnF")<<f2sd->first<<"_st"<<shriStep++ ); // debug
+      dumpFunction(SMESH_Comment("moveBoundaryOnF")<<f2sd->first<<"_st"<<shriStep ); // debug
       shrinked = false;
-      for ( unsigned i = 0; i < lEdges.size(); ++i )
+      for ( size_t i = 0; i < lEdges.size(); ++i )
       {
         shrinked |= lEdges[i]->SetNewLength2d( surface,F,helper );
       }
       dumpFunctionEnd();
 
       // Move nodes on EDGE's
+      // (XYZ is set as soon as a needed length reached in SetNewLength2d())
       set< _Shrinker1D* >::iterator shr = eShri1D.begin();
       for ( ; shr != eShri1D.end(); ++shr )
         (*shr)->Compute( /*set3D=*/false, helper );
@@ -3632,7 +4115,7 @@ bool _ViscousBuilder::shrink()
       // Smoothing in 2D
       // -----------------
       int nbNoImpSteps = 0;
-      bool moved = true;
+      bool       moved = true;
       badNb = 1;
       while (( nbNoImpSteps < 5 && badNb > 0) && moved)
       {
@@ -3641,11 +4124,10 @@ bool _ViscousBuilder::shrink()
         int oldBadNb = badNb;
         badNb = 0;
         moved = false;
-        for ( unsigned i = 0; i < nodesToSmooth.size(); ++i )
+        for ( size_t i = 0; i < nodesToSmooth.size(); ++i )
         {
           moved |= nodesToSmooth[i].Smooth( badNb,surface,helper,refSign,
-                                            /*isCentroidal=*/isConcaveFace,
-                                            /*set3D=*/isConcaveFace);
+                                            smoothType, /*set3D=*/isConcaveFace);
         }
         if ( badNb < oldBadNb )
           nbNoImpSteps = 0;
@@ -3656,39 +4138,75 @@ bool _ViscousBuilder::shrink()
       }
       if ( badNb > 0 )
         return error(SMESH_Comment("Can't shrink 2D mesh on face ") << f2sd->first );
-    }
+      if ( shriStep > 200 )
+        return error(SMESH_Comment("Infinite loop at shrinking 2D mesh on face ") << f2sd->first );
+
+      // Fix narrow triangles by swapping diagonals
+      // ---------------------------------------
+      if ( toFixTria )
+      {
+        set<const SMDS_MeshNode*> usedNodes;
+        fixBadFaces( F, helper, /*is2D=*/true, shriStep, & usedNodes); // swap diagonals
+
+        // update working data
+        set<const SMDS_MeshNode*>::iterator n;
+        for ( size_t i = 0; i < nodesToSmooth.size() && !usedNodes.empty(); ++i )
+        {
+          n = usedNodes.find( nodesToSmooth[ i ]._node );
+          if ( n != usedNodes.end())
+          {
+            getSimplices( nodesToSmooth[ i ]._node,
+                          nodesToSmooth[ i ]._simplices,
+                          ignoreShapes, NULL,
+                          /*sortSimplices=*/ smoothType == _SmoothNode::ANGULAR );
+            usedNodes.erase( n );
+          }
+        }
+        for ( size_t i = 0; i < lEdges.size() && !usedNodes.empty(); ++i )
+        {
+          n = usedNodes.find( /*tgtNode=*/ lEdges[i]->_nodes.back() );
+          if ( n != usedNodes.end())
+          {
+            getSimplices( lEdges[i]->_nodes.back(),
+                          lEdges[i]->_simplices,
+                          ignoreShapes );
+            usedNodes.erase( n );
+          }
+        }
+      }
+    } // while ( shrinked )
+
     // No wrongly shaped faces remain; final smooth. Set node XYZ.
-    // First, find out a needed quality of smoothing (high for quadrangles only)
-    bool highQuality;
+    bool isStructuredFixed = false;
+    if ( SMESH_2D_Algo* algo = dynamic_cast<SMESH_2D_Algo*>( sm->GetAlgo() ))
+      isStructuredFixed = algo->FixInternalNodes( *data._proxyMesh, F );
+    if ( !isStructuredFixed )
     {
-      const bool hasTria = _mesh->NbTriangles(), hasQuad = _mesh->NbQuadrangles();
-      if ( hasTria != hasQuad ) {
-        highQuality = hasQuad;
-      }
-      else {
-        set<int> nbNodesSet;
-        SMDS_ElemIteratorPtr fIt = smDS->GetElements();
-        while ( fIt->more() && nbNodesSet.size() < 2 )
-          nbNodesSet.insert( fIt->next()->NbCornerNodes() );
-        highQuality = ( *nbNodesSet.begin() == 4 );
+      if ( isConcaveFace ) // fix narrow faces by swapping diagonals
+        fixBadFaces( F, helper, /*is2D=*/false, ++shriStep );
+
+      for ( int st = 3; st; --st )
+      {
+        switch( st ) {
+        case 1: smoothType = _SmoothNode::LAPLACIAN; break;
+        case 2: smoothType = _SmoothNode::LAPLACIAN; break;
+        case 3: smoothType = _SmoothNode::ANGULAR; break;
+        }
+        dumpFunction(SMESH_Comment("shrinkFace")<<f2sd->first<<"_st"<<++smooStep); // debug
+        for ( size_t i = 0; i < nodesToSmooth.size(); ++i )
+        {
+          nodesToSmooth[i].Smooth( badNb,surface,helper,refSign,
+                                   smoothType,/*set3D=*/st==1 );
+        }
+        dumpFunctionEnd();
       }
     }
-    if ( !highQuality && isConcaveFace )
-      fixBadFaces( F, helper ); // fix narrow faces by swaping diagonals
-    for ( int st = highQuality ? 10 : 3; st; --st )
-    {
-      dumpFunction(SMESH_Comment("shrinkFace")<<f2sd->first<<"_st"<<++smooStep); // debug
-      for ( unsigned i = 0; i < nodesToSmooth.size(); ++i )
-        nodesToSmooth[i].Smooth( badNb,surface,helper,refSign,
-                                 /*isCentroidal=*/isConcaveFace,/*set3D=*/st==1 );
-      dumpFunctionEnd();
-    }
     // Set an event listener to clear FACE sub-mesh together with SOLID sub-mesh
-    _SrinkShapeListener::ToClearSubMeshWithSolid( sm, data._solid );
+    VISCOUS_3D::ToClearSubWithMain( sm, data._solid );
 
     if ( !getMeshDS()->IsEmbeddedMode() )
       // Log node movement
-      for ( unsigned i = 0; i < nodesToSmooth.size(); ++i )
+      for ( size_t i = 0; i < nodesToSmooth.size(); ++i )
       {
         SMESH_TNodeXYZ p ( nodesToSmooth[i]._node );
         getMeshDS()->MoveNode( nodesToSmooth[i]._node, p.X(), p.Y(), p.Z() );
@@ -3730,55 +4248,56 @@ bool _ViscousBuilder::prepareEdgeToShrink( _LayerEdge&            edge,
     double uvLen = uvDir.Magnitude();
     uvDir /= uvLen;
     edge._normal.SetCoord( uvDir.X(),uvDir.Y(), 0);
-
-    // IMPORTANT to have src nodes NOT yet REPLACED by tgt nodes in shrinked faces
-    vector<const SMDS_MeshElement*> faces;
-    multimap< double, const SMDS_MeshNode* > proj2node;
-    SMDS_ElemIteratorPtr fIt = srcNode->GetInverseElementIterator(SMDSAbs_Face);
-    while ( fIt->more() )
-    {
-      const SMDS_MeshElement* f = fIt->next();
-      if ( faceSubMesh->Contains( f ))
-        faces.push_back( f );
-    }
-    for ( unsigned i = 0; i < faces.size(); ++i )
-    {
-      const int nbNodes = faces[i]->NbCornerNodes();
-      for ( int j = 0; j < nbNodes; ++j )
-      {
-        const SMDS_MeshNode* n = faces[i]->GetNode(j);
-        if ( n == srcNode ) continue;
-        if ( n->GetPosition()->GetTypeOfPosition() != SMDS_TOP_FACE &&
-             ( faces.size() > 1 || nbNodes > 3 ))
-          continue;
-        gp_Pnt2d uv = helper.GetNodeUV( F, n );
-        gp_Vec2d uvDirN( srcUV, uv );
-        double proj = uvDirN * uvDir;
-        proj2node.insert( make_pair( proj, n ));
-      }
-    }
-
-    multimap< double, const SMDS_MeshNode* >::iterator p2n = proj2node.begin(), p2nEnd;
-    const double minProj = p2n->first;
-    const double projThreshold = 1.1 * uvLen;
-    if ( minProj > projThreshold )
-    {
-      // tgtNode is located so that it does not make faces with wrong orientation
-      return true;
-    }
+    edge._len = uvLen;
+
+    // // IMPORTANT to have src nodes NOT yet REPLACED by tgt nodes in shrinked faces
+    // vector<const SMDS_MeshElement*> faces;
+    // multimap< double, const SMDS_MeshNode* > proj2node;
+    // SMDS_ElemIteratorPtr fIt = srcNode->GetInverseElementIterator(SMDSAbs_Face);
+    // while ( fIt->more() )
+    // {
+    //   const SMDS_MeshElement* f = fIt->next();
+    //   if ( faceSubMesh->Contains( f ))
+    //     faces.push_back( f );
+    // }
+    // for ( size_t i = 0; i < faces.size(); ++i )
+    // {
+    //   const int nbNodes = faces[i]->NbCornerNodes();
+    //   for ( int j = 0; j < nbNodes; ++j )
+    //   {
+    //     const SMDS_MeshNode* n = faces[i]->GetNode(j);
+    //     if ( n == srcNode ) continue;
+    //     if ( n->GetPosition()->GetTypeOfPosition() != SMDS_TOP_FACE &&
+    //          ( faces.size() > 1 || nbNodes > 3 ))
+    //       continue;
+    //     gp_Pnt2d uv = helper.GetNodeUV( F, n );
+    //     gp_Vec2d uvDirN( srcUV, uv );
+    //     double proj = uvDirN * uvDir;
+    //     proj2node.insert( make_pair( proj, n ));
+    //   }
+    // }
+
+    // multimap< double, const SMDS_MeshNode* >::iterator p2n = proj2node.begin(), p2nEnd;
+    // const double       minProj = p2n->first;
+    // const double projThreshold = 1.1 * uvLen;
+    // if ( minProj > projThreshold )
+    // {
+    //   // tgtNode is located so that it does not make faces with wrong orientation
+    //   return true;
+    // }
     edge._pos.resize(1);
     edge._pos[0].SetCoord( tgtUV.X(), tgtUV.Y(), 0 );
 
     // store most risky nodes in _simplices
-    p2nEnd = proj2node.lower_bound( projThreshold );
-    int nbSimpl = ( std::distance( p2n, p2nEnd ) + 1) / 2;
-    edge._simplices.resize( nbSimpl );
-    for ( int i = 0; i < nbSimpl; ++i )
-    {
-      edge._simplices[i]._nPrev = p2n->second;
-      if ( ++p2n != p2nEnd )
-        edge._simplices[i]._nNext = p2n->second;
-    }
+    // p2nEnd = proj2node.lower_bound( projThreshold );
+    // int nbSimpl = ( std::distance( p2n, p2nEnd ) + 1) / 2;
+    // edge._simplices.resize( nbSimpl );
+    // for ( int i = 0; i < nbSimpl; ++i )
+    // {
+    //   edge._simplices[i]._nPrev = p2n->second;
+    //   if ( ++p2n != p2nEnd )
+    //     edge._simplices[i]._nNext = p2n->second;
+    // }
     // set UV of source node to target node
     SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( tgtNode->GetPosition() );
     pos->SetUParameter( srcUV.X() );
@@ -3886,23 +4405,28 @@ bool _ViscousBuilder::prepareEdgeToShrink( _LayerEdge&            edge,
  */
 //================================================================================
 
-void _ViscousBuilder::fixBadFaces(const TopoDS_Face& F, SMESH_MesherHelper& helper)
+void _ViscousBuilder::fixBadFaces(const TopoDS_Face&          F,
+                                  SMESH_MesherHelper&         helper,
+                                  const bool                  is2D,
+                                  const int                   step,
+                                  set<const SMDS_MeshNode*> * involvedNodes)
 {
   SMESH::Controls::AspectRatio qualifier;
   SMESH::Controls::TSequenceOfXYZ points(3), points1(3), points2(3);
-  const double maxAspectRatio = 4.;
+  const double maxAspectRatio = is2D ? 4. : 2;
+  NodeCoordHelper xyz( F, helper, is2D );
 
   // find bad triangles
 
   vector< const SMDS_MeshElement* > badTrias;
   vector< double >                  badAspects;
-  SMESHDS_SubMesh* sm = helper.GetMeshDS()->MeshElements( F );
+  SMESHDS_SubMesh*      sm = helper.GetMeshDS()->MeshElements( F );
   SMDS_ElemIteratorPtr fIt = sm->GetElements();
   while ( fIt->more() )
   {
     const SMDS_MeshElement * f = fIt->next();
     if ( f->NbCornerNodes() != 3 ) continue;
-    for ( int iP = 0; iP < 3; ++iP ) points(iP+1) = SMESH_TNodeXYZ( f->GetNode(iP));
+    for ( int iP = 0; iP < 3; ++iP ) points(iP+1) = xyz( f->GetNode(iP));
     double aspect = qualifier.GetValue( points );
     if ( aspect > maxAspectRatio )
     {
@@ -3910,6 +4434,18 @@ void _ViscousBuilder::fixBadFaces(const TopoDS_Face& F, SMESH_MesherHelper& help
       badAspects.push_back( aspect );
     }
   }
+  if ( step == 1 )
+  {
+    dumpFunction(SMESH_Comment("beforeSwapDiagonals_F")<<helper.GetSubShapeID());
+    SMDS_ElemIteratorPtr fIt = sm->GetElements();
+    while ( fIt->more() )
+    {
+      const SMDS_MeshElement * f = fIt->next();
+      if ( f->NbCornerNodes() == 3 )
+        dumpChangeNodes( f );
+    }
+    dumpFunctionEnd();
+  }
   if ( badTrias.empty() )
     return;
 
@@ -3925,9 +4461,10 @@ void _ViscousBuilder::fixBadFaces(const TopoDS_Face& F, SMESH_MesherHelper& help
     double  aspRatio [3];
     int i1, i2, i3;
 
-    involvedFaces.insert( badTrias[iTia] );
+    if ( !involvedFaces.insert( badTrias[iTia] ).second )
+      continue;
     for ( int iP = 0; iP < 3; ++iP )
-      points(iP+1) = SMESH_TNodeXYZ( badTrias[iTia]->GetNode(iP));
+      points(iP+1) = xyz( badTrias[iTia]->GetNode(iP));
 
     // find triangles adjacent to badTrias[iTia] with better aspect ratio after diag-swaping
     int bestCouple = -1;
@@ -3936,14 +4473,16 @@ void _ViscousBuilder::fixBadFaces(const TopoDS_Face& F, SMESH_MesherHelper& help
       const SMDS_MeshNode* n1 = badTrias[iTia]->GetNode( iSide );
       const SMDS_MeshNode* n2 = badTrias[iTia]->GetNode(( iSide+1 ) % 3 );
       trias [iSide].first  = badTrias[iTia];
-      trias [iSide].second = SMESH_MeshEditor::FindFaceInSet( n1, n2, emptySet, involvedFaces,
-                                                              & i1, & i2 );
-      if ( ! trias[iSide].second || trias[iSide].second->NbCornerNodes() != 3 )
+      trias [iSide].second = SMESH_MeshAlgos::FindFaceInSet( n1, n2, emptySet, involvedFaces,
+                                                             & i1, & i2 );
+      if (( ! trias[iSide].second ) ||
+          ( trias[iSide].second->NbCornerNodes() != 3 ) ||
+          ( ! sm->Contains( trias[iSide].second )))
         continue;
 
       // aspect ratio of an adjacent tria
       for ( int iP = 0; iP < 3; ++iP )
-        points2(iP+1) = SMESH_TNodeXYZ( trias[iSide].second->GetNode(iP));
+        points2(iP+1) = xyz( trias[iSide].second->GetNode(iP));
       double aspectInit = qualifier.GetValue( points2 );
 
       // arrange nodes as after diag-swaping
@@ -3960,6 +4499,12 @@ void _ViscousBuilder::fixBadFaces(const TopoDS_Face& F, SMESH_MesherHelper& help
       if ( aspRatio[ iSide ] > aspectInit + badAspects[ iTia ] )
         continue;
 
+      // prevent inversion of a triangle
+      gp_Vec norm1 = gp_Vec( points1(1), points1(3) ) ^ gp_Vec( points1(1), points1(2) );
+      gp_Vec norm2 = gp_Vec( points2(1), points2(3) ) ^ gp_Vec( points2(1), points2(2) );
+      if ( norm1 * norm2 < 0. && norm1.Angle( norm2 ) > 70./180.*M_PI )
+        continue;
+
       if ( bestCouple < 0 || aspRatio[ bestCouple ] > aspRatio[ iSide ] )
         bestCouple = iSide;
     }
@@ -3980,17 +4525,25 @@ void _ViscousBuilder::fixBadFaces(const TopoDS_Face& F, SMESH_MesherHelper& help
   // swap diagonals
 
   SMESH_MeshEditor editor( helper.GetMesh() );
-  dumpFunction(SMESH_Comment("beforeSwapDiagonals_F")<<helper.GetSubShapeID());
+  dumpFunction(SMESH_Comment("beforeSwapDiagonals_F")<<helper.GetSubShapeID()<<"_"<<step);
   for ( size_t i = 0; i < triaCouples.size(); ++i )
   {
     dumpChangeNodes( triaCouples[i].first );
     dumpChangeNodes( triaCouples[i].second );
     editor.InverseDiag( triaCouples[i].first, triaCouples[i].second );
   }
-  dumpFunctionEnd();
+
+  if ( involvedNodes )
+    for ( size_t i = 0; i < triaCouples.size(); ++i )
+    {
+      involvedNodes->insert( triaCouples[i].first->begin_nodes(),
+                             triaCouples[i].first->end_nodes() );
+      involvedNodes->insert( triaCouples[i].second->begin_nodes(),
+                             triaCouples[i].second->end_nodes() );
+    }
 
   // just for debug dump resulting triangles
-  dumpFunction(SMESH_Comment("swapDiagonals_F")<<helper.GetSubShapeID());
+  dumpFunction(SMESH_Comment("swapDiagonals_F")<<helper.GetSubShapeID()<<"_"<<step);
   for ( size_t i = 0; i < triaCouples.size(); ++i )
   {
     dumpChangeNodes( triaCouples[i].first );
@@ -4016,39 +4569,40 @@ bool _LayerEdge::SetNewLength2d( Handle(Geom_Surface)& surface,
   if ( _sWOL.ShapeType() == TopAbs_FACE )
   {
     gp_XY    curUV = helper.GetNodeUV( F, tgtNode );
-    gp_Pnt2d tgtUV( _pos[0].X(), _pos[0].Y());
+    gp_Pnt2d tgtUV( _pos[0].X(), _pos[0].Y() );
     gp_Vec2d uvDir( _normal.X(), _normal.Y() );
     const double uvLen = tgtUV.Distance( curUV );
+    const double kSafe = Max( 0.5, 1. - 0.1 * _simplices.size() );
 
     // Select shrinking step such that not to make faces with wrong orientation.
-    const double kSafe = 0.8;
-    const double minStepSize = uvLen / 10;
     double stepSize = uvLen;
-    for ( unsigned i = 0; i < _simplices.size(); ++i )
+    for ( size_t i = 0; i < _simplices.size(); ++i )
     {
-      const SMDS_MeshNode* nn[2] = { _simplices[i]._nPrev, _simplices[i]._nNext };
-      for ( int j = 0; j < 2; ++j )
-        if ( const SMDS_MeshNode* n = nn[j] )
-        {
-          gp_XY uv = helper.GetNodeUV( F, n );
-          gp_Vec2d uvDirN( curUV, uv );
-          double proj = uvDirN * uvDir * kSafe;
-          if ( proj < stepSize && proj > minStepSize )
-            stepSize = proj;
-        }
+      // find intersection of 2 lines: curUV-tgtUV and that connecting simplex nodes
+      gp_XY uvN1 = helper.GetNodeUV( F, _simplices[i]._nPrev );
+      gp_XY uvN2 = helper.GetNodeUV( F, _simplices[i]._nNext );
+      gp_XY dirN = uvN2 - uvN1;
+      double det = uvDir.Crossed( dirN );
+      if ( Abs( det )  < std::numeric_limits<double>::min() ) continue;
+      gp_XY dirN2Cur = curUV - uvN1;
+      double step = dirN.Crossed( dirN2Cur ) / det;
+      if ( step > 0 )
+        stepSize = Min( step, stepSize );
     }
-
     gp_Pnt2d newUV;
-    if ( stepSize == uvLen )
+    if ( uvLen - stepSize < _len / 200. )
     {
       newUV = tgtUV;
       _pos.clear();
     }
+    else if ( stepSize > 0 )
+    {
+      newUV = curUV + uvDir.XY() * stepSize * kSafe;
+    }
     else
     {
-      newUV = curUV + uvDir.XY() * stepSize;
+      return true;
     }
-
     SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( tgtNode->GetPosition() );
     pos->SetUParameter( newUV.X() );
     pos->SetVParameter( newUV.Y() );
@@ -4063,22 +4617,22 @@ bool _LayerEdge::SetNewLength2d( Handle(Geom_Surface)& surface,
   {
     TopoDS_Edge E = TopoDS::Edge( _sWOL );
     const SMDS_MeshNode* n2 = _simplices[0]._nPrev;
+    SMDS_EdgePosition* tgtPos = static_cast<SMDS_EdgePosition*>( tgtNode->GetPosition() );
 
     const double u2 = helper.GetNodeU( E, n2, tgtNode );
     const double uSrc   = _pos[0].Coord( U_SRC );
     const double lenTgt = _pos[0].Coord( LEN_TGT );
 
     double newU = _pos[0].Coord( U_TGT );
-    if ( lenTgt < 0.99 * fabs( uSrc-u2 ))
+    if ( lenTgt < 0.99 * fabs( uSrc-u2 )) // n2 got out of src-tgt range
     {
       _pos.clear();
     }
     else
     {
-      newU = 0.1 * uSrc + 0.9 * u2;
+      newU = 0.1 * tgtPos->GetUParameter() + 0.9 * u2;
     }
-    SMDS_EdgePosition* pos = static_cast<SMDS_EdgePosition*>( tgtNode->GetPosition() );
-    pos->SetUParameter( newU );
+    tgtPos->SetUParameter( newU );
 #ifdef __myDEBUG
     gp_XY newUV = helper.GetNodeUV( F, tgtNode, _nodes[0]);
     gp_Pnt p = surface->Value( newUV.X(), newUV.Y() );
@@ -4100,7 +4654,7 @@ bool _SmoothNode::Smooth(int&                  badNb,
                          Handle(Geom_Surface)& surface,
                          SMESH_MesherHelper&   helper,
                          const double          refSign,
-                         bool                  isCentroidal,
+                         SmoothType            how,
                          bool                  set3D)
 {
   const TopoDS_Face& face = TopoDS::Face( helper.GetSubShape() );
@@ -4112,7 +4666,24 @@ bool _SmoothNode::Smooth(int&                  badNb,
 
   // compute new UV for the node
   gp_XY newPos (0,0);
-  if ( isCentroidal && _simplices.size() > 3 )
+  if ( how == TFI && _simplices.size() == 4 )
+  {
+    gp_XY corners[4];
+    for ( size_t i = 0; i < _simplices.size(); ++i )
+      if ( _simplices[i]._nOpp )
+        corners[i] = helper.GetNodeUV( face, _simplices[i]._nOpp, _node );
+      else
+        throw SALOME_Exception(LOCALIZED("TFI smoothing: _Simplex::_nOpp not set!"));
+
+    newPos = helper.calcTFI ( 0.5, 0.5,
+                              corners[0], corners[1], corners[2], corners[3],
+                              uv[1], uv[2], uv[3], uv[0] );
+  }
+  else if ( how == ANGULAR )
+  {
+    newPos = computeAngularPos( uv, helper.GetNodeUV( face, _node ), refSign );
+  }
+  else if ( how == CENTROIDAL && _simplices.size() > 3 )
   {
     // average centers of diagonals wieghted with their reciprocal lengths
     if ( _simplices.size() == 4 )
@@ -4137,13 +4708,12 @@ bool _SmoothNode::Smooth(int&                  badNb,
           newPos += w * ( uv[i]+uv[i2] );
         }
       }
-      newPos /= 2 * sumWeight;
+      newPos /= 2 * sumWeight; // 2 is to get a middle between uv's
     }
   }
   else
   {
     // Laplacian smooth
-    isCentroidal = false;
     for ( size_t i = 0; i < _simplices.size(); ++i )
       newPos += uv[i];
     newPos /= _simplices.size();
@@ -4152,17 +4722,15 @@ bool _SmoothNode::Smooth(int&                  badNb,
   // count quality metrics (orientation) of triangles around the node
   int nbOkBefore = 0;
   gp_XY tgtUV = helper.GetNodeUV( face, _node );
-  for ( unsigned i = 0; i < _simplices.size(); ++i )
+  for ( size_t i = 0; i < _simplices.size(); ++i )
     nbOkBefore += _simplices[i].IsForward( tgtUV, _node, face, helper, refSign );
 
   int nbOkAfter = 0;
-  for ( unsigned i = 0; i < _simplices.size(); ++i )
+  for ( size_t i = 0; i < _simplices.size(); ++i )
     nbOkAfter += _simplices[i].IsForward( newPos, _node, face, helper, refSign );
 
   if ( nbOkAfter < nbOkBefore )
   {
-    // if ( isCentroidal )
-    //   return Smooth( badNb, surface, helper, refSign, !isCentroidal, set3D );
     badNb += _simplices.size() - nbOkBefore;
     return false;
   }
@@ -4185,6 +4753,66 @@ bool _SmoothNode::Smooth(int&                  badNb,
   return ( (tgtUV-newPos).SquareModulus() > 1e-10 );
 }
 
+//================================================================================
+/*!
+ * \brief Computes new UV using angle based smoothing technic
+ */
+//================================================================================
+
+gp_XY _SmoothNode::computeAngularPos(vector<gp_XY>& uv,
+                                     const gp_XY&   uvToFix,
+                                     const double   refSign)
+{
+  uv.push_back( uv.front() );
+
+  vector< gp_XY >  edgeDir ( uv.size() );
+  vector< double > edgeSize( uv.size() );
+  for ( size_t i = 1; i < edgeDir.size(); ++i )
+  {
+    edgeDir [i-1] = uv[i] - uv[i-1];
+    edgeSize[i-1] = edgeDir[i-1].Modulus();
+    if ( edgeSize[i-1] < numeric_limits<double>::min() )
+      edgeDir[i-1].SetX( 100 );
+    else
+      edgeDir[i-1] /= edgeSize[i-1] * refSign;
+  }
+  edgeDir.back()  = edgeDir.front();
+  edgeSize.back() = edgeSize.front();
+
+  gp_XY  newPos(0,0);
+  int    nbEdges = 0;
+  double sumSize = 0;
+  for ( size_t i = 1; i < edgeDir.size(); ++i )
+  {
+    if ( edgeDir[i-1].X() > 1. ) continue;
+    int i1 = i-1;
+    while ( edgeDir[i].X() > 1. && ++i < edgeDir.size() );
+    if ( i == edgeDir.size() ) break;
+    gp_XY p = uv[i];
+    gp_XY norm1( -edgeDir[i1].Y(), edgeDir[i1].X() );
+    gp_XY norm2( -edgeDir[i].Y(),  edgeDir[i].X() );
+    gp_XY bisec = norm1 + norm2;
+    double bisecSize = bisec.Modulus();
+    if ( bisecSize < numeric_limits<double>::min() )
+    {
+      bisec = -edgeDir[i1] + edgeDir[i];
+      bisecSize = bisec.Modulus();
+    }
+    bisec /= bisecSize;
+
+    gp_XY  dirToN  = uvToFix - p;
+    double distToN = dirToN.Modulus();
+    if ( bisec * dirToN < 0 )
+      distToN = -distToN;
+
+    newPos += ( p + bisec * distToN ) * ( edgeSize[i1] + edgeSize[i] );
+    ++nbEdges;
+    sumSize += edgeSize[i1] + edgeSize[i];
+  }
+  newPos /= /*nbEdges * */sumSize;
+  return newPos;
+}
+
 //================================================================================
 /*!
  * \brief Delete _SolidData
@@ -4193,7 +4821,7 @@ bool _SmoothNode::Smooth(int&                  badNb,
 
 _SolidData::~_SolidData()
 {
-  for ( unsigned i = 0; i < _edges.size(); ++i )
+  for ( size_t i = 0; i < _edges.size(); ++i )
   {
     if ( _edges[i] && _edges[i]->_2neibors )
       delete _edges[i]->_2neibors;
@@ -4266,7 +4894,7 @@ void _Shrinker1D::AddEdge( const _LayerEdge* e, SMESH_MesherHelper& helper )
   {
     // remove target node of the _LayerEdge from _nodes
     int nbFound = 0;
-    for ( unsigned i = 0; i < _nodes.size(); ++i )
+    for ( size_t i = 0; i < _nodes.size(); ++i )
       if ( !_nodes[i] || _nodes[i] == tgtNode0 || _nodes[i] == tgtNode1 )
         _nodes[i] = 0, nbFound++;
     if ( nbFound == _nodes.size() )
@@ -4304,7 +4932,7 @@ void _Shrinker1D::Compute(bool set3D, SMESH_MesherHelper& helper)
       l = helper.GetNodeU( E, _edges[1]->_nodes.back(), _nodes.back() );
     double totLen = GCPnts_AbscissaPoint::Length( aCurve, f, l );
 
-    for ( unsigned i = 0; i < _nodes.size(); ++i )
+    for ( size_t i = 0; i < _nodes.size(); ++i )
     {
       if ( !_nodes[i] ) continue;
       double len = totLen * _normPar[i];
@@ -4326,7 +4954,7 @@ void _Shrinker1D::Compute(bool set3D, SMESH_MesherHelper& helper)
     if ( _edges[1] )
       l = helper.GetNodeU( E, _edges[1]->_nodes.back(), _nodes.back() );
     
-    for ( unsigned i = 0; i < _nodes.size(); ++i )
+    for ( size_t i = 0; i < _nodes.size(); ++i )
     {
       if ( !_nodes[i] ) continue;
       double u = f * ( 1-_normPar[i] ) + l * _normPar[i];
@@ -4345,7 +4973,7 @@ void _Shrinker1D::Compute(bool set3D, SMESH_MesherHelper& helper)
 void _Shrinker1D::RestoreParams()
 {
   if ( _done )
-    for ( unsigned i = 0; i < _nodes.size(); ++i )
+    for ( size_t i = 0; i < _nodes.size(); ++i )
     {
       if ( !_nodes[i] ) continue;
       SMDS_EdgePosition* pos = static_cast<SMDS_EdgePosition*>( _nodes[i]->GetPosition() );
@@ -4398,7 +5026,7 @@ bool _ViscousBuilder::addBoundaryElements()
 {
   SMESH_MesherHelper helper( *_mesh );
 
-  for ( unsigned i = 0; i < _sdVec.size(); ++i )
+  for ( size_t i = 0; i < _sdVec.size(); ++i )
   {
     _SolidData& data = _sdVec[i];
     TopTools_IndexedMapOfShape geomEdges;
@@ -4458,7 +5086,7 @@ bool _ViscousBuilder::addBoundaryElements()
         reverse = ( helper.GetSubShapeOri( F, E ) == TopAbs_REVERSED );
         if ( helper.GetSubShapeOri( data._solid, F ) == TopAbs_REVERSED )
           reverse = !reverse, F.Reverse();
-        if ( SMESH_Algo::IsReversedSubMesh( TopoDS::Face(F), getMeshDS() ))
+        if ( helper.IsReversedSubMesh( TopoDS::Face(F) ))
           reverse = !reverse;
       }
       else
@@ -4469,7 +5097,7 @@ bool _ViscousBuilder::addBoundaryElements()
         {
           const TopoDS_Shape* pF = fIt->next();
           if ( helper.IsSubShape( *pF, data._solid) &&
-               !_ignoreShapeIds.count( e2f->first ))
+               !data._ignoreFaceIds.count( e2f->first ))
             F = *pF;
         }
       }
@@ -4485,7 +5113,7 @@ bool _ViscousBuilder::addBoundaryElements()
       // Make faces
       const int dj1 = reverse ? 0 : 1;
       const int dj2 = reverse ? 1 : 0;
-      for ( unsigned j = 1; j < ledges.size(); ++j )
+      for ( size_t j = 1; j < ledges.size(); ++j )
       {
         vector< const SMDS_MeshNode*>&  nn1 = ledges[j-dj1]->_nodes;
         vector< const SMDS_MeshNode*>&  nn2 = ledges[j-dj2]->_nodes;