Salome HOME
52459: Viscous layers are not normal to the surface.
[modules/smesh.git] / src / StdMeshers / StdMeshers_ViscousLayers.cxx
index 251cfc37778e749674a22070909850e79707d560..a2a44a23829f391426f09303ed72a60a62c9baef 100644 (file)
@@ -95,6 +95,12 @@ namespace VISCOUS_3D
   enum UIndex { U_TGT = 1, U_SRC, LEN_TGT };
 
   const double theMinSmoothCosin = 0.1;
+  const double theSmoothThickToElemSizeRatio = 0.3;
+
+  bool needSmoothing( double cosin, double tgtThick, double elemSize )
+  {
+    return cosin * tgtThick > theSmoothThickToElemSizeRatio * elemSize;
+  }
 
   /*!
    * \brief SMESH_ProxyMesh computed by _ViscousBuilder for a SOLID.
@@ -286,6 +292,7 @@ namespace VISCOUS_3D
         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 );
       }
@@ -294,30 +301,7 @@ namespace VISCOUS_3D
     double lenDelta(double len) const { return _k * ( _r + len ); }
     double lenDeltaByDist(double dist) const { return dist * _h2lenRatio; }
   };
-  struct _LayerEdge;
-  //--------------------------------------------------------------------------------
-  /*!
-   * Structure used to smooth a _LayerEdge (master) based on an EDGE.
-   */
-  struct _2NearEdges
-  {
-    // target nodes of 2 neighbour _LayerEdge's based on the same EDGE
-    const SMDS_MeshNode* _nodes[2];
-    // vectors from source nodes of 2 _LayerEdge's to the source node of master _LayerEdge
-    //gp_XYZ               _vec[2];
-    double               _wgt[2]; // weights of _nodes
-    _LayerEdge*          _edges[2];
-
-     // normal to plane passing through _LayerEdge._normal and tangent of EDGE
-    gp_XYZ*              _plnNorm;
-
-    _2NearEdges() { _nodes[0]=_nodes[1]=0; _plnNorm = 0; }
-    void reverse() {
-      std::swap( _nodes[0], _nodes[1] );
-      std::swap( _wgt  [0], _wgt  [1] );
-      std::swap( _edges[0], _edges[1] );
-    }
-  };
+  struct _2NearEdges;
   //--------------------------------------------------------------------------------
   /*!
    * \brief Edge normal to surface, connecting a node on solid surface (_nodes[0])
@@ -367,6 +351,7 @@ namespace VISCOUS_3D
                        double&              dist,
                        const double&        epsilon) const;
     gp_Ax1 LastSegment(double& segLen) const;
+    gp_XY  LastUV( const TopoDS_Face& F ) const;
     bool   IsOnEdge() const { return _2neibors; }
     gp_XYZ Copy( _LayerEdge& other, SMESH_MesherHelper& helper );
     void   SetCosin( double cosin );
@@ -379,6 +364,31 @@ namespace VISCOUS_3D
       return cmpNodes ? ( e1->_nodes[0]->GetID() < e2->_nodes[0]->GetID()) : ( e1 < e2 );
     }
   };
+  struct _LayerEdge;
+  //--------------------------------------------------------------------------------
+  /*!
+   * Structure used to smooth a _LayerEdge based on an EDGE.
+   */
+  struct _2NearEdges
+  {
+    double               _wgt  [2]; // weights of _nodes
+    _LayerEdge*          _edges[2];
+
+     // normal to plane passing through _LayerEdge._normal and tangent of EDGE
+    gp_XYZ*              _plnNorm;
+
+    _2NearEdges() { _edges[0]=_edges[1]=0; _plnNorm = 0; }
+    const SMDS_MeshNode* tgtNode(bool is2nd) {
+      return _edges[is2nd] ? _edges[is2nd]->_nodes.back() : 0;
+    }
+    const SMDS_MeshNode* srcNode(bool is2nd) {
+      return _edges[is2nd] ? _edges[is2nd]->_nodes[0] : 0;
+    }
+    void reverse() {
+      std::swap( _wgt  [0], _wgt  [1] );
+      std::swap( _edges[0], _edges[1] );
+    }
+  };
   //--------------------------------------------------------------------------------
   /*!
    * \brief Convex FACE whose radius of curvature is less than the thickness of 
@@ -419,12 +429,13 @@ namespace VISCOUS_3D
     TopoDS_Shape                    _hypShape;
     _MeshOfSolid*                   _proxyMesh;
     set<TGeomID>                    _reversedFaceIds;
-    set<TGeomID>                    _ignoreFaceIds;
+    set<TGeomID>                    _ignoreFaceIds; // WOL FACEs and FACEs of other SOLIDS
 
-    double                          _stepSize, _stepSizeCoeff;
+    double                          _stepSize, _stepSizeCoeff, _geomSize;
     const SMDS_MeshNode*            _stepSizeNodes[2];
 
-    TNode2Edge                      _n2eMap;
+    TNode2Edge                      _n2eMap; // nodes and _LayerEdge's based on them
+
     // map to find _n2eMap of another _SolidData by a shrink shape shared by two _SolidData's
     map< TGeomID, TNode2Edge* >     _s2neMap;
     // edges of _n2eMap. We keep same data in two containers because
@@ -440,8 +451,9 @@ namespace VISCOUS_3D
     // Convex FACEs whose radius of curvature is less than the thickness of layers
     map< TGeomID, _ConvexFace >      _convexFaces;
 
-    // FACE's WOL, srink on which is forbiden due to algo on the adjacent SOLID
-    set< TGeomID >                   _noShrinkFaces;
+    // shapes (EDGEs and VERTEXes) srink from which is forbiden due to collisions with
+    // the adjacent SOLID
+    set< TGeomID >                   _noShrinkShapes;
 
     // <EDGE to smooth on> to <it's curve> -- for analytic smooth
     map< TGeomID,Handle(Geom_Curve)> _edge2curve;
@@ -452,7 +464,7 @@ namespace VISCOUS_3D
 
     double                           _epsilon; // precision for SegTriaInter()
 
-    int                              _index; // for debug
+    TGeomID                          _index; // SOLID id, for debug
 
     _SolidData(const TopoDS_Shape&             s=TopoDS_Shape(),
                const StdMeshers_ViscousLayers* h=0,
@@ -486,7 +498,7 @@ namespace VISCOUS_3D
 
     bool GetShapeEdges(const TGeomID shapeID, size_t& edgeEnd, int* iBeg=0, int* iEnd=0 ) const;
 
-    void AddFacesToSmooth( const set< TGeomID >& faceIDs );
+    void AddShapesToSmooth( const set< TGeomID >& shapeIDs );
   };
   //--------------------------------------------------------------------------------
   /*!
@@ -583,12 +595,13 @@ namespace VISCOUS_3D
                        const bool          toSort = false);
     void findSimplexTestEdges( _SolidData&                    data,
                                vector< vector<_LayerEdge*> >& edgesByGeom);
+    void computeGeomSize( _SolidData& data );
     bool sortEdges( _SolidData&                    data,
                     vector< vector<_LayerEdge*> >& edgesByGeom);
     void limitStepSizeByCurvature( _SolidData&  data );
     void limitStepSize( _SolidData&             data,
                         const SMDS_MeshElement* face,
-                        const double            cosin);
+                        const _LayerEdge*       maxCosinEdge );
     void limitStepSize( _SolidData& data, const double minSize);
     bool inflate(_SolidData& data);
     bool smoothAndCheck(_SolidData& data, const int nbSteps, double & distToIntersection);
@@ -607,6 +620,7 @@ namespace VISCOUS_3D
     bool prepareEdgeToShrink( _LayerEdge& edge, const TopoDS_Face& F,
                               SMESH_MesherHelper& helper,
                               const SMESHDS_SubMesh* faceSubMesh );
+    void restoreNoShrink( _LayerEdge& edge ) const;
     void fixBadFaces(const TopoDS_Face&          F,
                      SMESH_MesherHelper&         helper,
                      const bool                  is2D,
@@ -615,7 +629,7 @@ namespace VISCOUS_3D
     bool addBoundaryElements();
 
     bool error( const string& text, int solidID=-1 );
-    SMESHDS_Mesh* getMeshDS() { return _mesh->GetMeshDS(); }
+    SMESHDS_Mesh* getMeshDS() const { return _mesh->GetMeshDS(); }
 
     // debug
     void makeGroupOfLE();
@@ -658,7 +672,8 @@ namespace VISCOUS_3D
     virtual SMDSAbs_ElementType  GetType() const              { return SMDSAbs_Face; }
     virtual vtkIdType GetVtkType() const                      { return -1; }
     virtual SMDSAbs_EntityType   GetEntityType() const        { return SMDSEntity_Last; }
-    virtual SMDSAbs_GeometryType GetGeomType() const          { return SMDSGeom_TRIANGLE; }
+    virtual SMDSAbs_GeometryType GetGeomType() const
+    { return _nn.size() == 3 ? SMDSGeom_TRIANGLE : SMDSGeom_QUADRANGLE; }
     virtual SMDS_ElemIteratorPtr elementsIterator(SMDSAbs_ElementType) const
     { return SMDS_ElemIteratorPtr( new SMDS_NodeVectorElemIterator( _nn.begin(), _nn.end()));}
   };
@@ -716,6 +731,7 @@ namespace VISCOUS_3D
       return _surface->Value( uv.X(), uv.Y() ).XYZ();
     }
   };
+
 } // namespace VISCOUS_3D
 
 
@@ -841,22 +857,32 @@ namespace
     gp_Vec dir;
     double f,l; gp_Pnt p;
     Handle(Geom_Curve) c = BRep_Tool::Curve( E, f, l );
+    if ( c.IsNull() ) return gp_XYZ( 1e100, 1e100, 1e100 );
     double u = helper.GetNodeU( E, atNode );
     c->D1( u, p, dir );
     return dir.XYZ();
   }
   //--------------------------------------------------------------------------------
+  gp_XYZ getFaceDir( const TopoDS_Face& F, const TopoDS_Vertex& fromV,
+                     const SMDS_MeshNode* node, SMESH_MesherHelper& helper, bool& ok,
+                     double* cosin=0);
+  //--------------------------------------------------------------------------------
   gp_XYZ getFaceDir( const TopoDS_Face& F, const TopoDS_Edge& fromE,
                      const SMDS_MeshNode* node, SMESH_MesherHelper& helper, bool& ok)
   {
+    double f,l;
+    Handle(Geom_Curve) c = BRep_Tool::Curve( fromE, f, l );
+    if ( c.IsNull() )
+    {
+      TopoDS_Vertex v = helper.IthVertex( 0, fromE );
+      return getFaceDir( F, v, node, helper, ok );
+    }
     gp_XY uv = helper.GetNodeUV( F, node, 0, &ok );
     Handle(Geom_Surface) surface = BRep_Tool::Surface( F );
     gp_Pnt p; gp_Vec du, dv, norm;
     surface->D1( uv.X(),uv.Y(), p, du,dv );
     norm = du ^ dv;
 
-    double f,l;
-    Handle(Geom_Curve) c = BRep_Tool::Curve( fromE, f, l );
     double u = helper.GetNodeU( fromE, node, 0, &ok );
     c->D1( u, p, du );
     TopAbs_Orientation o = helper.GetSubShapeOri( F.Oriented(TopAbs_FORWARD), fromE);
@@ -880,7 +906,7 @@ namespace
   //--------------------------------------------------------------------------------
   gp_XYZ getFaceDir( const TopoDS_Face& F, const TopoDS_Vertex& fromV,
                      const SMDS_MeshNode* node, SMESH_MesherHelper& helper,
-                     bool& ok, double* cosin=0)
+                     bool& ok, double* cosin)
   {
     TopoDS_Face faceFrw = F;
     faceFrw.Orientation( TopAbs_FORWARD );
@@ -896,12 +922,12 @@ namespace
         if ( SMESH_Algo::isDegenerated( e )) continue;
         TopExp::Vertices( e, VV[0], VV[1], /*CumOri=*/true );
         if ( VV[1].IsSame( fromV )) {
+          nbEdges += edges[ 0 ].IsNull();
           edges[ 0 ] = e;
-          nbEdges++;
         }
         else if ( VV[0].IsSame( fromV )) {
+          nbEdges += edges[ 1 ].IsNull();
           edges[ 1 ] = e;
-          nbEdges++;
         }
       }
     }
@@ -921,7 +947,7 @@ namespace
 
       // get angle between the 2 edges
       gp_Vec faceNormal;
-      double angle = helper.GetAngle( edges[0], edges[1], faceFrw, &faceNormal );
+      double angle = helper.GetAngle( edges[0], edges[1], faceFrw, fromV, &faceNormal );
       if ( Abs( angle ) < 5 * M_PI/180 )
       {
         dir = ( faceNormal.XYZ() ^ edgeDir[0].Reversed()) + ( faceNormal.XYZ() ^ edgeDir[1] );
@@ -939,7 +965,7 @@ namespace
     }
     else if ( nbEdges == 1 )
     {
-      dir = getFaceDir( faceFrw, edges[0], node, helper, ok );
+      dir = getFaceDir( faceFrw, edges[ edges[0].IsNull() ], node, helper, ok );
       if ( cosin ) *cosin = 1.;
     }
     else
@@ -1003,17 +1029,74 @@ namespace
         while ( SMESH_Algo::isDegenerated( wires[iW]->Edge( iE2 )))
           iE2 = ( iE2 + 1 ) % nbEdges;
         double angle = helper.GetAngle( wires[iW]->Edge( iE1 ),
-                                        wires[iW]->Edge( iE2 ), F );
+                                        wires[iW]->Edge( iE2 ), F,
+                                        wires[iW]->FirstVertex( iE2 ));
         if ( angle < -5. * M_PI / 180. )
           return true;
       }
     }
     return false;
   }
+  //================================================================================
+  /*!
+   * \brief Computes mimimal 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
+   *  \return bool - true if faceSize computed
+   */
+  //================================================================================
+
+  bool getDistFromEdge( const SMDS_MeshElement* face,
+                        const SMDS_MeshNode*    nodeOnEdge,
+                        double &                faceSize )
+  {
+    faceSize = Precision::Infinite();
+    bool done = false;
+
+    int nbN  = face->NbCornerNodes();
+    int iOnE = face->GetNodeIndex( nodeOnEdge );
+    int iNext[2] = { SMESH_MesherHelper::WrapIndex( iOnE+1, nbN ),
+                     SMESH_MesherHelper::WrapIndex( iOnE-1, nbN ) };
+    const SMDS_MeshNode* nNext[2] = { face->GetNode( iNext[0] ),
+                                      face->GetNode( iNext[1] ) };
+    gp_XYZ segVec, segEnd = SMESH_TNodeXYZ( nodeOnEdge ); // segment on EDGE
+    double segLen = -1.;
+    // look for two neighbor not in-FACE nodes of face
+    for ( int i = 0; i < 2; ++i )
+    {
+      if ( nNext[i]->GetPosition()->GetDim() != 2 &&
+           nNext[i]->GetID() < nodeOnEdge->GetID() )
+      {
+        // look for an in-FACE node
+        for ( int iN = 0; iN < nbN; ++iN )
+        {
+          if ( iN == iOnE || iN == iNext[i] )
+            continue;
+          SMESH_TNodeXYZ pInFace = face->GetNode( iN );
+          gp_XYZ v = pInFace - segEnd;
+          if ( segLen < 0 )
+          {
+            segVec = SMESH_TNodeXYZ( nNext[i] ) - segEnd;
+            segLen = segVec.Modulus();
+          }
+          double distToSeg = v.Crossed( segVec ).Modulus() / segLen;
+          faceSize = Min( faceSize, distToSeg );
+          done = true;
+        }
+        segLen = -1;
+      }
+    }
+    return done;
+  }
+
   //--------------------------------------------------------------------------------
   // DEBUG. Dump intermediate node positions into a python script
+  // HOWTO use: run python commands written in a console to see
+  //  construction steps of viscous layers
 #ifdef __myDEBUG
   ofstream* py;
+  int       theNbFunc;
   struct PyDump {
     PyDump() {
       const char* fname = "/tmp/viscous.py";
@@ -1024,19 +1107,24 @@ namespace
           << "smesh  = smeshBuilder.New(salome.myStudy)" << endl
           << "meshSO = smesh.GetCurrentStudy().FindObjectID('0:1:2:3')" << endl
           << "mesh   = smesh.Mesh( meshSO.GetObject() )"<<endl;
+      theNbFunc = 0;
     }
     void Finish() {
-      if (py)
-        *py << "mesh.MakeGroup('Viscous Prisms',SMESH.VOLUME,SMESH.FT_ElemGeomType,'=',SMESH.Geom_PENTA)"<<endl;
+      if (py) {
+        *py << "mesh.GroupOnFilter(SMESH.VOLUME,'Viscous Prisms',"
+          "smesh.GetFilter(SMESH.VOLUME,SMESH.FT_ElemGeomType,'=',SMESH.Geom_PENTA))"<<endl;
+        *py << "mesh.GroupOnFilter(SMESH.VOLUME,'Neg Volumes',"
+          "smesh.GetFilter(SMESH.VOLUME,SMESH.FT_Volume3D,'<',0))"<<endl;
+      }
       delete py; py=0;
     }
-    ~PyDump() { Finish(); }
+    ~PyDump() { Finish(); cout << "NB FUNCTIONS: " << theNbFunc << endl; }
   };
 #define dumpFunction(f) { _dumpFunction(f, __LINE__);}
 #define dumpMove(n)     { _dumpMove(n, __LINE__);}
 #define dumpCmd(txt)    { _dumpCmd(txt, __LINE__);}
   void _dumpFunction(const string& fun, int ln)
-  { if (py) *py<< "def "<<fun<<"(): # "<< ln <<endl; cout<<fun<<"()"<<endl;}
+  { if (py) *py<< "def "<<fun<<"(): # "<< ln <<endl; cout<<fun<<"()"<<endl; ++theNbFunc; }
   void _dumpMove(const SMDS_MeshNode* n, int ln)
   { if (py) *py<< "  mesh.MoveNode( "<<n->GetID()<< ", "<< n->X()
                << ", "<<n->Y()<<", "<< n->Z()<< ")\t\t # "<< ln <<endl; }
@@ -1082,19 +1170,36 @@ _ViscousBuilder::_ViscousBuilder()
 
 bool _ViscousBuilder::error(const string& text, int solidId )
 {
+  const string prefix = string("Viscous layers builder: ");
   _error->myName    = COMPERR_ALGO_FAILED;
-  _error->myComment = string("Viscous layers builder: ") + text;
+  _error->myComment = prefix + text;
   if ( _mesh )
   {
     SMESH_subMesh* sm = _mesh->GetSubMeshContaining( solidId );
     if ( !sm && !_sdVec.empty() )
-      sm = _mesh->GetSubMeshContaining( _sdVec[0]._index );
+      sm = _mesh->GetSubMeshContaining( solidId = _sdVec[0]._index );
     if ( sm && sm->GetSubShape().ShapeType() == TopAbs_SOLID )
     {
       SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
       if ( smError && smError->myAlgo )
         _error->myAlgo = smError->myAlgo;
       smError = _error;
+      sm->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE );
+    }
+    // set KO to all solids
+    for ( size_t i = 0; i < _sdVec.size(); ++i )
+    {
+      if ( _sdVec[i]._index == solidId )
+        continue;
+      sm = _mesh->GetSubMesh( _sdVec[i]._solid );
+      if ( !sm->IsEmpty() )
+        continue;
+      SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
+      if ( !smError || smError->IsOK() )
+      {
+        smError = SMESH_ComputeError::New( COMPERR_ALGO_FAILED, prefix + "failed");
+        sm->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE );
+      }
     }
   }
   makeGroupOfLE(); // debug
@@ -1350,6 +1455,7 @@ bool _ViscousBuilder::findFacesWithLayers()
         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 ))
@@ -1358,8 +1464,9 @@ bool _ViscousBuilder::findFacesWithLayers()
           int  iFace = getMeshDS()->ShapeToIndex( fWOL );
           if ( iSolid > 0 && !_sdVec[ iSolid-1 ]._ignoreFaceIds.count( iFace ))
           {
-            _sdVec[i]._noShrinkFaces.insert( iFace );
-            fWOL.Nullify();
+            //_sdVec[i]._noShrinkShapes.insert( iFace );
+            //fWOL.Nullify();
+            collision = true;
           }
         }
       // add edge to maps
@@ -1367,6 +1474,12 @@ bool _ViscousBuilder::findFacesWithLayers()
       {
         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 );
+        }
       }
     }
   }
@@ -1375,54 +1488,70 @@ bool _ViscousBuilder::findFacesWithLayers()
   set< string > notSupportAlgos; notSupportAlgos.insert("Hexa_3D");
   for ( size_t i = 0; i < _sdVec.size(); ++i )
   {
-    TopTools_MapOfShape noShrinkVertices;
     map< TGeomID, TopoDS_Shape >::iterator e2f = _sdVec[i]._shrinkShape2Shape.begin();
     for ( ; e2f != _sdVec[i]._shrinkShape2Shape.end(); ++e2f )
     {
       const TopoDS_Shape& fWOL = e2f->second;
-      TGeomID           edgeID = e2f->first;
+      const TGeomID     edgeID = e2f->first;
       bool notShrinkFace = false;
       PShapeIteratorPtr soIt = helper.GetAncestors(fWOL, *_mesh, TopAbs_SOLID);
-      while ( soIt->more())
+      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;
-        for ( size_t j = 0; j < _sdVec.size(); ++j )
+        size_t iSolid = 0;
+        for ( ; iSolid < _sdVec.size(); ++iSolid )
         {
-          if ( _sdVec[j]._solid.IsSame( *solid ) )
-            if ( _sdVec[j]._shrinkShape2Shape.count( edgeID ))
+          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);
+              }
+            }
+          }
         }
-      }
-      if ( notShrinkFace )
-      {
-        _sdVec[i]._noShrinkFaces.insert( getMeshDS()->ShapeToIndex( fWOL ));
-        for ( TopExp_Explorer vExp( fWOL, TopAbs_VERTEX ); vExp.More(); vExp.Next() )
-          noShrinkVertices.Add( vExp.Current() );
-      }
-    }
-    // erase from _shrinkShape2Shape all srink EDGE's of a SOLID connected
-    // to the found not shrinked fWOL's
-    e2f = _sdVec[i]._shrinkShape2Shape.begin();
-    for ( ; e2f != _sdVec[i]._shrinkShape2Shape.end(); )
-    {
-      TGeomID edgeID = e2f->first;
-      TopoDS_Vertex VV[2];
-      TopExp::Vertices( TopoDS::Edge( getMeshDS()->IndexToShape( edgeID )),VV[0],VV[1]);
-      if ( noShrinkVertices.Contains( VV[0] ) || noShrinkVertices.Contains( VV[1] ))
-      {
-        _sdVec[i]._noShrinkFaces.insert( getMeshDS()->ShapeToIndex( e2f->second ));
-        _sdVec[i]._shrinkShape2Shape.erase( e2f++ );
-      }
-      else
-      {
-        e2f++;
-      }
-    }
-  }
+
+      } // 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
 
@@ -1444,8 +1573,8 @@ bool _ViscousBuilder::findFacesWithLayers()
         {
           totalNbFaces++;
           const int fID = getMeshDS()->ShapeToIndex( *f );
-          if ( _sdVec[i]._ignoreFaceIds.count ( fID ) &&
-               !_sdVec[i]._noShrinkFaces.count( fID ))
+          if ( _sdVec[i]._ignoreFaceIds.count ( fID ) /*&&
+               !_sdVec[i]._noShrinkShapes.count( fID )*/)
             facesWOL.push_back( *f );
         }
       }
@@ -1523,15 +1652,14 @@ bool _ViscousBuilder::makeLayer(_SolidData& data)
 {
   // get all sub-shapes to make layers on
   set<TGeomID> subIds, faceIds;
-  subIds = data._noShrinkFaces;
+  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( getMeshDS()->ShapeToIndex( exp.Current() )))
+      if ( ! data._ignoreFaceIds.count( fSubM->GetId() ))
         faceIds.insert( fSubM->GetId() );
-      SMESH_subMeshIteratorPtr subIt =
-        fSubM->getDependsOnIterator(/*includeSelf=*/true, /*complexShapeFirst=*/false);
+      SMESH_subMeshIteratorPtr subIt = fSubM->getDependsOnIterator(/*includeSelf=*/true);
       while ( subIt->more() )
         subIds.insert( subIt->next()->GetId() );
     }
@@ -1564,7 +1692,7 @@ bool _ViscousBuilder::makeLayer(_SolidData& data)
 
   SMESH_MesherHelper helper( *_mesh );
   helper.SetSubShape( data._solid );
-  helper.SetElementsOnShape(true);
+  helper.SetElementsOnShape( true );
 
   vector< const SMDS_MeshNode*> newNodes; // of a mesh face
   TNode2Edge::iterator n2e2;
@@ -1586,27 +1714,51 @@ bool _ViscousBuilder::makeLayer(_SolidData& data)
     while ( eIt->more() )
     {
       const SMDS_MeshElement* face = eIt->next();
+      double          faceMaxCosin = -1;
+      _LayerEdge*     maxCosinEdge = 0;
+      int             nbDegenNodes = 0;
+
       newNodes.resize( face->NbCornerNodes() );
-      double faceMaxCosin = -1;
-      for ( int i = 0 ; i < face->NbCornerNodes(); ++i )
+      for ( size_t i = 0 ; i < newNodes.size(); ++i )
       {
-        const SMDS_MeshNode* n = face->GetNode(i);
+        const SMDS_MeshNode* n = face->GetNode( i );
+        const int      shapeID = n->getshapeId();
+        const bool onDegenShap = helper.IsDegenShape( shapeID );
+        const bool onDegenEdge = ( onDegenShap && n->GetPosition()->GetDim() == 1 );
+        if ( onDegenShap )
+        {
+          if ( onDegenEdge )
+          {
+            // substitute n on a degenerated EDGE with a node on a corresponding VERTEX
+            const TopoDS_Shape& E = getMeshDS()->IndexToShape( shapeID );
+            TopoDS_Vertex       V = helper.IthVertex( 0, TopoDS::Edge( E ));
+            if ( const SMDS_MeshNode* vN = SMESH_Algo::VertexNode( V, getMeshDS() )) {
+              n = vN;
+              nbDegenNodes++;
+            }
+          }
+          else
+          {
+            nbDegenNodes++;
+          }
+        }
         TNode2Edge::iterator n2e = data._n2eMap.insert( make_pair( n, (_LayerEdge*)0 )).first;
         if ( !(*n2e).second )
         {
           // add a _LayerEdge
           _LayerEdge* edge = new _LayerEdge();
-          n2e->second = edge;
           edge->_nodes.push_back( n );
-          const int shapeID = n->getshapeId();
+          n2e->second = edge;
           edgesByGeom[ shapeID ].push_back( edge );
+          const bool noShrink = data._noShrinkShapes.count( shapeID );
 
           SMESH_TNodeXYZ xyz( n );
 
           // set edge data or find already refined _LayerEdge and get data from it
-          if ( n->GetPosition()->GetTypeOfPosition() != SMDS_TOP_FACE &&
-               ( s2ne = data._s2neMap.find( shapeID )) != data._s2neMap.end() &&
-               ( n2e2 = (*s2ne).second->find( n )) != s2ne->second->end())
+          if (( !noShrink                                                     ) &&
+              ( n->GetPosition()->GetTypeOfPosition() != SMDS_TOP_FACE        ) &&
+              (( s2ne = data._s2neMap.find( shapeID )) != data._s2neMap.end() ) &&
+              (( n2e2 = (*s2ne).second->find( n )) != s2ne->second->end()     ))
           {
             _LayerEdge* foundEdge = (*n2e2).second;
             gp_XYZ        lastPos = edge->Copy( *foundEdge, helper );
@@ -1617,19 +1769,29 @@ bool _ViscousBuilder::makeLayer(_SolidData& data)
           }
           else
           {
-            edge->_nodes.push_back( helper.AddNode( xyz.X(), xyz.Y(), xyz.Z() ));
+            if ( !noShrink )
+            {
+              edge->_nodes.push_back( helper.AddNode( xyz.X(), xyz.Y(), xyz.Z() ));
+            }
             if ( !setEdgeData( *edge, subIds, helper, data ))
               return false;
           }
           dumpMove(edge->_nodes.back());
-          if ( edge->_cosin > 0.01 )
+
+          if ( edge->_cosin > faceMaxCosin )
           {
-            if ( edge->_cosin > faceMaxCosin )
-              faceMaxCosin = edge->_cosin;
+            faceMaxCosin = edge->_cosin;
+            maxCosinEdge = edge;
           }
         }
         newNodes[ i ] = n2e->second->_nodes.back();
+
+        if ( onDegenEdge )
+          data._n2eMap.insert( make_pair( face->GetNode( i ), n2e->second ));
       }
+      if ( newNodes.size() - nbDegenNodes < 2 )
+        continue;
+
       // create a temporary face
       const SMDS_MeshElement* newFace =
         new _TmpMeshFace( newNodes, --_tmpFaceID, face->getshapeId() );
@@ -1637,7 +1799,8 @@ bool _ViscousBuilder::makeLayer(_SolidData& data)
 
       // compute inflation step size by min size of element on a convex surface
       if ( faceMaxCosin > theMinSmoothCosin )
-        limitStepSize( data, face, faceMaxCosin );
+        limitStepSize( data, face, maxCosinEdge );
+
     } // loop on 2D elements on a FACE
   } // loop on FACEs of a SOLID
 
@@ -1652,28 +1815,61 @@ bool _ViscousBuilder::makeLayer(_SolidData& data)
   // limit data._stepSize depending on surface curvature and fill data._convexFaces
   limitStepSizeByCurvature( data ); // !!! it must be before node substitution in _Simplex
 
-  // Set target nodes into _Simplex and _2NearEdges of _LayerEdge's
+  // Set target nodes into _Simplex and _LayerEdge's to _2NearEdges
   TNode2Edge::iterator n2e;
+  const SMDS_MeshNode* nn[2];
   for ( size_t i = 0; i < data._edges.size(); ++i )
   {
-    if ( data._edges[i]->IsOnEdge())
+    _LayerEdge* edge = data._edges[i];
+    if ( edge->IsOnEdge() )
+    {
+      // get neighbor nodes
+      bool hasData = ( edge->_2neibors->_edges[0] );
+      if ( hasData ) // _LayerEdge is a copy of another one
+      {
+        nn[0] = edge->_2neibors->srcNode(0);
+        nn[1] = edge->_2neibors->srcNode(1);
+      }
+      else if ( !findNeiborsOnEdge( edge, nn[0],nn[1], data ))
+      {
+        return false;
+      }
+      // set neighbor _LayerEdge's
       for ( int j = 0; j < 2; ++j )
       {
-        if ( data._edges[i]->_nodes.back()->NbInverseElements(SMDSAbs_Volume) > 0 )
-          break; // _LayerEdge is shared by two _SolidData's
-        const SMDS_MeshNode* & n = data._edges[i]->_2neibors->_nodes[j];
-        if (( n2e = data._n2eMap.find( n )) == data._n2eMap.end() )
+        if (( n2e = data._n2eMap.find( nn[j] )) == data._n2eMap.end() )
           return error("_LayerEdge not found by src node", data._index);
-        n = (*n2e).second->_nodes.back();
-        data._edges[i]->_2neibors->_edges[j] = n2e->second;
+        edge->_2neibors->_edges[j] = n2e->second;
       }
-    //else
-    for ( size_t j = 0; j < data._edges[i]->_simplices.size(); ++j )
+      if ( !hasData )
+        edge->SetDataByNeighbors( nn[0], nn[1], helper);
+    }
+
+    for ( size_t j = 0; j < edge->_simplices.size(); ++j )
     {
-      _Simplex& s = data._edges[i]->_simplices[j];
+      _Simplex& s = edge->_simplices[j];
       s._nNext = data._n2eMap[ s._nNext ]->_nodes.back();
       s._nPrev = data._n2eMap[ s._nPrev ]->_nodes.back();
     }
+
+    // For an _LayerEdge on a degenerated EDGE, copy some data from
+    // a corresponding _LayerEdge on a VERTEX
+    // (issue 52453, pb on a downloaded SampleCase2-Tet-netgen-mephisto.hdf)
+    if ( helper.IsDegenShape( edge->_nodes[0]->getshapeId() ))
+    {
+      // Generally we should not get here
+      const TopoDS_Shape& E = getMeshDS()->IndexToShape( edge->_nodes[0]->getshapeId() );
+      if ( E.ShapeType() != TopAbs_EDGE )
+        continue;
+      TopoDS_Vertex V = helper.IthVertex( 0, TopoDS::Edge( E ));
+      const SMDS_MeshNode* vN = SMESH_Algo::VertexNode( V, getMeshDS() );
+      if (( n2e = data._n2eMap.find( vN )) == data._n2eMap.end() )
+        continue;
+      const _LayerEdge* vEdge = n2e->second;
+      edge->_normal    = vEdge->_normal;
+      edge->_lenFactor = vEdge->_lenFactor;
+      edge->_cosin     = vEdge->_cosin;
+    }
   }
 
   dumpFunctionEnd();
@@ -1688,7 +1884,7 @@ bool _ViscousBuilder::makeLayer(_SolidData& data)
 
 void _ViscousBuilder::limitStepSize( _SolidData&             data,
                                      const SMDS_MeshElement* face,
-                                     const double            cosin)
+                                     const _LayerEdge*       maxCosinEdge )
 {
   int iN = 0;
   double minSize = 10 * data._stepSize;
@@ -1696,20 +1892,20 @@ void _ViscousBuilder::limitStepSize( _SolidData&             data,
   for ( int i = 0; i < nbNodes; ++i )
   {
     const SMDS_MeshNode* nextN = face->GetNode( SMESH_MesherHelper::WrapIndex( i+1, nbNodes ));
-    const SMDS_MeshNode* curN = face->GetNode( i );
+    const SMDS_MeshNode*  curN = face->GetNode( i );
     if ( nextN->GetPosition()->GetTypeOfPosition() == SMDS_TOP_FACE ||
-         curN->GetPosition()->GetTypeOfPosition() == SMDS_TOP_FACE )
+         curN-> GetPosition()->GetTypeOfPosition() == SMDS_TOP_FACE )
     {
-      double dist = SMESH_TNodeXYZ( face->GetNode(i)).Distance( nextN );
+      double dist = SMESH_TNodeXYZ( curN ).Distance( nextN );
       if ( dist < minSize )
         minSize = dist, iN = i;
     }
   }
-  double newStep = 0.8 * minSize / cosin;
+  double newStep = 0.8 * minSize / maxCosinEdge->_lenFactor;
   if ( newStep < data._stepSize )
   {
     data._stepSize = newStep;
-    data._stepSizeCoeff = 0.8 / cosin;
+    data._stepSizeCoeff = 0.8 / maxCosinEdge->_lenFactor;
     data._stepSizeNodes[0] = face->GetNode( iN );
     data._stepSizeNodes[1] = face->GetNode( SMESH_MesherHelper::WrapIndex( iN+1, nbNodes ));
   }
@@ -1755,6 +1951,10 @@ void _ViscousBuilder::limitStepSizeByCurvature( _SolidData& data )
   for ( ; face.More(); face.Next() )
   {
     const TopoDS_Face& F = TopoDS::Face( face.Current() );
+    SMESH_subMesh *   sm = _mesh->GetSubMesh( F );
+    const TGeomID faceID = sm->GetId();
+    if ( data._ignoreFaceIds.count( faceID )) continue;
+
     BRepAdaptor_Surface surface( F, false );
     surfProp.SetSurface( surface );
 
@@ -1762,9 +1962,6 @@ void _ViscousBuilder::limitStepSizeByCurvature( _SolidData& data )
     int iBeg, iEnd;
 
     _ConvexFace cnvFace;
-    SMESH_subMesh *            sm = _mesh->GetSubMesh( F );
-    const TGeomID          faceID = sm->GetId();
-    if ( data._ignoreFaceIds.count( faceID )) continue;
     const double        oriFactor = ( F.Orientation() == TopAbs_REVERSED ? +1. : -1. );
     SMESH_subMeshIteratorPtr smIt = sm->getDependsOnIterator(/*includeSelf=*/true);
     while ( smIt->more() )
@@ -1779,8 +1976,8 @@ void _ViscousBuilder::limitStepSizeByCurvature( _SolidData& data )
         continue;
       // check concavity and curvature and limit data._stepSize
       int nbLEdges = iEnd - iBeg;
-      int step = Max( 1, nbLEdges / nbTestPnt );
-      for ( ; iBeg < iEnd; iBeg += step )
+      int iStep     = Max( 1, nbLEdges / nbTestPnt );
+      for ( ; iBeg < iEnd; iBeg += iStep )
       {
         gp_XY uv = helper.GetNodeUV( F, data._edges[ iBeg ]->_nodes[0] );
         surfProp.SetParameters( uv.X(), uv.Y() );
@@ -1870,8 +2067,12 @@ void _ViscousBuilder::limitStepSizeByCurvature( _SolidData& data )
 bool _ViscousBuilder::sortEdges( _SolidData&                    data,
                                  vector< vector<_LayerEdge*> >& edgesByGeom)
 {
+  // define allowed thickness
+  computeGeomSize( data ); // compute data._geomSize
+  const double tgtThick = Min( 0.5 * data._geomSize, data._hyp->GetTotalThickness() );
+
   // Find shapes needing smoothing; such a shape has _LayerEdge._normal on it's
-  // boundry inclined at a sharp angle to the shape
+  // boundry inclined to the shape at a sharp angle
 
   list< TGeomID > shapesToSmooth;
   
@@ -1888,28 +2089,33 @@ bool _ViscousBuilder::sortEdges( _SolidData&                    data,
     {
     case TopAbs_EDGE: {
 
-      bool isShrinkEdge = !eS[0]->_sWOL.IsNull();
+      if ( SMESH_Algo::isDegenerated( TopoDS::Edge( S )))
+        break;
+      //bool isShrinkEdge = !eS[0]->_sWOL.IsNull();
       for ( TopoDS_Iterator vIt( S ); vIt.More() && !needSmooth; vIt.Next() )
       {
         TGeomID iV = getMeshDS()->ShapeToIndex( vIt.Value() );
         vector<_LayerEdge*>& eV = edgesByGeom[ iV ];
         if ( eV.empty() ) continue;
-        double cosin = eV[0]->_cosin;
-        bool badCosin =
-          ( !eV[0]->_sWOL.IsNull() && ( eV[0]->_sWOL.ShapeType() == TopAbs_EDGE || !isShrinkEdge));
-        if ( badCosin )
+        gp_Vec eDir = getEdgeDir( TopoDS::Edge( S ), TopoDS::Vertex( vIt.Value() ));
+        double angle = eDir.Angle( eV[0]->_normal );
+        double cosin = Cos( angle );
+        if ( cosin > theMinSmoothCosin )
         {
-          gp_Vec dir1, dir2;
-          if ( eV[0]->_sWOL.ShapeType() == TopAbs_EDGE )
-            dir1 = getEdgeDir( TopoDS::Edge( eV[0]->_sWOL ), TopoDS::Vertex( vIt.Value() ));
-          else
-            dir1 = getFaceDir( TopoDS::Face( eV[0]->_sWOL ), TopoDS::Vertex( vIt.Value() ),
-                               eV[0]->_nodes[0], helper, ok);
-          dir2 = getEdgeDir( TopoDS::Edge( S ), TopoDS::Vertex( vIt.Value() ));
-          double angle = dir1.Angle( dir2 );
-          cosin = cos( angle );
+          // compare tgtThick with the length of an end segment
+          SMDS_ElemIteratorPtr eIt = eV[0]->_nodes[0]->GetInverseElementIterator(SMDSAbs_Edge);
+          while ( eIt->more() )
+          {
+            const SMDS_MeshElement* endSeg = eIt->next();
+            if ( endSeg->getshapeId() == iS )
+            {
+              double segLen =
+                SMESH_TNodeXYZ( endSeg->GetNode(0) ).Distance( endSeg->GetNode(1 ));
+              needSmooth = needSmoothing( cosin, tgtThick, segLen );
+              break;
+            }
+          }
         }
-        needSmooth = ( cosin > theMinSmoothCosin );
       }
       break;
     }
@@ -1920,25 +2126,38 @@ bool _ViscousBuilder::sortEdges( _SolidData&                    data,
         TGeomID iE = getMeshDS()->ShapeToIndex( eExp.Current() );
         vector<_LayerEdge*>& eE = edgesByGeom[ iE ];
         if ( eE.empty() ) continue;
-        if ( eE[0]->_sWOL.IsNull() )
+        // TopLoc_Location loc;
+        // Handle(Geom_Surface) surface = BRep_Tool::Surface( TopoDS::Face( S ), loc );
+        // bool isPlane = GeomLib_IsPlanarSurface( surface ).IsPlanar();
+        //if ( eE[0]->_sWOL.IsNull() )
         {
+          double faceSize;
           for ( size_t i = 0; i < eE.size() && !needSmooth; ++i )
-            needSmooth = ( eE[i]->_cosin > theMinSmoothCosin );
-        }
-        else
-        {
-          const TopoDS_Face& F1 = TopoDS::Face( S );
-          const TopoDS_Face& F2 = TopoDS::Face( eE[0]->_sWOL );
-          const TopoDS_Edge& E  = TopoDS::Edge( eExp.Current() );
-          for ( size_t i = 0; i < eE.size() && !needSmooth; ++i )
-          {
-            gp_Vec dir1 = getFaceDir( F1, E, eE[i]->_nodes[0], helper, ok );
-            gp_Vec dir2 = getFaceDir( F2, E, eE[i]->_nodes[0], helper, ok );
-            double angle = dir1.Angle( dir2 );
-            double cosin = cos( angle );
-            needSmooth = ( cosin > theMinSmoothCosin );
-          }
+            if ( eE[i]->_cosin > theMinSmoothCosin )
+            {
+              SMDS_ElemIteratorPtr fIt = eE[i]->_nodes[0]->GetInverseElementIterator(SMDSAbs_Face);
+              while ( fIt->more() && !needSmooth )
+              {
+                const SMDS_MeshElement* face = fIt->next();
+                if ( getDistFromEdge( face, eE[i]->_nodes[0], faceSize ))
+                  needSmooth = needSmoothing( eE[i]->_cosin, tgtThick, faceSize );
+              }
+            }
         }
+        // else
+        // {
+        //   const TopoDS_Face& F1 = TopoDS::Face( S );
+        //   const TopoDS_Face& F2 = TopoDS::Face( eE[0]->_sWOL );
+        //   const TopoDS_Edge& E  = TopoDS::Edge( eExp.Current() );
+        //   for ( size_t i = 0; i < eE.size() && !needSmooth; ++i )
+        //   {
+        //     gp_Vec dir1 = getFaceDir( F1, E, eE[i]->_nodes[0], helper, ok );
+        //     gp_Vec dir2 = getFaceDir( F2, E, eE[i]->_nodes[0], helper, ok );
+        //     double angle = dir1.Angle(  );
+        //     double cosin = cos( angle );
+        //     needSmooth = ( cosin > theMinSmoothCosin );
+        //   }
+        // }
       }
       break;
     }
@@ -1946,6 +2165,7 @@ bool _ViscousBuilder::sortEdges( _SolidData&                    data,
       continue;
     default:;
     }
+
     if ( needSmooth )
     {
       if ( S.ShapeType() == TopAbs_EDGE ) shapesToSmooth.push_front( iS );
@@ -2000,8 +2220,8 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge&         edge,
   const SMDS_MeshNode* node = edge._nodes[0]; // source node
   SMDS_TypeOfPosition posType = node->GetPosition()->GetTypeOfPosition();
 
-  edge._len = 0;
-  edge._2neibors = 0;
+  edge._len       = 0;
+  edge._2neibors  = 0;
   edge._curvature = 0;
 
   // --------------------------
@@ -2012,18 +2232,16 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge&         edge,
   edge._normal.SetCoord(0,0,0);
 
   int totalNbFaces = 0;
-  gp_Pnt p;
-  gp_Vec du, dv, geomNorm;
+  gp_Vec geomNorm;
   bool normOK = true;
 
-  TGeomID shapeInd = node->getshapeId();
+  const TGeomID shapeInd = node->getshapeId();
   map< TGeomID, TopoDS_Shape >::const_iterator s2s = data._shrinkShape2Shape.find( shapeInd );
-  bool onShrinkShape ( s2s != data._shrinkShape2Shape.end() );
-  TopoDS_Shape vertEdge;
+  const bool onShrinkShape ( s2s != data._shrinkShape2Shape.end() );
 
   if ( onShrinkShape ) // one of faces the node is on has no layers
   {
-    vertEdge = getMeshDS()->IndexToShape( s2s->first ); // vertex or edge
+    TopoDS_Shape vertEdge = getMeshDS()->IndexToShape( s2s->first ); // vertex or edge
     if ( s2s->second.ShapeType() == TopAbs_EDGE )
     {
       // inflate from VERTEX along EDGE
@@ -2152,14 +2370,16 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge&         edge,
     if ( edge._sWOL.ShapeType() == TopAbs_EDGE )
     {
       double u = helper.GetNodeU( TopoDS::Edge( edge._sWOL ), node, 0, &normOK );
-      edge._pos.push_back( gp_XYZ( u, 0, 0));
-      getMeshDS()->SetNodeOnEdge( tgtNode, TopoDS::Edge( edge._sWOL ), u );
+      edge._pos.push_back( gp_XYZ( u, 0, 0 ));
+      if ( edge._nodes.size() > 1 )
+        getMeshDS()->SetNodeOnEdge( tgtNode, TopoDS::Edge( edge._sWOL ), u );
     }
     else // TopAbs_FACE
     {
       gp_XY uv = helper.GetNodeUV( TopoDS::Face( edge._sWOL ), node, 0, &normOK );
       edge._pos.push_back( gp_XYZ( uv.X(), uv.Y(), 0));
-      getMeshDS()->SetNodeOnFace( tgtNode, TopoDS::Face( edge._sWOL ), uv.X(), uv.Y() );
+      if ( edge._nodes.size() > 1 )
+        getMeshDS()->SetNodeOnFace( tgtNode, TopoDS::Face( edge._sWOL ), uv.X(), uv.Y() );
     }
   }
   else
@@ -2189,14 +2409,14 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge&         edge,
   {
     edge._2neibors = new _2NearEdges;
     // target node instead of source ones will be set later
-    if ( ! findNeiborsOnEdge( &edge,
-                              edge._2neibors->_nodes[0],
-                              edge._2neibors->_nodes[1],
-                              data))
-      return false;
-    edge.SetDataByNeighbors( edge._2neibors->_nodes[0],
-                             edge._2neibors->_nodes[1],
-                             helper);
+    // if ( ! findNeiborsOnEdge( &edge,
+    //                           edge._2neibors->_nodes[0],
+    //                           edge._2neibors->_nodes[1],
+    //                           data))
+    //   return false;
+    // edge.SetDataByNeighbors( edge._2neibors->_nodes[0],
+    //                          edge._2neibors->_nodes[1],
+    //                          helper);
   }
 
   edge.SetCosin( edge._cosin ); // to update edge._lenFactor
@@ -2271,12 +2491,30 @@ gp_XYZ _ViscousBuilder::getFaceNormal(const SMDS_MeshNode* node,
   isOK = false;
 
   Handle(Geom_Surface) surface = BRep_Tool::Surface( face );
-  if ( GeomLib::NormEstim( surface, uv, 1e-10, normal ) < 3 )
+  int pointKind = GeomLib::NormEstim( surface, uv, 1e-5, normal );
+  enum { REGULAR = 0, QUASYSINGULAR, CONICAL, IMPOSSIBLE };
+
+  if ( pointKind == IMPOSSIBLE &&
+       node->GetPosition()->GetDim() == 2 ) // node inside the FACE
   {
-    normal;
+    // probably NormEstim() failed due to a too high tolerance
+    pointKind = GeomLib::NormEstim( surface, uv, 1e-20, normal );
+    isOK = ( pointKind < IMPOSSIBLE );
+  }
+  if ( pointKind < IMPOSSIBLE )
+  {
+    if ( pointKind != REGULAR &&
+         !shiftInside &&
+         node->GetPosition()->GetDim() < 2 ) // FACE boundary
+    {
+      gp_XYZ normShift = getFaceNormal( node, face, helper, isOK, /*shiftInside=*/true );
+      if ( normShift * normal.XYZ() < 0. )
+        normal = normShift;
+    }
     isOK = true;
   }
-  else // hard singularity
+
+  if ( !isOK ) // hard singularity, to call with shiftInside=true ?
   {
     const TGeomID faceID = helper.GetMeshDS()->ShapeToIndex( face );
 
@@ -2289,7 +2527,9 @@ gp_XYZ _ViscousBuilder::getFaceNormal(const SMDS_MeshNode* node,
         isOK = SMESH_MeshAlgos::FaceNormal( f, (gp_XYZ&) normal.XYZ(), /*normalized=*/true );
         if ( isOK )
         {
-          if ( helper.IsReversedSubMesh( face ))
+          TopoDS_Face ff = face;
+          ff.Orientation( TopAbs_FORWARD );
+          if ( helper.IsReversedSubMesh( ff ))
             normal.Reverse();
           break;
         }
@@ -2359,12 +2599,10 @@ gp_XYZ _ViscousBuilder::getWeigthedNormal( const SMDS_MeshNode*         n,
     }
     else
     {
-      TopoDS_Vertex v10 = SMESH_MesherHelper::IthVertex( 1, ee[ 0 ]);
-      TopoDS_Vertex v01 = SMESH_MesherHelper::IthVertex( 0, ee[ 1 ]);
-      if ( !v10.IsSame( v01 ))
+      if ( !V.IsSame( SMESH_MesherHelper::IthVertex( 0, ee[ 1 ] )))
         std::swap( ee[0], ee[1] );
     }
-    angles[i] = SMESH_MesherHelper::GetAngle( ee[0], ee[1], F );
+    angles[i] = SMESH_MesherHelper::GetAngle( ee[0], ee[1], F, TopoDS::Vertex( V ));
   }
 
   // compute a weighted normal
@@ -2392,11 +2630,10 @@ bool _ViscousBuilder::findNeiborsOnEdge(const _LayerEdge*     edge,
                                         _SolidData&           data)
 {
   const SMDS_MeshNode* node = edge->_nodes[0];
-  const int shapeInd = node->getshapeId();
-  SMESHDS_SubMesh* edgeSM = 0;
+  const int        shapeInd = node->getshapeId();
+  SMESHDS_SubMesh*   edgeSM = 0;
   if ( node->GetPosition()->GetTypeOfPosition() == SMDS_TOP_EDGE )
   {
-    
     edgeSM = getMeshDS()->MeshElements( shapeInd );
     if ( !edgeSM || edgeSM->NbElements() == 0 )
       return error(SMESH_Comment("Not meshed EDGE ") << shapeInd, data._index);
@@ -2444,11 +2681,11 @@ void _LayerEdge::SetDataByNeighbors( const SMDS_MeshNode* n1,
 
   // Set _curvature
 
-  double sumLen = vec1.Modulus() + vec2.Modulus();
+  double      sumLen = vec1.Modulus() + vec2.Modulus();
   _2neibors->_wgt[0] = 1 - vec1.Modulus() / sumLen;
   _2neibors->_wgt[1] = 1 - vec2.Modulus() / sumLen;
   double avgNormProj = 0.5 * ( _normal * vec1 + _normal * vec2 );
-  double avgLen = 0.5 * ( vec1.Modulus() + vec2.Modulus() );
+  double      avgLen = 0.5 * ( vec1.Modulus() + vec2.Modulus() );
   if ( _curvature ) delete _curvature;
   _curvature = _Curvature::New( avgNormProj, avgLen );
   // if ( _curvature )
@@ -2462,10 +2699,13 @@ void _LayerEdge::SetDataByNeighbors( const SMDS_MeshNode* n1,
   if ( _sWOL.IsNull() )
   {
     TopoDS_Shape S = helper.GetSubShapeByNode( _nodes[0], helper.GetMeshDS() );
-    gp_XYZ dirE = getEdgeDir( TopoDS::Edge( S ), _nodes[0], helper );
+    TopoDS_Edge  E = TopoDS::Edge( S );
+    // if ( SMESH_Algo::isDegenerated( E ))
+    //   return;
+    gp_XYZ dirE    = getEdgeDir( E, _nodes[0], helper );
     gp_XYZ plnNorm = dirE ^ _normal;
-    double proj0 = plnNorm * vec1;
-    double proj1 = plnNorm * vec2;
+    double proj0   = plnNorm * vec1;
+    double proj1   = plnNorm * vec2;
     if ( fabs( proj0 ) > 1e-10 || fabs( proj1 ) > 1e-10 )
     {
       if ( _2neibors->_plnNorm ) delete _2neibors->_plnNorm;
@@ -2632,6 +2872,31 @@ void _ViscousBuilder::makeGroupOfLE()
 #endif
 }
 
+//================================================================================
+/*!
+ * \brief Find maximal _LayerEdge length (layer thickness) limited by geometry
+ */
+//================================================================================
+
+void _ViscousBuilder::computeGeomSize( _SolidData& data )
+{
+  data._geomSize = Precision::Infinite();
+  double intersecDist;
+  auto_ptr<SMESH_ElementSearcher> searcher
+    ( SMESH_MeshAlgos::GetElementSearcher( *getMeshDS(),
+                                           data._proxyMesh->GetFaces( data._solid )) );
+
+  TNode2Edge::iterator n2e = data._n2eMap.begin(), n2eEnd = data._n2eMap.end();
+  for ( ; n2e != n2eEnd; ++n2e )
+  {
+    _LayerEdge* edge = n2e->second;
+    if ( edge->IsOnEdge() ) continue;
+    edge->FindIntersection( *searcher, intersecDist, data._epsilon );
+    if ( data._geomSize > intersecDist && intersecDist > 0 )
+      data._geomSize = intersecDist;
+  }
+}
+
 //================================================================================
 /*!
  * \brief Increase length of _LayerEdge's to reach the required thickness of layers
@@ -2644,19 +2909,8 @@ bool _ViscousBuilder::inflate(_SolidData& data)
 
   // Limit inflation step size by geometry size found by itersecting
   // normals of _LayerEdge's with mesh faces
-  double geomSize = Precision::Infinite(), intersecDist;
-  auto_ptr<SMESH_ElementSearcher> searcher
-    ( SMESH_MeshAlgos::GetElementSearcher( *getMeshDS(),
-                                           data._proxyMesh->GetFaces( data._solid )) );
-  for ( size_t i = 0; i < data._edges.size(); ++i )
-  {
-    if ( data._edges[i]->IsOnEdge() ) continue;
-    data._edges[i]->FindIntersection( *searcher, intersecDist, data._epsilon );
-    if ( geomSize > intersecDist && intersecDist > 0 )
-      geomSize = intersecDist;
-  }
-  if ( data._stepSize > 0.3 * geomSize )
-    limitStepSize( data, 0.3 * geomSize );
+  if ( data._stepSize > 0.3 * data._geomSize )
+    limitStepSize( data, 0.3 * data._geomSize );
 
   const double tgtThick = data._hyp->GetTotalThickness();
   if ( data._stepSize > tgtThick )
@@ -2665,7 +2919,7 @@ bool _ViscousBuilder::inflate(_SolidData& data)
   if ( data._stepSize < 1. )
     data._epsilon = data._stepSize * 1e-7;
 
-  debugMsg( "-- geomSize = " << geomSize << ", stepSize = " << data._stepSize );
+  debugMsg( "-- geomSize = " << data._geomSize << ", stepSize = " << data._stepSize );
 
   double avgThick = 0, curThick = 0, distToIntersection = Precision::Infinite();
   int nbSteps = 0, nbRepeats = 0;
@@ -2743,6 +2997,22 @@ bool _ViscousBuilder::inflate(_SolidData& data)
                                     " average reached thickness is " << avgThick ));
     }
 
+
+  // Restore position of src nodes moved by infaltion on _noShrinkShapes
+  dumpFunction(SMESH_Comment("restoNoShrink_So")<<data._index); // debug
+  int iBeg, iEnd = 0;
+  for ( int iS = 0; iS < data._endEdgeOnShape.size(); ++iS )
+  {
+    iBeg = iEnd;
+    iEnd = data._endEdgeOnShape[ iS ];
+    if ( data._edges[ iBeg ]->_nodes.size() == 1 )
+      for ( ; iBeg < iEnd; ++iBeg )
+      {
+        restoreNoShrink( *data._edges[ iBeg ] );
+      }
+  }
+  dumpFunctionEnd();
+
   return true;
 }
 
@@ -2818,11 +3088,15 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data,
         int oldBadNb = badNb;
         badNb = 0;
         moved = false;
-        for ( int i = iBeg; i < iEnd; ++i )
-          moved |= data._edges[i]->Smooth(badNb);
+        if ( step % 2 )
+          for ( int i = iBeg; i < iEnd; ++i )
+            moved |= data._edges[i]->Smooth(badNb);
+        else
+          for ( int i = iEnd-1; i >= iBeg; --i )
+            moved |= data._edges[i]->Smooth(badNb);
         improved = ( badNb < oldBadNb );
 
-        // issue 22576. no bad faces but still there are intersections to fix
+        // issue 22576 -- no bad faces but still there are intersections to fix
         if ( improved && badNb == 0 )
           stepLimit = step + 3;
 
@@ -2877,6 +3151,8 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data,
   int iLE = 0;
   for ( size_t i = 0; i < data._edges.size(); ++i )
   {
+    if ( !data._edges[i]->_sWOL.IsNull() )
+      continue;
     if ( data._edges[i]->FindIntersection( *searcher, dist, data._epsilon, &intFace ))
       return false;
     if ( distToIntersection > dist )
@@ -2956,8 +3232,8 @@ Handle(Geom_Curve) _SolidData::CurveForSmooth( const TopoDS_Edge&    E,
           bndBox.Add( SMESH_TNodeXYZ( nIt->next() ));
         gp_XYZ size = bndBox.CornerMax() - bndBox.CornerMin();
 
-        SMESH_TNodeXYZ p0( _edges[iFrom]->_2neibors->_nodes[0] );
-        SMESH_TNodeXYZ p1( _edges[iFrom]->_2neibors->_nodes[1] );
+        SMESH_TNodeXYZ p0( _edges[iFrom]->_2neibors->tgtNode(0) );
+        SMESH_TNodeXYZ p1( _edges[iFrom]->_2neibors->tgtNode(1) );
         const double lineTol = 1e-2 * ( p0 - p1 ).Modulus();
         for ( int i = 0; i < 3 && !isLine; ++i )
           isLine = ( size.Coord( i+1 ) <= lineTol );
@@ -3040,10 +3316,10 @@ void _SolidData::SortOnEdge( const TopoDS_Edge&  E,
 
   // set _2neibors according to the new order
   for ( int i = iFrom; i < iTo-1; ++i )
-    if ( _edges[i]->_2neibors->_nodes[1] != _edges[i+1]->_nodes.back() )
+    if ( _edges[i]->_2neibors->tgtNode(1) != _edges[i+1]->_nodes.back() )
       _edges[i]->_2neibors->reverse();
   if ( u2edge.size() > 1 &&
-       _edges[iTo-1]->_2neibors->_nodes[0] != _edges[iTo-2]->_nodes.back() )
+       _edges[iTo-1]->_2neibors->tgtNode(0) != _edges[iTo-2]->_nodes.back() )
     _edges[iTo-1]->_2neibors->reverse();
 }
 
@@ -3080,7 +3356,7 @@ bool _SolidData::GetShapeEdges(const TGeomID shapeID,
  */
 //================================================================================
 
-void _SolidData::AddFacesToSmooth( const set< TGeomID >& faceIDs )
+void _SolidData::AddShapesToSmooth( const set< TGeomID >& faceIDs )
 {
   // convert faceIDs to indices in _endEdgeOnShape
   set< size_t > iEnds;
@@ -3093,7 +3369,7 @@ void _SolidData::AddFacesToSmooth( const set< TGeomID >& faceIDs )
   set< size_t >::iterator endsIt = iEnds.begin();
 
   // "add" by move of _nbShapesToSmooth
-  int nbFacesToAdd = faceIDs.size();
+  int nbFacesToAdd = iEnds.size();
   while ( endsIt != iEnds.end() && *endsIt == _nbShapesToSmooth )
   {
     ++endsIt;
@@ -3169,8 +3445,8 @@ bool _ViscousBuilder::smoothAnalyticEdge( _SolidData&           data,
   {
     if ( F.IsNull() ) // 3D
     {
-      SMESH_TNodeXYZ p0( data._edges[iFrom]->_2neibors->_nodes[0]);
-      SMESH_TNodeXYZ p1( data._edges[iTo-1]->_2neibors->_nodes[1]);
+      SMESH_TNodeXYZ p0( data._edges[iFrom]->_2neibors->tgtNode(0));
+      SMESH_TNodeXYZ p1( data._edges[iTo-1]->_2neibors->tgtNode(1));
       for ( int i = iFrom; i < iTo; ++i )
       {
         double r = len[i-iFrom] / len.back();
@@ -3183,10 +3459,12 @@ bool _ViscousBuilder::smoothAnalyticEdge( _SolidData&           data,
     }
     else
     {
-      gp_XY uv0 = helper.GetNodeUV( F, data._edges[iFrom]->_2neibors->_nodes[0]);
-      gp_XY uv1 = helper.GetNodeUV( F, data._edges[iTo-1]->_2neibors->_nodes[1]);
-      if ( data._edges[iFrom]->_2neibors->_nodes[0] ==
-           data._edges[iTo-1]->_2neibors->_nodes[1] ) // closed edge
+      // gp_XY uv0 = helper.GetNodeUV( F, data._edges[iFrom]->_2neibors->tgtNode(0));
+      // gp_XY uv1 = helper.GetNodeUV( F, data._edges[iTo-1]->_2neibors->tgtNode(1));
+      gp_XY uv0 = data._edges[iFrom]->_2neibors->_edges[0]->LastUV( F );
+      gp_XY uv1 = data._edges[iTo-1]->_2neibors->_edges[1]->LastUV( F );
+      if ( data._edges[iFrom]->_2neibors->tgtNode(0) ==
+           data._edges[iTo-1]->_2neibors->tgtNode(1) ) // closed edge
       {
         int iPeriodic = helper.GetPeriodicIndex();
         if ( iPeriodic == 1 || iPeriodic == 2 )
@@ -3223,8 +3501,8 @@ bool _ViscousBuilder::smoothAnalyticEdge( _SolidData&           data,
 
     if ( F.IsNull() ) // 3D
     {
-      if ( data._edges[iFrom]->_2neibors->_nodes[0] ==
-           data._edges[iTo-1]->_2neibors->_nodes[1] )
+      if ( data._edges[iFrom]->_2neibors->tgtNode(0) ==
+           data._edges[iTo-1]->_2neibors->tgtNode(1) )
         return true; // closed EDGE - nothing to do
 
       return false; // TODO ???
@@ -3233,9 +3511,12 @@ bool _ViscousBuilder::smoothAnalyticEdge( _SolidData&           data,
     {
       const gp_XY center( center3D.X(), center3D.Y() );
 
-      gp_XY uv0 = helper.GetNodeUV( F, data._edges[iFrom]->_2neibors->_nodes[0]);
-      gp_XY uvM = helper.GetNodeUV( F, data._edges[iFrom]->_nodes.back());
-      gp_XY uv1 = helper.GetNodeUV( F, data._edges[iTo-1]->_2neibors->_nodes[1]);
+      gp_XY uv0 = data._edges[iFrom]->_2neibors->_edges[0]->LastUV( F );
+      gp_XY uvM = data._edges[iFrom]->LastUV( F );
+      gp_XY uv1 = data._edges[iTo-1]->_2neibors->_edges[1]->LastUV( F );
+      // gp_XY uv0 = helper.GetNodeUV( F, data._edges[iFrom]->_2neibors->tgtNode(0));
+      // gp_XY uvM = helper.GetNodeUV( F, data._edges[iFrom]->_nodes.back());
+      // gp_XY uv1 = helper.GetNodeUV( F, data._edges[iTo-1]->_2neibors->tgtNode(1));
       gp_Vec2d vec0( center, uv0 );
       gp_Vec2d vecM( center, uvM );
       gp_Vec2d vec1( center, uv1 );
@@ -3299,7 +3580,7 @@ bool _ViscousBuilder::updateNormals( _SolidData&         data,
       const SMDS_MeshNode* tgt1 = edge->_nodes.back();
       for ( int j = 0; j < 2; ++j ) // loop on _2NearEdges
       {
-        const SMDS_MeshNode* tgt2 = edge->_2neibors->_nodes[j];
+        const SMDS_MeshNode* tgt2 = edge->_2neibors->tgtNode(j);
         pair< set< SMESH_TLink >::iterator, bool > link_isnew =
           extrudedLinks.insert( SMESH_TLink( tgt1, tgt2 ));
         if ( !link_isnew.second )
@@ -3340,7 +3621,9 @@ bool _ViscousBuilder::updateNormals( _SolidData&         data,
   for ( size_t i = 0; i < data._edges.size(); ++i )
   {
     _LayerEdge* edge = data._edges[i];
-    if ( !edge->IsOnEdge() || !edge->_sWOL.IsNull() ) continue;
+    if (( !edge->IsOnEdge() ) &&
+        ( edge->_sWOL.IsNull() || edge->_sWOL.ShapeType() != TopAbs_FACE ))
+      continue;
     if ( edge->FindIntersection( *searcher, dist, eps, &face ))
     {
       const _TmpMeshFaceOnEdge* f = (const _TmpMeshFaceOnEdge*) face;
@@ -3360,107 +3643,163 @@ bool _ViscousBuilder::updateNormals( _SolidData&         data,
   {
     dumpFunction(SMESH_Comment("updateNormals")<<data._index);
 
+    set< TGeomID > shapesToSmooth;
+
+    // vector to store new _normal and _cosin for each edge in edge2CloseEdge
+    vector< pair< _LayerEdge*, _LayerEdge > > edge2newEdge( edge2CloseEdge.size() );
+
     TLEdge2LEdgeSet::iterator e2ee = edge2CloseEdge.begin();
-    for ( ; e2ee != edge2CloseEdge.end(); ++e2ee )
+    for ( size_t iE = 0; e2ee != edge2CloseEdge.end(); ++e2ee, ++iE )
     {
-      _LayerEdge* edge1       = e2ee->first;
-      _LayerEdge* edge2       = 0;
-      set< _LayerEdge*, _LayerEdgeCmp >& ee  = e2ee->second;
+      _LayerEdge* edge1 = e2ee->first;
+      _LayerEdge* edge2 = 0;
+      set< _LayerEdge*, _LayerEdgeCmp >& ee = e2ee->second;
+
+      edge2newEdge[ iE ].first = NULL;
 
       // find EDGEs the edges reside
-      TopoDS_Edge E1, E2;
-      TopoDS_Shape S = helper.GetSubShapeByNode( edge1->_nodes[0], getMeshDS() );
-      if ( S.ShapeType() != TopAbs_EDGE )
-        continue; // TODO: find EDGE by VERTEX
-      E1 = TopoDS::Edge( S );
+      // TopoDS_Edge E1, E2;
+      // TopoDS_Shape S = helper.GetSubShapeByNode( edge1->_nodes[0], getMeshDS() );
+      // if ( S.ShapeType() != TopAbs_EDGE )
+      //   continue; // TODO: find EDGE by VERTEX
+      // E1 = TopoDS::Edge( S );
       set< _LayerEdge*, _LayerEdgeCmp >::iterator eIt = ee.begin();
-      while ( E2.IsNull() && eIt != ee.end())
+      for ( ; !edge2 && eIt != ee.end(); ++eIt )
       {
-        _LayerEdge* e2 = *eIt++;
-        TopoDS_Shape S = helper.GetSubShapeByNode( e2->_nodes[0], getMeshDS() );
-        if ( S.ShapeType() == TopAbs_EDGE )
-          E2 = TopoDS::Edge( S ), edge2 = e2;
+        if ( edge1->_sWOL == (*eIt)->_sWOL )
+          edge2 = *eIt;
       }
-      if ( E2.IsNull() ) continue; // TODO: find EDGE by VERTEX
+      if ( !edge2 ) continue;
+
+      edge2newEdge[ iE ].first = edge1;
+      _LayerEdge& newEdge = edge2newEdge[ iE ].second;
+      // while ( E2.IsNull() && eIt != ee.end())
+      // {
+      //   _LayerEdge* e2 = *eIt++;
+      //   TopoDS_Shape S = helper.GetSubShapeByNode( e2->_nodes[0], getMeshDS() );
+      //   if ( S.ShapeType() == TopAbs_EDGE )
+      //     E2 = TopoDS::Edge( S ), edge2 = e2;
+      // }
+      // if ( E2.IsNull() ) continue; // TODO: find EDGE by VERTEX
 
       // find 3 FACEs sharing 2 EDGEs
 
-      TopoDS_Face FF1[2], FF2[2];
-      PShapeIteratorPtr fIt = helper.GetAncestors(E1, *_mesh, TopAbs_FACE);
-      while ( fIt->more() && FF1[1].IsNull())
+      // TopoDS_Face FF1[2], FF2[2];
+      // PShapeIteratorPtr fIt = helper.GetAncestors(E1, *_mesh, TopAbs_FACE);
+      // while ( fIt->more() && FF1[1].IsNull() )
+      // {
+      //   const TopoDS_Face *F = (const TopoDS_Face*) fIt->next();
+      //   if ( helper.IsSubShape( *F, data._solid))
+      //     FF1[ FF1[0].IsNull() ? 0 : 1 ] = *F;
+      // }
+      // fIt = helper.GetAncestors(E2, *_mesh, TopAbs_FACE);
+      // while ( fIt->more() && FF2[1].IsNull())
+      // {
+      //   const TopoDS_Face *F = (const TopoDS_Face*) fIt->next();
+      //   if ( helper.IsSubShape( *F, data._solid))
+      //     FF2[ FF2[0].IsNull() ? 0 : 1 ] = *F;
+      // }
+      // // exclude a FACE common to E1 and E2 (put it to FFn[1] )
+      // if ( FF1[0].IsSame( FF2[0]) || FF1[0].IsSame( FF2[1]))
+      //   std::swap( FF1[0], FF1[1] );
+      // if ( FF2[0].IsSame( FF1[0]) )
+      //   std::swap( FF2[0], FF2[1] );
+      // if ( FF1[0].IsNull() || FF2[0].IsNull() )
+      //   continue;
+
+      // get a new normal for edge1
+      //bool ok;
+      gp_Vec dir1 = edge1->_normal, dir2 = edge2->_normal;
+      // if ( edge1->_cosin < 0 )
+      //   dir1 = getFaceDir( FF1[0], E1, edge1->_nodes[0], helper, ok ).Normalized();
+      // if ( edge2->_cosin < 0 )
+      //   dir2 = getFaceDir( FF2[0], E2, edge2->_nodes[0], helper, ok ).Normalized();
+
+      double cos1 = Abs( edge1->_cosin ), cos2 = Abs( edge2->_cosin );
+      double wgt1 = ( cos1 + 0.001 ) / ( cos1 + cos2 + 0.002 );
+      double wgt2 = ( cos2 + 0.001 ) / ( cos1 + cos2 + 0.002 );
+      newEdge._normal = ( wgt1 * dir1 + wgt2 * dir2 ).XYZ();
+      newEdge._normal.Normalize();
+
+      // cout << edge1->_nodes[0]->GetID() << " "
+      //      << edge2->_nodes[0]->GetID() << " NORM: "
+      //      << newEdge._normal.X() << ", " << newEdge._normal.Y() << ", " << newEdge._normal.Z() << endl;
+
+      // get new cosin
+      if ( cos1 < theMinSmoothCosin )
       {
-        const TopoDS_Face *F = (const TopoDS_Face*) fIt->next();
-        if ( helper.IsSubShape( *F, data._solid))
-          FF1[ FF1[0].IsNull() ? 0 : 1 ] = *F;
+        newEdge._cosin = edge2->_cosin;
       }
-      fIt = helper.GetAncestors(E2, *_mesh, TopAbs_FACE);
-      while ( fIt->more() && FF2[1].IsNull())
+      else if ( cos2 > theMinSmoothCosin ) // both cos1 and cos2 > theMinSmoothCosin
       {
-        const TopoDS_Face *F = (const TopoDS_Face*) fIt->next();
-        if ( helper.IsSubShape( *F, data._solid))
-          FF2[ FF2[0].IsNull() ? 0 : 1 ] = *F;
+        // gp_Vec dirInFace;
+        // if ( edge1->_cosin < 0 )
+        //   dirInFace = dir1;
+        // else
+        //   dirInFace = getFaceDir( FF1[0], E1, edge1->_nodes[0], helper, ok );
+        // double angle = dirInFace.Angle( edge1->_normal ); // [0,PI]
+        // edge1->SetCosin( Cos( angle ));
+        //newEdge._cosin = 0; // ???????????
+        newEdge._cosin = ( wgt1 * cos1 + wgt2 * cos2 ) * edge1->_cosin / cos1;
       }
-      // exclude a FACE common to E1 and E2 (put it at [1] in FF* )
-      if ( FF1[0].IsSame( FF2[0]) || FF1[0].IsSame( FF2[1]))
-        std::swap( FF1[0], FF1[1] );
-      if ( FF2[0].IsSame( FF1[0]) )
-        std::swap( FF2[0], FF2[1] );
-      if ( FF1[0].IsNull() || FF2[0].IsNull() )
-        continue;
-
-      // get a new normal for edge1
-      bool ok;
-      gp_Vec dir1 = edge1->_normal, dir2 = edge2->_normal;
-      if ( edge1->_cosin < 0 )
-        dir1 = getFaceDir( FF1[0], E1, edge1->_nodes[0], helper, ok ).Normalized();
-      if ( edge2->_cosin < 0 )
-        dir2 = getFaceDir( FF2[0], E2, edge2->_nodes[0], helper, ok ).Normalized();
-      // gp_Vec dir1 = getFaceDir( FF1[0], E1, edge1->_nodes[0], helper, ok );
-      // gp_Vec dir2 = getFaceDir( FF2[0], E2, edge2->_nodes[0], helper, ok2 );
-      // double wgt1 = ( edge1->_cosin + 1 ) / ( edge1->_cosin + edge2->_cosin + 2 );
-      // double wgt2 = ( edge2->_cosin + 1 ) / ( edge1->_cosin + edge2->_cosin + 2 );
-      // gp_Vec newNorm = wgt1 * dir1 + wgt2 * dir2;
-      // newNorm.Normalize();
-
-      double wgt1 = ( edge1->_cosin + 1 ) / ( edge1->_cosin + edge2->_cosin + 2 );
-      double wgt2 = ( edge2->_cosin + 1 ) / ( edge1->_cosin + edge2->_cosin + 2 );
-      gp_Vec newNorm = wgt1 * dir1 + wgt2 * dir2;
-      newNorm.Normalize();
-
-      edge1->_normal = newNorm.XYZ();
-
-      // update data of edge1 depending on _normal
-      const SMDS_MeshNode *n1, *n2;
-      n1 = edge1->_2neibors->_edges[0]->_nodes[0];
-      n2 = edge1->_2neibors->_edges[1]->_nodes[0];
-      edge1->SetDataByNeighbors( n1, n2, helper );
-      gp_Vec dirInFace;
-      if ( edge1->_cosin < 0 )
-        dirInFace = dir1;
       else
-        dirInFace = getFaceDir( FF1[0], E1, edge1->_nodes[0], helper, ok );
-      double angle = dir1.Angle( edge1->_normal ); // [0,PI]
-      edge1->SetCosin( cos( angle ));
+      {
+        newEdge._cosin = edge1->_cosin;
+      }
 
-      // limit data._stepSize
-      if ( edge1->_cosin > theMinSmoothCosin )
+      // find shapes that need smoothing due to change of _normal
+      if ( edge1->_cosin  < theMinSmoothCosin &&
+           newEdge._cosin > theMinSmoothCosin )
       {
-        SMDS_ElemIteratorPtr fIt = edge1->_nodes[0]->GetInverseElementIterator(SMDSAbs_Face);
-        while ( fIt->more() )
-          limitStepSize( data, fIt->next(), edge1->_cosin );
+        if ( edge1->_sWOL.IsNull() )
+        {
+          SMDS_ElemIteratorPtr fIt = edge1->_nodes[0]->GetInverseElementIterator(SMDSAbs_Face);
+          while ( fIt->more() )
+            shapesToSmooth.insert( fIt->next()->getshapeId() );
+          //limitStepSize( data, fIt->next(), edge1->_cosin ); // too late
+        }
+        else // edge1 inflates along a FACE
+        {
+          TopoDS_Shape V = helper.GetSubShapeByNode( edge1->_nodes[0], getMeshDS() );
+          PShapeIteratorPtr eIt = helper.GetAncestors( V, *_mesh, TopAbs_EDGE );
+          while ( const TopoDS_Shape* E = eIt->next() )
+          {
+            if ( !helper.IsSubShape( *E, /*FACE=*/edge1->_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 )
+              shapesToSmooth.insert( getMeshDS()->ShapeToIndex( *E ));
+          }
+        }
       }
-      // set new XYZ of target node
+    }
+
+    data.AddShapesToSmooth( shapesToSmooth );
+
+    // Update data of edges depending on a new _normal
+
+    for ( size_t iE = 0; iE < edge2newEdge.size(); ++iE )
+    {
+      _LayerEdge*   edge1 = edge2newEdge[ iE ].first;
+      _LayerEdge& newEdge = edge2newEdge[ iE ].second;
+      if ( !edge1 ) continue;
+
+      edge1->_normal = newEdge._normal;
+      edge1->SetCosin( newEdge._cosin );
       edge1->InvalidateStep( 1 );
       edge1->_len = 0;
       edge1->SetNewLength( data._stepSize, helper );
-    }
+      if ( edge1->IsOnEdge() )
+      {
+        const SMDS_MeshNode * n1 = edge1->_2neibors->srcNode(0);
+        const SMDS_MeshNode * n2 = edge1->_2neibors->srcNode(1);
+        edge1->SetDataByNeighbors( n1, n2, helper );
+      }
 
-    // Update normals and other dependent data of not intersecting _LayerEdge's
-    // neighboring the intersecting ones
+      // Update normals and other dependent data of not intersecting _LayerEdge's
+      // neighboring the intersecting ones
 
-    for ( e2ee = edge2CloseEdge.begin(); e2ee != edge2CloseEdge.end(); ++e2ee )
-    {
-      _LayerEdge* edge1 = e2ee->first;
       if ( !edge1->_2neibors )
         continue;
       for ( int j = 0; j < 2; ++j ) // loop on 2 neighbors
@@ -3469,7 +3808,7 @@ bool _ViscousBuilder::updateNormals( _SolidData&         data,
         if ( edge2CloseEdge.count ( neighbor ))
           continue; // j-th neighbor is also intersected
         _LayerEdge* prevEdge = edge1;
-        const int nbSteps = 6;
+        const int nbSteps = 10;
         for ( int step = nbSteps; step; --step ) // step from edge1 in j-th direction
         {
           if ( !neighbor->_2neibors )
@@ -3598,7 +3937,13 @@ bool _ViscousBuilder::updateNormalsOfConvexFaces( _SolidData&         data,
         for ( ; iBeg < iEnd; ++iBeg )
           avgNormal += data._edges[ iBeg ]->_normal;
       }
-      avgNormal.Normalize();
+      double normSize = avgNormal.SquareModulus();
+      if ( normSize < 1e-200 )
+      {
+        debugMsg( "updateNormalsOfConvexFaces(): zero avgNormal" );
+        return false;
+      }
+      avgNormal /= Sqrt( normSize );
 
       // compute new _LayerEdge::_cosin on EDGEs
       double avgCosin = 0;
@@ -3857,7 +4202,7 @@ bool _ViscousBuilder::updateNormalsOfConvexFaces( _SolidData&         data,
         }
       }
     }
-    data.AddFacesToSmooth( adjFacesToSmooth );
+    data.AddShapesToSmooth( adjFacesToSmooth );
 
     dumpFunctionEnd();
 
@@ -3957,7 +4302,7 @@ bool _CentralCurveOnEdge::FindNewNormal( const gp_Pnt& center, gp_XYZ& newNormal
 
     newNormal += norm;
     double sz = newNormal.Modulus();
-    if ( Abs ( sz ) < 1e-200 )
+    if ( sz < 1e-200 )
       break;
     newNormal /= sz;
     return true;
@@ -4013,7 +4358,7 @@ bool _LayerEdge::FindIntersection( SMESH_ElementSearcher&   searcher,
   bool segmentIntersected = false;
   distance = Precision::Infinite();
   int iFace = -1; // intersected face
-  for ( size_t j = 0 ; j < suspectFaces.size() && !segmentIntersected; ++j )
+  for ( size_t j = 0 ; j < suspectFaces.size() /*&& !segmentIntersected*/; ++j )
   {
     const SMDS_MeshElement* face = suspectFaces[j];
     if ( face->GetNodeIndex( _nodes.back() ) >= 0 ||
@@ -4031,7 +4376,7 @@ bool _LayerEdge::FindIntersection( SMESH_ElementSearcher&   searcher,
     {
       const SMDS_MeshNode* tria[3];
       tria[0] = *nIt++;
-      tria[1] = *nIt++;;
+      tria[1] = *nIt++;
       for ( int n2 = 2; n2 < nbNodes && !intFound; ++n2 )
       {
         tria[2] = *nIt++;
@@ -4041,7 +4386,7 @@ bool _LayerEdge::FindIntersection( SMESH_ElementSearcher&   searcher,
     }
     if ( intFound )
     {
-      if ( dist < segLen*(1.01) && dist > -(_len-segLen) )
+      if ( dist < segLen*(1.01) && dist > -(_len*_lenFactor-segLen) )
         segmentIntersected = true;
       if ( distance > dist )
         distance = dist, iFace = j;
@@ -4123,6 +4468,31 @@ gp_Ax1 _LayerEdge::LastSegment(double& segLen) const
   return segDir;
 }
 
+//================================================================================
+/*!
+ * \brief Return the last position of the target node on a FACE. 
+ *  \param [in] F - the FACE this _LayerEdge is inflated along
+ *  \return gp_XY - result UV
+ */
+//================================================================================
+
+gp_XY _LayerEdge::LastUV( const TopoDS_Face& F ) const
+{
+  if ( F.IsSame( _sWOL )) // F is my FACE
+    return gp_XY( _pos.back().X(), _pos.back().Y() );
+
+  if ( _sWOL.IsNull() || _sWOL.ShapeType() != TopAbs_EDGE ) // wrong call
+    return gp_XY( 1e100, 1e100 );
+
+  // _sWOL is EDGE of F; _pos.back().X() is the last U on the EDGE
+  double f, l, u = _pos.back().X();
+  Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface( TopoDS::Edge(_sWOL), F, f,l);
+  if ( !C2d.IsNull() && f <= u && u <= l )
+    return C2d->Value( u ).XY();
+
+  return gp_XY( 1e100, 1e100 );
+}
+
 //================================================================================
 /*!
  * \brief Test intersection of the last segment with a given triangle
@@ -4150,9 +4520,9 @@ bool _LayerEdge::SegTriaInter( const gp_Ax1&        lastSegment,
   /* calculate distance from vert0 to ray origin */
   gp_XYZ tvec = orig - vert0;
 
-  if ( tvec * dir > EPSILON )
+  //if ( tvec * dir > EPSILON )
     // intersected face is at back side of the temporary face this _LayerEdge belongs to
-    return false;
+    //return false;
 
   gp_XYZ edge1 = vert1 - vert0;
   gp_XYZ edge2 = vert2 - vert0;
@@ -4164,26 +4534,29 @@ bool _LayerEdge::SegTriaInter( const gp_Ax1&        lastSegment,
   double det = edge1 * pvec;
 
   if (det > -EPSILON && det < EPSILON)
-    return 0;
+    return false;
   double inv_det = 1.0 / det;
 
   /* calculate U parameter and test bounds */
   double u = ( tvec * pvec ) * inv_det;
-  if (u < 0.0 || u > 1.0)
-    return 0;
+  //if (u < 0.0 || u > 1.0)
+  if (u < -EPSILON || u > 1.0 + EPSILON)
+    return false;
 
   /* prepare to test V parameter */
   gp_XYZ qvec = tvec ^ edge1;
 
   /* calculate V parameter and test bounds */
   double v = (dir * qvec) * inv_det;
-  if ( v < 0.0 || u + v > 1.0 )
-    return 0;
+  //if ( v < 0.0 || u + v > 1.0 )
+  if ( v < -EPSILON || u + v > 1.0 + EPSILON)
+    return false;
 
   /* calculate t, ray intersects triangle */
   t = (edge2 * qvec) * inv_det;
 
-  return true;
+  //return true;
+  return t > 0.;
 }
 
 //================================================================================
@@ -4203,9 +4576,9 @@ bool _LayerEdge::SmoothOnEdge(Handle(Geom_Surface)& surface,
   SMESH_TNodeXYZ oldPos( tgtNode );
   double dist01, distNewOld;
   
-  SMESH_TNodeXYZ p0( _2neibors->_nodes[0]);
-  SMESH_TNodeXYZ p1( _2neibors->_nodes[1]);
-  dist01 = p0.Distance( _2neibors->_nodes[1] );
+  SMESH_TNodeXYZ p0( _2neibors->tgtNode(0));
+  SMESH_TNodeXYZ p1( _2neibors->tgtNode(1));
+  dist01 = p0.Distance( _2neibors->tgtNode(1) );
 
   gp_Pnt newPos = p0 * _2neibors->_wgt[0] + p1 * _2neibors->_wgt[1];
   double lenDelta = 0;
@@ -4272,16 +4645,27 @@ bool _LayerEdge::Smooth(int& badNb)
     newPos += SMESH_TNodeXYZ( _simplices[i]._nPrev );
   newPos /= _simplices.size();
 
+  const gp_XYZ& curPos ( _pos.back() );
+  const gp_Pnt  prevPos( _pos[ _pos.size()-2 ]);
   if ( _curvature )
-    newPos += _normal * _curvature->lenDelta( _len );
-
-  gp_Pnt prevPos( _pos[ _pos.size()-2 ]);
+  {
+    double delta  = _curvature->lenDelta( _len );
+    if ( delta > 0 )
+      newPos += _normal * delta;
+    else
+    {
+      double segLen = _normal * ( newPos - prevPos.XYZ() );
+      if ( segLen + delta > 0 )
+        newPos += _normal * delta;
+    }
+    // double segLenChange = _normal * ( curPos - newPos );
+    // newPos += 0.5 * _normal * segLenChange;
+  }
 
   // count quality metrics (orientation) of tetras around _tgtNode
   int nbOkBefore = 0;
-  SMESH_TNodeXYZ tgtXYZ( _nodes.back() );
   for ( size_t i = 0; i < _simplices.size(); ++i )
-    nbOkBefore += _simplices[i].IsForward( _nodes[0], &tgtXYZ );
+    nbOkBefore += _simplices[i].IsForward( _nodes[0], &curPos );
 
   int nbOkAfter = 0;
   for ( size_t i = 0; i < _simplices.size(); ++i )
@@ -4334,17 +4718,23 @@ void _LayerEdge::SetNewLength( double len, SMESH_MesherHelper& helper )
       double u = Precision::Infinite(); // to force projection w/o distance check
       helper.CheckNodeU( TopoDS::Edge( _sWOL ), n, u, 1e-10, /*force=*/true, distXYZ );
       _pos.back().SetCoord( u, 0, 0 );
-      SMDS_EdgePosition* pos = static_cast<SMDS_EdgePosition*>( n->GetPosition() );
-      pos->SetUParameter( u );
+      if ( _nodes.size() > 1 )
+      {
+        SMDS_EdgePosition* pos = static_cast<SMDS_EdgePosition*>( n->GetPosition() );
+        pos->SetUParameter( u );
+      }
     }
     else //  TopAbs_FACE
     {
       gp_XY uv( Precision::Infinite(), 0 );
       helper.CheckNodeUV( TopoDS::Face( _sWOL ), n, uv, 1e-10, /*force=*/true, distXYZ );
       _pos.back().SetCoord( uv.X(), uv.Y(), 0 );
-      SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( n->GetPosition() );
-      pos->SetUParameter( uv.X() );
-      pos->SetVParameter( uv.Y() );
+      if ( _nodes.size() > 1 )
+      {
+        SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( n->GetPosition() );
+        pos->SetUParameter( uv.X() );
+        pos->SetVParameter( uv.Y() );
+      }
     }
     n->setXYZ( distXYZ[1], distXYZ[2], distXYZ[3]);
   }
@@ -4417,10 +4807,15 @@ bool _ViscousBuilder::refine(_SolidData& data)
   TNode2Edge* n2eMap = 0;
   TNode2Edge::iterator n2e;
 
+  // Create intermediate nodes on each _LayerEdge
+
   for ( size_t i = 0; i < data._edges.size(); ++i )
   {
     _LayerEdge& edge = *data._edges[i];
 
+    if ( edge._nodes.size() < 2 )
+      continue; // on _noShrinkShapes
+
     // get accumulated length of segments
     vector< double > segLen( edge._pos.size() );
     segLen[0] = 0.0;
@@ -4459,21 +4854,22 @@ bool _ViscousBuilder::refine(_SolidData& data)
       n2eMap = ( s2ne == data._s2neMap.end() ) ? 0 : n2eMap = s2ne->second;
       prevBaseId = baseShapeId;
     }
+    _LayerEdge* edgeOnSameNode = 0;
     if ( n2eMap && (( n2e = n2eMap->find( edge._nodes[0] )) != n2eMap->end() ))
     {
-      _LayerEdge*  foundEdge = n2e->second;
-      const gp_XYZ& foundPos = foundEdge->_pos.back();
-      SMDS_PositionPtr lastPos = tgtNode->GetPosition();
+      edgeOnSameNode = n2e->second;
+      const gp_XYZ& otherTgtPos = edgeOnSameNode->_pos.back();
+      SMDS_PositionPtr  lastPos = tgtNode->GetPosition();
       if ( isOnEdge )
       {
         SMDS_EdgePosition* epos = static_cast<SMDS_EdgePosition*>( lastPos );
-        epos->SetUParameter( foundPos.X() );
+        epos->SetUParameter( otherTgtPos.X() );
       }
       else
       {
         SMDS_FacePosition* fpos = static_cast<SMDS_FacePosition*>( lastPos );
-        fpos->SetUParameter( foundPos.X() );
-        fpos->SetVParameter( foundPos.Y() );
+        fpos->SetUParameter( otherTgtPos.X() );
+        fpos->SetVParameter( otherTgtPos.Y() );
       }
     }
     // calculate height of the first layer
@@ -4505,7 +4901,7 @@ bool _ViscousBuilder::refine(_SolidData& data)
       double r = ( segLen[iSeg] - hSum ) / ( segLen[iSeg] - segLen[iPrevSeg] );
       gp_Pnt pos = r * edge._pos[iPrevSeg] + (1-r) * edge._pos[iSeg];
 
-      SMDS_MeshNode*& node = const_cast< SMDS_MeshNode*& >(edge._nodes[ iStep ]);
+      SMDS_MeshNode*& node = const_cast< SMDS_MeshNode*& >( edge._nodes[ iStep ]);
       if ( !edge._sWOL.IsNull() )
       {
         // compute XYZ by parameters <pos>
@@ -4563,8 +4959,20 @@ bool _ViscousBuilder::refine(_SolidData& data)
         }
         node->setXYZ( pos.X(), pos.Y(), pos.Z() );
       }
+    } // loop on edge._nodes
+
+    if ( !edge._sWOL.IsNull() ) // prepare for shrink()
+    {
+      if ( isOnEdge )
+        edge._pos.back().SetCoord( u, 0,0);
+      else
+        edge._pos.back().SetCoord( uv.X(), uv.Y() ,0);
+
+      if ( edgeOnSameNode )
+        edgeOnSameNode->_pos.back() = edge._pos.back();
     }
-  }
+
+  } // loop on data._edges to create nodes
 
   if ( !getMeshDS()->IsEmbeddedMode() )
     // Log node movement
@@ -4575,50 +4983,151 @@ bool _ViscousBuilder::refine(_SolidData& data)
       getMeshDS()->MoveNode( p._node, p.X(), p.Y(), p.Z() );
     }
 
-  // TODO: make quadratic prisms and polyhedrons(?)
+  // Create volumes
 
   helper.SetElementsOnShape(true);
 
+  vector< vector<const SMDS_MeshNode*>* > nnVec;
+  set< vector<const SMDS_MeshNode*>* >    nnSet;
+  set< int > degenEdgeInd;
+  vector<const SMDS_MeshElement*> degenVols;
+
   TopExp_Explorer exp( data._solid, TopAbs_FACE );
   for ( ; exp.More(); exp.Next() )
   {
     if ( data._ignoreFaceIds.count( getMeshDS()->ShapeToIndex( exp.Current() )))
       continue;
-    SMESHDS_SubMesh* fSubM = getMeshDS()->MeshElements( exp.Current() );
+    SMESHDS_SubMesh*   fSubM = getMeshDS()->MeshElements( exp.Current() );
     SMDS_ElemIteratorPtr fIt = fSubM->GetElements();
-    vector< vector<const SMDS_MeshNode*>* > nnVec;
     while ( fIt->more() )
     {
       const SMDS_MeshElement* face = fIt->next();
-      int nbNodes = face->NbCornerNodes();
+      const int            nbNodes = face->NbCornerNodes();
       nnVec.resize( nbNodes );
-      SMDS_ElemIteratorPtr nIt = face->nodesIterator();
+      nnSet.clear();
+      degenEdgeInd.clear();
+      int nbZ = 0;
+      SMDS_NodeIteratorPtr nIt = face->nodeIterator();
       for ( int iN = 0; iN < nbNodes; ++iN )
       {
-        const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nIt->next() );
+        const SMDS_MeshNode* n = nIt->next();
         nnVec[ iN ] = & data._n2eMap[ n ]->_nodes;
+        if ( nnVec[ iN ]->size() < 2 )
+          degenEdgeInd.insert( iN );
+        else
+          nbZ = nnVec[ iN ]->size();
+
+        if ( helper.HasDegeneratedEdges() )
+          nnSet.insert( nnVec[ iN ]);
       }
+      if ( nbZ == 0 )
+        continue;
+      if ( 0 < nnSet.size() && nnSet.size() < 3 )
+        continue;
 
-      int nbZ = nnVec[0]->size();
       switch ( nbNodes )
       {
       case 3:
-        for ( int iZ = 1; iZ < nbZ; ++iZ )
-          helper.AddVolume( (*nnVec[0])[iZ-1], (*nnVec[1])[iZ-1], (*nnVec[2])[iZ-1],
-                            (*nnVec[0])[iZ],   (*nnVec[1])[iZ],   (*nnVec[2])[iZ]);
+        switch ( degenEdgeInd.size() )
+        {
+        case 0: // PENTA
+        {
+          for ( int iZ = 1; iZ < nbZ; ++iZ )
+            helper.AddVolume( (*nnVec[0])[iZ-1], (*nnVec[1])[iZ-1], (*nnVec[2])[iZ-1],
+                              (*nnVec[0])[iZ],   (*nnVec[1])[iZ],   (*nnVec[2])[iZ]);
+          break;
+        }
+        case 1: // PYRAM
+        {
+          int i2 = *degenEdgeInd.begin();
+          int i0 = helper.WrapIndex( i2 - 1, nbNodes );
+          int i1 = helper.WrapIndex( i2 + 1, nbNodes );
+          for ( int iZ = 1; iZ < nbZ; ++iZ )
+            helper.AddVolume( (*nnVec[i0])[iZ-1], (*nnVec[i1])[iZ-1],
+                              (*nnVec[i1])[iZ],   (*nnVec[i0])[iZ],   (*nnVec[i2])[0]);
+          break;
+        }
+        case 2: // TETRA
+        {
+          int i3 = !degenEdgeInd.count(0) ? 0 : !degenEdgeInd.count(1) ? 1 : 2;
+          for ( int iZ = 1; iZ < nbZ; ++iZ )
+            helper.AddVolume( (*nnVec[0])[iZ-1], (*nnVec[1])[iZ-1], (*nnVec[2])[iZ-1],
+                              (*nnVec[i3])[iZ]);
+          break;
+        }
+        }
         break;
+
       case 4:
-        for ( int iZ = 1; iZ < nbZ; ++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]);
+        switch ( degenEdgeInd.size() )
+        {
+        case 0: // HEX
+        {
+          for ( int iZ = 1; iZ < nbZ; ++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]);
+          break;
+        }
+        case 2: // PENTA?
+        {
+          int i2 = *degenEdgeInd.begin();
+          int i3 = *degenEdgeInd.rbegin();
+          bool ok = ( i3 - i2 == 1 );
+          if ( i2 == 0 && i3 == 3 ) { i2 = 3; i3 = 0; ok = true; }
+          int i0 = helper.WrapIndex( i3 + 1, nbNodes );
+          int i1 = helper.WrapIndex( i0 + 1, nbNodes );
+          for ( int iZ = 1; iZ < nbZ; ++iZ )
+          {
+            const SMDS_MeshElement* vol =
+              helper.AddVolume( (*nnVec[i3])[0], (*nnVec[i0])[iZ], (*nnVec[i0])[iZ-1],
+                                (*nnVec[i2])[0], (*nnVec[i1])[iZ], (*nnVec[i1])[iZ-1]);
+            if ( !ok && vol )
+              degenVols.push_back( vol );
+          }
+          break;
+        }
+        case 3: // degen HEX
+        {
+          const SMDS_MeshNode* nn[8];
+          for ( int iZ = 1; iZ < nbZ; ++iZ )
+          {
+            const SMDS_MeshElement* vol =
+              helper.AddVolume( nnVec[0]->size() > 1 ? (*nnVec[0])[iZ-1] : (*nnVec[0])[0],
+                                nnVec[1]->size() > 1 ? (*nnVec[1])[iZ-1] : (*nnVec[1])[0],
+                                nnVec[2]->size() > 1 ? (*nnVec[2])[iZ-1] : (*nnVec[2])[0],
+                                nnVec[3]->size() > 1 ? (*nnVec[3])[iZ-1] : (*nnVec[3])[0],
+                                nnVec[0]->size() > 1 ? (*nnVec[0])[iZ]   : (*nnVec[0])[0],
+                                nnVec[1]->size() > 1 ? (*nnVec[1])[iZ]   : (*nnVec[1])[0],
+                                nnVec[2]->size() > 1 ? (*nnVec[2])[iZ]   : (*nnVec[2])[0],
+                                nnVec[3]->size() > 1 ? (*nnVec[3])[iZ]   : (*nnVec[3])[0]);
+            degenVols.push_back( vol );
+          }
+        }
         break;
+        }
+        break;
+
       default:
         return error("Not supported type of element", data._index);
-      }
+
+      } // switch ( nbNodes )
+    } // while ( fIt->more() )
+  } // loop on FACEs
+
+  if ( !degenVols.empty() )
+  {
+    SMESH_ComputeErrorPtr& err = _mesh->GetSubMesh( data._solid )->GetComputeError();
+    if ( !err || err->IsOK() )
+    {
+      err.reset( new SMESH_ComputeError( COMPERR_WARNING,
+                                         "Degenerated volumes created" ));
+      err->myBadElements.insert( err->myBadElements.end(),
+                                 degenVols.begin(),degenVols.end() );
     }
   }
+
   return true;
 }
 
@@ -4664,20 +5173,19 @@ bool _ViscousBuilder::shrink()
 
   // EDGE's to shrink
   map< TGeomID, _Shrinker1D > e2shrMap;
+  vector< _LayerEdge* > lEdges;
 
   // loop on FACES to srink mesh on
   map< TGeomID, _SolidData* >::iterator f2sd = f2sdMap.begin();
   for ( ; f2sd != f2sdMap.end(); ++f2sd )
   {
-    _SolidData&     data = *f2sd->second;
-    TNode2Edge&   n2eMap = data._n2eMap;
-    const TopoDS_Face& F = TopoDS::Face( getMeshDS()->IndexToShape( f2sd->first ));
-
-    Handle(Geom_Surface) surface = BRep_Tool::Surface(F);
-
+    _SolidData&      data = *f2sd->second;
+    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);
+
     helper.SetSubShape(F);
 
     // ===========================
@@ -4712,23 +5220,30 @@ bool _ViscousBuilder::shrink()
     }
 
     // Find _LayerEdge's inflated along F
-    vector< _LayerEdge* > lEdges;
+    lEdges.clear();
     {
-      SMESH_subMeshIteratorPtr subIt =
-        sm->getDependsOnIterator(/*includeSelf=*/false, /*complexShapeFirst=*/false);
+      set< TGeomID > subIDs;
+      SMESH_subMeshIteratorPtr subIt = sm->getDependsOnIterator(/*includeSelf=*/false);
       while ( subIt->more() )
+        subIDs.insert( subIt->next()->GetId() );
+
+      int iBeg, iEnd = 0;
+      for ( int iS = 0; iS < data._endEdgeOnShape.size() && !subIDs.empty(); ++iS )
       {
-        SMESH_subMesh*     sub = subIt->next();
-        SMESHDS_SubMesh* subDS = sub->GetSubMeshDS();
-        if ( subDS->NbNodes() == 0 || !n2eMap.count( subDS->GetNodes()->next() ))
-          continue;
-        SMDS_NodeIteratorPtr nIt = subDS->GetNodes();
-        while ( nIt->more() )
-        {
-          _LayerEdge* edge = n2eMap[ nIt->next() ];
-          lEdges.push_back( edge );
-          prepareEdgeToShrink( *edge, F, helper, smDS );
-        }
+        iBeg = iEnd;
+        iEnd = data._endEdgeOnShape[ iS ];
+        TGeomID shapeID = data._edges[ iBeg ]->_nodes[0]->getshapeId();
+        set< TGeomID >::iterator idIt = subIDs.find( shapeID );
+        if ( idIt == subIDs.end() ||
+             data._edges[ iBeg ]->_sWOL.IsNull() ) continue;
+        subIDs.erase( idIt );
+
+        if ( !data._noShrinkShapes.count( shapeID ))
+          for ( ; iBeg < iEnd; ++iBeg )
+          {
+            lEdges.push_back( data._edges[ iBeg ] );
+            prepareEdgeToShrink( *data._edges[ iBeg ], F, helper, smDS );
+          }
       }
     }
 
@@ -4739,6 +5254,7 @@ bool _ViscousBuilder::shrink()
         dumpChangeNodes( f );
 
     // Replace source nodes by target nodes in mesh faces to shrink
+    dumpFunction(SMESH_Comment("replNodesOnFace")<<f2sd->first); // debug
     const SMDS_MeshNode* nodes[20];
     for ( size_t i = 0; i < lEdges.size(); ++i )
     {
@@ -4758,6 +5274,7 @@ bool _ViscousBuilder::shrink()
           nodes[iN] = ( n == srcNode ? tgtNode : n );
         }
         helper.GetMeshDS()->ChangeElementNodes( f, nodes, f->NbNodes() );
+        dumpChangeNodes( f );
       }
     }
 
@@ -4767,6 +5284,7 @@ bool _ViscousBuilder::shrink()
     // Create _SmoothNode's on face F
     vector< _SmoothNode > nodesToSmooth( smoothNodes.size() );
     {
+      dumpFunction(SMESH_Comment("fixUVOnFace")<<f2sd->first); // debug
       const bool sortSimplices = isConcaveFace;
       for ( size_t i = 0; i < smoothNodes.size(); ++i )
       {
@@ -4907,6 +5425,17 @@ bool _ViscousBuilder::shrink()
           }
         }
       }
+      // TODO: check effect of this additional smooth
+      // additional laplacian smooth to increase allowed shrink step
+      // for ( int st = 1; st; --st )
+      // {
+      //   dumpFunction(SMESH_Comment("shrinkFace")<<f2sd->first<<"_st"<<++smooStep); // debug
+      //   for ( size_t i = 0; i < nodesToSmooth.size(); ++i )
+      //   {
+      //     nodesToSmooth[i].Smooth( badNb,surface,helper,refSign,
+      //                              _SmoothNode::LAPLACIAN,/*set3D=*/false);
+      //   }
+      // }
     } // while ( shrinked )
 
     // No wrongly shaped faces remain; final smooth. Set node XYZ.
@@ -4971,16 +5500,14 @@ bool _ViscousBuilder::prepareEdgeToShrink( _LayerEdge&            edge,
   const SMDS_MeshNode* srcNode = edge._nodes[0];
   const SMDS_MeshNode* tgtNode = edge._nodes.back();
 
-  edge._pos.clear();
-
   if ( edge._sWOL.ShapeType() == TopAbs_FACE )
   {
-    gp_XY srcUV = helper.GetNodeUV( F, srcNode );
-    gp_XY tgtUV = helper.GetNodeUV( F, tgtNode );
+    gp_XY srcUV( edge._pos[0].X(), edge._pos[0].Y() );//helper.GetNodeUV( F, srcNode );
+    gp_XY tgtUV = edge.LastUV( F );                   //helper.GetNodeUV( F, tgtNode );
     gp_Vec2d uvDir( srcUV, tgtUV );
     double uvLen = uvDir.Magnitude();
     uvDir /= uvLen;
-    edge._normal.SetCoord( uvDir.X(),uvDir.Y(), 0);
+    edge._normal.SetCoord( uvDir.X(),uvDir.Y(), 0 );
     edge._len = uvLen;
 
     edge._pos.resize(1);
@@ -4993,7 +5520,7 @@ bool _ViscousBuilder::prepareEdgeToShrink( _LayerEdge&            edge,
   }
   else // _sWOL is TopAbs_EDGE
   {
-    TopoDS_Edge E = TopoDS::Edge( edge._sWOL);
+    const TopoDS_Edge&    E = TopoDS::Edge( edge._sWOL );
     SMESHDS_SubMesh* edgeSM = getMeshDS()->MeshElements( E );
     if ( !edgeSM || edgeSM->NbElements() == 0 )
       return error(SMESH_Comment("Not meshed EDGE ") << getMeshDS()->ShapeToIndex( E ));
@@ -5012,7 +5539,9 @@ 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();
 
     if ( fabs( uSrc-uTgt ) < 0.99 * fabs( uSrc-u2 ))
     {
@@ -5027,13 +5556,56 @@ bool _ViscousBuilder::prepareEdgeToShrink( _LayerEdge&            edge,
     edge._simplices.resize( 1 );
     edge._simplices[0]._nPrev = n2;
 
-    // set UV of source node to target node
+    // set U of source node to the target node
     SMDS_EdgePosition* pos = static_cast<SMDS_EdgePosition*>( tgtNode->GetPosition() );
     pos->SetUParameter( uSrc );
   }
   return true;
 }
 
+//================================================================================
+/*!
+ * \brief Restore position of a sole node of a _LayerEdge based on _noShrinkShapes
+ */
+//================================================================================
+
+void _ViscousBuilder::restoreNoShrink( _LayerEdge& edge ) const
+{
+  if ( edge._nodes.size() == 1 )
+  {
+    edge._pos.clear();
+    edge._len = 0;
+
+    const SMDS_MeshNode* srcNode = edge._nodes[0];
+    TopoDS_Shape S = SMESH_MesherHelper::GetSubShapeByNode( srcNode, getMeshDS() );
+    if ( S.IsNull() ) return;
+
+    gp_Pnt p;
+
+    switch ( S.ShapeType() )
+    {
+    case TopAbs_EDGE:
+    {
+      double f,l;
+      TopLoc_Location loc;
+      Handle(Geom_Curve) curve = BRep_Tool::Curve( TopoDS::Edge( S ), loc, f, l );
+      if ( curve.IsNull() ) return;
+      SMDS_EdgePosition* ePos = static_cast<SMDS_EdgePosition*>( srcNode->GetPosition() );
+      p = curve->Value( ePos->GetUParameter() );
+      break;
+    }
+    case TopAbs_VERTEX:
+    {
+      p = BRep_Tool::Pnt( TopoDS::Vertex( S ));
+      break;
+    }
+    default: return;
+    }
+    getMeshDS()->MoveNode( srcNode, p.X(), p.Y(), p.Z() );
+    dumpMove( srcNode );
+  }
+}
+
 //================================================================================
 /*!
  * \brief Try to fix triangles with high aspect ratio by swaping diagonals
@@ -5210,7 +5782,7 @@ bool _LayerEdge::SetNewLength2d( Handle(Geom_Surface)& surface,
     const double kSafe = Max( 0.5, 1. - 0.1 * _simplices.size() );
 
     // Select shrinking step such that not to make faces with wrong orientation.
-    double stepSize = uvLen;
+    double stepSize = 1e100;
     for ( size_t i = 0; i < _simplices.size(); ++i )
     {
       // find intersection of 2 lines: curUV-tgtUV and that connecting simplex nodes
@@ -5225,7 +5797,7 @@ bool _LayerEdge::SetNewLength2d( Handle(Geom_Surface)& surface,
         stepSize = Min( step, stepSize );
     }
     gp_Pnt2d newUV;
-    if ( uvLen - stepSize < _len / 200. )
+    if ( uvLen <= stepSize )
     {
       newUV = tgtUV;
       _pos.clear();
@@ -5250,11 +5822,11 @@ bool _LayerEdge::SetNewLength2d( Handle(Geom_Surface)& surface,
   }
   else // _sWOL is TopAbs_EDGE
   {
-    TopoDS_Edge E = TopoDS::Edge( _sWOL );
-    const SMDS_MeshNode* n2 = _simplices[0]._nPrev;
+    const TopoDS_Edge&      E = TopoDS::Edge( _sWOL );
+    const SMDS_MeshNode*   n2 = _simplices[0]._nPrev;
     SMDS_EdgePosition* tgtPos = static_cast<SMDS_EdgePosition*>( tgtNode->GetPosition() );
 
-    const double u2 = helper.GetNodeU( E, n2, tgtNode );
+    const double u2     = helper.GetNodeU( E, n2, tgtNode );
     const double uSrc   = _pos[0].Coord( U_SRC );
     const double lenTgt = _pos[0].Coord( LEN_TGT );
 
@@ -5661,6 +6233,8 @@ bool _ViscousBuilder::addBoundaryElements()
 {
   SMESH_MesherHelper helper( *_mesh );
 
+  vector< const SMDS_MeshNode* > faceNodes;
+
   for ( size_t i = 0; i < _sdVec.size(); ++i )
   {
     _SolidData& data = _sdVec[i];
@@ -5669,6 +6243,8 @@ bool _ViscousBuilder::addBoundaryElements()
     for ( int iE = 1; iE <= geomEdges.Extent(); ++iE )
     {
       const TopoDS_Edge& E = TopoDS::Edge( geomEdges(iE));
+      if ( data._noShrinkShapes.count( getMeshDS()->ShapeToIndex( E )))
+        continue;
 
       // Get _LayerEdge's based on E
 
@@ -5701,8 +6277,13 @@ bool _ViscousBuilder::addBoundaryElements()
         if ( nbSharedPyram > 1 )
           continue; // not free border of the pyramid
 
-        if ( getMeshDS()->FindFace( ledges[0]->_nodes[0], ledges[0]->_nodes[1],
-                                    ledges[1]->_nodes[0], ledges[1]->_nodes[1]))
+        faceNodes.clear();
+        faceNodes.push_back( ledges[0]->_nodes[0] );
+        faceNodes.push_back( ledges[1]->_nodes[0] );
+        if ( ledges[0]->_nodes.size() > 1 ) faceNodes.push_back( ledges[0]->_nodes[1] );
+        if ( ledges[1]->_nodes.size() > 1 ) faceNodes.push_back( ledges[1]->_nodes[1] );
+
+        if ( getMeshDS()->FindElement( faceNodes, SMDSAbs_Face, /*noMedium=*/true))
           continue; // faces already created
       }
       for ( ++u2n; u2n != u2nodes.end(); ++u2n )
@@ -5752,12 +6333,33 @@ bool _ViscousBuilder::addBoundaryElements()
       {
         vector< const SMDS_MeshNode*>&  nn1 = ledges[j-dj1]->_nodes;
         vector< const SMDS_MeshNode*>&  nn2 = ledges[j-dj2]->_nodes;
-        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] ));
+        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] ));
+          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] ));
+        }
+        else if ( nn1.size() == 1 )
+        {
+          if ( isOnFace )
+            for ( size_t z = 1; z < nn2.size(); ++z )
+              sm->AddElement( getMeshDS()->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] ));
+        }
         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]));
+        {
+          if ( isOnFace )
+            for ( size_t z = 1; z < nn1.size(); ++z )
+              sm->AddElement( getMeshDS()->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] ));
+        }
       }
 
       // Make edges
@@ -5767,7 +6369,7 @@ bool _ViscousBuilder::addBoundaryElements()
         if ( !edge->_sWOL.IsNull() && edge->_sWOL.ShapeType() == TopAbs_EDGE )
         {
           vector< const SMDS_MeshNode*>&  nn = edge->_nodes;
-          if ( nn[1]->GetInverseElementIterator( SMDSAbs_Edge )->more() )
+          if ( nn.size() < 2 || nn[1]->GetInverseElementIterator( SMDSAbs_Edge )->more() )
             continue;
           helper.SetSubShape( edge->_sWOL );
           helper.SetElementsOnShape( true );