Salome HOME
Fix compilation errors using gcc-5.X relating to explicit stream::operator bool()
[modules/smesh.git] / src / StdMeshers / StdMeshers_Adaptive1D.cxx
index 38b703539bc074e4fa5b2e4b8228e5b54430f64b..8ffcc0ce255b21818cd42c9fc905a41763123887 100644 (file)
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2013  CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2015  CEA/DEN, EDF R&D, OPEN CASCADE
 //
 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
@@ -6,7 +6,7 @@
 // This library is free software; you can redistribute it and/or
 // modify it under the terms of the GNU Lesser General Public
 // License as published by the Free Software Foundation; either
-// version 2.1 of the License.
+// version 2.1 of the License, or (at your option) any later version.
 //
 // This library is distributed in the hope that it will be useful,
 // but WITHOUT ANY WARRANTY; without even the implied warranty of
 
 #include <BRepAdaptor_Curve.hxx>
 #include <BRepAdaptor_Surface.hxx>
+#include <BRepBndLib.hxx>
 #include <BRepMesh_IncrementalMesh.hxx>
 #include <BRep_Tool.hxx>
 #include <Bnd_B3d.hxx>
+#include <Bnd_Box.hxx>
 #include <GCPnts_AbscissaPoint.hxx>
 #include <GeomAdaptor_Curve.hxx>
 #include <Geom_Curve.hxx>
 #include <Poly_Array1OfTriangle.hxx>
+#include <Poly_PolygonOnTriangulation.hxx>
 #include <Poly_Triangulation.hxx>
 #include <TColgp_Array1OfPnt.hxx>
 #include <TopExp.hxx>
@@ -167,6 +170,37 @@ namespace // internal utils
     SegSizeTree*                 mySizeTree;
   };
 
+  //================================================================================
+  /*!
+   * \brief Segment of Poly_PolygonOnTriangulation
+   */
+  struct Segment
+  {
+    gp_XYZ myPos, myDir;
+    double myLength;
+
+    void Init( const gp_Pnt& p1, const gp_Pnt& p2 )
+    {
+      myPos    = p1.XYZ();
+      myDir    = p2.XYZ() - p1.XYZ();
+      myLength = myDir.Modulus();
+      if ( myLength > std::numeric_limits<double>::min() )
+        myDir /= myLength;
+    }
+    bool Distance( const gp_Pnt& P, double& dist ) const // returns length of normal projection
+    {
+      gp_XYZ p = P.XYZ();
+      p.Subtract( myPos );
+      double proj = p.Dot( myDir );
+      if ( 0 < proj && proj < myLength )
+      {
+        p.Cross( myDir );
+        dist = p.Modulus();
+        return true;
+      }
+      return false;
+    }
+  };
   //================================================================================
   /*!
    * \brief Data of triangle used to locate it in an octree and to find distance
@@ -174,14 +208,17 @@ namespace // internal utils
    */
   struct Triangle
   {
-    Bnd_B3d myBox;
-    bool    myIsChecked; // to mark treated trias instead of using std::set
+    Bnd_B3d  myBox;
+    bool     myIsChecked; // to mark treated trias instead of using std::set
+    bool     myHasNodeOnVertex;
+    Segment* mySegments[3];
     // data for DistToProjection()
-    gp_XYZ myN0, myEdge1, myEdge2, myNorm, myPVec;
-    double myInvDet, myMaxSize2;
+    gp_XYZ   myN0, myEdge1, myEdge2, myNorm, myPVec;
+    double   myInvDet, myMaxSize2;
 
     void Init( const gp_Pnt& n1, const gp_Pnt& n2, const gp_Pnt& n3 ); 
     bool DistToProjection( const gp_Pnt& p, double& dist ) const;
+    bool DistToSegment   ( const gp_Pnt& p, double& dist ) const;
   };
   //================================================================================
   /*!
@@ -197,6 +234,7 @@ namespace // internal utils
   struct TriaTreeData : public ElemTreeData
   {
     vector< Triangle >           myTrias;
+    vector< Segment >            mySegments;
     double                       myFaceTol;
     double                       myTriasDeflection;
     BBox                         myBBox;
@@ -243,25 +281,25 @@ namespace // internal utils
   };
   //================================================================================
   /*!
-   * \brief BRepMesh_IncrementalMesh with access to its protected Bnd_Box
+   * \brief Link of two nodes
    */
-  struct IncrementalMesh : public BRepMesh_IncrementalMesh
+  struct NLink : public std::pair< int, int >
   {
-    IncrementalMesh(const TopoDS_Shape& shape,
-                    const Standard_Real deflection,
-                    const bool          relative):
-      BRepMesh_IncrementalMesh( shape, deflection, relative )
+    NLink( int n1, int n2 )
     {
+      if ( n1 < n2 )
+      {
+        first  = n1;
+        second = n2;
+      }
+      else
+      {
+        first  = n2;
+        second = n1;
+      }
     }
-    Bnd_B3d GetBox() const
-    {
-      Standard_Real TXmin, TYmin, TZmin, TXmax, TYmax, TZmax;
-      myBox.Get(TXmin, TYmin, TZmin, TXmax, TYmax, TZmax);
-      Bnd_B3d bb;
-      bb.Add( gp_XYZ( TXmin, TYmin, TZmin ));
-      bb.Add( gp_XYZ( TXmax, TYmax, TZmax ));
-      return bb;
-    }
+    int N1() const { return first; }
+    int N2() const { return second; }
   };
 
   //================================================================================
@@ -295,33 +333,82 @@ namespace // internal utils
       }
       for ( int i = myNodes->Lower(); i <= myNodes->Upper(); ++i )
         myBBox.Add( myNodes->Value(i).XYZ() );
-
     }
   }
+  //================================================================================
+  /*!
+   * \brief Prepare data for search of trinagles in GetMinDistInSphere()
+   */
+  //================================================================================
+
   void TriaTreeData::PrepareToTriaSearch()
   {
     if ( !myTrias.empty() ) return; // already done
     if ( !myPolyTrias ) return;
 
+    // get all boundary links and nodes on VERTEXes
+    map< NLink, Segment* > linkToSegMap;
+    map< NLink, Segment* >::iterator l2s;
+    set< int > vertexNodes;
+    TopLoc_Location loc;
+    Handle(Poly_Triangulation) tr = BRep_Tool::Triangulation( mySurface.Face(), loc );
+    if ( !tr.IsNull() )
+    {
+      TopTools_IndexedMapOfShape edgeMap;
+      TopExp::MapShapes( mySurface.Face(), TopAbs_EDGE, edgeMap );
+      for ( int iE = 1; iE <= edgeMap.Extent(); ++iE )
+      {
+        const TopoDS_Edge& edge = TopoDS::Edge( edgeMap( iE ));
+        Handle(Poly_PolygonOnTriangulation) polygon =
+          BRep_Tool::PolygonOnTriangulation( edge, tr, loc );
+        if ( polygon.IsNull()  )
+          continue;
+        const TColStd_Array1OfInteger& nodes = polygon->Nodes();
+        for ( int i = nodes.Lower(); i < nodes.Upper(); ++i )
+          linkToSegMap.insert( make_pair( NLink( nodes(i), nodes(i+1)), (Segment*)0 ));
+        vertexNodes.insert( nodes( nodes.Lower()));
+        vertexNodes.insert( nodes( nodes.Upper()));
+      }
+      // fill mySegments by boundary links
+      mySegments.resize( linkToSegMap.size() );
+      int iS = 0;
+      for ( l2s = linkToSegMap.begin(); l2s != linkToSegMap.end(); ++l2s, ++iS )
+      {
+        const NLink& link = (*l2s).first;
+        (*l2s).second = & mySegments[ iS ];
+        mySegments[ iS ].Init( myNodes->Value( link.N1() ),
+                               myNodes->Value( link.N2() ));
+      }
+    }
+
+    // initialize myTrias
     myTrias.resize( myPolyTrias->Length() );
     Standard_Integer n1,n2,n3;
     for ( int i = 1; i <= myPolyTrias->Upper(); ++i )
     {
+      Triangle & t = myTrias[ i-1 ];
       myPolyTrias->Value( i ).Get( n1,n2,n3 );
-      myTrias[ i-1 ].Init( myNodes->Value( n1 ),
-                           myNodes->Value( n2 ),
-                           myNodes->Value( n3 ));
+      t.Init( myNodes->Value( n1 ),
+              myNodes->Value( n2 ),
+              myNodes->Value( n3 ));
+      int nbSeg = 0;
+      if (( l2s = linkToSegMap.find( NLink( n1, n2 ))) != linkToSegMap.end())
+        t.mySegments[ nbSeg++ ] = l2s->second;
+      if (( l2s = linkToSegMap.find( NLink( n2, n3 ))) != linkToSegMap.end())
+        t.mySegments[ nbSeg++ ] = l2s->second;
+      if (( l2s = linkToSegMap.find( NLink( n3, n1 ))) != linkToSegMap.end())
+        t.mySegments[ nbSeg++ ] = l2s->second;
+      while ( nbSeg < 3 )
+        t.mySegments[ nbSeg++ ] = NULL;
+
+      t.myIsChecked = false;
+      t.myHasNodeOnVertex = ( vertexNodes.count( n1 ) ||
+                              vertexNodes.count( n2 ) ||
+                              vertexNodes.count( n3 ));
     }
-    myTree->FillIn();
 
-    // TODO: mark triangles with nodes on VERTEXes to
-    // less frequently compare with avoidPnt in GetMinDistInSphere()
-    //
-    // Handle(Poly_PolygonOnTriangulation) polygon =
-    //   BRep_Tool::PolygonOnTriangulation( edge, tr, loc );
-    // if ( polygon.IsNull() || !pologon.HasParameters() )
-    //   continue;
-    // Handle(TColStd_Array1OfInteger) nodeIDs = polygon->Nodes();
+    // fill the tree of triangles
+    myTree->FillIn();
   }
 
   //================================================================================
@@ -347,10 +434,8 @@ namespace // internal utils
       isConstSize = false;
     }
 
-    typedef std::pair<int,int> TLink;
-    TLink link;
-    map< TLink, double >           lenOfDoneLink;
-    map< TLink, double >::iterator link2len;
+    map< NLink, double >           lenOfDoneLink;
+    map< NLink, double >::iterator link2len;
 
     Standard_Integer n[4];
     gp_Pnt p[4];
@@ -373,11 +458,7 @@ namespace // internal utils
       maxLinkLen = 0;
       for ( int j = 0; j < 3; ++j )
       {
-        if ( n[j] < n[j+1] )
-          link = TLink( n[j], n[j+1] );
-        else
-          link = TLink( n[j+1], n[j] );
-        link2len  = lenOfDoneLink.insert( make_pair( link, -1. )).first;
+        link2len  = lenOfDoneLink.insert( make_pair( NLink( n[j], n[j+1] ), -1. )).first;
         isDone[j] = !((*link2len).second < 0 );
         a[j]      = isDone[j] ? (*link2len).second : (*link2len).second = p[j].Distance( p[j+1] );
         if ( isDone[j] )
@@ -399,6 +480,8 @@ namespace // internal utils
       if ( !isDone[ jLongest ] || !isConstSize )
       {
         //++nbLinks;
+        if ( size < numeric_limits<double>::min() )
+          continue;
         int nb = Max( 1, int( maxLinkLen / size / 2 ));
         for ( int k = 0; k <= nb; ++k )
         {
@@ -428,10 +511,13 @@ namespace // internal utils
   {
     double minDist2 = 1e100;
     const double tol2 = myFaceTol * myFaceTol;
+    const double dMin2 = myTriasDeflection * myTriasDeflection;
 
     TriaTreeData* me = const_cast<TriaTreeData*>( this );
     me->myFoundTriaIDs.clear();
     myTree->GetElementsInSphere( p.XYZ(), radius, me->myFoundTriaIDs );
+    if ( myFoundTriaIDs.empty() )
+      return minDist2;
 
     Standard_Integer n[ 3 ];
     for ( size_t i = 0; i < myFoundTriaIDs.size(); ++i )
@@ -442,19 +528,35 @@ namespace // internal utils
       t.myIsChecked = true;
 
       double d, minD2 = minDist2;
-      bool avoidTria = false;
       myPolyTrias->Value( myFoundTriaIDs[i]+1 ).Get( n[0],n[1],n[2] );
-      for ( int i = 0; i < 3; ++i )
+      if ( avoidPnt && t.myHasNodeOnVertex )
       {
-        const gp_Pnt& pn = myNodes->Value(n[i]);
-        if ( avoidTria = ( avoidPnt && pn.SquareDistance(*avoidPnt) <= tol2 ))
-          break;
-        if ( !projectedOnly )
-          minD2 = Min( minD2, pn.SquareDistance( p ));
+        bool avoidTria = false;
+        for ( int i = 0; i < 3; ++i )
+        {
+          const gp_Pnt& pn = myNodes->Value(n[i]);
+          if (( avoidTria = ( pn.SquareDistance( *avoidPnt ) <= tol2 )))
+            break;
+          if ( !projectedOnly )
+            minD2 = Min( minD2, pn.SquareDistance( p ));
+        }
+        if ( avoidTria )
+          continue;
+        if (( projectedOnly || minD2 < t.myMaxSize2 ) &&
+            ( t.DistToProjection( p, d ) || t.DistToSegment( p, d )))
+          minD2 = Min( minD2, d*d );
+        minDist2 = Min( minDist2, minD2 );
       }
-      if ( !avoidTria )
+      else if ( projectedOnly )
       {
-        if ( minD2 < t.myMaxSize2 && t.DistToProjection( p, d ))
+        if ( t.DistToProjection( p, d ) && d*d > dMin2 )
+          minDist2 = Min( minDist2, d*d );
+      }
+      else
+      {
+        for ( int i = 0; i < 3; ++i )
+          minD2 = Min( minD2, p.SquareDistance( myNodes->Value(n[i]) ));
+        if ( minD2 < t.myMaxSize2  && ( t.DistToProjection( p, d ) || t.DistToSegment( p, d )))
           minD2 = Min( minD2, d*d );
         minDist2 = Min( minDist2, minD2 );
       }
@@ -524,6 +626,31 @@ namespace // internal utils
     return true;
   }
 
+  //================================================================================
+  /*!
+   * \brief Compute distance from a point to either of mySegments. Return false if the point
+   *        is not projected on a segment
+   */
+  //================================================================================
+
+  bool Triangle::DistToSegment( const gp_Pnt& p, double& dist ) const
+  {
+    dist = 1e100;
+    bool res = false;
+    double d;
+    for ( int i = 0; i < 3; ++i )
+    {
+      if ( !mySegments[ i ])
+        break;
+      if ( mySegments[ i ]->Distance( p, d ))
+      {
+        res = true;
+        dist = Min( dist, d );
+      }
+    }
+    return res;
+  }
+
   //================================================================================
   /*!
    * \brief Consturct ElementBndBoxTree of Poly_Triangulation of a FACE
@@ -576,7 +703,7 @@ namespace // internal utils
   void ElementBndBoxTree::buildChildrenData()
   {
     ElemTreeData* data = GetElemData();
-    for ( int i = 0; i < _elementIDs.size(); ++i )
+    for ( size_t i = 0; i < _elementIDs.size(); ++i )
     {
       const Bnd_B3d* elemBox = data->GetBox( _elementIDs[i] );
       for (int j = 0; j < 8; j++)
@@ -591,7 +718,7 @@ namespace // internal utils
     {
       ElementBndBoxTree* child = static_cast<ElementBndBoxTree*>( myChildren[j] );
       child->_elementIDs = data->myWorkIDs[ j ];
-      if ( child->_elementIDs.size() <= theMaxNbElemsInLeaf )
+      if ((int) child->_elementIDs.size() <= theMaxNbElemsInLeaf )
         child->myIsLeaf = true;
       data->myWorkIDs[ j ].clear();
     }
@@ -614,7 +741,7 @@ namespace // internal utils
       if ( isLeaf() )
       {
         ElemTreeData* data = GetElemData();
-        for ( int i = 0; i < _elementIDs.size(); ++i )
+        for ( size_t i = 0; i < _elementIDs.size(); ++i )
           if ( !data->GetBox( _elementIDs[i] )->IsOut( center, radius ))
             foundElemIDs.push_back( _elementIDs[i] );
       }
@@ -865,7 +992,7 @@ ostream & StdMeshers_Adaptive1D::SaveTo(ostream & save)
 istream & StdMeshers_Adaptive1D::LoadFrom(istream & load)
 {
   int dummyParam;
-  bool isOK = (load >> myMinSize >> myMaxSize >> myDeflection >> dummyParam >> dummyParam);
+  bool isOK = static_cast<bool>(load >> myMinSize >> myMaxSize >> myDeflection >> dummyParam >> dummyParam);
   if (!isOK)
     load.clear(ios::badbit | load.rdstate());
   return load;
@@ -932,9 +1059,9 @@ bool StdMeshers_Adaptive1D::SetParametersByMesh(const SMESH_Mesh*   theMesh,
 bool StdMeshers_Adaptive1D::SetParametersByDefaults(const TDefaults&  dflts,
                                                     const SMESH_Mesh* /*theMesh*/)
 {
-  myMinSize = dflts._elemLength / 100;
+  myMinSize = dflts._elemLength / 10;
   myMaxSize = dflts._elemLength * 2;
-  myDeflection = myMinSize / 10;
+  myDeflection = myMinSize / 7;
   return true;
 }
 
@@ -991,7 +1118,7 @@ void AdaptiveAlgo::SetHypothesis( const StdMeshers_Adaptive1D* hyp )
 bool AdaptiveAlgo::Compute(SMESH_Mesh &         theMesh,
                            const TopoDS_Shape & theShape)
 {
-  //*theProgress = 0.01;
+  // *theProgress = 0.01;
 
   if ( myHyp->GetMinSize() > myHyp->GetMaxSize() )
     return error( "Bad parameters: min size > max size" );
@@ -1005,12 +1132,21 @@ bool AdaptiveAlgo::Compute(SMESH_Mesh &         theMesh,
   TopExp::MapShapes( theMesh.GetShapeToMesh(), TopAbs_FACE, faceMap );
 
   // Triangulate the shape with the given deflection ?????????
+  {
+    BRepMesh_IncrementalMesh im( theMesh.GetShapeToMesh(), myHyp->GetDeflection(), /*isRelatif=*/0);
+  }
+
+  // get a bnd box
   Bnd_B3d box;
   {
-    IncrementalMesh im( theMesh.GetShapeToMesh(), myHyp->GetDeflection(), /*Relatif=*/false);
-    box = im.GetBox();
+    Bnd_Box aBox;
+    BRepBndLib::Add( theMesh.GetShapeToMesh(), aBox);
+    Standard_Real TXmin, TYmin, TZmin, TXmax, TYmax, TZmax;
+    aBox.Get(TXmin, TYmin, TZmin, TXmax, TYmax, TZmax);
+    box.Add( gp_XYZ( TXmin, TYmin, TZmin ));
+    box.Add( gp_XYZ( TXmax, TYmax, TZmax ));
   }
-  //*theProgress = 0.3;
+  // *theProgress = 0.3;
 
   // holder of segment size at each point
   SegSizeTree sizeTree( box, grading, myHyp->GetMinSize(), myHyp->GetMaxSize() );
@@ -1043,6 +1179,7 @@ bool AdaptiveAlgo::Compute(SMESH_Mesh &         theMesh,
       eData.AddPoint( eData.myPoints.end(), eData.myC3d.LastParameter() );
     }
   }
+  if ( myEdges.empty() ) return true;
   if ( _computeCanceled ) return false;
 
   // Take into account size of already existing segments
@@ -1060,7 +1197,7 @@ bool AdaptiveAlgo::Compute(SMESH_Mesh &         theMesh,
   StdMeshers_Regular_1D::_value[ DEFLECTION_IND ] = myHyp->GetDeflection();
 
   list< double > params;
-  for ( int iE = 0; iE < myEdges.size(); ++iE )
+  for ( size_t iE = 0; iE < myEdges.size(); ++iE )
   {
     EdgeData& eData = myEdges[ iE ];
     //cout << "E " << theMesh.GetMeshDS()->ShapeToIndex( eData.Edge() ) << endl;
@@ -1106,7 +1243,7 @@ bool AdaptiveAlgo::Compute(SMESH_Mesh &         theMesh,
 
     triaSearcher->SetSizeByTrias( sizeTree, myHyp->GetDeflection() );
 
-    for ( int iE = 0; iE < myEdges.size(); ++iE )
+    for ( size_t iE = 0; iE < myEdges.size(); ++iE )
     {
       EdgeData& eData = myEdges[ iE ];
 
@@ -1153,31 +1290,34 @@ bool AdaptiveAlgo::Compute(SMESH_Mesh &         theMesh,
         //cout << "E " << theMesh.GetMeshDS()->ShapeToIndex( eData.Edge() ) << endl;
         sizeDecreased = false;
         const gp_Pnt* avoidPnt = & eData.First().myP;
+        EdgeData::TPntIter pItLast = --eData.myPoints.end(), pItFirst = eData.myPoints.begin();
         for ( pIt1 = eData.myPoints.begin(); pIt1 != eData.myPoints.end();  )
         {
           double distToFace =
             triaSearcher->GetMinDistInSphere( pIt1->myP, pIt1->mySegSize, isAdjFace, avoidPnt );
           double allowedSize = Max( minSize, distToFace*( 1. + grading ));
-          if ( 1.1 * allowedSize < pIt1->mySegSize  )
+          if ( allowedSize < pIt1->mySegSize  )
           {
-            sizeDecreased = true;
-            sizeTree.SetSize( pIt1->myP, allowedSize );
             // cout << "E " << theMesh.GetMeshDS()->ShapeToIndex( eData.Edge() )
-            //      << "\t SetSize " << allowedSize << " at "
-            //      << pIt1->myP.X() <<", "<< pIt1->myP.Y()<<", "<<pIt1->myP.Z() << endl;
-            pIt2 = pIt1;
-            if ( --pIt2 != eData.myPoints.end() && pIt2->mySegSize > allowedSize )
-              sizeTree.SetSize( eData.myC3d.Value( 0.6*pIt2->myU + 0.4*pIt1->myU ), allowedSize );
-            pIt2 = pIt1;
-            if ( ++pIt2 != eData.myPoints.end() && pIt2->mySegSize > allowedSize )
-              sizeTree.SetSize( eData.myC3d.Value( 0.6*pIt2->myU + 0.4*pIt1->myU ), allowedSize );
+            //      << "\t closure detected " << endl;
+            if ( 1.1 * allowedSize < pIt1->mySegSize  )
+            {
+              sizeDecreased = true;
+              sizeTree.SetSize( pIt1->myP, allowedSize );
+              // cout << "E " << theMesh.GetMeshDS()->ShapeToIndex( eData.Edge() )
+              //      << "\t SetSize " << allowedSize << " at "
+              //      << pIt1->myP.X() <<", "<< pIt1->myP.Y()<<", "<<pIt1->myP.Z() << endl;
+              pIt2 = pIt1;
+              if ( pIt1 != pItFirst && ( --pIt2 )->mySegSize > allowedSize )
+                sizeTree.SetSize( eData.myC3d.Value( 0.6*pIt2->myU + 0.4*pIt1->myU ), allowedSize );
+              pIt2 = pIt1;
+              if ( pIt1 != pItLast  && ( ++pIt2 )->mySegSize > allowedSize )
+                sizeTree.SetSize( eData.myC3d.Value( 0.6*pIt2->myU + 0.4*pIt1->myU ), allowedSize );
+            }
             pIt1->mySegSize = allowedSize;
           }
           ++pIt1;
-          if ( & (*pIt1) == & eData.Last() )
-            avoidPnt = & eData.Last().myP;
-          else
-            avoidPnt = NULL;
+          avoidPnt = ( pIt1 == pItLast ) ? & eData.Last().myP : NULL;
 
           if ( iLoop > 20 )
           {
@@ -1191,7 +1331,7 @@ bool AdaptiveAlgo::Compute(SMESH_Mesh &         theMesh,
       } // while ( sizeDecreased )
     } // loop on myEdges
 
-    //*theProgress = 0.3 + 0.3 * iF / double( faceMap.Extent() );
+    // *theProgress = 0.3 + 0.3 * iF / double( faceMap.Extent() );
 
   } // loop on faceMap
 
@@ -1214,24 +1354,26 @@ bool AdaptiveAlgo::makeSegments()
 
   vector< double > nbSegs, params;
 
-  for ( int iE = 0; iE < myEdges.size(); ++iE )
+  for ( size_t iE = 0; iE < myEdges.size(); ++iE )
   {
     EdgeData& eData = myEdges[ iE ];
 
-    // estimate roughly min segement size on the EDGE
+    // estimate roughly min segment size on the EDGE
     double edgeMinSize = myHyp->GetMaxSize();
     EdgeData::TPntIter pIt1 = eData.myPoints.begin();
     for ( ; pIt1 != eData.myPoints.end(); ++pIt1 )
-      edgeMinSize = Min( edgeMinSize, mySizeTree->GetSize( pIt1->myP ));
+      edgeMinSize = Min( edgeMinSize,
+                         Min( pIt1->mySegSize, mySizeTree->GetSize( pIt1->myP )));
 
-    const double f = eData.myC3d.FirstParameter(), l = eData.myC3d.LastParameter();
+    const double      f = eData.myC3d.FirstParameter(), l = eData.myC3d.LastParameter();
     const double parLen = l - f;
     const int  nbDivSeg = 5;
-    int           nbDiv = Max( 1, int ( eData.myLength / edgeMinSize * nbDivSeg ));
+    size_t        nbDiv = Max( 1, int ( eData.myLength / edgeMinSize * nbDivSeg ));
 
     // compute nb of segments
-    bool toRecompute = true;
+    bool  toRecompute = true;
     double maxSegSize = 0;
+    size_t i = 1, segCount;
     //cout << "E " << theMesh.GetMeshDS()->ShapeToIndex( eData.Edge() ) << endl;
     while ( toRecompute ) // recompute if segment size at some point is less than edgeMinSize/nbDivSeg
     {
@@ -1239,11 +1381,32 @@ bool AdaptiveAlgo::makeSegments()
       nbSegs[0] = 0;
       toRecompute = false;
 
+      // fill nbSegs with segment size stored in EdgeData::ProbePnt::mySegSize which can
+      // be less than size in mySizeTree
+      pIt1 = eData.myPoints.begin();
+      EdgeData::ProbePnt* pp1 = &(*pIt1), *pp2;
+      for ( ++pIt1; pIt1 != eData.myPoints.end(); ++pIt1 )
+      {
+        pp2 = &(*pIt1);
+        double size1 = Min( pp1->mySegSize, myHyp->GetMaxSize() );
+        double size2 = Min( pp2->mySegSize, myHyp->GetMaxSize() );
+        double r, u, du = pp2->myU - pp1->myU;
+        while(( u = f + parLen * i / nbDiv ) < pp2->myU )
+        {
+          r = ( u - pp1->myU ) / du;
+          nbSegs[i] = (1-r) * size1 + r * size2;
+          ++i;
+        }
+        if ( i < nbSegs.size() )
+          nbSegs[i] = size2;
+        pp1 = pp2;
+      }
+      // fill nbSegs with local nb of segments
       gp_Pnt p1 = eData.First().myP, p2, pDiv = p1;
-      for ( size_t i = 1, segCount = 1; i < nbSegs.size(); ++i )
+      for ( i = 1, segCount = 1; i < nbSegs.size(); ++i )
       {
         p2 = eData.myC3d.Value( f + parLen * i / nbDiv );
-        double locSize = Min( mySizeTree->GetSize( p2 ), myHyp->GetMaxSize() );
+        double locSize = Min( mySizeTree->GetSize( p2 ), nbSegs[i] );
         double nb      = p1.Distance( p2 ) / locSize;
         // if ( nbSegs.size() < 30 )
         //   cout << "locSize " << locSize << " nb " << nb << endl;
@@ -1266,13 +1429,13 @@ bool AdaptiveAlgo::makeSegments()
     }
 
     // compute parameters of nodes
-    int nbSegFinal = int(floor(nbSegs.back()+0.5));
+    size_t nbSegFinal = Max( 1, int(floor( nbSegs.back() + 0.5 )));
     double fact = nbSegFinal / nbSegs.back();
     if ( maxSegSize / fact > myHyp->GetMaxSize() )
       fact = ++nbSegFinal / nbSegs.back();
     //cout << "nbSegs.back() " << nbSegs.back() << " nbSegFinal " << nbSegFinal << endl;
     params.clear();
-    for ( int i = 0, segCount = 1; segCount < nbSegFinal; ++segCount )
+    for ( i = 0, segCount = 1; segCount < nbSegFinal; ++segCount )
     {
       while ( nbSegs[i] * fact < segCount )
         ++i;
@@ -1300,7 +1463,7 @@ bool AdaptiveAlgo::makeSegments()
     helper.SetElementsOnShape( true );
     const int ID = 0;
     const SMDS_MeshNode *n1 = nf, *n2;
-    for ( size_t i = 0; i < params.size(); ++i, n1 = n2 )
+    for ( i = 0; i < params.size(); ++i, n1 = n2 )
     {
       gp_Pnt p2 = eData.myC3d.Value( params[i] );
       n2 = helper.AddNode( p2.X(), p2.Y(), p2.Z(), ID, params[i] );
@@ -1339,7 +1502,7 @@ bool AdaptiveAlgo::Evaluate(SMESH_Mesh &         theMesh,
 
   for ( ; edExp.More(); edExp.Next() )
   {
-    const TopoDS_Edge & edge = TopoDS::Edge( edExp.Current() );
+    //const TopoDS_Edge & edge = TopoDS::Edge( edExp.Current() );
     StdMeshers_Regular_1D::Evaluate( theMesh, theShape, theResMap );
   }
   return true;