Salome HOME
Revert "Merge branch 'yan/parallel_mesh2'"
[modules/smesh.git] / src / StdMeshers / StdMeshers_ViscousLayers.cxx
index 075d1b46d68c1f32faaff7fac158250546923410..680708d7ea81b1578e264bf8870c206a848ef67a 100644 (file)
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2019  CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2022  CEA/DEN, EDF R&D, OPEN CASCADE
 //
 // This library is free software; you can redistribute it and/or
 // modify it under the terms of the GNU Lesser General Public
 #include "SMESH_subMeshEventListener.hxx"
 #include "StdMeshers_FaceSide.hxx"
 #include "StdMeshers_ProjectionUtils.hxx"
+#include "StdMeshers_Quadrangle_2D.hxx"
 #include "StdMeshers_ViscousLayers2D.hxx"
 
 #include <Adaptor3d_HSurface.hxx>
 #include <BRepAdaptor_Curve.hxx>
 #include <BRepAdaptor_Curve2d.hxx>
 #include <BRepAdaptor_Surface.hxx>
-//#include <BRepLProp_CLProps.hxx>
 #include <BRepLProp_SLProps.hxx>
 #include <BRepOffsetAPI_MakeOffsetShape.hxx>
 #include <BRep_Tool.hxx>
 #include <unordered_map>
 
 #ifdef _DEBUG_
-//#define __myDEBUG
+#ifndef WIN32
+#define __myDEBUG
+#endif
 //#define __NOT_INVALIDATE_BAD_SMOOTH
 //#define __NODES_AT_POS
 #endif
@@ -206,8 +208,8 @@ namespace VISCOUS_3D
     virtual void ProcessEvent(const int                       event,
                               const int                       eventType,
                               SMESH_subMesh*                  subMesh,
-                              SMESH_subMeshEventListenerData* data,
-                              const SMESH_Hypothesis*         hyp)
+                              SMESH_subMeshEventListenerData* /*data*/,
+                              const SMESH_Hypothesis*         /*hyp*/)
     {
       if (( SMESH_subMesh::COMPUTE_EVENT       == eventType ) &&
           ( SMESH_subMesh::CHECK_COMPUTE_STATE != event &&
@@ -385,6 +387,7 @@ namespace VISCOUS_3D
   struct _LayerEdge;
   struct _EdgesOnShape;
   struct _Smoother1D;
+  struct _Mapper2D;
   typedef map< const SMDS_MeshNode*, _LayerEdge*, TIDCompare > TNode2Edge;
 
   //--------------------------------------------------------------------------------
@@ -444,6 +447,10 @@ namespace VISCOUS_3D
                          const TopoDS_Face&    F,
                          _EdgesOnShape&        eos,
                          SMESH_MesherHelper&   helper );
+    bool UpdatePositionOnSWOL( SMDS_MeshNode*      n,
+                               double              tol,
+                               _EdgesOnShape&      eos,
+                               SMESH_MesherHelper& helper );
     void SetDataByNeighbors( const SMDS_MeshNode* n1,
                              const SMDS_MeshNode* n2,
                              const _EdgesOnShape& eos,
@@ -494,7 +501,7 @@ namespace VISCOUS_3D
     bool   IsOnFace() const { return ( _nodes[0]->GetPosition()->GetDim() == 2 ); }
     int    BaseShapeDim() const { return _nodes[0]->GetPosition()->GetDim(); }
     gp_XYZ Copy( _LayerEdge& other, _EdgesOnShape& eos, SMESH_MesherHelper& helper );
-    void   SetCosin( double cosin );
+    double SetCosin( double cosin );
     void   SetNormal( const gp_XYZ& n ) { _normal = n; }
     void   SetMaxLen( double l ) { _maxLen = l; }
     int    NbSteps() const { return _pos.size() - 1; } // nb inlation steps
@@ -612,12 +619,30 @@ namespace VISCOUS_3D
         _thickness      = Max( _thickness, hyp->GetTotalThickness() );
         _stretchFactor += hyp->GetStretchFactor();
         _method         = hyp->GetMethod();
+        if ( _groupName.empty() )
+          _groupName = hyp->GetGroupName();
       }
     }
     double GetTotalThickness() const { return _thickness; /*_nbHyps ? _thickness / _nbHyps : 0;*/ }
     double GetStretchFactor()  const { return _nbHyps ? _stretchFactor / _nbHyps : 0; }
     int    GetNumberLayers()   const { return _nbLayers; }
     int    GetMethod()         const { return _method; }
+    bool   ToCreateGroup()     const { return !_groupName.empty(); }
+    const std::string& GetGroupName() const { return _groupName; }
+
+    double Get1stLayerThickness( double realThickness = 0.) const
+    {
+      const double T = ( realThickness > 0 ) ? realThickness : GetTotalThickness();
+      const double f = GetStretchFactor();
+      const int    N = GetNumberLayers();
+      const double fPowN = pow( f, N );
+      double h0;
+      if ( fPowN - 1 <= numeric_limits<double>::min() )
+        h0 = T / N;
+      else
+        h0 = T * ( f - 1 )/( fPowN - 1 );
+      return h0;
+    }
 
     bool   UseSurfaceNormal()  const
     { return _method == StdMeshers_ViscousLayers::SURF_OFFSET_SMOOTH; }
@@ -636,8 +661,9 @@ namespace VISCOUS_3D
     static bool Equals( double v1, double v2 ) { return Abs( v1 - v2 ) < 0.01 * ( v1 + v2 ); }
 
   private:
-    int     _nbLayers, _nbHyps, _method;
-    double  _thickness, _stretchFactor;
+    int         _nbLayers, _nbHyps, _method;
+    double      _thickness, _stretchFactor;
+    std::string _groupName;
   };
 
   //--------------------------------------------------------------------------------
@@ -667,6 +693,8 @@ namespace VISCOUS_3D
 
     Handle(ShapeAnalysis_Surface) _offsetSurf;
     _LayerEdge*                   _edgeForOffset;
+    double                        _offsetValue;
+    _Mapper2D*                    _mapper2D;
 
     _SolidData*            _data; // parent SOLID
 
@@ -680,8 +708,11 @@ namespace VISCOUS_3D
     { return std::find( _eosC1.begin(), _eosC1.end(), other ) != _eosC1.end(); }
     bool             GetNormal( const SMDS_MeshElement* face, gp_Vec& norm );
     _SolidData&      GetData() const { return *_data; }
+    char             ShapeTypeLetter() const
+    { switch ( ShapeType() ) { case TopAbs_FACE: return 'F'; case TopAbs_EDGE: return 'E';
+      case TopAbs_VERTEX: return 'V'; default: return 'S'; }}
 
-    _EdgesOnShape(): _shapeID(-1), _subMesh(0), _toSmooth(false), _edgeSmoother(0) {}
+    _EdgesOnShape(): _shapeID(-1), _subMesh(0), _toSmooth(false), _edgeSmoother(0), _mapper2D(0) {}
     ~_EdgesOnShape();
   };
 
@@ -912,7 +943,7 @@ namespace VISCOUS_3D
 
   private:
 
-    bool findSolidsWithLayers();
+    bool findSolidsWithLayers(const bool checkFaceMesh=true);
     bool setBefore( _SolidData& solidBefore, _SolidData& solidAfter );
     bool findFacesWithLayers(const bool onlyWith=false);
     void findPeriodicFaces();
@@ -920,7 +951,7 @@ namespace VISCOUS_3D
                         const StdMeshers_ViscousLayers* hyp,
                         const TopoDS_Shape&             hypShape,
                         set<TGeomID>&                   ignoreFaces);
-    void makeEdgesOnShape();
+    int makeEdgesOnShape();
     bool makeLayer(_SolidData& data);
     void setShapeData( _EdgesOnShape& eos, SMESH_subMesh* sm, _SolidData& data );
     bool setEdgeData( _LayerEdge& edge, _EdgesOnShape& eos,
@@ -1086,7 +1117,7 @@ namespace VISCOUS_3D
                              Handle(ShapeAnalysis_Surface)& surface,
                              const TopoDS_Face&             F,
                              SMESH_MesherHelper&            helper);
-    bool smoothComplexEdge( _SolidData&                    data,
+    bool smoothComplexEdge( _SolidData&                     data,
                             Handle(ShapeAnalysis_Surface)& surface,
                             const TopoDS_Face&             F,
                             SMESH_MesherHelper&            helper);
@@ -1100,6 +1131,22 @@ namespace VISCOUS_3D
 
     void offPointsToPython() const; // debug
   };
+
+  //--------------------------------------------------------------------------------
+  /*!
+   * \brief Compute positions of nodes of 2D structured mesh using TFI
+   */
+  class _Mapper2D
+  {
+    FaceQuadStruct _quadPoints;
+
+    UVPtStruct& uvPnt( size_t i, size_t j ) { return _quadPoints.UVPt( i, j ); }
+
+  public:
+    _Mapper2D( const TParam2ColumnMap & param2ColumnMap, const TNode2Edge& n2eMap );
+    bool ComputeNodePositions();
+  };
+
   //--------------------------------------------------------------------------------
   /*!
    * \brief Class of temporary mesh face.
@@ -1247,7 +1294,8 @@ namespace VISCOUS_3D
 StdMeshers_ViscousLayers::StdMeshers_ViscousLayers(int hypId, SMESH_Gen* gen)
   :SMESH_Hypothesis(hypId, gen),
    _isToIgnoreShapes(1), _nbLayers(1), _thickness(1), _stretchFactor(1),
-   _method( SURF_OFFSET_SMOOTH )
+   _method( SURF_OFFSET_SMOOTH ),
+   _groupName("")
 {
   _name = StdMeshers_ViscousLayers::GetHypType();
   _param_algo_dim = -3; // auxiliary hyp used by 3D algos
@@ -1279,6 +1327,15 @@ void StdMeshers_ViscousLayers::SetMethod( ExtrusionMethod method )
   if ( _method != method )
     _method = method, NotifySubMeshesHypothesisModification();
 } // --------------------------------------------------------------------------------
+void StdMeshers_ViscousLayers::SetGroupName(const std::string& name)
+{
+  if ( _groupName != name )
+  {
+    _groupName = name;
+    if ( !_groupName.empty() )
+      NotifySubMeshesHypothesisModification();
+  }
+} // --------------------------------------------------------------------------------
 SMESH_ProxyMesh::Ptr
 StdMeshers_ViscousLayers::Compute(SMESH_Mesh&         theMesh,
                                   const TopoDS_Shape& theShape,
@@ -1333,6 +1390,9 @@ std::ostream & StdMeshers_ViscousLayers::SaveTo(std::ostream & save)
     save << " " << _shapeIds[i];
   save << " " << !_isToIgnoreShapes; // negate to keep the behavior in old studies.
   save << " " << _method;
+  save << " " << _groupName.size();
+  if ( !_groupName.empty() )
+    save << " " << _groupName;
   return save;
 } // --------------------------------------------------------------------------------
 std::istream & StdMeshers_ViscousLayers::LoadFrom(std::istream & load)
@@ -1345,14 +1405,21 @@ std::istream & StdMeshers_ViscousLayers::LoadFrom(std::istream & load)
     _isToIgnoreShapes = !shapeToTreat;
     if ( load >> method )
       _method = (ExtrusionMethod) method;
+    int nameSize = 0;
+    if ( load >> nameSize && nameSize > 0 )
+    {
+      _groupName.resize( nameSize );
+      load.get( _groupName[0] ); // remove a white-space
+      load.getline( &_groupName[0], nameSize + 1 );
+    }
   }
   else {
     _isToIgnoreShapes = true; // old behavior
   }
   return load;
 } // --------------------------------------------------------------------------------
-bool StdMeshers_ViscousLayers::SetParametersByMesh(const SMESH_Mesh*   theMesh,
-                                                   const TopoDS_Shape& theShape)
+bool StdMeshers_ViscousLayers::SetParametersByMesh(const SMESH_Mesh*   /*theMesh*/,
+                                                   const TopoDS_Shape& /*theShape*/)
 {
   // TODO
   return false;
@@ -1378,22 +1445,77 @@ bool StdMeshers_ViscousLayers::IsShapeWithLayers(int shapeIndex) const
     ( std::find( _shapeIds.begin(), _shapeIds.end(), shapeIndex ) != _shapeIds.end() );
   return IsToIgnoreShapes() ? !isIn : isIn;
 }
+
+// --------------------------------------------------------------------------------
+SMDS_MeshGroup* StdMeshers_ViscousLayers::CreateGroup( const std::string&  theName,
+                                                       SMESH_Mesh&         theMesh,
+                                                       SMDSAbs_ElementType theType)
+{
+  SMESH_Group*      group = 0;
+  SMDS_MeshGroup* groupDS = 0;
+
+  if ( theName.empty() )
+    return groupDS;
+       
+  if ( SMESH_Mesh::GroupIteratorPtr grIt = theMesh.GetGroups() )
+    while( grIt->more() && !group )
+    {
+      group = grIt->next();
+      if ( !group ||
+           group->GetGroupDS()->GetType() != theType ||
+           group->GetName()               != theName ||
+           !dynamic_cast< SMESHDS_Group* >( group->GetGroupDS() ))
+        group = 0;
+    }
+  if ( !group )
+    group = theMesh.AddGroup( theType, theName.c_str() );
+
+  groupDS = & dynamic_cast< SMESHDS_Group* >( group->GetGroupDS() )->SMDSGroup();
+
+  return groupDS;
+}
+
 // END StdMeshers_ViscousLayers hypothesis
 //================================================================================
 
 namespace VISCOUS_3D
 {
-  gp_XYZ getEdgeDir( const TopoDS_Edge& E, const TopoDS_Vertex& fromV )
+  gp_XYZ getEdgeDir( const TopoDS_Edge& E, const TopoDS_Vertex& fromV,
+                     const double h0, bool* isRegularEdge = nullptr )
   {
     gp_Vec dir;
     double f,l;
     Handle(Geom_Curve) c = BRep_Tool::Curve( E, f, l );
     if ( c.IsNull() ) return gp_XYZ( Precision::Infinite(), 1e100, 1e100 );
-    gp_Pnt p = BRep_Tool::Pnt( fromV );
-    double distF = p.SquareDistance( c->Value( f ));
-    double distL = p.SquareDistance( c->Value( l ));
+    gp_Pnt  p = BRep_Tool::Pnt( fromV );
+    gp_Pnt pf = c->Value( f ), pl = c->Value( l );
+    double distF = p.SquareDistance( pf );
+    double distL = p.SquareDistance( pl );
     c->D1(( distF < distL ? f : l), p, dir );
     if ( distL < distF ) dir.Reverse();
+    bool isDifficult = false;
+    if ( dir.SquareMagnitude() < h0 * h0 ) // check dir orientation
+    {
+      gp_Pnt& pClose = distF < distL ? pf : pl;
+      gp_Pnt&   pFar = distF < distL ? pl : pf;
+      gp_Pnt    pMid = 0.9 * pClose.XYZ() + 0.1 * pFar.XYZ();
+      gp_Vec vMid( p, pMid );
+      double     dot = vMid * dir;
+      double    cos2 = dot * dot / dir.SquareMagnitude() / vMid.SquareMagnitude();
+      if ( cos2 < 0.7 * 0.7 || dot < 0 ) // large angle between dir and vMid
+      {
+        double uClose = distF < distL ? f : l;
+        double   uFar = distF < distL ? l : f;
+        double      r = h0 / SMESH_Algo::EdgeLength( E );
+        double   uMid = ( 1 - r ) * uClose + r * uFar;
+        pMid = c->Value( uMid );
+        dir = gp_Vec( p, pMid );
+        isDifficult = true;
+      }
+    }
+    if ( isRegularEdge )
+      *isRegularEdge = !isDifficult;
+
     return dir.XYZ();
   }
   //--------------------------------------------------------------------------------
@@ -1410,8 +1532,8 @@ namespace VISCOUS_3D
   }
   //--------------------------------------------------------------------------------
   gp_XYZ getFaceDir( const TopoDS_Face& F, const TopoDS_Vertex& fromV,
-                     const SMDS_MeshNode* node, SMESH_MesherHelper& helper, bool& ok,
-                     double* cosin=0);
+                     const SMDS_MeshNode* node, SMESH_MesherHelper& helper, bool& ok/*,
+                     double* cosin=0*/);
   //--------------------------------------------------------------------------------
   gp_XYZ getFaceDir( const TopoDS_Face& F, const TopoDS_Edge& fromE,
                      const SMDS_MeshNode* node, SMESH_MesherHelper& helper, bool& ok)
@@ -1452,7 +1574,7 @@ namespace VISCOUS_3D
   //--------------------------------------------------------------------------------
   gp_XYZ getFaceDir( const TopoDS_Face& F, const TopoDS_Vertex& fromV,
                      const SMDS_MeshNode* node, SMESH_MesherHelper& helper,
-                     bool& ok, double* cosin)
+                     bool& ok/*, double* cosin*/)
   {
     TopoDS_Face faceFrw = F;
     faceFrw.Orientation( TopAbs_FORWARD );
@@ -1484,7 +1606,7 @@ namespace VISCOUS_3D
       ok = true;
       for ( size_t i = 0; i < nbEdges && ok; ++i )
       {
-        edgeDir[i] = getEdgeDir( edges[i], fromV );
+        edgeDir[i] = getEdgeDir( edges[i], fromV, 0.1 * SMESH_Algo::EdgeLength( edges[i] ));
         double size2 = edgeDir[i].SquareModulus();
         if (( ok = size2 > numeric_limits<double>::min() ))
           edgeDir[i] /= sqrt( size2 );
@@ -1504,15 +1626,15 @@ namespace VISCOUS_3D
         if ( angle < 0 )
           dir.Reverse();
       }
-      if ( cosin ) {
-        double angle = gp_Vec( edgeDir[0] ).Angle( dir );
-        *cosin = Cos( angle );
-      }
+      // if ( cosin ) {
+      //   double angle = faceNormal.Angle( dir );
+      //   *cosin = Cos( angle );
+      // }
     }
     else if ( nbEdges == 1 )
     {
       dir = getFaceDir( faceFrw, edges[ edges[0].IsNull() ], node, helper, ok );
-      if ( cosin ) *cosin = 1.;
+      //if ( cosin ) *cosin = 1.;
     }
     else
     {
@@ -1706,8 +1828,8 @@ namespace VISCOUS_3D
 
   //--------------------------------------------------------------------------------
   // DEBUG. Dump intermediate node positions into a python script
-  // HOWTO use: run python commands written in a console to see
-  //  construction steps of viscous layers
+  // HOWTO use: run python commands written in a console and defined in /tmp/viscous.py
+  // to see construction steps of viscous layers
 #ifdef __myDEBUG
   ostream* py;
   int      theNbPyFunc;
@@ -1716,7 +1838,7 @@ namespace VISCOUS_3D
     PyDump(SMESH_Mesh& m) {
       int tag = 3 + m.GetId();
       const char* fname = "/tmp/viscous.py";
-      cout << "execfile('"<<fname<<"')"<<endl;
+      cout << "exec(open('"<<fname<<"','rb').read() )"<<endl;
       py = _pyStream = new ofstream(fname);
       *py << "import SMESH" << endl
           << "from salome.smesh import smeshBuilder" << endl
@@ -1908,9 +2030,6 @@ SMESH_ComputeErrorPtr _ViscousBuilder::Compute(SMESH_Mesh&         theMesh,
   if ( _ViscousListener::GetSolidMesh( _mesh, exp.Current(), /*toCreate=*/false))
     return SMESH_ComputeErrorPtr(); // everything already computed
 
-  PyDump debugDump( theMesh );
-  _pyDump = &debugDump;
-
   // TODO: ignore already computed SOLIDs
   if ( !findSolidsWithLayers())
     return _error;
@@ -1918,16 +2037,15 @@ SMESH_ComputeErrorPtr _ViscousBuilder::Compute(SMESH_Mesh&         theMesh,
   if ( !findFacesWithLayers() )
     return _error;
 
-  // for ( size_t i = 0; i < _sdVec.size(); ++i )
-  // {
-  //   if ( ! makeLayer( _sdVec[ i ]))   // create _LayerEdge's
-  //     return _error;
-  // }
-
-  makeEdgesOnShape();
+  if ( !makeEdgesOnShape() )
+    return _error;
 
   findPeriodicFaces();
 
+  PyDump debugDump( theMesh );
+  _pyDump = &debugDump;
+
+
   for ( size_t i = 0; i < _sdVec.size(); ++i )
   {
     size_t iSD = 0;
@@ -1936,6 +2054,8 @@ SMESH_ComputeErrorPtr _ViscousBuilder::Compute(SMESH_Mesh&         theMesh,
            !_sdVec[iSD]._solid.IsNull() &&
            !_sdVec[iSD]._done )
         break;
+    if ( iSD == _sdVec.size() )
+      break; // all done
 
     if ( ! makeLayer(_sdVec[iSD]) )   // create _LayerEdge's
       return _error;
@@ -1985,7 +2105,7 @@ SMESH_ComputeErrorPtr _ViscousBuilder::CheckHypotheses( SMESH_Mesh&         mesh
     return SMESH_ComputeErrorPtr(); // everything already computed
 
 
-  findSolidsWithLayers();
+  findSolidsWithLayers( /*checkFaceMesh=*/false );
   bool ok = findFacesWithLayers( true );
 
   // remove _MeshOfSolid's of _SolidData's
@@ -2004,7 +2124,7 @@ SMESH_ComputeErrorPtr _ViscousBuilder::CheckHypotheses( SMESH_Mesh&         mesh
  */
 //================================================================================
 
-bool _ViscousBuilder::findSolidsWithLayers()
+bool _ViscousBuilder::findSolidsWithLayers(const bool checkFaceMesh)
 {
   // get all solids
   TopTools_IndexedMapOfShape allSolids;
@@ -2014,13 +2134,28 @@ bool _ViscousBuilder::findSolidsWithLayers()
   SMESH_HypoFilter filter;
   for ( int i = 1; i <= allSolids.Extent(); ++i )
   {
-    // find StdMeshers_ViscousLayers hyp assigned to the i-th solid
     SMESH_subMesh* sm = _mesh->GetSubMesh( allSolids(i) );
     if ( sm->GetSubMeshDS() && sm->GetSubMeshDS()->NbElements() > 0 )
       continue; // solid is already meshed
+    // TODO: check if algo is hidden
     SMESH_Algo* algo = sm->GetAlgo();
     if ( !algo ) continue;
-    // TODO: check if algo is hidden
+    // check if all FACEs are meshed, which can be false if Compute() a sub-shape
+    if ( checkFaceMesh )
+    {
+      bool facesMeshed = true;
+      SMESH_subMeshIteratorPtr smIt = sm->getDependsOnIterator(false,true);
+      while ( smIt->more() && facesMeshed )
+      {
+        SMESH_subMesh * faceSM = smIt->next();
+        if ( faceSM->GetSubShape().ShapeType() != TopAbs_FACE )
+          break;
+        facesMeshed = faceSM->IsMeshComputed();
+      }
+      if ( !facesMeshed )
+        continue;
+    }
+    // find StdMeshers_ViscousLayers hyp assigned to the i-th solid
     const list <const SMESHDS_Hypothesis *> & allHyps =
       algo->GetUsedHypothesis(*_mesh, allSolids(i), /*ignoreAuxiliary=*/false);
     _SolidData* soData = 0;
@@ -2547,6 +2682,7 @@ bool _ViscousBuilder::makeLayer(_SolidData& data)
 
   // Create temporary faces and _LayerEdge's
 
+  debugMsg( "######################" );
   dumpFunction(SMESH_Comment("makeLayers_")<<data._index);
 
   vector< _EdgesOnShape >& edgesByGeom = data._edgesOnShape;
@@ -2607,7 +2743,7 @@ bool _ViscousBuilder::makeLayer(_SolidData& data)
             nbDegenNodes++;
           }
         }
-        TNode2Edge::iterator n2e = data._n2eMap.insert( make_pair( n, (_LayerEdge*)0 )).first;
+        TNode2Edge::iterator n2e = data._n2eMap.insert({ n, nullptr }).first;
         if ( !(*n2e).second )
         {
           // add a _LayerEdge
@@ -2641,13 +2777,12 @@ bool _ViscousBuilder::makeLayer(_SolidData& data)
             if ( !setEdgeData( *edge, edgesByGeom[ shapeID ], helper, data ))
               return false;
 
-            if ( edge->_nodes.size() < 2 )
-              edge->Block( data );
-              //data._noShrinkShapes.insert( shapeID );
+            if ( edge->_nodes.size() < 2 && !noShrink )
+              edge->Block( data ); // a sole node is moved only if noShrink
           }
           dumpMove(edge->_nodes.back());
 
-          if ( edge->_cosin > faceMaxCosin )
+          if ( edge->_cosin > faceMaxCosin && edge->_nodes.size() > 1 )
           {
             faceMaxCosin = edge->_cosin;
             maxCosinEdge = edge;
@@ -2697,6 +2832,42 @@ bool _ViscousBuilder::makeLayer(_SolidData& data)
         if (( n2e = data._n2eMap.find( *node )) != data._n2eMap.end() )
           edge->_neibors.push_back( n2e->second );
     }
+
+    // Fix uv of nodes on periodic FACEs (bos #20643)
+
+    if ( eos.ShapeType() != TopAbs_EDGE ||
+         eos.SWOLType()  != TopAbs_FACE ||
+         eos.size() == 0 )
+      continue;
+
+    const TopoDS_Face& F = TopoDS::Face( eos._sWOL );
+    SMESH_MesherHelper faceHelper( *_mesh );
+    faceHelper.SetSubShape( F );
+    faceHelper.ToFixNodeParameters( true );
+    if ( faceHelper.GetPeriodicIndex() == 0 )
+      continue;
+
+    SMESHDS_SubMesh* smDS = getMeshDS()->MeshElements( F );
+    if ( !smDS || smDS->GetNodes() == 0 )
+      continue;
+
+    bool toCheck = true;
+    const double tol = 2 * helper.MaxTolerance( F );
+    for ( SMDS_NodeIteratorPtr nIt = smDS->GetNodes(); nIt->more(); )
+    {
+      const SMDS_MeshNode* node = nIt->next();
+      gp_XY uvNew( Precision::Infinite(), 0 );
+      if ( toCheck )
+      {
+        toCheck = false;
+        gp_XY uv = faceHelper.GetNodeUV( F, node );
+        if ( ! faceHelper.CheckNodeUV( F, node, uvNew, tol, /*force=*/true ))
+          break; // projection on F failed
+        if (( uv - uvNew ).Modulus() < Precision::Confusion() )
+          break; // current uv is OK
+      }
+      faceHelper.CheckNodeUV( F, node, uvNew, tol, /*force=*/true );
+    }
   }
 
   data._epsilon = 1e-7;
@@ -2867,6 +3038,7 @@ void _ViscousBuilder::limitStepSizeByCurvature( _SolidData& data )
     cnvFace._face = F;
     cnvFace._normalsFixed = false;
     cnvFace._isTooCurved = false;
+    cnvFace._normalsFixedOnBorders = false;
 
     double maxCurvature = cnvFace.GetMaxCurvature( data, eof, surfProp, helper );
     if ( maxCurvature > 0 )
@@ -3050,13 +3222,13 @@ bool _ViscousBuilder::findShapesToSmooth( _SolidData& data )
     if ( SMESH_Algo::isDegenerated( E ) || !edgesOfSmooFaces.Contains( E ))
       continue;
 
-    double tgtThick = eos._hyp.GetTotalThickness();
+    double tgtThick = eos._hyp.GetTotalThickness(), h0 = eos._hyp.Get1stLayerThickness();
     for ( TopoDS_Iterator vIt( E ); vIt.More() && !eos._toSmooth; vIt.Next() )
     {
       TGeomID iV = getMeshDS()->ShapeToIndex( vIt.Value() );
       vector<_LayerEdge*>& eV = edgesByGeom[ iV ]._edges;
       if ( eV.empty() || eV[0]->Is( _LayerEdge::MULTI_NORMAL )) continue;
-      gp_Vec  eDir    = getEdgeDir( E, TopoDS::Vertex( vIt.Value() ));
+      gp_Vec  eDir    = getEdgeDir( E, TopoDS::Vertex( vIt.Value() ), h0 );
       double angle    = eDir.Angle( eV[0]->_normal );
       double cosin    = Cos( angle );
       double cosinAbs = Abs( cosin );
@@ -3234,7 +3406,7 @@ bool _ViscousBuilder::findShapesToSmooth( _SolidData& data )
     {
       _EdgesOnShape* eoe = data.GetShapeEdges( *e );
       if ( !eoe ) continue; // other solid
-      gp_XYZ eDir = getEdgeDir( TopoDS::Edge( *e ), V );
+      gp_XYZ eDir = getEdgeDir( TopoDS::Edge( *e ), V, eoe->_hyp.Get1stLayerThickness() );
       if ( !Precision::IsInfinite( eDir.X() ))
         dirOfEdges.push_back( make_pair( eoe, eDir.Normalized() ));
     }
@@ -3249,12 +3421,22 @@ bool _ViscousBuilder::findShapesToSmooth( _SolidData& data )
           if ( isC1 )
           {
             double maxEdgeLen = 3 * Min( eov._edges[0]->_maxLen, eov._hyp.GetTotalThickness() );
-            double eLen1 = SMESH_Algo::EdgeLength( TopoDS::Edge( dirOfEdges[i].first->_shape ));
-            double eLen2 = SMESH_Algo::EdgeLength( TopoDS::Edge( dirOfEdges[j].first->_shape ));
-            if ( eLen1 < maxEdgeLen ) eov._eosC1.push_back( dirOfEdges[i].first );
-            if ( eLen2 < maxEdgeLen ) eov._eosC1.push_back( dirOfEdges[j].first );
-            dirOfEdges[i].first = 0;
-            dirOfEdges[j].first = 0;
+            for ( int isJ = 0; isJ < 2; ++isJ ) // loop on [i,j]
+            {
+              size_t k = isJ ? j : i;
+              const TopoDS_Edge& e = TopoDS::Edge( dirOfEdges[k].first->_shape );
+              double eLen = SMESH_Algo::EdgeLength( e );
+              if ( eLen < maxEdgeLen )
+              {
+                TopoDS_Shape oppV = SMESH_MesherHelper::IthVertex( 0, e );
+                if ( oppV.IsSame( V ))
+                  oppV = SMESH_MesherHelper::IthVertex( 1, e );
+                _EdgesOnShape* eovOpp = data.GetShapeEdges( oppV );
+                if ( dirOfEdges[k].second * eovOpp->_edges[0]->_normal < 0 )
+                  eov._eosC1.push_back( dirOfEdges[k].first );
+              }
+              dirOfEdges[k].first = 0;
+            }
           }
         }
   } // fill _eosC1 of VERTEXes
@@ -3270,9 +3452,10 @@ bool _ViscousBuilder::findShapesToSmooth( _SolidData& data )
  */
 //================================================================================
 
-void _ViscousBuilder::makeEdgesOnShape()
+int _ViscousBuilder::makeEdgesOnShape()
 {
   const int nbShapes = getMeshDS()->MaxShapeIndex();
+  int nbSolidsWL = 0;
 
   for ( size_t i = 0; i < _sdVec.size(); ++i )
   {
@@ -3281,6 +3464,7 @@ void _ViscousBuilder::makeEdgesOnShape()
     edgesByGeom.resize( nbShapes+1 );
 
     // set data of _EdgesOnShape's
+    int nbShapesWL = 0;
     if ( SMESH_subMesh* sm = _mesh->GetSubMesh( data._solid ))
     {
       SMESH_subMeshIteratorPtr smIt = sm->getDependsOnIterator(/*includeSelf=*/false);
@@ -3292,9 +3476,21 @@ void _ViscousBuilder::makeEdgesOnShape()
           continue;
 
         setShapeData( edgesByGeom[ sm->GetId() ], sm, data );
+
+        nbShapesWL += ( sm->GetSubShape().ShapeType() == TopAbs_FACE );
       }
     }
+    if ( nbShapesWL == 0 ) // no shapes with layers in a SOLID
+    {
+      data._done = true;
+      SMESHUtils::FreeVector( edgesByGeom );
+    }
+    else
+    {
+      ++nbSolidsWL;
+    }
   }
+  return nbSolidsWL;
 }
 
 //================================================================================
@@ -3320,6 +3516,7 @@ void _ViscousBuilder::setShapeData( _EdgesOnShape& eos,
     eos._shape.Orientation( helper.GetSubShapeOri( data._solid, eos._shape ));
   eos._toSmooth = false;
   eos._data = &data;
+  eos._mapper2D = nullptr;
 
   // set _SWOL
   map< TGeomID, TopoDS_Shape >::const_iterator s2s =
@@ -3440,6 +3637,7 @@ bool _EdgesOnShape::GetNormal( const SMDS_MeshElement* face, gp_Vec& norm )
 _EdgesOnShape::~_EdgesOnShape()
 {
   delete _edgeSmoother;
+  delete _mapper2D;
 }
 
 //================================================================================
@@ -3521,13 +3719,14 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge&         edge,
       if ( eos.SWOLType() == TopAbs_EDGE )
       {
         // inflate from VERTEX along EDGE
-        edge._normal = getEdgeDir( TopoDS::Edge( eos._sWOL ), TopoDS::Vertex( eos._shape ));
+        edge._normal = getEdgeDir( TopoDS::Edge( eos._sWOL ), TopoDS::Vertex( eos._shape ),
+                                   eos._hyp.Get1stLayerThickness(), &eos._isRegularSWOL );
       }
       else if ( eos.ShapeType() == TopAbs_VERTEX )
       {
         // inflate from VERTEX along FACE
         edge._normal = getFaceDir( TopoDS::Face( eos._sWOL ), TopoDS::Vertex( eos._shape ),
-                                   node, helper, normOK, &edge._cosin);
+                                   node, helper, normOK/*, &edge._cosin*/);
       }
       else
       {
@@ -3616,30 +3815,37 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge&         edge,
       break;
     }
     case TopAbs_VERTEX: {
+      TopoDS_Vertex V  = TopoDS::Vertex( eos._shape );
+      gp_Vec inFaceDir = getFaceDir( face2Norm[0].first , V, node, helper, normOK );
+      double angle     = inFaceDir.Angle( edge._normal ); // [0,PI]
+      edge._cosin      = Cos( angle );
       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 = 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 ) > Abs( edge._cosin ))
-                edge._cosin = cosin;
+        totalNbFaces--;
+      if ( totalNbFaces > 1 || helper.IsSeamShape( node->getshapeId() ))
+        for ( int iF = 1; iF < totalNbFaces; ++iF )
+        {
+          F = face2Norm[ iF ].first;
+          inFaceDir = getFaceDir( F, V, node, helper, normOK=true );
+          if ( normOK ) {
+            if ( onShrinkShape )
+            {
+              gp_Vec faceNorm = getFaceNormal( node, F, helper, normOK );
+              if ( !normOK ) continue;
+              if ( helper.GetSubShapeOri( data._solid, F ) != TopAbs_REVERSED )
+                faceNorm.Reverse();
+              angle = 0.5 * M_PI - faceNorm.Angle( edge._normal );
+              if ( inFaceDir * edge._normal < 0 )
+                angle = M_PI - angle;
+            }
+            else
+            {
+              angle = inFaceDir.Angle( edge._normal );
             }
+            double cosin = Cos( angle );
+            if ( Abs( cosin ) > Abs( edge._cosin ))
+              edge._cosin = cosin;
           }
-      }
+        }
       break;
     }
     default:
@@ -3664,7 +3870,20 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge&         edge,
   // Set the rest data
   // --------------------
 
-  edge.SetCosin( edge._cosin ); // to update edge._lenFactor
+  double realLenFactor = edge.SetCosin( edge._cosin ); // to update edge._lenFactor
+  // if ( realLenFactor > 3 )
+  // {
+  //   edge._cosin = 1;
+  //   if ( onShrinkShape )
+  //   {
+  //     edge.Set( _LayerEdge::RISKY_SWOL );
+  //     edge._lenFactor = 2;
+  //   }
+  //   else
+  //   {
+  //     edge._lenFactor = 1;
+  //   }
+  // }
 
   if ( onShrinkShape )
   {
@@ -3688,7 +3907,7 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge&         edge,
         getMeshDS()->SetNodeOnFace( tgtNode, TopoDS::Face( eos._sWOL ), uv.X(), uv.Y() );
     }
 
-    if ( edge._nodes.size() > 1 )
+    //if ( edge._nodes.size() > 1 ) -- allow RISKY_SWOL on noShrink shape
     {
       // check if an angle between a FACE with layers and SWOL is sharp,
       // else the edge should not inflate
@@ -3703,10 +3922,14 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge&         edge,
           geomNorm.Reverse(); // inside the SOLID
         if ( geomNorm * edge._normal < -0.001 )
         {
-          getMeshDS()->RemoveFreeNode( tgtNode, 0, /*fromGroups=*/false );
-          edge._nodes.resize( 1 );
+          if ( edge._nodes.size() > 1 )
+          {
+            getMeshDS()->RemoveFreeNode( tgtNode, 0, /*fromGroups=*/false );
+            edge._nodes.resize( 1 );
+          }
         }
-        else if ( edge._lenFactor > 3 )
+        else if ( realLenFactor > 3 ) ///  -- moved to SetCosin()
+          //else if ( edge._lenFactor > 3 )
         {
           edge._lenFactor = 2;
           edge.Set( _LayerEdge::RISKY_SWOL );
@@ -3873,7 +4096,7 @@ gp_XYZ _ViscousBuilder::getFaceNormal(const SMDS_MeshNode* node,
 
 bool _ViscousBuilder::getFaceNormalAtSingularity( const gp_XY&        uv,
                                                   const TopoDS_Face&  face,
-                                                  SMESH_MesherHelper& helper,
+                                                  SMESH_MesherHelper& /*helper*/,
                                                   gp_Dir&             normal )
 {
   BRepAdaptor_Surface surface( face );
@@ -4073,7 +4296,7 @@ void _OffsetPlane::ComputeIntersectionLine( _OffsetPlane&        pln,
     // 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 );
+    lineDir  = getEdgeDir( E, V, 0.1 * SMESH_Algo::EdgeLength( E ));
   }
   else
   {
@@ -4324,12 +4547,23 @@ gp_XYZ _LayerEdge::Copy( _LayerEdge&         other,
  */
 //================================================================================
 
-void _LayerEdge::SetCosin( double cosin )
+double _LayerEdge::SetCosin( double cosin )
 {
   _cosin = cosin;
   cosin = Abs( _cosin );
-  //_lenFactor = ( cosin < 1.-1e-12 ) ?  Min( 2., 1./sqrt(1-cosin*cosin )) : 1.0;
-  _lenFactor = ( cosin < 1.-1e-12 ) ?  1./sqrt(1-cosin*cosin ) : 1.0;
+  //_lenFactor = ( cosin < 1.-1e-12 ) ?  1./sqrt(1-cosin*cosin ) : 1.0;
+  double realLenFactor;
+  if ( cosin < 1.-1e-12 )
+  {
+    _lenFactor = realLenFactor = 1./sqrt(1-cosin*cosin );
+  }
+  else
+  {
+    _lenFactor = 1;
+    realLenFactor = Precision::Infinite();
+  }
+
+  return realLenFactor;
 }
 
 //================================================================================
@@ -4510,7 +4744,7 @@ void _ViscousBuilder::computeGeomSize( _SolidData& data )
     double thinkness = eos._hyp.GetTotalThickness();
     for ( size_t i = 0; i < eos._edges.size(); ++i )
     {
-      if ( eos._edges[i]->Is( _LayerEdge::BLOCKED )) continue;
+      if ( eos._edges[i]->_nodes.size() < 2 ) continue;
       eos._edges[i]->SetMaxLen( thinkness );
       eos._edges[i]->FindIntersection( *searcher, intersecDist, data._epsilon, eos, &face );
       if ( intersecDist > 0 && face )
@@ -4701,7 +4935,7 @@ bool _ViscousBuilder::inflate(_SolidData& data)
         if ( eos._edges[i]->_nodes.size() > 1 )
           avgThick    += Min( 1., eos._edges[i]->_len / shapeTgtThick );
         else
-          avgThick    += shapeTgtThick;
+          avgThick    += 1;
         nbActiveEdges += ( ! eos._edges[i]->Is( _LayerEdge::BLOCKED ));
       }
     }
@@ -4882,24 +5116,41 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data,
         for ( size_t iEOS = 0; iEOS < eosC1.size(); ++iEOS )
         {
           isConcaveFace[ iEOS ] = data._concaveFaces.count( eosC1[ iEOS ]->_shapeID  );
-          vector< _LayerEdge* > & edges = eosC1[ iEOS ]->_edges;
-          for ( size_t i = 0; i < edges.size(); ++i )
-            if ( edges[i]->Is( _LayerEdge::MOVED ) ||
-                 edges[i]->Is( _LayerEdge::NEAR_BOUNDARY ))
-              movedEdges.push_back( edges[i] );
 
+          if ( eosC1[ iEOS ]->_mapper2D )
+          {
+            // compute node position by boundary node position in structured mesh
+            dumpFunction(SMESH_Comment("map2dS")<<data._index<<"_Fa"<<eos._shapeID
+                         <<"_InfStep"<<infStep);
+
+            eosC1[ iEOS ]->_mapper2D->ComputeNodePositions();
+
+            for ( _LayerEdge* le : eosC1[ iEOS ]->_edges )
+              le->_pos.back() = SMESH_NodeXYZ( le->_nodes.back() );
+
+            dumpFunctionEnd();
+          }
+          else
+          {
+            for ( _LayerEdge* le : eosC1[ iEOS ]->_edges )
+              if ( le->Is( _LayerEdge::MOVED ) ||
+                   le->Is( _LayerEdge::NEAR_BOUNDARY ))
+                movedEdges.push_back( le );
+          }
           makeOffsetSurface( *eosC1[ iEOS ], helper );
         }
 
         int step = 0, stepLimit = 5, nbBad = 0;
         while (( ++step <= stepLimit ) || improved )
         {
-          dumpFunction(SMESH_Comment("smooth")<<data._index<<"_Fa"<<sInd
-                       <<"_InfStep"<<infStep<<"_"<<step); // debug
           int oldBadNb = nbBad;
           badEdges.clear();
 
 #ifdef INCREMENTAL_SMOOTH
+          // smooth moved only
+          if ( !movedEdges.empty() )
+            dumpFunction(SMESH_Comment("smooth")<<data._index<<"_Fa"<<sInd
+                         <<"_InfStep"<<infStep<<"_"<<step); // debug
           bool findBest = false; // ( step == stepLimit );
           for ( size_t i = 0; i < movedEdges.size(); ++i )
           {
@@ -4908,9 +5159,14 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data,
               badEdges.push_back( movedEdges[i] );
           }
 #else
+          // smooth all
+          dumpFunction(SMESH_Comment("smooth")<<data._index<<"_Fa"<<sInd
+                       <<"_InfStep"<<infStep<<"_"<<step); // debug
           bool findBest = ( step == stepLimit || isConcaveFace[ iEOS ]);
           for ( size_t iEOS = 0; iEOS < eosC1.size(); ++iEOS )
           {
+            if ( eosC1[ iEOS ]->_mapper2D )
+              continue;
             vector< _LayerEdge* > & edges = eosC1[ iEOS ]->_edges;
             for ( size_t i = 0; i < edges.size(); ++i )
             {
@@ -4981,11 +5237,11 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data,
 
         } // smoothing steps
 
-        // project -- to prevent intersections or fix bad simplices
+        // project -- to prevent intersections or to fix bad simplices
         for ( size_t iEOS = 0; iEOS < eosC1.size(); ++iEOS )
         {
           if ( ! eosC1[ iEOS ]->_eosConcaVer.empty() || nbBad > 0 )
-            putOnOffsetSurface( *eosC1[ iEOS ], infStep, eosC1 );
+            putOnOffsetSurface( *eosC1[ iEOS ], -infStep, eosC1 );
         }
 
         //if ( !badEdges.empty() )
@@ -5157,7 +5413,7 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data,
           continue;
 
         // ignore intersection with intFace of an adjacent FACE
-        if ( dist > 0.1 * eos._edges[i]->_len )
+        if ( dist > 0.01 * eos._edges[i]->_len )
         {
           bool toIgnore = false;
           if (  eos._toSmooth )
@@ -5193,8 +5449,11 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data,
 
         // intersection not ignored
 
-        if ( toBlockInfaltion &&
-             dist < ( eos._edges[i]->_len * theThickToIntersection ))
+        double minDist = 0;
+        if ( eos._edges[i]->_maxLen < 0.99 * eos._hyp.GetTotalThickness() ) // limited length
+          minDist = eos._edges[i]->_len * theThickToIntersection;
+
+        if ( toBlockInfaltion && dist < minDist  )
         {
           if ( is1stBlocked ) { is1stBlocked = false; // debug
             dumpFunction(SMESH_Comment("blockIntersected") <<data._index<<"_InfStep"<<infStep);
@@ -5232,7 +5491,9 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data,
   } // loop on data._edgesOnShape
 
   if ( !is1stBlocked )
+  {
     dumpFunctionEnd();
+  }
 
   if ( closestFace && le )
   {
@@ -5426,14 +5687,14 @@ 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() );
-  double offset = baseSurface->Gap();
+  eos._offsetValue = baseSurface->Gap();
 
   eos._offsetSurf.Nullify();
 
   try
   {
     BRepOffsetAPI_MakeOffsetShape offsetMaker;
-    offsetMaker.PerformByJoin( eos._shape, -offset, Precision::Confusion() );
+    offsetMaker.PerformByJoin( eos._shape, -eos._offsetValue, Precision::Confusion() );
     if ( !offsetMaker.IsDone() ) return;
 
     TopExp_Explorer fExp( offsetMaker.Shape(), TopAbs_FACE );
@@ -5445,7 +5706,7 @@ void _ViscousBuilder::makeOffsetSurface( _EdgesOnShape& eos, SMESH_MesherHelper&
 
     eos._offsetSurf = new ShapeAnalysis_Surface( surf );
   }
-  catch ( Standard_Failure )
+  catch ( Standard_Failure& )
   {
   }
 }
@@ -5485,6 +5746,7 @@ void _ViscousBuilder::putOnOffsetSurface( _EdgesOnShape&            eos,
     return;
 
   double preci = BRep_Tool::Tolerance( TopoDS::Face( eof->_shape )), vol;
+  bool neighborHasRiskySWOL = false;
   for ( size_t i = 0; i < eos._edges.size(); ++i )
   {
     _LayerEdge* edge = eos._edges[i];
@@ -5496,12 +5758,23 @@ void _ViscousBuilder::putOnOffsetSurface( _EdgesOnShape&            eos,
       if ( !edge->Is( _LayerEdge::UPD_NORMAL_CONV ))
         continue;
     }
+    else if ( moveAll == _LayerEdge::RISKY_SWOL )
+    {
+      if ( !edge->Is( _LayerEdge::RISKY_SWOL ) ||
+           edge->_cosin < 0 )
+        continue;
+    }
     else if ( !moveAll && !edge->Is( _LayerEdge::MOVED ))
       continue;
 
     int nbBlockedAround = 0;
     for ( size_t iN = 0; iN < edge->_neibors.size(); ++iN )
+    {
       nbBlockedAround += edge->_neibors[iN]->Is( _LayerEdge::BLOCKED );
+      if ( edge->_neibors[iN]->Is( _LayerEdge::RISKY_SWOL ) &&
+           edge->_neibors[iN]->_cosin > 0 )
+        neighborHasRiskySWOL = true;
+    }
     if ( nbBlockedAround > 1 )
       continue;
 
@@ -5530,6 +5803,35 @@ void _ViscousBuilder::putOnOffsetSurface( _EdgesOnShape&            eos,
       {
         edge->_normal = ( newP - prevP ).Normalized();
       }
+      // if ( edge->_len < eof->_offsetValue )
+      //   edge->_len = eof->_offsetValue;
+
+      if ( !eos._sWOL.IsNull() ) // RISKY_SWOL
+      {
+        double change = eof->_offsetSurf->Gap() / eof->_offsetValue;
+        if (( newP - tgtP.XYZ() ) * edge->_normal < 0 )
+          change = 1 - change;
+        else
+          change = 1 + change;
+        gp_XYZ shitfVec    = tgtP.XYZ() - SMESH_NodeXYZ( edge->_nodes[0] );
+        gp_XYZ newShiftVec = shitfVec * change;
+        double shift       = edge->_normal * shitfVec;
+        double newShift    = edge->_normal * newShiftVec;
+        newP = tgtP.XYZ() + edge->_normal * ( newShift - shift );
+
+        uv = eof->_offsetSurf->NextValueOfUV( edge->_curvature->_uv, newP, preci );
+        if ( eof->_offsetSurf->Gap() < edge->_len )
+        {
+          edge->_curvature->_uv = uv;
+          newP = eof->_offsetSurf->Value( uv ).XYZ();
+        }
+        n->setXYZ( newP.X(), newP.Y(), newP.Z());
+        if ( !edge->UpdatePositionOnSWOL( n, /*tol=*/10 * edge->_len / ( edge->NbSteps() + 1 ),
+                                          eos, eos.GetData().GetHelper() ))
+        {
+          debugMsg("UpdatePositionOnSWOL fails in putOnOffsetSurface()" );
+        }
+      }
     }
   }
 
@@ -5543,12 +5845,13 @@ void _ViscousBuilder::putOnOffsetSurface( _EdgesOnShape&            eos,
       break;
   if ( i < eos._edges.size() )
   {
-    dumpFunction(SMESH_Comment("putOnOffsetSurface_S") << eos._shapeID
-                 << "_InfStep" << infStep << "_" << smooStep );
+    dumpFunction(SMESH_Comment("putOnOffsetSurface_") << eos.ShapeTypeLetter() << eos._shapeID
+                 << "_InfStep" << infStep << "_" << Abs( smooStep ));
     for ( ; i < eos._edges.size(); ++i )
     {
-      if ( eos._edges[i]->Is( _LayerEdge::MARKED ))
+      if ( eos._edges[i]->Is( _LayerEdge::MARKED )) {
         dumpMove( eos._edges[i]->_nodes.back() );
+      }
     }
     dumpFunctionEnd();
   }
@@ -5564,7 +5867,7 @@ void _ViscousBuilder::putOnOffsetSurface( _EdgesOnShape&            eos,
     SMESH_subMeshIteratorPtr smIt = eos._subMesh->getDependsOnIterator(/*includeSelf=*/false);
     while ( smIt->more() )
     {
-      SMESH_subMesh* sm = smIt->next();
+      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;
@@ -5573,6 +5876,28 @@ void _ViscousBuilder::putOnOffsetSurface( _EdgesOnShape&            eos,
     }
     cnvFace->_normalsFixedOnBorders = true;
   }
+
+
+  // bos #20643
+  // negative smooStep means "final step", where we don't treat RISKY_SWOL edges
+  // as edges based on FACE are a bit late comparing with them
+  if ( smooStep >= 0 &&
+       neighborHasRiskySWOL &&
+       moveAll != _LayerEdge::RISKY_SWOL &&
+       eos.ShapeType() == TopAbs_FACE )
+  {
+    // 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::RISKY_SWOL );
+    }
+  }
 }
 
 //================================================================================
@@ -5875,9 +6200,11 @@ bool _Smoother1D::smoothAnalyticEdge( _SolidData&                    data,
           tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() );
           dumpMove( tgtNode );
 
-          SMDS_FacePositionPtr pos = tgtNode->GetPosition();
-          pos->SetUParameter( newUV.X() );
-          pos->SetVParameter( newUV.Y() );
+          if ( SMDS_FacePositionPtr pos = tgtNode->GetPosition() ) // NULL if F is noShrink
+          {
+            pos->SetUParameter( newUV.X() );
+            pos->SetVParameter( newUV.Y() );
+          }
 
           gp_XYZ newUV0( newUV.X(), newUV.Y(), 0 );
 
@@ -5987,10 +6314,11 @@ bool _Smoother1D::smoothAnalyticEdge( _SolidData&                    data,
         tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() );
         dumpMove( tgtNode );
 
-        SMDS_FacePositionPtr pos = tgtNode->GetPosition();
-        pos->SetUParameter( newUV.X() );
-        pos->SetVParameter( newUV.Y() );
-
+        if ( SMDS_FacePositionPtr pos = tgtNode->GetPosition() ) // NULL if F is noShrink
+        {
+          pos->SetUParameter( newUV.X() );
+          pos->SetVParameter( newUV.Y() );
+        }
         _eos[i]->Set( _LayerEdge::SMOOTHED ); // to check in refine() (IPAL54237)
       }
     }
@@ -6006,10 +6334,10 @@ bool _Smoother1D::smoothAnalyticEdge( _SolidData&                    data,
  */
 //================================================================================
 
-bool _Smoother1D::smoothComplexEdge( _SolidData&                    data,
+bool _Smoother1D::smoothComplexEdge( _SolidData&                    /*data*/,
                                      Handle(ShapeAnalysis_Surface)& surface,
                                      const TopoDS_Face&             F,
-                                     SMESH_MesherHelper&            helper)
+                                     SMESH_MesherHelper&            /*helper*/)
 {
   if ( _offPoints.empty() )
     return false;
@@ -6238,6 +6566,7 @@ void _Smoother1D::prepare(_SolidData& data)
   _edgeDir[1] = getEdgeDir( E, leOnV[1]->_nodes[0], data.GetHelper() );
   _leOnV[0]._cosin = Abs( leOnV[0]->_cosin );
   _leOnV[1]._cosin = Abs( leOnV[1]->_cosin );
+  _leOnV[0]._flags = _leOnV[1]._flags = 0;
   if ( _eos._sWOL.IsNull() ) // 3D
     for ( int iEnd = 0; iEnd < 2; ++iEnd )
       _leOnV[iEnd]._cosin = Abs( _edgeDir[iEnd].Normalized() * leOnV[iEnd]->_normal );
@@ -6247,7 +6576,7 @@ void _Smoother1D::prepare(_SolidData& data)
 
   // divide E to have offset segments with low deflection
   BRepAdaptor_Curve c3dAdaptor( E );
-  const double curDeflect = 0.1; //0.01; // Curvature deflection == |p1p2]*sin(p1p2,p1pM)
+  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 )
@@ -6385,6 +6714,16 @@ gp_XYZ _Smoother1D::getNormalNormal( const gp_XYZ & normal,
   // if ( size == 0 ) // MULTI_NORMAL _LayerEdge
   //   return gp_XYZ( 1e-100, 1e-100, 1e-100 );
 
+  if ( size < 1e-5 ) // normal || edgeDir (almost) at inflation along EDGE (bos #20643)
+  {
+    const _LayerEdge* le = _eos._edges[ _eos._edges.size() / 2 ];
+    const gp_XYZ& leNorm = le->_normal;
+
+    cross = leNorm ^ edgeDir;
+    norm = edgeDir ^ cross;
+    size = norm.Modulus();
+  }
+
   return norm / size;
 }
 
@@ -6397,7 +6736,7 @@ gp_XYZ _Smoother1D::getNormalNormal( const gp_XYZ & normal,
 void _Smoother1D::offPointsToPython() const
 {
   const char* fname = "/tmp/offPoints.py";
-  cout << "execfile('"<<fname<<"')"<<endl;
+  cout << "exec(open('"<<fname<<"','rb').read() )"<<endl;
   ofstream py(fname);
   py << "import SMESH" << endl
      << "from salome.smesh import smeshBuilder" << endl
@@ -6612,20 +6951,98 @@ void _SolidData::PrepareEdgesToSmoothOnFace( _EdgesOnShape* eos, bool substitute
     eos->_edgeForOffset = 0;
 
     double maxCosin = -1;
+    //bool hasNoShrink = false;
     for ( TopExp_Explorer eExp( eos->_shape, TopAbs_EDGE ); eExp.More(); eExp.Next() )
     {
       _EdgesOnShape* eoe = GetShapeEdges( eExp.Current() );
       if ( !eoe || eoe->_edges.empty() ) continue;
 
+      // if ( eos->GetData()._noShrinkShapes.count( eoe->_shapeID ))
+      //   hasNoShrink = true;
+
       vector<_LayerEdge*>& eE = eoe->_edges;
       _LayerEdge* e = eE[ eE.size() / 2 ];
-      if ( e->_cosin > maxCosin )
+      if ( !e->Is( _LayerEdge::RISKY_SWOL ) && e->_cosin > maxCosin )
       {
         eos->_edgeForOffset = e;
         maxCosin = e->_cosin;
       }
+
+      if ( !eoe->_sWOL.IsNull() )
+        for ( _LayerEdge* le : eoe->_edges )
+          if ( le->Is( _LayerEdge::RISKY_SWOL ) && e->_cosin > 0 )
+          {
+            // make _neibors on FACE be smoothed after le->Is( BLOCKED )
+            for ( _LayerEdge* neibor : le->_neibors )
+            {
+              int shapeDim =  neibor->BaseShapeDim();
+              if ( shapeDim == 2 )
+                neibor->Set( _LayerEdge::NEAR_BOUNDARY ); // on FACE
+              else if ( shapeDim == 0 )
+                neibor->Set( _LayerEdge::RISKY_SWOL );    // on VERTEX
+
+              if ( !neibor->_curvature )
+              {
+                gp_XY uv = helper.GetNodeUV( F, neibor->_nodes[0] );
+                neibor->_curvature = _Factory::NewCurvature();
+                neibor->_curvature->_r = 0;
+                neibor->_curvature->_k = 0;
+                neibor->_curvature->_h2lenRatio = 0;
+                neibor->_curvature->_uv = uv;
+              }
+            }
+          }
+    } // loop on EDGEs
+
+    // Try to initialize _Mapper2D
+
+    // if ( hasNoShrink )
+    //   return;
+
+    SMDS_ElemIteratorPtr fIt = eos->_subMesh->GetSubMeshDS()->GetElements();
+    if ( !fIt->more() || fIt->next()->NbCornerNodes() != 4 )
+      return;
+
+    // get EDGEs of quadrangle bottom
+    std::list< TopoDS_Edge > edges;
+    std::list< int > nbEdgesInWire;
+    int nbWire = SMESH_Block::GetOrderedEdges( F, edges, nbEdgesInWire );
+    if ( nbWire != 1 || nbEdgesInWire.front() < 4 )
+      return;
+    const SMDS_MeshNode* node;
+    while ( true ) // make edges start at a corner VERTEX
+    {
+      node = SMESH_Algo::VertexNode( helper.IthVertex( 0, edges.front() ), helper.GetMeshDS() );
+      if ( node && helper.IsCornerOfStructure( node, eos->_subMesh->GetSubMeshDS(), helper ))
+        break;
+      edges.pop_front();
+      if ( edges.empty() )
+        return;
     }
-  }
+    std::list< TopoDS_Edge >::iterator edgeIt = edges.begin();
+    while ( true ) // make edges finish at a corner VERTEX
+    {
+      node = SMESH_Algo::VertexNode( helper.IthVertex( 1, *edgeIt ), helper.GetMeshDS() );
+      ++edgeIt;
+      if ( node && helper.IsCornerOfStructure( node, eos->_subMesh->GetSubMeshDS(), helper ))
+      {
+        edges.erase( edgeIt, edges.end() );
+        break;
+      }
+      if ( edgeIt == edges.end() )
+        return;
+    }
+
+    // get structure of nodes
+    TParam2ColumnMap param2ColumnMap;
+    if ( !helper.LoadNodeColumns( param2ColumnMap, F, edges, helper.GetMeshDS() ))
+      return;
+
+    eos->_mapper2D = new _Mapper2D( param2ColumnMap, eos->GetData()._n2eMap );
+
+  } // if eos is of curved FACE
+
+  return;
 }
 
 //================================================================================
@@ -6666,7 +7083,7 @@ void _SolidData::AddShapesToSmooth( const set< _EdgesOnShape* >& eosToSmooth,
  */
 //================================================================================
 
-void _ViscousBuilder::limitMaxLenByCurvature( _SolidData& data, SMESH_MesherHelper& helper )
+void _ViscousBuilder::limitMaxLenByCurvature( _SolidData& data, SMESH_MesherHelper& /*helper*/ )
 {
   // find intersection of neighbor _LayerEdge's to limit _maxLen
   // according to local curvature (IPAL52648)
@@ -6718,9 +7135,9 @@ void _ViscousBuilder::limitMaxLenByCurvature( _SolidData& data, SMESH_MesherHelp
 
 void _ViscousBuilder::limitMaxLenByCurvature( _LayerEdge*    e1,
                                               _LayerEdge*    e2,
-                                              _EdgesOnShape& eos1,
-                                              _EdgesOnShape& eos2,
-                                              const bool     isSmoothable )
+                                              _EdgesOnShape& /*eos1*/,
+                                              _EdgesOnShape& /*eos2*/,
+                                              const bool     /*isSmoothable*/ )
 {
   if (( e1->_nodes[0]->GetPosition()->GetDim() !=
         e2->_nodes[0]->GetPosition()->GetDim() ) &&
@@ -7039,7 +7456,7 @@ void _ViscousBuilder::findEdgesToUpdateNormalNearConvexFace( _ConvexFace &
 bool _ViscousBuilder::updateNormals( _SolidData&         data,
                                      SMESH_MesherHelper& helper,
                                      int                 stepNb,
-                                     double              stepSize)
+                                     double              /*stepSize*/)
 {
   updateNormalsOfC1Vertices( data );
 
@@ -7204,7 +7621,7 @@ 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 ))
+      if ( edge->Is( _LayerEdge::BLOCKED ) && newEdge._maxLen > edge->_len )
         continue;
 
       // Check if a new _normal is OK:
@@ -7237,7 +7654,8 @@ bool _ViscousBuilder::updateNormals( _SolidData&         data,
           PShapeIteratorPtr eIt = helper.GetAncestors( V, *_mesh, TopAbs_EDGE, &eos->_sWOL );
           while ( const TopoDS_Shape* E = eIt->next() )
           {
-            gp_Vec edgeDir = getEdgeDir( TopoDS::Edge( *E ), TopoDS::Vertex( V ));
+            gp_Vec edgeDir = getEdgeDir( TopoDS::Edge( *E ), TopoDS::Vertex( V ),
+                                         eos->_hyp.Get1stLayerThickness() );
             double   angle = edgeDir.Angle( newEdge._normal ); // [0,PI]
             if ( angle < M_PI / 2 )
               shapesToSmooth.insert( data.GetShapeEdges( *E ));
@@ -7379,7 +7797,7 @@ bool _ViscousBuilder::isNewNormalOk( _SolidData&   data,
 //================================================================================
 
 bool _ViscousBuilder::updateNormalsOfSmoothed( _SolidData&         data,
-                                               SMESH_MesherHelper& helper,
+                                               SMESH_MesherHelper& /*helper*/,
                                                const int           nbSteps,
                                                const double        stepSize )
 {
@@ -8837,7 +9255,7 @@ int _LayerEdge::Smooth(const int step, const bool isConcaveFace, bool findBest )
 //================================================================================
 
 void _LayerEdge::ChooseSmooFunction( const set< TGeomID >& concaveVertices,
-                                     const TNode2Edge&     n2eMap)
+                                     const TNode2Edge&     /*n2eMap*/)
 {
   if ( _smooFunction ) return;
 
@@ -9030,7 +9448,7 @@ gp_XYZ _LayerEdge::smoothAngular()
       else
         norm += cross;
     }
-    catch (Standard_Failure) { // if |cross| == 0.
+    catch (Standard_Failure&) { // if |cross| == 0.
     }
   }
   gp_XYZ vec = newPos - pN;
@@ -9523,6 +9941,8 @@ void _LayerEdge::SetNewLength( double len, _EdgesOnShape& eos, SMESH_MesherHelpe
     Block( eos.GetData() );
   }
   const double lenDelta = len - _len;
+  // if ( lenDelta < 0 )
+  //   return;
   if ( lenDelta < len * 1e-3  )
   {
     Block( eos.GetData() );
@@ -9566,46 +9986,13 @@ void _LayerEdge::SetNewLength( double len, _EdgesOnShape& eos, SMESH_MesherHelpe
   _pos.push_back( newXYZ );
 
   if ( !eos._sWOL.IsNull() )
-  {
-    double distXYZ[4];
-    bool uvOK = false;
-    if ( eos.SWOLType() == TopAbs_EDGE )
-    {
-      double u = Precision::Infinite(); // to force projection w/o distance check
-      uvOK = helper.CheckNodeU( TopoDS::Edge( eos._sWOL ), n, u,
-                                /*tol=*/2*lenDelta, /*force=*/true, distXYZ );
-      _pos.back().SetCoord( u, 0, 0 );
-      if ( _nodes.size() > 1 && uvOK )
-      {
-        SMDS_EdgePositionPtr pos = n->GetPosition();
-        pos->SetUParameter( u );
-      }
-    }
-    else //  TopAbs_FACE
-    {
-      gp_XY uv( Precision::Infinite(), 0 );
-      uvOK = helper.CheckNodeUV( TopoDS::Face( eos._sWOL ), n, uv,
-                                 /*tol=*/2*lenDelta, /*force=*/true, distXYZ );
-      _pos.back().SetCoord( uv.X(), uv.Y(), 0 );
-      if ( _nodes.size() > 1 && uvOK )
-      {
-        SMDS_FacePositionPtr pos = n->GetPosition();
-        pos->SetUParameter( uv.X() );
-        pos->SetVParameter( uv.Y() );
-      }
-    }
-    if ( uvOK )
-    {
-      n->setXYZ( distXYZ[1], distXYZ[2], distXYZ[3]);
-    }
-    else
+    if ( !UpdatePositionOnSWOL( n, 2*lenDelta, eos, helper ))
     {
       n->setXYZ( oldXYZ.X(), oldXYZ.Y(), oldXYZ.Z() );
       _pos.pop_back();
       Block( eos.GetData() );
       return;
     }
-  }
 
   _len = len;
 
@@ -9614,13 +10001,57 @@ void _LayerEdge::SetNewLength( double len, _EdgesOnShape& eos, SMESH_MesherHelpe
   {
     for ( size_t i = 0; i < _neibors.size(); ++i )
       //if (  _len > _neibors[i]->GetSmooLen() )
-        _neibors[i]->Set( MOVED );
+      _neibors[i]->Set( MOVED );
 
     Set( MOVED );
   }
   dumpMove( n ); //debug
 }
 
+
+//================================================================================
+/*!
+ * \brief Update last position on SWOL by projecting node on SWOL
+*/
+//================================================================================
+
+bool _LayerEdge::UpdatePositionOnSWOL( SMDS_MeshNode*      n,
+                                       double              tol,
+                                       _EdgesOnShape&      eos,
+                                       SMESH_MesherHelper& helper )
+{
+  double distXYZ[4];
+  bool uvOK = false;
+  if ( eos.SWOLType() == TopAbs_EDGE )
+  {
+    double u = Precision::Infinite(); // to force projection w/o distance check
+    uvOK = helper.CheckNodeU( TopoDS::Edge( eos._sWOL ), n, u, tol, /*force=*/true, distXYZ );
+    _pos.back().SetCoord( u, 0, 0 );
+    if ( _nodes.size() > 1 && uvOK )
+    {
+      SMDS_EdgePositionPtr pos = n->GetPosition();
+      pos->SetUParameter( u );
+    }
+  }
+  else //  TopAbs_FACE
+  {
+    gp_XY uv( Precision::Infinite(), 0 );
+    uvOK = helper.CheckNodeUV( TopoDS::Face( eos._sWOL ), n, uv, tol, /*force=*/true, distXYZ );
+    _pos.back().SetCoord( uv.X(), uv.Y(), 0 );
+    if ( _nodes.size() > 1 && uvOK )
+    {
+      SMDS_FacePositionPtr pos = n->GetPosition();
+      pos->SetUParameter( uv.X() );
+      pos->SetVParameter( uv.Y() );
+    }
+  }
+  if ( uvOK )
+  {
+    n->setXYZ( distXYZ[1], distXYZ[2], distXYZ[3]);
+  }
+  return uvOK;
+}
+
 //================================================================================
 /*!
  * \brief Set BLOCKED flag and propagate limited _maxLen to _neibors
@@ -10037,11 +10468,16 @@ bool _ViscousBuilder::refine(_SolidData& data)
             segLen[j] = segLen[j-1] + (edge._pos[j-1] - edge._pos[j] ).Modulus();
         }
       }
-      else if ( !surface.IsNull() ) // SWOL surface with singularities
+      else // SWOL is surface with singularities or irregularly parametrized curve
       {
         pos3D.resize( edge._pos.size() );
-        for ( size_t j = 0; j < edge._pos.size(); ++j )
-          pos3D[j] = surface->Value( edge._pos[j].X(), edge._pos[j].Y() ).XYZ();
+
+        if ( !surface.IsNull() )
+          for ( size_t j = 0; j < edge._pos.size(); ++j )
+            pos3D[j] = surface->Value( edge._pos[j].X(), edge._pos[j].Y() ).XYZ();
+        else if ( !curve.IsNull() )
+          for ( size_t j = 0; j < edge._pos.size(); ++j )
+            pos3D[j] = curve->Value( edge._pos[j].X() ).XYZ();
 
         for ( size_t j = 1; j < edge._pos.size(); ++j )
           segLen[j] = segLen[j-1] + ( pos3D[j-1] - pos3D[j] ).Modulus();
@@ -10073,7 +10509,8 @@ bool _ViscousBuilder::refine(_SolidData& data)
       if ( n2eMap && (( n2e = n2eMap->find( edge._nodes[0] )) != n2eMap->end() ))
       {
         edgeOnSameNode = n2e->second;
-        useExistingPos = ( edgeOnSameNode->_len < edge._len );
+        useExistingPos = ( edgeOnSameNode->_len < edge._len ||
+                           segLen[0] == segLen.back() ); // too short inflation step (bos #20643)
         const gp_XYZ& otherTgtPos = edgeOnSameNode->_pos.back();
         SMDS_PositionPtr  lastPos = tgtNode->GetPosition();
         if ( isOnEdge )
@@ -10088,26 +10525,16 @@ bool _ViscousBuilder::refine(_SolidData& data)
           fpos->SetVParameter( otherTgtPos.Y() );
         }
       }
-      // calculate height of the first layer
-      double h0;
-      const double T = segLen.back(); //data._hyp.GetTotalThickness();
-      const double f = eos._hyp.GetStretchFactor();
-      const int    N = eos._hyp.GetNumberLayers();
-      const double fPowN = pow( f, N );
-      if ( fPowN - 1 <= numeric_limits<double>::min() )
-        h0 = T / N;
-      else
-        h0 = T * ( f - 1 )/( fPowN - 1 );
-
-      const double zeroLen = std::numeric_limits<double>::min();
 
       // create intermediate nodes
-      double hSum = 0, hi = h0/f;
+      const double      h0 = eos._hyp.Get1stLayerThickness( segLen.back() );
+      const double zeroLen = std::numeric_limits<double>::min();
+      double hSum = 0, hi = h0/eos._hyp.GetStretchFactor();
       size_t iSeg = 1;
       for ( size_t iStep = 1; iStep < edge._nodes.size(); ++iStep )
       {
         // compute an intermediate position
-        hi *= f;
+        hi *= eos._hyp.GetStretchFactor();
         hSum += hi;
         while ( hSum > segLen[iSeg] && iSeg < segLen.size()-1 )
           ++iSeg;
@@ -10230,6 +10657,11 @@ bool _ViscousBuilder::refine(_SolidData& data)
     const TGeomID faceID = getMeshDS()->ShapeToIndex( exp.Current() );
     if ( data._ignoreFaceIds.count( faceID ))
       continue;
+    _EdgesOnShape*    eos = data.GetShapeEdges( faceID );
+    SMDS_MeshGroup* group = StdMeshers_ViscousLayers::CreateGroup( eos->_hyp.GetGroupName(),
+                                                                   *helper.GetMesh(),
+                                                                   SMDSAbs_Volume );
+    std::vector< const SMDS_MeshElement* > vols;
     const bool isReversedFace = data._reversedFaceIds.count( faceID );
     SMESHDS_SubMesh*    fSubM = getMeshDS()->MeshElements( exp.Current() );
     SMDS_ElemIteratorPtr  fIt = fSubM->GetElements();
@@ -10260,14 +10692,20 @@ bool _ViscousBuilder::refine(_SolidData& data)
       if ( 0 < nnSet.size() && nnSet.size() < 3 )
         continue;
 
+      vols.clear();
+      const SMDS_MeshElement* vol;
+
       switch ( nbNodes )
       {
       case 3: // TRIA
       {
         // PENTA
         for ( size_t iZ = 1; iZ < minZ; ++iZ )
-          helper.AddVolume( (*nnVec[0])[iZ-1], (*nnVec[1])[iZ-1], (*nnVec[2])[iZ-1],
-                            (*nnVec[0])[iZ],   (*nnVec[1])[iZ],   (*nnVec[2])[iZ]);
+        {
+          vol = helper.AddVolume( (*nnVec[0])[iZ-1], (*nnVec[1])[iZ-1], (*nnVec[2])[iZ-1],
+                                  (*nnVec[0])[iZ],   (*nnVec[1])[iZ],   (*nnVec[2])[iZ]);
+          vols.push_back( vol );
+        }
 
         for ( size_t iZ = minZ; iZ < maxZ; ++iZ )
         {
@@ -10280,16 +10718,18 @@ bool _ViscousBuilder::refine(_SolidData& data)
             int i2 = *degenEdgeInd.begin();
             int i0 = helper.WrapIndex( i2 - 1, nbNodes );
             int i1 = helper.WrapIndex( i2 + 1, nbNodes );
-            helper.AddVolume( (*nnVec[i0])[iZ-1], (*nnVec[i1])[iZ-1],
-                              (*nnVec[i1])[iZ  ], (*nnVec[i0])[iZ  ], (*nnVec[i2]).back());
+            vol = helper.AddVolume( (*nnVec[i0])[iZ-1], (*nnVec[i1])[iZ-1],
+                                    (*nnVec[i1])[iZ  ], (*nnVec[i0])[iZ  ], (*nnVec[i2]).back());
+            vols.push_back( vol );
           }
           else  // TETRA
           {
             int i3 = !degenEdgeInd.count(0) ? 0 : !degenEdgeInd.count(1) ? 1 : 2;
-            helper.AddVolume( (*nnVec[  0 ])[ i3 == 0 ? iZ-1 : nnVec[0]->size()-1 ],
-                              (*nnVec[  1 ])[ i3 == 1 ? iZ-1 : nnVec[1]->size()-1 ],
-                              (*nnVec[  2 ])[ i3 == 2 ? iZ-1 : nnVec[2]->size()-1 ],
-                              (*nnVec[ i3 ])[ iZ ]);
+            vol = helper.AddVolume( (*nnVec[  0 ])[ i3 == 0 ? iZ-1 : nnVec[0]->size()-1 ],
+                                    (*nnVec[  1 ])[ i3 == 1 ? iZ-1 : nnVec[1]->size()-1 ],
+                                    (*nnVec[  2 ])[ i3 == 2 ? iZ-1 : nnVec[2]->size()-1 ],
+                                    (*nnVec[ i3 ])[ iZ ]);
+            vols.push_back( vol );
           }
         }
         break; // TRIA
@@ -10298,10 +10738,13 @@ bool _ViscousBuilder::refine(_SolidData& data)
       {
         // HEX
         for ( size_t iZ = 1; iZ < minZ; ++iZ )
-          helper.AddVolume( (*nnVec[0])[iZ-1], (*nnVec[1])[iZ-1],
-                            (*nnVec[2])[iZ-1], (*nnVec[3])[iZ-1],
-                            (*nnVec[0])[iZ],   (*nnVec[1])[iZ],
-                            (*nnVec[2])[iZ],   (*nnVec[3])[iZ]);
+        {
+          vol = helper.AddVolume( (*nnVec[0])[iZ-1], (*nnVec[1])[iZ-1],
+                                  (*nnVec[2])[iZ-1], (*nnVec[3])[iZ-1],
+                                  (*nnVec[0])[iZ],   (*nnVec[1])[iZ],
+                                  (*nnVec[2])[iZ],   (*nnVec[3])[iZ]);
+          vols.push_back( vol );
+        }
 
         for ( size_t iZ = minZ; iZ < maxZ; ++iZ )
         {
@@ -10320,9 +10763,9 @@ bool _ViscousBuilder::refine(_SolidData& data)
             int i0 = helper.WrapIndex( i3 + 1, nbNodes );
             int i1 = helper.WrapIndex( i0 + 1, nbNodes );
 
-            const SMDS_MeshElement* vol =
-              helper.AddVolume( nnVec[i3]->back(), (*nnVec[i0])[iZ], (*nnVec[i0])[iZ-1],
-                                nnVec[i2]->back(), (*nnVec[i1])[iZ], (*nnVec[i1])[iZ-1]);
+            vol = helper.AddVolume( nnVec[i3]->back(), (*nnVec[i0])[iZ], (*nnVec[i0])[iZ-1],
+                                    nnVec[i2]->back(), (*nnVec[i1])[iZ], (*nnVec[i1])[iZ-1]);
+            vols.push_back( vol );
             if ( !ok && vol )
               degenVols.push_back( vol );
           }
@@ -10330,15 +10773,15 @@ bool _ViscousBuilder::refine(_SolidData& data)
 
           default: // degen HEX
           {
-            const SMDS_MeshElement* vol =
-              helper.AddVolume( nnVec[0]->size() > iZ-1 ? (*nnVec[0])[iZ-1] : nnVec[0]->back(),
-                                nnVec[1]->size() > iZ-1 ? (*nnVec[1])[iZ-1] : nnVec[1]->back(),
-                                nnVec[2]->size() > iZ-1 ? (*nnVec[2])[iZ-1] : nnVec[2]->back(),
-                                nnVec[3]->size() > iZ-1 ? (*nnVec[3])[iZ-1] : nnVec[3]->back(),
-                                nnVec[0]->size() > iZ   ? (*nnVec[0])[iZ]   : nnVec[0]->back(),
-                                nnVec[1]->size() > iZ   ? (*nnVec[1])[iZ]   : nnVec[1]->back(),
-                                nnVec[2]->size() > iZ   ? (*nnVec[2])[iZ]   : nnVec[2]->back(),
-                                nnVec[3]->size() > iZ   ? (*nnVec[3])[iZ]   : nnVec[3]->back());
+            vol = helper.AddVolume( nnVec[0]->size() > iZ-1 ? (*nnVec[0])[iZ-1] : nnVec[0]->back(),
+                                    nnVec[1]->size() > iZ-1 ? (*nnVec[1])[iZ-1] : nnVec[1]->back(),
+                                    nnVec[2]->size() > iZ-1 ? (*nnVec[2])[iZ-1] : nnVec[2]->back(),
+                                    nnVec[3]->size() > iZ-1 ? (*nnVec[3])[iZ-1] : nnVec[3]->back(),
+                                    nnVec[0]->size() > iZ   ? (*nnVec[0])[iZ]   : nnVec[0]->back(),
+                                    nnVec[1]->size() > iZ   ? (*nnVec[1])[iZ]   : nnVec[1]->back(),
+                                    nnVec[2]->size() > iZ   ? (*nnVec[2])[iZ]   : nnVec[2]->back(),
+                                    nnVec[3]->size() > iZ   ? (*nnVec[3])[iZ]   : nnVec[3]->back());
+            vols.push_back( vol );
             degenVols.push_back( vol );
           }
           }
@@ -10349,6 +10792,11 @@ bool _ViscousBuilder::refine(_SolidData& data)
         return error("Not supported type of element", data._index);
 
       } // switch ( nbNodes )
+
+      if ( group )
+        for ( size_t i = 0; i < vols.size(); ++i )
+          group->Add( vols[ i ]);
+
     } // while ( fIt->more() )
   } // loop on FACEs
 
@@ -10386,6 +10834,8 @@ namespace VISCOUS_3D
     PeriodicFaces( ShrinkFace* sf1, ShrinkFace* sf2 ): _shriFace{ sf1, sf2 } {}
     bool IncludeShrunk( const TopoDS_Face& face, const TopTools_MapOfShape& shrunkFaces ) const;
     bool MoveNodes( const TopoDS_Face& tgtFace );
+    void Clear() { _nnMap.clear(); }
+    bool IsEmpty() const { return _nnMap.empty(); }
   };
 
   //--------------------------------------------------------------------------------
@@ -10403,12 +10853,20 @@ namespace VISCOUS_3D
       std::vector< SMESH_NodeXYZ > _nodes;
       TopAbs_ShapeEnum             _vertSWOLType[2]; // shrink part includes VERTEXes
       AverageHyp*                  _vertHyp[2];
+      double                       _edgeWOLLen[2]; // length of wol EDGE
+      double                       _tol; // to compare _edgeWOLLen's
 
       BndPart():
         _isShrink(0), _isReverse(0), _nbSegments(0), _hyp(0),
-        _vertSWOLType{ TopAbs_WIRE, TopAbs_WIRE }, _vertHyp{ 0, 0 }
+        _vertSWOLType{ TopAbs_WIRE, TopAbs_WIRE }, _vertHyp{ 0, 0 }, _edgeWOLLen{ 0., 0.}
       {}
 
+      bool IsEqualLengthEWOL( const BndPart& other ) const
+      {
+        return ( std::abs( _edgeWOLLen[0] - other._edgeWOLLen[0] ) < _tol &&
+                 std::abs( _edgeWOLLen[1] - other._edgeWOLLen[1] ) < _tol );
+      }
+
       bool operator==( const BndPart& other ) const
       {
         return ( _isShrink       == other._isShrink &&
@@ -10419,7 +10877,8 @@ namespace VISCOUS_3D
                  (( !_isShrink ) ||
                   ( *_hyp        == *other._hyp &&
                     vertHyp1()   == other.vertHyp1() &&
-                    vertHyp2()   == other.vertHyp2() ))
+                    vertHyp2()   == other.vertHyp2() &&
+                    IsEqualLengthEWOL( other )))
                  );
       }
       bool CanAppend( const BndPart& other )
@@ -10437,10 +10896,12 @@ namespace VISCOUS_3D
         bool hasCommonNode = ( _nodes.back()->GetID() == other._nodes.front()->GetID() );
         _nodes.insert( _nodes.end(), other._nodes.begin() + hasCommonNode, other._nodes.end() );
         _vertSWOLType[1] = other._vertSWOLType[1];
-        if ( _isShrink )
-          _vertHyp[1] = other._vertHyp[1];
+        if ( _isShrink ) {
+          _vertHyp[1]    = other._vertHyp[1];
+          _edgeWOLLen[1] = other._edgeWOLLen[1];
+        }
       }
-      const SMDS_MeshNode*    Node(size_t i)  const
+      const SMDS_MeshNode* Node(size_t i)  const
       {
         return _nodes[ _isReverse ? ( _nodes.size() - 1 - i ) : i ]._node;
       }
@@ -10455,14 +10916,13 @@ namespace VISCOUS_3D
     SMESH_subMesh*       _subMesh;
     _SolidData*          _data1;
     _SolidData*          _data2;
-    //bool                 _isPeriodic;
 
     std::list< BndPart > _boundary;
     int                  _boundarySize, _nbBoundaryParts;
 
     void Init( SMESH_subMesh* sm, _SolidData* sd1, _SolidData* sd2 )
     {
-      _subMesh = sm; _data1 = sd1; _data2 = sd2; //_isPeriodic = false;
+      _subMesh = sm; _data1 = sd1; _data2 = sd2;
     }
     bool IsSame( const TopoDS_Face& face ) const
     {
@@ -10525,11 +10985,15 @@ namespace VISCOUS_3D
         if ( !periodic._trsf.Solve( srcPnts, tgtPnts )) {
           continue;
         }
-        double tol = std::numeric_limits<double>::max();
+        double tol = std::numeric_limits<double>::max(); // tolerance by segment size
         for ( size_t i = 1; i < srcPnts.size(); ++i ) {
           tol = Min( tol, ( srcPnts[i-1] - srcPnts[i] ).SquareModulus() );
         }
         tol = 0.01 * Sqrt( tol );
+        for ( BndPart& boundary : _boundary ) { // tolerance by VL thickness
+          if ( boundary._isShrink )
+            tol = Min( tol, boundary._hyp->Get1stLayerThickness() / 50. );
+        }
         bool nodeCoincide = true;
         TNodeNodeMap::iterator n2n = periodic._nnMap.begin();
         for ( ; n2n != periodic._nnMap.end() &&  nodeCoincide; ++n2n )
@@ -10537,7 +11001,7 @@ namespace VISCOUS_3D
           SMESH_NodeXYZ nSrc = n2n->first;
           SMESH_NodeXYZ nTgt = n2n->second;
           gp_XYZ pTgt = periodic._trsf.Transform( nSrc );
-          nodeCoincide = (( pTgt - nTgt ).SquareModulus() < tol );
+          nodeCoincide = (( pTgt - nTgt ).SquareModulus() < tol * tol );
         }
         if ( nodeCoincide )
           return true;
@@ -10585,6 +11049,11 @@ namespace VISCOUS_3D
       for ( int iE = 0; iE < nbEdgesInWire.front(); ++iE )
       {
         BndPart bndPart;
+
+        std::vector<const SMDS_MeshNode*> nodes = fSide.GetOrderedNodes( iE );
+        bndPart._nodes.assign( nodes.begin(), nodes.end() );
+        bndPart._nbSegments = bndPart._nodes.size() - 1;
+
         _EdgesOnShape*  eos = _data1->GetShapeEdges( fSide.EdgeID( iE ));
 
         bndPart._isShrink = ( eos->SWOLType() == TopAbs_FACE );
@@ -10613,10 +11082,14 @@ namespace VISCOUS_3D
                 bndPart._vertSWOLType[iV] = eov[iV]->SWOLType();
             }
           }
+          bndPart._edgeWOLLen[0] = fSide.EdgeLength( iE - 1 );
+          bndPart._edgeWOLLen[1] = fSide.EdgeLength( iE + 1 );
+
+          bndPart._tol = std::numeric_limits<double>::max(); // tolerance by segment size
+          for ( size_t i = 1; i < bndPart._nodes.size(); ++i )
+            bndPart._tol = Min( bndPart._tol,
+                                ( bndPart._nodes[i-1] - bndPart._nodes[i] ).SquareModulus() );
         }
-        std::vector<const SMDS_MeshNode*> nodes = fSide.GetOrderedNodes( iE );
-        bndPart._nodes.assign( nodes.begin(), nodes.end() );
-        bndPart._nbSegments = bndPart._nodes.size() - 1;
 
         if ( _boundary.empty() || ! _boundary.back().CanAppend( bndPart ))
           _boundary.push_back( bndPart );
@@ -10646,7 +11119,7 @@ namespace VISCOUS_3D
     {
       points.reserve( _boundarySize );
       size_t  nb = _boundary.rbegin()->_nodes.size();
-      int lastID = _boundary.rbegin()->Node( nb - 1 )->GetID();
+      smIdType lastID = _boundary.rbegin()->Node( nb - 1 )->GetID();
       std::list< BndPart >::const_iterator part = _boundary.begin();
       for ( ; part != _boundary.end(); ++part )
       {
@@ -10680,6 +11153,13 @@ namespace VISCOUS_3D
           return & _periodicFaces[ i ];
       return 0;
     }
+    void ClearPeriodic( const TopoDS_Face& face )
+    {
+      for ( size_t i = 0; i < _periodicFaces.size(); ++i )
+        if ( _periodicFaces[ i ]._shriFace[0]->IsSame( face ) ||
+             _periodicFaces[ i ]._shriFace[1]->IsSame( face ))
+          _periodicFaces[ i ].Clear();
+    }
   };
 
   //================================================================================
@@ -10689,6 +11169,7 @@ namespace VISCOUS_3D
   bool PeriodicFaces::IncludeShrunk( const TopoDS_Face&         face,
                                      const TopTools_MapOfShape& shrunkFaces ) const
   {
+    if ( IsEmpty() ) return false;
     return (( _shriFace[0]->IsSame( face ) && _shriFace[1]->IsShrunk( shrunkFaces )) ||
             ( _shriFace[1]->IsSame( face ) && _shriFace[0]->IsShrunk( shrunkFaces )));
   }
@@ -10709,11 +11190,15 @@ namespace VISCOUS_3D
     if ( iSrc != 0 )
     {
       trsfInverse = _trsf;
-      trsfInverse.Invert();
+      if ( !trsfInverse.Invert())
+        return false;
       trsf = &trsfInverse;
     }
     SMESHDS_Mesh* meshDS = dataSrc->GetHelper().GetMeshDS();
 
+    dumpFunction(SMESH_Comment("periodicMoveNodes_F")
+                               << _shriFace[iSrc]->_subMesh->GetId() << "_F"
+                               << _shriFace[iTgt]->_subMesh->GetId() );
     TNode2Edge::iterator n2e;
     TNodeNodeMap::iterator n2n = _nnMap.begin();
     for ( ; n2n != _nnMap.end(); ++n2n )
@@ -10743,14 +11228,17 @@ namespace VISCOUS_3D
           SMESH_NodeXYZ pSrc = leSrc->_nodes[ iN ];
           gp_XYZ pTgt = trsf->Transform( pSrc );
           meshDS->MoveNode( leTgt->_nodes[ iN ], pTgt.X(), pTgt.Y(), pTgt.Z() );
+
+          dumpMove( leTgt->_nodes[ iN ]);
         }
       }
     }
     bool done = ( n2n == _nnMap.end() );
-    // cout << "MMMMMMMOOOOOOOOOOVVVVVVVVVVVEEEEEEEE "
-    //      << _shriFace[iSrc]->_subMesh->GetId() << " -> "
-    //      << _shriFace[iTgt]->_subMesh->GetId() << " -- "
-    //      << ( done ? "DONE" : "FAIL") << endl;
+    debugMsg( "PeriodicFaces::MoveNodes "
+              << _shriFace[iSrc]->_subMesh->GetId() << " -> "
+              << _shriFace[iTgt]->_subMesh->GetId() << " -- "
+              << ( done ? "DONE" : "FAIL"));
+    dumpFunctionEnd();
 
     return done;
   }
@@ -11210,6 +11698,8 @@ bool _ViscousBuilder::shrink(_SolidData& theData)
             getMeshDS()->RemoveFreeNode( n, smDS, /*fromGroups=*/false );
         }
       }
+      _periodicity->ClearPeriodic( F );
+
       // restore position and UV of target nodes
       gp_Pnt p;
       for ( size_t iS = 0; iS < subEOS.size(); ++iS )
@@ -11283,6 +11773,8 @@ bool _ViscousBuilder::shrink(_SolidData& theData)
       {
         _EdgesOnShape& eos = * subEOS[ iS ];
         if ( eos.ShapeType() != TopAbs_EDGE ) continue;
+        if ( eos.size() == 0 )
+          continue;
 
         const TopoDS_Edge& E = TopoDS::Edge( eos._shape );
         data.SortOnEdge( E, eos._edges );
@@ -11305,6 +11797,19 @@ bool _ViscousBuilder::shrink(_SolidData& theData)
           uvPtVec[ i ].param = helper.GetNodeU( E, edges[i]->_nodes[0] );
           uvPtVec[ i ].SetUV( helper.GetNodeUV( F, edges[i]->_nodes.back() ));
         }
+        if ( uvPtVec[ 0 ].node == uvPtVec.back().node &&            // closed
+             helper.IsSeamShape( uvPtVec[ 0 ].node->GetShapeID() ))
+        {
+          uvPtVec[ 0 ].SetUV( helper.GetNodeUV( F,
+                                                edges[0]->_nodes.back(),
+                                                edges[1]->_nodes.back() ));
+          size_t i = edges.size() - 1;
+          uvPtVec[ i ].SetUV( helper.GetNodeUV( F,
+                                                edges[i  ]->_nodes.back(),
+                                                edges[i-1]->_nodes.back() ));
+        }
+        // if ( edges.empty() )
+        //   continue;
         BRep_Tool::Range( E, uvPtVec[0].param, uvPtVec.back().param );
         StdMeshers_FaceSide fSide( uvPtVec, F, E, _mesh );
         StdMeshers_ViscousLayers2D::SetProxyMeshOfEdge( fSide );
@@ -11319,8 +11824,12 @@ bool _ViscousBuilder::shrink(_SolidData& theData)
       smDS->Clear();
 
       // compute the mesh on the FACE
+      TopTools_IndexedMapOfShape allowed(1);
+      allowed.Add( sm->GetSubShape() );
+      sm->SetAllowedSubShapes( &allowed );
       sm->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE );
       sm->ComputeStateEngine( SMESH_subMesh::COMPUTE_SUBMESH );
+      sm->SetAllowedSubShapes( nullptr );
 
       // re-fill proxy sub-meshes of the FACE
       for ( size_t i = 0 ; i < _sdVec.size(); ++i )
@@ -11398,7 +11907,7 @@ bool _ViscousBuilder::shrink(_SolidData& theData)
 bool _ViscousBuilder::prepareEdgeToShrink( _LayerEdge&            edge,
                                            _EdgesOnShape&         eos,
                                            SMESH_MesherHelper&    helper,
-                                           const SMESHDS_SubMesh* faceSubMesh)
+                                           const SMESHDS_SubMesh* /*faceSubMesh*/)
 {
   const SMDS_MeshNode* srcNode = edge._nodes[0];
   const SMDS_MeshNode* tgtNode = edge._nodes.back();
@@ -11452,6 +11961,14 @@ bool _ViscousBuilder::prepareEdgeToShrink( _LayerEdge&            edge,
     if ( !n2 )
       return error(SMESH_Comment("Wrongly meshed EDGE ") << getMeshDS()->ShapeToIndex( E ));
 
+    if ( n2 == tgtNode       || // for 3D_mesh_GHS3D_01/B1
+         n2 == edge._nodes[1] ) // bos #20643
+    {
+      // shrunk by other SOLID
+      edge.Set( _LayerEdge::SHRUNK ); // ???
+      return true;
+    }
+
     double uSrc = helper.GetNodeU( E, srcNode, n2 );
     double uTgt = helper.GetNodeU( E, tgtNode, srcNode );
     double u2   = helper.GetNodeU( E, n2,      srcNode );
@@ -11707,8 +12224,8 @@ bool _LayerEdge::SetNewLength2d( Handle(Geom_Surface)& surface,
         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 );
+      gp_XY uvN1 = helper.GetNodeUV( F, _simplices[i]._nPrev, tgtNode );
+      gp_XY uvN2 = helper.GetNodeUV( F, _simplices[i]._nNext, tgtNode );
       gp_XY dirN = uvN2 - uvN1;
       double det = uvDir.Crossed( dirN );
       if ( Abs( det )  < std::numeric_limits<double>::min() ) continue;
@@ -11740,6 +12257,8 @@ bool _LayerEdge::SetNewLength2d( Handle(Geom_Surface)& surface,
     gp_Pnt p = surface->Value( newUV.X(), newUV.Y() );
     tgtNode->setXYZ( p.X(), p.Y(), p.Z() );
     dumpMove( tgtNode );
+#else
+    if ( surface.IsNull() ) {}
 #endif
   }
   else // _sWOL is TopAbs_EDGE
@@ -11911,17 +12430,16 @@ gp_XY _SmoothNode::computeAngularPos(vector<gp_XY>& uv,
   edgeSize.back() = edgeSize.front();
 
   gp_XY  newPos(0,0);
-  //int    nbEdges = 0;
-  double sumSize = 0;
+  double sumWgt = 0;
   for ( size_t i = 1; i < edgeDir.size(); ++i )
   {
-    if ( edgeDir[i-1].X() > 1. ) continue;
-    int i1 = i-1;
+    const int i1 = i-1;
+    if ( edgeDir[i1].X() > 1. ) continue;
     while ( edgeDir[i].X() > 1. && ++i < edgeDir.size() );
     if ( i == edgeDir.size() ) break;
     gp_XY p = uv[i];
     gp_XY norm1( -edgeDir[i1].Y(), edgeDir[i1].X() );
-    gp_XY norm2( -edgeDir[i].Y(),  edgeDir[i].X() );
+    gp_XY norm2( -edgeDir[i ].Y(), edgeDir[i ].X() );
     gp_XY bisec = norm1 + norm2;
     double bisecSize = bisec.Modulus();
     if ( bisecSize < numeric_limits<double>::min() )
@@ -11931,16 +12449,16 @@ gp_XY _SmoothNode::computeAngularPos(vector<gp_XY>& uv,
     }
     bisec /= bisecSize;
 
-    gp_XY  dirToN  = uvToFix - p;
-    double distToN = dirToN.Modulus();
+    gp_XY   dirToN = uvToFix - p;
+    double distToN = bisec * dirToN;
     if ( bisec * dirToN < 0 )
       distToN = -distToN;
 
-    newPos += ( p + bisec * distToN ) * ( edgeSize[i1] + edgeSize[i] );
-    //++nbEdges;
-    sumSize += edgeSize[i1] + edgeSize[i];
+    double wgt = edgeSize[i1] + edgeSize[i];
+    newPos += ( p + bisec * distToN ) * wgt;
+    sumWgt += wgt;
   }
-  newPos /= /*nbEdges * */sumSize;
+  newPos /= sumWgt;
   return newPos;
 }
 
@@ -11975,11 +12493,18 @@ void _Shrinker1D::AddEdge( const _LayerEdge*   e,
   double u = helper.GetNodeU( _geomEdge, e->_nodes[0], e->_nodes.back());
   _edges[ u < 0.5*(f+l) ? 0 : 1 ] = e;
 
-  // Update _nodes
+  // Check if the nodes are already shrunk by another SOLID
 
   const SMDS_MeshNode* tgtNode0 = TgtNode( 0 );
   const SMDS_MeshNode* tgtNode1 = TgtNode( 1 );
 
+  _done = (( tgtNode0 && tgtNode0->NbInverseElements( SMDSAbs_Edge ) == 2 ) ||
+           ( tgtNode1 && tgtNode1->NbInverseElements( SMDSAbs_Edge ) == 2 ));
+  if ( _done )
+    _nodes.resize( 1, nullptr );
+
+  // Update _nodes
+
   if ( _nodes.empty() )
   {
     SMESHDS_SubMesh * eSubMesh = helper.GetMeshDS()->MeshElements( _geomEdge );
@@ -11990,7 +12515,7 @@ void _Shrinker1D::AddEdge( const _LayerEdge*   e,
     GeomAdaptor_Curve aCurve(C, f,l);
     const double totLen = GCPnts_AbscissaPoint::Length(aCurve, f, l);
 
-    int nbExpectNodes = eSubMesh->NbNodes();
+    smIdType nbExpectNodes = eSubMesh->NbNodes();
     _initU  .reserve( nbExpectNodes );
     _normPar.reserve( nbExpectNodes );
     _nodes  .reserve( nbExpectNodes );
@@ -12048,6 +12573,8 @@ void _Shrinker1D::Compute(bool set3D, SMESH_MesherHelper& helper)
   double f,l;
   if ( set3D || _done )
   {
+    dumpFunction(SMESH_Comment("shrink1D_E") << helper.GetMeshDS()->ShapeToIndex( _geomEdge )<<
+                 "_F" << helper.GetSubShapeID() );
     Handle(Geom_Curve) C = BRep_Tool::Curve(_geomEdge, f,l);
     GeomAdaptor_Curve aCurve(C, f,l);
 
@@ -12069,7 +12596,9 @@ void _Shrinker1D::Compute(bool set3D, SMESH_MesherHelper& helper)
       pos->SetUParameter( u );
       gp_Pnt p = C->Value( u );
       const_cast< SMDS_MeshNode*>( _nodes[i] )->setXYZ( p.X(), p.Y(), p.Z() );
+      dumpMove( _nodes[i] );
     }
+    dumpFunctionEnd();
   }
   else
   {
@@ -12142,6 +12671,140 @@ void _Shrinker1D::SwapSrcTgtNodes( SMESHDS_Mesh* mesh )
   }
 }
 
+//================================================================================
+/*!
+ * \brief Setup quadPoints
+ */
+//================================================================================
+
+_Mapper2D::_Mapper2D( const TParam2ColumnMap & param2ColumnMap, const TNode2Edge& n2eMap )
+{
+  size_t i, iSize = _quadPoints.iSize = param2ColumnMap.size();
+  size_t j, jSize = _quadPoints.jSize = param2ColumnMap.begin()->second.size();
+  if ( _quadPoints.iSize < 3 ||
+       _quadPoints.jSize < 3 )
+    return;
+  _quadPoints.uv_grid.resize( iSize * jSize );
+
+  // set nodes
+  i = 0;
+  for ( auto & u_columnNodes : param2ColumnMap )
+  {
+    for ( j = 0; j < u_columnNodes.second.size(); ++j )
+      _quadPoints.UVPt( i, j ).node = u_columnNodes.second[ j ];
+    ++i;
+  }
+
+  // compute x parameter on borders
+  uvPnt( 0, 0       ).x = 0;
+  uvPnt( 0, jSize-1 ).x = 0;
+  gp_Pnt p0, pPrev0 = SMESH_NodeXYZ( uvPnt( 0, 0       ).node );
+  gp_Pnt p1, pPrev1 = SMESH_NodeXYZ( uvPnt( 0, jSize-1 ).node );
+  for ( i = 1; i < iSize; ++i )
+  {
+    p0 = SMESH_NodeXYZ( uvPnt( i, 0       ).node );
+    p1 = SMESH_NodeXYZ( uvPnt( i, jSize-1 ).node );
+    uvPnt( i, 0       ).x = uvPnt( i-1, 0       ).x + p0.Distance( pPrev0 );
+    uvPnt( i, jSize-1 ).x = uvPnt( i-1, jSize-1 ).x + p1.Distance( pPrev1 );
+    pPrev0 = p0;
+    pPrev1 = p1;
+  }
+  for ( i = 1; i < iSize-1; ++i )
+  {
+    uvPnt( i, 0       ).x /= uvPnt( iSize-1, 0       ).x;
+    uvPnt( i, jSize-1 ).x /= uvPnt( iSize-1, jSize-1 ).x;
+    uvPnt( i, 0       ).y = 0;
+    uvPnt( i, jSize-1 ).y = 1;
+  }
+
+  // compute y parameter on borders
+  uvPnt( 0,       0 ).y = 0;
+  uvPnt( iSize-1, 0 ).y = 0;
+  pPrev0 = SMESH_NodeXYZ( uvPnt( 0,       0 ).node );
+  pPrev1 = SMESH_NodeXYZ( uvPnt( iSize-1, 0 ).node );
+  for ( j = 1; j < jSize; ++j )
+  {
+    p0 = SMESH_NodeXYZ( uvPnt( 0,       j ).node );
+    p1 = SMESH_NodeXYZ( uvPnt( iSize-1, j ).node );
+    uvPnt( 0,       j ).y = uvPnt( 0,       j-1 ).y + p0.Distance( pPrev0 );
+    uvPnt( iSize-1, j ).y = uvPnt( iSize-1, j-1 ).y + p1.Distance( pPrev1 );
+    pPrev0 = p0;
+    pPrev1 = p1;
+  }
+  for ( j = 1; j < jSize-1; ++j )
+  {
+    uvPnt( 0,       j ).y /= uvPnt( 0,       jSize-1 ).y;
+    uvPnt( iSize-1, j ).y /= uvPnt( iSize-1, jSize-1 ).y;
+    uvPnt( 0,       j ).x = 0;
+    uvPnt( iSize-1, j ).x = 1;
+  }
+
+  // compute xy of internal nodes
+  for ( i = 1; i < iSize-1; ++i )
+  {
+    const double x0 = uvPnt( i, 0       ).x;
+    const double x1 = uvPnt( i, jSize-1 ).x;
+    for ( j = 1; j < jSize-1; ++j )
+    {
+      const double y0 = uvPnt( 0,       j ).y;
+      const double y1 = uvPnt( iSize-1, j ).y;
+      double x = (x0 + y0 * (x1 - x0)) / (1 - (y1 - y0) * (x1 - x0));
+      double y = y0 + x * (y1 - y0);
+      uvPnt( i, j ).x = x;
+      uvPnt( i, j ).y = y;
+    }
+  }
+
+  // replace base nodes with target ones
+  for ( i = 0; i < iSize; ++i )
+    for ( j = 0; j < jSize; ++j )
+    {
+      auto n2e = n2eMap.find( uvPnt( i, j ).node );
+      uvPnt( i, j ).node = n2e->second->_nodes.back();
+    }
+
+  return;
+}
+
+//================================================================================
+/*!
+ * \brief Compute positions of nodes of 2D structured mesh using TFI
+ */
+//================================================================================
+
+bool _Mapper2D::ComputeNodePositions()
+{
+  if ( _quadPoints.uv_grid.empty() )
+    return true;
+
+  size_t i, iSize = _quadPoints.iSize;
+  size_t j, jSize = _quadPoints.jSize;
+
+  SMESH_NodeXYZ a0 ( uvPnt( 0,       0       ).node );
+  SMESH_NodeXYZ a1 ( uvPnt( iSize-1, 0       ).node );
+  SMESH_NodeXYZ a2 ( uvPnt( iSize-1, jSize-1 ).node );
+  SMESH_NodeXYZ a3 ( uvPnt( 0,       jSize-1 ).node );
+
+  for ( i = 1; i < iSize-1; ++i )
+  {
+    SMESH_NodeXYZ p0 ( uvPnt( i, 0       ).node );
+    SMESH_NodeXYZ p2 ( uvPnt( i, jSize-1 ).node );
+    for ( j = 1; j < jSize-1; ++j )
+    {
+      SMESH_NodeXYZ p1 ( uvPnt( iSize-1, j ).node );
+      SMESH_NodeXYZ p3 ( uvPnt( 0,       j ).node );
+      double x = uvPnt( i, j ).x;
+      double y = uvPnt( i, j ).y;
+
+      gp_XYZ p = SMESH_MesherHelper::calcTFI( x, y, a0,a1,a2,a3, p0,p1,p2,p3 );
+      const_cast< SMDS_MeshNode* >( uvPnt( i, j ).node )->setXYZ( p.X(), p.Y(), p.Z() );
+
+      dumpMove( uvPnt( i, j ).node );
+    }
+  }
+  return true;
+}
+
 //================================================================================
 /*!
  * \brief Creates 2D and 1D elements on boundaries of new prisms