Salome HOME
Copyright update 2020
[modules/smesh.git] / src / StdMeshers / StdMeshers_ViscousLayers.cxx
index a384191b2be9903701f64001acb577d3e5007b28..310b164d5add0f235547b13ddf6c5fed30941035 100644 (file)
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2016  CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2020  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
@@ -23,6 +23,7 @@
 
 #include "StdMeshers_ViscousLayers.hxx"
 
+#include "ObjectPool.hxx"
 #include "SMDS_EdgePosition.hxx"
 #include "SMDS_FaceOfNodes.hxx"
 #include "SMDS_FacePosition.hxx"
@@ -33,6 +34,7 @@
 #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>
@@ -371,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; }
   };
@@ -582,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;
     }
@@ -622,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; }
@@ -636,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;
   };
 
   //--------------------------------------------------------------------------------
@@ -683,6 +687,7 @@ namespace VISCOUS_3D
     _SolidData&      GetData() const { return *_data; }
 
     _EdgesOnShape(): _shapeID(-1), _subMesh(0), _toSmooth(false), _edgeSmoother(0) {}
+    ~_EdgesOnShape();
   };
 
   //--------------------------------------------------------------------------------
@@ -742,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
@@ -758,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
@@ -785,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 );
@@ -887,6 +893,7 @@ namespace VISCOUS_3D
                             const double   refSign );
   };
   struct PyDump;
+  struct Periodicity;
   //--------------------------------------------------------------------------------
   /*!
    * \brief Builder of viscous layers
@@ -910,13 +917,15 @@ namespace VISCOUS_3D
 
   private:
 
-    bool findSolidsWithLayers();
+    bool findSolidsWithLayers(const bool checkFaceMesh=true);
     bool setBefore( _SolidData& solidBefore, _SolidData& solidAfter );
     bool findFacesWithLayers(const bool onlyWith=false);
+    void findPeriodicFaces();
     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,
@@ -999,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;
   };
   //--------------------------------------------------------------------------------
   /*!
@@ -1211,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
 
 
@@ -1221,7 +1252,8 @@ namespace VISCOUS_3D
 StdMeshers_ViscousLayers::StdMeshers_ViscousLayers(int hypId, SMESH_Gen* gen)
   :SMESH_Hypothesis(hypId, gen),
    _isToIgnoreShapes(1), _nbLayers(1), _thickness(1), _stretchFactor(1),
-   _method( SURF_OFFSET_SMOOTH )
+   _method( SURF_OFFSET_SMOOTH ),
+   _groupName("")
 {
   _name = StdMeshers_ViscousLayers::GetHypType();
   _param_algo_dim = -3; // auxiliary hyp used by 3D algos
@@ -1253,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,
@@ -1307,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)
@@ -1319,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
@@ -1352,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
 //================================================================================
 
@@ -1590,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
@@ -1690,7 +1771,7 @@ namespace VISCOUS_3D
     PyDump(SMESH_Mesh& m) {
       int tag = 3 + m.GetId();
       const char* fname = "/tmp/viscous.py";
-      cout << "execfile('"<<fname<<"')"<<endl;
+      cout << "exec(open('"<<fname<<"','rb').read() )"<<endl;
       py = _pyStream = new ofstream(fname);
       *py << "import SMESH" << endl
           << "from salome.smesh import smeshBuilder" << endl
@@ -1872,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() )
@@ -1883,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
@@ -1919,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 );
@@ -1945,7 +2040,7 @@ SMESH_ComputeErrorPtr _ViscousBuilder::CheckHypotheses( SMESH_Mesh&         mesh
     return SMESH_ComputeErrorPtr(); // everything already computed
 
 
-  findSolidsWithLayers();
+  findSolidsWithLayers( /*checkFaceMesh=*/false );
   bool ok = findFacesWithLayers( true );
 
   // remove _MeshOfSolid's of _SolidData's
@@ -1964,7 +2059,7 @@ SMESH_ComputeErrorPtr _ViscousBuilder::CheckHypotheses( SMESH_Mesh&         mesh
  */
 //================================================================================
 
-bool _ViscousBuilder::findSolidsWithLayers()
+bool _ViscousBuilder::findSolidsWithLayers(const bool checkFaceMesh)
 {
   // get all solids
   TopTools_IndexedMapOfShape allSolids;
@@ -1974,13 +2069,28 @@ bool _ViscousBuilder::findSolidsWithLayers()
   SMESH_HypoFilter filter;
   for ( int i = 1; i <= allSolids.Extent(); ++i )
   {
-    // find StdMeshers_ViscousLayers hyp assigned to the i-th solid
     SMESH_subMesh* sm = _mesh->GetSubMesh( allSolids(i) );
     if ( sm->GetSubMeshDS() && sm->GetSubMeshDS()->NbElements() > 0 )
       continue; // solid is already meshed
+    // TODO: check if algo is hidden
     SMESH_Algo* algo = sm->GetAlgo();
     if ( !algo ) continue;
-    // TODO: check if algo is hidden
+    // check if all FACEs are meshed, which can be false if Compute() a sub-shape
+    if ( checkFaceMesh )
+    {
+      bool facesMeshed = true;
+      SMESH_subMeshIteratorPtr smIt = sm->getDependsOnIterator(false,true);
+      while ( smIt->more() && facesMeshed )
+      {
+        SMESH_subMesh * faceSM = smIt->next();
+        if ( faceSM->GetSubShape().ShapeType() != TopAbs_FACE )
+          break;
+        facesMeshed = faceSM->IsMeshComputed();
+      }
+      if ( !facesMeshed )
+        continue;
+    }
+    // find StdMeshers_ViscousLayers hyp assigned to the i-th solid
     const list <const SMESHDS_Hypothesis *> & allHyps =
       algo->GetUsedHypothesis(*_mesh, allSolids(i), /*ignoreAuxiliary=*/false);
     _SolidData* soData = 0;
@@ -2145,7 +2255,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 ))
@@ -2181,7 +2291,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())
@@ -2261,7 +2371,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 )
@@ -2486,17 +2596,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();
@@ -2520,6 +2619,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;
 
@@ -2530,34 +2631,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() )
@@ -2595,7 +2681,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 );
@@ -3233,12 +3319,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
@@ -3248,6 +3344,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
@@ -3382,6 +3511,16 @@ bool _EdgesOnShape::GetNormal( const SMDS_MeshElement* face, gp_Vec& norm )
   return ok;
 }
 
+//================================================================================
+/*!
+ * \brief EdgesOnShape destructor
+ */
+//================================================================================
+
+_EdgesOnShape::~_EdgesOnShape()
+{
+  delete _edgeSmoother;
+}
 
 //================================================================================
 /*!
@@ -3396,12 +3535,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
@@ -3674,7 +3814,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
   }
 
@@ -4090,7 +4230,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
  */
 //================================================================================
 
@@ -4135,7 +4275,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
  */
 //================================================================================
 
@@ -4160,7 +4328,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()
@@ -4204,8 +4371,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 )
@@ -6307,7 +6477,7 @@ gp_XYZ _Smoother1D::getNormalNormal( const gp_XYZ & normal,
 void _Smoother1D::offPointsToPython() const
 {
   const char* fname = "/tmp/offPoints.py";
-  cout << "execfile('"<<fname<<"')"<<endl;
+  cout << "exec(open('"<<fname<<"','rb').read() )"<<endl;
   ofstream py(fname);
   py << "import SMESH" << endl
      << "from salome.smesh import smeshBuilder" << endl
@@ -6926,7 +7096,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;
       }
@@ -8268,7 +8438,7 @@ void _LayerEdge::MoveNearConcaVer( const _EdgesOnShape*    eov,
     if ( !edgeF->_curvature )
       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;
@@ -8952,7 +9122,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
  */
 //================================================================================
 
@@ -10140,6 +10310,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();
@@ -10170,14 +10345,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 )
         {
@@ -10190,16 +10371,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
@@ -10208,10 +10391,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 )
         {
@@ -10230,9 +10416,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 );
           }
@@ -10240,15 +10426,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 );
           }
           }
@@ -10259,6 +10445,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
 
@@ -10278,6 +10469,451 @@ bool _ViscousBuilder::refine(_SolidData& data)
   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
@@ -10294,11 +10930,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 ))
@@ -10357,15 +10993,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() )
@@ -10411,11 +11056,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 );
+          }
       }
     }
 
@@ -10487,13 +11133,16 @@ bool _ViscousBuilder::shrink(_SolidData& theData)
         if ( eos.SWOLType() == TopAbs_EDGE )
         {
           SMESH_subMesh* edgeSM = _mesh->GetSubMesh( eos._sWOL );
-          _Shrinker1D& shrinker  = e2shrMap[ edgeSM->GetId() ];
-          eShri1D.insert( & shrinker );
-          shrinker.AddEdge( eos._edges[0], eos, helper );
           VISCOUS_3D::ToClearSubWithMain( edgeSM, data._solid );
-          // restore params of nodes on EDGE if the EDGE has been already
-          // shrinked while shrinking other FACE
-          shrinker.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 )
         {
@@ -10508,7 +11157,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 ) {
@@ -10527,24 +11176,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();
@@ -10635,7 +11284,7 @@ bool _ViscousBuilder::shrink(_SolidData& theData)
       //   }
       // }
 
-    } // while ( shrinked )
+    } // while ( shrunk )
 
     if ( !errMsg.empty() ) // Try to re-compute the shrink FACE
     {
@@ -10673,6 +11322,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 )
@@ -10799,7 +11450,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;
@@ -10843,7 +11494,7 @@ bool _ViscousBuilder::shrink(_SolidData& theData)
   } // 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 )
@@ -10915,6 +11566,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 );
@@ -11407,34 +12065,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
@@ -11600,7 +12230,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
  */
 //================================================================================