Salome HOME
decoupage insereFissureGenerale, renommée construitFissureGenerale, premiere phase...
[modules/smesh.git] / src / StdMeshers / StdMeshers_ViscousLayers.cxx
index fa4b13639e741badaa1804b465692c73f198f060..db5469ea228a823a843cc75993ed49714f634883 100644 (file)
@@ -83,7 +83,7 @@
 #include <cmath>
 #include <limits>
 
-//#define __myDEBUG
+#define __myDEBUG
 
 using namespace std;
 
@@ -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.
@@ -349,6 +355,7 @@ namespace VISCOUS_3D
     bool   IsOnEdge() const { return _2neibors; }
     gp_XYZ Copy( _LayerEdge& other, SMESH_MesherHelper& helper );
     void   SetCosin( double cosin );
+    int    NbSteps() const { return _pos.size() - 1; } // nb inlation steps
   };
   struct _LayerEdgeCmp
   {
@@ -408,6 +415,36 @@ namespace VISCOUS_3D
     bool CheckPrisms() const;
   };
 
+  //--------------------------------------------------------------------------------
+  /*!
+   * \brief Layers parameters got by averaging several hypotheses
+   */
+  struct AverageHyp
+  {
+    AverageHyp( const StdMeshers_ViscousLayers* hyp = 0 )
+      :_nbLayers(0), _nbHyps(0), _thickness(0), _stretchFactor(0)
+    {
+      Add( hyp );
+    }
+    void Add( const StdMeshers_ViscousLayers* hyp )
+    {
+      if ( hyp )
+      {
+        _nbHyps++;
+        _nbLayers       = hyp->GetNumberLayers();
+        //_thickness     += hyp->GetTotalThickness();
+        _thickness     = Max( _thickness, hyp->GetTotalThickness() );
+        _stretchFactor += hyp->GetStretchFactor();
+      }
+    }
+    double GetTotalThickness() const { return _thickness; /*_nbHyps ? _thickness / _nbHyps : 0;*/ }
+    double GetStretchFactor()  const { return _nbHyps ? _stretchFactor / _nbHyps : 0; }
+    int    GetNumberLayers()   const { return _nbLayers; }
+  private:
+    int     _nbLayers, _nbHyps;
+    double  _thickness, _stretchFactor;
+  };
+
   //--------------------------------------------------------------------------------
 
   typedef map< const SMDS_MeshNode*, _LayerEdge*, TIDCompare > TNode2Edge;
@@ -418,14 +455,17 @@ namespace VISCOUS_3D
    */
   struct _SolidData
   {
+    typedef const StdMeshers_ViscousLayers* THyp;
     TopoDS_Shape                    _solid;
-    const StdMeshers_ViscousLayers* _hyp;
-    TopoDS_Shape                    _hypShape;
+    TGeomID                         _index; // SOLID id
     _MeshOfSolid*                   _proxyMesh;
-    set<TGeomID>                    _reversedFaceIds;
-    set<TGeomID>                    _ignoreFaceIds; // WOL FACEs and FACEs of other SOLIDS
+    list< THyp >                    _hyps;
+    list< TopoDS_Shape >            _hypShapes;
+    map< TGeomID, THyp >            _face2hyp; // filled if _hyps.size() > 1
+    set< TGeomID >                  _reversedFaceIds;
+    set< TGeomID >                  _ignoreFaceIds; // WOL FACEs and FACEs of other SOLIDs
 
-    double                          _stepSize, _stepSizeCoeff;
+    double                          _stepSize, _stepSizeCoeff, _geomSize;
     const SMDS_MeshNode*            _stepSizeNodes[2];
 
     TNode2Edge                      _n2eMap; // nodes and _LayerEdge's based on them
@@ -433,7 +473,7 @@ namespace VISCOUS_3D
     // 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
-    // iteration over the map is 5 time longer than over the vector
+    // iteration over the map is 5 times longer than over the vector
     vector< _LayerEdge* >           _edges;
 
     // key:   an id of shape (EDGE or VERTEX) shared by a FACE with
@@ -456,15 +496,16 @@ namespace VISCOUS_3D
     vector< int >                    _endEdgeOnShape;
     int                              _nbShapesToSmooth;
 
-    double                           _epsilon; // precision for SegTriaInter()
+    // data of averaged StdMeshers_ViscousLayers parameters for each shape with _LayerEdge's
+    vector< AverageHyp >             _hypOnShape;
+    double                           _maxThickness; // of all _hyps
+    double                           _minThickness; // of all _hyps
 
-    TGeomID                          _index; // SOLID id, for debug
+    double                           _epsilon; // precision for SegTriaInter()
 
-    _SolidData(const TopoDS_Shape&             s=TopoDS_Shape(),
-               const StdMeshers_ViscousLayers* h=0,
-               const TopoDS_Shape&             hs=TopoDS_Shape(),
-               _MeshOfSolid*                   m=0)
-      :_solid(s), _hyp(h), _hypShape(hs), _proxyMesh(m) {}
+    _SolidData(const TopoDS_Shape& s=TopoDS_Shape(),
+               _MeshOfSolid*       m=0)
+      :_solid(s), _proxyMesh(m) {}
     ~_SolidData();
 
     Handle(Geom_Curve) CurveForSmooth( const TopoDS_Edge&    E,
@@ -490,7 +531,7 @@ namespace VISCOUS_3D
       iEnd = _endEdgeOnShape[ end ];
     }
 
-    bool GetShapeEdges(const TGeomID shapeID, size_t& edgeEnd, int* iBeg=0, int* iEnd=0 ) const;
+    bool GetShapeEdges(const TGeomID shapeID, size_t& iEdgeEnd, int* iBeg=0, int* iEnd=0 ) const;
 
     void AddShapesToSmooth( const set< TGeomID >& shapeIDs );
   };
@@ -557,6 +598,9 @@ namespace VISCOUS_3D
     // does it's job
     SMESH_ComputeErrorPtr Compute(SMESH_Mesh&         mesh,
                                   const TopoDS_Shape& shape);
+    // check validity of hypotheses
+    SMESH_ComputeErrorPtr CheckHypotheses( SMESH_Mesh&         mesh,
+                                           const TopoDS_Shape& shape );
 
     // restore event listeners used to clear an inferior dim sub-mesh modified by viscous layers
     void RestoreListeners();
@@ -567,7 +611,11 @@ namespace VISCOUS_3D
   private:
 
     bool findSolidsWithLayers();
-    bool findFacesWithLayers();
+    bool findFacesWithLayers(const bool onlyWith=false);
+    void getIgnoreFaces(const TopoDS_Shape&             solid,
+                        const StdMeshers_ViscousLayers* hyp,
+                        const TopoDS_Shape&             hypShape,
+                        set<TGeomID>&                   ignoreFaces);
     bool makeLayer(_SolidData& data);
     bool setEdgeData(_LayerEdge& edge, const set<TGeomID>& subIds,
                      SMESH_MesherHelper& helper, _SolidData& data);
@@ -589,6 +637,7 @@ 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 );
@@ -724,6 +773,7 @@ namespace VISCOUS_3D
       return _surface->Value( uv.X(), uv.Y() ).XYZ();
     }
   };
+
 } // namespace VISCOUS_3D
 
 
@@ -824,6 +874,27 @@ bool StdMeshers_ViscousLayers::SetParametersByMesh(const SMESH_Mesh*   theMesh,
 {
   // TODO
   return false;
+} // --------------------------------------------------------------------------------
+SMESH_ComputeErrorPtr
+StdMeshers_ViscousLayers::CheckHypothesis(SMESH_Mesh&                          theMesh,
+                                          const TopoDS_Shape&                  theShape,
+                                          SMESH_Hypothesis::Hypothesis_Status& theStatus)
+{
+  VISCOUS_3D::_ViscousBuilder bulder;
+  SMESH_ComputeErrorPtr err = bulder.CheckHypotheses( theMesh, theShape );
+  if ( err && !err->IsOK() )
+    theStatus = SMESH_Hypothesis::HYP_INCOMPAT_HYPS;
+  else
+    theStatus = SMESH_Hypothesis::HYP_OK;
+
+  return err;
+}
+// --------------------------------------------------------------------------------
+bool StdMeshers_ViscousLayers::IsShapeWithLayers(int shapeIndex) const
+{
+  bool isIn =
+    ( std::find( _shapeIds.begin(), _shapeIds.end(), shapeIndex ) != _shapeIds.end() );
+  return IsToIgnoreShapes() ? !isIn : isIn;
 }
 // END StdMeshers_ViscousLayers hypothesis
 //================================================================================
@@ -849,6 +920,7 @@ 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();
@@ -1028,6 +1100,59 @@ namespace
     }
     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
@@ -1251,6 +1376,34 @@ SMESH_ComputeErrorPtr _ViscousBuilder::Compute(SMESH_Mesh&         theMesh,
   return _error;
 }
 
+//================================================================================
+/*!
+ * \brief Check validity of hypotheses
+ */
+//================================================================================
+
+SMESH_ComputeErrorPtr _ViscousBuilder::CheckHypotheses( SMESH_Mesh&         mesh,
+                                                        const TopoDS_Shape& shape )
+{
+  _mesh = & mesh;
+
+  if ( _ViscousListener::GetSolidMesh( _mesh, shape, /*toCreate=*/false))
+    return SMESH_ComputeErrorPtr(); // everything already computed
+
+
+  findSolidsWithLayers();
+  bool ok = findFacesWithLayers();
+
+  // remove _MeshOfSolid's of _SolidData's
+  for ( size_t i = 0; i < _sdVec.size(); ++i )
+    _ViscousListener::RemoveSolidMesh( _mesh, _sdVec[i]._solid );
+
+  if ( !ok )
+    return _error;
+
+  return SMESH_ComputeErrorPtr();
+}
+
 //================================================================================
 /*!
  * \brief Finds SOLIDs to compute using viscous layers. Fills _sdVec
@@ -1274,22 +1427,28 @@ bool _ViscousBuilder::findSolidsWithLayers()
     // TODO: check if algo is hidden
     const list <const SMESHDS_Hypothesis *> & allHyps =
       algo->GetUsedHypothesis(*_mesh, allSolids(i), /*ignoreAuxiliary=*/false);
+    _SolidData* soData = 0;
     list< const SMESHDS_Hypothesis *>::const_iterator hyp = allHyps.begin();
     const StdMeshers_ViscousLayers* viscHyp = 0;
-    for ( ; hyp != allHyps.end() && !viscHyp; ++hyp )
-      viscHyp = dynamic_cast<const StdMeshers_ViscousLayers*>( *hyp );
-    if ( viscHyp )
-    {
-      TopoDS_Shape hypShape;
-      filter.Init( filter.Is( viscHyp ));
-      _mesh->GetHypothesis( allSolids(i), filter, true, &hypShape );
+    for ( ; hyp != allHyps.end(); ++hyp )
+      if ( viscHyp = dynamic_cast<const StdMeshers_ViscousLayers*>( *hyp ))
+      {
+        TopoDS_Shape hypShape;
+        filter.Init( filter.Is( viscHyp ));
+        _mesh->GetHypothesis( allSolids(i), filter, true, &hypShape );
 
-      _MeshOfSolid* proxyMesh = _ViscousListener::GetSolidMesh( _mesh,
-                                                                allSolids(i),
-                                                                /*toCreate=*/true);
-      _sdVec.push_back( _SolidData( allSolids(i), viscHyp, hypShape, proxyMesh ));
-      _sdVec.back()._index = getMeshDS()->ShapeToIndex( allSolids(i));
-    }
+        if ( !soData )
+        {
+          _MeshOfSolid* proxyMesh = _ViscousListener::GetSolidMesh( _mesh,
+                                                                    allSolids(i),
+                                                                    /*toCreate=*/true);
+          _sdVec.push_back( _SolidData( allSolids(i), proxyMesh ));
+          soData = & _sdVec.back();
+          soData->_index = getMeshDS()->ShapeToIndex( allSolids(i));
+        }
+        soData->_hyps.push_back( viscHyp );
+        soData->_hypShapes.push_back( hypShape );
+      }
   }
   if ( _sdVec.empty() )
     return error
@@ -1304,7 +1463,7 @@ bool _ViscousBuilder::findSolidsWithLayers()
  */
 //================================================================================
 
-bool _ViscousBuilder::findFacesWithLayers()
+bool _ViscousBuilder::findFacesWithLayers(const bool onlyWith)
 {
   SMESH_MesherHelper helper( *_mesh );
   TopExp_Explorer exp;
@@ -1315,56 +1474,99 @@ bool _ViscousBuilder::findFacesWithLayers()
   {
     solids.Add( _sdVec[i]._solid );
 
-    vector<TGeomID> ids = _sdVec[i]._hyp->GetBndShapes();
-    if ( _sdVec[i]._hyp->IsToIgnoreShapes() ) // FACEs to ignore are given
+    // get faces to ignore defined by each hyp
+    typedef const StdMeshers_ViscousLayers* THyp;
+    typedef std::pair< set<TGeomID>, THyp > TFacesOfHyp;
+    list< TFacesOfHyp > ignoreFacesOfHyps;
+    list< THyp >::iterator              hyp = _sdVec[i]._hyps.begin();
+    list< TopoDS_Shape >::iterator hypShape = _sdVec[i]._hypShapes.begin();
+    for ( ; hyp != _sdVec[i]._hyps.end(); ++hyp, ++hypShape )
     {
-      for ( size_t ii = 0; ii < ids.size(); ++ii )
-      {
-        const TopoDS_Shape& s = getMeshDS()->IndexToShape( ids[ii] );
-        if ( !s.IsNull() && s.ShapeType() == TopAbs_FACE )
-          _sdVec[i]._ignoreFaceIds.insert( ids[ii] );
-      }
+      ignoreFacesOfHyps.push_back( TFacesOfHyp( set<TGeomID>(), *hyp ));
+      getIgnoreFaces( _sdVec[i]._solid, *hyp, *hypShape, ignoreFacesOfHyps.back().first );
     }
-    else // FACEs with layers are given
+
+    // fill _SolidData::_face2hyp and check compatibility of hypotheses
+    const int nbHyps = _sdVec[i]._hyps.size();
+    if ( nbHyps > 1 )
     {
-      exp.Init( _sdVec[i]._solid, TopAbs_FACE );
-      for ( ; exp.More(); exp.Next() )
+      // check if two hypotheses define different parameters for the same FACE
+      list< TFacesOfHyp >::iterator igFacesOfHyp;
+      for ( exp.Init( _sdVec[i]._solid, TopAbs_FACE ); exp.More(); exp.Next() )
       {
-        TGeomID faceInd = getMeshDS()->ShapeToIndex( exp.Current() );
-        if ( find( ids.begin(), ids.end(), faceInd ) == ids.end() )
-          _sdVec[i]._ignoreFaceIds.insert( faceInd );
+        const TGeomID faceID = getMeshDS()->ShapeToIndex( exp.Current() );
+        THyp hyp = 0;
+        igFacesOfHyp = ignoreFacesOfHyps.begin();
+        for ( ; igFacesOfHyp != ignoreFacesOfHyps.end(); ++igFacesOfHyp )
+          if ( ! igFacesOfHyp->first.count( faceID ))
+          {
+            if ( hyp )
+              return error(SMESH_Comment("Several hypotheses define "
+                                         "Viscous Layers on the face #") << faceID );
+            hyp = igFacesOfHyp->second;
+          }
+        if ( hyp )
+          _sdVec[i]._face2hyp.insert( make_pair( faceID, hyp ));
+        else
+          _sdVec[i]._ignoreFaceIds.insert( faceID );
       }
-    }
 
-    // ignore internal FACEs if inlets and outlets are specified
+      // check if two hypotheses define different number of viscous layers for
+      // adjacent faces of a solid
+      set< int > nbLayersSet;
+      igFacesOfHyp = ignoreFacesOfHyps.begin();
+      for ( ; igFacesOfHyp != ignoreFacesOfHyps.end(); ++igFacesOfHyp )
+      {
+        nbLayersSet.insert( igFacesOfHyp->second->GetNumberLayers() );
+      }
+      if ( nbLayersSet.size() > 1 )
+      {
+        for ( exp.Init( _sdVec[i]._solid, TopAbs_EDGE ); exp.More(); exp.Next() )
+        {
+          PShapeIteratorPtr fIt = helper.GetAncestors( exp.Current(), *_mesh, TopAbs_FACE );
+          THyp hyp1 = 0, hyp2 = 0;
+          while( const TopoDS_Shape* face = fIt->next() )
+          {
+            const TGeomID faceID = getMeshDS()->ShapeToIndex( *face );
+            map< TGeomID, THyp >::iterator f2h = _sdVec[i]._face2hyp.find( faceID );
+            if ( f2h != _sdVec[i]._face2hyp.end() )
+            {
+              ( hyp1 ? hyp2 : hyp1 ) = f2h->second;
+            }
+          }
+          if ( hyp1 && hyp2 &&
+               hyp1->GetNumberLayers() != hyp2->GetNumberLayers() )
+          {
+            return error("Two hypotheses define different number of "
+                         "viscous layers on adjacent faces");
+          }
+        }
+      }
+    } // if ( nbHyps > 1 )
+    else
     {
-      TopTools_IndexedDataMapOfShapeListOfShape solidsOfFace;
-      if ( _sdVec[i]._hyp->IsToIgnoreShapes() )
-        TopExp::MapShapesAndAncestors( _sdVec[i]._hypShape,
-                                       TopAbs_FACE, TopAbs_SOLID, solidsOfFace);
+      _sdVec[i]._ignoreFaceIds.swap( ignoreFacesOfHyps.back().first );
+    }
 
+    // fill _SolidData::_reversedFaceIds
+    {
       exp.Init( _sdVec[i]._solid.Oriented( TopAbs_FORWARD ), TopAbs_FACE );
       for ( ; exp.More(); exp.Next() )
       {
         const TopoDS_Face& face = TopoDS::Face( exp.Current() );
-        if ( helper.NbAncestors( face, *_mesh, TopAbs_SOLID ) < 2 )
-          continue;
-
-        const TGeomID faceInd = getMeshDS()->ShapeToIndex( face );
-        if ( _sdVec[i]._hyp->IsToIgnoreShapes() )
-        {
-          int nbSolids = solidsOfFace.FindFromKey( face ).Extent();
-          if ( nbSolids > 1 )
-            _sdVec[i]._ignoreFaceIds.insert( faceInd );
-        }
-
-        if ( helper.IsReversedSubMesh( face ))
+        const TGeomID faceID = getMeshDS()->ShapeToIndex( face );
+        if ( //!sdVec[i]._ignoreFaceIds.count( faceID ) && ???????
+             helper.NbAncestors( face, *_mesh, TopAbs_SOLID ) > 1 &&
+             helper.IsReversedSubMesh( face ))
         {
-          _sdVec[i]._reversedFaceIds.insert( faceInd );
+          _sdVec[i]._reversedFaceIds.insert( faceID );
         }
       }
     }
-  }
+  } // loop on _sdVec
+
+  if ( onlyWith ) // is called to check hypotheses compatibility only
+    return true;
 
   // Find faces to shrink mesh on (solution 2 in issue 0020832);
   TopTools_IndexedMapOfShape shapes;
@@ -1580,6 +1782,60 @@ bool _ViscousBuilder::findFacesWithLayers()
   return true;
 }
 
+//================================================================================
+/*!
+ * \brief Finds FACEs w/o layers for a given SOLID by an hypothesis
+ */
+//================================================================================
+
+void _ViscousBuilder::getIgnoreFaces(const TopoDS_Shape&             solid,
+                                     const StdMeshers_ViscousLayers* hyp,
+                                     const TopoDS_Shape&             hypShape,
+                                     set<TGeomID>&                   ignoreFaceIds)
+{
+  TopExp_Explorer exp;
+
+  vector<TGeomID> ids = hyp->GetBndShapes();
+  if ( hyp->IsToIgnoreShapes() ) // FACEs to ignore are given
+  {
+    for ( size_t ii = 0; ii < ids.size(); ++ii )
+    {
+      const TopoDS_Shape& s = getMeshDS()->IndexToShape( ids[ii] );
+      if ( !s.IsNull() && s.ShapeType() == TopAbs_FACE )
+        ignoreFaceIds.insert( ids[ii] );
+    }
+  }
+  else // FACEs with layers are given
+  {
+    exp.Init( solid, TopAbs_FACE );
+    for ( ; exp.More(); exp.Next() )
+    {
+      TGeomID faceInd = getMeshDS()->ShapeToIndex( exp.Current() );
+      if ( find( ids.begin(), ids.end(), faceInd ) == ids.end() )
+        ignoreFaceIds.insert( faceInd );
+    }
+  }
+
+  // ignore internal FACEs if inlets and outlets are specified
+  if ( hyp->IsToIgnoreShapes() )
+  {
+    TopTools_IndexedDataMapOfShapeListOfShape solidsOfFace;
+    TopExp::MapShapesAndAncestors( hypShape,
+                                   TopAbs_FACE, TopAbs_SOLID, solidsOfFace);
+
+    for ( exp.Init( solid, TopAbs_FACE ); exp.More(); exp.Next() )
+    {
+      const TopoDS_Face& face = TopoDS::Face( exp.Current() );
+      if ( SMESH_MesherHelper::NbAncestors( face, *_mesh, TopAbs_SOLID ) < 2 )
+        continue;
+
+      int nbSolids = solidsOfFace.FindFromKey( face ).Extent();
+      if ( nbSolids > 1 )
+        ignoreFaceIds.insert( getMeshDS()->ShapeToIndex( face ));
+    }
+  }
+}
+
 //================================================================================
 /*!
  * \brief Create the inner surface of the viscous layer and prepare data for infation
@@ -1652,22 +1908,43 @@ 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;
-      _LayerEdge* maxCosinEdge = 0;
-      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();
-          const bool noShrink = data._noShrinkShapes.count( shapeID );
+          n2e->second = edge;
           edgesByGeom[ shapeID ].push_back( edge );
+          const bool noShrink = data._noShrinkShapes.count( shapeID );
 
           SMESH_TNodeXYZ xyz( n );
 
@@ -1702,7 +1979,13 @@ bool _ViscousBuilder::makeLayer(_SolidData& data)
           }
         }
         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() );
@@ -1711,6 +1994,7 @@ bool _ViscousBuilder::makeLayer(_SolidData& data)
       // compute inflation step size by min size of element on a convex surface
       if ( faceMaxCosin > theMinSmoothCosin )
         limitStepSize( data, face, maxCosinEdge );
+
     } // loop on 2D elements on a FACE
   } // loop on FACEs of a SOLID
 
@@ -1730,16 +2014,17 @@ bool _ViscousBuilder::makeLayer(_SolidData& data)
   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 = ( data._edges[i]->_2neibors->_edges[0] );
+      bool hasData = ( edge->_2neibors->_edges[0] );
       if ( hasData ) // _LayerEdge is a copy of another one
       {
-        nn[0] = data._edges[i]->_2neibors->srcNode(0);
-        nn[1] = data._edges[i]->_2neibors->srcNode(1);
+        nn[0] = edge->_2neibors->srcNode(0);
+        nn[1] = edge->_2neibors->srcNode(1);
       }
-      else if ( !findNeiborsOnEdge( data._edges[i], nn[0],nn[1], data ))
+      else if ( !findNeiborsOnEdge( edge, nn[0],nn[1], data ))
       {
         return false;
       }
@@ -1748,18 +2033,37 @@ bool _ViscousBuilder::makeLayer(_SolidData& data)
       {
         if (( n2e = data._n2eMap.find( nn[j] )) == data._n2eMap.end() )
           return error("_LayerEdge not found by src node", data._index);
-        data._edges[i]->_2neibors->_edges[j] = n2e->second;
+        edge->_2neibors->_edges[j] = n2e->second;
       }
       if ( !hasData )
-        data._edges[i]->SetDataByNeighbors( nn[0], nn[1], helper);
+        edge->SetDataByNeighbors( nn[0], nn[1], helper);
     }
 
-    for ( size_t j = 0; j < data._edges[i]->_simplices.size(); ++j )
+    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();
@@ -1830,7 +2134,6 @@ void _ViscousBuilder::limitStepSize( _SolidData& data, const double minSize )
 void _ViscousBuilder::limitStepSizeByCurvature( _SolidData& data )
 {
   const int nbTestPnt = 5; // on a FACE sub-shape
-  const double minCurvature = 0.9 / data._hyp->GetTotalThickness();
 
   BRepLProp_SLProps surfProp( 2, 1e-6 );
   SMESH_MesherHelper helper( *_mesh );
@@ -1865,8 +2168,9 @@ void _ViscousBuilder::limitStepSizeByCurvature( _SolidData& data )
       else
         continue;
       // check concavity and curvature and limit data._stepSize
+      const double minCurvature = 0.9 / data._hypOnShape[ edgesEnd ].GetTotalThickness();
       int nbLEdges = iEnd - iBeg;
-      int iStep     = Max( 1, nbLEdges / nbTestPnt );
+      int iStep    = Max( 1, nbLEdges / nbTestPnt );
       for ( ; iBeg < iEnd; iBeg += iStep )
       {
         gp_XY uv = helper.GetNodeUV( F, data._edges[ iBeg ]->_nodes[0] );
@@ -1957,8 +2261,21 @@ void _ViscousBuilder::limitStepSizeByCurvature( _SolidData& data )
 bool _ViscousBuilder::sortEdges( _SolidData&                    data,
                                  vector< vector<_LayerEdge*> >& edgesByGeom)
 {
+  // define allowed thickness
+  computeGeomSize( data ); // compute data._geomSize
+
+  data._maxThickness = 0;
+  data._minThickness = 1e100;
+  list< const StdMeshers_ViscousLayers* >::iterator hyp = data._hyps.begin();
+  for ( ; hyp != data._hyps.end(); ++hyp )
+  {
+    data._maxThickness = Max( data._maxThickness, (*hyp)->GetTotalThickness() );
+    data._minThickness = Min( data._minThickness, (*hyp)->GetTotalThickness() );
+  }
+  const double tgtThick = /*Min( 0.5 * data._geomSize, */data._maxThickness;
+
   // Find shapes needing smoothing; such a shape has _LayerEdge._normal on it's
-  // boundry inclined at a sharp angle to the shape
+  // boundry inclined to the shape at a sharp angle
 
   list< TGeomID > shapesToSmooth;
   
@@ -1975,31 +2292,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 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 );
-        // }
         gp_Vec eDir = getEdgeDir( TopoDS::Edge( S ), TopoDS::Vertex( vIt.Value() ));
         double angle = eDir.Angle( eV[0]->_normal );
         double cosin = Cos( angle );
-        needSmooth = ( cosin > theMinSmoothCosin );
+        if ( cosin > theMinSmoothCosin )
+        {
+          // 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;
+            }
+          }
+        }
       }
       break;
     }
@@ -2010,25 +2329,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;
     }
@@ -2036,6 +2368,7 @@ bool _ViscousBuilder::sortEdges( _SolidData&                    data,
       continue;
     default:;
     }
+
     if ( needSmooth )
     {
       if ( S.ShapeType() == TopAbs_EDGE ) shapesToSmooth.push_front( iS );
@@ -2070,6 +2403,45 @@ bool _ViscousBuilder::sortEdges( _SolidData&                    data,
     //eVec.clear();
   }
 
+  // compute average StdMeshers_ViscousLayers parameters for each shape
+
+  data._hypOnShape.clear();
+  if ( data._hyps.size() == 1 )
+  {
+    data._hypOnShape.resize( data._endEdgeOnShape.size(), AverageHyp( data._hyps.back() ));
+  }
+  else
+  {
+    data._hypOnShape.resize( data._endEdgeOnShape.size() );
+    map< TGeomID, const StdMeshers_ViscousLayers* >::iterator f2hyp;
+    for ( size_t i = 0; i < data._endEdgeOnShape.size(); ++i )
+    {
+      int       iEnd = data._endEdgeOnShape[i];
+      _LayerEdge* LE = data._edges[ iEnd-1 ];
+      TGeomID iShape = LE->_nodes[0]->getshapeId();
+      const TopoDS_Shape& S = getMeshDS()->IndexToShape( iShape );
+      if ( S.ShapeType() == TopAbs_FACE )
+      {
+        if (( f2hyp = data._face2hyp.find( iShape )) != data._face2hyp.end() )
+        {
+          data._hypOnShape[ i ].Add( f2hyp->second );
+        }
+      }
+      else
+      {
+        PShapeIteratorPtr fIt = SMESH_MesherHelper::GetAncestors( S, *_mesh, TopAbs_FACE );
+        while ( const TopoDS_Shape* face = fIt->next() )
+        {
+          TGeomID faceID = getMeshDS()->ShapeToIndex( *face );
+          if (( f2hyp = data._face2hyp.find( faceID )) != data._face2hyp.end() )
+          {
+            data._hypOnShape[ i ].Add( f2hyp->second );
+          }
+        }
+      }
+    }
+  }
+
   return ok;
 }
 
@@ -2551,11 +2923,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 )
@@ -2569,10 +2941,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;
@@ -2739,6 +3114,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
@@ -2751,46 +3151,40 @@ 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 )
-    limitStepSize( data, tgtThick );
+  const double tgtThick = data._maxThickness;
+  if ( data._stepSize > data._minThickness )
+    limitStepSize( data, data._minThickness );
 
   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;
-  while ( 1.01 * avgThick < tgtThick )
+  int iBeg, iEnd, iS;
+  while ( avgThick < 0.99 )
   {
     // new target length
     curThick += data._stepSize;
     if ( curThick > tgtThick )
     {
-      curThick = tgtThick + ( tgtThick-avgThick ) * nbRepeats;
+      curThick = tgtThick + tgtThick*( 1.-avgThick ) * nbRepeats;
       nbRepeats++;
     }
 
     // Elongate _LayerEdge's
     dumpFunction(SMESH_Comment("inflate")<<data._index<<"_step"<<nbSteps); // debug
-    for ( size_t i = 0; i < data._edges.size(); ++i )
+    for ( iBeg = 0, iS = 0; iS < data._endEdgeOnShape.size(); ++iS )
     {
-      data._edges[i]->SetNewLength( curThick, helper );
+      const double shapeCurThick = Min( curThick, data._hypOnShape[ iS ].GetTotalThickness() );
+      for ( iEnd = data._endEdgeOnShape[ iS ]; iBeg < iEnd; ++iBeg )
+      {
+        data._edges[iBeg]->SetNewLength( shapeCurThick, helper );
+      }
     }
     dumpFunctionEnd();
 
@@ -2815,16 +3209,22 @@ bool _ViscousBuilder::inflate(_SolidData& data)
 
     // Evaluate achieved thickness
     avgThick = 0;
-    for ( size_t i = 0; i < data._edges.size(); ++i )
-      avgThick += data._edges[i]->_len;
+    for ( iBeg = 0, iS = 0; iS < data._endEdgeOnShape.size(); ++iS )
+    {
+      const double shapeTgtThick = data._hypOnShape[ iS ].GetTotalThickness();
+      for ( iEnd = data._endEdgeOnShape[ iS ]; iBeg < iEnd; ++iBeg )
+      {
+        avgThick += Min( 1., data._edges[iBeg]->_len / shapeTgtThick );
+      }
+    }
     avgThick /= data._edges.size();
-    debugMsg( "-- Thickness " << avgThick << " reached" );
+    debugMsg( "-- Thickness " << curThick << " ("<< avgThick*100 << "%) reached" );
 
-    if ( distToIntersection < avgThick*1.5 )
+    if ( distToIntersection < tgtThick*avgThick*1.5 )
     {
       debugMsg( "-- Stop inflation since "
                 << " distToIntersection( "<<distToIntersection<<" ) < avgThick( "
-                << avgThick << " ) * 1.5" );
+                << tgtThick*avgThick << " ) * 1.5" );
       break;
     }
     // new step size
@@ -2833,12 +3233,12 @@ bool _ViscousBuilder::inflate(_SolidData& data)
       data._stepSize = data._stepSizeCoeff *
         SMESH_TNodeXYZ(data._stepSizeNodes[0]).Distance(data._stepSizeNodes[1]);
 
-  } // while ( 1.01 * avgThick < tgtThick )
+  } // while ( avgThick < 0.99 )
 
   if (nbSteps == 0 )
     return error("failed at the very first inflation step", data._index);
 
-  if ( 1.01 * avgThick < tgtThick )
+  if ( avgThick < 0.99 )
     if ( SMESH_subMesh* sm = _mesh->GetSubMeshContaining( data._index ))
     {
       SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
@@ -2847,14 +3247,13 @@ bool _ViscousBuilder::inflate(_SolidData& data)
           ( new SMESH_ComputeError (COMPERR_WARNING,
                                     SMESH_Comment("Thickness ") << tgtThick <<
                                     " of viscous layers not reached,"
-                                    " average reached thickness is " << avgThick ));
+                                    " average reached thickness is " << avgThick*tgtThick));
     }
 
 
   // 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 )
+  for ( iEnd = iS = 0; iS < data._endEdgeOnShape.size(); ++iS )
   {
     iBeg = iEnd;
     iEnd = data._endEdgeOnShape[ iS ];
@@ -2894,6 +3293,16 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data,
     iBeg = iEnd;
     iEnd = data._endEdgeOnShape[ iS ];
 
+    // bool toSmooth = false;
+    // for ( int i = iBeg; i < iEnd; ++i )
+    //   toSmooth = data._edges[ iBeg ]->NbSteps() >= nbSteps+1;
+    // if ( !toSmooth )
+    // {
+    //   if ( iS+1 == data._nbShapesToSmooth )
+    //     data._nbShapesToSmooth--;
+    //   continue; // target length reached some steps before
+    // }
+
     if ( !data._edges[ iBeg ]->_sWOL.IsNull() &&
          data._edges[ iBeg ]->_sWOL.ShapeType() == TopAbs_FACE )
     {
@@ -2942,10 +3351,10 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data,
         badNb = 0;
         moved = false;
         if ( step % 2 )
-          for ( int i = iBeg; i < iEnd; ++i )
+          for ( int i = iBeg; i < iEnd; ++i ) // iterate forward
             moved |= data._edges[i]->Smooth(badNb);
         else
-          for ( int i = iEnd-1; i >= iBeg; --i )
+          for ( int i = iEnd-1; i >= iBeg; --i ) // iterate backward
             moved |= data._edges[i]->Smooth(badNb);
         improved = ( badNb < oldBadNb );
 
@@ -3183,14 +3592,14 @@ void _SolidData::SortOnEdge( const TopoDS_Edge&  E,
 //================================================================================
 
 bool _SolidData::GetShapeEdges(const TGeomID shapeID,
-                               size_t &      edgesEnd,
+                               size_t &      iEdgesEnd,
                                int*          iBeg,
                                int*          iEnd ) const
 {
   int beg = 0, end = 0;
-  for ( edgesEnd = 0; edgesEnd < _endEdgeOnShape.size(); ++edgesEnd )
+  for ( iEdgesEnd = 0; iEdgesEnd < _endEdgeOnShape.size(); ++iEdgesEnd )
   {
-    end = _endEdgeOnShape[ edgesEnd ];
+    end = _endEdgeOnShape[ iEdgesEnd ];
     TGeomID sID = _edges[ beg ]->_nodes[0]->getshapeId();
     if ( sID == shapeID )
     {
@@ -4552,7 +4961,7 @@ void _LayerEdge::SetNewLength( double len, SMESH_MesherHelper& helper )
 {
   if ( _len - len > -1e-6 )
   {
-    _pos.push_back( _pos.back() );
+    //_pos.push_back( _pos.back() );
     return;
   }
 
@@ -4662,6 +5071,8 @@ bool _ViscousBuilder::refine(_SolidData& data)
 
   // Create intermediate nodes on each _LayerEdge
 
+  int iS = 0, iEnd = data._endEdgeOnShape[ iS ];
+
   for ( size_t i = 0; i < data._edges.size(); ++i )
   {
     _LayerEdge& edge = *data._edges[i];
@@ -4669,6 +5080,11 @@ bool _ViscousBuilder::refine(_SolidData& data)
     if ( edge._nodes.size() < 2 )
       continue; // on _noShrinkShapes
 
+    // get parameters of layers for the edge
+    if ( i == iEnd )
+      iEnd = data._endEdgeOnShape[ ++iS ];
+    const AverageHyp& hyp = data._hypOnShape[ iS ];
+
     // get accumulated length of segments
     vector< double > segLen( edge._pos.size() );
     segLen[0] = 0.0;
@@ -4679,7 +5095,7 @@ bool _ViscousBuilder::refine(_SolidData& data)
     const SMDS_MeshNode* tgtNode = edge._nodes.back();
     if ( edge._nodes.size() == 2 )
     {
-      edge._nodes.resize( data._hyp->GetNumberLayers() + 1, 0 );
+      edge._nodes.resize( hyp.GetNumberLayers() + 1, 0 );
       edge._nodes[1] = 0;
       edge._nodes.back() = tgtNode;
     }
@@ -4728,8 +5144,8 @@ bool _ViscousBuilder::refine(_SolidData& data)
     // calculate height of the first layer
     double h0;
     const double T = segLen.back(); //data._hyp.GetTotalThickness();
-    const double f = data._hyp->GetStretchFactor();
-    const int    N = data._hyp->GetNumberLayers();
+    const double f = hyp.GetStretchFactor();
+    const int    N = hyp.GetNumberLayers();
     const double fPowN = pow( f, N );
     if ( fPowN - 1 <= numeric_limits<double>::min() )
       h0 = T / N;
@@ -4841,35 +5257,45 @@ bool _ViscousBuilder::refine(_SolidData& data)
   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() )))
+    const TGeomID faceID = getMeshDS()->ShapeToIndex( exp.Current() );
+    if ( data._ignoreFaceIds.count( faceID ))
       continue;
-    SMESHDS_SubMesh*   fSubM = getMeshDS()->MeshElements( exp.Current() );
-    SMDS_ElemIteratorPtr fIt = fSubM->GetElements();
+    const bool isReversedFace = data._reversedFaceIds.count( faceID );
+    SMESHDS_SubMesh*    fSubM = getMeshDS()->MeshElements( exp.Current() );
+    SMDS_ElemIteratorPtr  fIt = fSubM->GetElements();
     while ( fIt->more() )
     {
       const SMDS_MeshElement* face = fIt->next();
       const int            nbNodes = face->NbCornerNodes();
       nnVec.resize( nbNodes );
+      nnSet.clear();
       degenEdgeInd.clear();
       int nbZ = 0;
-      SMDS_ElemIteratorPtr nIt = face->nodesIterator();
+      SMDS_NodeIteratorPtr nIt = face->nodeIterator();
       for ( int iN = 0; iN < nbNodes; ++iN )
       {
-        const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nIt->next() );
-        nnVec[ iN ] = & data._n2eMap[ n ]->_nodes;
-        if ( nnVec[ iN ]->size() < 2 )
+        const SMDS_MeshNode* n = nIt->next();
+        const int i = isReversedFace ? nbNodes-1-iN : iN;
+        nnVec[ i ] = & data._n2eMap[ n ]->_nodes;
+        if ( nnVec[ i ]->size() < 2 )
           degenEdgeInd.insert( iN );
         else
-          nbZ = nnVec[ iN ]->size();
+          nbZ = nnVec[ i ]->size();
+
+        if ( helper.HasDegeneratedEdges() )
+          nnSet.insert( nnVec[ i ]);
       }
       if ( nbZ == 0 )
         continue;
+      if ( 0 < nnSet.size() && nnSet.size() < 3 )
+        continue;
 
       switch ( nbNodes )
       {
@@ -6079,6 +6505,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];
@@ -6121,8 +6549,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 )