Salome HOME
Fix regressions
[modules/smesh.git] / src / StdMeshers / StdMeshers_ViscousLayers.cxx
index 1a117aff222f58264a7e1d2ff30eb78281070213..a6530b9244a5782abe3a820e82075e2aab951077 100644 (file)
@@ -95,7 +95,7 @@
 #include <string>
 
 #ifdef _DEBUG_
-#define __myDEBUG
+//#define __myDEBUG
 //#define __NOT_INVALIDATE_BAD_SMOOTH
 //#define __NODES_AT_POS
 #endif
@@ -409,7 +409,7 @@ namespace VISCOUS_3D
 
     gp_XYZ              _normal;    // to boundary of solid
     vector<gp_XYZ>      _pos;       // points computed during inflation
-    double              _len;       // length achived with the last inflation step
+    double              _len;       // length achieved with the last inflation step
     double              _maxLen;    // maximal possible length
     double              _cosin;     // of angle (_normal ^ surface)
     double              _minAngle;  // of _simplices
@@ -439,11 +439,13 @@ namespace VISCOUS_3D
                   SMOOTHED_C1     = 0x0000800, // is on _eosC1
                   DISTORTED       = 0x0001000, // was bad before smoothing
                   RISKY_SWOL      = 0x0002000, // SWOL is parallel to a source FACE
-                  UNUSED_FLAG     = 0x0100000
+                  SHRUNK          = 0x0004000, // target node reached a tgt position while shrink()
+                  UNUSED_FLAG     = 0x0100000  // to add use flags after
     };
     bool Is   ( int flag ) const { return _flags & flag; }
     void Set  ( int flag ) { _flags |= flag; }
     void Unset( int flag ) { _flags &= ~flag; }
+    std::string DumpFlags() const; // debug
 
     void SetNewLength( double len, _EdgesOnShape& eos, SMESH_MesherHelper& helper );
     bool SetNewLength2d( Handle(Geom_Surface)& surface,
@@ -721,6 +723,7 @@ namespace VISCOUS_3D
   {
     typedef const StdMeshers_ViscousLayers* THyp;
     TopoDS_Shape                    _solid;
+    TopTools_MapOfShape             _before; // SOLIDs to be computed before _solid
     TGeomID                         _index; // SOLID id
     _MeshOfSolid*                   _proxyMesh;
     list< THyp >                    _hyps;
@@ -891,6 +894,7 @@ namespace VISCOUS_3D
   private:
 
     bool findSolidsWithLayers();
+    bool setBefore( _SolidData& solidBefore, _SolidData& solidAfter );
     bool findFacesWithLayers(const bool onlyWith=false);
     void getIgnoreFaces(const TopoDS_Shape&             solid,
                         const StdMeshers_ViscousLayers* hyp,
@@ -957,7 +961,7 @@ namespace VISCOUS_3D
                         _LayerEdge&   edge,
                         const gp_XYZ& newNormal);
     bool refine(_SolidData& data);
-    bool shrink();
+    bool shrink(_SolidData& data);
     bool prepareEdgeToShrink( _LayerEdge& edge, _EdgesOnShape& eos,
                               SMESH_MesherHelper& helper,
                               const SMESHDS_SubMesh* faceSubMesh );
@@ -967,7 +971,7 @@ namespace VISCOUS_3D
                      const bool                  is2D,
                      const int                   step,
                      set<const SMDS_MeshNode*> * involvedNodes=NULL);
-    bool addBoundaryElements();
+    bool addBoundaryElements(_SolidData& data);
 
     bool error( const string& text, int solidID=-1 );
     SMESHDS_Mesh* getMeshDS() const { return _mesh->GetMeshDS(); }
@@ -975,11 +979,14 @@ namespace VISCOUS_3D
     // debug
     void makeGroupOfLE();
 
-    SMESH_Mesh*           _mesh;
-    SMESH_ComputeErrorPtr _error;
+    SMESH_Mesh*                _mesh;
+    SMESH_ComputeErrorPtr      _error;
 
-    vector< _SolidData >  _sdVec;
-    int                   _tmpFaceID;
+    vector<                    _SolidData >  _sdVec;
+    TopTools_IndexedMapOfShape _solids; // to find _SolidData by a solid
+    TopTools_MapOfShape        _shrinkedFaces;
+
+    int                        _tmpFaceID;
   };
   //--------------------------------------------------------------------------------
   /*!
@@ -1232,8 +1239,8 @@ StdMeshers_ViscousLayers::Compute(SMESH_Mesh&         theMesh,
                                   const bool          toMakeN2NMap) const
 {
   using namespace VISCOUS_3D;
-  _ViscousBuilder bulder;
-  SMESH_ComputeErrorPtr err = bulder.Compute( theMesh, theShape );
+  _ViscousBuilder builder;
+  SMESH_ComputeErrorPtr err = builder.Compute( theMesh, theShape );
   if ( err && !err->IsOK() )
     return SMESH_ProxyMesh::Ptr();
 
@@ -1245,7 +1252,7 @@ StdMeshers_ViscousLayers::Compute(SMESH_Mesh&         theMesh,
          _ViscousListener::GetSolidMesh( &theMesh, exp.Current(), /*toCreate=*/false))
     {
       if ( toMakeN2NMap && !pm->_n2nMapComputed )
-        if ( !bulder.MakeN2NMap( pm ))
+        if ( !builder.MakeN2NMap( pm ))
           return SMESH_ProxyMesh::Ptr();
       components.push_back( SMESH_ProxyMesh::Ptr( pm ));
       pm->myIsDeletable = false; // it will de deleted by boost::shared_ptr
@@ -1309,8 +1316,8 @@ StdMeshers_ViscousLayers::CheckHypothesis(SMESH_Mesh&                          t
                                           const TopoDS_Shape&                  theShape,
                                           SMESH_Hypothesis::Hypothesis_Status& theStatus)
 {
-  VISCOUS_3D::_ViscousBuilder bulder;
-  SMESH_ComputeErrorPtr err = bulder.CheckHypotheses( theMesh, theShape );
+  VISCOUS_3D::_ViscousBuilder builder;
+  SMESH_ComputeErrorPtr err = builder.CheckHypotheses( theMesh, theShape );
   if ( err && !err->IsOK() )
     theStatus = SMESH_Hypothesis::HYP_INCOMPAT_HYPS;
   else
@@ -1834,8 +1841,6 @@ bool _ViscousBuilder::MakeN2NMap( _MeshOfSolid* pm )
 SMESH_ComputeErrorPtr _ViscousBuilder::Compute(SMESH_Mesh&         theMesh,
                                                const TopoDS_Shape& theShape)
 {
-  // TODO: set priority of solids during Gen::Compute()
-
   _mesh = & theMesh;
 
   // check if proxy mesh already computed
@@ -1857,22 +1862,37 @@ SMESH_ComputeErrorPtr _ViscousBuilder::Compute(SMESH_Mesh&         theMesh,
 
   for ( size_t i = 0; i < _sdVec.size(); ++i )
   {
-    if ( ! makeLayer(_sdVec[i]) )
+    size_t iSD = 0;
+    for ( iSD = 0; iSD < _sdVec.size(); ++iSD ) // find next SOLID to compute
+      if ( _sdVec[iSD]._before.IsEmpty() &&
+           !_sdVec[iSD]._solid.IsNull() &&
+           _sdVec[iSD]._n2eMap.empty() )
+        break;
+
+    if ( ! makeLayer(_sdVec[iSD]) )   // create _LayerEdge's
       return _error;
 
-    if ( _sdVec[i]._n2eMap.size() == 0 )
+    if ( _sdVec[iSD]._n2eMap.size() == 0 ) // no layers in a SOLID
+    {
+      _sdVec[iSD]._solid.Nullify();
       continue;
-    
-    if ( ! inflate(_sdVec[i]) )
+    }
+
+    if ( ! inflate(_sdVec[iSD]) )     // increase length of _LayerEdge's
       return _error;
 
-    if ( ! refine(_sdVec[i]) )
+    if ( ! refine(_sdVec[iSD]) )      // create nodes and prisms
       return _error;
-  }
-  if ( !shrink() )
-    return _error;
 
-  addBoundaryElements();
+    if ( ! shrink(_sdVec[iSD]) )      // shrink 2D mesh on FACEs w/o layer
+      return _error;
+
+    addBoundaryElements(_sdVec[iSD]); // create quadrangles on prism bare sides
+
+    const TopoDS_Shape& solid = _sdVec[iSD]._solid;
+    for ( iSD = 0; iSD < _sdVec.size(); ++iSD )
+      _sdVec[iSD]._before.Remove( solid );
+  }
 
   makeGroupOfLE(); // debug
   debugDump.Finish();
@@ -1921,12 +1941,14 @@ bool _ViscousBuilder::findSolidsWithLayers()
   TopExp::MapShapes( _mesh->GetShapeToMesh(), TopAbs_SOLID, allSolids );
   _sdVec.reserve( allSolids.Extent());
 
-  SMESH_Gen* gen = _mesh->GetGen();
   SMESH_HypoFilter filter;
   for ( int i = 1; i <= allSolids.Extent(); ++i )
   {
     // find StdMeshers_ViscousLayers hyp assigned to the i-th solid
-    SMESH_Algo* algo = gen->GetAlgo( *_mesh, allSolids(i) );
+    SMESH_subMesh* sm = _mesh->GetSubMesh( allSolids(i) );
+    if ( sm->GetSubMeshDS() && sm->GetSubMeshDS()->NbElements() > 0 )
+      continue; // solid is already meshed
+    SMESH_Algo* algo = sm->GetAlgo();
     if ( !algo ) continue;
     // TODO: check if algo is hidden
     const list <const SMESHDS_Hypothesis *> & allHyps =
@@ -1951,6 +1973,7 @@ bool _ViscousBuilder::findSolidsWithLayers()
           soData->_index = getMeshDS()->ShapeToIndex( allSolids(i));
           soData->_helper = new SMESH_MesherHelper( *_mesh );
           soData->_helper->SetSubShape( allSolids(i) );
+          _solids.Add( allSolids(i) );
         }
         soData->_hyps.push_back( viscHyp );
         soData->_hypShapes.push_back( hypShape );
@@ -1965,7 +1988,37 @@ bool _ViscousBuilder::findSolidsWithLayers()
 
 //================================================================================
 /*!
- * \brief 
+ * \brief Set a _SolidData to be computed before another
+ */
+//================================================================================
+
+bool _ViscousBuilder::setBefore( _SolidData& solidBefore, _SolidData& solidAfter )
+{
+  // check possibility to set this order; get all solids before solidBefore
+  TopTools_IndexedMapOfShape allSolidsBefore;
+  allSolidsBefore.Add( solidBefore._solid );
+  for ( int i = 1; i <= allSolidsBefore.Extent(); ++i )
+  {
+    int iSD = _solids.FindIndex( allSolidsBefore(i) );
+    if ( iSD )
+    {
+      TopTools_MapIteratorOfMapOfShape soIt( _sdVec[ iSD-1 ]._before );
+      for ( ; soIt.More(); soIt.Next() )
+        allSolidsBefore.Add( soIt.Value() );
+    }
+  }
+  if ( allSolidsBefore.Contains( solidAfter._solid ))
+    return false;
+
+  for ( int i = 1; i <= allSolidsBefore.Extent(); ++i )
+    solidAfter._before.Add( allSolidsBefore(i) );
+
+  return true;
+}
+
+//================================================================================
+/*!
+ * \brief
  */
 //================================================================================
 
@@ -1973,13 +2026,10 @@ bool _ViscousBuilder::findFacesWithLayers(const bool onlyWith)
 {
   SMESH_MesherHelper helper( *_mesh );
   TopExp_Explorer exp;
-  TopTools_IndexedMapOfShape solids;
 
   // collect all faces-to-ignore defined by hyp
   for ( size_t i = 0; i < _sdVec.size(); ++i )
   {
-    solids.Add( _sdVec[i]._solid );
-
     // get faces-to-ignore defined by each hyp
     typedef const StdMeshers_ViscousLayers* THyp;
     typedef std::pair< set<TGeomID>, THyp > TFacesOfHyp;
@@ -2075,8 +2125,9 @@ bool _ViscousBuilder::findFacesWithLayers(const bool onlyWith)
     }
   }
 
-  // Find faces to shrink mesh on (solution 2 in issue 0020832);
+  // Find FACEs to shrink mesh on (solution 2 in issue 0020832): fill in _shrinkShape2Shape
   TopTools_IndexedMapOfShape shapes;
+  std::string structAlgoName = "Hexa_3D";
   for ( size_t i = 0; i < _sdVec.size(); ++i )
   {
     shapes.Clear();
@@ -2084,121 +2135,32 @@ bool _ViscousBuilder::findFacesWithLayers(const bool onlyWith)
     for ( int iE = 1; iE <= shapes.Extent(); ++iE )
     {
       const TopoDS_Shape& edge = shapes(iE);
-      // find 2 faces sharing an edge
+      // find 2 FACEs sharing an EDGE
       TopoDS_Shape FF[2];
-      PShapeIteratorPtr fIt = helper.GetAncestors(edge, *_mesh, TopAbs_FACE);
+      PShapeIteratorPtr fIt = helper.GetAncestors(edge, *_mesh, TopAbs_FACE, &_sdVec[i]._solid);
       while ( fIt->more())
       {
         const TopoDS_Shape* f = fIt->next();
-        if ( helper.IsSubShape( *f, _sdVec[i]._solid))
-          FF[ int( !FF[0].IsNull()) ] = *f;
+        FF[ int( !FF[0].IsNull()) ] = *f;
       }
       if( FF[1].IsNull() ) continue; // seam edge can be shared by 1 FACE only
+
       // check presence of layers on them
       int ignore[2];
       for ( int j = 0; j < 2; ++j )
-        ignore[j] = _sdVec[i]._ignoreFaceIds.count ( getMeshDS()->ShapeToIndex( FF[j] ));
+        ignore[j] = _sdVec[i]._ignoreFaceIds.count( getMeshDS()->ShapeToIndex( FF[j] ));
       if ( ignore[0] == ignore[1] )
         continue; // nothing interesting
       TopoDS_Shape fWOL = FF[ ignore[0] ? 0 : 1 ];
-      // check presence of layers on fWOL within an adjacent SOLID
-      bool collision = false;
-      PShapeIteratorPtr sIt = helper.GetAncestors( fWOL, *_mesh, TopAbs_SOLID );
-      while ( const TopoDS_Shape* solid = sIt->next() )
-        if ( !solid->IsSame( _sdVec[i]._solid ))
-        {
-          int iSolid = solids.FindIndex( *solid );
-          int  iFace = getMeshDS()->ShapeToIndex( fWOL );
-          if ( iSolid > 0 && !_sdVec[ iSolid-1 ]._ignoreFaceIds.count( iFace ))
-          {
-            //_sdVec[i]._noShrinkShapes.insert( iFace );
-            //fWOL.Nullify();
-            collision = true;
-          }
-        }
-      // add edge to maps
+
+      // add EDGE to maps
       if ( !fWOL.IsNull())
       {
         TGeomID edgeInd = getMeshDS()->ShapeToIndex( edge );
         _sdVec[i]._shrinkShape2Shape.insert( make_pair( edgeInd, fWOL ));
-        if ( collision )
-        {
-          // _shrinkShape2Shape will be used to temporary inflate _LayerEdge's based
-          // on the edge but shrink won't be performed
-          _sdVec[i]._noShrinkShapes.insert( edgeInd );
-        }
       }
     }
   }
-  // Exclude from _shrinkShape2Shape FACE's that can't be shrinked since
-  // the algo of the SOLID sharing the FACE does not support it
-  set< string > notSupportAlgos; notSupportAlgos.insert("Hexa_3D");
-  for ( size_t i = 0; i < _sdVec.size(); ++i )
-  {
-    map< TGeomID, TopoDS_Shape >::iterator e2f = _sdVec[i]._shrinkShape2Shape.begin();
-    for ( ; e2f != _sdVec[i]._shrinkShape2Shape.end(); ++e2f )
-    {
-      const TopoDS_Shape& fWOL = e2f->second;
-      const TGeomID     edgeID = e2f->first;
-      bool notShrinkFace = false;
-      PShapeIteratorPtr soIt = helper.GetAncestors(fWOL, *_mesh, TopAbs_SOLID);
-      while ( soIt->more() )
-      {
-        const TopoDS_Shape* solid = soIt->next();
-        if ( _sdVec[i]._solid.IsSame( *solid )) continue;
-        SMESH_Algo* algo = _mesh->GetGen()->GetAlgo( *_mesh, *solid );
-        if ( !algo || !notSupportAlgos.count( algo->GetName() )) continue;
-        notShrinkFace = true;
-        size_t iSolid = 0;
-        for ( ; iSolid < _sdVec.size(); ++iSolid )
-        {
-          if ( _sdVec[iSolid]._solid.IsSame( *solid ) ) {
-            if ( _sdVec[iSolid]._shrinkShape2Shape.count( edgeID ))
-              notShrinkFace = false;
-            break;
-          }
-        }
-        if ( notShrinkFace )
-        {
-          _sdVec[i]._noShrinkShapes.insert( edgeID );
-
-          // add VERTEXes of the edge in _noShrinkShapes
-          TopoDS_Shape edge = getMeshDS()->IndexToShape( edgeID );
-          for ( TopoDS_Iterator vIt( edge ); vIt.More(); vIt.Next() )
-            _sdVec[i]._noShrinkShapes.insert( getMeshDS()->ShapeToIndex( vIt.Value() ));
-
-          // check if there is a collision with to-shrink-from EDGEs in iSolid
-          if ( iSolid == _sdVec.size() )
-            continue; // no VL in the solid
-          shapes.Clear();
-          TopExp::MapShapes( fWOL, TopAbs_EDGE, shapes);
-          for ( int iE = 1; iE <= shapes.Extent(); ++iE )
-          {
-            const TopoDS_Edge& E = TopoDS::Edge( shapes( iE ));
-            const TGeomID    eID = getMeshDS()->ShapeToIndex( E );
-            if ( eID == edgeID ||
-                 !_sdVec[iSolid]._shrinkShape2Shape.count( eID ) ||
-                 _sdVec[i]._noShrinkShapes.count( eID ))
-              continue;
-            for ( int is1st = 0; is1st < 2; ++is1st )
-            {
-              TopoDS_Vertex V = helper.IthVertex( is1st, E );
-              if ( _sdVec[i]._noShrinkShapes.count( getMeshDS()->ShapeToIndex( V ) ))
-              {
-                // _sdVec[i]._noShrinkShapes.insert( eID );
-                // V = helper.IthVertex( !is1st, E );
-                // _sdVec[i]._noShrinkShapes.insert( getMeshDS()->ShapeToIndex( V ));
-                //iE = 0; // re-start the loop on EDGEs of fWOL
-                return error("No way to make a conformal mesh with "
-                             "the given set of faces with layers", _sdVec[i]._index);
-              }
-            }
-          }
-        }
-
-      } // while ( soIt->more() )
-    } // loop on _sdVec[i]._shrinkShape2Shape
-  } // loop on _sdVec to fill in _SolidData::_noShrinkShapes
 
   // Find the SHAPE along which to inflate _LayerEdge based on VERTEX
 
@@ -2212,18 +2174,14 @@ bool _ViscousBuilder::findFacesWithLayers(const bool onlyWith)
       // find faces WOL sharing the vertex
       vector< TopoDS_Shape > facesWOL;
       size_t totalNbFaces = 0;
-      PShapeIteratorPtr fIt = helper.GetAncestors(vertex, *_mesh, TopAbs_FACE);
+      PShapeIteratorPtr fIt = helper.GetAncestors(vertex, *_mesh, TopAbs_FACE, &_sdVec[i]._solid );
       while ( fIt->more())
       {
         const TopoDS_Shape* f = fIt->next();
-        if ( helper.IsSubShape( *f, _sdVec[i]._solid ) )
-        {
-          totalNbFaces++;
-          const int fID = getMeshDS()->ShapeToIndex( *f );
-          if ( _sdVec[i]._ignoreFaceIds.count ( fID ) /*&&
-               !_sdVec[i]._noShrinkShapes.count( fID )*/)
-            facesWOL.push_back( *f );
-        }
+        totalNbFaces++;
+        const int fID = getMeshDS()->ShapeToIndex( *f );
+        if ( _sdVec[i]._ignoreFaceIds.count ( fID ) /*&& !_sdVec[i]._noShrinkShapes.count( fID )*/)
+          facesWOL.push_back( *f );
       }
       if ( facesWOL.size() == totalNbFaces || facesWOL.empty() )
         continue; // no layers at this vertex or no WOL
@@ -2273,7 +2231,130 @@ bool _ViscousBuilder::findFacesWithLayers(const bool onlyWith)
     }
   }
 
-  // add FACEs of other SOLIDs to _ignoreFaceIds
+  // Add to _noShrinkShapes sub-shapes of FACE's that can't be shrinked since
+  // the algo of the SOLID sharing the FACE does not support it or for other reasons
+  set< string > notSupportAlgos; notSupportAlgos.insert( structAlgoName );
+  for ( size_t i = 0; i < _sdVec.size(); ++i )
+  {
+    map< TGeomID, TopoDS_Shape >::iterator e2f = _sdVec[i]._shrinkShape2Shape.begin();
+    for ( ; e2f != _sdVec[i]._shrinkShape2Shape.end(); ++e2f )
+    {
+      const TopoDS_Shape& fWOL = e2f->second;
+      const TGeomID     edgeID = e2f->first;
+      TGeomID           faceID = getMeshDS()->ShapeToIndex( fWOL );
+      TopoDS_Shape        edge = getMeshDS()->IndexToShape( edgeID );
+      if ( edge.ShapeType() != TopAbs_EDGE )
+        continue; // shrink shape is VERTEX
+
+      TopoDS_Shape solid;
+      PShapeIteratorPtr soIt = helper.GetAncestors(fWOL, *_mesh, TopAbs_SOLID);
+      while ( soIt->more() && solid.IsNull() )
+      {
+        const TopoDS_Shape* so = soIt->next();
+        if ( !so->IsSame( _sdVec[i]._solid ))
+          solid = *so;
+      }
+      if ( solid.IsNull() )
+        continue;
+
+      bool noShrinkE = false;
+      SMESH_Algo*  algo = _mesh->GetSubMesh( solid )->GetAlgo();
+      bool isStructured = ( algo && algo->GetName() == structAlgoName );
+      size_t     iSolid = _solids.FindIndex( solid ) - 1;
+      if ( iSolid < _sdVec.size() && _sdVec[ iSolid ]._ignoreFaceIds.count( faceID ))
+      {
+        // the adjacent SOLID has NO layers on fWOL;
+        // shrink allowed if
+        // - there are layers on the EDGE in the adjacent SOLID
+        // - there are NO layers in the adjacent SOLID && algo is unstructured and computed later
+        bool hasWLAdj = (_sdVec[iSolid]._shrinkShape2Shape.count( edgeID ));
+        bool shrinkAllowed = (( hasWLAdj ) ||
+                              ( !isStructured && setBefore( _sdVec[ i ], _sdVec[ iSolid ] )));
+        noShrinkE = !shrinkAllowed;
+      }
+      else if ( iSolid < _sdVec.size() )
+      {
+        // the adjacent SOLID has layers on fWOL;
+        // check if SOLID's mesh is unstructured and then try to set it
+        // to be computed after the i-th solid
+        if ( isStructured || !setBefore( _sdVec[ i ], _sdVec[ iSolid ] ))
+          noShrinkE = true; // don't shrink fWOL
+      }
+      else
+      {
+        // the adjacent SOLID has NO layers at all
+        noShrinkE = isStructured;
+      }
+
+      if ( noShrinkE )
+      {
+        _sdVec[i]._noShrinkShapes.insert( edgeID );
+
+        // check if there is a collision with to-shrink-from EDGEs in iSolid
+        // if ( iSolid < _sdVec.size() )
+        // {
+        //   shapes.Clear();
+        //   TopExp::MapShapes( fWOL, TopAbs_EDGE, shapes);
+        //   for ( int iE = 1; iE <= shapes.Extent(); ++iE )
+        //   {
+        //     const TopoDS_Edge& E = TopoDS::Edge( shapes( iE ));
+        //     const TGeomID    eID = getMeshDS()->ShapeToIndex( E );
+        //     if ( eID == edgeID ||
+        //          !_sdVec[iSolid]._shrinkShape2Shape.count( eID ) ||
+        //          _sdVec[i]._noShrinkShapes.count( eID ))
+        //       continue;
+        //     for ( int is1st = 0; is1st < 2; ++is1st )
+        //     {
+        //       TopoDS_Vertex V = helper.IthVertex( is1st, E );
+        //       if ( _sdVec[i]._noShrinkShapes.count( getMeshDS()->ShapeToIndex( V ) ))
+        //       {
+        //         return error("No way to make a conformal mesh with "
+        //                      "the given set of faces with layers", _sdVec[i]._index);
+        //       }
+        //     }
+        //   }
+        // }
+      }
+
+      // add VERTEXes of the edge in _noShrinkShapes, which is necessary if
+      // _shrinkShape2Shape is different in the adjacent SOLID
+      for ( TopoDS_Iterator vIt( edge ); vIt.More(); vIt.Next() )
+      {
+        TGeomID vID = getMeshDS()->ShapeToIndex( vIt.Value() );
+        bool noShrinkV = false;
+
+        if ( iSolid < _sdVec.size() )
+        {
+          if ( _sdVec[ iSolid ]._ignoreFaceIds.count( faceID ))
+          {
+            map< TGeomID, TopoDS_Shape >::iterator i2S, i2SAdj;
+            i2S    = _sdVec[i     ]._shrinkShape2Shape.find( vID );
+            i2SAdj = _sdVec[iSolid]._shrinkShape2Shape.find( vID );
+            if ( i2SAdj == _sdVec[iSolid]._shrinkShape2Shape.end() )
+              noShrinkV = ( i2S->second.ShapeType() == TopAbs_EDGE || isStructured );
+            else
+              noShrinkV = ( ! i2S->second.IsSame( i2SAdj->second ));
+          }
+          else
+          {
+            noShrinkV = noShrinkE;
+          }
+        }
+        else
+        {
+          // the adjacent SOLID has NO layers at all
+          noShrinkV = ( isStructured ||
+                        _sdVec[i]._shrinkShape2Shape[ vID ].ShapeType() == TopAbs_EDGE );
+        }
+        if ( noShrinkV )
+          _sdVec[i]._noShrinkShapes.insert( vID );
+      }
+
+    } // loop on _sdVec[i]._shrinkShape2Shape
+  } // loop on _sdVec to fill in _SolidData::_noShrinkShapes
+
+
+    // add FACEs of other SOLIDs to _ignoreFaceIds
   for ( size_t i = 0; i < _sdVec.size(); ++i )
   {
     shapes.Clear();
@@ -2872,7 +2953,7 @@ bool _ViscousBuilder::findShapesToSmooth( _SolidData& data )
   //const double tgtThick = /*Min( 0.5 * data._geomSize, */data._maxThickness;
 
   // Find shapes needing smoothing; such a shape has _LayerEdge._normal on it's
-  // boundry inclined to the shape at a sharp angle
+  // boundary inclined to the shape at a sharp angle
 
   //list< TGeomID > shapesToSmooth;
   TopTools_MapOfShape edgesOfSmooFaces;
@@ -4081,7 +4162,7 @@ void _LayerEdge::SetDataByNeighbors( const SMDS_MeshNode* n1,
 //================================================================================
 /*!
  * \brief Copy data from a _LayerEdge of other SOLID and based on the same node;
- * this and other _LayerEdge's are inflated along a FACE or an EDGE
+ * this and the other _LayerEdge are inflated along a FACE or an EDGE
  */
 //================================================================================
 
@@ -4431,7 +4512,10 @@ bool _ViscousBuilder::inflate(_SolidData& data)
       const double shapeTgtThick = eos._hyp.GetTotalThickness();
       for ( size_t i = 0; i < eos._edges.size(); ++i )
       {
-        avgThick      += Min( 1., eos._edges[i]->_len / shapeTgtThick );
+        if ( eos._edges[i]->_nodes.size() > 1 )
+          avgThick    += Min( 1., eos._edges[i]->_len / shapeTgtThick );
+        else
+          avgThick    += shapeTgtThick;
         nbActiveEdges += ( ! eos._edges[i]->Is( _LayerEdge::BLOCKED ));
       }
     }
@@ -6605,11 +6689,9 @@ bool _ViscousBuilder::updateNormals( _SolidData&         data,
         else // edge inflates along a FACE
         {
           TopoDS_Shape V = helper.GetSubShapeByNode( edge->_nodes[0], getMeshDS() );
-          PShapeIteratorPtr eIt = helper.GetAncestors( V, *_mesh, TopAbs_EDGE );
+          PShapeIteratorPtr eIt = helper.GetAncestors( V, *_mesh, TopAbs_EDGE, &eos->_sWOL );
           while ( const TopoDS_Shape* E = eIt->next() )
           {
-            if ( !helper.IsSubShape( *E, /*FACE=*/eos->_sWOL ))
-              continue;
             gp_Vec edgeDir = getEdgeDir( TopoDS::Edge( *E ), TopoDS::Vertex( V ));
             double   angle = edgeDir.Angle( newEdge._normal ); // [0,PI]
             if ( angle < M_PI / 2 )
@@ -8433,7 +8515,7 @@ gp_XYZ _LayerEdge::smoothNefPolygon()
 {
   gp_XYZ newPos(0,0,0);
 
-  // get a plane to seach a solution on
+  // get a plane to search a solution on
 
   vector< gp_XYZ > vecs( _simplices.size() + 1 );
   size_t i;
@@ -8624,7 +8706,7 @@ gp_XYZ _LayerEdge::smoothNefPolygon()
 { ////////////////////////////////// NEW
   gp_XYZ newPos(0,0,0);
 
-  // get a plane to seach a solution on
+  // get a plane to search a solution on
 
   size_t i;
   gp_XYZ center(0,0,0);
@@ -8971,7 +9053,7 @@ void _LayerEdge::Block( _SolidData& data )
       minDist   = Min( pSrc.SquareDistance( pTgtN ), minDist );
       minDist   = Min( pTgt.SquareDistance( pSrcN ), minDist );
       double newMaxLen = edge->_maxLen + 0.5 * Sqrt( minDist );
-      if ( edge->_nodes[0]->getshapeId() == neibor->_nodes[0]->getshapeId() )
+      //if ( edge->_nodes[0]->getshapeId() == neibor->_nodes[0]->getshapeId() ) viscous_layers_00/A3
       {
         newMaxLen *= edge->_lenFactor / neibor->_lenFactor;
       }
@@ -9125,10 +9207,48 @@ void _LayerEdge::SmoothPos( const vector< double >& segLen, const double tol )
 
 //================================================================================
 /*!
- * \brief Create layers of prisms
+ * \brief Print flags
  */
 //================================================================================
 
+std::string _LayerEdge::DumpFlags() const
+{
+  SMESH_Comment dump;
+  for ( int flag = 1; flag < 0x1000000; flag *= 2 )
+    if ( _flags & flag )
+    {
+      EFlags f = (EFlags) flag;
+      switch ( f ) {
+      case TO_SMOOTH:       dump << "TO_SMOOTH";       break;
+      case MOVED:           dump << "MOVED";           break;
+      case SMOOTHED:        dump << "SMOOTHED";        break;
+      case DIFFICULT:       dump << "DIFFICULT";       break;
+      case ON_CONCAVE_FACE: dump << "ON_CONCAVE_FACE"; break;
+      case BLOCKED:         dump << "BLOCKED";         break;
+      case INTERSECTED:     dump << "INTERSECTED";     break;
+      case NORMAL_UPDATED:  dump << "NORMAL_UPDATED";  break;
+      case MARKED:          dump << "MARKED";          break;
+      case MULTI_NORMAL:    dump << "MULTI_NORMAL";    break;
+      case NEAR_BOUNDARY:   dump << "NEAR_BOUNDARY";   break;
+      case SMOOTHED_C1:     dump << "SMOOTHED_C1";     break;
+      case DISTORTED:       dump << "DISTORTED";       break;
+      case RISKY_SWOL:      dump << "RISKY_SWOL";      break;
+      case SHRUNK:          dump << "SHRUNK";          break;
+      case UNUSED_FLAG:     dump << "UNUSED_FLAG";     break;
+      }
+      dump << " ";
+    }
+  cout << dump << endl;
+  return dump;
+}
+
+//================================================================================
+/*!
+  case brief:
+  default:
+*/
+//================================================================================
+
 bool _ViscousBuilder::refine(_SolidData& data)
 {
   SMESH_MesherHelper& helper = data.GetHelper();
@@ -9475,7 +9595,6 @@ bool _ViscousBuilder::refine(_SolidData& data)
   set< vector<const SMDS_MeshNode*>* >    nnSet;
   set< int >                       degenEdgeInd;
   vector<const SMDS_MeshElement*>     degenVols;
-  vector<int>                       isRiskySWOL;
 
   TopExp_Explorer exp( data._solid, TopAbs_FACE );
   for ( ; exp.More(); exp.Next() )
@@ -9493,7 +9612,6 @@ bool _ViscousBuilder::refine(_SolidData& data)
       nnVec.resize( nbNodes );
       nnSet.clear();
       degenEdgeInd.clear();
-      isRiskySWOL.resize( nbNodes );
       size_t maxZ = 0, minZ = std::numeric_limits<size_t>::max();
       SMDS_NodeIteratorPtr nIt = face->nodeIterator();
       for ( int iN = 0; iN < nbNodes; ++iN )
@@ -9504,7 +9622,6 @@ bool _ViscousBuilder::refine(_SolidData& data)
         nnVec[ i ] = & edge->_nodes;
         maxZ = std::max( maxZ, nnVec[ i ]->size() );
         minZ = std::min( minZ, nnVec[ i ]->size() );
-        //isRiskySWOL[ i ] = edge->Is( _LayerEdge::RISKY_SWOL );
 
         if ( helper.HasDegeneratedEdges() )
           nnSet.insert( nnVec[ i ]);
@@ -9613,7 +9730,7 @@ bool _ViscousBuilder::refine(_SolidData& data)
     if ( !err || err->IsOK() )
     {
       err.reset( new SMESH_ComputeError( COMPERR_WARNING,
-                                         "Degenerated volumes created" ));
+                                         "Bad quality volumes created" ));
       err->myBadElements.insert( err->myBadElements.end(),
                                  degenVols.begin(),degenVols.end() );
     }
@@ -9628,71 +9745,94 @@ bool _ViscousBuilder::refine(_SolidData& data)
  */
 //================================================================================
 
-bool _ViscousBuilder::shrink()
+bool _ViscousBuilder::shrink(_SolidData& theData)
 {
-  // make map of (ids of FACEs to shrink mesh on) to (_SolidData containing _LayerEdge's
-  // inflated along FACE or EDGE)
-  map< TGeomID, _SolidData* > f2sdMap;
+  // make map of (ids of FACEs to shrink mesh on) to (list of _SolidData containing
+  // _LayerEdge's inflated along FACE or EDGE)
+  map< TGeomID, list< _SolidData* > > f2sdMap;
   for ( size_t i = 0 ; i < _sdVec.size(); ++i )
   {
     _SolidData& data = _sdVec[i];
-    TopTools_MapOfShape FFMap;
     map< TGeomID, TopoDS_Shape >::iterator s2s = data._shrinkShape2Shape.begin();
     for (; s2s != data._shrinkShape2Shape.end(); ++s2s )
-      if ( s2s->second.ShapeType() == TopAbs_FACE )
+      if ( s2s->second.ShapeType() == TopAbs_FACE && !_shrinkedFaces.Contains( s2s->second ))
       {
-        f2sdMap.insert( make_pair( getMeshDS()->ShapeToIndex( s2s->second ), &data ));
+        f2sdMap[ getMeshDS()->ShapeToIndex( s2s->second )].push_back( &data );
 
-        if ( FFMap.Add( (*s2s).second ))
-          // Put mesh faces on the shrinked FACE to the proxy sub-mesh to avoid
-          // usage of mesh faces made in addBoundaryElements() by the 3D algo or
-          // by StdMeshers_QuadToTriaAdaptor
-          if ( SMESHDS_SubMesh* smDS = getMeshDS()->MeshElements( s2s->second ))
+        // Put mesh faces on the shrinked FACE to the proxy sub-mesh to avoid
+        // usage of mesh faces made in addBoundaryElements() by the 3D algo or
+        // by StdMeshers_QuadToTriaAdaptor
+        if ( SMESHDS_SubMesh* smDS = getMeshDS()->MeshElements( s2s->second ))
+        {
+          SMESH_ProxyMesh::SubMesh* proxySub =
+            data._proxyMesh->getFaceSubM( TopoDS::Face( s2s->second ), /*create=*/true);
+          if ( proxySub->NbElements() == 0 )
           {
-            SMESH_ProxyMesh::SubMesh* proxySub =
-              data._proxyMesh->getFaceSubM( TopoDS::Face( s2s->second ), /*create=*/true);
             SMDS_ElemIteratorPtr fIt = smDS->GetElements();
             while ( fIt->more() )
-              proxySub->AddElement( fIt->next() );
-            // as a result 3D algo will use elements from proxySub and not from smDS
+            {
+              const SMDS_MeshElement* f = fIt->next();
+              // as a result 3D algo will use elements from proxySub and not from smDS
+              proxySub->AddElement( f );
+              f->setIsMarked( true );
+
+              // Mark nodes on the FACE to discriminate them from nodes
+              // added by addBoundaryElements(); marked nodes are to be smoothed while shrink()
+              for ( int iN = 0, nbN = f->NbNodes(); iN < nbN; ++iN )
+              {
+                const SMDS_MeshNode* n = f->GetNode( iN );
+                if ( n->GetPosition()->GetDim() == 2 )
+                  n->setIsMarked( true );
+              }
+            }
           }
+        }
       }
   }
 
   SMESH_MesherHelper helper( *_mesh );
   helper.ToFixNodeParameters( true );
 
-  // EDGE's to shrink
+  // EDGEs to shrink
   map< TGeomID, _Shrinker1D > e2shrMap;
   vector< _EdgesOnShape* > subEOS;
   vector< _LayerEdge* > lEdges;
 
-  // loop on FACES to srink mesh on
-  map< TGeomID, _SolidData* >::iterator f2sd = f2sdMap.begin();
+  // loop on FACEs to srink mesh on
+  map< TGeomID, list< _SolidData* > >::iterator f2sd = f2sdMap.begin();
   for ( ; f2sd != f2sdMap.end(); ++f2sd )
   {
-    _SolidData&      data = *f2sd->second;
+    list< _SolidData* > & dataList = f2sd->second;
+    if ( dataList.front()->_n2eMap.empty() ||
+         dataList.back() ->_n2eMap.empty() )
+      continue; // not yet computed
+    if ( dataList.front() != &theData &&
+         dataList.back()  != &theData )
+      continue;
+
+    _SolidData&      data = *dataList.front();
+    _SolidData*     data2 = dataList.size() > 1 ? dataList.back() : 0;
     const TopoDS_Face&  F = TopoDS::Face( getMeshDS()->IndexToShape( f2sd->first ));
     SMESH_subMesh*     sm = _mesh->GetSubMesh( F );
     SMESHDS_SubMesh* smDS = sm->GetSubMeshDS();
 
-    Handle(Geom_Surface) surface = BRep_Tool::Surface(F);
+    Handle(Geom_Surface) surface = BRep_Tool::Surface( F );
 
-    helper.SetSubShape(F);
+    _shrinkedFaces.Add( F );
+    helper.SetSubShape( F );
 
     // ===========================
     // Prepare data for shrinking
     // ===========================
 
-    // Collect nodes to smooth, as src nodes are not yet replaced by tgt ones
-    // and hence all nodes on a FACE connected to 2d elements are to be smoothed
+    // Collect nodes to smooth (they are marked at the beginning of this method)
     vector < const SMDS_MeshNode* > smoothNodes;
     {
       SMDS_NodeIteratorPtr nIt = smDS->GetNodes();
       while ( nIt->more() )
       {
         const SMDS_MeshNode* n = nIt->next();
-        if ( n->NbInverseElements( SMDSAbs_Face ) > 0 )
+        if ( n->isMarked() )
           smoothNodes.push_back( n );
       }
     }
@@ -9704,10 +9844,10 @@ bool _ViscousBuilder::shrink()
     {
       vector<_Simplex> simplices;
       _Simplex::GetSimplices( smoothNodes[0], simplices, ignoreShapes );
-      helper.GetNodeUV( F, simplices[0]._nPrev, 0, &isOkUV ); // fix UV of silpmex nodes
+      helper.GetNodeUV( F, simplices[0]._nPrev, 0, &isOkUV ); // fix UV of simplex nodes
       helper.GetNodeUV( F, simplices[0]._nNext, 0, &isOkUV );
       gp_XY uv = helper.GetNodeUV( F, smoothNodes[0], 0, &isOkUV );
-      if ( !simplices[0].IsForward(uv, smoothNodes[0], F, helper,refSign) )
+      if ( !simplices[0].IsForward(uv, smoothNodes[0], F, helper, refSign ))
         refSign = -1;
     }
 
@@ -9723,8 +9863,13 @@ bool _ViscousBuilder::shrink()
         if ( data._noShrinkShapes.count( subID ))
           continue;
         _EdgesOnShape* eos = data.GetShapeEdges( subID );
-        if ( !eos || eos->_sWOL.IsNull() ) continue;
-
+        if ( !eos || eos->_sWOL.IsNull() )
+          if ( data2 ) // check in adjacent SOLID
+          {
+            eos = data2->GetShapeEdges( subID );
+            if ( !eos || eos->_sWOL.IsNull() )
+              continue;
+          }
         subEOS.push_back( eos );
 
         for ( size_t i = 0; i < eos->_edges.size(); ++i )
@@ -9757,7 +9902,7 @@ bool _ViscousBuilder::shrink()
         while ( fIt->more() )
         {
           const SMDS_MeshElement* f = fIt->next();
-          if ( !smDS->Contains( f ))
+          if ( !smDS->Contains( f ) || !f->isMarked() )
             continue;
           SMDS_NodeIteratorPtr nIt = f->nodeIterator();
           for ( int iN = 0; nIt->more(); ++iN )
@@ -9784,7 +9929,7 @@ bool _ViscousBuilder::shrink()
       {
         const SMDS_MeshNode* n = smoothNodes[i];
         nodesToSmooth[ i ]._node = n;
-        // src nodes must be replaced by tgt nodes to have tgt nodes in _simplices
+        // src nodes must be already replaced by tgt nodes to have tgt nodes in _simplices
         _Simplex::GetSimplices( n, nodesToSmooth[ i ]._simplices, ignoreShapes, 0, sortSimplices);
         // fix up incorrect uv of nodes on the FACE
         helper.GetNodeUV( F, n, 0, &isOkUV);
@@ -9815,6 +9960,10 @@ bool _ViscousBuilder::shrink()
         {
           _LayerEdge& edge = * eos._edges[i];
           _Simplex::GetSimplices( /*tgtNode=*/edge._nodes.back(), edge._simplices, ignoreShapes );
+
+          // additionally mark tgt node; only marked nodes will be used in SetNewLength2d()
+          // not-marked nodes are those added by refine()
+          edge._nodes.back()->setIsMarked( true );
         }
       }
     }
@@ -9951,6 +10100,8 @@ bool _ViscousBuilder::shrink()
 
     if ( !errMsg.empty() ) // Try to re-compute the shrink FACE
     {
+      debugMsg( "Re-compute FACE " << f2sd->first << " because " << errMsg );
+
       // remove faces
       SMESHDS_SubMesh* psm = data._proxyMesh->getFaceSubM( F );
       {
@@ -9990,7 +10141,8 @@ bool _ViscousBuilder::shrink()
         {
           _LayerEdge*       edge = subEOS[iS]->_edges[i];
           SMDS_MeshNode* tgtNode = const_cast< SMDS_MeshNode*& >( edge->_nodes.back() );
-          if ( edge->_pos.empty() ) continue;
+          if ( edge->_pos.empty() ||
+               edge->Is( _LayerEdge::SHRUNK )) continue;
           if ( subEOS[iS]->SWOLType() == TopAbs_FACE )
           {
             SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( tgtNode->GetPosition() );
@@ -10146,6 +10298,8 @@ bool _ViscousBuilder::shrink()
 
     // Set an event listener to clear FACE sub-mesh together with SOLID sub-mesh
     VISCOUS_3D::ToClearSubWithMain( sm, data._solid );
+    if ( data2 )
+      VISCOUS_3D::ToClearSubWithMain( sm, data2->_solid );
 
   } // loop on FACES to srink mesh on
 
@@ -10178,6 +10332,7 @@ bool _ViscousBuilder::prepareEdgeToShrink( _LayerEdge&            edge,
     if ( tgtNode->GetPosition()->GetDim() != 2 ) // not inflated edge
     {
       edge._pos.clear();
+      edge.Set( _LayerEdge::SHRUNK );
       return srcNode == tgtNode;
     }
     gp_XY srcUV ( edge._pos[0].X(), edge._pos[0].Y() );          //helper.GetNodeUV( F, srcNode );
@@ -10188,7 +10343,7 @@ bool _ViscousBuilder::prepareEdgeToShrink( _LayerEdge&            edge,
     edge._normal.SetCoord( uvDir.X(),uvDir.Y(), 0 );
     edge._len = uvLen;
 
-    edge._pos.resize(1);
+    //edge._pos.resize(1);
     edge._pos[0].SetCoord( tgtUV.X(), tgtUV.Y(), 0 );
 
     // set UV of source node to target node
@@ -10201,6 +10356,7 @@ bool _ViscousBuilder::prepareEdgeToShrink( _LayerEdge&            edge,
     if ( tgtNode->GetPosition()->GetDim() != 1 ) // not inflated edge
     {
       edge._pos.clear();
+      edge.Set( _LayerEdge::SHRUNK );
       return srcNode == tgtNode;
     }
     const TopoDS_Edge&    E = TopoDS::Edge( eos._sWOL );
@@ -10222,16 +10378,17 @@ bool _ViscousBuilder::prepareEdgeToShrink( _LayerEdge&            edge,
 
     double uSrc = helper.GetNodeU( E, srcNode, n2 );
     double uTgt = helper.GetNodeU( E, tgtNode, srcNode );
-    double u2   = helper.GetNodeU( E, n2, srcNode );
+    double u2   = helper.GetNodeU( E, n2,      srcNode );
 
-    edge._pos.clear();
+    //edge._pos.clear();
 
     if ( fabs( uSrc-uTgt ) < 0.99 * fabs( uSrc-u2 ))
     {
       // tgtNode is located so that it does not make faces with wrong orientation
+      edge.Set( _LayerEdge::SHRUNK );
       return true;
     }
-    edge._pos.resize(1);
+    //edge._pos.resize(1);
     edge._pos[0].SetCoord( U_TGT, uTgt );
     edge._pos[0].SetCoord( U_SRC, uSrc );
     edge._pos[0].SetCoord( LEN_TGT, fabs( uSrc-uTgt ));
@@ -10452,7 +10609,7 @@ bool _LayerEdge::SetNewLength2d( Handle(Geom_Surface)& surface,
                                  _EdgesOnShape&        eos,
                                  SMESH_MesherHelper&   helper )
 {
-  if ( _pos.empty() )
+  if ( Is( SHRUNK ))
     return false; // already at the target position
 
   SMDS_MeshNode* tgtNode = const_cast< SMDS_MeshNode*& >( _nodes.back() );
@@ -10469,6 +10626,10 @@ bool _LayerEdge::SetNewLength2d( Handle(Geom_Surface)& surface,
     double stepSize = 1e100;
     for ( size_t i = 0; i < _simplices.size(); ++i )
     {
+      if ( !_simplices[i]._nPrev->isMarked() ||
+           !_simplices[i]._nNext->isMarked() )
+        continue; // simplex of quadrangle created by addBoundaryElements()
+
       // find intersection of 2 lines: curUV-tgtUV and that connecting simplex nodes
       gp_XY uvN1 = helper.GetNodeUV( F, _simplices[i]._nPrev );
       gp_XY uvN2 = helper.GetNodeUV( F, _simplices[i]._nNext );
@@ -10484,7 +10645,8 @@ bool _LayerEdge::SetNewLength2d( Handle(Geom_Surface)& surface,
     if ( uvLen <= stepSize )
     {
       newUV = tgtUV;
-      _pos.clear();
+      Set( SHRUNK );
+      //_pos.clear();
     }
     else if ( stepSize > 0 )
     {
@@ -10517,7 +10679,8 @@ bool _LayerEdge::SetNewLength2d( Handle(Geom_Surface)& surface,
     double newU = _pos[0].Coord( U_TGT );
     if ( lenTgt < 0.99 * fabs( uSrc-u2 )) // n2 got out of src-tgt range
     {
-      _pos.clear();
+      Set( _LayerEdge::SHRUNK );
+      //_pos.clear();
     }
     else
     {
@@ -10750,7 +10913,7 @@ void _Shrinker1D::AddEdge( const _LayerEdge*   e,
     _done = false;
   }
   // check _LayerEdge
-  if ( e == _edges[0] || e == _edges[1] )
+  if ( e == _edges[0] || e == _edges[1] || e->_nodes.size() < 2 )
     return;
   if ( eos.SWOLType() != TopAbs_EDGE )
     throw SALOME_Exception(LOCALIZED("Wrong _LayerEdge is added"));
@@ -10787,9 +10950,18 @@ void _Shrinker1D::AddEdge( const _LayerEdge*   e,
     while ( nIt->more() )
     {
       const SMDS_MeshNode* node = nIt->next();
+
+      // skip refinement nodes
       if ( node->NbInverseElements(SMDSAbs_Edge) == 0 ||
            node == tgtNode0 || node == tgtNode1 )
-        continue; // refinement nodes
+        continue;
+      bool hasMarkedFace = false;
+      SMDS_ElemIteratorPtr fIt = node->GetInverseElementIterator(SMDSAbs_Face);
+      while ( fIt->more() && !hasMarkedFace )
+        hasMarkedFace = fIt->next()->isMarked();
+      if ( !hasMarkedFace )
+        continue;
+
       _nodes.push_back( node );
       _initU.push_back( helper.GetNodeU( _geomEdge, node ));
       double len = GCPnts_AbscissaPoint::Length(aCurve, f, _initU.back());
@@ -10822,8 +10994,8 @@ void _Shrinker1D::Compute(bool set3D, SMESH_MesherHelper& helper)
   if ( !e ) e = _edges[1];
   if ( !e ) return;
 
-  _done =  (( !_edges[0] || _edges[0]->_pos.empty() ) &&
-            ( !_edges[1] || _edges[1]->_pos.empty() ));
+  _done =  (( !_edges[0] || _edges[0]->Is( _LayerEdge::SHRUNK )) &&
+            ( !_edges[1] || _edges[1]->Is( _LayerEdge::SHRUNK )));
 
   double f,l;
   if ( set3D || _done )
@@ -10904,11 +11076,12 @@ void _Shrinker1D::SwapSrcTgtNodes( SMESHDS_Mesh* mesh )
     if ( !eSubMesh ) return;
     const SMDS_MeshNode* srcNode = _edges[i]->_nodes[0];
     const SMDS_MeshNode* tgtNode = _edges[i]->_nodes.back();
+    const SMDS_MeshNode* scdNode = _edges[i]->_nodes[1];
     SMDS_ElemIteratorPtr eIt = srcNode->GetInverseElementIterator(SMDSAbs_Edge);
     while ( eIt->more() )
     {
       const SMDS_MeshElement* e = eIt->next();
-      if ( !eSubMesh->Contains( e ))
+      if ( !eSubMesh->Contains( e ) || e->GetNodeIndex( scdNode ) >= 0 )
           continue;
       SMDS_ElemIteratorPtr nIt = e->nodesIterator();
       for ( int iN = 0; iN < e->NbNodes(); ++iN )
@@ -10927,21 +11100,22 @@ void _Shrinker1D::SwapSrcTgtNodes( SMESHDS_Mesh* mesh )
  */
 //================================================================================
 
-bool _ViscousBuilder::addBoundaryElements()
+bool _ViscousBuilder::addBoundaryElements(_SolidData& data)
 {
   SMESH_MesherHelper helper( *_mesh );
 
   vector< const SMDS_MeshNode* > faceNodes;
 
-  for ( size_t i = 0; i < _sdVec.size(); ++i )
+  //for ( size_t i = 0; i < _sdVec.size(); ++i )
   {
-    _SolidData& data = _sdVec[i];
+    //_SolidData& data = _sdVec[i];
     TopTools_IndexedMapOfShape geomEdges;
     TopExp::MapShapes( data._solid, TopAbs_EDGE, geomEdges );
     for ( int iE = 1; iE <= geomEdges.Extent(); ++iE )
     {
       const TopoDS_Edge& E = TopoDS::Edge( geomEdges(iE));
-      if ( data._noShrinkShapes.count( getMeshDS()->ShapeToIndex( E )))
+      const TGeomID edgeID = getMeshDS()->ShapeToIndex( E );
+      if ( data._noShrinkShapes.count( edgeID ))
         continue;
 
       // Get _LayerEdge's based on E
@@ -10990,10 +11164,9 @@ bool _ViscousBuilder::addBoundaryElements()
       // Find out orientation and type of face to create
 
       bool reverse = false, isOnFace;
-      
-      map< TGeomID, TopoDS_Shape >::iterator e2f =
-        data._shrinkShape2Shape.find( getMeshDS()->ShapeToIndex( E ));
       TopoDS_Shape F;
+
+      map< TGeomID, TopoDS_Shape >::iterator e2f = data._shrinkShape2Shape.find( edgeID );
       if (( isOnFace = ( e2f != data._shrinkShape2Shape.end() )))
       {
         F = e2f->second.Oriented( TopAbs_FORWARD );
@@ -11003,17 +11176,12 @@ bool _ViscousBuilder::addBoundaryElements()
         if ( helper.IsReversedSubMesh( TopoDS::Face(F) ))
           reverse = !reverse;
       }
-      else
+      else if ( !data._ignoreFaceIds.count( e2f->first ))
       {
         // find FACE with layers sharing E
-        PShapeIteratorPtr fIt = helper.GetAncestors( E, *_mesh, TopAbs_FACE );
-        while ( fIt->more() && F.IsNull() )
-        {
-          const TopoDS_Shape* pF = fIt->next();
-          if ( helper.IsSubShape( *pF, data._solid) &&
-               !data._ignoreFaceIds.count( e2f->first ))
-            F = *pF;
-        }
+        PShapeIteratorPtr fIt = helper.GetAncestors( E, *_mesh, TopAbs_FACE, &data._solid );
+        if ( fIt->more() )
+          F = *( fIt->next() );
       }
       // Find the sub-mesh to add new faces
       SMESHDS_SubMesh* sm = 0;
@@ -11024,18 +11192,44 @@ bool _ViscousBuilder::addBoundaryElements()
       if ( !sm )
         return error("error in addBoundaryElements()", data._index);
 
+      // Find a proxy sub-mesh of the FACE of an adjacent SOLID, which will use the new boundary
+      // faces for 3D meshing (PAL23414)
+      SMESHDS_SubMesh* adjSM = 0;
+      if ( isOnFace )
+      {
+        const TGeomID   faceID = sm->GetID();
+        PShapeIteratorPtr soIt = helper.GetAncestors( F, *_mesh, TopAbs_SOLID );
+        while ( const TopoDS_Shape* solid = soIt->next() )
+          if ( !solid->IsSame( data._solid ))
+          {
+            size_t iData = _solids.FindIndex( *solid ) - 1;
+            if ( iData < _sdVec.size() &&
+                 _sdVec[ iData ]._ignoreFaceIds.count( faceID ) &&
+                 _sdVec[ iData ]._shrinkShape2Shape.count( edgeID ) == 0 )
+            {
+              SMESH_ProxyMesh::SubMesh* proxySub =
+                _sdVec[ iData ]._proxyMesh->getFaceSubM( TopoDS::Face( F ), /*create=*/false);
+              if ( proxySub && proxySub->NbElements() > 0 )
+                adjSM = proxySub;
+            }
+          }
+      }
+
       // Make faces
       const int dj1 = reverse ? 0 : 1;
       const int dj2 = reverse ? 1 : 0;
+      vector< const SMDS_MeshElement*> ff; // new faces row
+      SMESHDS_Mesh* m = getMeshDS();
       for ( size_t j = 1; j < ledges.size(); ++j )
       {
         vector< const SMDS_MeshNode*>&  nn1 = ledges[j-dj1]->_nodes;
         vector< const SMDS_MeshNode*>&  nn2 = ledges[j-dj2]->_nodes;
+        ff.resize( std::max( nn1.size(), nn2.size() ), NULL );
         if ( nn1.size() == nn2.size() )
         {
           if ( isOnFace )
             for ( size_t z = 1; z < nn1.size(); ++z )
-              sm->AddElement( getMeshDS()->AddFace( nn1[z-1], nn2[z-1], nn2[z], nn1[z] ));
+              sm->AddElement( ff[z-1] = m->AddFace( nn1[z-1], nn2[z-1], nn2[z], nn1[z] ));
           else
             for ( size_t z = 1; z < nn1.size(); ++z )
               sm->AddElement( new SMDS_FaceOfNodes( nn1[z-1], nn2[z-1], nn2[z], nn1[z] ));
@@ -11044,7 +11238,7 @@ bool _ViscousBuilder::addBoundaryElements()
         {
           if ( isOnFace )
             for ( size_t z = 1; z < nn2.size(); ++z )
-              sm->AddElement( getMeshDS()->AddFace( nn1[0], nn2[z-1], nn2[z] ));
+              sm->AddElement( ff[z-1] = m->AddFace( nn1[0], nn2[z-1], nn2[z] ));
           else
             for ( size_t z = 1; z < nn2.size(); ++z )
               sm->AddElement( new SMDS_FaceOfNodes( nn1[0], nn2[z-1], nn2[z] ));
@@ -11053,11 +11247,19 @@ bool _ViscousBuilder::addBoundaryElements()
         {
           if ( isOnFace )
             for ( size_t z = 1; z < nn1.size(); ++z )
-              sm->AddElement( getMeshDS()->AddFace( nn1[z-1], nn2[0], nn1[z] ));
+              sm->AddElement( ff[z-1] = m->AddFace( nn1[z-1], nn2[0], nn1[z] ));
           else
             for ( size_t z = 1; z < nn1.size(); ++z )
               sm->AddElement( new SMDS_FaceOfNodes( nn1[z-1], nn2[0], nn2[z] ));
         }
+
+        if ( adjSM ) // add faces to a proxy SM of the adjacent SOLID
+        {
+          for ( size_t z = 0; z < ff.size(); ++z )
+            if ( ff[ z ])
+              adjSM->AddElement( ff[ z ]);
+          ff.clear();
+        }
       }
 
       // Make edges
@@ -11068,7 +11270,7 @@ bool _ViscousBuilder::addBoundaryElements()
         if ( eos && eos->SWOLType() == TopAbs_EDGE )
         {
           vector< const SMDS_MeshNode*>&  nn = edge->_nodes;
-          if ( nn.size() < 2 || nn[1]->GetInverseElementIterator( SMDSAbs_Edge )->more() )
+          if ( nn.size() < 2 || nn[1]->NbInverseElements( SMDSAbs_Edge ) >= 2 )
             continue;
           helper.SetSubShape( eos->_sWOL );
           helper.SetElementsOnShape( true );