Salome HOME
52459: Viscous layers are not normal to the surface.
[modules/smesh.git] / src / StdMeshers / StdMeshers_ViscousLayers.cxx
index b5a076c5e3edc5ca459963484f9c044197a9d31d..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.
@@ -425,7 +431,7 @@ namespace VISCOUS_3D
     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
@@ -589,6 +595,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 +731,7 @@ namespace VISCOUS_3D
       return _surface->Value( uv.X(), uv.Y() ).XYZ();
     }
   };
+
 } // namespace VISCOUS_3D
 
 
@@ -849,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);
@@ -888,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 );
@@ -1019,6 +1037,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
@@ -1643,22 +1714,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 );
 
@@ -1693,7 +1785,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() );
@@ -1702,6 +1800,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
 
@@ -1721,16 +1820,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;
       }
@@ -1739,18 +1839,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();
@@ -1948,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;
   
@@ -1966,31 +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 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;
     }
@@ -2001,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;
     }
@@ -2027,6 +2165,7 @@ bool _ViscousBuilder::sortEdges( _SolidData&                    data,
       continue;
     default:;
     }
+
     if ( needSmooth )
     {
       if ( S.ShapeType() == TopAbs_EDGE ) shapesToSmooth.push_front( iS );
@@ -2354,9 +2493,19 @@ gp_XYZ _ViscousBuilder::getFaceNormal(const SMDS_MeshNode* node,
   Handle(Geom_Surface) surface = BRep_Tool::Surface( face );
   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
+  {
+    // 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 )
+    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. )
@@ -2364,7 +2513,8 @@ gp_XYZ _ViscousBuilder::getFaceNormal(const SMDS_MeshNode* node,
     }
     isOK = true;
   }
-  else // hard singularity, to call with shiftInside=true ?
+
+  if ( !isOK ) // hard singularity, to call with shiftInside=true ?
   {
     const TGeomID faceID = helper.GetMeshDS()->ShapeToIndex( face );
 
@@ -2377,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;
         }
@@ -2529,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 )
@@ -2547,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;
@@ -2717,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
@@ -2729,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 )
@@ -2750,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;
@@ -4819,6 +4988,7 @@ 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;
 
@@ -4834,20 +5004,26 @@ bool _ViscousBuilder::refine(_SolidData& data)
       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() );
+        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;
 
       switch ( nbNodes )
       {
@@ -6057,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];
@@ -6099,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 )