Salome HOME
http://www.salome-platform.org/forum/forum_10/450300019
[modules/smesh.git] / src / StdMeshers / StdMeshers_ViscousLayers.cxx
index 1f6d46263460186259fe57f14cafee07b51facb4..b2f6855b01c065ec0245b23b6a550c94e1552f25 100644 (file)
@@ -45,6 +45,8 @@
 #include "StdMeshers_FaceSide.hxx"
 
 #include <BRepAdaptor_Curve2d.hxx>
+#include <BRepAdaptor_Surface.hxx>
+#include <BRepLProp_SLProps.hxx>
 #include <BRep_Tool.hxx>
 #include <Bnd_B2d.hxx>
 #include <Bnd_B3d.hxx>
@@ -61,6 +63,7 @@
 #include <Geom_TrimmedCurve.hxx>
 #include <Precision.hxx>
 #include <Standard_ErrorHandler.hxx>
+#include <Standard_Failure.hxx>
 #include <TColStd_Array1OfReal.hxx>
 #include <TopExp.hxx>
 #include <TopExp_Explorer.hxx>
@@ -323,7 +326,7 @@ namespace VISCOUS_3D
 
     gp_XYZ              _normal; // to solid surface
     vector<gp_XYZ>      _pos; // points computed during inflation
-    double              _len; // length achived with the last step
+    double              _len; // length achived with the last inflation step
     double              _cosin; // of angle (_normal ^ surface)
     double              _lenFactor; // to compute _len taking _cosin into account
 
@@ -362,7 +365,7 @@ namespace VISCOUS_3D
                        const double&        epsilon) const;
     gp_Ax1 LastSegment(double& segLen) const;
     bool IsOnEdge() const { return _2neibors; }
-    void Copy( _LayerEdge& other, SMESH_MesherHelper& helper );
+    gp_XYZ Copy( _LayerEdge& other, SMESH_MesherHelper& helper );
     void SetCosin( double cosin );
   };
   struct _LayerEdgeCmp
@@ -394,9 +397,13 @@ namespace VISCOUS_3D
     const SMDS_MeshNode*            _stepSizeNodes[2];
 
     TNode2Edge                      _n2eMap;
+    // 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
     vector< _LayerEdge* >           _edges;
+    // edges on EDGE's with null _sWOL, whose _simplices are used to stop inflation
+    vector< _LayerEdge* >           _simplexTestEdges;
 
     // key:   an id of shape (EDGE or VERTEX) shared by a FACE with
     //        layers and a FACE w/o layers
@@ -479,6 +486,9 @@ namespace VISCOUS_3D
     bool makeLayer(_SolidData& data);
     bool setEdgeData(_LayerEdge& edge, const set<TGeomID>& subIds,
                      SMESH_MesherHelper& helper, _SolidData& data);
+    gp_XYZ getWeigthedNormal( const SMDS_MeshNode*         n,
+                              std::pair< TGeomID, gp_XYZ > fId2Normal[],
+                              const int                    nbFaces );
     bool findNeiborsOnEdge(const _LayerEdge*     edge,
                            const SMDS_MeshNode*& n1,
                            const SMDS_MeshNode*& n2,
@@ -487,6 +497,8 @@ namespace VISCOUS_3D
                        const set<TGeomID>& ingnoreShapes,
                        const _SolidData*   dataToCheckOri = 0,
                        const bool          toSort = false);
+    void findSimplexTestEdges( _SolidData&                    data,
+                               vector< vector<_LayerEdge*> >& edgesByGeom);
     bool sortEdges( _SolidData&                    data,
                     vector< vector<_LayerEdge*> >& edgesByGeom);
     void limitStepSize( _SolidData&             data,
@@ -794,23 +806,29 @@ namespace
     // get average dir of edges going fromV
     gp_XYZ edgeDir;
     //if ( edges.size() > 1 )
-      for ( size_t i = 0; i < edges.size(); ++i )
-      {
-        edgeDir = getEdgeDir( edges[i], fromV );
-        double size2 = edgeDir.SquareModulus();
-        if ( size2 > numeric_limits<double>::min() )
-          edgeDir /= sqrt( size2 );
-        else
-          ok = false;
-        dir += edgeDir;
-      }
+    for ( size_t i = 0; i < edges.size(); ++i )
+    {
+      edgeDir = getEdgeDir( edges[i], fromV );
+      double size2 = edgeDir.SquareModulus();
+      if ( size2 > numeric_limits<double>::min() )
+        edgeDir /= sqrt( size2 );
+      else
+        ok = false;
+      dir += edgeDir;
+    }
     gp_XYZ fromEdgeDir = getFaceDir( F, edges[0], node, helper, ok );
-    if ( edges.size() == 1 )
-      dir = fromEdgeDir;
-    else if ( dir.SquareModulus() < 0.1 ) // ~< 20 degrees
-      dir = fromEdgeDir + getFaceDir( F, edges[1], node, helper, ok );
-    else if ( dir * fromEdgeDir < 0 )
-      dir *= -1;
+    try {
+      if ( edges.size() == 1 )
+        dir = fromEdgeDir;
+      else if ( dir.SquareModulus() < 0.1 ) // ~< 20 degrees
+        dir = fromEdgeDir.Normalized() + getFaceDir( F, edges[1], node, helper, ok ).Normalized();
+      else if ( dir * fromEdgeDir < 0 )
+        dir *= -1;
+    }
+    catch ( Standard_Failure )
+    {
+      ok = false;
+    }
     if ( ok )
     {
       //dir /= edges.size();
@@ -1389,6 +1407,19 @@ bool _ViscousBuilder::findFacesWithLayers()
     }
   }
 
+  // add FACEs of other SOLIDs to _ignoreFaceIds
+  for ( size_t i = 0; i < _sdVec.size(); ++i )
+  {
+    shapes.Clear();
+    TopExp::MapShapes(_sdVec[i]._solid, TopAbs_FACE, shapes);
+
+    for ( exp.Init( _mesh->GetShapeToMesh(), TopAbs_FACE ); exp.More(); exp.Next() )
+    {
+      if ( !shapes.Contains( exp.Current() ))
+        _sdVec[i]._ignoreFaceIds.insert( getMeshDS()->ShapeToIndex( exp.Current() ));
+    }
+  }
+
   return true;
 }
 
@@ -1416,7 +1447,6 @@ bool _ViscousBuilder::makeLayer(_SolidData& data)
     }
 
   // make a map to find new nodes on sub-shapes shared with other SOLID
-  map< TGeomID, TNode2Edge* > s2neMap;
   map< TGeomID, TNode2Edge* >::iterator s2ne;
   map< TGeomID, TopoDS_Shape >::iterator s2s = data._shrinkShape2Shape.begin();
   for (; s2s != data._shrinkShape2Shape.end(); ++s2s )
@@ -1429,7 +1459,7 @@ bool _ViscousBuilder::makeLayer(_SolidData& data)
       if ( s2s2 != _sdVec[i]._shrinkShape2Shape.end() &&
            *s2s == *s2s2 && !_sdVec[i]._n2eMap.empty() )
       {
-        s2neMap.insert( make_pair( shapeInd, &_sdVec[i]._n2eMap ));
+        data._s2neMap.insert( make_pair( shapeInd, &_sdVec[i]._n2eMap ));
         break;
       }
     }
@@ -1481,21 +1511,23 @@ bool _ViscousBuilder::makeLayer(_SolidData& data)
           const int shapeID = n->getshapeId();
           edgesByGeom[ shapeID ].push_back( edge );
 
+          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 = s2neMap.find( shapeID )) != s2neMap.end() &&
+               ( s2ne = data._s2neMap.find( shapeID )) != data._s2neMap.end() &&
                ( n2e2 = (*s2ne).second->find( n )) != s2ne->second->end())
           {
             _LayerEdge* foundEdge = (*n2e2).second;
-            edge->Copy( *foundEdge, helper );
-            // location of the last node is modified but we can restore
-            // it by node position on _sWOL stored by the node
+            gp_XYZ        lastPos = edge->Copy( *foundEdge, helper );
+            foundEdge->_pos.push_back( lastPos );
+            // location of the last node is modified and we restore it by foundEdge->_pos.back()
             const_cast< SMDS_MeshNode* >
-              ( edge->_nodes.back() )->setXYZ( n->X(), n->Y(), n->Z() );
+              ( edge->_nodes.back() )->setXYZ( xyz.X(), xyz.Y(), xyz.Z() );
           }
           else
           {
-            edge->_nodes.push_back( helper.AddNode( n->X(), n->Y(), n->Z() ));
+            edge->_nodes.push_back( helper.AddNode( xyz.X(), xyz.Y(), xyz.Z() ));
             if ( !setEdgeData( *edge, subIds, helper, data ))
               return false;
           }
@@ -1522,12 +1554,14 @@ bool _ViscousBuilder::makeLayer(_SolidData& data)
   if ( data._stepSize < 1. )
     data._epsilon *= data._stepSize;
 
-  // Put _LayerEdge's into the vector data._edges
+  // fill data._simplexTestEdges
+  findSimplexTestEdges( data, edgesByGeom );
 
+  // Put _LayerEdge's into the vector data._edges
   if ( !sortEdges( data, edgesByGeom ))
     return false;
 
-  // Set target nodes into _Simplex and _2NearEdges
+  // Set target nodes into _Simplex and _2NearEdges of _LayerEdge's
   TNode2Edge::iterator n2e;
   for ( size_t i = 0; i < data._edges.size(); ++i )
   {
@@ -1542,13 +1576,13 @@ bool _ViscousBuilder::makeLayer(_SolidData& data)
         n = (*n2e).second->_nodes.back();
         data._edges[i]->_2neibors->_edges[j] = n2e->second;
       }
-    else
-      for ( size_t j = 0; j < data._edges[i]->_simplices.size(); ++j )
-      {
-        _Simplex& s = data._edges[i]->_simplices[j];
-        s._nNext = data._n2eMap[ s._nNext ]->_nodes.back();
-        s._nPrev = data._n2eMap[ s._nPrev ]->_nodes.back();
-      }
+    //else
+    for ( size_t j = 0; j < data._edges[i]->_simplices.size(); ++j )
+    {
+      _Simplex& s = data._edges[i]->_simplices[j];
+      s._nNext = data._n2eMap[ s._nNext ]->_nodes.back();
+      s._nPrev = data._n2eMap[ s._nPrev ]->_nodes.back();
+    }
   }
 
   dumpFunctionEnd();
@@ -1596,7 +1630,7 @@ void _ViscousBuilder::limitStepSize( _SolidData&             data,
  */
 //================================================================================
 
-void _ViscousBuilder::limitStepSize( _SolidData& data, const double minSize)
+void _ViscousBuilder::limitStepSize( _SolidData& data, const double minSize )
 {
   if ( minSize < data._stepSize )
   {
@@ -1610,6 +1644,96 @@ void _ViscousBuilder::limitStepSize( _SolidData& data, const double minSize)
   }
 }
 
+//================================================================================
+/*!
+ * Fill data._simplexTestEdges. These _LayerEdge's are used to stop inflation
+ * in the case where there are no _LayerEdge's on a curved convex FACE,
+ * as e.g. on a fillet surface with no internal nodes - issue 22580,
+ * so that collision of viscous internal faces is not detected by check of
+ * intersection of _LayerEdge's with the viscous internal faces.
+ */
+//================================================================================
+
+void _ViscousBuilder::findSimplexTestEdges( _SolidData&                    data,
+                                            vector< vector<_LayerEdge*> >& edgesByGeom)
+{
+  data._simplexTestEdges.clear();
+
+  SMESH_MesherHelper helper( *_mesh );
+
+  vector< vector<_LayerEdge*> * > ledgesOnEdges;
+  set< const SMDS_MeshNode* > usedNodes;
+
+  const double minCurvature = 1. / data._hyp->GetTotalThickness();
+
+  for ( size_t iS = 0; iS < edgesByGeom.size(); ++iS )
+  {
+    // look for a FACE with layers and w/o _LayerEdge's
+    const vector<_LayerEdge*>& eS = edgesByGeom[iS];
+    if ( !eS.empty() ) continue;
+    const TopoDS_Shape& S = getMeshDS()->IndexToShape( iS );
+    if ( S.IsNull() || S.ShapeType() != TopAbs_FACE ) continue;
+    if ( data._ignoreFaceIds.count( iS )) continue;
+
+    const TopoDS_Face& F = TopoDS::Face( S );
+
+    // look for _LayerEdge's on EDGEs with null _sWOL
+    ledgesOnEdges.clear();
+    TopExp_Explorer eExp( F, TopAbs_EDGE );
+    for ( ; eExp.More(); eExp.Next() )
+    {
+      TGeomID iE = getMeshDS()->ShapeToIndex( eExp.Current() );
+      vector<_LayerEdge*>& eE = edgesByGeom[iE];
+      if ( !eE.empty() && eE[0]->_sWOL.IsNull() )
+        ledgesOnEdges.push_back( & eE );
+    }
+    if ( ledgesOnEdges.empty() ) continue;
+
+    // check FACE convexity
+    const _LayerEdge* le = ledgesOnEdges[0]->back();
+    gp_XY uv = helper.GetNodeUV( F, le->_nodes[0] );
+    BRepAdaptor_Surface surf( F );
+    BRepLProp_SLProps surfProp( surf, uv.X(), uv.Y(), 2, 1e-6 );
+    if ( !surfProp.IsCurvatureDefined() )
+      continue;
+    double surfCurvature = Max( Abs( surfProp.MaxCurvature() ),
+                                Abs( surfProp.MinCurvature() ));
+    if ( surfCurvature < minCurvature )
+      continue;
+    gp_Dir minDir, maxDir;
+    surfProp.CurvatureDirections( maxDir, minDir );
+    if ( F.Orientation() == TopAbs_REVERSED ) {
+      maxDir.Reverse(); minDir.Reverse();
+    }
+    const gp_XYZ& inDir = le->_normal;
+    if ( inDir * maxDir.XYZ() < 0 &&
+         inDir * minDir.XYZ() < 0 )
+      continue;
+
+    limitStepSize( data, 0.9 / surfCurvature );
+
+    // add _simplices to the _LayerEdge's
+    for ( size_t iE = 0; iE < ledgesOnEdges.size(); ++iE )
+    {
+      const vector<_LayerEdge*>& ledges = *ledgesOnEdges[iE];
+      for ( size_t iLE = 0; iLE < ledges.size(); ++iLE )
+      {
+        _LayerEdge* ledge = ledges[iLE];
+        const SMDS_MeshNode* srcNode = ledge->_nodes[0];
+        if ( !usedNodes.insert( srcNode ).second ) continue;
+
+        getSimplices( srcNode, ledge->_simplices, data._ignoreFaceIds, &data );
+        for ( size_t i = 0; i < ledge->_simplices.size(); ++i )
+        {
+          usedNodes.insert( ledge->_simplices[i]._nPrev );
+          usedNodes.insert( ledge->_simplices[i]._nNext );
+        }
+        data._simplexTestEdges.push_back( ledge );
+      }
+    }
+  }
+}
+
 //================================================================================
 /*!
  * \brief Separate shapes (and _LayerEdge's on them) to smooth from the rest ones
@@ -1631,7 +1755,7 @@ bool _ViscousBuilder::sortEdges( _SolidData&                    data,
   {
     vector<_LayerEdge*>& eS = edgesByGeom[iS];
     if ( eS.empty() ) continue;
-    TopoDS_Shape S = getMeshDS()->IndexToShape( iS );
+    const TopoDS_Shape& S = getMeshDS()->IndexToShape( iS );
     bool needSmooth = false;
     switch ( S.ShapeType() )
     {
@@ -1804,12 +1928,12 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge&         edge,
 
     set<TGeomID>::iterator id = faceIds.begin();
     TopoDS_Face F;
+    std::pair< TGeomID, gp_XYZ > id2Norm[20];
     for ( ; id != faceIds.end(); ++id )
     {
       const TopoDS_Shape& s = getMeshDS()->IndexToShape( *id );
       if ( s.IsNull() || s.ShapeType() != TopAbs_FACE || !subIds.count( *id ))
         continue;
-      totalNbFaces++;
       F = TopoDS::Face( s );
 
       gp_XY uv = helper.GetNodeUV( F, node, 0, &normOK );
@@ -1844,12 +1968,22 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge&         edge,
       }
       if ( helper.GetSubShapeOri( data._solid, F ) != TopAbs_REVERSED )
         geomNorm.Reverse();
+      id2Norm[ totalNbFaces ].first  = *id;
+      id2Norm[ totalNbFaces ].second = geomNorm.XYZ();
+      totalNbFaces++;
       edge._normal += geomNorm.XYZ();
     }
     if ( totalNbFaces == 0 )
       return error(SMESH_Comment("Can't get normal to node ") << node->GetID(), data._index);
 
-    edge._normal /= totalNbFaces;
+    if ( totalNbFaces < 3 )
+    {
+      //edge._normal /= totalNbFaces;
+    }
+    else
+    {
+      edge._normal = getWeigthedNormal( node, id2Norm, totalNbFaces );
+    }
 
     switch ( posType )
     {
@@ -1951,6 +2085,87 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge&         edge,
   return true;
 }
 
+//================================================================================
+/*!
+ * \brief Return normal at a node weighted with angles taken by FACEs
+ *  \param [in] n - the node
+ *  \param [in] fId2Normal - FACE ids and normals
+ *  \param [in] nbFaces - nb of FACEs meeting at the node
+ *  \return gp_XYZ - computed normal
+ */
+//================================================================================
+
+gp_XYZ _ViscousBuilder::getWeigthedNormal( const SMDS_MeshNode*         n,
+                                           std::pair< TGeomID, gp_XYZ > fId2Normal[],
+                                           const int                    nbFaces )
+{
+  gp_XYZ resNorm(0,0,0);
+  TopoDS_Shape V = SMESH_MesherHelper::GetSubShapeByNode( n, getMeshDS() );
+  if ( V.ShapeType() != TopAbs_VERTEX )
+  {
+    for ( int i = 0; i < nbFaces; ++i )
+      resNorm += fId2Normal[i].second / nbFaces ;
+    return resNorm;
+  }
+
+  double angles[30];
+  for ( int i = 0; i < nbFaces; ++i )
+  {
+    const TopoDS_Face& F = TopoDS::Face( getMeshDS()->IndexToShape( fId2Normal[i].first ));
+
+    // look for two EDGEs shared by F and other FACEs within fId2Normal
+    TopoDS_Edge ee[2];
+    int nbE = 0;
+    PShapeIteratorPtr eIt = SMESH_MesherHelper::GetAncestors( V, *_mesh, TopAbs_EDGE );
+    while ( const TopoDS_Shape* E = eIt->next() )
+    {
+      if ( !SMESH_MesherHelper::IsSubShape( *E, F ))
+        continue;
+      bool isSharedEdge = false;
+      for ( int j = 0; j < nbFaces && !isSharedEdge; ++j )
+      {
+        if ( i == j ) continue;
+        const TopoDS_Shape& otherF = getMeshDS()->IndexToShape( fId2Normal[j].first );
+        isSharedEdge = SMESH_MesherHelper::IsSubShape( *E, otherF );
+      }
+      if ( !isSharedEdge )
+        continue;
+      ee[ nbE ] = TopoDS::Edge( *E );
+      ee[ nbE ].Orientation( SMESH_MesherHelper::GetSubShapeOri( F, *E ));
+      if ( ++nbE == 2 )
+        break;
+    }
+
+    // get an angle between the two EDGEs
+    angles[i] = 0;
+    if ( nbE < 1 ) continue;
+    if ( nbE == 1 )
+    {
+      ee[ 1 ] == ee[ 0 ];
+    }
+    else
+    {
+      TopoDS_Vertex v10 = SMESH_MesherHelper::IthVertex( 1, ee[ 0 ]);
+      TopoDS_Vertex v01 = SMESH_MesherHelper::IthVertex( 0, ee[ 1 ]);
+      if ( !v10.IsSame( v01 ))
+        std::swap( ee[0], ee[1] );
+    }
+    angles[i] = SMESH_MesherHelper::GetAngle( ee[0], ee[1], F );
+  }
+
+  // compute a weighted normal
+  double sumAngle = 0;
+  for ( int i = 0; i < nbFaces; ++i )
+  {
+    angles[i] = ( angles[i] > 2*M_PI )  ?  0  :  M_PI - angles[i];
+    sumAngle += angles[i];
+  }
+  for ( int i = 0; i < nbFaces; ++i )
+    resNorm += angles[i] / sumAngle * fId2Normal[i].second;
+
+  return resNorm;
+}
+
 //================================================================================
 /*!
  * \brief Find 2 neigbor nodes of a node on EDGE
@@ -2054,7 +2269,7 @@ void _LayerEdge::SetDataByNeighbors( const SMDS_MeshNode* n1,
  */
 //================================================================================
 
-void _LayerEdge::Copy( _LayerEdge& other, SMESH_MesherHelper& helper )
+gp_XYZ _LayerEdge::Copy( _LayerEdge& other, SMESH_MesherHelper& helper )
 {
   _nodes     = other._nodes;
   _normal    = other._normal;
@@ -2066,16 +2281,25 @@ void _LayerEdge::Copy( _LayerEdge& other, SMESH_MesherHelper& helper )
   _curvature = 0; std::swap( _curvature, other._curvature );
   _2neibors  = 0; std::swap( _2neibors,  other._2neibors );
 
+  gp_XYZ lastPos( 0,0,0 );
   if ( _sWOL.ShapeType() == TopAbs_EDGE )
   {
     double u = helper.GetNodeU( TopoDS::Edge( _sWOL ), _nodes[0] );
     _pos.push_back( gp_XYZ( u, 0, 0));
+
+    u = helper.GetNodeU( TopoDS::Edge( _sWOL ), _nodes.back() );
+    lastPos.SetX( u );
   }
   else // TopAbs_FACE
   {
     gp_XY uv = helper.GetNodeUV( TopoDS::Face( _sWOL ), _nodes[0]);
     _pos.push_back( gp_XYZ( uv.X(), uv.Y(), 0));
+
+    uv = helper.GetNodeUV( TopoDS::Face( _sWOL ), _nodes.back() );
+    lastPos.SetX( uv.X() );
+    lastPos.SetY( uv.Y() );
   }
+  return lastPos;
 }
 
 //================================================================================
@@ -2087,7 +2311,7 @@ void _LayerEdge::Copy( _LayerEdge& other, SMESH_MesherHelper& helper )
 void _LayerEdge::SetCosin( double cosin )
 {
   _cosin = cosin;
-  _lenFactor = ( _cosin > 0.1 ) ?  1./sqrt(1-_cosin*_cosin) : 1.0;
+  _lenFactor = ( Abs( _cosin ) > 0.1 ) ?  1./sqrt(1-_cosin*_cosin) : 1.0;
 }
 
 //================================================================================
@@ -2379,8 +2603,8 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data,
     else
     {
       // smooth on FACE's
-      int step = 0, badNb = 0; moved = true;
-      while (( ++step <= 5 && moved ) || improved )
+      int step = 0, stepLimit = 5, badNb = 0; moved = true;
+      while (( ++step <= stepLimit && moved ) || improved )
       {
         dumpFunction(SMESH_Comment("smooth")<<data._index<<"_Fa"<<sInd
                      <<"_InfStep"<<nbSteps<<"_"<<step); // debug
@@ -2391,6 +2615,10 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data,
           moved |= data._edges[i]->Smooth(badNb);
         improved = ( badNb < oldBadNb );
 
+        // issue 22576. no bad faces but still there are intersections to fix
+        if ( improved && badNb == 0 )
+          stepLimit = step + 3;
+
         dumpFunctionEnd();
       }
       if ( badNb > 0 )
@@ -2415,6 +2643,24 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data,
     }
   } // loop on shapes to smooth
 
+  // Check orientation of simplices of _simplexTestEdges
+  for ( size_t i = 0; i < data._simplexTestEdges.size(); ++i )
+  {
+    const _LayerEdge* edge = data._simplexTestEdges[i];
+    SMESH_TNodeXYZ tgtXYZ( edge->_nodes.back() );
+    for ( size_t j = 0; j < edge->_simplices.size(); ++j )
+      if ( !edge->_simplices[j].IsForward( edge->_nodes[0], &tgtXYZ ))
+      {
+#ifdef __myDEBUG
+        cout << "Bad simplex of _simplexTestEdges ("
+             << " "<< edge->_nodes[0]->GetID()<< " "<< tgtXYZ._node->GetID()
+             << " "<< edge->_simplices[j]._nPrev->GetID()
+             << " "<< edge->_simplices[j]._nNext->GetID() << " )" << endl;
+#endif
+        return false;
+      }
+  }
+
   // Check if the last segments of _LayerEdge intersects 2D elements;
   // checked elements are either temporary faces or faces on surfaces w/o the layers
 
@@ -2867,19 +3113,19 @@ bool _ViscousBuilder::updateNormals( _SolidData&         data,
       if ( FF1[0].IsNull() || FF2[0].IsNull() )
         continue;
 
-//       // get a new normal for edge1
+      // 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();
+      // 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 );
@@ -2893,7 +3139,7 @@ bool _ViscousBuilder::updateNormals( _SolidData&         data,
       n1 = edge1->_2neibors->_edges[0]->_nodes[0];
       n2 = edge1->_2neibors->_edges[1]->_nodes[0];
       //if ( !findNeiborsOnEdge( edge1, n1, n2, data ))
-      //continue;
+      //  continue;
       edge1->SetDataByNeighbors( n1, n2, helper );
       gp_Vec dirInFace;
       if ( edge1->_cosin < 0 )
@@ -3429,10 +3675,14 @@ bool _ViscousBuilder::refine(_SolidData& data)
   Handle(Geom_Surface) surface;
   TopoDS_Edge geomEdge;
   TopoDS_Face geomFace;
+  TopoDS_Shape prevSWOL;
   TopLoc_Location loc;
-  double f,l, u/*, distXYZ[4]*/;
+  double f,l, u;
   gp_XY uv;
   bool isOnEdge;
+  TGeomID prevBaseId = -1;
+  TNode2Edge* n2eMap = 0;
+  TNode2Edge::iterator n2e;
 
   for ( size_t i = 0; i < data._edges.size(); ++i )
   {
@@ -3452,27 +3702,46 @@ bool _ViscousBuilder::refine(_SolidData& data)
       edge._nodes[1] = 0;
       edge._nodes.back() = tgtNode;
     }
-    if ( !edge._sWOL.IsNull() )
+    // get data of a shrink shape
+    if ( !edge._sWOL.IsNull() && edge._sWOL != prevSWOL )
     {
       isOnEdge = ( edge._sWOL.ShapeType() == TopAbs_EDGE );
-      // restore position of the last node
-//       gp_Pnt p;
       if ( isOnEdge )
       {
         geomEdge = TopoDS::Edge( edge._sWOL );
-        curve = BRep_Tool::Curve( geomEdge, loc, f,l);
-//         double u = helper.GetNodeU( tgtNode );
-//         p = curve->Value( u );
+        curve    = BRep_Tool::Curve( geomEdge, loc, f,l);
       }
       else
       {
         geomFace = TopoDS::Face( edge._sWOL );
-        surface = BRep_Tool::Surface( geomFace, loc );
-//         gp_XY uv = helper.GetNodeUV( tgtNode );
-//         p = surface->Value( uv.X(), uv.Y() );
+        surface  = BRep_Tool::Surface( geomFace, loc );
+      }
+      prevSWOL = edge._sWOL;
+    }
+    // restore shapePos of the last node by already treated _LayerEdge of another _SolidData
+    const TGeomID baseShapeId = edge._nodes[0]->getshapeId();
+    if ( baseShapeId != prevBaseId )
+    {
+      map< TGeomID, TNode2Edge* >::iterator s2ne = data._s2neMap.find( baseShapeId );
+      n2eMap = ( s2ne == data._s2neMap.end() ) ? 0 : n2eMap = s2ne->second;
+      prevBaseId = baseShapeId;
+    }
+    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();
+      if ( isOnEdge )
+      {
+        SMDS_EdgePosition* epos = static_cast<SMDS_EdgePosition*>( lastPos );
+        epos->SetUParameter( foundPos.X() );
+      }
+      else
+      {
+        SMDS_FacePosition* fpos = static_cast<SMDS_FacePosition*>( lastPos );
+        fpos->SetUParameter( foundPos.X() );
+        fpos->SetVParameter( foundPos.Y() );
       }
-//       p.Transform( loc );
-//       const_cast< SMDS_MeshNode* >( tgtNode )->setXYZ( p.X(), p.Y(), p.Z() );
     }
     // calculate height of the first layer
     double h0;
@@ -3510,12 +3779,14 @@ bool _ViscousBuilder::refine(_SolidData& data)
         if ( isOnEdge )
         {
           u = pos.X();
-          pos = curve->Value( u ).Transformed(loc);
+          if ( !node )
+            pos = curve->Value( u ).Transformed(loc);
         }
         else
         {
           uv.SetCoord( pos.X(), pos.Y() );
-          pos = surface->Value( pos.X(), pos.Y() ).Transformed(loc);
+          if ( !node )
+            pos = surface->Value( pos.X(), pos.Y() ).Transformed(loc);
         }
       }
       // create or update the node
@@ -3543,11 +3814,18 @@ bool _ViscousBuilder::refine(_SolidData& data)
           {
             u = 0.5 * ( u + helper.GetNodeU( geomEdge, node ));
             pos = curve->Value( u ).Transformed(loc);
+
+            SMDS_EdgePosition* epos = static_cast<SMDS_EdgePosition*>( node->GetPosition() );
+            epos->SetUParameter( u );
           }
           else
           {
             uv = 0.5 * ( uv + helper.GetNodeUV( geomFace, node ));
             pos = surface->Value( uv.X(), uv.Y()).Transformed(loc);
+
+            SMDS_FacePosition* fpos = static_cast<SMDS_FacePosition*>( node->GetPosition() );
+            fpos->SetUParameter( uv.X() );
+            fpos->SetVParameter( uv.Y() );
           }
         }
         node->setXYZ( pos.X(), pos.Y(), pos.Z() );