Salome HOME
PR: link problem detection on ubuntu 12.04
[modules/smesh.git] / src / StdMeshers / StdMeshers_ViscousLayers2D.cxx
index 6d9bc78747ff83bd5738fdc1bc5fb87a68656add..2fbab5f95cb38fd2e781db67fda89752bcf15a6a 100644 (file)
@@ -53,6 +53,8 @@
 #include <Bnd_B3d.hxx>
 #include <ElCLib.hxx>
 #include <GCPnts_AbscissaPoint.hxx>
+#include <Geom2dAdaptor_Curve.hxx>
+#include <Geom2dInt_GInter.hxx>
 #include <Geom2d_Circle.hxx>
 #include <Geom2d_Line.hxx>
 #include <Geom2d_TrimmedCurve.hxx>
@@ -61,6 +63,7 @@
 #include <Geom_Curve.hxx>
 #include <Geom_Line.hxx>
 #include <Geom_TrimmedCurve.hxx>
+#include <IntRes2d_IntersectionPoint.hxx>
 #include <Precision.hxx>
 #include <Standard_ErrorHandler.hxx>
 #include <TColStd_Array1OfReal.hxx>
@@ -239,6 +242,7 @@ namespace VISCOUS_2D
     StdMeshers_FaceSide* _wire;
     int                  _edgeInd;     // index of my EDGE in _wire
     bool                 _advancable;  // true if there is a viscous layer on my EDGE
+    bool                 _isStraight2D;// pcurve type
     _PolyLine*           _leftLine;    // lines of neighbour EDGE's
     _PolyLine*           _rightLine;
     int                  _firstPntInd; // index in vector<UVPtStruct> of _wire
@@ -257,10 +261,19 @@ namespace VISCOUS_2D
     typedef vector< _Segment >::iterator   TSegIterator;
     typedef vector< _LayerEdge >::iterator TEdgeIterator;
 
+    TIDSortedElemSet     _newFaces; // faces generated from this line
+
     bool IsCommonEdgeShared( const _PolyLine& other );
-    size_t FirstLEdge() const { return _leftLine->_advancable ? 1 : 0; }
-    bool IsAdjacent( const _Segment& seg ) const
+    size_t FirstLEdge() const
+    {
+      return ( _leftLine->_advancable && _lEdges.size() > 2 ) ? 1 : 0;
+    }
+    bool IsAdjacent( const _Segment& seg, const _LayerEdge* LE=0 ) const
     {
+      if ( LE && seg._indexInLine < _lEdges.size() &&
+           ( seg._uv[0] == & LE->_uvIn ||
+             seg._uv[1] == & LE->_uvIn ))
+        return true;
       return ( & seg == &_leftLine->_segments.back() ||
                & seg == &_rightLine->_segments[0] );
     }
@@ -278,6 +291,7 @@ namespace VISCOUS_2D
 
     bool Compute(const _Segment& seg1, const _Segment& seg2, bool seg2IsRay = false )
     {
+      const double eps = 1e-10;
       _vec1  = seg1.p2() - seg1.p1(); 
       _vec2  = seg2.p2() - seg2.p1(); 
       _vec21 = seg1.p1() - seg2.p1(); 
@@ -285,10 +299,10 @@ namespace VISCOUS_2D
       if ( fabs(_D) < std::numeric_limits<double>::min())
         return false;
       _param1 = _vec2.Crossed(_vec21) / _D; 
-      if (_param1 < 0 || _param1 > 1 )
+      if (_param1 < -eps || _param1 > 1 + eps )
         return false;
       _param2 = _vec1.Crossed(_vec21) / _D; 
-      if (_param2 < 0 || ( !seg2IsRay && _param2 > 1 ))
+      if (_param2 < -eps || ( !seg2IsRay && _param2 > 1 + eps ))
         return false;
       return true;
     }
@@ -326,10 +340,15 @@ namespace VISCOUS_2D
     bool fixCollisions();
     bool refine();
     bool shrink();
+    bool improve();
     bool toShrinkForAdjacent( const TopoDS_Face& adjFace,
                               const TopoDS_Edge& E,
                               const TopoDS_Vertex& V);
     void setLenRatio( _LayerEdge& LE, const gp_Pnt& pOut );
+    void setLayerEdgeData( _LayerEdge&           lEdge,
+                           const double          u,
+                           Handle(Geom2d_Curve)& pcurve,
+                           const bool            reverse);
     void adjustCommonEdge( _PolyLine& LL, _PolyLine& LR );
     void calcLayersHeight(const double    totalThick,
                           vector<double>& heights);
@@ -421,6 +440,8 @@ StdMeshers_ViscousLayers2D::Compute(SMESH_Mesh&        theMesh,
       theMesh.GetSubMesh( theFace )->GetComputeError() = error;
     else if ( !pm )
       pm.reset( new SMESH_ProxyMesh( theMesh ));
+    if ( getenv("ONLY_VL2D"))
+      pm.reset();
   }
   else
   {
@@ -516,8 +537,6 @@ SMESH_ProxyMesh::Ptr _ViscousBuilder2D::Compute()
   if ( !_error->IsOK() )
     return _proxyMesh;
 
-  //PyDump debugDump;
-
   if ( !findEdgesWithLayers() ) // analysis of a shape
     return _proxyMesh;
 
@@ -533,11 +552,7 @@ SMESH_ProxyMesh::Ptr _ViscousBuilder2D::Compute()
   if ( ! refine() ) // make faces
     return _proxyMesh;
 
-  // for ( size_t i = 0; i < _facesToRecompute.size(); ++i )
-  //   _mesh->GetSubMesh( _facesToRecompute[i] )->ComputeStateEngine( SMESH_subMesh::COMPUTE );
-
-  //makeGroupOfLE(); // debug
-  //debugDump.Finish();
+  //improve();
 
   return _proxyMesh;
 }
@@ -658,30 +673,29 @@ bool _ViscousBuilder2D::makePolyLines()
       while ( points[ iPnt ].normParam < lastNormPar )
         ++iPnt;
       L._lastPntInd = iPnt;
-      L._lEdges.resize( L._lastPntInd - L._firstPntInd + 1 );
+      L._lEdges.resize( Max( 3, L._lastPntInd - L._firstPntInd + 1 )); // 3 edges minimum
 
       // TODO: add more _LayerEdge's to strongly curved EDGEs
       // in order not to miss collisions
 
       Handle(Geom2d_Curve) pcurve = L._wire->Curve2d( L._edgeInd );
-      gp_Pnt2d uv; gp_Vec2d tangent;
+      const bool reverse = ( L._wire->Edge( iE ).Orientation() == TopAbs_REVERSED );
       for ( int i = L._firstPntInd; i <= L._lastPntInd; ++i )
       {
         _LayerEdge& lEdge = L._lEdges[ i - L._firstPntInd ];
         const double u = ( i == L._firstPntInd ? wire->FirstU(iE) : points[ i ].param );
-        pcurve->D1( u , uv, tangent );
-        tangent.Normalize();
-        if ( L._wire->Edge( iE ).Orientation() == TopAbs_REVERSED )
-          tangent.Reverse();
-        lEdge._uvOut = lEdge._uvIn = uv.XY();
-        lEdge._normal2D.SetCoord( -tangent.Y(), tangent.X() );
-        lEdge._ray.SetLocation( lEdge._uvOut );
-        lEdge._ray.SetDirection( lEdge._normal2D );
-        lEdge._isBlocked = false;
-        lEdge._length2D  = 0;
-
+        setLayerEdgeData( lEdge, u, pcurve, reverse );
         setLenRatio( lEdge, SMESH_TNodeXYZ( points[ i ].node ) );
       }
+      if ( L._lastPntInd - L._firstPntInd + 1 < 3 ) // add 3d _LayerEdge in the middle
+      {
+        L._lEdges[2] = L._lEdges[1];
+        const double u = 0.5 * ( wire->FirstU(iE) + wire->LastU(iE) );
+        setLayerEdgeData( L._lEdges[1], u, pcurve, reverse );
+        gp_Pnt p = 0.5 * ( SMESH_TNodeXYZ( points[ L._firstPntInd ].node ) +
+                           SMESH_TNodeXYZ( points[ L._lastPntInd ].node ));
+        setLenRatio( L._lEdges[1], p );
+      }
     }
   }
 
@@ -710,13 +724,14 @@ bool _ViscousBuilder2D::makePolyLines()
     L._segTree.reset( new _SegmentTree( L._segments ));
   }
 
-  // Evaluate possible _thickness if required layers thickness seems too high
-  // -------------------------------------------------------------------------
+  // Evaluate max possible _thickness if required layers thickness seems too high
+  // ----------------------------------------------------------------------------
 
   _thickness = _hyp->GetTotalThickness();
   _SegmentTree::box_type faceBndBox2D;
   for ( iPoLine = 0; iPoLine < _polyLineVec.size(); ++iPoLine )
     faceBndBox2D.Add( *_polyLineVec[ iPoLine]._segTree->getBox() );
+  double boxTol = 1e-3 * sqrt( faceBndBox2D.SquareExtent() );
   //
   if ( _thickness * maxLen2dTo3dRatio > sqrt( faceBndBox2D.SquareExtent() ) / 10 )
   {
@@ -726,9 +741,15 @@ bool _ViscousBuilder2D::makePolyLines()
     for ( size_t iL1 = 0; iL1 < _polyLineVec.size(); ++iL1 )
     {
       _PolyLine& L1 = _polyLineVec[ iL1 ];
+      _SegmentTree::box_type boxL1 = * L1._segTree->getBox();
+      boxL1.Enlarge( boxTol );
       for ( size_t iL2 = iL1+1; iL2 < _polyLineVec.size(); ++iL2 )
       {
         _PolyLine& L2 = _polyLineVec[ iL2 ];
+        _SegmentTree::box_type boxL2 = * L2._segTree->getBox();
+        boxL2.Enlarge( boxTol );
+        if ( boxL1.IsOut( boxL2 ))
+          continue;
         for ( size_t iLE = 1; iLE < L1._lEdges.size(); ++iLE )
         {
           foundSegs.clear();
@@ -790,16 +811,16 @@ bool _ViscousBuilder2D::makePolyLines()
   {
     lineBoxes[ iPoLine ] = *_polyLineVec[ iPoLine ]._segTree->getBox();
     if ( _polyLineVec[ iPoLine ]._advancable )
-      lineBoxes[ iPoLine ].Enlarge( maxLen2dTo3dRatio * _thickness );
+      lineBoxes[ iPoLine ].Enlarge( maxLen2dTo3dRatio * _thickness * 2 );
   }
   // _reachableLines
   for ( iPoLine = 0; iPoLine < _polyLineVec.size(); ++iPoLine )
   {
     _PolyLine& L1 = _polyLineVec[ iPoLine ];
-    for ( size_t i = 0; i < _polyLineVec.size(); ++i )
+    for ( size_t iL2 = 0; iL2 < _polyLineVec.size(); ++iL2 )
     {
-      _PolyLine& L2 = _polyLineVec[ i ];
-      if ( iPoLine == i || lineBoxes[ iPoLine ].IsOut( lineBoxes[ i ]))
+      _PolyLine& L2 = _polyLineVec[ iL2 ];
+      if ( iPoLine == iL2 || lineBoxes[ iPoLine ].IsOut( lineBoxes[ iL2 ]))
         continue;
       if ( !L1._advancable && ( L1._leftLine == &L2 || L1._rightLine == &L2 ))
         continue;
@@ -808,10 +829,8 @@ bool _ViscousBuilder2D::makePolyLines()
       for ( size_t iLE = 1; iLE < L1._lEdges.size(); iLE += iDelta )
       {
         _LayerEdge& LE = L1._lEdges[iLE];
-        if ( !lineBoxes[ i ].IsOut ( LE._uvOut,
-                                     LE._uvOut + LE._normal2D * _thickness * LE._len2dTo3dRatio )
-             &&
-             !L1.IsAdjacent( L2._segments[0] ))
+        if ( !lineBoxes[ iL2 ].IsOut ( LE._uvOut,
+                                       LE._uvOut + LE._normal2D *_thickness * LE._len2dTo3dRatio ))
         {
           L1._reachableLines.push_back( & L2 );
           break;
@@ -820,7 +839,8 @@ bool _ViscousBuilder2D::makePolyLines()
     }
     // add self to _reachableLines
     Geom2dAdaptor_Curve pcurve( L1._wire->Curve2d( L1._edgeInd ));
-    if ( pcurve.GetType() != GeomAbs_Line )
+    L1._isStraight2D = ( pcurve.GetType() == GeomAbs_Line );
+    if ( !L1._isStraight2D )
     {
       // TODO: check carefully
       L1._reachableLines.push_back( & L1 );
@@ -849,14 +869,17 @@ void _ViscousBuilder2D::adjustCommonEdge( _PolyLine& LL, _PolyLine& LR )
   gp_XY normL    = EL._normal2D;
   gp_XY normR    = ER._normal2D;
   gp_XY tangL ( normL.Y(), -normL.X() );
-  //gp_XY tangR ( normR.Y(), -normR.X() );
-
-  gp_XY normCommon = ( normL + normR ).Normalized(); // average normal at VERTEX
 
+  // set common direction to a VERTEX _LayerEdge shared by two _PolyLine's
+  gp_XY normCommon = ( normL * int( LL._advancable ) +
+                       normR * int( LR._advancable )).Normalized();
   EL._normal2D = normCommon;
   EL._ray.SetLocation ( EL._uvOut );
   EL._ray.SetDirection( EL._normal2D );
-
+  if ( nbAdvancableL == 1 ) { // _normal2D is true normal (not average)
+    EL._isBlocked = true; // prevent intersecting with _Segments of _advancable line
+    EL._length2D  = 0;
+  }
   // update _LayerEdge::_len2dTo3dRatio according to a new direction
   const vector<UVPtStruct>& points = LL._wire->GetUVPtStruct();
   setLenRatio( EL, SMESH_TNodeXYZ( points[ LL._lastPntInd ].node ));
@@ -865,48 +888,85 @@ void _ViscousBuilder2D::adjustCommonEdge( _PolyLine& LL, _PolyLine& LR )
 
   const double dotNormTang = normR * tangL;
   const bool    largeAngle = Abs( dotNormTang ) > 0.2;
-  if ( largeAngle )
+  if ( largeAngle ) // not 180 degrees
   {
     // recompute _len2dTo3dRatio to take into account angle between EDGEs
     gp_Vec2d oldNorm( LL._advancable ? normL : normR );
-    double fact = 1. / Max( 0.3, Cos( oldNorm.Angle(  normCommon )));
-    EL._len2dTo3dRatio *= fact;
+    double angleFactor  = 1. / Max( 0.3, Cos( oldNorm.Angle( normCommon )));
+    EL._len2dTo3dRatio *= angleFactor;
     ER._len2dTo3dRatio  = EL._len2dTo3dRatio;
 
+    gp_XY normAvg = ( normL + normR ).Normalized(); // average normal at VERTEX
+
     if ( dotNormTang < 0. ) // ---------------------------- CONVEX ANGLE
     {
-      // Remove _LayerEdge's intersecting the normCommon
+      // Remove _LayerEdge's intersecting the normAvg to avoid collisions
+      // during inflate().
       //
+      // find max length of the VERTEX based _LayerEdge whose direction is normAvg
+      double maxLen2D       = _thickness * EL._len2dTo3dRatio;
       const gp_XY& pCommOut = ER._uvOut;
-      gp_XY pCommIn( pCommOut + normCommon * _thickness * EL._len2dTo3dRatio );
+      gp_XY        pCommIn  = pCommOut + normAvg * maxLen2D;
       _Segment segCommon( pCommOut, pCommIn );
       _SegmentIntersection intersection;
+      vector< const _Segment* > foundSegs;
+      for ( size_t iL1 = 0; iL1 < _polyLineVec.size(); ++iL1 )
+      {
+        _PolyLine& L1 = _polyLineVec[ iL1 ];
+        const _SegmentTree::box_type* boxL1 = L1._segTree->getBox();
+        if ( boxL1->IsOut ( pCommOut, pCommIn ))
+          continue;
+        for ( size_t iLE = 1; iLE < L1._lEdges.size(); ++iLE )
+        {
+          foundSegs.clear();
+          L1._segTree->GetSegmentsNear( segCommon, foundSegs );
+          for ( size_t i = 0; i < foundSegs.size(); ++i )
+            if ( intersection.Compute( *foundSegs[i], segCommon ) &&
+                 intersection._param2 > 1e-10 )
+            {
+              double len2D = intersection._param2 * maxLen2D / ( 2 + L1._advancable );
+              if ( len2D < maxLen2D ) {
+                maxLen2D = len2D;
+                pCommIn  = pCommOut + normAvg * maxLen2D; // here length of segCommon changes
+              }
+            }
+        }
+      }
+
+      // remove _LayerEdge's intersecting segCommon
       for ( int isR = 0; isR < 2; ++isR ) // loop on [ LL, LR ]
       {
         _PolyLine&                 L = isR ? LR : LL;
         _PolyLine::TEdgeIterator eIt = isR ? L._lEdges.begin()+1 : L._lEdges.end()-2;
         int                      dIt = isR ? +1 : -1;
-        // at least 2 _LayerEdge's should remain in a _PolyLine (if _advancable)
-        if ( L._lEdges.size() < 3 ) continue;
+        if ( nbAdvancableL == 1 && L._advancable && normL * normR > -0.01 )
+          continue;  // obtuse internal angle
+        // at least 3 _LayerEdge's should remain in a _PolyLine
+        if ( L._lEdges.size() < 4 ) continue;
         size_t iLE = 1;
+        _SegmentIntersection lastIntersection;
         for ( ; iLE < L._lEdges.size(); ++iLE, eIt += dIt )
         {
           gp_XY uvIn = eIt->_uvOut + eIt->_normal2D * _thickness * eIt->_len2dTo3dRatio;
           _Segment segOfEdge( eIt->_uvOut, uvIn );
           if ( !intersection.Compute( segCommon, segOfEdge ))
             break;
+          lastIntersection._param1 = intersection._param1;
+          lastIntersection._param2 = intersection._param2;
         }
         if ( iLE >= L._lEdges.size () - 1 )
         {
           // all _LayerEdge's intersect the segCommon, limit inflation
           // of remaining 2 _LayerEdge's
-          vector< _LayerEdge > newEdgeVec( );
+          vector< _LayerEdge > newEdgeVec( Min( 3, L._lEdges.size() ));
           newEdgeVec.front() = L._lEdges.front();
           newEdgeVec.back()  = L._lEdges.back();
+          if ( newEdgeVec.size() == 3 )
+            newEdgeVec[1] = L._lEdges[ L._lEdges.size() / 2 ];
           L._lEdges.swap( newEdgeVec );
-          if ( !isR ) std::swap( intersection._param1 , intersection._param2 );
-          L._lEdges.front()._len2dTo3dRatio *= intersection._param1;
-          L._lEdges.back ()._len2dTo3dRatio *= intersection._param2;
+          if ( !isR ) std::swap( lastIntersection._param1 , lastIntersection._param2 );
+          L._lEdges.front()._len2dTo3dRatio *= lastIntersection._param1; // ??
+          L._lEdges.back ()._len2dTo3dRatio *= lastIntersection._param2;
         }
         else if ( iLE != 1 )
         {
@@ -922,14 +982,44 @@ void _ViscousBuilder2D::adjustCommonEdge( _PolyLine& LL, _PolyLine& LR )
     {
       if ( nbAdvancableL == 1 )
       {
-        // make that the _LayerEdge at VERTEX is not shared by LL and LR
+        // make that the _LayerEdge at VERTEX is not shared by LL and LR:
+        // different normals is a sign that they are not shared
         _LayerEdge& notSharedEdge = LL._advancable ? LR._lEdges[0] : LL._lEdges.back();
+        _LayerEdge&    sharedEdge = LR._advancable ? LR._lEdges[0] : LL._lEdges.back();
+
         notSharedEdge._normal2D.SetCoord( 0.,0. );
+        sharedEdge._normal2D     = normAvg;
+        sharedEdge._isBlocked    = false;
+        notSharedEdge._isBlocked = true;
       }
     }
   }
 }
 
+//================================================================================
+/*!
+ * \brief initialize data of a _LayerEdge
+ */
+//================================================================================
+
+void _ViscousBuilder2D::setLayerEdgeData( _LayerEdge&           lEdge,
+                                          const double          u,
+                                          Handle(Geom2d_Curve)& pcurve,
+                                          const bool            reverse)
+{
+  gp_Pnt2d uv; gp_Vec2d tangent;
+  pcurve->D1( u, uv, tangent );
+  tangent.Normalize();
+  if ( reverse )
+    tangent.Reverse();
+  lEdge._uvOut = lEdge._uvIn = uv.XY();
+  lEdge._normal2D.SetCoord( -tangent.Y(), tangent.X() );
+  lEdge._ray.SetLocation( lEdge._uvOut );
+  lEdge._ray.SetDirection( lEdge._normal2D );
+  lEdge._isBlocked = false;
+  lEdge._length2D  = 0;
+}
+
 //================================================================================
 /*!
  * \brief Compute and set _LayerEdge::_len2dTo3dRatio
@@ -973,7 +1063,7 @@ bool _ViscousBuilder2D::inflate()
         foundSegs.clear();
         L2._segTree->GetSegmentsNear( L1._lEdges[iLE]._ray, foundSegs );
         for ( size_t i = 0; i < foundSegs.size(); ++i )
-          if ( ! L1.IsAdjacent( *foundSegs[i] ) &&
+          if ( ! L1.IsAdjacent( *foundSegs[i], & L1._lEdges[iLE] ) &&
                intersection.Compute( *foundSegs[i], L1._lEdges[iLE]._ray ))
           {
             double distToL2 = intersection._param2 / L1._lEdges[iLE]._len2dTo3dRatio;
@@ -1031,6 +1121,57 @@ bool _ViscousBuilder2D::inflate()
   // if (nbSteps == 0 )
   //   return error("failed at the very first inflation step");
 
+
+  // remove _LayerEdge's of one line intersecting with each other
+  for ( size_t iL = 0; iL < _polyLineVec.size(); ++iL )
+  {
+    _PolyLine& L = _polyLineVec[ iL ];
+    if ( !L._advancable ) continue;
+
+    // replace an inactive (1st) _LayerEdge with an active one of a neighbour _PolyLine
+    if ( /*!L._leftLine->_advancable &&*/ L.IsCommonEdgeShared( *L._leftLine ) ) {
+      L._lEdges[0] = L._leftLine->_lEdges.back();
+    }
+    if ( !L._rightLine->_advancable && L.IsCommonEdgeShared( *L._rightLine ) ) {
+      L._lEdges.back() = L._rightLine->_lEdges[0];
+    }
+
+    _SegmentIntersection intersection;
+    for ( int isR = 0; ( isR < 2 && L._lEdges.size() > 2 ); ++isR )
+    {
+      int nbRemove = 0, deltaIt = isR ? -1 : +1;
+      _PolyLine::TEdgeIterator eIt = isR ? L._lEdges.end()-1 : L._lEdges.begin();
+      if ( eIt->_length2D == 0 ) continue;
+      _Segment seg1( eIt->_uvOut, eIt->_uvIn );
+      for ( eIt += deltaIt; nbRemove < L._lEdges.size()-1; eIt += deltaIt )
+      {
+        _Segment seg2( eIt->_uvOut, eIt->_uvIn );
+        if ( !intersection.Compute( seg1, seg2 ))
+          break;
+        ++nbRemove;
+      }
+      if ( nbRemove > 0 ) {
+        if ( nbRemove == L._lEdges.size()-1 ) // 1st and last _LayerEdge's intersect
+        {
+          --nbRemove;
+          _LayerEdge& L0 = L._lEdges.front();
+          _LayerEdge& L1 = L._lEdges.back();
+          L0._length2D *= intersection._param1 * 0.5;
+          L1._length2D *= intersection._param2 * 0.5;
+          L0._uvIn = L0._uvOut + L0._normal2D * L0._length2D;
+          L1._uvIn = L1._uvOut + L1._normal2D * L1._length2D;
+          if ( L.IsCommonEdgeShared( *L._leftLine ))
+            L._leftLine->_lEdges.back() = L0;
+        }
+        if ( isR )
+          L._lEdges.erase( L._lEdges.end()-nbRemove-1,
+                           L._lEdges.end()-nbRemove );
+        else
+          L._lEdges.erase( L._lEdges.begin()+1,
+                           L._lEdges.begin()+1+nbRemove );
+      }
+    }
+  }
   return true;
 }
 
@@ -1049,6 +1190,7 @@ bool _ViscousBuilder2D::fixCollisions()
   _SegmentIntersection intersection;
 
   list< pair< _LayerEdge*, double > > edgeLenLimitList;
+  list< _LayerEdge* >                 blockedEdgesList;
 
   for ( size_t iL1 = 0; iL1 < _polyLineVec.size(); ++iL1 )
   {
@@ -1057,8 +1199,7 @@ bool _ViscousBuilder2D::fixCollisions()
     for ( size_t iL2 = 0; iL2 < L1._reachableLines.size(); ++iL2 )
     {
       _PolyLine& L2 = * L1._reachableLines[ iL2 ];
-      //for ( size_t iLE = L1.FirstLEdge(); iLE < L1._lEdges.size(); ++iLE )
-      for ( size_t iLE = 1; iLE < L1._lEdges.size()-1; ++iLE )
+      for ( size_t iLE = L1.FirstLEdge(); iLE < L1._lEdges.size(); ++iLE )
       {
         _LayerEdge& LE1 = L1._lEdges[iLE];
         if ( LE1._isBlocked ) continue;
@@ -1066,7 +1207,7 @@ bool _ViscousBuilder2D::fixCollisions()
         L2._segTree->GetSegmentsNear( LE1._ray, foundSegs );
         for ( size_t i = 0; i < foundSegs.size(); ++i )
         {
-          if ( ! L1.IsAdjacent( *foundSegs[i] ) &&
+          if ( ! L1.IsAdjacent( *foundSegs[i], &LE1 ) &&
                intersection.Compute( *foundSegs[i], LE1._ray ))
           {
             const double dist2DToL2 = intersection._param2;
@@ -1075,11 +1216,12 @@ bool _ViscousBuilder2D::fixCollisions()
             {
               if ( newLen2D < LE1._length2D )
               {
+                blockedEdgesList.push_back( &LE1 );
                 if ( L1._advancable )
                 {
                   edgeLenLimitList.push_back( make_pair( &LE1, newLen2D ));
-                  L2._lEdges[ foundSegs[i]->_indexInLine     ]._isBlocked = true;
-                  L2._lEdges[ foundSegs[i]->_indexInLine + 1 ]._isBlocked = true;
+                  blockedEdgesList.push_back( &L2._lEdges[ foundSegs[i]->_indexInLine     ]);
+                  blockedEdgesList.push_back( &L2._lEdges[ foundSegs[i]->_indexInLine + 1 ]);
                 }
                 else // here dist2DToL2 < 0 and LE1._length2D == 0
                 {
@@ -1091,21 +1233,9 @@ bool _ViscousBuilder2D::fixCollisions()
 
                   edgeLenLimitList.push_back( make_pair( &LE2[0], newLen2D ));
                   edgeLenLimitList.push_back( make_pair( &LE2[1], newLen2D ));
-                  LE2[0]._isBlocked = true;
-                  LE2[1]._isBlocked = true;
                 }
               }
-              LE1._isBlocked = true; // !! after SetNewLength()
             }
-            // else
-            // {
-            //   double step2D = newLen2D - LE1._length2D;
-            //   double step   = step2D / LE1._len2dTo3dRatio;
-            //   if ( step > maxStep )
-            //     maxStep = step;
-            //   if ( step < minStep )
-            //     minStep = step;
-            // }
           }
         }
       }
@@ -1118,8 +1248,15 @@ bool _ViscousBuilder2D::fixCollisions()
   {
     _LayerEdge* LE = edge2Len->first;
     LE->SetNewLength( edge2Len->second / LE->_len2dTo3dRatio );
+    LE->_isBlocked = true;
   }
 
+  // block inflation of _LayerEdge's
+  list< _LayerEdge* >::iterator edge = blockedEdgesList.begin();
+  for ( ; edge != blockedEdgesList.end(); ++edge )
+    (*edge)->_isBlocked = true;
+
+  // find a not blocked _LayerEdge
   for ( size_t iL = 0; iL < _polyLineVec.size(); ++iL )
   {
     _PolyLine& L = _polyLineVec[ iL ];
@@ -1141,7 +1278,7 @@ bool _ViscousBuilder2D::fixCollisions()
 
 bool _ViscousBuilder2D::shrink()
 {
-  gp_Pnt2d uv; gp_Vec2d tangent;
+  gp_Pnt2d uv; //gp_Vec2d tangent;
   _SegmentIntersection intersection;
   double sign;
 
@@ -1150,7 +1287,8 @@ bool _ViscousBuilder2D::shrink()
     _PolyLine& L = _polyLineVec[ iL1 ]; // line with no layers
     if ( L._advancable )
       continue;
-    if ( !L._rightLine->_advancable && !L._leftLine->_advancable )
+    const int nbAdvancable = ( L._rightLine->_advancable + L._leftLine->_advancable );
+    if ( nbAdvancable == 0 )
       continue;
 
     const TopoDS_Edge&        E = L._wire->Edge      ( L._edgeInd );
@@ -1234,7 +1372,12 @@ bool _ViscousBuilder2D::shrink()
     double u1 = L._wire->FirstU( L._edgeInd ), uf = u1;
     double u2 = L._wire->LastU ( L._edgeInd ), ul = u2;
 
-    // Get length of existing segments (from edge start to node) and their nodes
+    // a ratio to pass 2D <--> 1D
+    const double len1D = 1e-3;
+    const double len2D = pcurve->Value(uf).Distance( pcurve->Value(uf+len1D));
+    double len1dTo2dRatio = len1D / len2D;
+
+    // Get length of existing segments (from an edge start to a node) and their nodes
     const vector<UVPtStruct>& points = L._wire->GetUVPtStruct();
     UVPtStructVec nodeDataVec( & points[ L._firstPntInd ],
                                & points[ L._lastPntInd + 1 ]);
@@ -1275,47 +1418,73 @@ bool _ViscousBuilder2D::shrink()
       double  u0 = isR ? ul : uf; // init value of the param to move
       int  iPEnd = isR ? nodeDataVec.size() - 1 : 0;
 
+      _LayerEdge& nearLE = isR ? L._lEdges.back() : L._lEdges.front();
+      _LayerEdge&  farLE = isR ? L._lEdges.front() : L._lEdges.back();
+
       // try to find length of advancement along L by intersecting L with
       // an adjacent _Segment of L2
 
-      double & length2D = ( isR ? L._lEdges.back() : L._lEdges.front() )._length2D;
+      double& length2D = nearLE._length2D;
+      double  length1D = 0;
       sign = ( isR ^ edgeReversed ) ? -1. : 1.;
-      pcurve->D1( u, uv, tangent );
+      //pcurve->D1( u, uv, tangent );
 
+      bool isConvex = false;
       if ( L2->_advancable )
       {
-        gp_Ax2d      edgeRay( uv, tangent * sign );
-        const _Segment& seg2( isR ? L2->_segments.front() : L2->_segments.back() );
-        // make an elongated seg2
-        gp_XY seg2Vec( seg2.p2() - seg2.p1() );
-        gp_XY longSeg2p1 = seg2.p1() - 1000 * seg2Vec;
-        gp_XY longSeg2p2 = seg2.p2() + 1000 * seg2Vec;
-        _Segment longSeg2( longSeg2p1, longSeg2p2 );
-        if ( intersection.Compute( longSeg2, edgeRay )) { // convex VERTEX
-
-          length2D = intersection._param2; /*  |L  seg2     
-                                            *  |  o---o--- 
-                                            *  | /    |    
-                                            *  |/     |  L2
-                                            *  x------x---      */
+        int iFSeg2 = isR ? 0 : L2->_segments.size() - 1;
+        int iLSeg2 = isR ? 1 : L2->_segments.size() - 2;
+        gp_XY uvLSeg2In  = L2->_lEdges[ iLSeg2 ]._uvIn;
+        gp_XY uvLSeg2Out = L2->_lEdges[ iLSeg2 ]._uvOut;
+        gp_XY uvFSeg2Out = L2->_lEdges[ iFSeg2 ]._uvOut;
+        Handle(Geom2d_Line) seg2Line = new Geom2d_Line( uvLSeg2In, uvFSeg2Out - uvLSeg2Out );
+
+        Geom2dAdaptor_Curve edgeCurve( pcurve, Min( uf, ul ), Max( uf, ul ));
+        Geom2dAdaptor_Curve seg2Curve( seg2Line );
+        Geom2dInt_GInter     curveInt( edgeCurve, seg2Curve, 1e-7, 1e-7 );
+        isConvex = ( curveInt.IsDone() && !curveInt.IsEmpty() );
+        if ( isConvex ) {
+          /*                   convex VERTEX */
+          length1D = Abs( u - curveInt.Point( 1 ).ParamOnFirst() );
+          double maxDist2d = 2 * L2->_lEdges[ iLSeg2 ]._length2D;
+          isConvex = ( length1D < maxDist2d * len1dTo2dRatio );
+                                                  /*  |L  seg2     
+                                                   *  |  o---o--- 
+                                                   *  | /    |    
+                                                   *  |/     |  L2
+                                                   *  x------x---      */
         }
-        else { /* concave VERTEX */        /*  o-----o--- 
-                                            *   \    |    
-                                            *    \   |  L2
-                                            *     x--x--- 
-                                            *    /        
-                                            * L /               */
-          length2D = ( isR ? L2->_lEdges.front() : L2->_lEdges.back() )._length2D;
+        if ( !isConvex ) { /* concave VERTEX */   /*  o-----o--- 
+                                                   *   \    |    
+                                                   *    \   |  L2
+                                                   *     x--x--- 
+                                                   *    /        
+                                                   * L /               */
+          length2D = L2->_lEdges[ iFSeg2 ]._length2D;
+          //if ( L2->_advancable ) continue;
         }
       }
       else // L2 is advancable but in the face adjacent by L
       {
-        length2D = ( isR ? L._lEdges.front() : L._lEdges.back() )._length2D;
-        if ( length2D == 0 )
-          length2D = ( isR ? L._leftLine->_lEdges.back() : L._rightLine->_lEdges.front() )._length2D;
+        length2D = farLE._length2D;
+        if ( length2D == 0 ) {
+          _LayerEdge& neighborLE =
+            ( isR ? L._leftLine->_lEdges.back() : L._rightLine->_lEdges.front() );
+          length2D = neighborLE._length2D;
+          if ( length2D == 0 )
+            length2D = _thickness * nearLE._len2dTo3dRatio;
+        }
       }
+
       // move u to the internal boundary of layers
-      u += length2D * sign;
+      double maxLen3D = Min( _thickness, edgeLen / ( 1 + nbAdvancable ));
+      double maxLen2D = maxLen3D * nearLE._len2dTo3dRatio;
+      if ( !length2D ) length2D = length1D / len1dTo2dRatio;
+      if ( Abs( length2D ) > maxLen2D )
+        length2D = maxLen2D;
+      nearLE._uvIn = nearLE._uvOut + nearLE._normal2D * length2D;
+
+      u += length2D * len1dTo2dRatio * sign;
       nodeDataVec[ iPEnd ].param = u;
 
       gp_Pnt2d newUV = pcurve->Value( u );
@@ -1508,20 +1677,59 @@ bool _ViscousBuilder2D::refine()
     _PolyLine& L = _polyLineVec[ iL ];
     if ( !L._advancable ) continue;
 
-    //if ( L._leftLine->_advancable ) L._lEdges[0] = L._leftLine->_lEdges.back();
-
-    // calculate intermediate UV on _LayerEdge's ( _LayerEdge::_uvRefined )
+    // replace an inactive (1st) _LayerEdge with an active one of a neighbour _PolyLine
     size_t iLE = 0, nbLE = L._lEdges.size();
-    if ( /*!L._leftLine->_advancable &&*/ L.IsCommonEdgeShared( *L._leftLine ))
+    const bool leftEdgeShared  = L.IsCommonEdgeShared( *L._leftLine );
+    const bool rightEdgeShared = L.IsCommonEdgeShared( *L._rightLine );
+    if ( /*!L._leftLine->_advancable &&*/ leftEdgeShared )
     {
       L._lEdges[0] = L._leftLine->_lEdges.back();
       iLE += int( !L._leftLine->_advancable );
     }
-    if ( !L._rightLine->_advancable && L.IsCommonEdgeShared( *L._rightLine ))
+    if ( !L._rightLine->_advancable && rightEdgeShared )
     {
       L._lEdges.back() = L._rightLine->_lEdges[0];
       --nbLE;
     }
+
+    // limit length of neighbour _LayerEdge's to avoid sharp change of layers thickness
+    vector< double > segLen( L._lEdges.size() );
+    segLen[0] = 0.0;
+    for ( size_t i = 1; i < segLen.size(); ++i )
+    {
+      // accumulate length of segments
+      double sLen = (L._lEdges[i-1]._uvOut - L._lEdges[i]._uvOut ).Modulus();
+      segLen[i] = segLen[i-1] + sLen;
+    }
+    for ( int isR = 0; isR < 2; ++isR )
+    {
+      size_t iF = 0, iL = L._lEdges.size()-1;
+      size_t *i = isR ? &iL : &iF;
+      //size_t iRef = *i;
+      _LayerEdge* prevLE = & L._lEdges[ *i ];
+      double weight = 0;
+      for ( ++iF, --iL; iF < L._lEdges.size()-1; ++iF, --iL )
+      {
+        _LayerEdge& LE = L._lEdges[*i];
+        if ( prevLE->_length2D > 0 ) {
+          gp_XY tangent ( LE._normal2D.Y(), -LE._normal2D.X() );
+          weight += Abs( tangent * ( prevLE->_uvIn - LE._uvIn )) / segLen.back();
+          gp_XY prevTang = ( LE._uvOut - prevLE->_uvOut );
+          gp_XY prevNorm    = gp_XY( -prevTang.Y(), prevTang.X() );
+          double prevProj   = prevNorm * ( prevLE->_uvIn - prevLE->_uvOut );
+          if ( prevProj > 0 ) {
+            prevProj /= prevTang.Modulus();
+            if ( LE._length2D < prevProj )
+              weight += 0.75 * ( 1 - weight ); // length decrease is more preferable
+            LE._length2D  = weight * LE._length2D + ( 1 - weight ) * prevProj;
+            LE._uvIn = LE._uvOut + LE._normal2D * LE._length2D;
+          }
+        }
+        prevLE = & LE;
+      }
+    }
+
+    // calculate intermediate UV on _LayerEdge's ( _LayerEdge::_uvRefined )
     for ( ; iLE < nbLE; ++iLE )
     {
       _LayerEdge& LE = L._lEdges[iLE];
@@ -1554,13 +1762,11 @@ bool _ViscousBuilder2D::refine()
 
     // Create layers of faces
 
-    int hasLeftNode  = ( !L._leftLine->_rightNodes.empty() );
-    int hasRightNode = ( !L._rightLine->_leftNodes.empty() );
+    int hasLeftNode  = ( !L._leftLine->_rightNodes.empty() && leftEdgeShared );
+    int hasRightNode = ( !L._rightLine->_leftNodes.empty() && rightEdgeShared );
     size_t iS, iN0 = hasLeftNode, nbN = innerNodes.size() - hasRightNode;
     L._leftNodes .resize( _hyp->GetNumberLayers() );
     L._rightNodes.resize( _hyp->GetNumberLayers() );
-    vector< double > segLen( L._lEdges.size() );
-    segLen[0] = 0.0;
     for ( int iF = 0; iF < _hyp->GetNumberLayers(); ++iF ) // loop on layers of faces
     {
       // get accumulated length of intermediate segments
@@ -1592,11 +1798,55 @@ bool _ViscousBuilder2D::refine()
       // create faces
       // TODO care of orientation
       for ( size_t i = 1; i < innerNodes.size(); ++i )
-        _helper.AddFace( outerNodes[ i-1 ], outerNodes[ i ],
-                         innerNodes[ i ],   innerNodes[ i-1 ]);
+        if ( SMDS_MeshElement* f = _helper.AddFace( outerNodes[ i-1 ], outerNodes[ i ],
+                                                    innerNodes[ i ],   innerNodes[ i-1 ]))
+          L._newFaces.insert( L._newFaces.end(), f );
 
       outerNodes.swap( innerNodes );
     }
+    // faces between not shared _LayerEdge's (at concave VERTEX)
+    for ( int isR = 0; isR < 2; ++isR )
+    {
+      if ( isR ? rightEdgeShared : leftEdgeShared)
+        continue;
+      vector< const SMDS_MeshNode* > &
+        lNodes = (isR ? L._rightNodes : L._leftLine->_rightNodes ),
+        rNodes = (isR ? L._rightLine->_leftNodes : L._leftNodes );
+      if ( lNodes.empty() || rNodes.empty() || lNodes.size() != rNodes.size() )
+        continue;
+
+      for ( size_t i = 1; i < lNodes.size(); ++i )
+        _helper.AddFace( lNodes[ i-1 ], rNodes[ i-1 ],
+                         rNodes[ i ],   lNodes[ i ]);
+
+      const UVPtStruct& ptOnVertex = points[ isR ? L._lastPntInd : L._firstPntInd ];
+      _helper.AddFace( ptOnVertex.node, rNodes[ 0 ], lNodes[ 0 ]);
+
+      // update nodeDataVec of an adjacent _PolyLine
+      // int iAdjEdge = isR ? L._rightLine->_edgeInd : L._leftLine->_edgeInd;
+      // _ProxyMeshOfFace::_EdgeSubMesh* adjEdgeSM
+      //   = getProxyMesh()->GetEdgeSubMesh( L._wire->EdgeID( iAdjEdge ));
+      // const UVPtStructVec& nodeDataVec = adjEdgeSM->GetUVPtStructVec();
+      // if ( !nodeDataVec.empty() )
+      // {
+      //   UVPtStruct ptOnVertex;
+      //   _LayerEdge& LE = isR ? L._lEdges.back() : L._lEdges.front();
+      //   ptOnVertex.u         = LE._uvRefined.back().X();
+      //   ptOnVertex.v         = LE._uvRefined.back().Y();
+      //   ptOnVertex.node      = isR ? L._rightNodes.back() : L._leftNodes.back();
+      //   ptOnVertex.param     = isR ? L._wire->FirstU( iAdjEdge ) :L._wire->LastU( iAdjEdge );
+      //   ptOnVertex.normParam = isR ? 1 : 0;
+      //   ptOnVertex.x         = ptOnVertex.normParam;
+      //   ptOnVertex.y         = ptOnVertex.normParam;
+
+      //   int iN  = isR ? _hyp->GetNumberLayers() : 0;
+      //   int nbN = nodeDataVec.size() - ( isR ? 0 : _hyp->GetNumberLayers() );
+      //   UVPtStructVec newNodeData( nodeDataVec.begin() + iN,
+      //                              nodeDataVec.begin() + nbN );
+      //   newNodeData.insert( isR ? newNodeData.begin() : newNodeData.end(), ptOnVertex );
+      //   adjEdgeSM->SetUVPtStructVec( newNodeData );
+      // }
+    }
 
     // Fill the _ProxyMeshOfFace
 
@@ -1624,6 +1874,56 @@ bool _ViscousBuilder2D::refine()
   return true;
 }
 
+//================================================================================
+/*!
+ * \brief Improve quality of the created mesh elements
+ */
+//================================================================================
+
+bool _ViscousBuilder2D::improve()
+{
+  if ( !_proxyMesh )
+    return false;
+
+  // fixed nodes on EDGE's
+  std::set<const SMDS_MeshNode*> fixedNodes;
+  for ( size_t iWire = 0; iWire < _faceSideVec.size(); ++iWire )
+  {
+    StdMeshers_FaceSidePtr      wire = _faceSideVec[ iWire ];
+    const vector<UVPtStruct>& points = wire->GetUVPtStruct();
+    for ( size_t i = 0; i < points.size(); ++i )
+      fixedNodes.insert( fixedNodes.end(), points[i].node );
+  }
+  // fixed proxy nodes
+  for ( size_t iL = 0; iL < _polyLineVec.size(); ++iL )
+  {
+    _PolyLine&         L = _polyLineVec[ iL ];
+    const TopoDS_Edge& E = L._wire->Edge( L._edgeInd );
+    if ( const SMESH_ProxyMesh::SubMesh* sm = _proxyMesh->GetProxySubMesh( E ))
+    {
+      const UVPtStructVec& points = sm->GetUVPtStructVec();
+      for ( size_t i = 0; i < points.size(); ++i )
+        fixedNodes.insert( fixedNodes.end(), points[i].node );
+    }
+    for ( size_t i = 0; i < L._rightNodes.size(); ++i )
+      fixedNodes.insert( fixedNodes.end(), L._rightNodes[i] );
+  }
+
+  // smoothing
+  SMESH_MeshEditor editor( _mesh );
+  for ( size_t iL = 0; iL < _polyLineVec.size(); ++iL )
+  {
+    _PolyLine& L = _polyLineVec[ iL ];
+    if ( L._isStraight2D ) continue;
+    // SMESH_MeshEditor::SmoothMethod how =
+    //   L._isStraight2D ? SMESH_MeshEditor::LAPLACIAN : SMESH_MeshEditor::CENTROIDAL;
+    //editor.Smooth( L._newFaces, fixedNodes, how, /*nbIt = */3 );
+    //editor.Smooth( L._newFaces, fixedNodes, SMESH_MeshEditor::LAPLACIAN, /*nbIt = */1 );
+    editor.Smooth( L._newFaces, fixedNodes, SMESH_MeshEditor::CENTROIDAL, /*nbIt = */3 );
+  }
+  return true;
+}
+
 //================================================================================
 /*!
  * \brief Remove elements and nodes from a face