Salome HOME
23444: EDF 14244 - Viscous Layer
[modules/smesh.git] / src / StdMeshers / StdMeshers_ViscousLayers.cxx
index acb0a56404894fe3b4587200eac20c788f8689c5..76f96eae201a68043f9eea3c03060b53273102af 100644 (file)
@@ -30,6 +30,7 @@
 #include "SMDS_SetIterator.hxx"
 #include "SMESHDS_Group.hxx"
 #include "SMESHDS_Hypothesis.hxx"
+#include "SMESHDS_Mesh.hxx"
 #include "SMESH_Algo.hxx"
 #include "SMESH_ComputeError.hxx"
 #include "SMESH_ControlsDef.hxx"
@@ -94,8 +95,9 @@
 #include <string>
 
 #ifdef _DEBUG_
-//#define __myDEBUG
+#define __myDEBUG
 //#define __NOT_INVALIDATE_BAD_SMOOTH
+//#define __NODES_AT_POS
 #endif
 
 #define INCREMENTAL_SMOOTH // smooth only if min angle is too small
@@ -407,7 +409,7 @@ namespace VISCOUS_3D
 
     gp_XYZ              _normal;    // to boundary of solid
     vector<gp_XYZ>      _pos;       // points computed during inflation
-    double              _len;       // length achived with the last inflation step
+    double              _len;       // length achieved with the last inflation step
     double              _maxLen;    // maximal possible length
     double              _cosin;     // of angle (_normal ^ surface)
     double              _minAngle;  // of _simplices
@@ -423,24 +425,28 @@ namespace VISCOUS_3D
     // data for smoothing of _LayerEdge's based on the EDGE
     _2NearEdges*        _2neibors;
 
-    enum EFlags { TO_SMOOTH       = 1,
-                  MOVED           = 2,   // set by _neibors[i]->SetNewLength()
-                  SMOOTHED        = 4,   // set by this->Smooth()
-                  DIFFICULT       = 8,   // near concave VERTEX
-                  ON_CONCAVE_FACE = 16,
-                  BLOCKED         = 32,  // not to inflate any more
-                  INTERSECTED     = 64,  // close intersection with a face found
-                  NORMAL_UPDATED  = 128,
-                  MARKED          = 256, // local usage
-                  MULTI_NORMAL    = 512, // a normal is invisible by some of surrounding faces
-                  NEAR_BOUNDARY   = 1024,// is near FACE boundary forcing smooth
-                  SMOOTHED_C1     = 2048,// is on _eosC1
-                  DISTORTED       = 4096,// was bad before smoothing
-                  RISKY_SWOL      = 8192 // SWOL is parallel to a source FACE
+    enum EFlags { TO_SMOOTH       = 0x0000001,
+                  MOVED           = 0x0000002, // set by _neibors[i]->SetNewLength()
+                  SMOOTHED        = 0x0000004, // set by _LayerEdge::Smooth()
+                  DIFFICULT       = 0x0000008, // near concave VERTEX
+                  ON_CONCAVE_FACE = 0x0000010,
+                  BLOCKED         = 0x0000020, // not to inflate any more
+                  INTERSECTED     = 0x0000040, // close intersection with a face found
+                  NORMAL_UPDATED  = 0x0000080,
+                  UPD_NORMAL_CONV = 0x0000100, // to update normal on boundary of concave FACE
+                  MARKED          = 0x0000200, // local usage
+                  MULTI_NORMAL    = 0x0000400, // a normal is invisible by some of surrounding faces
+                  NEAR_BOUNDARY   = 0x0000800, // is near FACE boundary forcing smooth
+                  SMOOTHED_C1     = 0x0001000, // is on _eosC1
+                  DISTORTED       = 0x0002000, // was bad before smoothing
+                  RISKY_SWOL      = 0x0004000, // SWOL is parallel to a source FACE
+                  SHRUNK          = 0x0008000, // target node reached a tgt position while shrink()
+                  UNUSED_FLAG     = 0x0100000  // to add user flags after
     };
-    bool Is   ( EFlags f ) const { return _flags & f; }
-    void Set  ( EFlags f ) { _flags |= f; }
-    void Unset( EFlags f ) { _flags &= ~f; }
+    bool Is   ( int flag ) const { return _flags & flag; }
+    void Set  ( int flag ) { _flags |= flag; }
+    void Unset( int flag ) { _flags &= ~flag; }
+    std::string DumpFlags() const; // debug
 
     void SetNewLength( double len, _EdgesOnShape& eos, SMESH_MesherHelper& helper );
     bool SetNewLength2d( Handle(Geom_Surface)& surface,
@@ -456,6 +462,7 @@ namespace VISCOUS_3D
     void ChooseSmooFunction(const set< TGeomID >& concaveVertices,
                             const TNode2Edge&     n2eMap);
     void SmoothPos( const vector< double >& segLen, const double tol );
+    int  GetSmoothedPos( const double tol );
     int  Smooth(const int step, const bool isConcaveFace, bool findBest);
     int  Smooth(const int step, bool findBest, vector< _LayerEdge* >& toSmooth );
     int  CheckNeiborsOnBoundary(vector< _LayerEdge* >* badNeibors = 0, bool * needSmooth = 0 );
@@ -489,10 +496,11 @@ namespace VISCOUS_3D
                            dist, epsilon );
     }
     const gp_XYZ& PrevPos() const { return _pos[ _pos.size() - 2 ]; }
-    const gp_XYZ& PrevCheckPos() const { return _pos[ Is( NORMAL_UPDATED ) ? _pos.size()-2 : 0 ]; }
+    gp_XYZ PrevCheckPos( _EdgesOnShape* eos=0 ) const;
     gp_Ax1 LastSegment(double& segLen, _EdgesOnShape& eos) const;
-    gp_XY  LastUV( const TopoDS_Face& F, _EdgesOnShape& eos ) const;
+    gp_XY  LastUV( const TopoDS_Face& F, _EdgesOnShape& eos, int which=-1 ) const;
     bool   IsOnEdge() const { return _2neibors; }
+    bool   IsOnFace() const { return ( _nodes[0]->GetPosition()->GetDim() == 2 ); }
     gp_XYZ Copy( _LayerEdge& other, _EdgesOnShape& eos, SMESH_MesherHelper& helper );
     void   SetCosin( double cosin );
     void   SetNormal( const gp_XYZ& n ) { _normal = n; }
@@ -658,6 +666,8 @@ namespace VISCOUS_3D
 
     _SolidData*            _data; // parent SOLID
 
+    _LayerEdge*      operator[](size_t i) const { return (_LayerEdge*) _edges[i]; }
+    size_t           size() const { return _edges.size(); }
     TopAbs_ShapeEnum ShapeType() const
     { return _shape.IsNull() ? TopAbs_SHAPE : _shape.ShapeType(); }
     TopAbs_ShapeEnum SWOLType() const
@@ -672,7 +682,7 @@ namespace VISCOUS_3D
 
   //--------------------------------------------------------------------------------
   /*!
-   * \brief Convex FACE whose radius of curvature is less than the thickness of 
+   * \brief Convex FACE whose radius of curvature is less than the thickness of
    *        layers. It is used to detect distortion of prisms based on a convex
    *        FACE and to update normals to enable further increasing the thickness
    */
@@ -686,7 +696,14 @@ namespace VISCOUS_3D
     // map a sub-shape to _SolidData::_edgesOnShape
     map< TGeomID, _EdgesOnShape* >  _subIdToEOS;
 
+    bool                            _isTooCurved;
     bool                            _normalsFixed;
+    bool                            _normalsFixedOnBorders; // used in putOnOffsetSurface()
+
+    double GetMaxCurvature( _SolidData&         data,
+                            _EdgesOnShape&      eof,
+                            BRepLProp_SLProps&  surfProp,
+                            SMESH_MesherHelper& helper);
 
     bool GetCenterOfCurvature( _LayerEdge*         ledge,
                                BRepLProp_SLProps&  surfProp,
@@ -717,6 +734,7 @@ namespace VISCOUS_3D
   {
     typedef const StdMeshers_ViscousLayers* THyp;
     TopoDS_Shape                    _solid;
+    TopTools_MapOfShape             _before; // SOLIDs to be computed before _solid
     TGeomID                         _index; // SOLID id
     _MeshOfSolid*                   _proxyMesh;
     list< THyp >                    _hyps;
@@ -750,8 +768,6 @@ namespace VISCOUS_3D
 
     int                              _nbShapesToSmooth;
 
-    //map< TGeomID,Handle(Geom_Curve)> _edge2curve;
-
     vector< _CollisionEdges >        _collisionEdges;
     set< TGeomID >                   _concaveFaces;
 
@@ -781,10 +797,10 @@ namespace VISCOUS_3D
 
     SMESH_MesherHelper& GetHelper() const { return *_helper; }
 
-    void UnmarkEdges() {
+    void UnmarkEdges( int flag = _LayerEdge::MARKED ) {
       for ( size_t i = 0; i < _edgesOnShape.size(); ++i )
         for ( size_t j = 0; j < _edgesOnShape[i]._edges.size(); ++j )
-          _edgesOnShape[i]._edges[j]->Unset( _LayerEdge::MARKED );
+          _edgesOnShape[i]._edges[j]->Unset( flag );
     }
     void AddShapesToSmooth( const set< _EdgesOnShape* >& shape,
                             const set< _EdgesOnShape* >* edgesNoAnaSmooth=0 );
@@ -805,8 +821,10 @@ namespace VISCOUS_3D
     _OffsetPlane() {
       _isLineOK[0] = _isLineOK[1] = false; _faceIndexNext[0] = _faceIndexNext[1] = -1;
     }
-    void   ComputeIntersectionLine( _OffsetPlane& pln );
-    gp_XYZ GetCommonPoint(bool& isFound) const;
+    void   ComputeIntersectionLine( _OffsetPlane&        pln, 
+                                    const TopoDS_Edge&   E,
+                                    const TopoDS_Vertex& V );
+    gp_XYZ GetCommonPoint(bool& isFound, const TopoDS_Vertex& V) const;
     int    NbLines() const { return _isLineOK[0] + _isLineOK[1]; }
   };
   //--------------------------------------------------------------------------------
@@ -861,6 +879,7 @@ namespace VISCOUS_3D
                             const gp_XY&   uvToFix,
                             const double   refSign );
   };
+  struct PyDump;
   //--------------------------------------------------------------------------------
   /*!
    * \brief Builder of viscous layers
@@ -885,6 +904,7 @@ namespace VISCOUS_3D
   private:
 
     bool findSolidsWithLayers();
+    bool setBefore( _SolidData& solidBefore, _SolidData& solidAfter );
     bool findFacesWithLayers(const bool onlyWith=false);
     void getIgnoreFaces(const TopoDS_Shape&             solid,
                         const StdMeshers_ViscousLayers* hyp,
@@ -906,7 +926,8 @@ namespace VISCOUS_3D
     gp_XYZ getWeigthedNormal( const _LayerEdge*                edge );
     gp_XYZ getNormalByOffset( _LayerEdge*                      edge,
                               std::pair< TopoDS_Face, gp_XYZ > fId2Normal[],
-                              int                              nbFaces );
+                              int                              nbFaces,
+                              bool                             lastNoOffset = false);
     bool findNeiborsOnEdge(const _LayerEdge*     edge,
                            const SMDS_MeshNode*& n1,
                            const SMDS_MeshNode*& n2,
@@ -929,8 +950,17 @@ namespace VISCOUS_3D
                               vector< _EdgesOnShape* >& eosC1,
                               const int                 infStep );
     void makeOffsetSurface( _EdgesOnShape& eos, SMESH_MesherHelper& );
-    void putOnOffsetSurface( _EdgesOnShape& eos, int infStep, int smooStep=0, bool moveAll=false );
+    void putOnOffsetSurface( _EdgesOnShape& eos, int infStep,
+                             vector< _EdgesOnShape* >& eosC1,
+                             int smooStep=0, int moveAll=false );
     void findCollisionEdges( _SolidData& data, SMESH_MesherHelper& helper );
+    void findEdgesToUpdateNormalNearConvexFace( _ConvexFace &       convFace,
+                                                _SolidData&         data,
+                                                SMESH_MesherHelper& helper );
+    void limitMaxLenByCurvature( _SolidData& data, SMESH_MesherHelper& helper );
+    void limitMaxLenByCurvature( _LayerEdge* e1, _LayerEdge* e2,
+                                 _EdgesOnShape& eos1, _EdgesOnShape& eos2,
+                                 SMESH_MesherHelper& helper );
     bool updateNormals( _SolidData& data, SMESH_MesherHelper& helper, int stepNb, double stepSize );
     bool updateNormalsOfConvexFaces( _SolidData&         data,
                                      SMESH_MesherHelper& helper,
@@ -944,7 +974,7 @@ namespace VISCOUS_3D
                         _LayerEdge&   edge,
                         const gp_XYZ& newNormal);
     bool refine(_SolidData& data);
-    bool shrink();
+    bool shrink(_SolidData& data);
     bool prepareEdgeToShrink( _LayerEdge& edge, _EdgesOnShape& eos,
                               SMESH_MesherHelper& helper,
                               const SMESHDS_SubMesh* faceSubMesh );
@@ -954,7 +984,7 @@ namespace VISCOUS_3D
                      const bool                  is2D,
                      const int                   step,
                      set<const SMDS_MeshNode*> * involvedNodes=NULL);
-    bool addBoundaryElements();
+    bool addBoundaryElements(_SolidData& data);
 
     bool error( const string& text, int solidID=-1 );
     SMESHDS_Mesh* getMeshDS() const { return _mesh->GetMeshDS(); }
@@ -962,11 +992,15 @@ namespace VISCOUS_3D
     // debug
     void makeGroupOfLE();
 
-    SMESH_Mesh*           _mesh;
-    SMESH_ComputeErrorPtr _error;
+    SMESH_Mesh*                _mesh;
+    SMESH_ComputeErrorPtr      _error;
 
-    vector< _SolidData >  _sdVec;
-    int                   _tmpFaceID;
+    vector<                    _SolidData >  _sdVec;
+    TopTools_IndexedMapOfShape _solids; // to find _SolidData by a solid
+    TopTools_MapOfShape        _shrinkedFaces;
+
+    int                        _tmpFaceID;
+    PyDump*                    _pyDump;
   };
   //--------------------------------------------------------------------------------
   /*!
@@ -1003,15 +1037,18 @@ namespace VISCOUS_3D
       double      _len;    // length reached at previous inflation step
       double      _param;  // on EDGE
       _2NearEdges _2edges; // 2 neighbor _LayerEdge's
+      gp_XYZ      _edgeDir;// EDGE tangent at _param
       double Distance( const OffPnt& p ) const { return ( _xyz - p._xyz ).Modulus(); }
     };
     vector< OffPnt >   _offPoints;
     vector< double >   _leParams; // normalized param of _eos._edges on EDGE
     Handle(Geom_Curve) _anaCurve; // for analytic smooth
     _LayerEdge         _leOnV[2]; // _LayerEdge's holding normal to the EDGE at VERTEXes
+    gp_XYZ             _edgeDir[2]; // tangent at VERTEXes
     size_t             _iSeg[2];  // index of segment where extreme tgt node is projected
     _EdgesOnShape&     _eos;
     double             _curveLen; // length of the EDGE
+    std::pair<int,int> _eToSmooth[2]; // <from,to> indices of _LayerEdge's in _eos
 
     static Handle(Geom_Curve) CurveForSmooth( const TopoDS_Edge&  E,
                                               _EdgesOnShape&      eos,
@@ -1025,31 +1062,24 @@ namespace VISCOUS_3D
     bool Perform(_SolidData&                    data,
                  Handle(ShapeAnalysis_Surface)& surface,
                  const TopoDS_Face&             F,
-                 SMESH_MesherHelper&            helper )
-    {
-      if ( _leParams.empty() || ( !isAnalytic() && _offPoints.empty() ))
-        prepare( data );
+                 SMESH_MesherHelper&            helper );
 
-      if ( isAnalytic() )
-        return smoothAnalyticEdge( data, surface, F, helper );
-      else
-        return smoothComplexEdge ( data, surface, F, helper );
-    }
     void prepare(_SolidData& data );
 
+    void findEdgesToSmooth();
+
+    bool isToSmooth( int iE );
+
     bool smoothAnalyticEdge( _SolidData&                    data,
                              Handle(ShapeAnalysis_Surface)& surface,
                              const TopoDS_Face&             F,
                              SMESH_MesherHelper&            helper);
-
     bool smoothComplexEdge( _SolidData&                    data,
                             Handle(ShapeAnalysis_Surface)& surface,
                             const TopoDS_Face&             F,
                             SMESH_MesherHelper&            helper);
-
-    void setNormalOnV( const bool          is2nd,
-                       SMESH_MesherHelper& helper);
-
+    gp_XYZ getNormalNormal( const gp_XYZ & normal,
+                            const gp_XYZ&  edgeDir);
     _LayerEdge* getLEdgeOnV( bool is2nd )
     {
       return _eos._edges[ is2nd ? _eos._edges.size()-1 : 0 ]->_2neibors->_edges[ is2nd ];
@@ -1217,8 +1247,8 @@ StdMeshers_ViscousLayers::Compute(SMESH_Mesh&         theMesh,
                                   const bool          toMakeN2NMap) const
 {
   using namespace VISCOUS_3D;
-  _ViscousBuilder bulder;
-  SMESH_ComputeErrorPtr err = bulder.Compute( theMesh, theShape );
+  _ViscousBuilder builder;
+  SMESH_ComputeErrorPtr err = builder.Compute( theMesh, theShape );
   if ( err && !err->IsOK() )
     return SMESH_ProxyMesh::Ptr();
 
@@ -1230,7 +1260,7 @@ StdMeshers_ViscousLayers::Compute(SMESH_Mesh&         theMesh,
          _ViscousListener::GetSolidMesh( &theMesh, exp.Current(), /*toCreate=*/false))
     {
       if ( toMakeN2NMap && !pm->_n2nMapComputed )
-        if ( !bulder.MakeN2NMap( pm ))
+        if ( !builder.MakeN2NMap( pm ))
           return SMESH_ProxyMesh::Ptr();
       components.push_back( SMESH_ProxyMesh::Ptr( pm ));
       pm->myIsDeletable = false; // it will de deleted by boost::shared_ptr
@@ -1294,8 +1324,8 @@ StdMeshers_ViscousLayers::CheckHypothesis(SMESH_Mesh&                          t
                                           const TopoDS_Shape&                  theShape,
                                           SMESH_Hypothesis::Hypothesis_Status& theStatus)
 {
-  VISCOUS_3D::_ViscousBuilder bulder;
-  SMESH_ComputeErrorPtr err = bulder.CheckHypotheses( theMesh, theShape );
+  VISCOUS_3D::_ViscousBuilder builder;
+  SMESH_ComputeErrorPtr err = builder.CheckHypotheses( theMesh, theShape );
   if ( err && !err->IsOK() )
     theStatus = SMESH_Hypothesis::HYP_INCOMPAT_HYPS;
   else
@@ -1641,14 +1671,15 @@ namespace VISCOUS_3D
   // HOWTO use: run python commands written in a console to see
   //  construction steps of viscous layers
 #ifdef __myDEBUG
-  ofstream* py;
-  int       theNbPyFunc;
-  struct PyDump {
+  ostream* py;
+  int      theNbPyFunc;
+  struct PyDump
+  {
     PyDump(SMESH_Mesh& m) {
       int tag = 3 + m.GetId();
       const char* fname = "/tmp/viscous.py";
       cout << "execfile('"<<fname<<"')"<<endl;
-      py = new ofstream(fname);
+      py = _pyStream = new ofstream(fname);
       *py << "import SMESH" << endl
           << "from salome.smesh import smeshBuilder" << endl
           << "smesh  = smeshBuilder.New(salome.myStudy)" << endl
@@ -1666,6 +1697,14 @@ namespace VISCOUS_3D
       delete py; py=0;
     }
     ~PyDump() { Finish(); cout << "NB FUNCTIONS: " << theNbPyFunc << endl; }
+    struct MyStream : public ostream
+    {
+      template <class T> ostream & operator<<( const T &anything ) { return *this ; }
+    };
+    void Pause() { py = &_mystream; }
+    void Resume() { py = _pyStream; }
+    MyStream _mystream;
+    ostream* _pyStream;
   };
 #define dumpFunction(f) { _dumpFunction(f, __LINE__);}
 #define dumpMove(n)     { _dumpMove(n, __LINE__);}
@@ -1688,7 +1727,7 @@ namespace VISCOUS_3D
 
 #else
 
-  struct PyDump { PyDump(SMESH_Mesh&) {} void Finish() {} };
+  struct PyDump { PyDump(SMESH_Mesh&) {} void Finish() {} void Pause() {} void Resume() {} };
 #define dumpFunction(f) f
 #define dumpMove(n)
 #define dumpMoveComm(n,txt)
@@ -1819,8 +1858,6 @@ bool _ViscousBuilder::MakeN2NMap( _MeshOfSolid* pm )
 SMESH_ComputeErrorPtr _ViscousBuilder::Compute(SMESH_Mesh&         theMesh,
                                                const TopoDS_Shape& theShape)
 {
-  // TODO: set priority of solids during Gen::Compute()
-
   _mesh = & theMesh;
 
   // check if proxy mesh already computed
@@ -1832,6 +1869,7 @@ SMESH_ComputeErrorPtr _ViscousBuilder::Compute(SMESH_Mesh&         theMesh,
     return SMESH_ComputeErrorPtr(); // everything already computed
 
   PyDump debugDump( theMesh );
+  _pyDump = &debugDump;
 
   // TODO: ignore already computed SOLIDs 
   if ( !findSolidsWithLayers())
@@ -1842,22 +1880,37 @@ SMESH_ComputeErrorPtr _ViscousBuilder::Compute(SMESH_Mesh&         theMesh,
 
   for ( size_t i = 0; i < _sdVec.size(); ++i )
   {
-    if ( ! makeLayer(_sdVec[i]) )
+    size_t iSD = 0;
+    for ( iSD = 0; iSD < _sdVec.size(); ++iSD ) // find next SOLID to compute
+      if ( _sdVec[iSD]._before.IsEmpty() &&
+           !_sdVec[iSD]._solid.IsNull() &&
+           _sdVec[iSD]._n2eMap.empty() )
+        break;
+
+    if ( ! makeLayer(_sdVec[iSD]) )   // create _LayerEdge's
       return _error;
 
-    if ( _sdVec[i]._n2eMap.size() == 0 )
+    if ( _sdVec[iSD]._n2eMap.size() == 0 ) // no layers in a SOLID
+    {
+      _sdVec[iSD]._solid.Nullify();
       continue;
-    
-    if ( ! inflate(_sdVec[i]) )
+    }
+
+    if ( ! inflate(_sdVec[iSD]) )     // increase length of _LayerEdge's
       return _error;
 
-    if ( ! refine(_sdVec[i]) )
+    if ( ! refine(_sdVec[iSD]) )      // create nodes and prisms
       return _error;
-  }
-  if ( !shrink() )
-    return _error;
 
-  addBoundaryElements();
+    if ( ! shrink(_sdVec[iSD]) )      // shrink 2D mesh on FACEs w/o layer
+      return _error;
+
+    addBoundaryElements(_sdVec[iSD]); // create quadrangles on prism bare sides
+
+    const TopoDS_Shape& solid = _sdVec[iSD]._solid;
+    for ( iSD = 0; iSD < _sdVec.size(); ++iSD )
+      _sdVec[iSD]._before.Remove( solid );
+  }
 
   makeGroupOfLE(); // debug
   debugDump.Finish();
@@ -1906,12 +1959,14 @@ bool _ViscousBuilder::findSolidsWithLayers()
   TopExp::MapShapes( _mesh->GetShapeToMesh(), TopAbs_SOLID, allSolids );
   _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
-    SMESH_Algo* algo = gen->GetAlgo( *_mesh, allSolids(i) );
+    SMESH_subMesh* sm = _mesh->GetSubMesh( allSolids(i) );
+    if ( sm->GetSubMeshDS() && sm->GetSubMeshDS()->NbElements() > 0 )
+      continue; // solid is already meshed
+    SMESH_Algo* algo = sm->GetAlgo();
     if ( !algo ) continue;
     // TODO: check if algo is hidden
     const list <const SMESHDS_Hypothesis *> & allHyps =
@@ -1936,6 +1991,7 @@ bool _ViscousBuilder::findSolidsWithLayers()
           soData->_index = getMeshDS()->ShapeToIndex( allSolids(i));
           soData->_helper = new SMESH_MesherHelper( *_mesh );
           soData->_helper->SetSubShape( allSolids(i) );
+          _solids.Add( allSolids(i) );
         }
         soData->_hyps.push_back( viscHyp );
         soData->_hypShapes.push_back( hypShape );
@@ -1950,7 +2006,37 @@ bool _ViscousBuilder::findSolidsWithLayers()
 
 //================================================================================
 /*!
- * \brief 
+ * \brief Set a _SolidData to be computed before another
+ */
+//================================================================================
+
+bool _ViscousBuilder::setBefore( _SolidData& solidBefore, _SolidData& solidAfter )
+{
+  // check possibility to set this order; get all solids before solidBefore
+  TopTools_IndexedMapOfShape allSolidsBefore;
+  allSolidsBefore.Add( solidBefore._solid );
+  for ( int i = 1; i <= allSolidsBefore.Extent(); ++i )
+  {
+    int iSD = _solids.FindIndex( allSolidsBefore(i) );
+    if ( iSD )
+    {
+      TopTools_MapIteratorOfMapOfShape soIt( _sdVec[ iSD-1 ]._before );
+      for ( ; soIt.More(); soIt.Next() )
+        allSolidsBefore.Add( soIt.Value() );
+    }
+  }
+  if ( allSolidsBefore.Contains( solidAfter._solid ))
+    return false;
+
+  for ( int i = 1; i <= allSolidsBefore.Extent(); ++i )
+    solidAfter._before.Add( allSolidsBefore(i) );
+
+  return true;
+}
+
+//================================================================================
+/*!
+ * \brief
  */
 //================================================================================
 
@@ -1958,13 +2044,10 @@ bool _ViscousBuilder::findFacesWithLayers(const bool onlyWith)
 {
   SMESH_MesherHelper helper( *_mesh );
   TopExp_Explorer exp;
-  TopTools_IndexedMapOfShape solids;
 
   // collect all faces-to-ignore defined by hyp
   for ( size_t i = 0; i < _sdVec.size(); ++i )
   {
-    solids.Add( _sdVec[i]._solid );
-
     // get faces-to-ignore defined by each hyp
     typedef const StdMeshers_ViscousLayers* THyp;
     typedef std::pair< set<TGeomID>, THyp > TFacesOfHyp;
@@ -2060,8 +2143,9 @@ bool _ViscousBuilder::findFacesWithLayers(const bool onlyWith)
     }
   }
 
-  // Find faces to shrink mesh on (solution 2 in issue 0020832);
+  // Find FACEs to shrink mesh on (solution 2 in issue 0020832): fill in _shrinkShape2Shape
   TopTools_IndexedMapOfShape shapes;
+  std::string structAlgoName = "Hexa_3D";
   for ( size_t i = 0; i < _sdVec.size(); ++i )
   {
     shapes.Clear();
@@ -2069,121 +2153,32 @@ bool _ViscousBuilder::findFacesWithLayers(const bool onlyWith)
     for ( int iE = 1; iE <= shapes.Extent(); ++iE )
     {
       const TopoDS_Shape& edge = shapes(iE);
-      // find 2 faces sharing an edge
+      // find 2 FACEs sharing an EDGE
       TopoDS_Shape FF[2];
-      PShapeIteratorPtr fIt = helper.GetAncestors(edge, *_mesh, TopAbs_FACE);
+      PShapeIteratorPtr fIt = helper.GetAncestors(edge, *_mesh, TopAbs_FACE, &_sdVec[i]._solid);
       while ( fIt->more())
       {
         const TopoDS_Shape* f = fIt->next();
-        if ( helper.IsSubShape( *f, _sdVec[i]._solid))
-          FF[ int( !FF[0].IsNull()) ] = *f;
+        FF[ int( !FF[0].IsNull()) ] = *f;
       }
       if( FF[1].IsNull() ) continue; // seam edge can be shared by 1 FACE only
+
       // check presence of layers on them
       int ignore[2];
       for ( int j = 0; j < 2; ++j )
-        ignore[j] = _sdVec[i]._ignoreFaceIds.count ( getMeshDS()->ShapeToIndex( FF[j] ));
+        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
-      bool collision = false;
-      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]._noShrinkShapes.insert( iFace );
-            //fWOL.Nullify();
-            collision = true;
-          }
-        }
-      // add edge to maps
+
+      // add EDGE to maps
       if ( !fWOL.IsNull())
       {
         TGeomID edgeInd = getMeshDS()->ShapeToIndex( edge );
         _sdVec[i]._shrinkShape2Shape.insert( make_pair( edgeInd, fWOL ));
-        if ( collision )
-        {
-          // _shrinkShape2Shape will be used to temporary inflate _LayerEdge's based
-          // on the edge but shrink won't be performed
-          _sdVec[i]._noShrinkShapes.insert( edgeInd );
-        }
       }
     }
   }
-  // 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 ( size_t i = 0; i < _sdVec.size(); ++i )
-  {
-    map< TGeomID, TopoDS_Shape >::iterator e2f = _sdVec[i]._shrinkShape2Shape.begin();
-    for ( ; e2f != _sdVec[i]._shrinkShape2Shape.end(); ++e2f )
-    {
-      const TopoDS_Shape& fWOL = e2f->second;
-      const TGeomID     edgeID = e2f->first;
-      bool notShrinkFace = false;
-      PShapeIteratorPtr soIt = helper.GetAncestors(fWOL, *_mesh, TopAbs_SOLID);
-      while ( soIt->more() )
-      {
-        const TopoDS_Shape* solid = soIt->next();
-        if ( _sdVec[i]._solid.IsSame( *solid )) continue;
-        SMESH_Algo* algo = _mesh->GetGen()->GetAlgo( *_mesh, *solid );
-        if ( !algo || !notSupportAlgos.count( algo->GetName() )) continue;
-        notShrinkFace = true;
-        size_t iSolid = 0;
-        for ( ; iSolid < _sdVec.size(); ++iSolid )
-        {
-          if ( _sdVec[iSolid]._solid.IsSame( *solid ) ) {
-            if ( _sdVec[iSolid]._shrinkShape2Shape.count( edgeID ))
-              notShrinkFace = false;
-            break;
-          }
-        }
-        if ( notShrinkFace )
-        {
-          _sdVec[i]._noShrinkShapes.insert( edgeID );
-
-          // add VERTEXes of the edge in _noShrinkShapes
-          TopoDS_Shape edge = getMeshDS()->IndexToShape( edgeID );
-          for ( TopoDS_Iterator vIt( edge ); vIt.More(); vIt.Next() )
-            _sdVec[i]._noShrinkShapes.insert( getMeshDS()->ShapeToIndex( vIt.Value() ));
-
-          // check if there is a collision with to-shrink-from EDGEs in iSolid
-          if ( iSolid == _sdVec.size() )
-            continue; // no VL in the solid
-          shapes.Clear();
-          TopExp::MapShapes( fWOL, TopAbs_EDGE, shapes);
-          for ( int iE = 1; iE <= shapes.Extent(); ++iE )
-          {
-            const TopoDS_Edge& E = TopoDS::Edge( shapes( iE ));
-            const TGeomID    eID = getMeshDS()->ShapeToIndex( E );
-            if ( eID == edgeID ||
-                 !_sdVec[iSolid]._shrinkShape2Shape.count( eID ) ||
-                 _sdVec[i]._noShrinkShapes.count( eID ))
-              continue;
-            for ( int is1st = 0; is1st < 2; ++is1st )
-            {
-              TopoDS_Vertex V = helper.IthVertex( is1st, E );
-              if ( _sdVec[i]._noShrinkShapes.count( getMeshDS()->ShapeToIndex( V ) ))
-              {
-                // _sdVec[i]._noShrinkShapes.insert( eID );
-                // V = helper.IthVertex( !is1st, E );
-                // _sdVec[i]._noShrinkShapes.insert( getMeshDS()->ShapeToIndex( V ));
-                //iE = 0; // re-start the loop on EDGEs of fWOL
-                return error("No way to make a conformal mesh with "
-                             "the given set of faces with layers", _sdVec[i]._index);
-              }
-            }
-          }
-        }
-
-      } // while ( soIt->more() )
-    } // loop on _sdVec[i]._shrinkShape2Shape
-  } // loop on _sdVec to fill in _SolidData::_noShrinkShapes
 
   // Find the SHAPE along which to inflate _LayerEdge based on VERTEX
 
@@ -2197,18 +2192,14 @@ bool _ViscousBuilder::findFacesWithLayers(const bool onlyWith)
       // find faces WOL sharing the vertex
       vector< TopoDS_Shape > facesWOL;
       size_t totalNbFaces = 0;
-      PShapeIteratorPtr fIt = helper.GetAncestors(vertex, *_mesh, TopAbs_FACE);
+      PShapeIteratorPtr fIt = helper.GetAncestors(vertex, *_mesh, TopAbs_FACE, &_sdVec[i]._solid );
       while ( fIt->more())
       {
         const TopoDS_Shape* f = fIt->next();
-        if ( helper.IsSubShape( *f, _sdVec[i]._solid ) )
-        {
-          totalNbFaces++;
-          const int fID = getMeshDS()->ShapeToIndex( *f );
-          if ( _sdVec[i]._ignoreFaceIds.count ( fID ) /*&&
-               !_sdVec[i]._noShrinkShapes.count( fID )*/)
-            facesWOL.push_back( *f );
-        }
+        totalNbFaces++;
+        const int fID = getMeshDS()->ShapeToIndex( *f );
+        if ( _sdVec[i]._ignoreFaceIds.count ( fID ) /*&& !_sdVec[i]._noShrinkShapes.count( fID )*/)
+          facesWOL.push_back( *f );
       }
       if ( facesWOL.size() == totalNbFaces || facesWOL.empty() )
         continue; // no layers at this vertex or no WOL
@@ -2258,7 +2249,130 @@ bool _ViscousBuilder::findFacesWithLayers(const bool onlyWith)
     }
   }
 
-  // add FACEs of other SOLIDs to _ignoreFaceIds
+  // Add to _noShrinkShapes sub-shapes of FACE's that can't be shrinked since
+  // the algo of the SOLID sharing the FACE does not support it or for other reasons
+  set< string > notSupportAlgos; notSupportAlgos.insert( structAlgoName );
+  for ( size_t i = 0; i < _sdVec.size(); ++i )
+  {
+    map< TGeomID, TopoDS_Shape >::iterator e2f = _sdVec[i]._shrinkShape2Shape.begin();
+    for ( ; e2f != _sdVec[i]._shrinkShape2Shape.end(); ++e2f )
+    {
+      const TopoDS_Shape& fWOL = e2f->second;
+      const TGeomID     edgeID = e2f->first;
+      TGeomID           faceID = getMeshDS()->ShapeToIndex( fWOL );
+      TopoDS_Shape        edge = getMeshDS()->IndexToShape( edgeID );
+      if ( edge.ShapeType() != TopAbs_EDGE )
+        continue; // shrink shape is VERTEX
+
+      TopoDS_Shape solid;
+      PShapeIteratorPtr soIt = helper.GetAncestors(fWOL, *_mesh, TopAbs_SOLID);
+      while ( soIt->more() && solid.IsNull() )
+      {
+        const TopoDS_Shape* so = soIt->next();
+        if ( !so->IsSame( _sdVec[i]._solid ))
+          solid = *so;
+      }
+      if ( solid.IsNull() )
+        continue;
+
+      bool noShrinkE = false;
+      SMESH_Algo*  algo = _mesh->GetSubMesh( solid )->GetAlgo();
+      bool isStructured = ( algo && algo->GetName() == structAlgoName );
+      size_t     iSolid = _solids.FindIndex( solid ) - 1;
+      if ( iSolid < _sdVec.size() && _sdVec[ iSolid ]._ignoreFaceIds.count( faceID ))
+      {
+        // the adjacent SOLID has NO layers on fWOL;
+        // shrink allowed if
+        // - there are layers on the EDGE in the adjacent SOLID
+        // - there are NO layers in the adjacent SOLID && algo is unstructured and computed later
+        bool hasWLAdj = (_sdVec[iSolid]._shrinkShape2Shape.count( edgeID ));
+        bool shrinkAllowed = (( hasWLAdj ) ||
+                              ( !isStructured && setBefore( _sdVec[ i ], _sdVec[ iSolid ] )));
+        noShrinkE = !shrinkAllowed;
+      }
+      else if ( iSolid < _sdVec.size() )
+      {
+        // the adjacent SOLID has layers on fWOL;
+        // check if SOLID's mesh is unstructured and then try to set it
+        // to be computed after the i-th solid
+        if ( isStructured || !setBefore( _sdVec[ i ], _sdVec[ iSolid ] ))
+          noShrinkE = true; // don't shrink fWOL
+      }
+      else
+      {
+        // the adjacent SOLID has NO layers at all
+        noShrinkE = isStructured;
+      }
+
+      if ( noShrinkE )
+      {
+        _sdVec[i]._noShrinkShapes.insert( edgeID );
+
+        // check if there is a collision with to-shrink-from EDGEs in iSolid
+        // if ( iSolid < _sdVec.size() )
+        // {
+        //   shapes.Clear();
+        //   TopExp::MapShapes( fWOL, TopAbs_EDGE, shapes);
+        //   for ( int iE = 1; iE <= shapes.Extent(); ++iE )
+        //   {
+        //     const TopoDS_Edge& E = TopoDS::Edge( shapes( iE ));
+        //     const TGeomID    eID = getMeshDS()->ShapeToIndex( E );
+        //     if ( eID == edgeID ||
+        //          !_sdVec[iSolid]._shrinkShape2Shape.count( eID ) ||
+        //          _sdVec[i]._noShrinkShapes.count( eID ))
+        //       continue;
+        //     for ( int is1st = 0; is1st < 2; ++is1st )
+        //     {
+        //       TopoDS_Vertex V = helper.IthVertex( is1st, E );
+        //       if ( _sdVec[i]._noShrinkShapes.count( getMeshDS()->ShapeToIndex( V ) ))
+        //       {
+        //         return error("No way to make a conformal mesh with "
+        //                      "the given set of faces with layers", _sdVec[i]._index);
+        //       }
+        //     }
+        //   }
+        // }
+      }
+
+      // add VERTEXes of the edge in _noShrinkShapes, which is necessary if
+      // _shrinkShape2Shape is different in the adjacent SOLID
+      for ( TopoDS_Iterator vIt( edge ); vIt.More(); vIt.Next() )
+      {
+        TGeomID vID = getMeshDS()->ShapeToIndex( vIt.Value() );
+        bool noShrinkV = false;
+
+        if ( iSolid < _sdVec.size() )
+        {
+          if ( _sdVec[ iSolid ]._ignoreFaceIds.count( faceID ))
+          {
+            map< TGeomID, TopoDS_Shape >::iterator i2S, i2SAdj;
+            i2S    = _sdVec[i     ]._shrinkShape2Shape.find( vID );
+            i2SAdj = _sdVec[iSolid]._shrinkShape2Shape.find( vID );
+            if ( i2SAdj == _sdVec[iSolid]._shrinkShape2Shape.end() )
+              noShrinkV = ( i2S->second.ShapeType() == TopAbs_EDGE || isStructured );
+            else
+              noShrinkV = ( ! i2S->second.IsSame( i2SAdj->second ));
+          }
+          else
+          {
+            noShrinkV = noShrinkE;
+          }
+        }
+        else
+        {
+          // the adjacent SOLID has NO layers at all
+          noShrinkV = ( isStructured ||
+                        _sdVec[i]._shrinkShape2Shape[ vID ].ShapeType() == TopAbs_EDGE );
+        }
+        if ( noShrinkV )
+          _sdVec[i]._noShrinkShapes.insert( vID );
+      }
+
+    } // loop on _sdVec[i]._shrinkShape2Shape
+  } // loop on _sdVec to fill in _SolidData::_noShrinkShapes
+
+
+    // add FACEs of other SOLIDs to _ignoreFaceIds
   for ( size_t i = 0; i < _sdVec.size(); ++i )
   {
     shapes.Clear();
@@ -2537,7 +2651,7 @@ bool _ViscousBuilder::makeLayer(_SolidData& data)
   if ( data._stepSize < 1. )
     data._epsilon *= data._stepSize;
 
-  if ( !findShapesToSmooth( data ))
+  if ( !findShapesToSmooth( data )) // _LayerEdge::_maxLen is computed here
     return false;
 
   // limit data._stepSize depending on surface curvature and fill data._convexFaces
@@ -2679,11 +2793,9 @@ void _ViscousBuilder::limitStepSize( _SolidData& data, const double minSize )
 
 void _ViscousBuilder::limitStepSizeByCurvature( _SolidData& data )
 {
-  const int nbTestPnt = 5; // on a FACE sub-shape
-
-  BRepLProp_SLProps surfProp( 2, 1e-6 );
   SMESH_MesherHelper helper( *_mesh );
 
+  BRepLProp_SLProps surfProp( 2, 1e-6 );
   data._convexFaces.clear();
 
   for ( size_t iS = 0; iS < data._edgesOnShape.size(); ++iS )
@@ -2694,57 +2806,62 @@ void _ViscousBuilder::limitStepSizeByCurvature( _SolidData& data )
       continue;
 
     TopoDS_Face        F = TopoDS::Face( eof._shape );
-    SMESH_subMesh *   sm = eof._subMesh;
     const TGeomID faceID = eof._shapeID;
 
     BRepAdaptor_Surface surface( F, false );
     surfProp.SetSurface( surface );
 
-    bool isTooCurved = false;
-
     _ConvexFace cnvFace;
-    const double        oriFactor = ( F.Orientation() == TopAbs_REVERSED ? +1. : -1. );
-    SMESH_subMeshIteratorPtr smIt = sm->getDependsOnIterator(/*includeSelf=*/true);
-    while ( smIt->more() )
+    cnvFace._face = F;
+    cnvFace._normalsFixed = false;
+    cnvFace._isTooCurved = false;
+
+    double maxCurvature = cnvFace.GetMaxCurvature( data, eof, surfProp, helper );
+    if ( maxCurvature > 0 )
     {
-      sm = smIt->next();
-      const TGeomID subID = sm->GetId();
-      // find _LayerEdge's of a sub-shape
-      _EdgesOnShape* eos;
-      if (( eos = data.GetShapeEdges( subID )))
-        cnvFace._subIdToEOS.insert( make_pair( subID, eos ));
-      else
-        continue;
-      // check concavity and curvature and limit data._stepSize
+      limitStepSize( data, 0.9 / maxCurvature );
+      findEdgesToUpdateNormalNearConvexFace( cnvFace, data, helper );
+    }
+    if ( !cnvFace._isTooCurved ) continue;
+
+    _ConvexFace & convFace =
+      data._convexFaces.insert( make_pair( faceID, cnvFace )).first->second;
+
+    // skip a closed surface (data._convexFaces is useful anyway)
+    bool isClosedF = false;
+    helper.SetSubShape( F );
+    if ( helper.HasRealSeam() )
+    {
+      // in the closed surface there must be a closed EDGE
+      for ( TopExp_Explorer eIt( F, TopAbs_EDGE ); eIt.More() && !isClosedF; eIt.Next() )
+        isClosedF = helper.IsClosedEdge( TopoDS::Edge( eIt.Current() ));
+    }
+    if ( isClosedF )
+    {
+      // limit _LayerEdge::_maxLen on the FACE
+      const double oriFactor    = ( F.Orientation() == TopAbs_REVERSED ? +1. : -1. );
       const double minCurvature =
-        1. / ( eos->_hyp.GetTotalThickness() * ( 1 + theThickToIntersection ));
-      size_t iStep = Max( 1, eos->_edges.size() / nbTestPnt );
-      for ( size_t i = 0; i < eos->_edges.size(); i += iStep )
+        1. / ( eof._hyp.GetTotalThickness() * ( 1 + theThickToIntersection ));
+      map< TGeomID, _EdgesOnShape* >::iterator id2eos = cnvFace._subIdToEOS.find( faceID );
+      if ( id2eos != cnvFace._subIdToEOS.end() )
       {
-        gp_XY uv = helper.GetNodeUV( F, eos->_edges[ i ]->_nodes[0] );
-        surfProp.SetParameters( uv.X(), uv.Y() );
-        if ( !surfProp.IsCurvatureDefined() )
-          continue;
-        if ( surfProp.MaxCurvature() * oriFactor > minCurvature )
-        {
-          limitStepSize( data, 0.9 / surfProp.MaxCurvature() * oriFactor );
-          isTooCurved = true;
-        }
-        if ( surfProp.MinCurvature() * oriFactor > minCurvature )
+        _EdgesOnShape& eos = * id2eos->second;
+        for ( size_t i = 0; i < eos._edges.size(); ++i )
         {
-          limitStepSize( data, 0.9 / surfProp.MinCurvature() * oriFactor );
-          isTooCurved = true;
+          _LayerEdge* ledge = eos._edges[ i ];
+          gp_XY uv = helper.GetNodeUV( F, ledge->_nodes[0] );
+          surfProp.SetParameters( uv.X(), uv.Y() );
+          if ( surfProp.IsCurvatureDefined() )
+          {
+            double curvature = Max( surfProp.MaxCurvature() * oriFactor,
+                                    surfProp.MinCurvature() * oriFactor );
+            if ( curvature > minCurvature )
+              ledge->_maxLen = Min( ledge->_maxLen, 1. / curvature );
+          }
         }
       }
-    } // loop on sub-shapes of the FACE
-
-    if ( !isTooCurved ) continue;
-
-    _ConvexFace & convFace =
-      data._convexFaces.insert( make_pair( faceID, cnvFace )).first->second;
-
-    convFace._face = F;
-    convFace._normalsFixed = false;
+      continue;
+    }
 
     // Fill _ConvexFace::_simplexTestEdges. These _LayerEdge's are used to detect
     // prism distortion.
@@ -2810,18 +2927,9 @@ bool _ViscousBuilder::findShapesToSmooth( _SolidData& data )
   // define allowed thickness
   computeGeomSize( data ); // compute data._geomSize and _LayerEdge::_maxLen
 
-  data._maxThickness = 0;
-  data._minThickness = 1e100;
-  list< const StdMeshers_ViscousLayers* >::iterator hyp = data._hyps.begin();
-  for ( ; hyp != data._hyps.end(); ++hyp )
-  {
-    data._maxThickness = Max( data._maxThickness, (*hyp)->GetTotalThickness() );
-    data._minThickness = Min( data._minThickness, (*hyp)->GetTotalThickness() );
-  }
-  //const double tgtThick = /*Min( 0.5 * data._geomSize, */data._maxThickness;
 
   // Find shapes needing smoothing; such a shape has _LayerEdge._normal on it's
-  // boundry inclined to the shape at a sharp angle
+  // boundary inclined to the shape at a sharp angle
 
   //list< TGeomID > shapesToSmooth;
   TopTools_MapOfShape edgesOfSmooFaces;
@@ -2890,7 +2998,7 @@ bool _ViscousBuilder::findShapesToSmooth( _SolidData& data )
     {
       TGeomID iV = getMeshDS()->ShapeToIndex( vIt.Value() );
       vector<_LayerEdge*>& eV = edgesByGeom[ iV ]._edges;
-      if ( eV.empty() ) continue;
+      if ( eV.empty() || eV[0]->Is( _LayerEdge::MULTI_NORMAL )) continue;
       gp_Vec  eDir    = getEdgeDir( E, TopoDS::Vertex( vIt.Value() ));
       double angle    = eDir.Angle( eV[0]->_normal );
       double cosin    = Cos( angle );
@@ -2917,8 +3025,8 @@ bool _ViscousBuilder::findShapesToSmooth( _SolidData& data )
         {
           eos._edgeSmoother = new _Smoother1D( curve, eos );
 
-          for ( size_t i = 0; i < eos._edges.size(); ++i )
-            eos._edges[i]->Set( _LayerEdge::TO_SMOOTH );
+          // for ( size_t i = 0; i < eos._edges.size(); ++i )
+          //   eos._edges[i]->Set( _LayerEdge::TO_SMOOTH );
         }
       }
     }
@@ -3297,9 +3405,14 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge&         edge,
   }
 
   // find _normal
+  bool fromVonF = false;
   if ( useGeometry )
   {
-    if ( onShrinkShape ) // one of faces the node is on has no layers
+    fromVonF = ( eos.ShapeType() == TopAbs_VERTEX &&
+                 eos.SWOLType()  == TopAbs_FACE  &&
+                 totalNbFaces > 1 );
+
+    if ( onShrinkShape && !fromVonF ) // one of faces the node is on has no layers
     {
       if ( eos.SWOLType() == TopAbs_EDGE )
       {
@@ -3319,12 +3432,15 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge&         edge,
                                    node, helper, normOK);
       }
     }
-    else // layers are on all FACEs of SOLID the node is on
+    else // layers are on all FACEs of SOLID the node is on (or fromVonF)
     {
+      if ( fromVonF )
+        face2Norm[ totalNbFaces++ ].first = TopoDS::Face( eos._sWOL );
+
       int nbOkNorms = 0;
-      for ( int iF = 0; iF < totalNbFaces; ++iF )
+      for ( int iF = totalNbFaces - 1; iF >= 0; --iF )
       {
-        F = TopoDS::Face( face2Norm[ iF ].first );
+        F = face2Norm[ iF ].first;
         geomNorm = getFaceNormal( node, F, helper, normOK );
         if ( !normOK ) continue;
         nbOkNorms++;
@@ -3337,11 +3453,16 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge&         edge,
       if ( nbOkNorms == 0 )
         return error(SMESH_Comment("Can't get normal to node ") << node->GetID(), data._index);
 
+      if ( totalNbFaces >= 3 )
+      {
+        edge._normal = getNormalByOffset( &edge, face2Norm, totalNbFaces, fromVonF );
+      }
+
       if ( edge._normal.Modulus() < 1e-3 && nbOkNorms > 1 )
       {
         // opposite normals, re-get normals at shifted positions (IPAL 52426)
         edge._normal.SetCoord( 0,0,0 );
-        for ( int iF = 0; iF < totalNbFaces; ++iF )
+        for ( int iF = 0; iF < totalNbFaces - fromVonF; ++iF )
         {
           const TopoDS_Face& F = face2Norm[iF].first;
           geomNorm = getFaceNormal( node, F, helper, normOK, /*shiftInside=*/true );
@@ -3352,11 +3473,6 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge&         edge,
           edge._normal += face2Norm[ iF ].second;
         }
       }
-
-      if ( totalNbFaces >= 3 )
-      {
-        edge._normal = getNormalByOffset( &edge, face2Norm, totalNbFaces );
-      }
     }
   }
   else // !useGeometry - get _normal using surrounding mesh faces
@@ -3396,20 +3512,26 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge&         edge,
       break;
     }
     case TopAbs_VERTEX: {
-      if ( eos.SWOLType() != TopAbs_FACE ) { // else _cosin is set by getFaceDir()
+      if ( fromVonF )
+      {
+        getFaceDir( TopoDS::Face( eos._sWOL ), TopoDS::Vertex( eos._shape ),
+                    node, helper, normOK, &edge._cosin );
+      }
+      else if ( eos.SWOLType() != TopAbs_FACE ) // else _cosin is set by getFaceDir()
+      {
         TopoDS_Vertex V  = TopoDS::Vertex( eos._shape );
         gp_Vec inFaceDir = getFaceDir( F, V, node, helper, normOK );
         double angle     = inFaceDir.Angle( edge._normal ); // [0,PI]
         edge._cosin      = Cos( angle );
         if ( totalNbFaces > 2 || helper.IsSeamShape( node->getshapeId() ))
-          for ( int iF = totalNbFaces-2; iF >=0; --iF )
+          for ( int iF = 1; iF < totalNbFaces; ++iF )
           {
             F = face2Norm[ iF ].first;
             inFaceDir = getFaceDir( F, V, node, helper, normOK=true );
             if ( normOK ) {
               double angle = inFaceDir.Angle( edge._normal );
               double cosin = Cos( angle );
-              if ( Abs( cosin ) > edge._cosin )
+              if ( Abs( cosin ) > Abs( edge._cosin ))
                 edge._cosin = cosin;
             }
           }
@@ -3732,7 +3854,8 @@ gp_XYZ _ViscousBuilder::getWeigthedNormal( const _LayerEdge* edge )
 
 gp_XYZ _ViscousBuilder::getNormalByOffset( _LayerEdge*                      edge,
                                            std::pair< TopoDS_Face, gp_XYZ > f2Normal[],
-                                           int                              nbFaces )
+                                           int                              nbFaces,
+                                           bool                             lastNoOffset)
 {
   SMESH_TNodeXYZ p0 = edge->_nodes[0];
 
@@ -3747,11 +3870,16 @@ gp_XYZ _ViscousBuilder::getNormalByOffset( _LayerEdge*                      edge
 
   // prepare _OffsetPlane's
   vector< _OffsetPlane > pln( nbFaces );
-  for ( int i = 0; i < nbFaces; ++i )
+  for ( int i = 0; i < nbFaces - lastNoOffset; ++i )
   {
     pln[i]._faceIndex = i;
     pln[i]._plane = gp_Pln( p0 + f2Normal[i].second, f2Normal[i].second );
   }
+  if ( lastNoOffset )
+  {
+    pln[ nbFaces - 1 ]._faceIndex = nbFaces - 1;
+    pln[ nbFaces - 1 ]._plane = gp_Pln( p0, f2Normal[ nbFaces - 1 ].second );
+  }
 
   // intersect neighboring OffsetPlane's
   PShapeIteratorPtr edgeIt = SMESH_MesherHelper::GetAncestors( V, *_mesh, TopAbs_EDGE );
@@ -3763,7 +3891,7 @@ gp_XYZ _ViscousBuilder::getNormalByOffset( _LayerEdge*                      edge
         (( f1 < 0 ) ? f1 : f2 ) = i;
 
     if ( f2 >= 0 )
-      pln[ f1 ].ComputeIntersectionLine( pln[ f2 ]);
+      pln[ f1 ].ComputeIntersectionLine( pln[ f2 ], TopoDS::Edge( *edge ), TopoDS::Vertex( V ));
   }
 
   // get a common point
@@ -3772,7 +3900,7 @@ gp_XYZ _ViscousBuilder::getNormalByOffset( _LayerEdge*                      edge
   bool isPointFound;
   for ( int i = 0; i < nbFaces; ++i )
   {
-    commonPnt += pln[ i ].GetCommonPoint( isPointFound );
+    commonPnt += pln[ i ].GetCommonPoint( isPointFound, TopoDS::Vertex( V ));
     nbPoints  += isPointFound;
   }
   gp_XYZ wgtNorm = getWeigthedNormal( edge );
@@ -3781,13 +3909,15 @@ gp_XYZ _ViscousBuilder::getNormalByOffset( _LayerEdge*                      edge
 
   commonPnt /= nbPoints;
   resNorm = commonPnt - p0;
+  if ( lastNoOffset )
+    return resNorm;
 
   // choose the best among resNorm and wgtNorm
   resNorm.Normalize();
   wgtNorm.Normalize();
   double resMinDot = std::numeric_limits<double>::max();
   double wgtMinDot = std::numeric_limits<double>::max();
-  for ( int i = 0; i < nbFaces; ++i )
+  for ( int i = 0; i < nbFaces - lastNoOffset; ++i )
   {
     resMinDot = Min( resMinDot, resNorm * f2Normal[i].second );
     wgtMinDot = Min( wgtMinDot, wgtNorm * f2Normal[i].second );
@@ -3807,7 +3937,9 @@ gp_XYZ _ViscousBuilder::getNormalByOffset( _LayerEdge*                      edge
  */
 //================================================================================
 
-void _OffsetPlane::ComputeIntersectionLine( _OffsetPlane& pln )
+void _OffsetPlane::ComputeIntersectionLine( _OffsetPlane&        pln,
+                                            const TopoDS_Edge&   E,
+                                            const TopoDS_Vertex& V )
 {
   int iNext = bool( _faceIndexNext[0] >= 0 );
   _faceIndexNext[ iNext ] = pln._faceIndex;
@@ -3831,31 +3963,37 @@ void _OffsetPlane::ComputeIntersectionLine( _OffsetPlane& pln )
     else       cooMax = 3;
   }
 
-  if ( Abs( lineDir.Coord( cooMax )) < 0.05 )
-    return;
-
   gp_Pnt linePos;
-  // the constants in the 2 plane equations
-  double d1 = - ( _plane.Axis().Direction().XYZ()     * _plane.Location().XYZ() );
-  double d2 = - ( pln._plane.Axis().Direction().XYZ() * pln._plane.Location().XYZ() );
-
-  switch ( cooMax ) {
-  case 1:
-    linePos.SetX(  0 );
-    linePos.SetY(( d2*n1.Z() - d1*n2.Z()) / lineDir.X() );
-    linePos.SetZ(( d1*n2.Y() - d2*n1.Y()) / lineDir.X() );
-    break;
-  case 2:
-    linePos.SetX(( d1*n2.Z() - d2*n1.Z()) / lineDir.Y() );
-    linePos.SetY(  0 );
-    linePos.SetZ(( d2*n1.X() - d1*n2.X()) / lineDir.Y() );
-    break;
-  case 3:
-    linePos.SetX(( d2*n1.Y() - d1*n2.Y()) / lineDir.Z() );
-    linePos.SetY(( d1*n2.X() - d2*n1.X()) / lineDir.Z() );
-    linePos.SetZ(  0 );
+  if ( Abs( lineDir.Coord( cooMax )) < 0.05 )
+  {
+    // parallel planes - intersection is an offset of the common EDGE
+    gp_Pnt p = BRep_Tool::Pnt( V );
+    linePos  = 0.5 * (( p.XYZ() + n1 ) + ( p.XYZ() + n2 ));
+    lineDir  = getEdgeDir( E, V );
   }
+  else
+  {
+    // the constants in the 2 plane equations
+    double d1 = - ( _plane.Axis().Direction().XYZ()     * _plane.Location().XYZ() );
+    double d2 = - ( pln._plane.Axis().Direction().XYZ() * pln._plane.Location().XYZ() );
 
+    switch ( cooMax ) {
+    case 1:
+      linePos.SetX(  0 );
+      linePos.SetY(( d2*n1.Z() - d1*n2.Z()) / lineDir.X() );
+      linePos.SetZ(( d1*n2.Y() - d2*n1.Y()) / lineDir.X() );
+      break;
+    case 2:
+      linePos.SetX(( d1*n2.Z() - d2*n1.Z()) / lineDir.Y() );
+      linePos.SetY(  0 );
+      linePos.SetZ(( d2*n1.X() - d1*n2.X()) / lineDir.Y() );
+      break;
+    case 3:
+      linePos.SetX(( d2*n1.Y() - d1*n2.Y()) / lineDir.Z() );
+      linePos.SetY(( d1*n2.X() - d2*n1.X()) / lineDir.Z() );
+      linePos.SetZ(  0 );
+    }
+  }
   gp_Lin& line = _lines[ iNext ];
   line.SetDirection( lineDir );
   line.SetLocation ( linePos );
@@ -3875,7 +4013,8 @@ void _OffsetPlane::ComputeIntersectionLine( _OffsetPlane& pln )
  */
 //================================================================================
 
-gp_XYZ _OffsetPlane::GetCommonPoint(bool& isFound) const
+gp_XYZ _OffsetPlane::GetCommonPoint(bool&                 isFound,
+                                    const TopoDS_Vertex & V) const
 {
   gp_XYZ p( 0,0,0 );
   isFound = false;
@@ -3884,13 +4023,22 @@ gp_XYZ _OffsetPlane::GetCommonPoint(bool& isFound) const
   {
     gp_Vec lPerp0 = _lines[0].Direction().XYZ() ^ _plane.Axis().Direction().XYZ();
     double  dot01 = lPerp0 * _lines[1].Direction().XYZ();
-    if ( Abs( dot01 ) > std::numeric_limits<double>::min() )
+    if ( Abs( dot01 ) > 0.05 )
     {
       gp_Vec l0l1 = _lines[1].Location().XYZ() - _lines[0].Location().XYZ();
       double   u1 = - ( lPerp0 * l0l1 ) / dot01;
       p = ( _lines[1].Location().XYZ() + _lines[1].Direction().XYZ() * u1 );
       isFound = true;
     }
+    else
+    {
+      gp_Pnt  pV ( BRep_Tool::Pnt( V ));
+      gp_Vec  lv0( _lines[0].Location(), pV    ),  lv1(_lines[1].Location(), pV     );
+      double dot0( lv0 * _lines[0].Direction() ), dot1( lv1 * _lines[1].Direction() );
+      p += 0.5 * ( _lines[0].Location().XYZ() + _lines[0].Direction().XYZ() * dot0 );
+      p += 0.5 * ( _lines[1].Location().XYZ() + _lines[1].Direction().XYZ() * dot1 );
+      isFound = true;
+    }
   }
 
   return p;
@@ -3996,7 +4144,7 @@ void _LayerEdge::SetDataByNeighbors( const SMDS_MeshNode* n1,
 //================================================================================
 /*!
  * \brief Copy data from a _LayerEdge of other SOLID and based on the same node;
- * this and other _LayerEdge's are inflated along a FACE or an EDGE
+ * this and the other _LayerEdge are inflated along a FACE or an EDGE
  */
 //================================================================================
 
@@ -4207,7 +4355,7 @@ void _ViscousBuilder::computeGeomSize( _SolidData& data )
     _EdgesOnShape& eos = data._edgesOnShape[ iS ];
     if ( eos._edges.empty() )
       continue;
-    // get neighbor faces intersection with which should not be considered since
+    // get neighbor faces, intersection with which should not be considered since
     // collisions are avoided by means of smoothing
     set< TGeomID > neighborFaces;
     if ( eos._hyp.ToSmooth() )
@@ -4237,6 +4385,78 @@ void _ViscousBuilder::computeGeomSize( _SolidData& data )
       }
     }
   }
+
+  data._maxThickness = 0;
+  data._minThickness = 1e100;
+  list< const StdMeshers_ViscousLayers* >::iterator hyp = data._hyps.begin();
+  for ( ; hyp != data._hyps.end(); ++hyp )
+  {
+    data._maxThickness = Max( data._maxThickness, (*hyp)->GetTotalThickness() );
+    data._minThickness = Min( data._minThickness, (*hyp)->GetTotalThickness() );
+  }
+
+  // Limit inflation step size by geometry size found by intersecting
+  // normals of _LayerEdge's with mesh faces
+  if ( data._stepSize > 0.3 * data._geomSize )
+    limitStepSize( data, 0.3 * data._geomSize );
+
+  if ( data._stepSize > data._minThickness )
+    limitStepSize( data, data._minThickness );
+
+
+  // -------------------------------------------------------------------------
+  // Detect _LayerEdge which can't intersect with opposite or neighbor layer,
+  // so no need in detecting intersection at each inflation step
+  // -------------------------------------------------------------------------
+
+  int nbSteps = data._maxThickness / data._stepSize;
+  if ( nbSteps < 3 || nbSteps * data._n2eMap.size() < 100000 )
+    return;
+
+  vector< const SMDS_MeshElement* > closeFaces;
+  int nbDetected = 0;
+
+  for ( size_t iS = 0; iS < data._edgesOnShape.size(); ++iS )
+  {
+    _EdgesOnShape& eos = data._edgesOnShape[ iS ];
+    if ( eos._edges.empty() || eos.ShapeType() != TopAbs_FACE )
+      continue;
+
+    for ( size_t i = 0; i < eos.size(); ++i )
+    {
+      SMESH_NodeXYZ p( eos[i]->_nodes[0] );
+      double radius = data._maxThickness + 2 * eos[i]->_maxLen;
+      closeFaces.clear();
+      searcher->GetElementsInSphere( p, radius, SMDSAbs_Face, closeFaces );
+
+      bool toIgnore = true;
+      for ( size_t iF = 0; iF < closeFaces.size()  && toIgnore; ++iF )
+        if ( !( toIgnore = ( closeFaces[ iF ]->getshapeId() == eos._shapeID ||
+                             data._ignoreFaceIds.count( closeFaces[ iF ]->getshapeId() ))))
+        {
+          // check if a _LayerEdge will inflate in a direction opposite to a direction
+          // toward a close face
+          bool allBehind = true;
+          for ( int iN = 0; iN < closeFaces[ iF ]->NbCornerNodes()  && allBehind; ++iN )
+          {
+            SMESH_NodeXYZ pi( closeFaces[ iF ]->GetNode( iN ));
+            allBehind = (( pi - p ) * eos[i]->_normal < 0.1 * data._stepSize );
+          }
+          toIgnore = allBehind;
+        }
+
+
+      if ( toIgnore ) // no need to detect intersection
+      {
+        eos[i]->Set( _LayerEdge::INTERSECTED );
+        ++nbDetected;
+      }
+    }
+  }
+
+  debugMsg( "Nb LE to intersect " << data._n2eMap.size()-nbDetected << ", ignore " << nbDetected );
+
+  return;
 }
 
 //================================================================================
@@ -4249,22 +4469,20 @@ bool _ViscousBuilder::inflate(_SolidData& data)
 {
   SMESH_MesherHelper helper( *_mesh );
 
-  // Limit inflation step size by geometry size found by itersecting
-  // normals of _LayerEdge's with mesh faces
-  if ( data._stepSize > 0.3 * data._geomSize )
-    limitStepSize( data, 0.3 * data._geomSize );
-
   const double tgtThick = data._maxThickness;
-  if ( data._stepSize > data._minThickness )
-    limitStepSize( data, data._minThickness );
 
   if ( data._stepSize < 1. )
     data._epsilon = data._stepSize * 1e-7;
 
   debugMsg( "-- geomSize = " << data._geomSize << ", stepSize = " << data._stepSize );
+  _pyDump->Pause();
 
   findCollisionEdges( data, helper );
 
+  limitMaxLenByCurvature( data, helper );
+
+  _pyDump->Resume();
+
   // limit length of _LayerEdge's around MULTI_NORMAL _LayerEdge's
   for ( size_t i = 0; i < data._edgesOnShape.size(); ++i )
     if ( data._edgesOnShape[i].ShapeType() == TopAbs_VERTEX &&
@@ -4344,7 +4562,10 @@ bool _ViscousBuilder::inflate(_SolidData& data)
       const double shapeTgtThick = eos._hyp.GetTotalThickness();
       for ( size_t i = 0; i < eos._edges.size(); ++i )
       {
-        avgThick      += Min( 1., eos._edges[i]->_len / shapeTgtThick );
+        if ( eos._edges[i]->_nodes.size() > 1 )
+          avgThick    += Min( 1., eos._edges[i]->_len / shapeTgtThick );
+        else
+          avgThick    += shapeTgtThick;
         nbActiveEdges += ( ! eos._edges[i]->Is( _LayerEdge::BLOCKED ));
       }
     }
@@ -4366,6 +4587,7 @@ bool _ViscousBuilder::inflate(_SolidData& data)
       break;
     }
 #endif
+
     // new step size
     limitStepSize( data, 0.25 * distToIntersection );
     if ( data._stepSizeNodes[0] )
@@ -4448,20 +4670,30 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data,
       {
         // smooth disabled by the user; check validy only
         if ( !isFace ) continue;
+        badEdges.clear();
         for ( size_t i = 0; i < eos._edges.size(); ++i )
         {
           _LayerEdge* edge = eos._edges[i];
           for ( size_t iF = 0; iF < edge->_simplices.size(); ++iF )
             if ( !edge->_simplices[iF].IsForward( edge->_nodes[0], edge->_pos.back(), vol ))
             {
-              debugMsg( "-- Stop inflation. Bad simplex ("
-                        << " "<< edge->_nodes[0]->GetID()
-                        << " "<< edge->_nodes.back()->GetID()
-                        << " "<< edge->_simplices[iF]._nPrev->GetID()
-                        << " "<< edge->_simplices[iF]._nNext->GetID() << " ) ");
-              return false;
+              // debugMsg( "-- Stop inflation. Bad simplex ("
+              //           << " "<< edge->_nodes[0]->GetID()
+              //           << " "<< edge->_nodes.back()->GetID()
+              //           << " "<< edge->_simplices[iF]._nPrev->GetID()
+              //           << " "<< edge->_simplices[iF]._nNext->GetID() << " ) ");
+              // return false;
+              badEdges.push_back( edge );
             }
         }
+        if ( !badEdges.empty() )
+        {
+          eosC1.resize(1);
+          eosC1[0] = &eos;
+          int nbBad = invalidateBadSmooth( data, helper, badEdges, eosC1, infStep );
+          if ( nbBad > 0 )
+            return false;
+        }
         continue; // goto the next EDGE or FACE
       }
 
@@ -4608,7 +4840,7 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data,
           if (( step % 3 == 1 ) || ( nbBad > 0 && step >= stepLimit / 2 ))
             for ( size_t iEOS = 0; iEOS < eosC1.size(); ++iEOS )
             {
-              putOnOffsetSurface( *eosC1[ iEOS ], infStep, step, /*moveAll=*/step == 1 );
+              putOnOffsetSurface( *eosC1[ iEOS ], infStep, eosC1, step, /*moveAll=*/step == 1 );
             }
 
         } // smoothing steps
@@ -4617,10 +4849,10 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data,
         for ( size_t iEOS = 0; iEOS < eosC1.size(); ++iEOS )
         {
           if ( ! eosC1[ iEOS ]->_eosConcaVer.empty() || nbBad > 0 )
-            putOnOffsetSurface( *eosC1[ iEOS ], infStep );
+            putOnOffsetSurface( *eosC1[ iEOS ], infStep, eosC1 );
         }
 
-        if ( !badEdges.empty() )
+        //if ( !badEdges.empty() )
         {
           badEdges.clear();
           for ( size_t iEOS = 0; iEOS < eosC1.size(); ++iEOS )
@@ -4631,10 +4863,11 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data,
 
               _LayerEdge* edge = eosC1[ iEOS ]->_edges[i];
               edge->CheckNeiborsOnBoundary( & badEdges );
-              if ( nbBad > 0 )
+              if (( nbBad > 0 ) ||
+                  ( edge->Is( _LayerEdge::BLOCKED ) && edge->Is( _LayerEdge::NEAR_BOUNDARY )))
               {
                 SMESH_TNodeXYZ tgtXYZ = edge->_nodes.back();
-                const gp_XYZ& prevXYZ = edge->PrevCheckPos();
+                gp_XYZ        prevXYZ = edge->PrevCheckPos();
                 for ( size_t j = 0; j < edge->_simplices.size(); ++j )
                   if ( !edge->_simplices[j].IsForward( &prevXYZ, &tgtXYZ, vol ))
                   {
@@ -4676,7 +4909,8 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data,
       _LayerEdge*      edge = eos._edges[i];
       if ( edge->_nodes.size() < 2 ) continue;
       SMESH_TNodeXYZ tgtXYZ = edge->_nodes.back();
-      const gp_XYZ& prevXYZ = edge->PrevCheckPos();
+      //SMESH_TNodeXYZ prevXYZ = edge->_nodes[0];
+      gp_XYZ        prevXYZ = edge->PrevCheckPos( &eos );
       //const gp_XYZ& prevXYZ = edge->PrevPos();
       for ( size_t j = 0; j < edge->_simplices.size(); ++j )
         if ( !edge->_simplices[j].IsForward( &prevXYZ, &tgtXYZ, vol ))
@@ -4727,39 +4961,42 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data,
         continue;
       if ( eos._edges[i]->FindIntersection( *searcher, dist, data._epsilon, eos, &intFace ))
       {
+        return false;
+        // commented due to "Illegal hash-positionPosition" error in NETGEN
+        // on Debian60 on viscous_layers_01/B2 case
         // Collision; try to deflate _LayerEdge's causing it
-        badEdges.clear();
-        badEdges.push_back( eos._edges[i] );
-        eosC1[0] = & eos;
-        int nbBad = invalidateBadSmooth( data, helper, badEdges, eosC1, infStep );
-        if ( nbBad > 0 )
-          return false;
+        // badEdges.clear();
+        // badEdges.push_back( eos._edges[i] );
+        // eosC1[0] = & eos;
+        // int nbBad = invalidateBadSmooth( data, helper, badEdges, eosC1, infStep );
+        // if ( nbBad > 0 )
+        //   return false;
 
-        badEdges.clear();
-        if ( _EdgesOnShape* eof = data.GetShapeEdges( intFace->getshapeId() ))
-        {
-          if ( const _TmpMeshFace* f = dynamic_cast< const _TmpMeshFace*>( intFace ))
-          {
-            const SMDS_MeshElement* srcFace =
-              eof->_subMesh->GetSubMeshDS()->GetElement( f->getIdInShape() );
-            SMDS_ElemIteratorPtr nIt = srcFace->nodesIterator();
-            while ( nIt->more() )
-            {
-              const SMDS_MeshNode* srcNode = static_cast<const SMDS_MeshNode*>( nIt->next() );
-              TNode2Edge::iterator n2e = data._n2eMap.find( srcNode );
-              if ( n2e != data._n2eMap.end() )
-                badEdges.push_back( n2e->second );
-            }
-            eosC1[0] = eof;
-            nbBad = invalidateBadSmooth( data, helper, badEdges, eosC1, infStep );
-            if ( nbBad > 0 )
-              return false;
-          }
-        }
-        if ( eos._edges[i]->FindIntersection( *searcher, dist, data._epsilon, eos, &intFace ))
-          return false;
-        else
-          continue;
+        // badEdges.clear();
+        // if ( _EdgesOnShape* eof = data.GetShapeEdges( intFace->getshapeId() ))
+        // {
+        //   if ( const _TmpMeshFace* f = dynamic_cast< const _TmpMeshFace*>( intFace ))
+        //   {
+        //     const SMDS_MeshElement* srcFace =
+        //       eof->_subMesh->GetSubMeshDS()->GetElement( f->getIdInShape() );
+        //     SMDS_ElemIteratorPtr nIt = srcFace->nodesIterator();
+        //     while ( nIt->more() )
+        //     {
+        //       const SMDS_MeshNode* srcNode = static_cast<const SMDS_MeshNode*>( nIt->next() );
+        //       TNode2Edge::iterator n2e = data._n2eMap.find( srcNode );
+        //       if ( n2e != data._n2eMap.end() )
+        //         badEdges.push_back( n2e->second );
+        //     }
+        //     eosC1[0] = eof;
+        //     nbBad = invalidateBadSmooth( data, helper, badEdges, eosC1, infStep );
+        //     if ( nbBad > 0 )
+        //       return false;
+        //   }
+        // }
+        // if ( eos._edges[i]->FindIntersection( *searcher, dist, data._epsilon, eos, &intFace ))
+        //   return false;
+        // else
+        //   continue;
       }
       if ( !intFace )
       {
@@ -4774,7 +5011,7 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data,
         // ignore intersection of a _LayerEdge based on a _ConvexFace with a face
         // lying on this _ConvexFace
         if ( _ConvexFace* convFace = data.GetConvexFace( intFace->getshapeId() ))
-          if ( convFace->_subIdToEOS.count ( eos._shapeID ))
+          if ( convFace->_isTooCurved && convFace->_subIdToEOS.count ( eos._shapeID ))
             continue;
 
         // ignore intersection of a _LayerEdge based on a FACE with an element on this FACE
@@ -4786,15 +5023,16 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data,
         if ( dist > 0 )
         {
           bool toIgnore = false;
-          if (  eos._edges[i]->Is( _LayerEdge::TO_SMOOTH ))
+          if (  eos._toSmooth )
           {
             const TopoDS_Shape& S = getMeshDS()->IndexToShape( intFace->getshapeId() );
             if ( !S.IsNull() && S.ShapeType() == TopAbs_FACE )
             {
-              TopExp_Explorer edge( eos._shape, TopAbs_EDGE );
-              for ( ; !toIgnore && edge.More(); edge.Next() )
-                // is adjacent - has a common EDGE
-                toIgnore = ( helper.IsSubShape( edge.Current(), S ));
+              TopExp_Explorer sub( eos._shape,
+                                   eos.ShapeType() == TopAbs_FACE ? TopAbs_EDGE : TopAbs_VERTEX );
+              for ( ; !toIgnore && sub.More(); sub.Next() )
+                // is adjacent - has a common EDGE or VERTEX
+                toIgnore = ( helper.IsSubShape( sub.Current(), S ));
 
               if ( toIgnore ) // check angle between normals
               {
@@ -4854,16 +5092,16 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data,
     } // loop on eos._edges
   } // loop on data._edgesOnShape
 
-#ifdef __myDEBUG
-  if ( closestFace )
+  if ( closestFace && le )
   {
+#ifdef __myDEBUG
     SMDS_MeshElement::iterator nIt = closestFace->begin_nodes();
     cout << "Shortest distance: _LayerEdge nodes: tgt " << le->_nodes.back()->GetID()
          << " src " << le->_nodes[0]->GetID()<< ", intersection with face ("
          << (*nIt++)->GetID()<<" "<< (*nIt++)->GetID()<<" "<< (*nIt++)->GetID()
          << ") distance = " << distToIntersection<< endl;
-  }
 #endif
+  }
 
   return true;
 }
@@ -4886,53 +5124,119 @@ int _ViscousBuilder::invalidateBadSmooth( _SolidData&               data,
 
   dumpFunction(SMESH_Comment("invalidateBadSmooth")<<"_S"<<eosC1[0]->_shapeID<<"_InfStep"<<infStep);
 
-  //data.UnmarkEdges();
+  enum {
+    INVALIDATED   = _LayerEdge::UNUSED_FLAG,
+    TO_INVALIDATE = _LayerEdge::UNUSED_FLAG * 2,
+    ADDED         = _LayerEdge::UNUSED_FLAG * 4
+  };
+  data.UnmarkEdges( TO_INVALIDATE & INVALIDATED & ADDED );
 
   double vol;
-  //size_t iniNbBad = badSmooEdges.size();
+  bool haveInvalidated = true;
+  while ( haveInvalidated )
+  {
+    haveInvalidated = false;
+    for ( size_t i = 0; i < badSmooEdges.size(); ++i )
+    {
+      _LayerEdge*   edge = badSmooEdges[i];
+      _EdgesOnShape* eos = data.GetShapeEdges( edge );
+      edge->Set( ADDED );
+      bool invalidated = false;
+      if ( edge->Is( TO_INVALIDATE ) && edge->NbSteps() > 1 )
+      {
+        edge->InvalidateStep( edge->NbSteps(), *eos, /*restoreLength=*/true );
+        edge->Block( data );
+        edge->Set( INVALIDATED );
+        edge->Unset( TO_INVALIDATE );
+        invalidated = true;
+        haveInvalidated = true;
+      }
+
+      // look for _LayerEdge's of bad _simplices
+      int nbBad = 0;
+      SMESH_TNodeXYZ tgtXYZ  = edge->_nodes.back();
+      gp_XYZ        prevXYZ1 = edge->PrevCheckPos( eos );
+      //const gp_XYZ& prevXYZ2 = edge->PrevPos();
+      for ( size_t j = 0; j < edge->_simplices.size(); ++j )
+      {
+        if (( edge->_simplices[j].IsForward( &prevXYZ1, &tgtXYZ, vol ))/* &&
+            ( &prevXYZ1 == &prevXYZ2 || edge->_simplices[j].IsForward( &prevXYZ2, &tgtXYZ, vol ))*/)
+          continue;
+
+        bool isBad = true;
+        _LayerEdge* ee[2] = { 0,0 };
+        for ( size_t iN = 0; iN < edge->_neibors.size() &&   !ee[1]  ; ++iN )
+          if ( edge->_simplices[j].Includes( edge->_neibors[iN]->_nodes.back() ))
+            ee[ ee[0] != 0 ] = edge->_neibors[iN];
+
+        int maxNbSteps = Max( ee[0]->NbSteps(), ee[1]->NbSteps() );
+        while ( maxNbSteps > edge->NbSteps() && isBad )
+        {
+          --maxNbSteps;
+          for ( int iE = 0; iE < 2; ++iE )
+          {
+            if ( ee[ iE ]->NbSteps() > maxNbSteps &&
+                 ee[ iE ]->NbSteps() > 1 )
+            {
+              _EdgesOnShape* eos = data.GetShapeEdges( ee[ iE ] );
+              ee[ iE ]->InvalidateStep( ee[ iE ]->NbSteps(), *eos, /*restoreLength=*/true );
+              ee[ iE ]->Block( data );
+              ee[ iE ]->Set( INVALIDATED );
+              haveInvalidated = true;
+            }
+          }
+          if (( edge->_simplices[j].IsForward( &prevXYZ1, &tgtXYZ, vol )) /*&&
+              ( &prevXYZ1 == &prevXYZ2 || edge->_simplices[j].IsForward( &prevXYZ2, &tgtXYZ, vol ))*/)
+            isBad = false;
+        }
+        nbBad += isBad;
+        if ( !ee[0]->Is( ADDED )) badSmooEdges.push_back( ee[0] );
+        if ( !ee[1]->Is( ADDED )) badSmooEdges.push_back( ee[1] );
+        ee[0]->Set( ADDED );
+        ee[1]->Set( ADDED );
+        if ( isBad )
+        {
+          ee[0]->Set( TO_INVALIDATE );
+          ee[1]->Set( TO_INVALIDATE );
+        }
+      }
+
+      if ( !invalidated &&  nbBad > 0  &&  edge->NbSteps() > 1 )
+      {
+        _EdgesOnShape* eos = data.GetShapeEdges( edge );
+        edge->InvalidateStep( edge->NbSteps(), *eos, /*restoreLength=*/true );
+        edge->Block( data );
+        edge->Set( INVALIDATED );
+        edge->Unset( TO_INVALIDATE );
+        haveInvalidated = true;
+      }
+    } // loop on badSmooEdges
+  } // while ( haveInvalidated )
+
+  // re-smooth on analytical EDGEs
   for ( size_t i = 0; i < badSmooEdges.size(); ++i )
   {
     _LayerEdge* edge = badSmooEdges[i];
-    if ( edge->NbSteps() < 2 /*|| edge->Is( _LayerEdge::MARKED )*/)
-      continue;
+    if ( !edge->Is( INVALIDATED )) continue;
 
     _EdgesOnShape* eos = data.GetShapeEdges( edge );
-    edge->InvalidateStep( edge->NbSteps(), *eos, /*restoreLength=*/true );
-    edge->Block( data );
-    //edge->Set( _LayerEdge::MARKED );
-
-    // look for _LayerEdge's of bad _simplices
-    SMESH_TNodeXYZ tgtXYZ  = edge->_nodes.back();
-    const gp_XYZ& prevXYZ1 = edge->PrevCheckPos();
-    const gp_XYZ& prevXYZ2 = edge->PrevPos();
-    for ( size_t j = 0; j < edge->_simplices.size(); ++j )
-    {
-      if (( edge->_simplices[j].IsForward( &prevXYZ1, &tgtXYZ, vol )) &&
-          ( &prevXYZ1 == &prevXYZ2 || edge->_simplices[j].IsForward( &prevXYZ2, &tgtXYZ, vol )))
-        continue;
-      for ( size_t iN = 0; iN < edge->_neibors.size(); ++iN )
-        if ( edge->_simplices[j].Includes( edge->_neibors[iN]->_nodes.back() ))
-          badSmooEdges.push_back( edge->_neibors[iN] );
-    }
-
     if ( eos->ShapeType() == TopAbs_VERTEX )
     {
-      // re-smooth on analytical EDGEs
       PShapeIteratorPtr eIt = helper.GetAncestors( eos->_shape, *_mesh, TopAbs_EDGE );
       while ( const TopoDS_Shape* e = eIt->next() )
         if ( _EdgesOnShape* eoe = data.GetShapeEdges( *e ))
           if ( eoe->_edgeSmoother && eoe->_edgeSmoother->isAnalytic() )
           {
-            TopoDS_Face F; Handle(ShapeAnalysis_Surface) surface;
-            if ( eoe->SWOLType() == TopAbs_FACE ) {
-              F       = TopoDS::Face( eoe->_sWOL );
-              surface = helper.GetSurface( F );
-            }
-            eoe->_edgeSmoother->Perform( data, surface, F, helper );
+            // TopoDS_Face F; Handle(ShapeAnalysis_Surface) surface;
+            // if ( eoe->SWOLType() == TopAbs_FACE ) {
+            //   F       = TopoDS::Face( eoe->_sWOL );
+            //   surface = helper.GetSurface( F );
+            // }
+            // eoe->_edgeSmoother->Perform( data, surface, F, helper );
+            eoe->_edgeSmoother->_anaCurve.Nullify();
           }
-
     }
-  } // loop on badSmooEdges
+  }
 
 
   // check result of invalidation
@@ -4945,7 +5249,7 @@ int _ViscousBuilder::invalidateBadSmooth( _SolidData&               data,
       if ( !eosC1[ iEOS ]->_sWOL.IsNull() ) continue;
       _LayerEdge*      edge = eosC1[ iEOS ]->_edges[i];
       SMESH_TNodeXYZ tgtXYZ = edge->_nodes.back();
-      const gp_XYZ& prevXYZ = edge->PrevCheckPos();
+      gp_XYZ        prevXYZ = edge->PrevCheckPos( eosC1[ iEOS ]);
       for ( size_t j = 0; j < edge->_simplices.size(); ++j )
         if ( !edge->_simplices[j].IsForward( &prevXYZ, &tgtXYZ, vol ))
         {
@@ -4979,7 +5283,7 @@ void _ViscousBuilder::makeOffsetSurface( _EdgesOnShape& eos, SMESH_MesherHelper&
 
   // find offset
   gp_Pnt   tgtP = SMESH_TNodeXYZ( eos._edgeForOffset->_nodes.back() );
-  gp_Pnt2d   uv = baseSurface->ValueOfUV( tgtP, Precision::Confusion() );
+  /*gp_Pnt2d uv=*/baseSurface->ValueOfUV( tgtP, Precision::Confusion() );
   double offset = baseSurface->Gap();
 
   eos._offsetSurf.Nullify();
@@ -5009,26 +5313,48 @@ void _ViscousBuilder::makeOffsetSurface( _EdgesOnShape& eos, SMESH_MesherHelper&
  */
 //================================================================================
 
-void _ViscousBuilder::putOnOffsetSurface( _EdgesOnShape& eos,
-                                          int            infStep,
-                                          int            smooStep,
-                                          bool           moveAll )
+void _ViscousBuilder::putOnOffsetSurface( _EdgesOnShape&            eos,
+                                          int                       infStep,
+                                          vector< _EdgesOnShape* >& eosC1,
+                                          int                       smooStep,
+                                          int                       moveAll )
 {
-  if ( eos._offsetSurf.IsNull() ||
-       eos.ShapeType() != TopAbs_FACE ||
-       eos._edgeForOffset == 0 ||
-       eos._edgeForOffset->Is( _LayerEdge::BLOCKED ))
+  _EdgesOnShape * eof = & eos;
+  if ( eos.ShapeType() != TopAbs_FACE ) // eos is a boundary of C1 FACE, look for the FACE eos
+  {
+    eof = 0;
+    for ( size_t i = 0; i < eosC1.size() && !eof; ++i )
+    {
+      if ( eosC1[i]->_offsetSurf.IsNull() ||
+           eosC1[i]->ShapeType() != TopAbs_FACE ||
+           eosC1[i]->_edgeForOffset == 0 ||
+           eosC1[i]->_edgeForOffset->Is( _LayerEdge::BLOCKED ))
+        continue;
+      if ( SMESH_MesherHelper::IsSubShape( eos._shape, eosC1[i]->_shape ))
+        eof = eosC1[i];
+    }
+  }
+  if ( !eof ||
+       eof->_offsetSurf.IsNull() ||
+       eof->ShapeType() != TopAbs_FACE ||
+       eof->_edgeForOffset == 0 ||
+       eof->_edgeForOffset->Is( _LayerEdge::BLOCKED ))
     return;
 
-  double preci = BRep_Tool::Tolerance( TopoDS::Face( eos._shape )), vol;
+  double preci = BRep_Tool::Tolerance( TopoDS::Face( eof->_shape )), vol;
   for ( size_t i = 0; i < eos._edges.size(); ++i )
   {
     _LayerEdge* edge = eos._edges[i];
     edge->Unset( _LayerEdge::MARKED );
     if ( edge->Is( _LayerEdge::BLOCKED ) || !edge->_curvature )
       continue;
-    if ( !moveAll && !edge->Is( _LayerEdge::MOVED ))
+    if ( moveAll == _LayerEdge::UPD_NORMAL_CONV )
+    {
+      if ( !edge->Is( _LayerEdge::UPD_NORMAL_CONV ))
         continue;
+    }
+    else if ( !moveAll && !edge->Is( _LayerEdge::MOVED ))
+      continue;
 
     int nbBlockedAround = 0;
     for ( size_t iN = 0; iN < edge->_neibors.size(); ++iN )
@@ -5037,12 +5363,12 @@ void _ViscousBuilder::putOnOffsetSurface( _EdgesOnShape& eos,
       continue;
 
     gp_Pnt tgtP = SMESH_TNodeXYZ( edge->_nodes.back() );
-    gp_Pnt2d uv = eos._offsetSurf->NextValueOfUV( edge->_curvature->_uv, tgtP, preci );
-    if ( eos._offsetSurf->Gap() > edge->_len ) continue; // NextValueOfUV() bug 
+    gp_Pnt2d uv = eof->_offsetSurf->NextValueOfUV( edge->_curvature->_uv, tgtP, preci );
+    if ( eof->_offsetSurf->Gap() > edge->_len ) continue; // NextValueOfUV() bug
     edge->_curvature->_uv = uv;
-    if ( eos._offsetSurf->Gap() < 10 * preci ) continue; // same pos
+    if ( eof->_offsetSurf->Gap() < 10 * preci ) continue; // same pos
 
-    gp_XYZ  newP = eos._offsetSurf->Value( uv ).XYZ();
+    gp_XYZ  newP = eof->_offsetSurf->Value( uv ).XYZ();
     gp_XYZ prevP = edge->PrevCheckPos();
     bool      ok = true;
     if ( !moveAll )
@@ -5057,9 +5383,15 @@ void _ViscousBuilder::putOnOffsetSurface( _EdgesOnShape& eos,
       edge->_pos.back() = newP;
 
       edge->Set( _LayerEdge::MARKED );
+      if ( moveAll == _LayerEdge::UPD_NORMAL_CONV )
+      {
+        edge->_normal = ( newP - prevP ).Normalized();
+      }
     }
   }
 
+
+
 #ifdef _DEBUG_
   // dumpMove() for debug
   size_t i = 0;
@@ -5068,7 +5400,7 @@ void _ViscousBuilder::putOnOffsetSurface( _EdgesOnShape& eos,
       break;
   if ( i < eos._edges.size() )
   {
-    dumpFunction(SMESH_Comment("putOnOffsetSurface_F") << eos._shapeID
+    dumpFunction(SMESH_Comment("putOnOffsetSurface_S") << eos._shapeID
                  << "_InfStep" << infStep << "_" << smooStep );
     for ( ; i < eos._edges.size(); ++i )
     {
@@ -5078,6 +5410,26 @@ void _ViscousBuilder::putOnOffsetSurface( _EdgesOnShape& eos,
     dumpFunctionEnd();
   }
 #endif
+
+  _ConvexFace* cnvFace;
+  if ( moveAll != _LayerEdge::UPD_NORMAL_CONV &&
+       eos.ShapeType() == TopAbs_FACE &&
+       (cnvFace = eos.GetData().GetConvexFace( eos._shapeID )) &&
+       !cnvFace->_normalsFixedOnBorders )
+  {
+    // put on the surface nodes built on FACE boundaries
+    SMESH_subMeshIteratorPtr smIt = eos._subMesh->getDependsOnIterator(/*includeSelf=*/false);
+    while ( smIt->more() )
+    {
+      SMESH_subMesh* sm = smIt->next();
+      _EdgesOnShape* subEOS = eos.GetData().GetShapeEdges( sm->GetId() );
+      if ( !subEOS->_sWOL.IsNull() ) continue;
+      if ( std::find( eosC1.begin(), eosC1.end(), subEOS ) != eosC1.end() ) continue;
+
+      putOnOffsetSurface( *subEOS, infStep, eosC1, smooStep, _LayerEdge::UPD_NORMAL_CONV );
+    }
+    cnvFace->_normalsFixedOnBorders = true;
+  }
 }
 
 //================================================================================
@@ -5175,6 +5527,106 @@ Handle(Geom_Curve) _Smoother1D::CurveForSmooth( const TopoDS_Edge&  E,
   return Handle(Geom_Curve)();
 }
 
+//================================================================================
+/*!
+ * \brief Smooth edges on EDGE
+ */
+//================================================================================
+
+bool _Smoother1D::Perform(_SolidData&                    data,
+                          Handle(ShapeAnalysis_Surface)& surface,
+                          const TopoDS_Face&             F,
+                          SMESH_MesherHelper&            helper )
+{
+  if ( _leParams.empty() || ( !isAnalytic() && _offPoints.empty() ))
+    prepare( data );
+
+  findEdgesToSmooth();
+  if ( isAnalytic() )
+    return smoothAnalyticEdge( data, surface, F, helper );
+  else
+    return smoothComplexEdge ( data, surface, F, helper );
+}
+
+//================================================================================
+/*!
+ * \brief Find edges to smooth
+ */
+//================================================================================
+
+void _Smoother1D::findEdgesToSmooth()
+{
+  _LayerEdge* leOnV[2] = { getLEdgeOnV(0), getLEdgeOnV(1) };
+  for ( int iEnd = 0; iEnd < 2; ++iEnd )
+    if ( leOnV[iEnd]->Is( _LayerEdge::NORMAL_UPDATED ))
+      _leOnV[iEnd]._cosin = Abs( _edgeDir[iEnd].Normalized() * leOnV[iEnd]->_normal );
+
+  _eToSmooth[0].first = _eToSmooth[0].second = 0;
+
+  for ( size_t i = 0; i < _eos.size(); ++i )
+  {
+    if ( !_eos[i]->Is( _LayerEdge::TO_SMOOTH ))
+    {
+      if ( needSmoothing( _leOnV[0]._cosin, _eos[i]->_len, _curveLen * _leParams[i] ) ||
+           isToSmooth( i ))
+        _eos[i]->Set( _LayerEdge::TO_SMOOTH );
+      else
+        break;
+    }
+    _eToSmooth[0].second = i+1;
+  }
+
+  _eToSmooth[1].first = _eToSmooth[1].second = _eos.size();
+
+  for ( int i = _eos.size() - 1; i >= _eToSmooth[0].second; --i )
+  {
+    if ( !_eos[i]->Is( _LayerEdge::TO_SMOOTH ))
+    {
+      if ( needSmoothing( _leOnV[1]._cosin, _eos[i]->_len, _curveLen * ( 1.-_leParams[i] )) ||
+           isToSmooth( i ))
+        _eos[i]->Set( _LayerEdge::TO_SMOOTH );
+      else
+        break;
+    }
+    _eToSmooth[1].first = i;
+  }
+}
+
+//================================================================================
+/*!
+ * \brief Check if iE-th _LayerEdge needs smoothing
+ */
+//================================================================================
+
+bool _Smoother1D::isToSmooth( int iE )
+{
+  SMESH_NodeXYZ pi( _eos[iE]->_nodes[0] );
+  SMESH_NodeXYZ p0( _eos[iE]->_2neibors->srcNode(0) );
+  SMESH_NodeXYZ p1( _eos[iE]->_2neibors->srcNode(1) );
+  gp_XYZ       seg0 = pi - p0;
+  gp_XYZ       seg1 = p1 - pi;
+  gp_XYZ    tangent =  seg0 + seg1;
+  double tangentLen = tangent.Modulus();
+  double  segMinLen = Min( seg0.Modulus(), seg1.Modulus() );
+  if ( tangentLen < std::numeric_limits<double>::min() )
+    return false;
+  tangent /= tangentLen;
+
+  for ( size_t i = 0; i < _eos[iE]->_neibors.size(); ++i )
+  {
+    _LayerEdge* ne = _eos[iE]->_neibors[i];
+    if ( !ne->Is( _LayerEdge::TO_SMOOTH ) ||
+         ne->_nodes.size() < 2 ||
+         ne->_nodes[0]->GetPosition()->GetDim() != 2 )
+      continue;
+    gp_XYZ edgeVec = SMESH_NodeXYZ( ne->_nodes.back() ) - SMESH_NodeXYZ( ne->_nodes[0] );
+    double    proj = edgeVec * tangent;
+    if ( needSmoothing( 1., proj, segMinLen ))
+      return true;
+  }
+  return false;
+}
+
 //================================================================================
 /*!
  * \brief smooth _LayerEdge's on a staight EDGE or circular EDGE
@@ -5188,75 +5640,121 @@ bool _Smoother1D::smoothAnalyticEdge( _SolidData&                    data,
 {
   if ( !isAnalytic() ) return false;
 
-  const size_t iFrom = 0, iTo = _eos._edges.size();
+  size_t iFrom = 0, iTo = _eos._edges.size();
 
   if ( _anaCurve->IsKind( STANDARD_TYPE( Geom_Line )))
   {
     if ( F.IsNull() ) // 3D
     {
-      SMESH_TNodeXYZ p0   ( _eos._edges[iFrom]->_2neibors->tgtNode(0) );
-      SMESH_TNodeXYZ p1   ( _eos._edges[iTo-1]->_2neibors->tgtNode(1) );
       SMESH_TNodeXYZ pSrc0( _eos._edges[iFrom]->_2neibors->srcNode(0) );
       SMESH_TNodeXYZ pSrc1( _eos._edges[iTo-1]->_2neibors->srcNode(1) );
-      gp_XYZ newPos;
-      for ( size_t i = iFrom; i < iTo; ++i )
-      {
-        _LayerEdge*       edge = _eos._edges[i];
-        SMDS_MeshNode* tgtNode = const_cast<SMDS_MeshNode*>( edge->_nodes.back() );
-        newPos = p0 * ( 1. - _leParams[i] ) + p1 * _leParams[i];
-
-        if ( _eos._edges[i]->Is( _LayerEdge::NORMAL_UPDATED ))
-        {
-          gp_XYZ curPos  = SMESH_TNodeXYZ ( tgtNode );
-          gp_XYZ lineDir = pSrc1 - pSrc0;
-          double   shift = ( lineDir * ( newPos - pSrc0 ) -
-                             lineDir * ( curPos - pSrc0 ));
-          newPos = curPos + lineDir * shift / lineDir.SquareModulus();
-        }
-        if ( _eos._edges[i]->Is( _LayerEdge::BLOCKED ))
+      //const   gp_XYZ lineDir = pSrc1 - pSrc0;
+      //_LayerEdge* vLE0 = getLEdgeOnV( 0 );
+      //_LayerEdge* vLE1 = getLEdgeOnV( 1 );
+      // bool shiftOnly = ( vLE0->Is( _LayerEdge::NORMAL_UPDATED ) ||
+      //                    vLE0->Is( _LayerEdge::BLOCKED ) ||
+      //                    vLE1->Is( _LayerEdge::NORMAL_UPDATED ) ||
+      //                    vLE1->Is( _LayerEdge::BLOCKED ));
+      for ( int iEnd = 0; iEnd < 2; ++iEnd )
+      {
+        iFrom = _eToSmooth[ iEnd ].first, iTo = _eToSmooth[ iEnd ].second;
+        if ( iFrom >= iTo ) continue;
+        SMESH_TNodeXYZ p0( _eos[iFrom]->_2neibors->tgtNode(0) );
+        SMESH_TNodeXYZ p1( _eos[iTo-1]->_2neibors->tgtNode(1) );
+        double param0 = ( iFrom == 0 ) ? 0. : _leParams[ iFrom-1 ];
+        double param1 = _leParams[ iTo ];
+        for ( size_t i = iFrom; i < iTo; ++i )
         {
-          SMESH_TNodeXYZ pSrc( edge->_nodes[0] );
-          double curThick = pSrc.SquareDistance( tgtNode );
-          double newThink = ( pSrc - newPos ).SquareModulus();
-          if ( newThink > curThick )
-            continue;
+          _LayerEdge*       edge = _eos[i];
+          SMDS_MeshNode* tgtNode = const_cast<SMDS_MeshNode*>( edge->_nodes.back() );
+          double           param = ( _leParams[i] - param0 ) / ( param1 - param0 );
+          gp_XYZ          newPos = p0 * ( 1. - param ) + p1 * param;
+
+          // if ( shiftOnly || edge->Is( _LayerEdge::NORMAL_UPDATED ))
+          // {
+          //   gp_XYZ curPos = SMESH_TNodeXYZ ( tgtNode );
+          //   double  shift = ( lineDir * ( newPos - pSrc0 ) -
+          //                     lineDir * ( curPos - pSrc0 ));
+          //   newPos = curPos + lineDir * shift / lineDir.SquareModulus();
+          // }
+          if ( edge->Is( _LayerEdge::BLOCKED ))
+          {
+            SMESH_TNodeXYZ pSrc( edge->_nodes[0] );
+            double curThick = pSrc.SquareDistance( tgtNode );
+            double newThink = ( pSrc - newPos ).SquareModulus();
+            if ( newThink > curThick )
+              continue;
+          }
+          edge->_pos.back() = newPos;
+          tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() );
+          dumpMove( tgtNode );
         }
-        edge->_pos.back() = newPos;
-        tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() );
-        dumpMove( tgtNode );
       }
     }
     else // 2D
     {
-      _LayerEdge* e0 = getLEdgeOnV( 0 );
-      _LayerEdge* e1 = getLEdgeOnV( 1 );
-      gp_XY uv0 = e0->LastUV( F, *data.GetShapeEdges( e0 ));
-      gp_XY uv1 = e1->LastUV( F, *data.GetShapeEdges( e1 ));
-      if ( e0->_nodes.back() == e1->_nodes.back() ) // closed edge
+      _LayerEdge* eV0 = getLEdgeOnV( 0 );
+      _LayerEdge* eV1 = getLEdgeOnV( 1 );
+      gp_XY      uvV0 = eV0->LastUV( F, *data.GetShapeEdges( eV0 ));
+      gp_XY      uvV1 = eV1->LastUV( F, *data.GetShapeEdges( eV1 ));
+      if ( eV0->_nodes.back() == eV1->_nodes.back() ) // closed edge
       {
         int iPeriodic = helper.GetPeriodicIndex();
         if ( iPeriodic == 1 || iPeriodic == 2 )
         {
-          uv1.SetCoord( iPeriodic, helper.GetOtherParam( uv1.Coord( iPeriodic )));
-          if ( uv0.Coord( iPeriodic ) > uv1.Coord( iPeriodic ))
-            std::swap( uv0, uv1 );
+          uvV1.SetCoord( iPeriodic, helper.GetOtherParam( uvV1.Coord( iPeriodic )));
+          if ( uvV0.Coord( iPeriodic ) > uvV1.Coord( iPeriodic ))
+            std::swap( uvV0, uvV1 );
         }
       }
-      const gp_XY rangeUV = uv1 - uv0;
-      for ( size_t i = iFrom; i < iTo; ++i )
-      {
-        if ( _eos._edges[i]->Is( _LayerEdge::BLOCKED )) continue;
-        gp_XY newUV = uv0 + _leParams[i] * rangeUV;
-        _eos._edges[i]->_pos.back().SetCoord( newUV.X(), newUV.Y(), 0 );
+      for ( int iEnd = 0; iEnd < 2; ++iEnd )
+      {
+        iFrom = _eToSmooth[ iEnd ].first, iTo = _eToSmooth[ iEnd ].second;
+        if ( iFrom >= iTo ) continue;
+        _LayerEdge* e0 = _eos[iFrom]->_2neibors->_edges[0];
+        _LayerEdge* e1 = _eos[iTo-1]->_2neibors->_edges[1];
+        gp_XY      uv0 = ( e0 == eV0 ) ? uvV0 : e0->LastUV( F, _eos );
+        gp_XY      uv1 = ( e1 == eV1 ) ? uvV1 : e1->LastUV( F, _eos );
+        double  param0 = ( iFrom == 0 ) ? 0. : _leParams[ iFrom-1 ];
+        double  param1 = _leParams[ iTo ];
+        gp_XY  rangeUV = uv1 - uv0;
+        for ( size_t i = iFrom; i < iTo; ++i )
+        {
+          if ( _eos[i]->Is( _LayerEdge::BLOCKED )) continue;
+          double param = ( _leParams[i] - param0 ) / ( param1 - param0 );
+          gp_XY newUV = uv0 + param * rangeUV;
 
-        gp_Pnt newPos = surface->Value( newUV.X(), newUV.Y() );
-        SMDS_MeshNode* tgtNode = const_cast<SMDS_MeshNode*>( _eos._edges[i]->_nodes.back() );
-        tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() );
-        dumpMove( tgtNode );
+          gp_Pnt newPos = surface->Value( newUV.X(), newUV.Y() );
+          SMDS_MeshNode* tgtNode = const_cast<SMDS_MeshNode*>( _eos[i]->_nodes.back() );
+          tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() );
+          dumpMove( tgtNode );
 
-        SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( tgtNode->GetPosition() );
-        pos->SetUParameter( newUV.X() );
-        pos->SetVParameter( newUV.Y() );
+          SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( tgtNode->GetPosition() );
+          pos->SetUParameter( newUV.X() );
+          pos->SetVParameter( newUV.Y() );
+
+          gp_XYZ newUV0( newUV.X(), newUV.Y(), 0 );
+
+          if ( !_eos[i]->Is( _LayerEdge::SMOOTHED ))
+          {
+            _eos[i]->Set( _LayerEdge::SMOOTHED ); // to check in refine() (IPAL54237)
+            if ( _eos[i]->_pos.size() > 2 )
+            {
+              // modify previous positions to make _LayerEdge less sharply bent
+              vector<gp_XYZ>& uvVec = _eos[i]->_pos;
+              const gp_XYZ  uvShift = newUV0 - uvVec.back();
+              const double     len2 = ( uvVec.back() - uvVec[ 0 ] ).SquareModulus();
+              int iPrev = uvVec.size() - 2;
+              while ( iPrev > 0 )
+              {
+                double r = ( uvVec[ iPrev ] - uvVec[0] ).SquareModulus() / len2;
+                uvVec[ iPrev ] += uvShift * r;
+                --iPrev;
+              }
+            }
+          }
+          _eos[i]->_pos.back() = newUV0;
+        }
       }
     }
     return true;
@@ -5295,9 +5793,10 @@ bool _Smoother1D::smoothAnalyticEdge( _SolidData&                    data,
       if ( uLast < 0 )
         uLast += 2 * M_PI;
       
-      for ( size_t i = iFrom; i < iTo; ++i )
+      for ( size_t i = 0; i < _eos.size(); ++i )
       {
-        if ( _eos._edges[i]->Is( _LayerEdge::BLOCKED )) continue;
+        if ( _eos[i]->Is( _LayerEdge::BLOCKED )) continue;
+        //if ( !_eos[i]->Is( _LayerEdge::TO_SMOOTH )) continue;
         double u = uLast * _leParams[i];
         gp_Pnt p = ElCLib::Value( u, newCirc );
         _eos._edges[i]->_pos.back() = p.XYZ();
@@ -5329,9 +5828,10 @@ bool _Smoother1D::smoothAnalyticEdge( _SolidData&                    data,
 
       gp_Ax2d   axis( center, vec0 );
       gp_Circ2d circ( axis, radius );
-      for ( size_t i = iFrom; i < iTo; ++i )
+      for ( size_t i = 0; i < _eos.size(); ++i )
       {
-        if ( _eos._edges[i]->Is( _LayerEdge::BLOCKED )) continue;
+        if ( _eos[i]->Is( _LayerEdge::BLOCKED )) continue;
+        //if ( !_eos[i]->Is( _LayerEdge::TO_SMOOTH )) continue;
         double    newU = uLast * _leParams[i];
         gp_Pnt2d newUV = ElCLib::Value( newU, circ );
         _eos._edges[i]->_pos.back().SetCoord( newUV.X(), newUV.Y(), 0 );
@@ -5344,6 +5844,8 @@ bool _Smoother1D::smoothAnalyticEdge( _SolidData&                    data,
         SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( tgtNode->GetPosition() );
         pos->SetUParameter( newUV.X() );
         pos->SetVParameter( newUV.Y() );
+
+        _eos[i]->Set( _LayerEdge::SMOOTHED ); // to check in refine() (IPAL54237)
       }
     }
     return true;
@@ -5366,11 +5868,15 @@ bool _Smoother1D::smoothComplexEdge( _SolidData&                    data,
   if ( _offPoints.empty() )
     return false;
 
-  // move _offPoints to positions along normals of _LayerEdge's
+  // ----------------------------------------------
+  // move _offPoints along normals of _LayerEdge's
+  // ----------------------------------------------
 
   _LayerEdge* e[2] = { getLEdgeOnV(0), getLEdgeOnV(1) };
-  if ( e[0]->Is( _LayerEdge::NORMAL_UPDATED )) setNormalOnV( 0, helper );
-  if ( e[1]->Is( _LayerEdge::NORMAL_UPDATED )) setNormalOnV( 1, helper );
+  if ( e[0]->Is( _LayerEdge::NORMAL_UPDATED ))
+    _leOnV[0]._normal = getNormalNormal( e[0]->_normal, _edgeDir[0] );
+  if ( e[1]->Is( _LayerEdge::NORMAL_UPDATED )) 
+    _leOnV[1]._normal = getNormalNormal( e[1]->_normal, _edgeDir[1] );
   _leOnV[0]._len = e[0]->_len;
   _leOnV[1]._len = e[1]->_len;
   for ( size_t i = 0; i < _offPoints.size(); i++ )
@@ -5382,12 +5888,15 @@ bool _Smoother1D::smoothComplexEdge( _SolidData&                    data,
     gp_XYZ  avgNorm = ( e0->_normal    * w0 + e1->_normal    * w1 ).Normalized();
     double  avgLen  = ( e0->_len       * w0 + e1->_len       * w1 );
     double  avgFact = ( e0->_lenFactor * w0 + e1->_lenFactor * w1 );
+    if ( e0->Is( _LayerEdge::NORMAL_UPDATED ) ||
+         e1->Is( _LayerEdge::NORMAL_UPDATED ))
+      avgNorm = getNormalNormal( avgNorm, _offPoints[i]._edgeDir );
 
     _offPoints[i]._xyz += avgNorm * ( avgLen - _offPoints[i]._len ) * avgFact;
     _offPoints[i]._len  = avgLen;
   }
 
-  double fTol;
+  double fTol = 0;
   if ( !surface.IsNull() ) // project _offPoints to the FACE
   {
     fTol = 100 * BRep_Tool::Tolerance( F );
@@ -5405,7 +5914,13 @@ bool _Smoother1D::smoothComplexEdge( _SolidData&                    data,
     }
   }
 
+  // -----------------------------------------------------------------
   // project tgt nodes of extreme _LayerEdge's to the offset segments
+  // -----------------------------------------------------------------
+
+  const int updatedOrBlocked = _LayerEdge::NORMAL_UPDATED | _LayerEdge::BLOCKED;
+  if ( e[0]->Is( updatedOrBlocked )) _iSeg[0] = 0;
+  if ( e[1]->Is( updatedOrBlocked )) _iSeg[1] = _offPoints.size()-2;
 
   gp_Pnt pExtreme[2], pProj[2];
   for ( int is2nd = 0; is2nd < 2; ++is2nd )
@@ -5414,26 +5929,26 @@ bool _Smoother1D::smoothComplexEdge( _SolidData&                    data,
     int  i = _iSeg[ is2nd ];
     int di = is2nd ? -1 : +1;
     bool projected = false;
-    double uOnSeg, uOnSegDiff, uOnSegBestDiff = Precision::Infinite(), uOnSegPrevDiff;
+    double uOnSeg, distMin = Precision::Infinite(), dist, distPrev = 0;
     int nbWorse = 0;
     do {
       gp_Vec v0p( _offPoints[i]._xyz, pExtreme[ is2nd ]    );
       gp_Vec v01( _offPoints[i]._xyz, _offPoints[i+1]._xyz );
-      uOnSeg     = ( v0p * v01 ) / v01.SquareMagnitude();
-      uOnSegDiff = Abs( uOnSeg - 0.5 );
-      projected  = ( uOnSegDiff <= 0.5 );
-      if ( uOnSegDiff < uOnSegBestDiff )
+      uOnSeg     = ( v0p * v01 ) / v01.SquareMagnitude();  // param [0,1] along v01
+      projected  = ( Abs( uOnSeg - 0.5 ) <= 0.5 );
+      dist       =  pExtreme[ is2nd ].SquareDistance( _offPoints[ i + ( uOnSeg > 0.5 )]._xyz );
+      if ( dist < distMin || projected )
       {
         _iSeg[ is2nd ] = i;
         pProj[ is2nd ] = _offPoints[i]._xyz + ( v01 * uOnSeg ).XYZ();
-        uOnSegBestDiff = uOnSegDiff;
+        distMin = dist;
       }
-      else if ( uOnSegDiff > uOnSegPrevDiff )
+      else if ( dist > distPrev )
       {
         if ( ++nbWorse > 3 ) // avoid projection to the middle of a closed EDGE
           break;
       }
-      uOnSegPrevDiff = uOnSegDiff;
+      distPrev = dist;
       i += di;
     }
     while ( !projected &&
@@ -5456,20 +5971,36 @@ bool _Smoother1D::smoothComplexEdge( _SolidData&                    data,
     return false;
   }
 
+  // adjust length of extreme LE (test viscous_layers_01/B7)
+  gp_Vec vDiv0( pExtreme[0], pProj[0] );
+  gp_Vec vDiv1( pExtreme[1], pProj[1] );
+  double d0 = vDiv0.Magnitude();
+  double d1 = vDiv1.Magnitude();
+  if ( e[0]->Is( _LayerEdge::BLOCKED )) {
+    if ( e[0]->_normal * vDiv0.XYZ() < 0 ) e[0]->_len += d0;
+    else                                   e[0]->_len -= d0;
+  }
+  if ( e[1]->Is( _LayerEdge::BLOCKED )) {
+    if ( e[1]->_normal * vDiv1.XYZ() < 0 ) e[1]->_len += d1;
+    else                                   e[1]->_len -= d1;
+  }
+
+  // ---------------------------------------------------------------------------------
   // compute normalized length of the offset segments located between the projections
+  // ---------------------------------------------------------------------------------
 
   size_t iSeg = 0, nbSeg = _iSeg[1] - _iSeg[0] + 1;
   vector< double > len( nbSeg + 1 );
   len[ iSeg++ ] = 0;
-  len[ iSeg++ ] = pProj[ 0 ].Distance( _offPoints[ _iSeg[0]+1 ]._xyz );
+  len[ iSeg++ ] = pProj[ 0 ].Distance( _offPoints[ _iSeg[0]+1 ]._xyz )/* * e[0]->_lenFactor*/;
   for ( size_t i = _iSeg[0]+1; i <= _iSeg[1]; ++i, ++iSeg )
   {
     len[ iSeg ] = len[ iSeg-1 ] + _offPoints[i].Distance( _offPoints[i+1] );
   }
-  len[ nbSeg ] -= pProj[ 1 ].Distance( _offPoints[ _iSeg[1]+1 ]._xyz );
+  len[ nbSeg ] -= pProj[ 1 ].Distance( _offPoints[ _iSeg[1]+1 ]._xyz )/* * e[1]->_lenFactor*/;
 
-  double d0 = pProj[0].Distance( pExtreme[0]);
-  double d1 = pProj[1].Distance( pExtreme[1]);
+  // d0 *= e[0]->_lenFactor;
+  // d1 *= e[1]->_lenFactor;
   double fullLen = len.back() - d0 - d1;
   for ( iSeg = 0; iSeg < len.size(); ++iSeg )
     len[iSeg] = ( len[iSeg] - d0 ) / fullLen;
@@ -5480,12 +6011,15 @@ bool _Smoother1D::smoothComplexEdge( _SolidData&                    data,
   _offPoints[ _iSeg[0]   ]._xyz = pExtreme[0].XYZ();
   _offPoints[ _iSeg[1]+ 1]._xyz = pExtreme[1].XYZ();
 
+  // -------------------------------------------------------------
   // distribute tgt nodes of _LayerEdge's between the projections
+  // -------------------------------------------------------------
 
   iSeg = 0;
-  for ( size_t i = 0; i < _eos._edges.size(); ++i )
+  for ( size_t i = 0; i < _eos.size(); ++i )
   {
-    if ( _eos._edges[i]->Is( _LayerEdge::BLOCKED )) continue;
+    if ( _eos[i]->Is( _LayerEdge::BLOCKED )) continue;
+    //if ( !_eos[i]->Is( _LayerEdge::TO_SMOOTH )) continue;
     while ( iSeg+2 < len.size() && _leParams[i] > len[ iSeg+1 ] )
       iSeg++;
     double r = ( _leParams[i] - len[ iSeg ]) / ( len[ iSeg+1 ] - len[ iSeg ]);
@@ -5494,17 +6028,17 @@ bool _Smoother1D::smoothComplexEdge( _SolidData&                    data,
 
     if ( surface.IsNull() )
     {
-      _eos._edges[i]->_pos.back() = p;
+      _eos[i]->_pos.back() = p;
     }
     else // project a new node position to a FACE
     {
-      gp_Pnt2d uv ( _eos._edges[i]->_pos.back().X(), _eos._edges[i]->_pos.back().Y() );
+      gp_Pnt2d uv ( _eos[i]->_pos.back().X(), _eos[i]->_pos.back().Y() );
       gp_Pnt2d uv2( surface->NextValueOfUV( uv, p, fTol ));
 
       p = surface->Value( uv2 ).XYZ();
-      _eos._edges[i]->_pos.back().SetCoord( uv2.X(), uv2.Y(), 0 );
+      _eos[i]->_pos.back().SetCoord( uv2.X(), uv2.Y(), 0 );
     }
-    SMDS_MeshNode* tgtNode = const_cast<SMDS_MeshNode*>( _eos._edges[i]->_nodes.back() );
+    SMDS_MeshNode* tgtNode = const_cast<SMDS_MeshNode*>( _eos[i]->_nodes.back() );
     tgtNode->setXYZ( p.X(), p.Y(), p.Z() );
     dumpMove( tgtNode );
   }
@@ -5529,71 +6063,43 @@ void _Smoother1D::prepare(_SolidData& data)
   // sort _LayerEdge's by position on the EDGE
   data.SortOnEdge( E, _eos._edges );
 
-  SMESH_MesherHelper& helper = data.GetHelper();
-
   // compute normalized param of _eos._edges on EDGE
   _leParams.resize( _eos._edges.size() + 1 );
   {
-    double curLen, prevLen = _leParams[0] = 1.0;
+    double curLen;
     gp_Pnt pPrev = SMESH_TNodeXYZ( getLEdgeOnV( 0 )->_nodes[0] );
     _leParams[0] = 0;
     for ( size_t i = 0; i < _eos._edges.size(); ++i )
     {
-      gp_Pnt p = SMESH_TNodeXYZ( _eos._edges[i]->_nodes[0] );
-      //curLen = prevLen * _eos._edges[i]->_2neibors->_wgt[1] / _eos._edges[i]->_2neibors->_wgt[0];
-      curLen = p.Distance( pPrev );
+      gp_Pnt p       = SMESH_TNodeXYZ( _eos._edges[i]->_nodes[0] );
+      curLen         = p.Distance( pPrev );
       _leParams[i+1] = _leParams[i] + curLen;
-      prevLen = curLen;
-      pPrev = p;
+      pPrev          = p;
     }
     double fullLen = _leParams.back() + pPrev.Distance( SMESH_TNodeXYZ( getLEdgeOnV(1)->_nodes[0]));
     for ( size_t i = 0; i < _leParams.size()-1; ++i )
       _leParams[i] = _leParams[i+1] / fullLen;
+    _leParams.back() = 1.;
   }
 
-  // find intersection of neighbor _LayerEdge's to limit _maxLen
-  // according to EDGE curvature (IPAL52648)
-  _LayerEdge* e0 = _eos._edges[0];
-  for ( size_t i = 1; i < _eos._edges.size(); ++i )
-  {
-    _LayerEdge* ei = _eos._edges[i];
-    gp_XYZ plnNorm = e0->_normal ^ ei->_normal;
-    gp_XYZ   perp0 = e0->_normal ^ plnNorm;
-    double   dot0i = perp0 * ei->_normal;
-    if ( Abs( dot0i ) > std::numeric_limits<double>::min() )
-    {
-      SMESH_TNodeXYZ srci( ei->_nodes[0] ), src0( e0->_nodes[0] );
-      double ui = ( perp0 * ( src0 - srci )) / dot0i;
-      if ( ui > 0 )
-      {
-        ei->_maxLen = Min(  ei->_maxLen, 0.75 * ui / ei->_lenFactor );
-        if ( ei->_maxLen < ei->_len )
-        {
-          ei->InvalidateStep( ei->NbSteps(), _eos, /*restoreLength=*/true  );
-          ei->SetNewLength( ei->_maxLen, _eos, helper );
-          ei->Block( data );
-        }
-        gp_Pnt pi = srci + ei->_normal * ui;
-        double u0 = pi.Distance( src0 );
-        e0->_maxLen = Min(  e0->_maxLen, 0.75 * u0 / e0->_lenFactor );
-        if ( e0->_maxLen < e0->_len )
-        {
-          e0->InvalidateStep( e0->NbSteps(), _eos, /*restoreLength=*/true  );
-          e0->SetNewLength( e0->_maxLen, _eos, helper );
-          e0->Block( data );
-        }
-      }
-    }
-    e0 = ei;
-  }
-    
+  _LayerEdge* leOnV[2] = { getLEdgeOnV(0), getLEdgeOnV(1) };
+
+  // get cosin to use in findEdgesToSmooth()
+  _edgeDir[0] = getEdgeDir( E, leOnV[0]->_nodes[0], data.GetHelper() );
+  _edgeDir[1] = getEdgeDir( E, leOnV[1]->_nodes[0], data.GetHelper() );
+  _leOnV[0]._cosin = Abs( leOnV[0]->_cosin );
+  _leOnV[1]._cosin = Abs( leOnV[1]->_cosin );
+  if ( _eos._sWOL.IsNull() ) // 3D
+    for ( int iEnd = 0; iEnd < 2; ++iEnd )
+      _leOnV[iEnd]._cosin = Abs( _edgeDir[iEnd].Normalized() * leOnV[iEnd]->_normal );
+
   if ( isAnalytic() )
     return;
 
   // divide E to have offset segments with low deflection
   BRepAdaptor_Curve c3dAdaptor( E );
-  const double curDeflect = 0.1; //0.3; // 0.01; // Curvature deflection
-  const double angDeflect = 0.1; //0.2; // 0.09; // Angular deflection
+  const double curDeflect = 0.1; //0.01; // Curvature deflection == |p1p2]*sin(p1p2,p1pM)
+  const double angDeflect = 0.1; //0.09; // Angular deflection == sin(p1pM,pMp2)
   GCPnts_TangentialDeflection discret(c3dAdaptor, angDeflect, curDeflect);
   if ( discret.NbPoints() <= 2 )
   {
@@ -5601,18 +6107,42 @@ void _Smoother1D::prepare(_SolidData& data)
     return;
   }
 
-  const double edgeLen = SMESH_Algo::EdgeLength( E );
-  const double u0      = c3dAdaptor.FirstParameter();
-  _offPoints.resize( discret.NbPoints() );
-  for ( size_t i = 0; i < _offPoints.size(); i++ )
+  const double u0 = c3dAdaptor.FirstParameter();
+  gp_Pnt p; gp_Vec tangent;
+  if ( discret.NbPoints() >= (int) _eos.size() + 2 )
   {
-    _offPoints[i]._xyz = discret.Value( i+1 ).XYZ();
-    // use OffPnt::_len to  TEMPORARY  store normalized param of an offset point
-    double u = discret.Parameter( i+1 );
-    _offPoints[i]._len = GCPnts_AbscissaPoint::Length( c3dAdaptor, u0, u ) / edgeLen;
+    _offPoints.resize( discret.NbPoints() );
+    for ( size_t i = 0; i < _offPoints.size(); i++ )
+    {
+      double u = discret.Parameter( i+1 );
+      c3dAdaptor.D1( u, p, tangent );
+      _offPoints[i]._xyz     = p.XYZ();
+      _offPoints[i]._edgeDir = tangent.XYZ();
+      _offPoints[i]._param = GCPnts_AbscissaPoint::Length( c3dAdaptor, u0, u ) / _curveLen;
+    }
   }
+  else
+  {
+    std::vector< double > params( _eos.size() + 2 );
 
-  _LayerEdge* leOnV[2] = { getLEdgeOnV(0), getLEdgeOnV(1) };
+    params[0]     = data.GetHelper().GetNodeU( E, leOnV[0]->_nodes[0] );
+    params.back() = data.GetHelper().GetNodeU( E, leOnV[1]->_nodes[0] );
+    for ( size_t i = 0; i < _eos.size(); i++ )
+      params[i+1] = data.GetHelper().GetNodeU( E, _eos[i]->_nodes[0] );
+
+    if ( params[1] > params[ _eos.size() ] )
+      std::reverse( params.begin() + 1, params.end() - 1 );
+
+    _offPoints.resize( _eos.size() + 2 );
+    for ( size_t i = 0; i < _offPoints.size(); i++ )
+    {
+      const double u = params[i];
+      c3dAdaptor.D1( u, p, tangent );
+      _offPoints[i]._xyz     = p.XYZ();
+      _offPoints[i]._edgeDir = tangent.XYZ();
+      _offPoints[i]._param = GCPnts_AbscissaPoint::Length( c3dAdaptor, u0, u ) / _curveLen;
+    }
+  }
 
   // set _2edges
   _offPoints    [0]._2edges.set( &_leOnV[0], &_leOnV[0], 0.5, 0.5 );
@@ -5627,7 +6157,7 @@ void _Smoother1D::prepare(_SolidData& data)
   {
     // find _LayerEdge's located before and after an offset point
     // (_eos._edges[ iLE ] is next after ePrev)
-    while ( iLE < _eos._edges.size() && _offPoints[i]._len > _leParams[ iLE ] )
+    while ( iLE < _eos._edges.size() && _offPoints[i]._param > _leParams[ iLE ] )
       ePrev = _eos._edges[ iLE++ ];
     eNext = ePrev->_2neibors->_edges[1];
 
@@ -5637,27 +6167,28 @@ void _Smoother1D::prepare(_SolidData& data)
     _offPoints[i]._2edges.set( ePrev, eNext, 1-r, r );
   }
 
-  int iLBO = _offPoints.size() - 2; // last but one
-  _offPoints[iLBO]._2edges._edges[1] = & _leOnV[1];
-
-  // {
-  //   TopoDS_Face face[2]; // FACEs sharing the EDGE
-  //   PShapeIteratorPtr fIt = helper.GetAncestors( _eos._shape, *helper.GetMesh(), TopAbs_FACE );
-  //   while ( const TopoDS_Shape* F = fIt->next() )
-  //   {
-  //     TGeomID fID = helper.GetMeshDS()->ShapeToIndex( *F );
-  //     if ( ! data._ignoreFaceIds.count( fID ))
-  //       face[ !face[0].IsNull() ] = *F;
-  //   }
-  //   if ( face[0].IsNull() ) return;
-  //   if ( face[1].IsNull() ) face[1] = face[0];
-  // }
-
+  // replace _LayerEdge's on VERTEX by _leOnV in _offPoints._2edges
+  for ( size_t i = 0; i < _offPoints.size(); i++ )
+    if ( _offPoints[i]._2edges._edges[0] == leOnV[0] )
+      _offPoints[i]._2edges._edges[0] = & _leOnV[0];
+    else break;
+  for ( size_t i = _offPoints.size()-1; i > 0; i-- )
+    if ( _offPoints[i]._2edges._edges[1] == leOnV[1] )
+      _offPoints[i]._2edges._edges[1] = & _leOnV[1];
+    else break;
 
   // set _normal of _leOnV[0] and _leOnV[1] to be normal to the EDGE
 
-  setNormalOnV( 0, data.GetHelper() );
-  setNormalOnV( 1, data.GetHelper() );
+  int iLBO = _offPoints.size() - 2; // last but one
+
+  if ( leOnV[ 0 ]->Is( _LayerEdge::MULTI_NORMAL ))
+    _leOnV[ 0 ]._normal = getNormalNormal( _eos._edges[1]->_normal, _edgeDir[0] );
+  else
+    _leOnV[ 0 ]._normal = getNormalNormal( leOnV[0]->_normal,       _edgeDir[0] );
+  if ( leOnV[ 1 ]->Is( _LayerEdge::MULTI_NORMAL ))
+    _leOnV[ 1 ]._normal = getNormalNormal( _eos._edges.back()->_normal, _edgeDir[1] );
+  else
+    _leOnV[ 1 ]._normal = getNormalNormal( leOnV[1]->_normal,           _edgeDir[1] );
   _leOnV[ 0 ]._len = 0;
   _leOnV[ 1 ]._len = 0;
   _leOnV[ 0 ]._lenFactor = _offPoints[1   ]._2edges._edges[1]->_lenFactor;
@@ -5691,22 +6222,21 @@ void _Smoother1D::prepare(_SolidData& data)
 
 //================================================================================
 /*!
- * \brief set _normal of _leOnV[is2nd] to be normal to the EDGE
+ * \brief return _normal of _leOnV[is2nd] normal to the EDGE
  */
 //================================================================================
 
-void _Smoother1D::setNormalOnV( const bool          is2nd,
-                                SMESH_MesherHelper& helper)
+gp_XYZ _Smoother1D::getNormalNormal( const gp_XYZ & normal,
+                                     const gp_XYZ&  edgeDir)
 {
-  _LayerEdge*    leOnV = getLEdgeOnV( is2nd );
-  const TopoDS_Edge& E = TopoDS::Edge( _eos._shape );
-  TopoDS_Shape       V = helper.GetSubShapeByNode( leOnV->_nodes[0], helper.GetMeshDS() );
-  gp_XYZ          eDir = getEdgeDir( E, TopoDS::Vertex( V ));
-  gp_XYZ         cross = leOnV->_normal ^ eDir;
-  gp_XYZ          norm = eDir ^ cross;
-  double          size = norm.Modulus();
+  gp_XYZ cross = normal ^ edgeDir;
+  gp_XYZ  norm = edgeDir ^ cross;
+  double  size = norm.Modulus();
+
+  // if ( size == 0 ) // MULTI_NORMAL _LayerEdge
+  //   return gp_XYZ( 1e-100, 1e-100, 1e-100 );
 
-  _leOnV[ is2nd ]._normal = norm / size;
+  return norm / size;
 }
 
 //================================================================================
@@ -5956,6 +6486,94 @@ void _SolidData::AddShapesToSmooth( const set< _EdgesOnShape* >& eosToSmooth,
     }
 }
 
+//================================================================================
+/*!
+ * \brief Limit _LayerEdge::_maxLen according to local curvature
+ */
+//================================================================================
+
+void _ViscousBuilder::limitMaxLenByCurvature( _SolidData& data, SMESH_MesherHelper& helper )
+{
+  // find intersection of neighbor _LayerEdge's to limit _maxLen
+  // according to local curvature (IPAL52648)
+
+  // This method must be called after findCollisionEdges() where _LayerEdge's
+  // get _lenFactor initialized in the case of eos._hyp.IsOffsetMethod()
+
+  for ( size_t iS = 0; iS < data._edgesOnShape.size(); ++iS )
+  {
+    _EdgesOnShape& eosI = data._edgesOnShape[iS];
+    if ( eosI._edges.empty() ) continue;
+    if ( !eosI._hyp.ToSmooth() )
+    {
+      for ( size_t i = 0; i < eosI._edges.size(); ++i )
+      {
+        _LayerEdge* eI = eosI._edges[i];
+        for ( size_t iN = 0; iN < eI->_neibors.size(); ++iN )
+        {
+          _LayerEdge* eN = eI->_neibors[iN];
+          if ( eI->_nodes[0]->GetID() < eN->_nodes[0]->GetID() ) // treat this pair once
+          {
+            _EdgesOnShape* eosN = data.GetShapeEdges( eN );
+            limitMaxLenByCurvature( eI, eN, eosI, *eosN, helper );
+          }
+        }
+      }
+    }
+    else if ( eosI.ShapeType() == TopAbs_EDGE )
+    {
+      const TopoDS_Edge& E = TopoDS::Edge( eosI._shape );
+      if ( SMESH_Algo::IsStraight( E, /*degenResult=*/true )) continue;
+
+      _LayerEdge* e0 = eosI._edges[0];
+      for ( size_t i = 1; i < eosI._edges.size(); ++i )
+      {
+        _LayerEdge* eI = eosI._edges[i];
+        limitMaxLenByCurvature( eI, e0, eosI, eosI, helper );
+        e0 = eI;
+      }
+    }
+  }
+}
+
+//================================================================================
+/*!
+ * \brief Limit _LayerEdge::_maxLen according to local curvature
+ */
+//================================================================================
+
+void _ViscousBuilder::limitMaxLenByCurvature( _LayerEdge*         e1,
+                                              _LayerEdge*         e2,
+                                              _EdgesOnShape&      eos1,
+                                              _EdgesOnShape&      eos2,
+                                              SMESH_MesherHelper& helper )
+{
+  gp_XYZ plnNorm = e1->_normal ^ e2->_normal;
+  double norSize = plnNorm.SquareModulus();
+  if ( norSize < std::numeric_limits<double>::min() )
+    return; // parallel normals
+
+  // find closest points of skew _LayerEdge's
+  SMESH_TNodeXYZ src1( e1->_nodes[0] ), src2( e2->_nodes[0] );
+  gp_XYZ dir12 = src2 - src1;
+  gp_XYZ perp1 = e1->_normal ^ plnNorm;
+  gp_XYZ perp2 = e2->_normal ^ plnNorm;
+  double  dot1 = perp2 * e1->_normal;
+  double  dot2 = perp1 * e2->_normal;
+  double    u1 =   ( perp2 * dir12 ) / dot1;
+  double    u2 = - ( perp1 * dir12 ) / dot2;
+  if ( u1 > 0 && u2 > 0 )
+  {
+    double ovl = ( u1 * e1->_normal * dir12 -
+                   u2 * e2->_normal * dir12 ) / dir12.SquareModulus();
+    if ( ovl > theSmoothThickToElemSizeRatio )
+    {    
+      e1->_maxLen = Min( e1->_maxLen, 0.75 * u1 / e1->_lenFactor );
+      e2->_maxLen = Min( e2->_maxLen, 0.75 * u2 / e2->_lenFactor );
+    }
+  }
+}
+
 //================================================================================
 /*!
  * \brief Fill data._collisionEdges
@@ -6034,10 +6652,10 @@ void _ViscousBuilder::findCollisionEdges( _SolidData& data, SMESH_MesherHelper&
   SMESHUtils::Deleter<SMESH_ElementSearcher> searcher
     ( SMESH_MeshAlgos::GetElementSearcher( *getMeshDS(), fIt ));
 
-  double dist1, dist2, segLen, eps;
+  double dist1, dist2, segLen, eps = 0.5;
   _CollisionEdges collEdges;
   vector< const SMDS_MeshElement* > suspectFaces;
-  const double angle30 = Cos( 30. * M_PI / 180. );
+  const double angle45 = Cos( 45. * M_PI / 180. );
 
   for ( size_t iS = 0; iS < data._edgesOnShape.size(); ++iS )
   {
@@ -6066,9 +6684,9 @@ void _ViscousBuilder::findCollisionEdges( _SolidData& data, SMESH_MesherHelper&
     // find intersecting _LayerEdge's
     for ( size_t i = 0; i < eos._edges.size(); ++i )
     {
+      if ( eos._edges[i]->Is( _LayerEdge::MULTI_NORMAL )) continue;
       _LayerEdge*   edge = eos._edges[i];
       gp_Ax1 lastSegment = edge->LastSegment( segLen, eos );
-      eps     = 0.5 * edge->_len;
       segLen *= 1.2;
 
       gp_Vec eSegDir0, eSegDir1;
@@ -6079,7 +6697,7 @@ void _ViscousBuilder::findCollisionEdges( _SolidData& data, SMESH_MesherHelper&
         eSegDir1 = SMESH_TNodeXYZ( edge->_2neibors->srcNode(1) ) - eP;
       }
       suspectFaces.clear();
-      searcher->GetElementsInSphere( SMESH_TNodeXYZ( edge->_nodes.back()), edge->_len,
+      searcher->GetElementsInSphere( SMESH_TNodeXYZ( edge->_nodes.back()), edge->_len * 2,
                                      SMDSAbs_Face, suspectFaces );
       collEdges._intEdges.clear();
       for ( size_t j = 0 ; j < suspectFaces.size(); ++j )
@@ -6108,10 +6726,10 @@ void _ViscousBuilder::findCollisionEdges( _SolidData& data, SMESH_MesherHelper&
         {
           // skip perpendicular EDGEs
           gp_Vec fSegDir  = SMESH_TNodeXYZ( f->_nn[0] ) - SMESH_TNodeXYZ( f->_nn[3] );
-          bool isParallel = ( isLessAngle( eSegDir0, fSegDir, angle30 ) ||
-                              isLessAngle( eSegDir1, fSegDir, angle30 ) ||
-                              isLessAngle( eSegDir0, fSegDir.Reversed(), angle30 ) ||
-                              isLessAngle( eSegDir1, fSegDir.Reversed(), angle30 ));
+          bool isParallel = ( isLessAngle( eSegDir0, fSegDir, angle45 ) ||
+                              isLessAngle( eSegDir1, fSegDir, angle45 ) ||
+                              isLessAngle( eSegDir0, fSegDir.Reversed(), angle45 ) ||
+                              isLessAngle( eSegDir1, fSegDir.Reversed(), angle45 ));
           if ( !isParallel )
             continue;
         }
@@ -6156,6 +6774,80 @@ void _ViscousBuilder::findCollisionEdges( _SolidData& data, SMESH_MesherHelper&
   }
 }
 
+//================================================================================
+/*!
+ * \brief Find _LayerEdge's located on boundary of a convex FACE whose normal
+ *        will be updated at each inflation step
+ */
+//================================================================================
+
+void _ViscousBuilder::findEdgesToUpdateNormalNearConvexFace( _ConvexFace &       convFace,
+                                                             _SolidData&         data,
+                                                             SMESH_MesherHelper& helper )
+{
+  const TGeomID convFaceID = getMeshDS()->ShapeToIndex( convFace._face );
+  const double       preci = BRep_Tool::Tolerance( convFace._face );
+  Handle(ShapeAnalysis_Surface) surface = helper.GetSurface( convFace._face );
+
+  bool edgesToUpdateFound = false;
+
+  map< TGeomID, _EdgesOnShape* >::iterator id2eos = convFace._subIdToEOS.begin();
+  for ( ; id2eos != convFace._subIdToEOS.end(); ++id2eos )
+  {
+    _EdgesOnShape& eos = * id2eos->second;
+    if ( !eos._sWOL.IsNull() ) continue;
+    if ( !eos._hyp.ToSmooth() ) continue;
+    for ( size_t i = 0; i < eos._edges.size(); ++i )
+    {
+      _LayerEdge* ledge = eos._edges[ i ];
+      if ( ledge->Is( _LayerEdge::UPD_NORMAL_CONV )) continue; // already checked
+      if ( ledge->Is( _LayerEdge::MULTI_NORMAL )) continue; // not inflatable
+
+      gp_XYZ tgtPos = ( SMESH_NodeXYZ( ledge->_nodes[0] ) +
+                        ledge->_normal * ledge->_lenFactor * ledge->_maxLen );
+
+      // the normal must be updated if distance from tgtPos to surface is less than
+      // target thickness
+
+      // find an initial UV for search of a projection of tgtPos to surface
+      const SMDS_MeshNode* nodeInFace = 0;
+      SMDS_ElemIteratorPtr fIt = ledge->_nodes[0]->GetInverseElementIterator(SMDSAbs_Face);
+      while ( fIt->more() && !nodeInFace )
+      {
+        const SMDS_MeshElement* f = fIt->next();
+        if ( convFaceID != f->getshapeId() ) continue;
+
+        SMDS_ElemIteratorPtr nIt = f->nodesIterator();
+        while ( nIt->more() && !nodeInFace )
+        {
+          const SMDS_MeshElement* n = nIt->next();
+          if ( n->getshapeId() == convFaceID )
+            nodeInFace = static_cast< const SMDS_MeshNode* >( n );
+        }
+      }
+      if ( !nodeInFace )
+        continue;
+      gp_XY uv = helper.GetNodeUV( convFace._face, nodeInFace );
+
+      // projection
+      surface->NextValueOfUV( uv, tgtPos, preci );
+      double  dist = surface->Gap();
+      if ( dist < 0.95 * ledge->_maxLen )
+      {
+        ledge->Set( _LayerEdge::UPD_NORMAL_CONV );
+        if ( !ledge->_curvature ) ledge->_curvature = new _Curvature;
+        ledge->_curvature->_uv.SetCoord( uv.X(), uv.Y() );
+        edgesToUpdateFound = true;
+      }
+    }
+  }
+
+  if ( !convFace._isTooCurved && edgesToUpdateFound )
+  {
+    data._convexFaces.insert( make_pair( convFaceID, convFace )).first->second;
+  }
+}
+
 //================================================================================
 /*!
  * \brief Modify normals of _LayerEdge's on EDGE's to avoid intersection with
@@ -6183,7 +6875,7 @@ bool _ViscousBuilder::updateNormals( _SolidData&         data,
 
   set< _EdgesOnShape* > shapesToSmooth, edgesNoAnaSmooth;
 
-  double segLen, dist1, dist2;
+  double segLen, dist1, dist2, dist;
   vector< pair< _LayerEdge*, double > > intEdgesDist;
   _TmpMeshFaceOnEdge quad( &zeroEdge, &zeroEdge, 0 );
 
@@ -6195,19 +6887,20 @@ bool _ViscousBuilder::updateNormals( _SolidData&         data,
     {
       _CollisionEdges& ce = data._collisionEdges[iE];
       _LayerEdge*   edge1 = ce._edge;
-      if ( !edge1 || edge1->Is( _LayerEdge::BLOCKED )) continue;
+      if ( !edge1 /*|| edge1->Is( _LayerEdge::BLOCKED )*/) continue;
       _EdgesOnShape* eos1 = data.GetShapeEdges( edge1 );
       if ( !eos1 ) continue;
 
       // detect intersections
       gp_Ax1 lastSeg = edge1->LastSegment( segLen, *eos1 );
-      double testLen = 1.5 * edge1->_maxLen; //2 + edge1->_len * edge1->_lenFactor;
-      double     eps = 0.5 * edge1->_len;
+      double testLen = 1.5 * edge1->_maxLen * edge1->_lenFactor;
+      double     eps = 0.5;
       intEdgesDist.clear();
       double minIntDist = Precision::Infinite();
       for ( size_t i = 0; i < ce._intEdges.size(); i += 2 )
       {
-        if ( ce._intEdges[i  ]->Is( _LayerEdge::BLOCKED ) ||
+        if ( edge1->Is( _LayerEdge::BLOCKED ) &&
+             ce._intEdges[i  ]->Is( _LayerEdge::BLOCKED ) &&
              ce._intEdges[i+1]->Is( _LayerEdge::BLOCKED ))
           continue;
         double dot  = edge1->_normal * quad.GetDir( ce._intEdges[i], ce._intEdges[i+1] );
@@ -6217,16 +6910,18 @@ bool _ViscousBuilder::updateNormals( _SolidData&         data,
         gp_XYZ pLast0 = pSrc0 + ( pTgt0 - pSrc0 ) * fact;
         gp_XYZ pLast1 = pSrc1 + ( pTgt1 - pSrc1 ) * fact;
         dist1 = dist2 = Precision::Infinite();
-        if ( !edge1->SegTriaInter( lastSeg, pSrc0, pTgt0, pSrc1, dist1, eps ) &&
-             !edge1->SegTriaInter( lastSeg, pSrc1, pTgt1, pTgt0, dist2, eps ))
-          continue;
-        if (( dist1 > testLen || dist1 < 0 ) &&
-            ( dist2 > testLen || dist2 < 0 ))
+        if ( !edge1->SegTriaInter( lastSeg, pSrc0, pLast0, pSrc1,  dist1, eps ) &&
+             !edge1->SegTriaInter( lastSeg, pSrc1, pLast1, pLast0, dist2, eps ))
           continue;
-
+        dist = dist1;
+        if ( dist > testLen || dist <= 0 )
+        {
+          dist = dist2;
+          if ( dist > testLen || dist <= 0 )
+            continue;
+        }
         // choose a closest edge
-        gp_Pnt intP( lastSeg.Location().XYZ() +
-                     lastSeg.Direction().XYZ() * ( Min( dist1, dist2 ) + segLen ));
+        gp_Pnt intP( lastSeg.Location().XYZ() + lastSeg.Direction().XYZ() * ( dist + segLen ));
         double d1 = intP.SquareDistance( pSrc0 );
         double d2 = intP.SquareDistance( pSrc1 );
         int iClose = i + ( d2 < d1 );
@@ -6245,15 +6940,14 @@ bool _ViscousBuilder::updateNormals( _SolidData&         data,
             ( d1 < d2 ? edgeJ : edge2 )->Set( _LayerEdge::MARKED );
           }
         }
-        intEdgesDist.push_back( make_pair( edge2, Min( dist1, dist2 )));
+        intEdgesDist.push_back( make_pair( edge2, dist ));
         // if ( Abs( d2 - d1 ) / Max( d2, d1 ) < 0.5 )
         // {
         //   iClose = i + !( d2 < d1 );
         //   intEdges.push_back( ce._intEdges[iClose] );
         //   ce._intEdges[iClose]->Unset( _LayerEdge::MARKED );
         // }
-        minIntDist = Min( edge1->_len * edge1->_lenFactor - segLen + dist1, minIntDist );
-        minIntDist = Min( edge1->_len * edge1->_lenFactor - segLen + dist2, minIntDist );
+        minIntDist = Min( edge1->_len * edge1->_lenFactor - segLen + dist, minIntDist );
       }
 
       //ce._edge = 0;
@@ -6261,8 +6955,10 @@ bool _ViscousBuilder::updateNormals( _SolidData&         data,
       // compute new _normals
       for ( size_t i = 0; i < intEdgesDist.size(); ++i )
       {
-        _LayerEdge* edge2    = intEdgesDist[i].first;
-        double       distWgt = edge1->_len / intEdgesDist[i].second;
+        _LayerEdge* edge2   = intEdgesDist[i].first;
+        double      distWgt = edge1->_len / intEdgesDist[i].second;
+        // if ( edge1->Is( _LayerEdge::BLOCKED ) &&
+        //      edge2->Is( _LayerEdge::BLOCKED )) continue;        
         if ( edge2->Is( _LayerEdge::MARKED )) continue;
         edge2->Set( _LayerEdge::MARKED );
 
@@ -6301,9 +6997,14 @@ bool _ViscousBuilder::updateNormals( _SolidData&         data,
         e2neIt->second._maxLen  = 0.7 * minIntDist / edge1->_lenFactor;
         if ( iter > 0 && sgn1 * sgn2 < 0 && edge1->_cosin < 0 )
           e2neIt->second._normal += dir2;
+
         e2neIt = edge2newEdge.insert( make_pair( edge2, zeroEdge )).first;
         e2neIt->second._normal += distWgt * newNormal;
-        e2neIt->second._cosin   = edge2->_cosin;
+        if ( Precision::IsInfinite( zeroEdge._maxLen ))
+        {
+          e2neIt->second._cosin  = edge2->_cosin;
+          e2neIt->second._maxLen = 1.3 * minIntDist / edge1->_lenFactor;
+        }
         if ( iter > 0 && sgn1 * sgn2 < 0 && edge2->_cosin < 0 )
           e2neIt->second._normal += dir1;
       }
@@ -6322,6 +7023,8 @@ bool _ViscousBuilder::updateNormals( _SolidData&         data,
       _LayerEdge*    edge = e2neIt->first;
       _LayerEdge& newEdge = e2neIt->second;
       _EdgesOnShape*  eos = data.GetShapeEdges( edge );
+      if ( edge->Is( _LayerEdge::BLOCKED && newEdge._maxLen > edge->_len ))
+        continue;
 
       // Check if a new _normal is OK:
       newEdge._normal.Normalize();
@@ -6350,11 +7053,9 @@ bool _ViscousBuilder::updateNormals( _SolidData&         data,
         else // edge inflates along a FACE
         {
           TopoDS_Shape V = helper.GetSubShapeByNode( edge->_nodes[0], getMeshDS() );
-          PShapeIteratorPtr eIt = helper.GetAncestors( V, *_mesh, TopAbs_EDGE );
+          PShapeIteratorPtr eIt = helper.GetAncestors( V, *_mesh, TopAbs_EDGE, &eos->_sWOL );
           while ( const TopoDS_Shape* E = eIt->next() )
           {
-            if ( !helper.IsSubShape( *E, /*FACE=*/eos->_sWOL ))
-              continue;
             gp_Vec edgeDir = getEdgeDir( TopoDS::Edge( *E ), TopoDS::Vertex( V ));
             double   angle = edgeDir.Angle( newEdge._normal ); // [0,PI]
             if ( angle < M_PI / 2 )
@@ -6390,7 +7091,7 @@ bool _ViscousBuilder::updateNormals( _SolidData&         data,
         edge1->SetDataByNeighbors( n1, n2, *eos1, helper );
       }
 
-      if ( !edge1->_2neibors )
+      if ( !edge1->_2neibors || !eos1->_sWOL.IsNull() )
         continue;
       for ( int j = 0; j < 2; ++j ) // loop on 2 neighbors
       {
@@ -6414,7 +7115,7 @@ bool _ViscousBuilder::updateNormals( _SolidData&         data,
             if ( nextEdge == prevEdge )
               nextEdge = neighbor->_2neibors->_edges[ ++iNext ];
           }
-          double r = double(step-1)/nbSteps;
+          double r = double(step-1)/nbSteps/(iter+1);
           if ( !nextEdge->_2neibors )
             r = Min( r, 0.5 );
 
@@ -6478,14 +7179,16 @@ bool _ViscousBuilder::isNewNormalOk( _SolidData&   data,
     newMinDot = Min( newNormal    * normFace, newMinDot );
     curMinDot = Min( edge._normal * normFace, curMinDot );
   }
+  bool ok = true;
   if ( newMinDot < 0.5 )
   {
-    return ( newMinDot >= curMinDot * 0.9 );
+    ok = ( newMinDot >= curMinDot * 0.9 );
     //return ( newMinDot >= ( curMinDot * ( 0.8 + 0.1 * edge.NbSteps() )));
     // double initMinDot2 = 1. - edge._cosin * edge._cosin;
     // return ( newMinDot * newMinDot ) >= ( 0.8 * initMinDot2 );
   }
-  return true;
+
+  return ok;
 }
 
 //================================================================================
@@ -6568,6 +7271,7 @@ void _ViscousBuilder::updateNormalsOfC1Vertices( _SolidData& data )
         oppV = SMESH_MesherHelper::IthVertex( 1, e );
       _EdgesOnShape* eovOpp = data.GetShapeEdges( oppV );
       if ( !eovOpp || eovOpp->_edges.empty() ) continue;
+      if ( eov._edges[0]->Is( _LayerEdge::BLOCKED )) continue;
 
       double curThickOpp = eovOpp->_edges[0]->_len * eovOpp->_edges[0]->_lenFactor;
       if ( curThickOpp + curThick < eLen )
@@ -6602,6 +7306,8 @@ bool _ViscousBuilder::updateNormalsOfConvexFaces( _SolidData&         data,
   for ( ; id2face != data._convexFaces.end(); ++id2face )
   {
     _ConvexFace & convFace = (*id2face).second;
+    convFace._normalsFixedOnBorders = false; // to update at each inflation step
+
     if ( convFace._normalsFixed )
       continue; // already fixed
     if ( convFace.CheckPrisms() )
@@ -6964,6 +7670,59 @@ bool _ViscousBuilder::updateNormalsOfConvexFaces( _SolidData&         data,
   return true;
 }
 
+//================================================================================
+/*!
+ * \brief Return max curvature of a FACE
+ */
+//================================================================================
+
+double _ConvexFace::GetMaxCurvature( _SolidData&         data,
+                                     _EdgesOnShape&      eof,
+                                     BRepLProp_SLProps&  surfProp,
+                                     SMESH_MesherHelper& helper)
+{
+  double maxCurvature = 0;
+
+  TopoDS_Face F = TopoDS::Face( eof._shape );
+
+  const int           nbTestPnt = 5;
+  const double        oriFactor = ( F.Orientation() == TopAbs_REVERSED ? +1. : -1. );
+  SMESH_subMeshIteratorPtr smIt = eof._subMesh->getDependsOnIterator(/*includeSelf=*/true);
+  while ( smIt->more() )
+  {
+    SMESH_subMesh* sm = smIt->next();
+    const TGeomID subID = sm->GetId();
+
+    // find _LayerEdge's of a sub-shape
+    _EdgesOnShape* eos;
+    if (( eos = data.GetShapeEdges( subID )))
+      this->_subIdToEOS.insert( make_pair( subID, eos ));
+    else
+      continue;
+
+    // check concavity and curvature and limit data._stepSize
+    const double minCurvature =
+      1. / ( eos->_hyp.GetTotalThickness() * ( 1 + theThickToIntersection ));
+    size_t iStep = Max( 1, eos->_edges.size() / nbTestPnt );
+    for ( size_t i = 0; i < eos->_edges.size(); i += iStep )
+    {
+      gp_XY uv = helper.GetNodeUV( F, eos->_edges[ i ]->_nodes[0] );
+      surfProp.SetParameters( uv.X(), uv.Y() );
+      if ( surfProp.IsCurvatureDefined() )
+      {
+        double curvature = Max( surfProp.MaxCurvature() * oriFactor,
+                                surfProp.MinCurvature() * oriFactor );
+        maxCurvature = Max( maxCurvature, curvature );
+
+        if ( curvature > minCurvature )
+          this->_isTooCurved = true;
+      }
+    }
+  } // loop on sub-shapes of the FACE
+
+  return maxCurvature;
+}
+
 //================================================================================
 /*!
  * \brief Finds a center of curvature of a surface at a _LayerEdge
@@ -7164,6 +7923,28 @@ bool _LayerEdge::FindIntersection( SMESH_ElementSearcher&   searcher,
   return segmentIntersected;
 }
 
+//================================================================================
+/*!
+ * \brief Returns a point used to check orientation of _simplices
+ */
+//================================================================================
+
+gp_XYZ _LayerEdge::PrevCheckPos( _EdgesOnShape* eos ) const
+{
+  size_t i = Is( NORMAL_UPDATED ) && IsOnFace() ? _pos.size()-2 : 0;
+
+  if ( !eos || eos->_sWOL.IsNull() )
+    return _pos[ i ];
+
+  if ( eos->SWOLType() == TopAbs_EDGE )
+  {
+    return BRepAdaptor_Curve( TopoDS::Edge( eos->_sWOL )).Value( _pos[i].X() ).XYZ();
+  }
+  //else //  TopAbs_FACE
+
+  return BRepAdaptor_Surface( TopoDS::Face( eos->_sWOL )).Value(_pos[i].X(), _pos[i].Y() ).XYZ();
+}
+
 //================================================================================
 /*!
  * \brief Returns size and direction of the last segment
@@ -7224,13 +8005,14 @@ gp_Ax1 _LayerEdge::LastSegment(double& segLen, _EdgesOnShape& eos) const
 
 //================================================================================
 /*!
- * \brief Return the last position of the target node on a FACE. 
+ * \brief Return the last (or \a which) position of the target node on a FACE. 
  *  \param [in] F - the FACE this _LayerEdge is inflated along
+ *  \param [in] which - index of position
  *  \return gp_XY - result UV
  */
 //================================================================================
 
-gp_XY _LayerEdge::LastUV( const TopoDS_Face& F, _EdgesOnShape& eos ) const
+gp_XY _LayerEdge::LastUV( const TopoDS_Face& F, _EdgesOnShape& eos, int which ) const
 {
   if ( F.IsSame( eos._sWOL )) // F is my FACE
     return gp_XY( _pos.back().X(), _pos.back().Y() );
@@ -7239,7 +8021,7 @@ gp_XY _LayerEdge::LastUV( const TopoDS_Face& F, _EdgesOnShape& eos ) const
     return gp_XY( 1e100, 1e100 );
 
   // _sWOL is EDGE of F; _pos.back().X() is the last U on the EDGE
-  double f, l, u = _pos.back().X();
+  double f, l, u = _pos[ which < 0 ? _pos.size()-1 : which ].X();
   Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface( TopoDS::Edge(eos._sWOL), F, f,l);
   if ( !C2d.IsNull() && f <= u && u <= l )
     return C2d->Value( u ).XY();
@@ -7266,7 +8048,7 @@ bool _LayerEdge::SegTriaInter( const gp_Ax1& lastSegment,
   const gp_Dir& dir  = lastSegment.Direction();
 
   /* calculate distance from vert0 to ray origin */
-  gp_XYZ tvec = orig.XYZ() - vert0;
+  //gp_XYZ tvec = orig.XYZ() - vert0;
 
   //if ( tvec * dir > EPSILON )
     // intersected face is at back side of the temporary face this _LayerEdge belongs to
@@ -7285,6 +8067,9 @@ bool _LayerEdge::SegTriaInter( const gp_Ax1& lastSegment,
   if ( det > -ANGL_EPSILON && det < ANGL_EPSILON )
     return false;
 
+  /* calculate distance from vert0 to ray origin */
+  gp_XYZ tvec = orig.XYZ() - vert0;
+
   /* calculate U parameter and test bounds */
   double u = ( tvec * pvec ) / det;
   //if (u < 0.0 || u > 1.0)
@@ -7314,7 +8099,7 @@ bool _LayerEdge::SegTriaInter( const gp_Ax1& lastSegment,
  *  \param [in] eov - EOS of the VERTEX
  *  \param [in] eos - EOS of the FACE
  *  \param [in] step - inflation step
- *  \param [in,out] badSmooEdges - not untangled _LayerEdge's
+ *  \param [in,out] badSmooEdges - tangled _LayerEdge's
  */
 //================================================================================
 
@@ -7573,7 +8358,9 @@ int _LayerEdge::CheckNeiborsOnBoundary( vector< _LayerEdge* >* badNeibors, bool
     if ( eN->_nodes[0]->getshapeId() == _nodes[0]->getshapeId() )
       continue;
     if ( needSmooth )
-      *needSmooth |= ( eN->Is( _LayerEdge::BLOCKED ) || eN->Is( _LayerEdge::NORMAL_UPDATED ));
+      *needSmooth |= ( eN->Is( _LayerEdge::BLOCKED ) ||
+                       eN->Is( _LayerEdge::NORMAL_UPDATED ) ||
+                       eN->_pos.size() != _pos.size() );
 
     SMESH_TNodeXYZ curPosN ( eN->_nodes.back() );
     SMESH_TNodeXYZ prevPosN( eN->_nodes[0] );
@@ -7619,7 +8406,7 @@ int _LayerEdge::Smooth(const int step, bool findBest, vector< _LayerEdge* >& toS
     findBest = true;
 
   const gp_XYZ& curPos  = _pos.back();
-  const gp_XYZ& prevPos = PrevCheckPos();
+  const gp_XYZ& prevPos = _pos[0]; //PrevPos();
 
   // quality metrics (orientation) of tetras around _tgtNode
   int nbOkBefore = 0;
@@ -7763,7 +8550,7 @@ int _LayerEdge::Smooth(const int step, const bool isConcaveFace, bool findBest )
     return 0; // not inflated
 
   const gp_XYZ& curPos  = _pos.back();
-  const gp_XYZ& prevPos = PrevCheckPos();
+  const gp_XYZ& prevPos = _pos[0]; //PrevCheckPos();
 
   // quality metrics (orientation) of tetras around _tgtNode
   int nbOkBefore = 0;
@@ -8148,7 +8935,7 @@ gp_XYZ _LayerEdge::smoothNefPolygon()
 {
   gp_XYZ newPos(0,0,0);
 
-  // get a plane to seach a solution on
+  // get a plane to search a solution on
 
   vector< gp_XYZ > vecs( _simplices.size() + 1 );
   size_t i;
@@ -8339,7 +9126,7 @@ gp_XYZ _LayerEdge::smoothNefPolygon()
 { ////////////////////////////////// NEW
   gp_XYZ newPos(0,0,0);
 
-  // get a plane to seach a solution on
+  // get a plane to search a solution on
 
   size_t i;
   gp_XYZ center(0,0,0);
@@ -8580,13 +9367,14 @@ void _LayerEdge::SetNewLength( double len, _EdgesOnShape& eos, SMESH_MesherHelpe
 
       // find point of intersection of the face plane located at baryCenter
       // and _normal located at newXYZ
-      double d    = -( faceNorm.XYZ() * baryCenter ); // d of plane equation ax+by+cz+d=0
-      double dot  = ( faceNorm.XYZ() * _normal );
+      double d   = -( faceNorm.XYZ() * baryCenter ); // d of plane equation ax+by+cz+d=0
+      double dot  ( faceNorm.XYZ() * _normal );
       if ( dot < std::numeric_limits<double>::min() )
         dot = lenDelta * 1e-3;
       double step = -( faceNorm.XYZ() * newXYZ + d ) / dot;
       newXYZ += step * _normal;
     }
+    _lenFactor = _normal * ( newXYZ - oldXYZ ) / lenDelta; // _lenFactor is used in InvalidateStep()
   }
   else
   {
@@ -8660,9 +9448,14 @@ void _LayerEdge::SetNewLength( double len, _EdgesOnShape& eos, SMESH_MesherHelpe
 
 void _LayerEdge::Block( _SolidData& data )
 {
-  if ( Is( BLOCKED )) return;
+  //if ( Is( BLOCKED )) return;
   Set( BLOCKED );
 
+  SMESH_Comment msg( "#BLOCK shape=");
+  msg << data.GetShapeEdges( this )->_shapeID
+      << ", nodes " << _nodes[0]->GetID() << ", " << _nodes.back()->GetID();
+  dumpCmd( msg + " -- BEGIN")
+
   _maxLen = _len;
   std::queue<_LayerEdge*> queue;
   queue.push( this );
@@ -8676,8 +9469,7 @@ void _LayerEdge::Block( _SolidData& data )
     for ( size_t iN = 0; iN < edge->_neibors.size(); ++iN )
     {
       _LayerEdge* neibor = edge->_neibors[iN];
-      if ( neibor->Is( BLOCKED ) ||
-           neibor->_maxLen < edge->_maxLen )
+      if ( neibor->_maxLen < edge->_maxLen * 1.01 )
         continue;
       pSrcN = SMESH_TNodeXYZ( neibor->_nodes[0] );
       pTgtN = SMESH_TNodeXYZ( neibor->_nodes.back() );
@@ -8686,9 +9478,11 @@ void _LayerEdge::Block( _SolidData& data )
       minDist   = Min( pSrc.SquareDistance( pTgtN ), minDist );
       minDist   = Min( pTgt.SquareDistance( pSrcN ), minDist );
       double newMaxLen = edge->_maxLen + 0.5 * Sqrt( minDist );
-      if ( edge->_nodes[0]->getshapeId() == neibor->_nodes[0]->getshapeId() )
+      //if ( edge->_nodes[0]->getshapeId() == neibor->_nodes[0]->getshapeId() ) viscous_layers_00/A3
       {
         newMaxLen *= edge->_lenFactor / neibor->_lenFactor;
+        // newMaxLen *= Min( edge->_lenFactor / neibor->_lenFactor,
+        //                   neibor->_lenFactor / edge->_lenFactor );
       }
       if ( neibor->_maxLen > newMaxLen )
       {
@@ -8697,14 +9491,16 @@ void _LayerEdge::Block( _SolidData& data )
         {
           _EdgesOnShape* eos = data.GetShapeEdges( neibor );
           while ( neibor->_len > neibor->_maxLen &&
-                  neibor->NbSteps() > 1 )
+                  neibor->NbSteps() > 0 )
             neibor->InvalidateStep( neibor->NbSteps(), *eos, /*restoreLength=*/true );
           neibor->SetNewLength( neibor->_maxLen, *eos, data.GetHelper() );
+          //neibor->Block( data );
         }
         queue.push( neibor );
       }
     }
   }
+  dumpCmd( msg + " -- END")
 }
 
 //================================================================================
@@ -8747,24 +9543,22 @@ void _LayerEdge::InvalidateStep( size_t curStep, const _EdgesOnShape& eos, bool
 
     if ( restoreLength )
     {
-      _len -= ( nXYZ.XYZ() - curXYZ ).Modulus() / _lenFactor;
+      if ( NbSteps() == 0 )
+        _len = 0.;
+      else
+        _len -= ( nXYZ.XYZ() - curXYZ ).Modulus() / _lenFactor;
     }
   }
 }
 
 //================================================================================
 /*!
- * \brief Smooth a path formed by _pos of a _LayerEdge smoothed on FACE
+ * \brief Return index of a _pos distant from _normal
  */
 //================================================================================
 
-void _LayerEdge::SmoothPos( const vector< double >& segLen, const double tol )
+int _LayerEdge::GetSmoothedPos( const double tol )
 {
-  //return;
-  if ( /*Is( NORMAL_UPDATED ) ||*/ _pos.size() <= 2 )
-    return;
-
-  // find the 1st smoothed _pos
   int iSmoothed = 0;
   for ( size_t i = 1; i < _pos.size() && !iSmoothed; ++i )
   {
@@ -8772,12 +9566,26 @@ void _LayerEdge::SmoothPos( const vector< double >& segLen, const double tol )
     if ( normDist > tol * tol )
       iSmoothed = i;
   }
+  return iSmoothed;
+}
+
+//================================================================================
+/*!
+ * \brief Smooth a path formed by _pos of a _LayerEdge smoothed on FACE
+ */
+//================================================================================
+
+void _LayerEdge::SmoothPos( const vector< double >& segLen, const double tol )
+{
+  if ( /*Is( NORMAL_UPDATED ) ||*/ _pos.size() <= 2 )
+    return;
+
+  // find the 1st smoothed _pos
+  int iSmoothed = GetSmoothedPos( tol );
   if ( !iSmoothed ) return;
 
-  if ( 1 || Is( DISTORTED ))
+  //if ( 1 || Is( DISTORTED ))
   {
-    // if ( segLen[ iSmoothed ] / segLen.back() < 0.5 )
-    //   return;
     gp_XYZ normal = _normal;
     if ( Is( NORMAL_UPDATED ))
       for ( size_t i = 1; i < _pos.size(); ++i )
@@ -8791,7 +9599,7 @@ void _LayerEdge::SmoothPos( const vector< double >& segLen, const double tol )
         }
       }
     const double r = 0.2;
-    for ( int iter = 0; iter < 3; ++iter )
+    for ( int iter = 0; iter < 50; ++iter )
     {
       double minDot = 1;
       for ( size_t i = Max( 1, iSmoothed-1-iter ); i < _pos.size()-1; ++i )
@@ -8804,11 +9612,11 @@ void _LayerEdge::SmoothPos( const vector< double >& segLen, const double tol )
         const_cast< double& >( segLen[i] ) = newLen;
         // check angle between normal and (_pos[i+1], _pos[i] )
         gp_XYZ posDir = _pos[i+1] - _pos[i];
-        double size   = posDir.Modulus();
+        double size   = posDir.SquareModulus();
         if ( size > RealSmall() )
-          minDot = Min( minDot, ( normal * posDir ) / size );
+          minDot = Min( minDot, ( normal * posDir ) * ( normal * posDir ) / size );
       }
-      if ( minDot > 0.5 )
+      if ( minDot > 0.5 * 0.5 )
         break;
     }
   }
@@ -8830,10 +9638,49 @@ void _LayerEdge::SmoothPos( const vector< double >& segLen, const double tol )
 
 //================================================================================
 /*!
- * \brief Create layers of prisms
+ * \brief Print flags
  */
 //================================================================================
 
+std::string _LayerEdge::DumpFlags() const
+{
+  SMESH_Comment dump;
+  for ( int flag = 1; flag < 0x1000000; flag *= 2 )
+    if ( _flags & flag )
+    {
+      EFlags f = (EFlags) flag;
+      switch ( f ) {
+      case TO_SMOOTH:       dump << "TO_SMOOTH";       break;
+      case MOVED:           dump << "MOVED";           break;
+      case SMOOTHED:        dump << "SMOOTHED";        break;
+      case DIFFICULT:       dump << "DIFFICULT";       break;
+      case ON_CONCAVE_FACE: dump << "ON_CONCAVE_FACE"; break;
+      case BLOCKED:         dump << "BLOCKED";         break;
+      case INTERSECTED:     dump << "INTERSECTED";     break;
+      case NORMAL_UPDATED:  dump << "NORMAL_UPDATED";  break;
+      case UPD_NORMAL_CONV: dump << "UPD_NORMAL_CONV"; break;
+      case MARKED:          dump << "MARKED";          break;
+      case MULTI_NORMAL:    dump << "MULTI_NORMAL";    break;
+      case NEAR_BOUNDARY:   dump << "NEAR_BOUNDARY";   break;
+      case SMOOTHED_C1:     dump << "SMOOTHED_C1";     break;
+      case DISTORTED:       dump << "DISTORTED";       break;
+      case RISKY_SWOL:      dump << "RISKY_SWOL";      break;
+      case SHRUNK:          dump << "SHRUNK";          break;
+      case UNUSED_FLAG:     dump << "UNUSED_FLAG";     break;
+      }
+      dump << " ";
+    }
+  cout << dump << endl;
+  return dump;
+}
+
+//================================================================================
+/*!
+  case brief:
+  default:
+*/
+//================================================================================
+
 bool _ViscousBuilder::refine(_SolidData& data)
 {
   SMESH_MesherHelper& helper = data.GetHelper();
@@ -8844,10 +9691,10 @@ bool _ViscousBuilder::refine(_SolidData& data)
   TopoDS_Edge geomEdge;
   TopoDS_Face geomFace;
   TopLoc_Location loc;
-  double f,l, u;
+  double f,l, u = 0;
   gp_XY uv;
   vector< gp_XYZ > pos3D;
-  bool isOnEdge;
+  bool isOnEdge, isTooConvexFace = false;
   TGeomID prevBaseId = -1;
   TNode2Edge* n2eMap = 0;
   TNode2Edge::iterator n2e;
@@ -8891,6 +9738,9 @@ bool _ViscousBuilder::refine(_SolidData& data)
         for ( size_t j = 0; j < eos._eosC1[i]->_edges.size(); ++j )
           eos._eosC1[i]->_edges[j]->Set( _LayerEdge::SMOOTHED_C1 );
       }
+      isTooConvexFace = false;
+      if ( _ConvexFace* cf = data.GetConvexFace( eos._shapeID ))
+        isTooConvexFace = cf->_isTooCurved;
     }
 
     vector< double > segLen;
@@ -8906,18 +9756,16 @@ bool _ViscousBuilder::refine(_SolidData& data)
       if ( eos._sWOL.IsNull() )
       {
         bool useNormal = true;
-        bool   usePos  = false;
-        bool smoothed  = false;
-        const double preci = 0.1 * edge._len;
-        if ( eos._toSmooth )
+        bool    usePos = false;
+        bool  smoothed = false;
+        double   preci = 0.1 * edge._len;
+        if ( eos._toSmooth && edge._pos.size() > 2 )
         {
-          gp_Pnt tgtExpected = edge._pos[0] + edge._normal * edge._len;
-          smoothed = tgtExpected.SquareDistance( edge._pos.back() ) > preci * preci;
+          smoothed = edge.GetSmoothedPos( preci );
         }
         if ( smoothed )
         {
-          if ( !surface.IsNull() &&
-               !data._convexFaces.count( eos._shapeID )) // edge smoothed on FACE
+          if ( !surface.IsNull() && !isTooConvexFace ) // edge smoothed on FACE
           {
             useNormal = usePos = false;
             gp_Pnt2d uv = helper.GetNodeUV( geomFace, edge._nodes[0] );
@@ -8931,13 +9779,15 @@ bool _ViscousBuilder::refine(_SolidData& data)
             }
           }
         }
-        else
+        else if ( !edge.Is( _LayerEdge::NORMAL_UPDATED ))
         {
+#ifndef __NODES_AT_POS
           useNormal = usePos = false;
           edge._pos[1] = edge._pos.back();
           edge._pos.resize( 2 );
           segLen.resize( 2 );
           segLen[ 1 ] = edge._len;
+#endif
         }
         if ( useNormal && edge.Is( _LayerEdge::NORMAL_UPDATED ))
         {
@@ -8966,11 +9816,12 @@ bool _ViscousBuilder::refine(_SolidData& data)
           while ( swapped )
           {
             swapped = false;
-            for ( size_t j = 1; j < edge._pos.size(); ++j )
+            for ( size_t j = 1; j < edge._pos.size()-1; ++j )
               if ( segLen[j] > segLen.back() )
               {
                 segLen.erase( segLen.begin() + j );
                 edge._pos.erase( edge._pos.begin() + j );
+                --j;
               }
               else if ( segLen[j] < segLen[j-1] )
               {
@@ -8981,14 +9832,28 @@ bool _ViscousBuilder::refine(_SolidData& data)
           }
         }
         // smooth a path formed by edge._pos
+#ifndef __NODES_AT_POS
         if (( smoothed ) /*&&
             ( eos.ShapeType() == TopAbs_FACE || edge.Is( _LayerEdge::SMOOTHED_C1 ))*/)
           edge.SmoothPos( segLen, preci );
+#endif
       }
       else if ( eos._isRegularSWOL ) // usual SWOL
       {
-        for ( size_t j = 1; j < edge._pos.size(); ++j )
-          segLen[j] = segLen[j-1] + (edge._pos[j-1] - edge._pos[j] ).Modulus();
+        if ( edge.Is( _LayerEdge::SMOOTHED ))
+        {
+          SMESH_NodeXYZ p0( edge._nodes[0] );
+          for ( size_t j = 1; j < edge._pos.size(); ++j )
+          {
+            gp_XYZ pj = surface->Value( edge._pos[j].X(), edge._pos[j].Y() ).XYZ();
+            segLen[j] = ( pj - p0 ) * edge._normal;
+          }
+        }
+        else
+        {
+          for ( size_t j = 1; j < edge._pos.size(); ++j )
+            segLen[j] = segLen[j-1] + (edge._pos[j-1] - edge._pos[j] ).Modulus();
+        }
       }
       else if ( !surface.IsNull() ) // SWOL surface with singularities
       {
@@ -9004,7 +9869,12 @@ bool _ViscousBuilder::refine(_SolidData& data)
       const SMDS_MeshNode* tgtNode = edge._nodes.back();
       if ( edge._nodes.size() == 2 )
       {
-        edge._nodes.resize( eos._hyp.GetNumberLayers() + 1, 0 );
+#ifdef __NODES_AT_POS
+        int nbNodes = edge._pos.size();
+#else
+        int nbNodes = eos._hyp.GetNumberLayers() + 1;
+#endif
+        edge._nodes.resize( nbNodes, 0 );
         edge._nodes[1] = 0;
         edge._nodes.back() = tgtNode;
       }
@@ -9064,7 +9934,9 @@ bool _ViscousBuilder::refine(_SolidData& data)
           --iPrevSeg;
         double   r = ( segLen[iSeg] - hSum ) / ( segLen[iSeg] - segLen[iPrevSeg] );
         gp_Pnt pos = r * edge._pos[iPrevSeg] + (1-r) * edge._pos[iSeg];
-
+#ifdef __NODES_AT_POS
+        pos = edge._pos[ iStep ];
+#endif
         SMDS_MeshNode*& node = const_cast< SMDS_MeshNode*& >( edge._nodes[ iStep ]);
         if ( !eos._sWOL.IsNull() )
         {
@@ -9169,7 +10041,6 @@ bool _ViscousBuilder::refine(_SolidData& data)
   set< vector<const SMDS_MeshNode*>* >    nnSet;
   set< int >                       degenEdgeInd;
   vector<const SMDS_MeshElement*>     degenVols;
-  vector<int>                       isRiskySWOL;
 
   TopExp_Explorer exp( data._solid, TopAbs_FACE );
   for ( ; exp.More(); exp.Next() )
@@ -9187,7 +10058,6 @@ bool _ViscousBuilder::refine(_SolidData& data)
       nnVec.resize( nbNodes );
       nnSet.clear();
       degenEdgeInd.clear();
-      isRiskySWOL.resize( nbNodes );
       size_t maxZ = 0, minZ = std::numeric_limits<size_t>::max();
       SMDS_NodeIteratorPtr nIt = face->nodeIterator();
       for ( int iN = 0; iN < nbNodes; ++iN )
@@ -9198,7 +10068,6 @@ bool _ViscousBuilder::refine(_SolidData& data)
         nnVec[ i ] = & edge->_nodes;
         maxZ = std::max( maxZ, nnVec[ i ]->size() );
         minZ = std::min( minZ, nnVec[ i ]->size() );
-        //isRiskySWOL[ i ] = edge->Is( _LayerEdge::RISKY_SWOL );
 
         if ( helper.HasDegeneratedEdges() )
           nnSet.insert( nnVec[ i ]);
@@ -9307,7 +10176,7 @@ bool _ViscousBuilder::refine(_SolidData& data)
     if ( !err || err->IsOK() )
     {
       err.reset( new SMESH_ComputeError( COMPERR_WARNING,
-                                         "Degenerated volumes created" ));
+                                         "Bad quality volumes created" ));
       err->myBadElements.insert( err->myBadElements.end(),
                                  degenVols.begin(),degenVols.end() );
     }
@@ -9322,71 +10191,94 @@ bool _ViscousBuilder::refine(_SolidData& data)
  */
 //================================================================================
 
-bool _ViscousBuilder::shrink()
+bool _ViscousBuilder::shrink(_SolidData& theData)
 {
-  // make map of (ids of FACEs to shrink mesh on) to (_SolidData containing _LayerEdge's
-  // inflated along FACE or EDGE)
-  map< TGeomID, _SolidData* > f2sdMap;
+  // make map of (ids of FACEs to shrink mesh on) to (list of _SolidData containing
+  // _LayerEdge's inflated along FACE or EDGE)
+  map< TGeomID, list< _SolidData* > > f2sdMap;
   for ( size_t i = 0 ; i < _sdVec.size(); ++i )
   {
     _SolidData& data = _sdVec[i];
-    TopTools_MapOfShape FFMap;
     map< TGeomID, TopoDS_Shape >::iterator s2s = data._shrinkShape2Shape.begin();
     for (; s2s != data._shrinkShape2Shape.end(); ++s2s )
-      if ( s2s->second.ShapeType() == TopAbs_FACE )
+      if ( s2s->second.ShapeType() == TopAbs_FACE && !_shrinkedFaces.Contains( s2s->second ))
       {
-        f2sdMap.insert( make_pair( getMeshDS()->ShapeToIndex( s2s->second ), &data ));
+        f2sdMap[ getMeshDS()->ShapeToIndex( s2s->second )].push_back( &data );
 
-        if ( FFMap.Add( (*s2s).second ))
-          // Put mesh faces on the shrinked FACE to the proxy sub-mesh to avoid
-          // usage of mesh faces made in addBoundaryElements() by the 3D algo or
-          // by StdMeshers_QuadToTriaAdaptor
-          if ( SMESHDS_SubMesh* smDS = getMeshDS()->MeshElements( s2s->second ))
+        // Put mesh faces on the shrinked FACE to the proxy sub-mesh to avoid
+        // usage of mesh faces made in addBoundaryElements() by the 3D algo or
+        // by StdMeshers_QuadToTriaAdaptor
+        if ( SMESHDS_SubMesh* smDS = getMeshDS()->MeshElements( s2s->second ))
+        {
+          SMESH_ProxyMesh::SubMesh* proxySub =
+            data._proxyMesh->getFaceSubM( TopoDS::Face( s2s->second ), /*create=*/true);
+          if ( proxySub->NbElements() == 0 )
           {
-            SMESH_ProxyMesh::SubMesh* proxySub =
-              data._proxyMesh->getFaceSubM( TopoDS::Face( s2s->second ), /*create=*/true);
             SMDS_ElemIteratorPtr fIt = smDS->GetElements();
             while ( fIt->more() )
-              proxySub->AddElement( fIt->next() );
-            // as a result 3D algo will use elements from proxySub and not from smDS
+            {
+              const SMDS_MeshElement* f = fIt->next();
+              // as a result 3D algo will use elements from proxySub and not from smDS
+              proxySub->AddElement( f );
+              f->setIsMarked( true );
+
+              // Mark nodes on the FACE to discriminate them from nodes
+              // added by addBoundaryElements(); marked nodes are to be smoothed while shrink()
+              for ( int iN = 0, nbN = f->NbNodes(); iN < nbN; ++iN )
+              {
+                const SMDS_MeshNode* n = f->GetNode( iN );
+                if ( n->GetPosition()->GetDim() == 2 )
+                  n->setIsMarked( true );
+              }
+            }
           }
+        }
       }
   }
 
   SMESH_MesherHelper helper( *_mesh );
   helper.ToFixNodeParameters( true );
 
-  // EDGE's to shrink
+  // EDGEs to shrink
   map< TGeomID, _Shrinker1D > e2shrMap;
   vector< _EdgesOnShape* > subEOS;
   vector< _LayerEdge* > lEdges;
 
-  // loop on FACES to srink mesh on
-  map< TGeomID, _SolidData* >::iterator f2sd = f2sdMap.begin();
+  // loop on FACEs to srink mesh on
+  map< TGeomID, list< _SolidData* > >::iterator f2sd = f2sdMap.begin();
   for ( ; f2sd != f2sdMap.end(); ++f2sd )
   {
-    _SolidData&      data = *f2sd->second;
+    list< _SolidData* > & dataList = f2sd->second;
+    if ( dataList.front()->_n2eMap.empty() ||
+         dataList.back() ->_n2eMap.empty() )
+      continue; // not yet computed
+    if ( dataList.front() != &theData &&
+         dataList.back()  != &theData )
+      continue;
+
+    _SolidData&      data = *dataList.front();
+    _SolidData*     data2 = dataList.size() > 1 ? dataList.back() : 0;
     const TopoDS_Face&  F = TopoDS::Face( getMeshDS()->IndexToShape( f2sd->first ));
     SMESH_subMesh*     sm = _mesh->GetSubMesh( F );
     SMESHDS_SubMesh* smDS = sm->GetSubMeshDS();
 
-    Handle(Geom_Surface) surface = BRep_Tool::Surface(F);
+    Handle(Geom_Surface) surface = BRep_Tool::Surface( F );
 
-    helper.SetSubShape(F);
+    _shrinkedFaces.Add( F );
+    helper.SetSubShape( F );
 
     // ===========================
     // Prepare data for shrinking
     // ===========================
 
-    // Collect nodes to smooth, as src nodes are not yet replaced by tgt ones
-    // and hence all nodes on a FACE connected to 2d elements are to be smoothed
+    // Collect nodes to smooth (they are marked at the beginning of this method)
     vector < const SMDS_MeshNode* > smoothNodes;
     {
       SMDS_NodeIteratorPtr nIt = smDS->GetNodes();
       while ( nIt->more() )
       {
         const SMDS_MeshNode* n = nIt->next();
-        if ( n->NbInverseElements( SMDSAbs_Face ) > 0 )
+        if ( n->isMarked() )
           smoothNodes.push_back( n );
       }
     }
@@ -9398,10 +10290,10 @@ bool _ViscousBuilder::shrink()
     {
       vector<_Simplex> simplices;
       _Simplex::GetSimplices( smoothNodes[0], simplices, ignoreShapes );
-      helper.GetNodeUV( F, simplices[0]._nPrev, 0, &isOkUV ); // fix UV of silpmex nodes
+      helper.GetNodeUV( F, simplices[0]._nPrev, 0, &isOkUV ); // fix UV of simplex nodes
       helper.GetNodeUV( F, simplices[0]._nNext, 0, &isOkUV );
       gp_XY uv = helper.GetNodeUV( F, smoothNodes[0], 0, &isOkUV );
-      if ( !simplices[0].IsForward(uv, smoothNodes[0], F, helper,refSign) )
+      if ( !simplices[0].IsForward(uv, smoothNodes[0], F, helper, refSign ))
         refSign = -1;
     }
 
@@ -9417,8 +10309,13 @@ bool _ViscousBuilder::shrink()
         if ( data._noShrinkShapes.count( subID ))
           continue;
         _EdgesOnShape* eos = data.GetShapeEdges( subID );
-        if ( !eos || eos->_sWOL.IsNull() ) continue;
-
+        if ( !eos || eos->_sWOL.IsNull() )
+          if ( data2 ) // check in adjacent SOLID
+          {
+            eos = data2->GetShapeEdges( subID );
+            if ( !eos || eos->_sWOL.IsNull() )
+              continue;
+          }
         subEOS.push_back( eos );
 
         for ( size_t i = 0; i < eos->_edges.size(); ++i )
@@ -9451,7 +10348,7 @@ bool _ViscousBuilder::shrink()
         while ( fIt->more() )
         {
           const SMDS_MeshElement* f = fIt->next();
-          if ( !smDS->Contains( f ))
+          if ( !smDS->Contains( f ) || !f->isMarked() )
             continue;
           SMDS_NodeIteratorPtr nIt = f->nodeIterator();
           for ( int iN = 0; nIt->more(); ++iN )
@@ -9478,7 +10375,7 @@ bool _ViscousBuilder::shrink()
       {
         const SMDS_MeshNode* n = smoothNodes[i];
         nodesToSmooth[ i ]._node = n;
-        // src nodes must be replaced by tgt nodes to have tgt nodes in _simplices
+        // src nodes must be already replaced by tgt nodes to have tgt nodes in _simplices
         _Simplex::GetSimplices( n, nodesToSmooth[ i ]._simplices, ignoreShapes, 0, sortSimplices);
         // fix up incorrect uv of nodes on the FACE
         helper.GetNodeUV( F, n, 0, &isOkUV);
@@ -9509,6 +10406,10 @@ bool _ViscousBuilder::shrink()
         {
           _LayerEdge& edge = * eos._edges[i];
           _Simplex::GetSimplices( /*tgtNode=*/edge._nodes.back(), edge._simplices, ignoreShapes );
+
+          // additionally mark tgt node; only marked nodes will be used in SetNewLength2d()
+          // not-marked nodes are those added by refine()
+          edge._nodes.back()->setIsMarked( true );
         }
       }
     }
@@ -9645,6 +10546,8 @@ bool _ViscousBuilder::shrink()
 
     if ( !errMsg.empty() ) // Try to re-compute the shrink FACE
     {
+      debugMsg( "Re-compute FACE " << f2sd->first << " because " << errMsg );
+
       // remove faces
       SMESHDS_SubMesh* psm = data._proxyMesh->getFaceSubM( F );
       {
@@ -9684,7 +10587,8 @@ bool _ViscousBuilder::shrink()
         {
           _LayerEdge*       edge = subEOS[iS]->_edges[i];
           SMDS_MeshNode* tgtNode = const_cast< SMDS_MeshNode*& >( edge->_nodes.back() );
-          if ( edge->_pos.empty() ) continue;
+          if ( edge->_pos.empty() ||
+               edge->Is( _LayerEdge::SHRUNK )) continue;
           if ( subEOS[iS]->SWOLType() == TopAbs_FACE )
           {
             SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( tgtNode->GetPosition() );
@@ -9840,6 +10744,8 @@ bool _ViscousBuilder::shrink()
 
     // Set an event listener to clear FACE sub-mesh together with SOLID sub-mesh
     VISCOUS_3D::ToClearSubWithMain( sm, data._solid );
+    if ( data2 )
+      VISCOUS_3D::ToClearSubWithMain( sm, data2->_solid );
 
   } // loop on FACES to srink mesh on
 
@@ -9872,6 +10778,7 @@ bool _ViscousBuilder::prepareEdgeToShrink( _LayerEdge&            edge,
     if ( tgtNode->GetPosition()->GetDim() != 2 ) // not inflated edge
     {
       edge._pos.clear();
+      edge.Set( _LayerEdge::SHRUNK );
       return srcNode == tgtNode;
     }
     gp_XY srcUV ( edge._pos[0].X(), edge._pos[0].Y() );          //helper.GetNodeUV( F, srcNode );
@@ -9882,7 +10789,7 @@ bool _ViscousBuilder::prepareEdgeToShrink( _LayerEdge&            edge,
     edge._normal.SetCoord( uvDir.X(),uvDir.Y(), 0 );
     edge._len = uvLen;
 
-    edge._pos.resize(1);
+    //edge._pos.resize(1);
     edge._pos[0].SetCoord( tgtUV.X(), tgtUV.Y(), 0 );
 
     // set UV of source node to target node
@@ -9895,6 +10802,7 @@ bool _ViscousBuilder::prepareEdgeToShrink( _LayerEdge&            edge,
     if ( tgtNode->GetPosition()->GetDim() != 1 ) // not inflated edge
     {
       edge._pos.clear();
+      edge.Set( _LayerEdge::SHRUNK );
       return srcNode == tgtNode;
     }
     const TopoDS_Edge&    E = TopoDS::Edge( eos._sWOL );
@@ -9916,16 +10824,17 @@ bool _ViscousBuilder::prepareEdgeToShrink( _LayerEdge&            edge,
 
     double uSrc = helper.GetNodeU( E, srcNode, n2 );
     double uTgt = helper.GetNodeU( E, tgtNode, srcNode );
-    double u2   = helper.GetNodeU( E, n2, srcNode );
+    double u2   = helper.GetNodeU( E, n2,      srcNode );
 
-    edge._pos.clear();
+    //edge._pos.clear();
 
     if ( fabs( uSrc-uTgt ) < 0.99 * fabs( uSrc-u2 ))
     {
       // tgtNode is located so that it does not make faces with wrong orientation
+      edge.Set( _LayerEdge::SHRUNK );
       return true;
     }
-    edge._pos.resize(1);
+    //edge._pos.resize(1);
     edge._pos[0].SetCoord( U_TGT, uTgt );
     edge._pos[0].SetCoord( U_SRC, uSrc );
     edge._pos[0].SetCoord( LEN_TGT, fabs( uSrc-uTgt ));
@@ -10146,7 +11055,7 @@ bool _LayerEdge::SetNewLength2d( Handle(Geom_Surface)& surface,
                                  _EdgesOnShape&        eos,
                                  SMESH_MesherHelper&   helper )
 {
-  if ( _pos.empty() )
+  if ( Is( SHRUNK ))
     return false; // already at the target position
 
   SMDS_MeshNode* tgtNode = const_cast< SMDS_MeshNode*& >( _nodes.back() );
@@ -10163,6 +11072,10 @@ bool _LayerEdge::SetNewLength2d( Handle(Geom_Surface)& surface,
     double stepSize = 1e100;
     for ( size_t i = 0; i < _simplices.size(); ++i )
     {
+      if ( !_simplices[i]._nPrev->isMarked() ||
+           !_simplices[i]._nNext->isMarked() )
+        continue; // simplex of quadrangle created by addBoundaryElements()
+
       // 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 );
@@ -10178,7 +11091,8 @@ bool _LayerEdge::SetNewLength2d( Handle(Geom_Surface)& surface,
     if ( uvLen <= stepSize )
     {
       newUV = tgtUV;
-      _pos.clear();
+      Set( SHRUNK );
+      //_pos.clear();
     }
     else if ( stepSize > 0 )
     {
@@ -10211,7 +11125,8 @@ bool _LayerEdge::SetNewLength2d( Handle(Geom_Surface)& surface,
     double newU = _pos[0].Coord( U_TGT );
     if ( lenTgt < 0.99 * fabs( uSrc-u2 )) // n2 got out of src-tgt range
     {
-      _pos.clear();
+      Set( _LayerEdge::SHRUNK );
+      //_pos.clear();
     }
     else
     {
@@ -10444,7 +11359,7 @@ void _Shrinker1D::AddEdge( const _LayerEdge*   e,
     _done = false;
   }
   // check _LayerEdge
-  if ( e == _edges[0] || e == _edges[1] )
+  if ( e == _edges[0] || e == _edges[1] || e->_nodes.size() < 2 )
     return;
   if ( eos.SWOLType() != TopAbs_EDGE )
     throw SALOME_Exception(LOCALIZED("Wrong _LayerEdge is added"));
@@ -10481,9 +11396,18 @@ void _Shrinker1D::AddEdge( const _LayerEdge*   e,
     while ( nIt->more() )
     {
       const SMDS_MeshNode* node = nIt->next();
+
+      // skip refinement nodes
       if ( node->NbInverseElements(SMDSAbs_Edge) == 0 ||
            node == tgtNode0 || node == tgtNode1 )
-        continue; // refinement nodes
+        continue;
+      bool hasMarkedFace = false;
+      SMDS_ElemIteratorPtr fIt = node->GetInverseElementIterator(SMDSAbs_Face);
+      while ( fIt->more() && !hasMarkedFace )
+        hasMarkedFace = fIt->next()->isMarked();
+      if ( !hasMarkedFace )
+        continue;
+
       _nodes.push_back( node );
       _initU.push_back( helper.GetNodeU( _geomEdge, node ));
       double len = GCPnts_AbscissaPoint::Length(aCurve, f, _initU.back());
@@ -10516,8 +11440,8 @@ void _Shrinker1D::Compute(bool set3D, SMESH_MesherHelper& helper)
   if ( !e ) e = _edges[1];
   if ( !e ) return;
 
-  _done =  (( !_edges[0] || _edges[0]->_pos.empty() ) &&
-            ( !_edges[1] || _edges[1]->_pos.empty() ));
+  _done =  (( !_edges[0] || _edges[0]->Is( _LayerEdge::SHRUNK )) &&
+            ( !_edges[1] || _edges[1]->Is( _LayerEdge::SHRUNK )));
 
   double f,l;
   if ( set3D || _done )
@@ -10598,11 +11522,12 @@ void _Shrinker1D::SwapSrcTgtNodes( SMESHDS_Mesh* mesh )
     if ( !eSubMesh ) return;
     const SMDS_MeshNode* srcNode = _edges[i]->_nodes[0];
     const SMDS_MeshNode* tgtNode = _edges[i]->_nodes.back();
+    const SMDS_MeshNode* scdNode = _edges[i]->_nodes[1];
     SMDS_ElemIteratorPtr eIt = srcNode->GetInverseElementIterator(SMDSAbs_Edge);
     while ( eIt->more() )
     {
       const SMDS_MeshElement* e = eIt->next();
-      if ( !eSubMesh->Contains( e ))
+      if ( !eSubMesh->Contains( e ) || e->GetNodeIndex( scdNode ) >= 0 )
           continue;
       SMDS_ElemIteratorPtr nIt = e->nodesIterator();
       for ( int iN = 0; iN < e->NbNodes(); ++iN )
@@ -10621,21 +11546,22 @@ void _Shrinker1D::SwapSrcTgtNodes( SMESHDS_Mesh* mesh )
  */
 //================================================================================
 
-bool _ViscousBuilder::addBoundaryElements()
+bool _ViscousBuilder::addBoundaryElements(_SolidData& data)
 {
   SMESH_MesherHelper helper( *_mesh );
 
   vector< const SMDS_MeshNode* > faceNodes;
 
-  for ( size_t i = 0; i < _sdVec.size(); ++i )
+  //for ( size_t i = 0; i < _sdVec.size(); ++i )
   {
-    _SolidData& data = _sdVec[i];
+    //_SolidData& data = _sdVec[i];
     TopTools_IndexedMapOfShape geomEdges;
     TopExp::MapShapes( data._solid, TopAbs_EDGE, geomEdges );
     for ( int iE = 1; iE <= geomEdges.Extent(); ++iE )
     {
       const TopoDS_Edge& E = TopoDS::Edge( geomEdges(iE));
-      if ( data._noShrinkShapes.count( getMeshDS()->ShapeToIndex( E )))
+      const TGeomID edgeID = getMeshDS()->ShapeToIndex( E );
+      if ( data._noShrinkShapes.count( edgeID ))
         continue;
 
       // Get _LayerEdge's based on E
@@ -10684,10 +11610,9 @@ bool _ViscousBuilder::addBoundaryElements()
       // Find out orientation and type of face to create
 
       bool reverse = false, isOnFace;
-      
-      map< TGeomID, TopoDS_Shape >::iterator e2f =
-        data._shrinkShape2Shape.find( getMeshDS()->ShapeToIndex( E ));
       TopoDS_Shape F;
+
+      map< TGeomID, TopoDS_Shape >::iterator e2f = data._shrinkShape2Shape.find( edgeID );
       if (( isOnFace = ( e2f != data._shrinkShape2Shape.end() )))
       {
         F = e2f->second.Oriented( TopAbs_FORWARD );
@@ -10697,17 +11622,12 @@ bool _ViscousBuilder::addBoundaryElements()
         if ( helper.IsReversedSubMesh( TopoDS::Face(F) ))
           reverse = !reverse;
       }
-      else
+      else if ( !data._ignoreFaceIds.count( e2f->first ))
       {
         // find FACE with layers sharing E
-        PShapeIteratorPtr fIt = helper.GetAncestors( E, *_mesh, TopAbs_FACE );
-        while ( fIt->more() && F.IsNull() )
-        {
-          const TopoDS_Shape* pF = fIt->next();
-          if ( helper.IsSubShape( *pF, data._solid) &&
-               !data._ignoreFaceIds.count( e2f->first ))
-            F = *pF;
-        }
+        PShapeIteratorPtr fIt = helper.GetAncestors( E, *_mesh, TopAbs_FACE, &data._solid );
+        if ( fIt->more() )
+          F = *( fIt->next() );
       }
       // Find the sub-mesh to add new faces
       SMESHDS_SubMesh* sm = 0;
@@ -10718,18 +11638,44 @@ bool _ViscousBuilder::addBoundaryElements()
       if ( !sm )
         return error("error in addBoundaryElements()", data._index);
 
+      // Find a proxy sub-mesh of the FACE of an adjacent SOLID, which will use the new boundary
+      // faces for 3D meshing (PAL23414)
+      SMESHDS_SubMesh* adjSM = 0;
+      if ( isOnFace )
+      {
+        const TGeomID   faceID = sm->GetID();
+        PShapeIteratorPtr soIt = helper.GetAncestors( F, *_mesh, TopAbs_SOLID );
+        while ( const TopoDS_Shape* solid = soIt->next() )
+          if ( !solid->IsSame( data._solid ))
+          {
+            size_t iData = _solids.FindIndex( *solid ) - 1;
+            if ( iData < _sdVec.size() &&
+                 _sdVec[ iData ]._ignoreFaceIds.count( faceID ) &&
+                 _sdVec[ iData ]._shrinkShape2Shape.count( edgeID ) == 0 )
+            {
+              SMESH_ProxyMesh::SubMesh* proxySub =
+                _sdVec[ iData ]._proxyMesh->getFaceSubM( TopoDS::Face( F ), /*create=*/false);
+              if ( proxySub && proxySub->NbElements() > 0 )
+                adjSM = proxySub;
+            }
+          }
+      }
+
       // Make faces
       const int dj1 = reverse ? 0 : 1;
       const int dj2 = reverse ? 1 : 0;
+      vector< const SMDS_MeshElement*> ff; // new faces row
+      SMESHDS_Mesh* m = getMeshDS();
       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;
+        ff.resize( std::max( nn1.size(), nn2.size() ), NULL );
         if ( nn1.size() == nn2.size() )
         {
           if ( isOnFace )
             for ( size_t z = 1; z < nn1.size(); ++z )
-              sm->AddElement( getMeshDS()->AddFace( nn1[z-1], nn2[z-1], nn2[z], nn1[z] ));
+              sm->AddElement( ff[z-1] = m->AddFace( nn1[z-1], nn2[z-1], nn2[z], nn1[z] ));
           else
             for ( size_t z = 1; z < nn1.size(); ++z )
               sm->AddElement( new SMDS_FaceOfNodes( nn1[z-1], nn2[z-1], nn2[z], nn1[z] ));
@@ -10738,7 +11684,7 @@ bool _ViscousBuilder::addBoundaryElements()
         {
           if ( isOnFace )
             for ( size_t z = 1; z < nn2.size(); ++z )
-              sm->AddElement( getMeshDS()->AddFace( nn1[0], nn2[z-1], nn2[z] ));
+              sm->AddElement( ff[z-1] = m->AddFace( nn1[0], nn2[z-1], nn2[z] ));
           else
             for ( size_t z = 1; z < nn2.size(); ++z )
               sm->AddElement( new SMDS_FaceOfNodes( nn1[0], nn2[z-1], nn2[z] ));
@@ -10747,11 +11693,19 @@ bool _ViscousBuilder::addBoundaryElements()
         {
           if ( isOnFace )
             for ( size_t z = 1; z < nn1.size(); ++z )
-              sm->AddElement( getMeshDS()->AddFace( nn1[z-1], nn2[0], nn1[z] ));
+              sm->AddElement( ff[z-1] = m->AddFace( nn1[z-1], nn2[0], nn1[z] ));
           else
             for ( size_t z = 1; z < nn1.size(); ++z )
               sm->AddElement( new SMDS_FaceOfNodes( nn1[z-1], nn2[0], nn2[z] ));
         }
+
+        if ( adjSM ) // add faces to a proxy SM of the adjacent SOLID
+        {
+          for ( size_t z = 0; z < ff.size(); ++z )
+            if ( ff[ z ])
+              adjSM->AddElement( ff[ z ]);
+          ff.clear();
+        }
       }
 
       // Make edges
@@ -10762,7 +11716,7 @@ bool _ViscousBuilder::addBoundaryElements()
         if ( eos && eos->SWOLType() == TopAbs_EDGE )
         {
           vector< const SMDS_MeshNode*>&  nn = edge->_nodes;
-          if ( nn.size() < 2 || nn[1]->GetInverseElementIterator( SMDSAbs_Edge )->more() )
+          if ( nn.size() < 2 || nn[1]->NbInverseElements( SMDSAbs_Edge ) >= 2 )
             continue;
           helper.SetSubShape( eos->_sWOL );
           helper.SetElementsOnShape( true );