Salome HOME
IPAL54633: Bad viscous layers
[modules/smesh.git] / src / StdMeshers / StdMeshers_ViscousLayers.cxx
index 8212b0a8b4d4f407e847ce6376be88d99f9a0162..ea9d2f55ca7496aac1b6156709ef174e8d976eb7 100644 (file)
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2016  CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2019  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 "StdMeshers_ViscousLayers.hxx"
 
+#include "ObjectPool.hxx"
 #include "SMDS_EdgePosition.hxx"
 #include "SMDS_FaceOfNodes.hxx"
 #include "SMDS_FacePosition.hxx"
 #include "SMDS_MeshNode.hxx"
+#include "SMDS_PolygonalFaceOfNodes.hxx"
 #include "SMDS_SetIterator.hxx"
 #include "SMESHDS_Group.hxx"
 #include "SMESHDS_Hypothesis.hxx"
 #include "SMESHDS_Mesh.hxx"
 #include "SMESH_Algo.hxx"
+#include "SMESH_Block.hxx"
 #include "SMESH_ComputeError.hxx"
 #include "SMESH_ControlsDef.hxx"
 #include "SMESH_Gen.hxx"
 #include "SMESH_HypoFilter.hxx"
 #include "SMESH_Mesh.hxx"
 #include "SMESH_MeshAlgos.hxx"
+#include "SMESH_MeshEditor.hxx"
 #include "SMESH_MesherHelper.hxx"
 #include "SMESH_ProxyMesh.hxx"
 #include "SMESH_subMesh.hxx"
 #include "SMESH_subMeshEventListener.hxx"
 #include "StdMeshers_FaceSide.hxx"
+#include "StdMeshers_ProjectionUtils.hxx"
 #include "StdMeshers_ViscousLayers2D.hxx"
 
 #include <Adaptor3d_HSurface.hxx>
 #include <list>
 #include <queue>
 #include <string>
+#include <unordered_map>
 
 #ifdef _DEBUG_
-#define __myDEBUG
+//#define __myDEBUG
 //#define __NOT_INVALIDATE_BAD_SMOOTH
 //#define __NODES_AT_POS
 #endif
@@ -114,7 +120,7 @@ namespace VISCOUS_3D
   enum UIndex { U_TGT = 1, U_SRC, LEN_TGT };
 
   const double theMinSmoothCosin = 0.1;
-  const double theSmoothThickToElemSizeRatio = 0.3;
+  const double theSmoothThickToElemSizeRatio = 0.6;
   const double theMinSmoothTriaAngle = 30;
   const double theMinSmoothQuadAngle = 45;
 
@@ -369,22 +375,7 @@ namespace VISCOUS_3D
     double   _h2lenRatio; // avgNormProj / (2*avgDist)
     gp_Pnt2d _uv; // UV used in putOnOffsetSurface()
   public:
-    static _Curvature* New( double avgNormProj, double avgDist )
-    {
-      _Curvature* c = 0;
-      if ( fabs( avgNormProj / avgDist ) > 1./200 )
-      {
-        c = new _Curvature;
-        c->_r = avgDist * avgDist / avgNormProj;
-        c->_k = avgDist * avgDist / c->_r / c->_r;
-        //c->_k = avgNormProj / c->_r;
-        c->_k *= ( c->_r < 0 ? 1/1.1 : 1.1 ); // not to be too restrictive
-        c->_h2lenRatio = avgNormProj / ( avgDist + avgDist );
-
-        c->_uv.SetCoord( 0., 0. );
-      }
-      return c;
-    }
+    static _Curvature* New( double avgNormProj, double avgDist );
     double lenDelta(double len) const { return _k * ( _r + len ); }
     double lenDeltaByDist(double dist) const { return dist * _h2lenRatio; }
   };
@@ -500,9 +491,12 @@ namespace VISCOUS_3D
     gp_Ax1 LastSegment(double& segLen, _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 ); }
+    int    BaseShapeDim() const { return _nodes[0]->GetPosition()->GetDim(); }
     gp_XYZ Copy( _LayerEdge& other, _EdgesOnShape& eos, SMESH_MesherHelper& helper );
     void   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
     bool   IsNeiborOnEdge( const _LayerEdge* edge ) const;
     void   SetSmooLen( double len ) { // set _len at which smoothing is needed
@@ -577,6 +571,7 @@ namespace VISCOUS_3D
     gp_XYZ*              _plnNorm;
 
     _2NearEdges() { _edges[0]=_edges[1]=0; _plnNorm = 0; }
+    ~_2NearEdges(){ delete _plnNorm; }
     const SMDS_MeshNode* tgtNode(bool is2nd) {
       return _edges[is2nd] ? _edges[is2nd]->_nodes.back() : 0;
     }
@@ -617,12 +612,16 @@ 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; }
 
     bool   UseSurfaceNormal()  const
     { return _method == StdMeshers_ViscousLayers::SURF_OFFSET_SMOOTH; }
@@ -631,9 +630,19 @@ namespace VISCOUS_3D
     bool   IsOffsetMethod()    const
     { return _method == StdMeshers_ViscousLayers::FACE_OFFSET; }
 
+    bool operator==( const AverageHyp& other ) const
+    {
+      return ( _nbLayers == other._nbLayers &&
+               _method   == other._method   &&
+               Equals( GetTotalThickness(), other.GetTotalThickness() ) &&
+               Equals( GetStretchFactor(), other.GetStretchFactor() ));
+    }
+    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;
   };
 
   //--------------------------------------------------------------------------------
@@ -657,7 +666,8 @@ namespace VISCOUS_3D
     vector< _EdgesOnShape* > _eosConcaVer; // edges at concave VERTEXes of a FACE
     vector< _EdgesOnShape* > _eosC1; // to smooth together several C1 continues shapes
 
-    vector< gp_XYZ >         _faceNormals; // if _shape is FACE
+    typedef std::unordered_map< const SMDS_MeshElement*, gp_XYZ > TFace2NormMap;
+    TFace2NormMap            _faceNormals; // if _shape is FACE
     vector< _EdgesOnShape* > _faceEOS; // to get _faceNormals of adjacent FACEs
 
     Handle(ShapeAnalysis_Surface) _offsetSurf;
@@ -677,6 +687,7 @@ namespace VISCOUS_3D
     _SolidData&      GetData() const { return *_data; }
 
     _EdgesOnShape(): _shapeID(-1), _subMesh(0), _toSmooth(false), _edgeSmoother(0) {}
+    ~_EdgesOnShape();
   };
 
   //--------------------------------------------------------------------------------
@@ -736,6 +747,7 @@ namespace VISCOUS_3D
     TopTools_MapOfShape             _before; // SOLIDs to be computed before _solid
     TGeomID                         _index; // SOLID id
     _MeshOfSolid*                   _proxyMesh;
+    bool                            _done;
     list< THyp >                    _hyps;
     list< TopoDS_Shape >            _hypShapes;
     map< TGeomID, THyp >            _face2hyp; // filled if _hyps.size() > 1
@@ -752,7 +764,7 @@ namespace VISCOUS_3D
     // _LayerEdge's with underlying shapes
     vector< _EdgesOnShape >         _edgesOnShape;
 
-    // key:   an id of shape (EDGE or VERTEX) shared by a FACE with
+    // key:   an ID of shape (EDGE or VERTEX) shared by a FACE with
     //        layers and a FACE w/o layers
     // value: the shape (FACE or EDGE) to shrink mesh on.
     //       _LayerEdge's basing on nodes on key shape are inflated along the value shape
@@ -761,7 +773,7 @@ namespace VISCOUS_3D
     // Convex FACEs whose radius of curvature is less than the thickness of layers
     map< TGeomID, _ConvexFace >      _convexFaces;
 
-    // shapes (EDGEs and VERTEXes) srink from which is forbidden due to collisions with
+    // shapes (EDGEs and VERTEXes) shrink from which is forbidden due to collisions with
     // the adjacent SOLID
     set< TGeomID >                   _noShrinkShapes;
 
@@ -779,8 +791,8 @@ namespace VISCOUS_3D
 
     _SolidData(const TopoDS_Shape& s=TopoDS_Shape(),
                _MeshOfSolid*       m=0)
-      :_solid(s), _proxyMesh(m), _helper(0) {}
-    ~_SolidData();
+      :_solid(s), _proxyMesh(m), _done(false),_helper(0) {}
+    ~_SolidData() { delete _helper; _helper = 0; }
 
     void SortOnEdge( const TopoDS_Edge& E, vector< _LayerEdge* >& edges);
     void Sort2NeiborsOnEdge( vector< _LayerEdge* >& edges );
@@ -844,6 +856,8 @@ namespace VISCOUS_3D
 
     void Append( const gp_Pnt& center, _LayerEdge* ledge )
     {
+      if ( ledge->Is( _LayerEdge::MULTI_NORMAL ))
+        return;
       if ( _curvaCenters.size() > 0 )
         _segLength2.push_back( center.SquareDistance( _curvaCenters.back() ));
       _curvaCenters.push_back( center );
@@ -879,6 +893,7 @@ namespace VISCOUS_3D
                             const double   refSign );
   };
   struct PyDump;
+  struct Periodicity;
   //--------------------------------------------------------------------------------
   /*!
    * \brief Builder of viscous layers
@@ -905,10 +920,12 @@ namespace VISCOUS_3D
     bool findSolidsWithLayers();
     bool setBefore( _SolidData& solidBefore, _SolidData& solidAfter );
     bool findFacesWithLayers(const bool onlyWith=false);
+    void findPeriodicFaces();
     void getIgnoreFaces(const TopoDS_Shape&             solid,
                         const StdMeshers_ViscousLayers* hyp,
                         const TopoDS_Shape&             hypShape,
                         set<TGeomID>&                   ignoreFaces);
+    void makeEdgesOnShape();
     bool makeLayer(_SolidData& data);
     void setShapeData( _EdgesOnShape& eos, SMESH_subMesh* sm, _SolidData& data );
     bool setEdgeData( _LayerEdge& edge, _EdgesOnShape& eos,
@@ -959,7 +976,7 @@ namespace VISCOUS_3D
     void limitMaxLenByCurvature( _SolidData& data, SMESH_MesherHelper& helper );
     void limitMaxLenByCurvature( _LayerEdge* e1, _LayerEdge* e2,
                                  _EdgesOnShape& eos1, _EdgesOnShape& eos2,
-                                 SMESH_MesherHelper& helper );
+                                 const bool isSmoothable );
     bool updateNormals( _SolidData& data, SMESH_MesherHelper& helper, int stepNb, double stepSize );
     bool updateNormalsOfConvexFaces( _SolidData&         data,
                                      SMESH_MesherHelper& helper,
@@ -991,15 +1008,16 @@ namespace VISCOUS_3D
     // debug
     void makeGroupOfLE();
 
-    SMESH_Mesh*                _mesh;
-    SMESH_ComputeErrorPtr      _error;
+    SMESH_Mesh*                  _mesh;
+    SMESH_ComputeErrorPtr        _error;
 
-    vector<                    _SolidData >  _sdVec;
-    TopTools_IndexedMapOfShape _solids; // to find _SolidData by a solid
-    TopTools_MapOfShape        _shrinkedFaces;
+    vector<                      _SolidData >  _sdVec;
+    TopTools_IndexedMapOfShape   _solids; // to find _SolidData by a solid
+    TopTools_MapOfShape          _shrunkFaces;
+    std::unique_ptr<Periodicity> _periodicity;
 
-    int                        _tmpFaceID;
-    PyDump*                    _pyDump;
+    int                          _tmpFaceID;
+    PyDump*                      _pyDump;
   };
   //--------------------------------------------------------------------------------
   /*!
@@ -1084,31 +1102,32 @@ namespace VISCOUS_3D
       return _eos._edges[ is2nd ? _eos._edges.size()-1 : 0 ]->_2neibors->_edges[ is2nd ];
     }
     bool isAnalytic() const { return !_anaCurve.IsNull(); }
+
+    void offPointsToPython() const; // debug
   };
   //--------------------------------------------------------------------------------
   /*!
    * \brief Class of temporary mesh face.
    * We can't use SMDS_FaceOfNodes since it's impossible to set it's ID which is
-   * needed because SMESH_ElementSearcher internaly uses set of elements sorted by ID
+   * needed because SMESH_ElementSearcher internally uses set of elements sorted by ID
    */
-  struct _TmpMeshFace : public SMDS_MeshElement
+  struct _TmpMeshFace : public SMDS_PolygonalFaceOfNodes
   {
-    vector<const SMDS_MeshNode* > _nn;
+    const SMDS_MeshElement* _srcFace;
+
     _TmpMeshFace( const vector<const SMDS_MeshNode*>& nodes,
-                  int id, int faceID=-1, int idInFace=-1):
-      SMDS_MeshElement(id), _nn(nodes) { setShapeId(faceID); setIdInShape(idInFace); }
-    virtual const SMDS_MeshNode* GetNode(const int ind) const { return _nn[ind]; }
-    virtual SMDSAbs_ElementType  GetType() const              { return SMDSAbs_Face; }
-    virtual vtkIdType GetVtkType() const                      { return -1; }
-    virtual SMDSAbs_EntityType   GetEntityType() const        { return SMDSEntity_Last; }
-    virtual SMDSAbs_GeometryType GetGeomType() const
-    { return _nn.size() == 3 ? SMDSGeom_TRIANGLE : SMDSGeom_QUADRANGLE; }
-    virtual SMDS_ElemIteratorPtr elementsIterator(SMDSAbs_ElementType) const
-    { return SMDS_ElemIteratorPtr( new SMDS_NodeVectorElemIterator( _nn.begin(), _nn.end()));}
+                  int                                 ID,
+                  int                                 faceID=-1,
+                  const SMDS_MeshElement*             srcFace=0 ):
+      SMDS_PolygonalFaceOfNodes(nodes), _srcFace( srcFace ) { setID( ID ); setShapeID( faceID ); }
+    virtual SMDSAbs_EntityType  GetEntityType() const
+    { return _srcFace ? _srcFace->GetEntityType() : SMDSEntity_Quadrangle; }
+    virtual SMDSAbs_GeometryType GetGeomType()  const
+    { return _srcFace ? _srcFace->GetGeomType() : SMDSGeom_QUADRANGLE; }
   };
   //--------------------------------------------------------------------------------
   /*!
-   * \brief Class of temporary mesh face storing _LayerEdge it's based on
+   * \brief Class of temporary mesh quadrangle face storing _LayerEdge it's based on
    */
   struct _TmpMeshFaceOnEdge : public _TmpMeshFace
   {
@@ -1116,17 +1135,21 @@ namespace VISCOUS_3D
     _TmpMeshFaceOnEdge( _LayerEdge* le1, _LayerEdge* le2, int ID ):
       _TmpMeshFace( vector<const SMDS_MeshNode*>(4), ID ), _le1(le1), _le2(le2)
     {
-      _nn[0]=_le1->_nodes[0];
-      _nn[1]=_le1->_nodes.back();
-      _nn[2]=_le2->_nodes.back();
-      _nn[3]=_le2->_nodes[0];
+      myNodes[0]=_le1->_nodes[0];
+      myNodes[1]=_le1->_nodes.back();
+      myNodes[2]=_le2->_nodes.back();
+      myNodes[3]=_le2->_nodes[0];
+    }
+    const SMDS_MeshNode* n( size_t i ) const
+    {
+      return myNodes[ i ];
     }
     gp_XYZ GetDir() const // return average direction of _LayerEdge's, normal to EDGE
     {
-      SMESH_TNodeXYZ p0s( _nn[0] );
-      SMESH_TNodeXYZ p0t( _nn[1] );
-      SMESH_TNodeXYZ p1t( _nn[2] );
-      SMESH_TNodeXYZ p1s( _nn[3] );
+      SMESH_TNodeXYZ p0s( myNodes[0] );
+      SMESH_TNodeXYZ p0t( myNodes[1] );
+      SMESH_TNodeXYZ p1t( myNodes[2] );
+      SMESH_TNodeXYZ p1s( myNodes[3] );
       gp_XYZ  v0 = p0t - p0s;
       gp_XYZ  v1 = p1t - p1s;
       gp_XYZ v01 = p1s - p0s;
@@ -1137,10 +1160,10 @@ namespace VISCOUS_3D
     }
     gp_XYZ GetDir(_LayerEdge* le1, _LayerEdge* le2) // return average direction of _LayerEdge's
     {
-      _nn[0]=le1->_nodes[0];
-      _nn[1]=le1->_nodes.back();
-      _nn[2]=le2->_nodes.back();
-      _nn[3]=le2->_nodes[0];
+      myNodes[0]=le1->_nodes[0];
+      myNodes[1]=le1->_nodes.back();
+      myNodes[2]=le2->_nodes.back();
+      myNodes[3]=le2->_nodes[0];
       return GetDir();
     }
   };
@@ -1198,6 +1221,27 @@ namespace VISCOUS_3D
             ( dot * dot ) / l1 / l2 >= ( cos * cos ));
   }
 
+  class _Factory
+  {
+    ObjectPool< _LayerEdge >  _edgePool;
+    ObjectPool< _Curvature >  _curvaturePool;
+    ObjectPool< _2NearEdges > _nearEdgesPool;
+
+    static _Factory* & me()
+    {
+      static _Factory* theFactory = 0;
+      return theFactory;
+    }
+  public:
+
+    _Factory()  { me() = this; }
+    ~_Factory() { me() = 0; }
+
+    static _LayerEdge*  NewLayerEdge() { return me()->_edgePool.getNew(); }
+    static _Curvature * NewCurvature() { return me()->_curvaturePool.getNew(); }
+    static _2NearEdges* NewNearEdges() { return me()->_nearEdgesPool.getNew(); }
+  };
+
 } // namespace VISCOUS_3D
 
 
@@ -1205,10 +1249,11 @@ namespace VISCOUS_3D
 //================================================================================
 // StdMeshers_ViscousLayers hypothesis
 //
-StdMeshers_ViscousLayers::StdMeshers_ViscousLayers(int hypId, int studyId, SMESH_Gen* gen)
-  :SMESH_Hypothesis(hypId, studyId, gen),
+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
@@ -1240,6 +1285,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,
@@ -1294,6 +1348,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)
@@ -1306,6 +1363,13 @@ 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
@@ -1339,6 +1403,36 @@ 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
 //================================================================================
 
@@ -1577,7 +1671,7 @@ namespace VISCOUS_3D
 
   //================================================================================
   /*!
-   * \brief Computes mimimal distance of face in-FACE nodes from an EDGE
+   * \brief Computes minimal distance of face in-FACE nodes from an EDGE
    *  \param [in] face - the mesh face to treat
    *  \param [in] nodeOnEdge - a node on the EDGE
    *  \param [out] faceSize - the computed distance
@@ -1603,8 +1697,8 @@ namespace VISCOUS_3D
     // look for two neighbor not in-FACE nodes of face
     for ( int i = 0; i < 2; ++i )
     {
-      if ( nNext[i]->GetPosition()->GetDim() != 2 &&
-           nNext[i]->GetID() < nodeOnEdge->GetID() )
+      if (( nNext[i]->GetPosition()->GetDim() != 2 ) &&
+          ( nodeOnEdge->GetPosition()->GetDim() == 0 || nNext[i]->GetID() < nodeOnEdge->GetID() ))
       {
         // look for an in-FACE node
         for ( int iN = 0; iN < nbN; ++iN )
@@ -1677,12 +1771,12 @@ 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
-          << "smesh  = smeshBuilder.New(salome.myStudy)" << endl
-          << "meshSO = smesh.GetCurrentStudy().FindObjectID('0:1:2:" << tag <<"')" << endl
+          << "smesh  = smeshBuilder.New()" << endl
+          << "meshSO = salome.myStudy.FindObjectID('0:1:2:" << tag <<"')" << endl
           << "mesh   = smesh.Mesh( meshSO.GetObject() )"<<endl;
       theNbPyFunc = 0;
     }
@@ -1859,6 +1953,8 @@ SMESH_ComputeErrorPtr _ViscousBuilder::Compute(SMESH_Mesh&         theMesh,
 {
   _mesh = & theMesh;
 
+  _Factory factory;
+
   // check if proxy mesh already computed
   TopExp_Explorer exp( theShape, TopAbs_SOLID );
   if ( !exp.More() )
@@ -1870,20 +1966,30 @@ SMESH_ComputeErrorPtr _ViscousBuilder::Compute(SMESH_Mesh&         theMesh,
   PyDump debugDump( theMesh );
   _pyDump = &debugDump;
 
-  // TODO: ignore already computed SOLIDs 
+  // TODO: ignore already computed SOLIDs
   if ( !findSolidsWithLayers())
     return _error;
 
   if ( !findFacesWithLayers() )
     return _error;
 
+  // for ( size_t i = 0; i < _sdVec.size(); ++i )
+  // {
+  //   if ( ! makeLayer( _sdVec[ i ]))   // create _LayerEdge's
+  //     return _error;
+  // }
+
+  makeEdgesOnShape();
+
+  findPeriodicFaces();
+
   for ( size_t i = 0; i < _sdVec.size(); ++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() )
+           !_sdVec[iSD]._done )
         break;
 
     if ( ! makeLayer(_sdVec[iSD]) )   // create _LayerEdge's
@@ -1906,6 +2012,8 @@ SMESH_ComputeErrorPtr _ViscousBuilder::Compute(SMESH_Mesh&         theMesh,
 
     addBoundaryElements(_sdVec[iSD]); // create quadrangles on prism bare sides
 
+    _sdVec[iSD]._done = true;
+
     const TopoDS_Shape& solid = _sdVec[iSD]._solid;
     for ( iSD = 0; iSD < _sdVec.size(); ++iSD )
       _sdVec[iSD]._before.Remove( solid );
@@ -2132,7 +2240,7 @@ bool _ViscousBuilder::findFacesWithLayers(const bool onlyWith)
     for ( ; exp.More(); exp.Next() )
     {
       const TopoDS_Face& face = TopoDS::Face( exp.Current() );
-      const TGeomID faceID = getMeshDS()->ShapeToIndex( face );
+      const TGeomID    faceID = getMeshDS()->ShapeToIndex( face );
       if ( //!sdVec[i]._ignoreFaceIds.count( faceID ) &&
           helper.NbAncestors( face, *_mesh, TopAbs_SOLID ) > 1 &&
           helper.IsReversedSubMesh( face ))
@@ -2168,7 +2276,7 @@ bool _ViscousBuilder::findFacesWithLayers(const bool onlyWith)
         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 ];
+      TopoDS_Shape fWOL = FF[ ignore[0] ? 0 : 1 ]; // FACE w/o layers
 
       // add EDGE to maps
       if ( !fWOL.IsNull())
@@ -2248,7 +2356,7 @@ bool _ViscousBuilder::findFacesWithLayers(const bool onlyWith)
     }
   }
 
-  // Add to _noShrinkShapes sub-shapes of FACE's that can't be shrinked since
+  // Add to _noShrinkShapes sub-shapes of FACE's that can't be shrunk 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 )
@@ -2338,7 +2446,7 @@ bool _ViscousBuilder::findFacesWithLayers(const bool onlyWith)
       for ( TopoDS_Iterator vIt( edge ); vIt.More(); vIt.Next() )
       {
         TGeomID vID = getMeshDS()->ShapeToIndex( vIt.Value() );
-        bool noShrinkV = false;
+        bool noShrinkV = false, noShrinkIfAdjMeshed = false;
 
         if ( iSolid < _sdVec.size() )
         {
@@ -2348,7 +2456,8 @@ bool _ViscousBuilder::findFacesWithLayers(const bool onlyWith)
             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 );
+              noShrinkV = (( isStructured ) ||
+                           ( noShrinkIfAdjMeshed = i2S->second.ShapeType() == TopAbs_EDGE ));
             else
               noShrinkV = ( ! i2S->second.IsSame( i2SAdj->second ));
           }
@@ -2360,8 +2469,31 @@ bool _ViscousBuilder::findFacesWithLayers(const bool onlyWith)
         else
         {
           // the adjacent SOLID has NO layers at all
-          noShrinkV = ( isStructured ||
-                        _sdVec[i]._shrinkShape2Shape[ vID ].ShapeType() == TopAbs_EDGE );
+          if ( isStructured )
+          {
+            noShrinkV = true;
+          }
+          else
+          {
+            noShrinkV = noShrinkIfAdjMeshed =
+              ( _sdVec[i]._shrinkShape2Shape[ vID ].ShapeType() == TopAbs_EDGE );
+          }
+        }
+
+        if ( noShrinkV && noShrinkIfAdjMeshed )
+        {
+          // noShrinkV if FACEs in the adjacent SOLID are meshed
+          PShapeIteratorPtr fIt = helper.GetAncestors( _sdVec[i]._shrinkShape2Shape[ vID ],
+                                                       *_mesh, TopAbs_FACE, &solid );
+          while ( fIt->more() )
+          {
+            const TopoDS_Shape* f = fIt->next();
+            if ( !f->IsSame( fWOL ))
+            {
+              noShrinkV = ! _mesh->GetSubMesh( *f )->IsEmpty();
+              break;
+            }
+          }
         }
         if ( noShrinkV )
           _sdVec[i]._noShrinkShapes.insert( vID );
@@ -2449,17 +2581,6 @@ void _ViscousBuilder::getIgnoreFaces(const TopoDS_Shape&             solid,
 
 bool _ViscousBuilder::makeLayer(_SolidData& data)
 {
-  // get all sub-shapes to make layers on
-  set<TGeomID> subIds, faceIds;
-  subIds = data._noShrinkShapes;
-  TopExp_Explorer exp( data._solid, TopAbs_FACE );
-  for ( ; exp.More(); exp.Next() )
-  {
-    SMESH_subMesh* fSubM = _mesh->GetSubMesh( exp.Current() );
-    if ( ! data._ignoreFaceIds.count( fSubM->GetId() ))
-      faceIds.insert( fSubM->GetId() );
-  }
-
   // make a map to find new nodes on sub-shapes shared with other SOLID
   map< TGeomID, TNode2Edge* >::iterator s2ne;
   map< TGeomID, TopoDS_Shape >::iterator s2s = data._shrinkShape2Shape.begin();
@@ -2483,6 +2604,8 @@ bool _ViscousBuilder::makeLayer(_SolidData& data)
 
   dumpFunction(SMESH_Comment("makeLayers_")<<data._index);
 
+  vector< _EdgesOnShape >& edgesByGeom = data._edgesOnShape;
+
   data._stepSize = Precision::Infinite();
   data._stepSizeNodes[0] = 0;
 
@@ -2493,34 +2616,19 @@ bool _ViscousBuilder::makeLayer(_SolidData& data)
   vector< const SMDS_MeshNode*> newNodes; // of a mesh face
   TNode2Edge::iterator n2e2;
 
-  // collect _LayerEdge's of shapes they are based on
-  vector< _EdgesOnShape >& edgesByGeom = data._edgesOnShape;
-  const int nbShapes = getMeshDS()->MaxShapeIndex();
-  edgesByGeom.resize( nbShapes+1 );
-
-  // set data of _EdgesOnShape's
-  if ( SMESH_subMesh* sm = _mesh->GetSubMesh( data._solid ))
-  {
-    SMESH_subMeshIteratorPtr smIt = sm->getDependsOnIterator(/*includeSelf=*/false);
-    while ( smIt->more() )
-    {
-      sm = smIt->next();
-      if ( sm->GetSubShape().ShapeType() == TopAbs_FACE &&
-           !faceIds.count( sm->GetId() ))
-        continue;
-      setShapeData( edgesByGeom[ sm->GetId() ], sm, data );
-    }
-  }
   // make _LayerEdge's
-  for ( set<TGeomID>::iterator id = faceIds.begin(); id != faceIds.end(); ++id )
+  for ( TopExp_Explorer exp( data._solid, TopAbs_FACE ); exp.More(); exp.Next() )
   {
-    const TopoDS_Face& F = TopoDS::Face( getMeshDS()->IndexToShape( *id ));
-    SMESH_subMesh* sm = _mesh->GetSubMesh( F );
+    const TopoDS_Face& F = TopoDS::Face( exp.Current() );
+    SMESH_subMesh*    sm = _mesh->GetSubMesh( F );
+    const TGeomID     id = sm->GetId();
+    if ( edgesByGeom[ id ]._shape.IsNull() )
+      continue; // no layers
     SMESH_ProxyMesh::SubMesh* proxySub =
       data._proxyMesh->getFaceSubM( F, /*create=*/true);
 
     SMESHDS_SubMesh* smDS = sm->GetSubMeshDS();
-    if ( !smDS ) return error(SMESH_Comment("Not meshed face ") << *id, data._index );
+    if ( !smDS ) return error(SMESH_Comment("Not meshed face ") << id, data._index );
 
     SMDS_ElemIteratorPtr eIt = smDS->GetElements();
     while ( eIt->more() )
@@ -2558,7 +2666,7 @@ bool _ViscousBuilder::makeLayer(_SolidData& data)
         if ( !(*n2e).second )
         {
           // add a _LayerEdge
-          _LayerEdge* edge = new _LayerEdge();
+          _LayerEdge* edge = _Factory::NewLayerEdge();
           edge->_nodes.push_back( n );
           n2e->second = edge;
           edgesByGeom[ shapeID ]._edges.push_back( edge );
@@ -2610,7 +2718,7 @@ bool _ViscousBuilder::makeLayer(_SolidData& data)
 
       // create a temporary face
       const SMDS_MeshElement* newFace =
-        new _TmpMeshFace( newNodes, --_tmpFaceID, face->getshapeId(), face->getIdInShape() );
+        new _TmpMeshFace( newNodes, --_tmpFaceID, face->GetShapeID(), face );
       proxySub->AddElement( newFace );
 
       // compute inflation step size by min size of element on a convex surface
@@ -2855,7 +2963,7 @@ void _ViscousBuilder::limitStepSizeByCurvature( _SolidData& data )
             double curvature = Max( surfProp.MaxCurvature() * oriFactor,
                                     surfProp.MinCurvature() * oriFactor );
             if ( curvature > minCurvature )
-              ledge->_maxLen = Min( ledge->_maxLen, 1. / curvature );
+              ledge->SetMaxLen( Min( ledge->_maxLen, 1. / curvature ));
           }
         }
       }
@@ -2876,7 +2984,12 @@ void _ViscousBuilder::limitStepSizeByCurvature( _SolidData& data )
         for ( size_t j = 0; j < ledge->_simplices.size(); ++j )
           if ( ledge->_simplices[j]._nNext->GetPosition()->GetDim() < 2 )
           {
-            convFace._simplexTestEdges.push_back( ledge );
+            // do not select _LayerEdge's neighboring sharp EDGEs
+            bool sharpNbr = false;
+            for ( size_t iN = 0; iN < ledge->_neibors.size()  && !sharpNbr; ++iN )
+              sharpNbr = ( ledge->_neibors[iN]->_cosin > theMinSmoothCosin );
+            if ( !sharpNbr )
+              convFace._simplexTestEdges.push_back( ledge );
             break;
           }
       }
@@ -2926,22 +3039,11 @@ 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
   // boundary inclined to the shape at a sharp angle
 
-  //list< TGeomID > shapesToSmooth;
   TopTools_MapOfShape edgesOfSmooFaces;
-
   SMESH_MesherHelper helper( *_mesh );
   bool ok = true;
 
@@ -2956,32 +3058,34 @@ bool _ViscousBuilder::findShapesToSmooth( _SolidData& data )
       continue;
 
     double tgtThick = eos._hyp.GetTotalThickness();
-    TopExp_Explorer eExp( edgesByGeom[iS]._shape, TopAbs_EDGE );
-    for ( ; eExp.More() && !eos._toSmooth; eExp.Next() )
+    SMESH_subMeshIteratorPtr subIt = eos._subMesh->getDependsOnIterator(/*includeSelf=*/false );
+    while ( subIt->more() && !eos._toSmooth )
     {
-      TGeomID iE = getMeshDS()->ShapeToIndex( eExp.Current() );
-      vector<_LayerEdge*>& eE = edgesByGeom[ iE ]._edges;
-      if ( eE.empty() ) continue;
+      TGeomID iSub = subIt->next()->GetId();
+      const vector<_LayerEdge*>& eSub = edgesByGeom[ iSub ]._edges;
+      if ( eSub.empty() ) continue;
 
       double faceSize;
-      for ( size_t i = 0; i < eE.size() && !eos._toSmooth; ++i )
-        if ( eE[i]->_cosin > theMinSmoothCosin )
+      for ( size_t i = 0; i < eSub.size() && !eos._toSmooth; ++i )
+        if ( eSub[i]->_cosin > theMinSmoothCosin )
         {
-          SMDS_ElemIteratorPtr fIt = eE[i]->_nodes[0]->GetInverseElementIterator(SMDSAbs_Face);
+          SMDS_ElemIteratorPtr fIt = eSub[i]->_nodes[0]->GetInverseElementIterator(SMDSAbs_Face);
           while ( fIt->more() && !eos._toSmooth )
           {
             const SMDS_MeshElement* face = fIt->next();
             if ( face->getshapeId() == eos._shapeID &&
-                 getDistFromEdge( face, eE[i]->_nodes[0], faceSize ))
+                 getDistFromEdge( face, eSub[i]->_nodes[0], faceSize ))
             {
-              eos._toSmooth = needSmoothing( eE[i]->_cosin, tgtThick, faceSize );
+              eos._toSmooth = needSmoothing( eSub[i]->_cosin,
+                                             tgtThick * eSub[i]->_lenFactor,
+                                             faceSize);
             }
           }
         }
     }
     if ( eos._toSmooth )
     {
-      for ( eExp.ReInit(); eExp.More(); eExp.Next() )
+      for ( TopExp_Explorer eExp( edgesByGeom[iS]._shape, TopAbs_EDGE ); eExp.More(); eExp.Next() )
         edgesOfSmooFaces.Add( eExp.Current() );
 
       data.PrepareEdgesToSmoothOnFace( &edgesByGeom[iS], /*substituteSrcNodes=*/false );
@@ -3025,8 +3129,8 @@ bool _ViscousBuilder::findShapesToSmooth( _SolidData& data )
           if ( endSeg->getshapeId() == (int) iS )
           {
             double segLen =
-              SMESH_TNodeXYZ( endSeg->GetNode(0) ).Distance( endSeg->GetNode(1 ));
-            eos._toSmooth = needSmoothing( cosinAbs, tgtThick, segLen );
+              SMESH_TNodeXYZ( endSeg->GetNode( 0 )).Distance( endSeg->GetNode( 1 ));
+            eos._toSmooth = needSmoothing( cosinAbs, tgtThick * eV[0]->_lenFactor, segLen );
           }
         }
         if ( eos._toSmooth )
@@ -3050,11 +3154,12 @@ bool _ViscousBuilder::findShapesToSmooth( _SolidData& data )
 
     if ( !eos._hyp.ToSmooth() )
       for ( size_t i = 0; i < eos._edges.size(); ++i )
-        eos._edges[i]->SetCosin( 0 );
+        //eos._edges[i]->SetCosin( 0 ); // keep _cosin to use in limitMaxLenByCurvature()
+        eos._edges[i]->_lenFactor = 1;
   }
 
 
-  // Fill _eosC1 to make that C1 FACEs and EGDEs between them to be smoothed as a whole
+  // Fill _eosC1 to make that C1 FACEs and EDGEs between them to be smoothed as a whole
 
   TopTools_MapOfShape c1VV;
 
@@ -3114,19 +3219,21 @@ bool _ViscousBuilder::findShapesToSmooth( _SolidData& data )
         {
           _EdgesOnShape* eof = data.GetShapeEdges( *face );
           if ( !eof ) continue; // other solid
-          if ( !eos.HasC1( eoe ))
-          {
-            eos._eosC1.push_back( eoe );
-            eoe->_toSmooth = false;
-            data.PrepareEdgesToSmoothOnFace( eoe, /*substituteSrcNodes=*/false );
-          }
-          if ( eos._shapeID != eof->_shapeID && !eos.HasC1( eof )) 
+          if ( eos._shapeID == eof->_shapeID ) continue;
+          if ( !eos.HasC1( eof ))
           {
+            // check the FACEs
             eos._eosC1.push_back( eof );
             eof->_toSmooth = false;
             data.PrepareEdgesToSmoothOnFace( eof, /*substituteSrcNodes=*/false );
             smQueue.push_back( eof->_subMesh );
           }
+          if ( !eos.HasC1( eoe ))
+          {
+            eos._eosC1.push_back( eoe );
+            eoe->_toSmooth = false;
+            data.PrepareEdgesToSmoothOnFace( eoe, /*substituteSrcNodes=*/false );
+          }
         }
       }
     }
@@ -3197,12 +3304,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
@@ -3212,6 +3329,39 @@ bool _ViscousBuilder::findShapesToSmooth( _SolidData& data )
   return ok;
 }
 
+//================================================================================
+/*!
+ * \brief Set up _SolidData::_edgesOnShape
+ */
+//================================================================================
+
+void _ViscousBuilder::makeEdgesOnShape()
+{
+  const int nbShapes = getMeshDS()->MaxShapeIndex();
+
+  for ( size_t i = 0; i < _sdVec.size(); ++i )
+  {
+    _SolidData& data = _sdVec[ i ];
+    vector< _EdgesOnShape >& edgesByGeom = data._edgesOnShape;
+    edgesByGeom.resize( nbShapes+1 );
+
+    // set data of _EdgesOnShape's
+    if ( SMESH_subMesh* sm = _mesh->GetSubMesh( data._solid ))
+    {
+      SMESH_subMeshIteratorPtr smIt = sm->getDependsOnIterator(/*includeSelf=*/false);
+      while ( smIt->more() )
+      {
+        sm = smIt->next();
+        if ( sm->GetSubShape().ShapeType() == TopAbs_FACE &&
+             data._ignoreFaceIds.count( sm->GetId() ))
+          continue;
+
+        setShapeData( edgesByGeom[ sm->GetId() ], sm, data );
+      }
+    }
+  }
+}
+
 //================================================================================
 /*!
  * \brief initialize data of _EdgesOnShape
@@ -3282,19 +3432,19 @@ void _ViscousBuilder::setShapeData( _EdgesOnShape& eos,
     if ( eos.ShapeType() == TopAbs_FACE ) // get normals to elements on a FACE
     {
       SMESHDS_SubMesh* smDS = sm->GetSubMeshDS();
-      eos._faceNormals.resize( smDS->NbElements() );
+      if ( !smDS ) return;
+      eos._faceNormals.reserve( smDS->NbElements() );
 
+      double oriFactor = helper.IsReversedSubMesh( TopoDS::Face( eos._shape )) ? 1.: -1.;
       SMDS_ElemIteratorPtr eIt = smDS->GetElements();
-      for ( int iF = 0; eIt->more(); ++iF )
+      for ( ; eIt->more(); )
       {
         const SMDS_MeshElement* face = eIt->next();
-        if ( !SMESH_MeshAlgos::FaceNormal( face, eos._faceNormals[iF], /*normalized=*/true ))
-          eos._faceNormals[iF].SetCoord( 0,0,0 );
+        gp_XYZ&                 norm = eos._faceNormals[face];
+        if ( !SMESH_MeshAlgos::FaceNormal( face, norm, /*normalized=*/true ))
+          norm.SetCoord( 0,0,0 );
+        norm *= oriFactor;
       }
-
-      if ( !helper.IsReversedSubMesh( TopoDS::Face( eos._shape )))
-        for ( size_t iF = 0; iF < eos._faceNormals.size(); ++iF )
-          eos._faceNormals[iF].Reverse();
     }
     else // find EOS of adjacent FACEs
     {
@@ -3320,7 +3470,7 @@ void _ViscousBuilder::setShapeData( _EdgesOnShape& eos,
 bool _EdgesOnShape::GetNormal( const SMDS_MeshElement* face, gp_Vec& norm )
 {
   bool ok = false;
-  const _EdgesOnShape* eos = 0;
+  _EdgesOnShape* eos = 0;
 
   if ( face->getshapeId() == _shapeID )
   {
@@ -3334,9 +3484,9 @@ bool _EdgesOnShape::GetNormal( const SMDS_MeshElement* face, gp_Vec& norm )
   }
 
   if (( eos ) &&
-      ( ok = ( face->getIdInShape() < (int) eos->_faceNormals.size() )))
+      ( ok = ( eos->_faceNormals.count( face ) )))
   {
-    norm = eos->_faceNormals[ face->getIdInShape() ];
+    norm = eos->_faceNormals[ face ];
   }
   else if ( !eos )
   {
@@ -3346,6 +3496,16 @@ bool _EdgesOnShape::GetNormal( const SMDS_MeshElement* face, gp_Vec& norm )
   return ok;
 }
 
+//================================================================================
+/*!
+ * \brief EdgesOnShape destructor
+ */
+//================================================================================
+
+_EdgesOnShape::~_EdgesOnShape()
+{
+  delete _edgeSmoother;
+}
 
 //================================================================================
 /*!
@@ -3360,12 +3520,13 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge&         edge,
 {
   const SMDS_MeshNode* node = edge._nodes[0]; // source node
 
-  edge._len       = 0;
-  edge._maxLen    = Precision::Infinite();
-  edge._minAngle  = 0;
-  edge._2neibors  = 0;
-  edge._curvature = 0;
-  edge._flags     = 0;
+  edge._len          = 0;
+  edge._maxLen       = Precision::Infinite();
+  edge._minAngle     = 0;
+  edge._2neibors     = 0;
+  edge._curvature    = 0;
+  edge._flags        = 0;
+  edge._smooFunction = 0;
 
   // --------------------------
   // Compute _normal and _cosin
@@ -3562,7 +3723,7 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge&         edge,
     getMeshDS()->RemoveFreeNode( edge._nodes.back(), 0, /*fromGroups=*/false );
     edge._nodes.resize( 1 );
     edge._normal.SetCoord( 0,0,0 );
-    edge._maxLen = 0;
+    edge.SetMaxLen( 0 );
   }
 
   // Set the rest data
@@ -3574,7 +3735,7 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge&         edge,
   {
     const SMDS_MeshNode* tgtNode = edge._nodes.back();
     if ( SMESHDS_SubMesh* sm = getMeshDS()->MeshElements( data._solid ))
-      sm->RemoveNode( tgtNode , /*isNodeDeleted=*/false );
+      sm->RemoveNode( tgtNode );
 
     // set initial position which is parameters on _sWOL in this case
     if ( eos.SWOLType() == TopAbs_EDGE )
@@ -3638,7 +3799,7 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge&         edge,
   if ( eos.ShapeType() == TopAbs_EDGE /*||
        ( onShrinkShape && posType == SMDS_TOP_VERTEX && fabs( edge._cosin ) < 1e-10 )*/)
   {
-    edge._2neibors = new _2NearEdges;
+    edge._2neibors = _Factory::NewNearEdges();
     // target nodes instead of source ones will be set later
   }
 
@@ -3699,7 +3860,7 @@ gp_XYZ _ViscousBuilder::getFaceNormal(const SMDS_MeshNode* node,
       isOK = false;
       return p.XYZ();
     }
-    Quantity_Parameter U,V;
+    Standard_Real U,V;
     projector.LowerDistanceParameters(U,V);
     uv.SetCoord( U,V );
   }
@@ -4054,7 +4215,7 @@ gp_XYZ _OffsetPlane::GetCommonPoint(bool&                 isFound,
 
 //================================================================================
 /*!
- * \brief Find 2 neigbor nodes of a node on EDGE
+ * \brief Find 2 neighbor nodes of a node on EDGE
  */
 //================================================================================
 
@@ -4099,7 +4260,35 @@ bool _ViscousBuilder::findNeiborsOnEdge(const _LayerEdge*     edge,
 
 //================================================================================
 /*!
- * \brief Set _curvature and _2neibors->_plnNorm by 2 neigbor nodes residing the same EDGE
+ * \brief Create _Curvature
+ */
+//================================================================================
+
+_Curvature* _Curvature::New( double avgNormProj, double avgDist )
+{
+  // double   _r; // radius
+  // double   _k; // factor to correct node smoothed position
+  // double   _h2lenRatio; // avgNormProj / (2*avgDist)
+  // gp_Pnt2d _uv; // UV used in putOnOffsetSurface()
+
+  _Curvature* c = 0;
+  if ( fabs( avgNormProj / avgDist ) > 1./200 )
+  {
+    c = _Factory::NewCurvature();
+    c->_r = avgDist * avgDist / avgNormProj;
+    c->_k = avgDist * avgDist / c->_r / c->_r;
+    //c->_k = avgNormProj / c->_r;
+    c->_k *= ( c->_r < 0 ? 1/1.1 : 1.1 ); // not to be too restrictive
+    c->_h2lenRatio = avgNormProj / ( avgDist + avgDist );
+
+    c->_uv.SetCoord( 0., 0. );
+  }
+  return c;
+}
+
+//================================================================================
+/*!
+ * \brief Set _curvature and _2neibors->_plnNorm by 2 neighbor nodes residing the same EDGE
  */
 //================================================================================
 
@@ -4110,6 +4299,8 @@ void _LayerEdge::SetDataByNeighbors( const SMDS_MeshNode* n1,
 {
   if ( eos.ShapeType() != TopAbs_EDGE )
     return;
+  if ( _curvature && Is( SMOOTHED_C1 ))
+    return;
 
   gp_XYZ  pos = SMESH_TNodeXYZ( _nodes[0] );
   gp_XYZ vec1 = pos - SMESH_TNodeXYZ( n1 );
@@ -4122,7 +4313,6 @@ void _LayerEdge::SetDataByNeighbors( const SMDS_MeshNode* n1,
   _2neibors->_wgt[1] = 1 - vec2.Modulus() / sumLen;
   double avgNormProj = 0.5 * ( _normal * vec1 + _normal * vec2 );
   double      avgLen = 0.5 * ( vec1.Modulus() + vec2.Modulus() );
-  if ( _curvature ) delete _curvature;
   _curvature = _Curvature::New( avgNormProj, avgLen );
   // if ( _curvature )
   //   debugMsg( _nodes[0]->GetID()
@@ -4166,8 +4356,11 @@ gp_XYZ _LayerEdge::Copy( _LayerEdge&         other,
   _lenFactor = other._lenFactor;
   _cosin     = other._cosin;
   _2neibors  = other._2neibors;
-  _curvature = 0; std::swap( _curvature, other._curvature );
-  _2neibors  = 0; std::swap( _2neibors,  other._2neibors );
+  _curvature = other._curvature;
+  _2neibors  = other._2neibors;
+  _maxLen    = Precision::Infinite();//other._maxLen;
+  _flags     = 0;
+  _smooFunction = 0;
 
   gp_XYZ lastPos( 0,0,0 );
   if ( eos.SWOLType() == TopAbs_EDGE )
@@ -4276,7 +4469,7 @@ void _Simplex::SortSimplices(vector<_Simplex>& simplices)
 
 //================================================================================
 /*!
- * \brief DEBUG. Create groups contating temorary data of _LayerEdge's
+ * \brief DEBUG. Create groups containing temporary data of _LayerEdge's
  */
 //================================================================================
 
@@ -4363,7 +4556,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() )
@@ -4383,16 +4576,88 @@ void _ViscousBuilder::computeGeomSize( _SolidData& data )
     for ( size_t i = 0; i < eos._edges.size(); ++i )
     {
       if ( eos._edges[i]->Is( _LayerEdge::BLOCKED )) continue;
-      eos._edges[i]->_maxLen = thinkness;
+      eos._edges[i]->SetMaxLen( thinkness );
       eos._edges[i]->FindIntersection( *searcher, intersecDist, data._epsilon, eos, &face );
       if ( intersecDist > 0 && face )
       {
         data._geomSize = Min( data._geomSize, intersecDist );
         if ( !neighborFaces.count( face->getshapeId() ))
-          eos._edges[i]->_maxLen = Min( thinkness, intersecDist / ( face->GetID() < 0 ? 3. : 2. ));
+          eos[i]->SetMaxLen( Min( thinkness, intersecDist / ( face->GetID() < 0 ? 3. : 2. )));
       }
     }
   }
+
+  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;
 }
 
 //================================================================================
@@ -4405,14 +4670,7 @@ 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;
@@ -4530,6 +4788,7 @@ bool _ViscousBuilder::inflate(_SolidData& data)
       break;
     }
 #endif
+
     // new step size
     limitStepSize( data, 0.25 * distToIntersection );
     if ( data._stepSizeNodes[0] )
@@ -4762,7 +5021,7 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data,
 
           if ( nbBad == oldBadNb  &&
                nbBad > 0 &&
-               step < stepLimit ) // smooth w/o chech of validity
+               step < stepLimit ) // smooth w/o check of validity
           {
             dumpFunctionEnd();
             dumpFunction(SMESH_Comment("smoothWoCheck")<<data._index<<"_Fa"<<sInd
@@ -4851,6 +5110,7 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data,
       _LayerEdge*      edge = eos._edges[i];
       if ( edge->_nodes.size() < 2 ) continue;
       SMESH_TNodeXYZ tgtXYZ = edge->_nodes.back();
+      //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 )
@@ -4890,6 +5150,7 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data,
   const SMDS_MeshElement* intFace = 0;
   const SMDS_MeshElement* closestFace = 0;
   _LayerEdge* le = 0;
+  bool is1stBlocked = true; // dbg
   for ( size_t iS = 0; iS < data._edgesOnShape.size(); ++iS )
   {
     _EdgesOnShape& eos = data._edgesOnShape[ iS ];
@@ -4961,7 +5222,7 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data,
           continue;
 
         // ignore intersection with intFace of an adjacent FACE
-        if ( dist > 0 )
+        if ( dist > 0.1 * eos._edges[i]->_len )
         {
           bool toIgnore = false;
           if (  eos._toSmooth )
@@ -5000,17 +5261,19 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data,
         if ( toBlockInfaltion &&
              dist < ( eos._edges[i]->_len * theThickToIntersection ))
         {
+          if ( is1stBlocked ) { is1stBlocked = false; // debug
+            dumpFunction(SMESH_Comment("blockIntersected") <<data._index<<"_InfStep"<<infStep);
+          }
           eos._edges[i]->Set( _LayerEdge::INTERSECTED ); // not to intersect
           eos._edges[i]->Block( data );                  // not to inflate
 
-          if ( _EdgesOnShape* eof = data.GetShapeEdges( intFace->getshapeId() ))
+          //if ( _EdgesOnShape* eof = data.GetShapeEdges( intFace->getshapeId() ))
           {
             // block _LayerEdge's, on top of which intFace is
             if ( const _TmpMeshFace* f = dynamic_cast< const _TmpMeshFace*>( intFace ))
             {
-              const SMDS_MeshElement* srcFace =
-                eof->_subMesh->GetSubMeshDS()->GetElement( f->getIdInShape() );
-              SMDS_ElemIteratorPtr nIt = srcFace->nodesIterator();
+              const SMDS_MeshElement* srcFace = f->_srcFace;
+              SMDS_ElemIteratorPtr        nIt = srcFace->nodesIterator();
               while ( nIt->more() )
               {
                 const SMDS_MeshNode* srcNode = static_cast<const SMDS_MeshNode*>( nIt->next() );
@@ -5033,11 +5296,14 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data,
     } // loop on eos._edges
   } // loop on data._edgesOnShape
 
+  if ( !is1stBlocked )
+    dumpFunctionEnd();
+
   if ( closestFace && le )
   {
 #ifdef __myDEBUG
     SMDS_MeshElement::iterator nIt = closestFace->begin_nodes();
-    cout << "Shortest distance: _LayerEdge nodes: tgt " << le->_nodes.back()->GetID()
+    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;
@@ -5231,7 +5497,8 @@ void _ViscousBuilder::makeOffsetSurface( _EdgesOnShape& eos, SMESH_MesherHelper&
 
   try
   {
-    BRepOffsetAPI_MakeOffsetShape offsetMaker( eos._shape, -offset, Precision::Confusion() );
+    BRepOffsetAPI_MakeOffsetShape offsetMaker;
+    offsetMaker.PerformByJoin( eos._shape, -offset, Precision::Confusion() );
     if ( !offsetMaker.IsDone() ) return;
 
     TopExp_Explorer fExp( offsetMaker.Shape(), TopAbs_FACE );
@@ -5508,8 +5775,10 @@ void _Smoother1D::findEdgesToSmooth()
   {
     if ( !_eos[i]->Is( _LayerEdge::TO_SMOOTH ))
     {
-      if ( needSmoothing( _leOnV[0]._cosin, _eos[i]->_len, _curveLen * _leParams[i] ) ||
-           isToSmooth( i ))
+      if ( needSmoothing( _leOnV[0]._cosin,
+                          _eos[i]->_len * leOnV[0]->_lenFactor, _curveLen * _leParams[i] ) ||
+           isToSmooth( i )
+           )
         _eos[i]->Set( _LayerEdge::TO_SMOOTH );
       else
         break;
@@ -5523,7 +5792,8 @@ void _Smoother1D::findEdgesToSmooth()
   {
     if ( !_eos[i]->Is( _LayerEdge::TO_SMOOTH ))
     {
-      if ( needSmoothing( _leOnV[1]._cosin, _eos[i]->_len, _curveLen * ( 1.-_leParams[i] )) ||
+      if ( needSmoothing( _leOnV[1]._cosin,
+                          _eos[i]->_len * leOnV[1]->_lenFactor, _curveLen * ( 1.-_leParams[i] )) ||
            isToSmooth( i ))
         _eos[i]->Set( _LayerEdge::TO_SMOOTH );
       else
@@ -5654,26 +5924,47 @@ bool _Smoother1D::smoothAnalyticEdge( _SolidData&                    data,
         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 ];
-        const gp_XY rangeUV = uv1 - uv0;
+        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;
-          _eos[i]->_pos.back().SetCoord( newUV.X(), newUV.Y(), 0 );
 
           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() );
+          SMDS_FacePositionPtr pos = 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;
         }
       }
     }
@@ -5761,9 +6052,11 @@ bool _Smoother1D::smoothAnalyticEdge( _SolidData&                    data,
         tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() );
         dumpMove( tgtNode );
 
-        SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( tgtNode->GetPosition() );
+        SMDS_FacePositionPtr pos = tgtNode->GetPosition();
         pos->SetUParameter( newUV.X() );
         pos->SetVParameter( newUV.Y() );
+
+        _eos[i]->Set( _LayerEdge::SMOOTHED ); // to check in refine() (IPAL54237)
       }
     }
     return true;
@@ -5836,16 +6129,19 @@ bool _Smoother1D::smoothComplexEdge( _SolidData&                    data,
   // project tgt nodes of extreme _LayerEdge's to the offset segments
   // -----------------------------------------------------------------
 
-  if ( e[0]->Is( _LayerEdge::NORMAL_UPDATED )) _iSeg[0] = 0;
-  if ( e[1]->Is( _LayerEdge::NORMAL_UPDATED )) _iSeg[1] = _offPoints.size()-2;
+  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];
+  bool isProjected[2];
   for ( int is2nd = 0; is2nd < 2; ++is2nd )
   {
     pExtreme[ is2nd ] = SMESH_TNodeXYZ( e[is2nd]->_nodes.back() );
     int  i = _iSeg[ is2nd ];
     int di = is2nd ? -1 : +1;
-    bool projected = false;
+    bool & projected = isProjected[ is2nd ];
+    projected = false;
     double uOnSeg, distMin = Precision::Infinite(), dist, distPrev = 0;
     int nbWorse = 0;
     do {
@@ -5892,38 +6188,43 @@ bool _Smoother1D::smoothComplexEdge( _SolidData&                    data,
   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]->_normal * vDiv0.XYZ() < 0 ) e[0]->_len += d0;
-  else                                   e[0]->_len -= d0;
-  if ( e[1]->_normal * vDiv1.XYZ() < 0 ) e[1]->_len += d1;
-  else                                   e[1]->_len -= d1;
+  double d1 = isProjected[1] ? vDiv1.Magnitude() : 0;
+  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
   // ---------------------------------------------------------------------------------
 
+  // temporary replace extreme _offPoints by pExtreme
+  gp_XYZ opXYZ[2] = { _offPoints[ _iSeg[0]   ]._xyz,
+                      _offPoints[ _iSeg[1]+1 ]._xyz };
+  _offPoints[ _iSeg[0]   ]._xyz = pExtreme[0].XYZ();
+  _offPoints[ _iSeg[1]+ 1]._xyz = pExtreme[1].XYZ();
+
   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 )/* * e[0]->_lenFactor*/;
+  len[ iSeg++ ] = pProj[ 0 ].Distance( _offPoints[ _iSeg[0]+1 ]._xyz );
   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 )/* * e[1]->_lenFactor*/;
+  // if ( isProjected[ 1 ])
+  //   len[ nbSeg ] -= pProj[ 1 ].Distance( _offPoints[ _iSeg[1]+1 ]._xyz );
+  // else
+  //   len[ nbSeg ] += pExtreme[ 1 ].Distance( _offPoints[ _iSeg[1]+1 ]._xyz );
 
-  // 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;
 
-  // temporary replace extreme _offPoints by pExtreme
-  gp_XYZ op[2] = { _offPoints[ _iSeg[0]   ]._xyz,
-                   _offPoints[ _iSeg[1]+1 ]._xyz };
-  _offPoints[ _iSeg[0]   ]._xyz = pExtreme[0].XYZ();
-  _offPoints[ _iSeg[1]+ 1]._xyz = pExtreme[1].XYZ();
-
   // -------------------------------------------------------------
   // distribute tgt nodes of _LayerEdge's between the projections
   // -------------------------------------------------------------
@@ -5956,8 +6257,8 @@ bool _Smoother1D::smoothComplexEdge( _SolidData&                    data,
     dumpMove( tgtNode );
   }
 
-  _offPoints[ _iSeg[0]   ]._xyz = op[0];
-  _offPoints[ _iSeg[1]+1 ]._xyz = op[1];
+  _offPoints[ _iSeg[0]   ]._xyz = opXYZ[0];
+  _offPoints[ _iSeg[1]+1 ]._xyz = opXYZ[1];
 
   return true;
 }
@@ -6011,8 +6312,8 @@ void _Smoother1D::prepare(_SolidData& data)
 
   // 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 )
   {
@@ -6022,14 +6323,39 @@ void _Smoother1D::prepare(_SolidData& data)
 
   const double u0 = c3dAdaptor.FirstParameter();
   gp_Pnt p; gp_Vec tangent;
-  _offPoints.resize( discret.NbPoints() );
-  for ( size_t i = 0; i < _offPoints.size(); i++ )
+  if ( discret.NbPoints() >= (int) _eos.size() + 2 )
+  {
+    _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
   {
-    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;
+    std::vector< double > params( _eos.size() + 2 );
+
+    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
@@ -6069,8 +6395,14 @@ void _Smoother1D::prepare(_SolidData& data)
 
   int iLBO = _offPoints.size() - 2; // last but one
 
-  _leOnV[ 0 ]._normal = getNormalNormal( leOnV[0]->_normal, _edgeDir[0] );
-  _leOnV[ 1 ]._normal = getNormalNormal( leOnV[1]->_normal, _edgeDir[1] );
+  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;
@@ -6104,7 +6436,7 @@ 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
  */
 //================================================================================
 
@@ -6115,9 +6447,36 @@ gp_XYZ _Smoother1D::getNormalNormal( const gp_XYZ & normal,
   gp_XYZ  norm = edgeDir ^ cross;
   double  size = norm.Modulus();
 
+  // if ( size == 0 ) // MULTI_NORMAL _LayerEdge
+  //   return gp_XYZ( 1e-100, 1e-100, 1e-100 );
+
   return norm / size;
 }
 
+//================================================================================
+/*!
+ * \brief Writes a script creating a mesh composed of _offPoints
+ */
+//================================================================================
+
+void _Smoother1D::offPointsToPython() const
+{
+  const char* fname = "/tmp/offPoints.py";
+  cout << "exec(open('"<<fname<<"','rb').read() )"<<endl;
+  ofstream py(fname);
+  py << "import SMESH" << endl
+     << "from salome.smesh import smeshBuilder" << endl
+     << "smesh  = smeshBuilder.New(salome.myStudy)" << endl
+     << "mesh   = smesh.Mesh( 'offPoints' )"<<endl;
+  for ( size_t i = 0; i < _offPoints.size(); i++ )
+  {
+    py << "mesh.AddNode( "
+       << _offPoints[i]._xyz.X() << ", "
+       << _offPoints[i]._xyz.Y() << ", "
+       << _offPoints[i]._xyz.Z() << " )" << endl;
+  }
+}
+
 //================================================================================
 /*!
  * \brief Sort _LayerEdge's by a parameter on a given EDGE
@@ -6225,43 +6584,43 @@ void _SolidData::PrepareEdgesToSmoothOnFace( _EdgesOnShape* eos, bool substitute
     }
 
     // SetSmooLen() to _LayerEdge's on FACE
-    for ( size_t i = 0; i < eos->_edges.size(); ++i )
-    {
-      eos->_edges[i]->SetSmooLen( Precision::Infinite() );
-    }
-    SMESH_subMeshIteratorPtr smIt = eos->_subMesh->getDependsOnIterator(/*includeSelf=*/false);
-    while ( smIt->more() ) // loop on sub-shapes of the FACE
-    {
-      _EdgesOnShape* eoe = GetShapeEdges( smIt->next()->GetId() );
-      if ( !eoe ) continue;
-
-      vector<_LayerEdge*>& eE = eoe->_edges;
-      for ( size_t iE = 0; iE < eE.size(); ++iE ) // loop on _LayerEdge's on EDGE or VERTEX
-      {
-        if ( eE[iE]->_cosin <= theMinSmoothCosin )
-          continue;
+    // for ( size_t i = 0; i < eos->_edges.size(); ++i )
+    // {
+    //   eos->_edges[i]->SetSmooLen( Precision::Infinite() );
+    // }
+    // SMESH_subMeshIteratorPtr smIt = eos->_subMesh->getDependsOnIterator(/*includeSelf=*/false);
+    // while ( smIt->more() ) // loop on sub-shapes of the FACE
+    // {
+    //   _EdgesOnShape* eoe = GetShapeEdges( smIt->next()->GetId() );
+    //   if ( !eoe ) continue;
 
-        SMDS_ElemIteratorPtr segIt = eE[iE]->_nodes[0]->GetInverseElementIterator(SMDSAbs_Edge);
-        while ( segIt->more() )
-        {
-          const SMDS_MeshElement* seg = segIt->next();
-          if ( !eos->_subMesh->DependsOn( seg->getshapeId() ))
-            continue;
-          if ( seg->GetNode(0) != eE[iE]->_nodes[0] )
-            continue; // not to check a seg twice
-          for ( size_t iN = 0; iN < eE[iE]->_neibors.size(); ++iN )
-          {
-            _LayerEdge* eN = eE[iE]->_neibors[iN];
-            if ( eN->_nodes[0]->getshapeId() != eos->_shapeID )
-              continue;
-            double dist    = SMESH_MeshAlgos::GetDistance( seg, SMESH_TNodeXYZ( eN->_nodes[0] ));
-            double smooLen = getSmoothingThickness( eE[iE]->_cosin, dist );
-            eN->SetSmooLen( Min( smooLen, eN->GetSmooLen() ));
-            eN->Set( _LayerEdge::NEAR_BOUNDARY );
-          }
-        }
-      }
-    }
+    //   vector<_LayerEdge*>& eE = eoe->_edges;
+    //   for ( size_t iE = 0; iE < eE.size(); ++iE ) // loop on _LayerEdge's on EDGE or VERTEX
+    //   {
+    //     if ( eE[iE]->_cosin <= theMinSmoothCosin )
+    //       continue;
+
+    //     SMDS_ElemIteratorPtr segIt = eE[iE]->_nodes[0]->GetInverseElementIterator(SMDSAbs_Edge);
+    //     while ( segIt->more() )
+    //     {
+    //       const SMDS_MeshElement* seg = segIt->next();
+    //       if ( !eos->_subMesh->DependsOn( seg->getshapeId() ))
+    //         continue;
+    //       if ( seg->GetNode(0) != eE[iE]->_nodes[0] )
+    //         continue; // not to check a seg twice
+    //       for ( size_t iN = 0; iN < eE[iE]->_neibors.size(); ++iN )
+    //       {
+    //         _LayerEdge* eN = eE[iE]->_neibors[iN];
+    //         if ( eN->_nodes[0]->getshapeId() != eos->_shapeID )
+    //           continue;
+    //         double dist    = SMESH_MeshAlgos::GetDistance( seg, SMESH_TNodeXYZ( eN->_nodes[0] ));
+    //         double smooLen = getSmoothingThickness( eE[iE]->_cosin, dist );
+    //         eN->SetSmooLen( Min( smooLen, eN->GetSmooLen() ));
+    //         eN->Set( _LayerEdge::NEAR_BOUNDARY );
+    //       }
+    //     }
+    //   }
+    // }
   } // if ( eos->ShapeType() == TopAbs_FACE )
 
   for ( size_t i = 0; i < eos->_edges.size(); ++i )
@@ -6299,11 +6658,12 @@ void _SolidData::PrepareEdgesToSmoothOnFace( _EdgesOnShape* eos, bool substitute
     avgLen      /= edge->_simplices.size();
     if (( edge->_curvature = _Curvature::New( avgNormProj, avgLen )))
     {
+      edge->Set( _LayerEdge::SMOOTHED_C1 );
       isCurved = true;
-      SMDS_FacePosition* fPos = dynamic_cast<SMDS_FacePosition*>( edge->_nodes[0]->GetPosition() );
+      SMDS_FacePositionPtr fPos = edge->_nodes[0]->GetPosition();
       if ( !fPos )
         for ( size_t iS = 0; iS < edge->_simplices.size()  &&  !fPos; ++iS )
-          fPos = dynamic_cast<SMDS_FacePosition*>( edge->_simplices[iS]._nPrev->GetPosition() );
+          fPos = edge->_simplices[iS]._nPrev->GetPosition();
       if ( fPos )
         edge->_curvature->_uv.SetCoord( fPos->GetUParameter(), fPos->GetVParameter() );
     }
@@ -6394,7 +6754,7 @@ void _ViscousBuilder::limitMaxLenByCurvature( _SolidData& data, SMESH_MesherHelp
           if ( eI->_nodes[0]->GetID() < eN->_nodes[0]->GetID() ) // treat this pair once
           {
             _EdgesOnShape* eosN = data.GetShapeEdges( eN );
-            limitMaxLenByCurvature( eI, eN, eosI, *eosN, helper );
+            limitMaxLenByCurvature( eI, eN, eosI, *eosN, eosI._hyp.ToSmooth() );
           }
         }
       }
@@ -6408,7 +6768,7 @@ void _ViscousBuilder::limitMaxLenByCurvature( _SolidData& data, SMESH_MesherHelp
       for ( size_t i = 1; i < eosI._edges.size(); ++i )
       {
         _LayerEdge* eI = eosI._edges[i];
-        limitMaxLenByCurvature( eI, e0, eosI, eosI, helper );
+        limitMaxLenByCurvature( eI, e0, eosI, eosI, eosI._hyp.ToSmooth() );
         e0 = eI;
       }
     }
@@ -6421,12 +6781,17 @@ void _ViscousBuilder::limitMaxLenByCurvature( _SolidData& data, SMESH_MesherHelp
  */
 //================================================================================
 
-void _ViscousBuilder::limitMaxLenByCurvature( _LayerEdge*         e1,
-                                              _LayerEdge*         e2,
-                                              _EdgesOnShape&      eos1,
-                                              _EdgesOnShape&      eos2,
-                                              SMESH_MesherHelper& helper )
+void _ViscousBuilder::limitMaxLenByCurvature( _LayerEdge*    e1,
+                                              _LayerEdge*    e2,
+                                              _EdgesOnShape& eos1,
+                                              _EdgesOnShape& eos2,
+                                              const bool     isSmoothable )
 {
+  if (( e1->_nodes[0]->GetPosition()->GetDim() !=
+        e2->_nodes[0]->GetPosition()->GetDim() ) &&
+      ( e1->_cosin < 0.75 ))
+    return; // angle > 90 deg at e1
+
   gp_XYZ plnNorm = e1->_normal ^ e2->_normal;
   double norSize = plnNorm.SquareModulus();
   if ( norSize < std::numeric_limits<double>::min() )
@@ -6446,9 +6811,10 @@ void _ViscousBuilder::limitMaxLenByCurvature( _LayerEdge*         e1,
     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 );
+    {
+      const double coef = 0.75;
+      e1->SetMaxLen( Min( e1->_maxLen, coef * u1 / e1->_lenFactor ));
+      e2->SetMaxLen( Min( e2->_maxLen, coef * u2 / e2->_lenFactor ));
     }
   }
 }
@@ -6469,6 +6835,7 @@ void _ViscousBuilder::findCollisionEdges( _SolidData& data, SMESH_MesherHelper&
     _EdgesOnShape& eos = data._edgesOnShape[iS];
     if ( eos._edges.empty() ) continue;
     if ( eos.ShapeType() != TopAbs_EDGE && eos.ShapeType() != TopAbs_VERTEX ) continue;
+    if ( !eos._sWOL.IsNull() ) continue; // PAL23566
 
     for ( size_t i = 0; i < eos._edges.size(); ++i )
     {
@@ -6515,7 +6882,7 @@ void _ViscousBuilder::findCollisionEdges( _SolidData& data, SMESH_MesherHelper&
              src2->GetID() < edge->_nodes[0]->GetID() )
           continue; // avoid using same segment twice
 
-        // a _LayerEdge containg tgt2
+        // a _LayerEdge containing tgt2
         _LayerEdge* neiborEdge = edge->_2neibors->_edges[j];
 
         _TmpMeshFaceOnEdge* f = new _TmpMeshFaceOnEdge( edge, neiborEdge, --_tmpFaceID );
@@ -6594,9 +6961,9 @@ void _ViscousBuilder::findCollisionEdges( _SolidData& data, SMESH_MesherHelper&
               ( f->_le2->IsOnEdge() && f->_le2->_2neibors->include( edge )))  continue;
         }
         dist1 = dist2 = Precision::Infinite();
-        if ( !edge->SegTriaInter( lastSegment, f->_nn[0], f->_nn[1], f->_nn[2], dist1, eps ))
+        if ( !edge->SegTriaInter( lastSegment, f->n(0), f->n(1), f->n(2), dist1, eps ))
           dist1 = Precision::Infinite();
-        if ( !edge->SegTriaInter( lastSegment, f->_nn[3], f->_nn[2], f->_nn[0], dist2, eps ))
+        if ( !edge->SegTriaInter( lastSegment, f->n(3), f->n(2), f->n(0), dist2, eps ))
           dist2 = Precision::Infinite();
         if (( dist1 > segLen ) && ( dist2 > segLen ))
           continue;
@@ -6604,7 +6971,7 @@ void _ViscousBuilder::findCollisionEdges( _SolidData& data, SMESH_MesherHelper&
         if ( edge->IsOnEdge() )
         {
           // skip perpendicular EDGEs
-          gp_Vec fSegDir  = SMESH_TNodeXYZ( f->_nn[0] ) - SMESH_TNodeXYZ( f->_nn[3] );
+          gp_Vec fSegDir  = SMESH_TNodeXYZ( f->n(0) ) - SMESH_TNodeXYZ( f->n(3) );
           bool isParallel = ( isLessAngle( eSegDir0, fSegDir, angle45 ) ||
                               isLessAngle( eSegDir1, fSegDir, angle45 ) ||
                               isLessAngle( eSegDir0, fSegDir.Reversed(), angle45 ) ||
@@ -6623,7 +6990,7 @@ void _ViscousBuilder::findCollisionEdges( _SolidData& data, SMESH_MesherHelper&
         // else
         // {
         //   double shortLen = 0.75 * ( Min( dist1, dist2 ) / edge->_lenFactor );
-        //   edge->_maxLen = Min( shortLen, edge->_maxLen );
+        //   edge->SetMaxLen( Min( shortLen, edge->_maxLen ));
         // }
       }
 
@@ -6714,7 +7081,7 @@ void _ViscousBuilder::findEdgesToUpdateNormalNearConvexFace( _ConvexFace &
       if ( dist < 0.95 * ledge->_maxLen )
       {
         ledge->Set( _LayerEdge::UPD_NORMAL_CONV );
-        if ( !ledge->_curvature ) ledge->_curvature = new _Curvature;
+        if ( !ledge->_curvature ) ledge->_curvature = _Factory::NewCurvature();
         ledge->_curvature->_uv.SetCoord( uv.X(), uv.Y() );
         edgesToUpdateFound = true;
       }
@@ -6834,8 +7201,8 @@ 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;
@@ -6873,12 +7240,17 @@ bool _ViscousBuilder::updateNormals( _SolidData&         data,
         e2neIt = edge2newEdge.insert( make_pair( edge1, zeroEdge )).first;
         e2neIt->second._normal += distWgt * newNormal;
         e2neIt->second._cosin   = newCos;
-        e2neIt->second._maxLen  = 0.7 * minIntDist / edge1->_lenFactor;
+        e2neIt->second.SetMaxLen( 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.SetMaxLen( 1.3 * minIntDist / edge1->_lenFactor );
+        }
         if ( iter > 0 && sgn1 * sgn2 < 0 && edge2->_cosin < 0 )
           e2neIt->second._normal += dir1;
       }
@@ -6895,9 +7267,10 @@ bool _ViscousBuilder::updateNormals( _SolidData&         data,
     for ( e2neIt = edge2newEdge.begin(); e2neIt != edge2newEdge.end(); ++e2neIt )
     {
       _LayerEdge*    edge = e2neIt->first;
-      if ( edge->Is( _LayerEdge::BLOCKED )) continue;
       _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();
@@ -6906,7 +7279,7 @@ bool _ViscousBuilder::updateNormals( _SolidData&         data,
         if ( newEdge._maxLen < edge->_len && iter > 0 ) // limit _maxLen
         {
           edge->InvalidateStep( stepNb + 1, *eos, /*restoreLength=*/true  );
-          edge->_maxLen = newEdge._maxLen;
+          edge->SetMaxLen( newEdge._maxLen );
           edge->SetNewLength( newEdge._maxLen, *eos, helper );
         }
         continue; // the new _normal is bad
@@ -7804,7 +8177,7 @@ bool _LayerEdge::FindIntersection( SMESH_ElementSearcher&   searcher,
 
 gp_XYZ _LayerEdge::PrevCheckPos( _EdgesOnShape* eos ) const
 {
-  size_t i = Is( NORMAL_UPDATED ) ? _pos.size()-2 : 0;
+  size_t i = Is( NORMAL_UPDATED ) && IsOnFace() ? _pos.size()-2 : 0;
 
   if ( !eos || eos->_sWOL.IsNull() )
     return _pos[ i ];
@@ -8034,7 +8407,7 @@ void _LayerEdge::MoveNearConcaVer( const _EdgesOnShape*    eov,
     prevPosV = surface.Value( prevPosV.X(), prevPosV.Y() ).XYZ();
   }
 
-  SMDS_FacePosition* fPos;
+  SMDS_FacePositionPtr fPos;
   //double r = 1. - Min( 0.9, step / 10. );
   for ( set< _LayerEdge* >::iterator e = edges.begin(); e != edges.end(); ++e )
   {
@@ -8048,9 +8421,9 @@ void _LayerEdge::MoveNearConcaVer( const _EdgesOnShape*    eov,
 
     // set _curvature to make edgeF updated by putOnOffsetSurface()
     if ( !edgeF->_curvature )
-      if (( fPos = dynamic_cast<SMDS_FacePosition*>( edgeF->_nodes[0]->GetPosition() )))
+      if (( fPos = edgeF->_nodes[0]->GetPosition() ))
       {
-        edgeF->_curvature = new _Curvature;
+        edgeF->_curvature = _Factory::NewCurvature();
         edgeF->_curvature->_r = 0;
         edgeF->_curvature->_k = 0;
         edgeF->_curvature->_h2lenRatio = 0;
@@ -8523,7 +8896,7 @@ int _LayerEdge::Smooth(const int step, const bool isConcaveFace, bool findBest )
 
 //================================================================================
 /*!
- * \brief Chooses a smoothing technic giving a position most close to an initial one.
+ * \brief Chooses a smoothing technique giving a position most close to an initial one.
  *        For a correct result, _simplices must contain nodes lying on geometry.
  */
 //================================================================================
@@ -8558,7 +8931,7 @@ void _LayerEdge::ChooseSmooFunction( const set< TGeomID >& concaveVertices,
       }
     }
 
-    // // this coice is done only if ( !concaveVertices.empty() ) for Grids/smesh/bugs_19/X1
+    // // this choice is done only if ( !concaveVertices.empty() ) for Grids/smesh/bugs_19/X1
     // // where the nodes are smoothed too far along a sphere thus creating
     // // inverted _simplices
     // double dist[theNbSmooFuns];
@@ -8734,7 +9107,7 @@ gp_XYZ _LayerEdge::smoothAngular()
 
 //================================================================================
 /*!
- * \brief Computes a new node position using weigthed node positions
+ * \brief Computes a new node position using weighted node positions
  */
 //================================================================================
 
@@ -9269,7 +9642,7 @@ void _LayerEdge::SetNewLength( double len, _EdgesOnShape& eos, SMESH_MesherHelpe
       _pos.back().SetCoord( u, 0, 0 );
       if ( _nodes.size() > 1 && uvOK )
       {
-        SMDS_EdgePosition* pos = static_cast<SMDS_EdgePosition*>( n->GetPosition() );
+        SMDS_EdgePositionPtr pos = n->GetPosition();
         pos->SetUParameter( u );
       }
     }
@@ -9281,7 +9654,7 @@ void _LayerEdge::SetNewLength( double len, _EdgesOnShape& eos, SMESH_MesherHelpe
       _pos.back().SetCoord( uv.X(), uv.Y(), 0 );
       if ( _nodes.size() > 1 && uvOK )
       {
-        SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( n->GetPosition() );
+        SMDS_FacePositionPtr pos = n->GetPosition();
         pos->SetUParameter( uv.X() );
         pos->SetVParameter( uv.Y() );
       }
@@ -9327,9 +9700,9 @@ void _LayerEdge::Block( _SolidData& data )
   SMESH_Comment msg( "#BLOCK shape=");
   msg << data.GetShapeEdges( this )->_shapeID
       << ", nodes " << _nodes[0]->GetID() << ", " << _nodes.back()->GetID();
-  dumpCmd( msg + " -- BEGIN")
+  dumpCmd( msg + " -- BEGIN");
 
-  _maxLen = _len;
+  SetMaxLen( _len );
   std::queue<_LayerEdge*> queue;
   queue.push( this );
 
@@ -9353,18 +9726,19 @@ void _LayerEdge::Block( _SolidData& data )
       double newMaxLen = edge->_maxLen + 0.5 * Sqrt( minDist );
       //if ( edge->_nodes[0]->getshapeId() == neibor->_nodes[0]->getshapeId() ) viscous_layers_00/A3
       {
-        newMaxLen *= edge->_lenFactor / neibor->_lenFactor;
+        //newMaxLen *= edge->_lenFactor / neibor->_lenFactor;
         // newMaxLen *= Min( edge->_lenFactor / neibor->_lenFactor,
         //                   neibor->_lenFactor / edge->_lenFactor );
       }
       if ( neibor->_maxLen > newMaxLen )
       {
-        neibor->_maxLen = newMaxLen;
+        neibor->SetMaxLen( newMaxLen );
         if ( neibor->_maxLen < neibor->_len )
         {
           _EdgesOnShape* eos = data.GetShapeEdges( neibor );
+          int       lastStep = neibor->Is( BLOCKED ) ? 1 : 0;
           while ( neibor->_len > neibor->_maxLen &&
-                  neibor->NbSteps() > 1 )
+                  neibor->NbSteps() > lastStep )
             neibor->InvalidateStep( neibor->NbSteps(), *eos, /*restoreLength=*/true );
           neibor->SetNewLength( neibor->_maxLen, *eos, data.GetHelper() );
           //neibor->Block( data );
@@ -9373,7 +9747,7 @@ void _LayerEdge::Block( _SolidData& data )
       }
     }
   }
-  dumpCmd( msg + " -- END")
+  dumpCmd( msg + " -- END");
 }
 
 //================================================================================
@@ -9396,7 +9770,7 @@ void _LayerEdge::InvalidateStep( size_t curStep, const _EdgesOnShape& eos, bool
       TopLoc_Location loc;
       if ( eos.SWOLType() == TopAbs_EDGE )
       {
-        SMDS_EdgePosition* pos = static_cast<SMDS_EdgePosition*>( n->GetPosition() );
+        SMDS_EdgePositionPtr pos = n->GetPosition();
         pos->SetUParameter( nXYZ.X() );
         double f,l;
         Handle(Geom_Curve) curve = BRep_Tool::Curve( TopoDS::Edge( eos._sWOL ), loc, f,l);
@@ -9404,7 +9778,7 @@ void _LayerEdge::InvalidateStep( size_t curStep, const _EdgesOnShape& eos, bool
       }
       else
       {
-        SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( n->GetPosition() );
+        SMDS_FacePositionPtr pos = n->GetPosition();
         pos->SetUParameter( nXYZ.X() );
         pos->SetVParameter( nXYZ.Y() );
         Handle(Geom_Surface) surface = BRep_Tool::Surface( TopoDS::Face(eos._sWOL), loc );
@@ -9416,9 +9790,15 @@ void _LayerEdge::InvalidateStep( size_t curStep, const _EdgesOnShape& eos, bool
 
     if ( restoreLength )
     {
-      _len -= ( nXYZ.XYZ() - curXYZ ).Modulus() / _lenFactor;
+      if ( NbSteps() == 0 )
+        _len = 0.;
+      else if ( IsOnFace() && Is( MOVED ))
+        _len = ( nXYZ.XYZ() - SMESH_NodeXYZ( _nodes[0] )) * _normal;
+      else
+        _len -= ( nXYZ.XYZ() - curXYZ ).Modulus() / _lenFactor;
     }
   }
+  return;
 }
 
 //================================================================================
@@ -9454,10 +9834,23 @@ void _LayerEdge::SmoothPos( const vector< double >& segLen, const double tol )
   int iSmoothed = GetSmoothedPos( tol );
   if ( !iSmoothed ) return;
 
-  //if ( 1 || Is( DISTORTED ))
+  gp_XYZ normal = _normal;
+  if ( Is( NORMAL_UPDATED ))
   {
-    gp_XYZ normal = _normal;
-    if ( Is( NORMAL_UPDATED ))
+    double minDot = 1;
+    for ( size_t i = 0; i < _neibors.size(); ++i )
+    {
+      if ( _neibors[i]->IsOnFace() )
+      {
+        double dot = _normal * _neibors[i]->_normal;
+        if ( dot < minDot )
+        {
+          normal = _neibors[i]->_normal;
+          minDot = dot;
+        }
+      }
+    }
+    if ( minDot == 1. )
       for ( size_t i = 1; i < _pos.size(); ++i )
       {
         normal = _pos[i] - _pos[0];
@@ -9468,42 +9861,29 @@ void _LayerEdge::SmoothPos( const vector< double >& segLen, const double tol )
           break;
         }
       }
-    const double r = 0.2;
-    for ( int iter = 0; iter < 50; ++iter )
+  }
+  const double r = 0.2;
+  for ( int iter = 0; iter < 50; ++iter )
+  {
+    double minDot = 1;
+    for ( size_t i = Max( 1, iSmoothed-1-iter ); i < _pos.size()-1; ++i )
     {
-      double minDot = 1;
-      for ( size_t i = Max( 1, iSmoothed-1-iter ); i < _pos.size()-1; ++i )
-      {
-        gp_XYZ midPos = 0.5 * ( _pos[i-1] + _pos[i+1] );
-        gp_XYZ newPos = ( 1-r ) * midPos + r * _pos[i];
-        _pos[i] = newPos;
-        double midLen = 0.5 * ( segLen[i-1] + segLen[i+1] );
-        double newLen = ( 1-r ) * midLen + r * segLen[i];
-        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.SquareModulus();
-        if ( size > RealSmall() )
-          minDot = Min( minDot, ( normal * posDir ) * ( normal * posDir ) / size );
-      }
-      if ( minDot > 0.5 * 0.5 )
-        break;
+      gp_XYZ midPos = 0.5 * ( _pos[i-1] + _pos[i+1] );
+      gp_XYZ newPos = ( 1-r ) * midPos + r * _pos[i];
+      _pos[i] = newPos;
+      double midLen = 0.5 * ( segLen[i-1] + segLen[i+1] );
+      double newLen = ( 1-r ) * midLen + r * segLen[i];
+      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.SquareModulus();
+      if ( size > RealSmall() )
+        minDot = Min( minDot, ( normal * posDir ) * ( normal * posDir ) / size );
     }
+    if ( minDot > 0.5 * 0.5 )
+      break;
   }
-  // else
-  // {
-  //   for ( size_t i = 1; i < _pos.size()-1; ++i )
-  //   {
-  //     if ((int) i < iSmoothed  &&  ( segLen[i] / segLen.back() < 0.5 ))
-  //       continue;
-
-  //     double     wgt = segLen[i] / segLen.back();
-  //     gp_XYZ normPos = _pos[0] + _normal * wgt * _len;
-  //     gp_XYZ tgtPos  = ( 1 - wgt ) * _pos[0] +  wgt * _pos.back();
-  //     gp_XYZ newPos  = ( 1 - wgt ) * normPos +  wgt * tgtPos;
-  //     _pos[i] = newPos;
-  //   }
-  // }
+  return;
 }
 
 //================================================================================
@@ -9544,11 +9924,11 @@ std::string _LayerEdge::DumpFlags() const
   return dump;
 }
 
+
 //================================================================================
 /*!
-  case brief:
-  default:
-*/
+ * \brief Create layers of prisms
+ */
 //================================================================================
 
 bool _ViscousBuilder::refine(_SolidData& data)
@@ -9603,11 +9983,8 @@ bool _ViscousBuilder::refine(_SolidData& data)
       surface  = helper.GetSurface( geomFace );
       // propagate _toSmooth back to _eosC1, which was unset in findShapesToSmooth()
       for ( size_t i = 0; i < eos._eosC1.size(); ++i )
-      {
         eos._eosC1[ i ]->_toSmooth = true;
-        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;
@@ -9710,8 +10087,20 @@ bool _ViscousBuilder::refine(_SolidData& data)
       }
       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
       {
@@ -9754,12 +10143,12 @@ bool _ViscousBuilder::refine(_SolidData& data)
         SMDS_PositionPtr  lastPos = tgtNode->GetPosition();
         if ( isOnEdge )
         {
-          SMDS_EdgePosition* epos = static_cast<SMDS_EdgePosition*>( lastPos );
+          SMDS_EdgePositionPtr epos = lastPos;
           epos->SetUParameter( otherTgtPos.X() );
         }
         else
         {
-          SMDS_FacePosition* fpos = static_cast<SMDS_FacePosition*>( lastPos );
+          SMDS_FacePositionPtr fpos = lastPos;
           fpos->SetUParameter( otherTgtPos.X() );
           fpos->SetVParameter( otherTgtPos.Y() );
         }
@@ -9848,7 +10237,7 @@ bool _ViscousBuilder::refine(_SolidData& data)
                 u = helper.GetNodeU( geomEdge, node );
               pos = curve->Value( u ).Transformed(loc);
 
-              SMDS_EdgePosition* epos = static_cast<SMDS_EdgePosition*>( node->GetPosition() );
+              SMDS_EdgePositionPtr epos = node->GetPosition();
               epos->SetUParameter( u );
             }
             else
@@ -9858,7 +10247,7 @@ bool _ViscousBuilder::refine(_SolidData& data)
                 uv = helper.GetNodeUV( geomFace, node );
               pos = surface->Value( uv );
 
-              SMDS_FacePosition* fpos = static_cast<SMDS_FacePosition*>( node->GetPosition() );
+              SMDS_FacePositionPtr fpos = node->GetPosition();
               fpos->SetUParameter( uv.X() );
               fpos->SetVParameter( uv.Y() );
             }
@@ -9906,6 +10295,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();
@@ -9936,14 +10330,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 )
         {
@@ -9956,16 +10356,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
@@ -9974,10 +10376,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 )
         {
@@ -9996,9 +10401,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 );
           }
@@ -10006,15 +10411,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 );
           }
           }
@@ -10025,6 +10430,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
 
@@ -10033,16 +10443,462 @@ bool _ViscousBuilder::refine(_SolidData& data)
     SMESH_ComputeErrorPtr& err = _mesh->GetSubMesh( data._solid )->GetComputeError();
     if ( !err || err->IsOK() )
     {
-      err.reset( new SMESH_ComputeError( COMPERR_WARNING,
-                                         "Bad quality volumes created" ));
-      err->myBadElements.insert( err->myBadElements.end(),
-                                 degenVols.begin(),degenVols.end() );
+      SMESH_BadInputElements* badElems =
+        new SMESH_BadInputElements( getMeshDS(), COMPERR_WARNING, "Bad quality volumes created" );
+      badElems->myBadElements.insert( badElems->myBadElements.end(),
+                                      degenVols.begin(),degenVols.end() );
+      err.reset( badElems );
     }
   }
 
   return true;
 }
 
+namespace VISCOUS_3D
+{
+  struct ShrinkFace;
+  //--------------------------------------------------------------------------------
+  /*!
+   * \brief Pair of periodic FACEs
+   */
+  struct PeriodicFaces
+  {
+    typedef StdMeshers_ProjectionUtils::TrsfFinder3D Trsf;
+
+    ShrinkFace*  _shriFace[2];
+    TNodeNodeMap _nnMap;
+    Trsf         _trsf;
+
+    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(); }
+  };
+
+  //--------------------------------------------------------------------------------
+  /*!
+   * \brief Shrink FACE data used to find periodic FACEs
+   */
+  struct ShrinkFace
+  {
+    // ................................................................................
+    struct BndPart //!< part of FACE boundary, either shrink or no-shrink
+    {
+      bool                         _isShrink, _isReverse;
+      int                          _nbSegments;
+      AverageHyp*                  _hyp;
+      std::vector< SMESH_NodeXYZ > _nodes;
+      TopAbs_ShapeEnum             _vertSWOLType[2]; // shrink part includes VERTEXes
+      AverageHyp*                  _vertHyp[2];
+
+      BndPart():
+        _isShrink(0), _isReverse(0), _nbSegments(0), _hyp(0),
+        _vertSWOLType{ TopAbs_WIRE, TopAbs_WIRE }, _vertHyp{ 0, 0 }
+      {}
+
+      bool operator==( const BndPart& other ) const
+      {
+        return ( _isShrink       == other._isShrink &&
+                 _nbSegments     == other._nbSegments &&
+                 _nodes.size()   == other._nodes.size() &&
+                 vertSWOLType1() == other.vertSWOLType1() &&
+                 vertSWOLType2() == other.vertSWOLType2() &&
+                 (( !_isShrink ) ||
+                  ( *_hyp        == *other._hyp &&
+                    vertHyp1()   == other.vertHyp1() &&
+                    vertHyp2()   == other.vertHyp2() ))
+                 );
+      }
+      bool CanAppend( const BndPart& other )
+      {
+        return ( _isShrink  == other._isShrink  &&
+                 (( !_isShrink ) ||
+                  ( *_hyp        == *other._hyp &&
+                    *_hyp        == vertHyp2()  &&
+                    vertHyp2()   == other.vertHyp1() ))
+                 );
+      }
+      void Append( const BndPart& other )
+      {
+        _nbSegments += other._nbSegments;
+        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];
+      }
+      const SMDS_MeshNode*    Node(size_t i)  const
+      {
+        return _nodes[ _isReverse ? ( _nodes.size() - 1 - i ) : i ]._node;
+      }
+      void Reverse() { _isReverse = !_isReverse; }
+      const TopAbs_ShapeEnum& vertSWOLType1() const { return _vertSWOLType[ _isReverse  ]; }
+      const TopAbs_ShapeEnum& vertSWOLType2() const { return _vertSWOLType[ !_isReverse ]; }
+      const AverageHyp&       vertHyp1()      const { return *(_vertHyp[ _isReverse  ]); }
+      const AverageHyp&       vertHyp2()      const { return *(_vertHyp[ !_isReverse ]); }
+    };
+    // ................................................................................
+
+    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;
+    }
+    bool IsSame( const TopoDS_Face& face ) const
+    {
+      return _subMesh->GetSubShape().IsSame( face );
+    }
+    bool IsShrunk( const TopTools_MapOfShape& shrunkFaces ) const
+    {
+      return shrunkFaces.Contains( _subMesh->GetSubShape() );
+    }
+
+    //================================================================================
+    /*!
+     * Check if meshes on two FACEs are equal
+     */
+    bool IsPeriodic( ShrinkFace& other, PeriodicFaces& periodic )
+    {
+      if ( !IsSameNbElements( other ))
+        return false;
+
+      this->SetBoundary();
+      other.SetBoundary();
+      if ( this->_boundarySize    != other._boundarySize ||
+           this->_nbBoundaryParts != other._nbBoundaryParts )
+        return false;
+
+      for ( int isReverse = 0; isReverse < 2; ++isReverse )
+      {
+        if ( isReverse )
+          Reverse( _boundary );
+
+        // check boundaries
+        bool equalBoundary = false;
+        for ( int iP = 0; iP < _nbBoundaryParts &&  !equalBoundary; ++iP )
+        {
+          if ( ! ( equalBoundary = ( this->_boundary == other._boundary )))
+            // set first part at end
+            _boundary.splice( _boundary.end(), _boundary, _boundary.begin() );
+        }
+        if ( !equalBoundary )
+          continue;
+
+        // check connectivity
+        std::set<const SMDS_MeshElement*> elemsThis, elemsOther;
+        this->GetElements( elemsThis  );
+        other.GetElements( elemsOther );
+        SMESH_MeshEditor::Sew_Error err =
+          SMESH_MeshEditor::FindMatchingNodes( elemsThis, elemsOther,
+                                               this->_boundary.front().Node(0),
+                                               other._boundary.front().Node(0),
+                                               this->_boundary.front().Node(1),
+                                               other._boundary.front().Node(1),
+                                               periodic._nnMap );
+        if ( err != SMESH_MeshEditor::SEW_OK )
+          continue;
+
+        // check node positions
+        std::vector< gp_XYZ > srcPnts, tgtPnts;
+        this->GetBoundaryPoints( srcPnts );
+        other.GetBoundaryPoints( tgtPnts );
+        if ( !periodic._trsf.Solve( srcPnts, tgtPnts )) {
+          continue;
+        }
+        double tol = std::numeric_limits<double>::max();
+        for ( size_t i = 1; i < srcPnts.size(); ++i ) {
+          tol = Min( tol, ( srcPnts[i-1] - srcPnts[i] ).SquareModulus() );
+        }
+        tol = 0.01 * Sqrt( tol );
+        bool nodeCoincide = true;
+        TNodeNodeMap::iterator n2n = periodic._nnMap.begin();
+        for ( ; n2n != periodic._nnMap.end() &&  nodeCoincide; ++n2n )
+        {
+          SMESH_NodeXYZ nSrc = n2n->first;
+          SMESH_NodeXYZ nTgt = n2n->second;
+          gp_XYZ pTgt = periodic._trsf.Transform( nSrc );
+          nodeCoincide = (( pTgt - nTgt ).SquareModulus() < tol );
+        }
+        if ( nodeCoincide )
+          return true;
+      }
+      return false;
+    }
+
+    bool IsSameNbElements( ShrinkFace& other ) // check number of mesh faces
+    {
+      SMESHDS_SubMesh* sm1 = this->_subMesh->GetSubMeshDS();
+      SMESHDS_SubMesh* sm2 = other._subMesh->GetSubMeshDS();
+      return ( sm1->NbElements() == sm2->NbElements() &&
+               sm1->NbNodes()    == sm2->NbNodes() );
+    }
+
+    void Reverse( std::list< BndPart >& boundary )
+    {
+      boundary.reverse();
+      for ( std::list< BndPart >::iterator part = boundary.begin(); part != boundary.end(); ++part )
+        part->Reverse();
+    }
+
+    void SetBoundary()
+    {
+      if ( !_boundary.empty() )
+        return;
+
+      TopoDS_Face F = TopoDS::Face( _subMesh->GetSubShape() );
+      if ( F.Orientation() >= TopAbs_INTERNAL ) F.Orientation( TopAbs_FORWARD );
+      std::list< TopoDS_Edge > edges;
+      std::list< int > nbEdgesInWire;
+      /*int nbWires =*/ SMESH_Block::GetOrderedEdges (F, edges, nbEdgesInWire);
+
+      // std::list< TopoDS_Edge >::iterator edgesEnd = edges.end();
+      // if ( nbWires > 1 ) {
+      //   edgesEnd = edges.begin();
+      //   std::advance( edgesEnd, nbEdgesInWire.front() );
+      // }
+      StdMeshers_FaceSide fSide( F, edges, _subMesh->GetFather(),
+                                 /*fwd=*/true, /*skipMedium=*/true );
+      _boundarySize = fSide.NbSegments();
+
+      //TopoDS_Vertex vv[2];
+      //std::list< TopoDS_Edge >::iterator edgeIt = edges.begin();
+      for ( int iE = 0; iE < nbEdgesInWire.front(); ++iE )
+      {
+        BndPart bndPart;
+        _EdgesOnShape*  eos = _data1->GetShapeEdges( fSide.EdgeID( iE ));
+
+        bndPart._isShrink = ( eos->SWOLType() == TopAbs_FACE );
+        if ( bndPart._isShrink )
+          if ((           _data1->_noShrinkShapes.count( eos->_shapeID )) ||
+              ( _data2 && _data2->_noShrinkShapes.count( eos->_shapeID )))
+            bndPart._isShrink = false;
+
+        if ( bndPart._isShrink )
+        {
+          bndPart._hyp = & eos->_hyp;
+          _EdgesOnShape* eov[2] = { _data1->GetShapeEdges( fSide.FirstVertex( iE )),
+                                    _data1->GetShapeEdges( fSide.LastVertex ( iE )) };
+          for ( int iV = 0; iV < 2; ++iV )
+          {
+            bndPart._vertHyp     [iV] = & eov[iV]->_hyp;
+            bndPart._vertSWOLType[iV] = eov[iV]->SWOLType();
+            if ( _data1->_noShrinkShapes.count( eov[iV]->_shapeID ))
+              bndPart._vertSWOLType[iV] = TopAbs_SHAPE;
+            if ( _data2 && bndPart._vertSWOLType[iV] != TopAbs_SHAPE )
+            {
+              eov[iV] = _data2->GetShapeEdges( iV ? fSide.LastVertex(iE) : fSide.FirstVertex(iE ));
+              if ( _data2->_noShrinkShapes.count( eov[iV]->_shapeID ))
+                bndPart._vertSWOLType[iV] = TopAbs_SHAPE;
+              else if ( eov[iV]->SWOLType() > bndPart._vertSWOLType[iV] )
+                bndPart._vertSWOLType[iV] = eov[iV]->SWOLType();
+            }
+          }
+        }
+        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 );
+        else
+          _boundary.back().Append( bndPart );
+      }
+
+      _nbBoundaryParts = _boundary.size();
+      if ( _nbBoundaryParts > 1 && _boundary.front()._isShrink == _boundary.back()._isShrink )
+      {
+        _boundary.back().Append( _boundary.front() );
+        _boundary.pop_front();
+        --_nbBoundaryParts;
+      }
+    }
+
+    void GetElements( std::set<const SMDS_MeshElement*>& theElems)
+    {
+      if ( SMESHDS_SubMesh* sm = _subMesh->GetSubMeshDS() )
+        for ( SMDS_ElemIteratorPtr fIt = sm->GetElements(); fIt->more(); )
+          theElems.insert( theElems.end(), fIt->next() );
+
+      return ;
+    }
+
+    void GetBoundaryPoints( std::vector< gp_XYZ >& points )
+    {
+      points.reserve( _boundarySize );
+      size_t  nb = _boundary.rbegin()->_nodes.size();
+      int lastID = _boundary.rbegin()->Node( nb - 1 )->GetID();
+      std::list< BndPart >::const_iterator part = _boundary.begin();
+      for ( ; part != _boundary.end(); ++part )
+      {
+        size_t nb = part->_nodes.size();
+        size_t iF = 0;
+        size_t iR = nb - 1;
+        size_t* i = part->_isReverse ? &iR : &iF;
+        if ( part->_nodes[ *i ]->GetID() == lastID )
+          ++iF, --iR;
+        for ( ; iF < nb; ++iF, --iR )
+          points.push_back( part->_nodes[ *i ]);
+        --iF, ++iR;
+        lastID = part->_nodes[ *i ]->GetID();
+      }
+    }
+  }; // struct ShrinkFace
+
+  //--------------------------------------------------------------------------------
+  /*!
+   * \brief Periodic FACEs
+   */
+  struct Periodicity
+  {
+    std::vector< ShrinkFace >    _shrinkFaces;
+    std::vector< PeriodicFaces > _periodicFaces;
+
+    PeriodicFaces* GetPeriodic( const TopoDS_Face& face, const TopTools_MapOfShape& shrunkFaces )
+    {
+      for ( size_t i = 0; i < _periodicFaces.size(); ++i )
+        if ( _periodicFaces[ i ].IncludeShrunk( face, shrunkFaces ))
+          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();
+    }
+  };
+
+  //================================================================================
+  /*!
+   * Check if a pair includes the given FACE and the other FACE is already shrunk
+   */
+  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 )));
+  }
+
+  //================================================================================
+  /*!
+   * Make equal meshes on periodic faces by moving corresponding nodes
+   */
+  bool PeriodicFaces::MoveNodes( const TopoDS_Face& tgtFace )
+  {
+    int iTgt = _shriFace[1]->IsSame( tgtFace );
+    int iSrc = 1 - iTgt;
+
+    _SolidData* dataSrc = _shriFace[iSrc]->_data1;
+    _SolidData* dataTgt = _shriFace[iTgt]->_data1;
+
+    Trsf * trsf = & _trsf, trsfInverse;
+    if ( iSrc != 0 )
+    {
+      trsfInverse = _trsf;
+      if ( !trsfInverse.Invert())
+        return false;
+      trsf = &trsfInverse;
+    }
+    SMESHDS_Mesh* meshDS = dataSrc->GetHelper().GetMeshDS();
+
+    TNode2Edge::iterator n2e;
+    TNodeNodeMap::iterator n2n = _nnMap.begin();
+    for ( ; n2n != _nnMap.end(); ++n2n )
+    {
+      const SMDS_MeshNode* const* nn = & n2n->first;
+      const SMDS_MeshNode*      nSrc = nn[ iSrc ];
+      const SMDS_MeshNode*      nTgt = nn[ iTgt ];
+
+      if (( nSrc->GetPosition()->GetDim() == 2 ) ||
+          (( n2e = dataSrc->_n2eMap.find( nSrc )) == dataSrc->_n2eMap.end() ))
+      {
+        SMESH_NodeXYZ pSrc = nSrc;
+        gp_XYZ pTgt = trsf->Transform( pSrc );
+        meshDS->MoveNode( nTgt, pTgt.X(), pTgt.Y(), pTgt.Z() );
+      }
+      else
+      {
+        _LayerEdge* leSrc = n2e->second;
+        n2e = dataTgt->_n2eMap.find( nTgt );
+        if ( n2e == dataTgt->_n2eMap.end() )
+          break;
+        _LayerEdge* leTgt = n2e->second;
+        if ( leSrc->_nodes.size() != leTgt->_nodes.size() )
+          break;
+        for ( size_t iN = 1; iN < leSrc->_nodes.size(); ++iN )
+        {
+          SMESH_NodeXYZ pSrc = leSrc->_nodes[ iN ];
+          gp_XYZ pTgt = trsf->Transform( pSrc );
+          meshDS->MoveNode( leTgt->_nodes[ iN ], pTgt.X(), pTgt.Y(), pTgt.Z() );
+        }
+      }
+    }
+    bool done = ( n2n == _nnMap.end() );
+    debugMsg( "PeriodicFaces::MoveNodes "
+              << _shriFace[iSrc]->_subMesh->GetId() << " -> "
+              << _shriFace[iTgt]->_subMesh->GetId() << " -- "
+              << ( done ? "DONE" : "FAIL"));
+
+    return done;
+  }
+} // namespace VISCOUS_3D; Periodicity part
+
+
+//================================================================================
+/*!
+ * \brief Find FACEs to shrink, that are equally meshed before shrink (i.e. periodic)
+ *        and should remain equal after shrink
+ */
+//================================================================================
+
+void _ViscousBuilder::findPeriodicFaces()
+{
+  // make map of (ids of FACEs to shrink mesh on) to (list of _SolidData containing
+  // _LayerEdge's inflated along FACE or EDGE)
+  std::map< TGeomID, std::list< _SolidData* > > id2sdMap;
+  for ( size_t i = 0 ; i < _sdVec.size(); ++i )
+  {
+    _SolidData& data = _sdVec[i];
+    std::map< TGeomID, TopoDS_Shape >::iterator s2s = data._shrinkShape2Shape.begin();
+    for (; s2s != data._shrinkShape2Shape.end(); ++s2s )
+      if ( s2s->second.ShapeType() == TopAbs_FACE )
+        id2sdMap[ getMeshDS()->ShapeToIndex( s2s->second )].push_back( &data );
+  }
+
+  _periodicity.reset( new Periodicity );
+  _periodicity->_shrinkFaces.resize( id2sdMap.size() );
+
+  std::map< TGeomID, std::list< _SolidData* > >::iterator id2sdIt = id2sdMap.begin();
+  for ( size_t i = 0; i < id2sdMap.size(); ++i, ++id2sdIt )
+  {
+    _SolidData* sd1 = id2sdIt->second.front();
+    _SolidData* sd2 = id2sdIt->second.back();
+    _periodicity->_shrinkFaces[ i ].Init( _mesh->GetSubMeshContaining( id2sdIt->first ), sd1, sd2 );
+  }
+
+  for (   size_t i1 = 0;      i1 < _periodicity->_shrinkFaces.size(); ++i1 )
+    for ( size_t i2 = i1 + 1; i2 < _periodicity->_shrinkFaces.size(); ++i2 )
+    {
+      PeriodicFaces pf( & _periodicity->_shrinkFaces[ i1 ],
+                        & _periodicity->_shrinkFaces[ i2 ]);
+      if ( pf._shriFace[0]->IsPeriodic( *pf._shriFace[1], pf ))
+      {
+        _periodicity->_periodicFaces.push_back( pf );
+      }
+    }
+  return;
+}
+
 //================================================================================
 /*!
  * \brief Shrink 2D mesh on faces to let space for inflated layers
@@ -10059,11 +10915,11 @@ bool _ViscousBuilder::shrink(_SolidData& theData)
     _SolidData& data = _sdVec[i];
     map< TGeomID, TopoDS_Shape >::iterator s2s = data._shrinkShape2Shape.begin();
     for (; s2s != data._shrinkShape2Shape.end(); ++s2s )
-      if ( s2s->second.ShapeType() == TopAbs_FACE && !_shrinkedFaces.Contains( s2s->second ))
+      if ( s2s->second.ShapeType() == TopAbs_FACE && !_shrunkFaces.Contains( s2s->second ))
       {
         f2sdMap[ getMeshDS()->ShapeToIndex( s2s->second )].push_back( &data );
 
-        // Put mesh faces on the shrinked FACE to the proxy sub-mesh to avoid
+        // Put mesh faces on the shrunk 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 ))
@@ -10102,7 +10958,7 @@ bool _ViscousBuilder::shrink(_SolidData& theData)
   vector< _EdgesOnShape* > subEOS;
   vector< _LayerEdge* > lEdges;
 
-  // loop on FACEs to srink mesh on
+  // loop on FACEs to shrink mesh on
   map< TGeomID, list< _SolidData* > >::iterator f2sd = f2sdMap.begin();
   for ( ; f2sd != f2sdMap.end(); ++f2sd )
   {
@@ -10122,15 +10978,24 @@ bool _ViscousBuilder::shrink(_SolidData& theData)
 
     Handle(Geom_Surface) surface = BRep_Tool::Surface( F );
 
-    _shrinkedFaces.Add( F );
+    _shrunkFaces.Add( F );
     helper.SetSubShape( F );
 
+    // ==============================
+    // Use periodicity to move nodes
+    // ==============================
+
+    PeriodicFaces* periodic = _periodicity->GetPeriodic( F, _shrunkFaces );
+    bool movedByPeriod = ( periodic && periodic->MoveNodes( F ));
+
     // ===========================
     // Prepare data for shrinking
     // ===========================
 
     // Collect nodes to smooth (they are marked at the beginning of this method)
     vector < const SMDS_MeshNode* > smoothNodes;
+
+    if ( !movedByPeriod )
     {
       SMDS_NodeIteratorPtr nIt = smDS->GetNodes();
       while ( nIt->more() )
@@ -10176,11 +11041,12 @@ bool _ViscousBuilder::shrink(_SolidData& theData)
           }
         subEOS.push_back( eos );
 
-        for ( size_t i = 0; i < eos->_edges.size(); ++i )
-        {
-          lEdges.push_back( eos->_edges[ i ] );
-          prepareEdgeToShrink( *eos->_edges[ i ], *eos, helper, smDS );
-        }
+        if ( !movedByPeriod )
+          for ( size_t i = 0; i < eos->_edges.size(); ++i )
+          {
+            lEdges.push_back( eos->_edges[ i ] );
+            prepareEdgeToShrink( *eos->_edges[ i ], *eos, helper, smDS );
+          }
       }
     }
 
@@ -10252,13 +11118,16 @@ bool _ViscousBuilder::shrink(_SolidData& theData)
         if ( eos.SWOLType() == TopAbs_EDGE )
         {
           SMESH_subMesh* edgeSM = _mesh->GetSubMesh( eos._sWOL );
-          _Shrinker1D& srinker  = e2shrMap[ edgeSM->GetId() ];
-          eShri1D.insert( & srinker );
-          srinker.AddEdge( eos._edges[0], eos, helper );
           VISCOUS_3D::ToClearSubWithMain( edgeSM, data._solid );
-          // restore params of nodes on EGDE if the EDGE has been already
-          // srinked while srinking other FACE
-          srinker.RestoreParams();
+          if ( !movedByPeriod )
+          {
+            _Shrinker1D& shrinker = e2shrMap[ edgeSM->GetId() ];
+            eShri1D.insert( & shrinker );
+            shrinker.AddEdge( eos._edges[0], eos, helper );
+            // restore params of nodes on EDGE if the EDGE has been already
+            // shrunk while shrinking other FACE
+            shrinker.RestoreParams();
+          }
         }
         for ( size_t i = 0; i < eos._edges.size(); ++i )
         {
@@ -10273,7 +11142,7 @@ bool _ViscousBuilder::shrink(_SolidData& theData)
     }
 
     bool toFixTria = false; // to improve quality of trias by diagonal swap
-    if ( isConcaveFace )
+    if ( isConcaveFace && !movedByPeriod )
     {
       const bool hasTria = _mesh->NbTriangles(), hasQuad = _mesh->NbQuadrangles();
       if ( hasTria != hasQuad ) {
@@ -10292,24 +11161,24 @@ bool _ViscousBuilder::shrink(_SolidData& theData)
     // Perform shrinking
     // ==================
 
-    bool shrinked = true;
+    bool shrunk = !movedByPeriod;
     int nbBad, shriStep=0, smooStep=0;
     _SmoothNode::SmoothType smoothType
       = isConcaveFace ? _SmoothNode::ANGULAR : _SmoothNode::LAPLACIAN;
     SMESH_Comment errMsg;
-    while ( shrinked )
+    while ( shrunk )
     {
       shriStep++;
       // Move boundary nodes (actually just set new UV)
       // -----------------------------------------------
       dumpFunction(SMESH_Comment("moveBoundaryOnF")<<f2sd->first<<"_st"<<shriStep ); // debug
-      shrinked = false;
+      shrunk = false;
       for ( size_t iS = 0; iS < subEOS.size(); ++iS )
       {
         _EdgesOnShape& eos = * subEOS[ iS ];
         for ( size_t i = 0; i < eos._edges.size(); ++i )
         {
-          shrinked |= eos._edges[i]->SetNewLength2d( surface, F, eos, helper );
+          shrunk |= eos._edges[i]->SetNewLength2d( surface, F, eos, helper );
         }
       }
       dumpFunctionEnd();
@@ -10400,7 +11269,7 @@ bool _ViscousBuilder::shrink(_SolidData& theData)
       //   }
       // }
 
-    } // while ( shrinked )
+    } // while ( shrunk )
 
     if ( !errMsg.empty() ) // Try to re-compute the shrink FACE
     {
@@ -10438,6 +11307,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 )
@@ -10449,14 +11320,14 @@ bool _ViscousBuilder::shrink(_SolidData& theData)
                edge->Is( _LayerEdge::SHRUNK )) continue;
           if ( subEOS[iS]->SWOLType() == TopAbs_FACE )
           {
-            SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( tgtNode->GetPosition() );
+            SMDS_FacePositionPtr pos = tgtNode->GetPosition();
             pos->SetUParameter( edge->_pos[0].X() );
             pos->SetVParameter( edge->_pos[0].Y() );
             p = surface->Value( edge->_pos[0].X(), edge->_pos[0].Y() );
           }
           else
           {
-            SMDS_EdgePosition* pos = static_cast<SMDS_EdgePosition*>( tgtNode->GetPosition() );
+            SMDS_EdgePositionPtr pos = tgtNode->GetPosition();
             pos->SetUParameter( edge->_pos[0].Coord( U_TGT ));
             p = BRepAdaptor_Curve( TopoDS::Edge( subEOS[iS]->_sWOL )).Value( pos->GetUParameter() );
           }
@@ -10564,7 +11435,7 @@ bool _ViscousBuilder::shrink(_SolidData& theData)
         return error( errMsg );
 
     } // end of re-meshing in case of failed smoothing
-    else
+    else if ( !movedByPeriod )
     {
       // No wrongly shaped faces remain; final smooth. Set node XYZ.
       bool isStructuredFixed = false;
@@ -10605,10 +11476,10 @@ bool _ViscousBuilder::shrink(_SolidData& theData)
     if ( data2 )
       VISCOUS_3D::ToClearSubWithMain( sm, data2->_solid );
 
-  } // loop on FACES to srink mesh on
+  } // loop on FACES to shrink mesh on
 
 
-  // Replace source nodes by target nodes in shrinked mesh edges
+  // Replace source nodes by target nodes in shrunk mesh edges
 
   map< int, _Shrinker1D >::iterator e2shr = e2shrMap.begin();
   for ( ; e2shr != e2shrMap.end(); ++e2shr )
@@ -10651,7 +11522,7 @@ bool _ViscousBuilder::prepareEdgeToShrink( _LayerEdge&            edge,
     edge._pos[0].SetCoord( tgtUV.X(), tgtUV.Y(), 0 );
 
     // set UV of source node to target node
-    SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( tgtNode->GetPosition() );
+    SMDS_FacePositionPtr pos = tgtNode->GetPosition();
     pos->SetUParameter( srcUV.X() );
     pos->SetVParameter( srcUV.Y() );
   }
@@ -10680,6 +11551,13 @@ 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
+    {
+      // 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 );
@@ -10701,7 +11579,7 @@ bool _ViscousBuilder::prepareEdgeToShrink( _LayerEdge&            edge,
     edge._simplices[0]._nPrev = n2;
 
     // set U of source node to the target node
-    SMDS_EdgePosition* pos = static_cast<SMDS_EdgePosition*>( tgtNode->GetPosition() );
+    SMDS_EdgePositionPtr pos = tgtNode->GetPosition();
     pos->SetUParameter( uSrc );
   }
   return true;
@@ -10734,7 +11612,7 @@ void _ViscousBuilder::restoreNoShrink( _LayerEdge& edge ) const
       TopLoc_Location loc;
       Handle(Geom_Curve) curve = BRep_Tool::Curve( TopoDS::Edge( S ), loc, f, l );
       if ( curve.IsNull() ) return;
-      SMDS_EdgePosition* ePos = static_cast<SMDS_EdgePosition*>( srcNode->GetPosition() );
+      SMDS_EdgePositionPtr ePos = srcNode->GetPosition();
       p = curve->Value( ePos->GetUParameter() );
       break;
     }
@@ -10752,7 +11630,7 @@ void _ViscousBuilder::restoreNoShrink( _LayerEdge& edge ) const
 
 //================================================================================
 /*!
- * \brief Try to fix triangles with high aspect ratio by swaping diagonals
+ * \brief Try to fix triangles with high aspect ratio by swapping diagonals
  */
 //================================================================================
 
@@ -10960,7 +11838,7 @@ bool _LayerEdge::SetNewLength2d( Handle(Geom_Surface)& surface,
     {
       return true;
     }
-    SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( tgtNode->GetPosition() );
+    SMDS_FacePositionPtr pos = tgtNode->GetPosition();
     pos->SetUParameter( newUV.X() );
     pos->SetVParameter( newUV.Y() );
 
@@ -10974,7 +11852,7 @@ bool _LayerEdge::SetNewLength2d( Handle(Geom_Surface)& surface,
   {
     const TopoDS_Edge&      E = TopoDS::Edge( eos._sWOL );
     const SMDS_MeshNode*   n2 = _simplices[0]._nPrev;
-    SMDS_EdgePosition* tgtPos = static_cast<SMDS_EdgePosition*>( tgtNode->GetPosition() );
+    SMDS_EdgePositionPtr tgtPos = tgtNode->GetPosition();
 
     const double u2     = helper.GetNodeU( E, n2, tgtNode );
     const double uSrc   = _pos[0].Coord( U_SRC );
@@ -11094,7 +11972,7 @@ bool _SmoothNode::Smooth(int&                  nbBad,
     return false;
   }
 
-  SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( _node->GetPosition() );
+  SMDS_FacePositionPtr pos = _node->GetPosition();
   pos->SetUParameter( newPos.X() );
   pos->SetVParameter( newPos.Y() );
 
@@ -11114,7 +11992,7 @@ bool _SmoothNode::Smooth(int&                  nbBad,
 
 //================================================================================
 /*!
- * \brief Computes new UV using angle based smoothing technic
+ * \brief Computes new UV using angle based smoothing technique
  */
 //================================================================================
 
@@ -11172,34 +12050,6 @@ gp_XY _SmoothNode::computeAngularPos(vector<gp_XY>& uv,
   return newPos;
 }
 
-//================================================================================
-/*!
- * \brief Delete _SolidData
- */
-//================================================================================
-
-_SolidData::~_SolidData()
-{
-  TNode2Edge::iterator n2e = _n2eMap.begin();
-  for ( ; n2e != _n2eMap.end(); ++n2e )
-  {
-    _LayerEdge* & e = n2e->second;
-    if ( e )
-    {
-      delete e->_curvature;
-      if ( e->_2neibors )
-        delete e->_2neibors->_plnNorm;
-      delete e->_2neibors;
-    }
-    delete e;
-    e = 0;
-  }
-  _n2eMap.clear();
-
-  delete _helper;
-  _helper = 0;
-}
-
 //================================================================================
 /*!
  * \brief Keep a _LayerEdge inflated along the EDGE
@@ -11321,7 +12171,7 @@ void _Shrinker1D::Compute(bool set3D, SMESH_MesherHelper& helper)
       if ( !discret.IsDone() )
         return throw SALOME_Exception(LOCALIZED("GCPnts_AbscissaPoint failed"));
       double u = discret.Parameter();
-      SMDS_EdgePosition* pos = static_cast<SMDS_EdgePosition*>( _nodes[i]->GetPosition() );
+      SMDS_EdgePositionPtr pos = _nodes[i]->GetPosition();
       pos->SetUParameter( u );
       gp_Pnt p = C->Value( u );
       const_cast< SMDS_MeshNode*>( _nodes[i] )->setXYZ( p.X(), p.Y(), p.Z() );
@@ -11339,7 +12189,7 @@ void _Shrinker1D::Compute(bool set3D, SMESH_MesherHelper& helper)
     {
       if ( !_nodes[i] ) continue;
       double u = f * ( 1-_normPar[i] ) + l * _normPar[i];
-      SMDS_EdgePosition* pos = static_cast<SMDS_EdgePosition*>( _nodes[i]->GetPosition() );
+      SMDS_EdgePositionPtr pos = _nodes[i]->GetPosition();
       pos->SetUParameter( u );
     }
   }
@@ -11357,7 +12207,7 @@ void _Shrinker1D::RestoreParams()
     for ( size_t i = 0; i < _nodes.size(); ++i )
     {
       if ( !_nodes[i] ) continue;
-      SMDS_EdgePosition* pos = static_cast<SMDS_EdgePosition*>( _nodes[i]->GetPosition() );
+      SMDS_EdgePositionPtr pos = _nodes[i]->GetPosition();
       pos->SetUParameter( _initU[i] );
     }
   _done = false;
@@ -11365,7 +12215,7 @@ void _Shrinker1D::RestoreParams()
 
 //================================================================================
 /*!
- * \brief Replace source nodes by target nodes in shrinked mesh edges
+ * \brief Replace source nodes by target nodes in shrunk mesh edges
  */
 //================================================================================
 
@@ -11444,11 +12294,11 @@ bool _ViscousBuilder::addBoundaryElements(_SolidData& data)
         const SMDS_MeshNode* tgtN0 = ledges[0]->_nodes.back();
         const SMDS_MeshNode* tgtN1 = ledges[1]->_nodes.back();
         int nbSharedPyram = 0;
-        SMDS_ElemIteratorPtr vIt = tgtN0->GetInverseElementIterator(SMDSAbs_Volume);
+        SMDS_ElemIteratorPtr vIt = tgtN1->GetInverseElementIterator(SMDSAbs_Volume);
         while ( vIt->more() )
         {
           const SMDS_MeshElement* v = vIt->next();
-          nbSharedPyram += int( v->GetNodeIndex( tgtN1 ) >= 0 );
+          nbSharedPyram += int( v->GetNodeIndex( tgtN0 ) >= 0 );
         }
         if ( nbSharedPyram > 1 )
           continue; // not free border of the pyramid