diff --git a/src/BOPAlgo/BOPAlgo_PaveFiller.cdl b/src/BOPAlgo/BOPAlgo_PaveFiller.cdl
index 63d569406b..9d943abfe1 100644
--- a/src/BOPAlgo/BOPAlgo_PaveFiller.cdl
+++ b/src/BOPAlgo/BOPAlgo_PaveFiller.cdl
@@ -418,7 +418,8 @@ is
     -- state information 
 
     UpdateFaceInfo(me:out;
-        theDME:out DataMapOfPaveBlockListOfPaveBlock from BOPDS) 
+        theDME:out DataMapOfPaveBlockListOfPaveBlock from BOPDS; 
+        theDMV: DataMapOfIntegerInteger from BOPCol) 
       is protected; 
     ---Purpose: 
     -- Updates the information about faces
diff --git a/src/BOPAlgo/BOPAlgo_PaveFiller_6.cxx b/src/BOPAlgo/BOPAlgo_PaveFiller_6.cxx
index 116160be9c..2edebeb5eb 100644
--- a/src/BOPAlgo/BOPAlgo_PaveFiller_6.cxx
+++ b/src/BOPAlgo/BOPAlgo_PaveFiller_6.cxx
@@ -32,11 +32,14 @@
 #include <TopoDS_Vertex.hxx>
 #include <TopoDS_Compound.hxx>
 
+#include <TopExp.hxx>
 #include <TopExp_Explorer.hxx>
 
 #include <BRep_Builder.hxx>
 #include <BRep_Tool.hxx>
 
+#include <BRepBuilderAPI_MakeVertex.hxx>
+
 #include <BRepBndLib.hxx>
 #include <BRepTools.hxx>
 
@@ -48,7 +51,10 @@
 #include <IntTools_SequenceOfPntOn2Faces.hxx>
 #include <IntTools_Curve.hxx>
 #include <IntTools_PntOn2Faces.hxx>
+#include <IntTools_ShrunkRange.hxx>
+#include <IntTools_Context.hxx>
 #include <IntTools_Tools.hxx>
+#include <IntTools_EdgeFace.hxx>
 
 #include <IntSurf_ListOfPntOn2S.hxx>
 #include <IntSurf_PntOn2S.hxx>
@@ -65,9 +71,6 @@
 #include <BOPCol_NCVector.hxx>
 #include <BOPCol_TBB.hxx>
 
-#include <IntTools_Context.hxx>
-#include <IntTools_Tools.hxx>
-
 #include <BOPDS_Interf.hxx>
 #include <BOPDS_Iterator.hxx>
 #include <BOPDS_Curve.hxx>
@@ -85,12 +88,9 @@
 #include <BOPDS_CoupleOfPaveBlocks.hxx>
 #include <BOPDS_FaceInfo.hxx>
 #include <BOPDS_CommonBlock.hxx>
+#include <BOPDS_DataMapOfPaveBlockListOfPaveBlock.hxx>
 
 #include <BOPAlgo_Tools.hxx>
-#include <BRepBuilderAPI_MakeVertex.hxx>
-#include <TopExp.hxx>
-#include <IntTools_ShrunkRange.hxx>
-#include <BOPDS_DataMapOfPaveBlockListOfPaveBlock.hxx>
 
 static void ToleranceFF(const BRepAdaptor_Surface& aBAS1,
                         const BRepAdaptor_Surface& aBAS2,
@@ -597,7 +597,7 @@ void BOPAlgo_PaveFiller::MakeBlocks()
   }
   //
   // update face info
-  UpdateFaceInfo(aDMExEdges);
+  UpdateFaceInfo(aDMExEdges, aDMI);
   //Update all pave blocks
   UpdatePaveBlocks(aDMI);
   //-----------------------------------------------------scope t
@@ -754,7 +754,8 @@ Standard_Integer BOPAlgo_PaveFiller::PostTreatFF
       if (!bHasPaveBlocks) {
         if (bOld) {
           aDMExEdges.ChangeFind(aPB1).Append(aPB1);
-        } else {
+        }
+        else {
           aSI.SetShapeType(aType);
           aSI.SetShape(aSx);
           iE=myDS->Append(aSI);
@@ -894,12 +895,13 @@ Standard_Integer BOPAlgo_PaveFiller::PostTreatFF
 //purpose  : 
 //=======================================================================
 void BOPAlgo_PaveFiller::UpdateFaceInfo
-  (BOPDS_DataMapOfPaveBlockListOfPaveBlock& theDME)
+  (BOPDS_DataMapOfPaveBlockListOfPaveBlock& theDME,
+   const BOPCol_DataMapOfIntegerInteger& theDMV)
 {
   Standard_Integer i, j, nV1, nF1, nF2, 
-                   aNbFF, aNbC, aNbP, aNbS, aNbPBIn;
-  BOPDS_IndexedMapOfPaveBlock aMPBCopy;
+                   aNbFF, aNbC, aNbP;
   BOPDS_ListIteratorOfListOfPaveBlock aItLPB;
+  BOPCol_MapOfInteger aMF;
   //
   BOPDS_VectorOfInterfFF& aFFs=myDS->InterfFF();
   aNbFF=aFFs.Extent();
@@ -912,6 +914,7 @@ void BOPAlgo_PaveFiller::UpdateFaceInfo
     BOPDS_FaceInfo& aFI1=myDS->ChangeFaceInfo(nF1);
     BOPDS_FaceInfo& aFI2=myDS->ChangeFaceInfo(nF2);
     //
+    // 1.1. Section edges
     BOPDS_VectorOfCurve& aVNC=aFF.ChangeCurves();
     aNbC=aVNC.Extent();
     for (j=0; j<aNbC; ++j) {
@@ -919,21 +922,24 @@ void BOPAlgo_PaveFiller::UpdateFaceInfo
       BOPDS_ListOfPaveBlock& aLPBC=aNC.ChangePaveBlocks();
       aItLPB.Initialize(aLPBC);
       //
+      // Treat existing pave blocks
       if (aItLPB.More() && theDME.IsBound(aLPBC.First())) {
         const Handle(BOPDS_PaveBlock)& aPB=aLPBC.First();
-        BOPDS_ListOfPaveBlock& aLPB = theDME.ChangeFind(aPB);
+        BOPDS_ListOfPaveBlock& aLPB=theDME.ChangeFind(aPB);
         UpdateExistingPaveBlocks(aPB, aLPB, nF1, nF2);
         aLPBC.Clear();
         continue;
       }
       //
-      for(; aItLPB.More(); aItLPB.Next()) {
+      // Add section edges to face info
+      for (; aItLPB.More(); aItLPB.Next()) {
         const Handle(BOPDS_PaveBlock)& aPB=aItLPB.Value();
         aFI1.ChangePaveBlocksSc().Add(aPB);
         aFI2.ChangePaveBlocksSc().Add(aPB);
       }
     }
-    // VerticesSc
+    //
+    // 1.2. Section vertices
     const BOPDS_VectorOfPoint& aVNP=aFF.Points();
     aNbP=aVNP.Extent();
     for (j=0; j<aNbP; ++j) {
@@ -945,43 +951,81 @@ void BOPAlgo_PaveFiller::UpdateFaceInfo
       aFI1.ChangeVerticesSc().Add(nV1);
       aFI2.ChangeVerticesSc().Add(nV1);
     }
+    //
+    aMF.Add(nF1);
+    aMF.Add(nF2);
   }
   //
-  //2. PaveBlocksIn
-  if (theDME.IsEmpty()) {
+  Standard_Boolean bVerts, bEdges;
+  //
+  bVerts = theDMV.Extent() > 0;
+  bEdges = theDME.Extent() > 0;
+  //
+  if (!bVerts && !bEdges) {
     return;
   }
   //
-  aNbS=myDS->NbSourceShapes();
-  for (i=0; i<aNbS; ++i) {
-    const BOPDS_ShapeInfo& aSI=myDS->ShapeInfo(i);
-    if (aSI.ShapeType()!=TopAbs_FACE) {
-      continue;
-    }
-    if(!myDS->HasFaceInfo(i)) {
-      continue;
-    }
-    BOPDS_FaceInfo& aFI=myDS->ChangeFaceInfo(i);
+  // 2. Update Face Info information with new vertices and new
+  //    pave blocks created in PostTreatFF from existing ones
+  Standard_Integer nV2, aNbPB;
+  BOPCol_MapIteratorOfMapOfInteger aItMF;
+  BOPCol_DataMapIteratorOfDataMapOfIntegerInteger aItMV;
+  //
+  aItMF.Initialize(aMF);
+  for (; aItMF.More(); aItMF.Next()) {
+    nF1 = aItMF.Value();
     //
-    BOPDS_IndexedMapOfPaveBlock& aMPBIn=aFI.ChangePaveBlocksIn();
-    aMPBCopy.Assign(aMPBIn);
-    aMPBIn.Clear();
+    BOPDS_FaceInfo& aFI = myDS->ChangeFaceInfo(nF1);
     //
-    aNbPBIn=aMPBCopy.Extent();
-    for (j=1; j<=aNbPBIn; ++j) {
-      const Handle(BOPDS_PaveBlock)& aPB = aMPBCopy(j);
-      if (theDME.IsBound(aPB)) {
-        const BOPDS_ListOfPaveBlock& aLPB = theDME.Find(aPB);
-        aItLPB.Initialize(aLPB);
-        for (; aItLPB.More(); aItLPB.Next()) {
-          const Handle(BOPDS_PaveBlock)& aPB1 = aItLPB.Value();
-          aMPBIn.Add(aPB1);
+    // 2.1. Update information about vertices
+    if (bVerts) {
+      BOPCol_MapOfInteger& aMVOn = aFI.ChangeVerticesOn();
+      BOPCol_MapOfInteger& aMVIn = aFI.ChangeVerticesIn();
+      //
+      aItMV.Initialize(theDMV);
+      for (; aItMV.More(); aItMV.Next()) {
+        nV1 = aItMV.Key();
+        nV2 = aItMV.Value();
+        //
+        if (aMVOn.Remove(nV1)) {
+          aMVOn.Add(nV2);
         }
-      } else {
-        aMPBIn.Add(aPB);
-      }
-    }//for (j=1; j<=aNbPBIn; ++j) {
-  }//for (i=0; i<aNbS; ++i) {
+        //
+        if (aMVIn.Remove(nV1)) {
+          aMVIn.Add(nV2);
+        }
+      } // for (; aItMV.More(); aItMV.Next()) {
+    } // if (bVerts) {
+    //
+    // 2.2. Update information about pave blocks
+    if (bEdges) {
+      BOPDS_IndexedMapOfPaveBlock& aMPBOn = aFI.ChangePaveBlocksOn();
+      BOPDS_IndexedMapOfPaveBlock& aMPBIn = aFI.ChangePaveBlocksIn();
+      //
+      BOPDS_IndexedMapOfPaveBlock aMPBCopy;
+      for (i = 0; i < 2; ++i) {
+        BOPDS_IndexedMapOfPaveBlock& aMPBOnIn = !i ? aMPBOn : aMPBIn;
+        aMPBCopy = aMPBOnIn;
+        aMPBOnIn.Clear();
+        //
+        aNbPB = aMPBCopy.Extent();
+        for (j = 1; j <= aNbPB; ++j) {
+          const Handle(BOPDS_PaveBlock)& aPB = aMPBCopy(j);
+          if (theDME.IsBound(aPB)) {
+            const BOPDS_ListOfPaveBlock& aLPB = theDME.Find(aPB);
+            aItLPB.Initialize(aLPB);
+            for (; aItLPB.More(); aItLPB.Next()) {
+              const Handle(BOPDS_PaveBlock)& aPB1 = aItLPB.Value();
+              aMPBOnIn.Add(aPB1);
+            }
+          }
+          else {
+            aMPBOnIn.Add(aPB);
+          }
+        } // for (j = 1; j <= aNbPB; ++j) {
+      } // for (i = 0; i < 2; ++i) {
+    } // if (bEdges) {
+  }
 }
 //=======================================================================
 //function : IsExistingVertex
@@ -1807,9 +1851,8 @@ void BOPAlgo_PaveFiller::RemoveUsedVertices(BOPDS_Curve& aNC,
   Handle(BOPDS_PaveBlock) aPB, aPB1, aPB2, aPB2n;
   Handle(BOPDS_CommonBlock) aCB;
   BOPDS_ListIteratorOfListOfPaveBlock aIt, aIt1, aIt2;
-  BOPDS_IndexedMapOfPaveBlock aMPB;
   //
-  //remove micro edges from aLPB
+  // 1. Remove micro edges from aLPB
   aIt.Initialize(aLPB);
   for (; aIt.More();) {
     aPB = aIt.Value();
@@ -1824,10 +1867,6 @@ void BOPAlgo_PaveFiller::RemoveUsedVertices(BOPDS_Curve& aNC,
   if (!aLPB.Extent()) {
     return;
   }
-  //update face info
-  myDS->UpdateFaceInfoOn(nF1);
-  //
-  myDS->UpdateFaceInfoOn(nF2);
   //
   BOPDS_FaceInfo& aFI1 = myDS->ChangeFaceInfo(nF1);
   BOPDS_FaceInfo& aFI2 = myDS->ChangeFaceInfo(nF2);
@@ -1837,7 +1876,7 @@ void BOPAlgo_PaveFiller::RemoveUsedVertices(BOPDS_Curve& aNC,
   BOPDS_IndexedMapOfPaveBlock& aMPBOn2 = aFI2.ChangePaveBlocksOn();
   BOPDS_IndexedMapOfPaveBlock& aMPBIn2 = aFI2.ChangePaveBlocksIn();
   //
-  // remove old pave blocks
+  // 2. Remove old pave blocks
   const Handle(BOPDS_CommonBlock)& aCB1 = myDS->CommonBlock(aPBf);
   bCB = !aCB1.IsNull();
   BOPDS_ListOfPaveBlock aLPB1;
@@ -1863,8 +1902,9 @@ void BOPAlgo_PaveFiller::RemoveUsedVertices(BOPDS_Curve& aNC,
     }
   }
   //
+  // 3. Update pave blocks
   if (bCB) {
-    //create new pave blocks
+    //create new common blocks
     const BOPCol_ListOfInteger& aFaces = aCB1->Faces();
     aIt.Initialize(aLPB);
     for (; aIt.More(); aIt.Next()) {
@@ -1891,55 +1931,66 @@ void BOPAlgo_PaveFiller::RemoveUsedVertices(BOPDS_Curve& aNC,
       aPB=aCB->PaveBlocks().First();
     }
   } 
-  //
-  aIt.Initialize(aLPB);
-  for (; aIt.More(); aIt.Next()) {
-    Handle(BOPDS_PaveBlock)& aPB = aIt.ChangeValue();
-    nE = aPB->OriginalEdge();
-    //
-    Standard_Integer nF = (aMPBOn1.Contains(aPBf) || 
-                           aMPBIn1.Contains(aPBf)) ? nF2 : nF1;
-    const TopoDS_Face& aF = *(TopoDS_Face*)&myDS->Shape(nF);
-    IntTools_Range aShrR(aPB->Pave1().Parameter(), 
-                         aPB->Pave2().Parameter());
-    const TopoDS_Edge& aE = *(TopoDS_Edge*)&myDS->Shape(aPB->Edge());
-    //
-    Standard_Boolean bCom = 
-      BOPTools_AlgoTools::IsBlockInOnFace(aShrR, aF, aE, myContext);
-    if (bCom) {
-      if (bCB) {
-        aCB = myDS->CommonBlock(aPB);
-        aCB->AddFace(nF);
-      } else {
-        aCB = new BOPDS_CommonBlock;
-        aCB->AddPaveBlock(aPB);
-        aCB->AddFace(nF1);
-        aCB->AddFace(nF2);
-        //
-        myDS->SetCommonBlock(aPB, aCB);
-      }
-      aMPB.Add(aPB);
-    }
-    if (!bCB) {
-      myDS->ChangePaveBlocks(nE).Append(aPB);
+  else {
+    nE = aPBf->OriginalEdge();
+    BOPDS_ListOfPaveBlock& aLPBE = myDS->ChangePaveBlocks(nE);
+    aIt.Initialize(aLPB);
+    for (; aIt.More(); aIt.Next()) {
+      aPB = aIt.Value();
+      aLPBE.Append(aPB);
     }
   }
   //
-  Standard_Integer i, aNbPB;
   Standard_Boolean bIn1, bIn2;
   //
   bIn1 = aMPBOn1.Contains(aPBf) || aMPBIn1.Contains(aPBf);
   bIn2 = aMPBOn2.Contains(aPBf) || aMPBIn2.Contains(aPBf);
   //
-  aNbPB=aMPB.Extent();
-  for (i=1; i<=aNbPB; ++i) {
-    aPB = aMPB(i);
-    if (!bIn1) {
-      aMPBIn1.Add(aPB);
-    }
+  if (bIn1 && bIn2) {
+    return;
+  }
+  //
+  // 4. Check new pave blocks for coincidence 
+  //    with the opposite face.
+  //    In case of coincidence create common blocks
+  Standard_Integer nF;
+  Standard_Real aTolE, aTolF;
+  //
+  nF = bIn1 ? nF2 : nF1;
+  const TopoDS_Face& aF = *(TopoDS_Face*)&myDS->Shape(nF);
+  BOPDS_IndexedMapOfPaveBlock& aMPBIn = bIn1 ? aMPBIn2 : aMPBIn1;
+  aTolF = BRep_Tool::Tolerance(aF);
+  //
+  aIt.Initialize(aLPB);
+  for (; aIt.More(); aIt.Next()) {
+    Handle(BOPDS_PaveBlock)& aPB = aIt.ChangeValue();
+    const TopoDS_Edge& aE = *(TopoDS_Edge*)&myDS->Shape(aPB->Edge());
+    aTolE = BRep_Tool::Tolerance(aE);
     //
-    if (!bIn2) {
-      aMPBIn2.Add(aPB);
+    IntTools_EdgeFace anEF;
+    anEF.SetEdge(aE);
+    anEF.SetFace(aF);
+    anEF.SetTolE(aTolE);
+    anEF.SetTolF(aTolF);
+    anEF.SetRange(aPB->Pave1().Parameter(), aPB->Pave2().Parameter());
+    anEF.SetContext(myContext);
+    anEF.Perform();
+    //
+    const IntTools_SequenceOfCommonPrts& aCPrts=anEF.CommonParts();
+    if (aCPrts.Length() == 1) {
+      Standard_Boolean bCoinc = (aCPrts(1).Type() == TopAbs_EDGE);
+      if (bCoinc) {
+        if (bCB) {
+          aCB = myDS->CommonBlock(aPB);
+        } else {
+          aCB = new BOPDS_CommonBlock;
+          aCB->AddPaveBlock(aPB);
+          myDS->SetCommonBlock(aPB, aCB);
+        }
+        aCB->AddFace(nF);
+        //
+        aMPBIn.Add(aPB);
+      }
     }
   }
 }
diff --git a/src/BOPTest/BOPTest_CheckCommands.cxx b/src/BOPTest/BOPTest_CheckCommands.cxx
index cfd5996c8b..7da70dc77f 100644
--- a/src/BOPTest/BOPTest_CheckCommands.cxx
+++ b/src/BOPTest/BOPTest_CheckCommands.cxx
@@ -904,11 +904,11 @@ Standard_Integer bopargcheck (Draw_Interpretor& di,
           di << "  Cases(" << S2_COnS << ")  Total shapes(" << S2_COnSAll << ")" << "\n";
         else
           di << "\n";
-
-        // warning
+      }
+      // warning
+      if(hasUnknown) {
         di << "\n";
-        if(hasUnknown)
-          di << "WARNING: The unexpected test break occurs!" << "\n";
+        di << "WARNING: The unexpected test break occurs!" << "\n";
       }
     } // full output
   } // has faulties
diff --git a/src/BOPTools/BOPTools_AlgoTools.cdl b/src/BOPTools/BOPTools_AlgoTools.cdl
index 923763a79d..77eebf826f 100644
--- a/src/BOPTools/BOPTools_AlgoTools.cdl
+++ b/src/BOPTools/BOPTools_AlgoTools.cdl
@@ -459,17 +459,6 @@ is
             theSolid:Solid from TopoDS)  
         returns Boolean from Standard;  
     ---Purpose: Returns true if the solid <theSolid> is inverted  
-     
-    ComputeTolerance(myclass;
-        theCurve3D : Curve from Geom;
-        theCurve2D : Curve from Geom2d;
-        theSurf    : Surface from Geom;
-        theMaxDist : out Real from Standard; 
-        theMaxPar  : out Real from Standard)
-    returns Boolean from Standard; 
-    ---Purpose:  
-    -- Computes the max distance between points  
-    -- taken from 3D and 2D curves by the same parameter
 
     ComputeTolerance(myclass;
         theFace : Face from TopoDS;
diff --git a/src/BOPTools/BOPTools_AlgoTools2D.cxx b/src/BOPTools/BOPTools_AlgoTools2D.cxx
index 78c0539dd8..61ad907712 100644
--- a/src/BOPTools/BOPTools_AlgoTools2D.cxx
+++ b/src/BOPTools/BOPTools_AlgoTools2D.cxx
@@ -579,7 +579,7 @@ void BOPTools_AlgoTools2D::MakePCurveOnFace
    Handle(Geom2d_Curve)& aC2D, 
    Standard_Real& TolReached2d)
 {
-  Standard_Real aTolR;
+  Standard_Real aTolR, aT;
   Handle(Geom2d_Curve) aC2DA;
   //
   Handle(Geom_Surface) aS=BRep_Tool::Surface(aF);
@@ -618,6 +618,14 @@ void BOPTools_AlgoTools2D::MakePCurveOnFace
   BOPTools_AlgoTools2D::AdjustPCurveOnFace (aF, aFirst, aLast, 
                                             aC2D, aC2DA);
   aC2D=aC2DA;
+  //
+  // compute the appropriate tolerance for the edge
+  if (IntTools_Tools::ComputeTolerance
+      (aC3D, aC2D, aS, aFirst, aLast, aTolR, aT)) {
+    if (aTolR > TolReached2d) {
+      TolReached2d = aTolR;
+    }
+  }
 }
 
 //=======================================================================
diff --git a/src/BOPTools/BOPTools_AlgoTools_1.cxx b/src/BOPTools/BOPTools_AlgoTools_1.cxx
index 46368c048e..cacb2276ac 100644
--- a/src/BOPTools/BOPTools_AlgoTools_1.cxx
+++ b/src/BOPTools/BOPTools_AlgoTools_1.cxx
@@ -77,8 +77,11 @@
 #include <TopTools_ListOfShape.hxx>
 #include <TopTools_ListIteratorOfListOfShape.hxx>
 //
+#include <IntTools_Tools.hxx>
+//
 #include <BOPCol_NCVector.hxx>
 #include <BOPCol_TBB.hxx>
+#include <BRepLib_CheckCurveOnSurface.hxx>
 
 static 
   void CheckEdge (const TopoDS_Edge& E,
@@ -313,131 +316,6 @@ typedef BOPCol_TBBCnt
 //
 //
 //=======================================================================
-//class    : BOPTools_CheckCurveOnSurface
-//purpose  : it is used to check the curve on the surface
-//=======================================================================
-#include <math_GlobOptMin.hxx>
-#include <math_MultipleVarFunctionWithHessian.hxx>
-#include <math_Matrix.hxx>
-#include <Geom2d_TrimmedCurve.hxx>
-
-class BOPTools_CheckCurveOnSurface : 
-  public math_MultipleVarFunctionWithHessian
-{
- public:
-  BOPTools_CheckCurveOnSurface(BOPTools_CheckCurveOnSurface&);
-  BOPTools_CheckCurveOnSurface(const Handle(Geom_Curve)& theC3D,
-                              const Handle(Geom2d_Curve)& theC2D,
-                              const Handle(Geom_Surface)& theSurf)
-    :
-      my3DCurve(theC3D),
-      my2DCurve(theC2D),
-      mySurf(theSurf)
-  {
-  }
-  //
-  virtual Standard_Integer NbVariables() const {
-    return 1;
-  }
-  //
-  virtual Standard_Boolean Value(const math_Vector& theX,
-                                 Standard_Real& theFVal) {
-    try {
-      const Standard_Real aPar = theX(1);
-      if (!CheckParameter(aPar))
-        return Standard_False;
-      gp_Pnt aP1, aP2;
-      gp_Pnt2d aP2d;
-      my3DCurve->D0(aPar, aP1);
-      my2DCurve->D0(aPar, aP2d);
-      mySurf->D0(aP2d.X(), aP2d.Y(), aP2);
-      //
-      theFVal = -1.0*aP1.SquareDistance(aP2);
-    }
-    catch(Standard_Failure) {
-      return Standard_False;
-    }
-    //
-    return Standard_True;
-  }
-  //
-  virtual Standard_Integer GetStateNumber() {
-    return 0;
-  }
-  //
-  virtual Standard_Boolean Gradient(const math_Vector& theX,
-                                    math_Vector& theGrad) {
-    try {
-      const Standard_Real aPar = theX(1);
-      if (!CheckParameter(aPar))
-        return Standard_False;
-      gp_Pnt aP1, aP2;
-      gp_Vec aDC3D, aDSU, aDSV;
-      gp_Pnt2d aP2d;
-      gp_Vec2d aDC2D;
-      
-      my3DCurve->D1(aPar, aP1, aDC3D);
-      my2DCurve->D1(aPar, aP2d, aDC2D);
-      mySurf->D1(aP2d.X(), aP2d.Y(), aP2, aDSU, aDSV);
-      
-      aP1.SetXYZ(aP1.XYZ() - aP2.XYZ());
-      aP2.SetXYZ(aDC3D.XYZ() - aDC2D.X()*aDSU.XYZ() - aDC2D.Y()*aDSV.XYZ());
-      
-      theGrad(1) = -2.0*aP1.XYZ().Dot(aP2.XYZ());
-    }
-    catch(Standard_Failure) {
-      return Standard_False;
-    }
-    
-    return Standard_True;
-  }
-  //   
-  virtual Standard_Boolean Values(const math_Vector& theX,
-                                  Standard_Real& theVal,
-                                  math_Vector& theGrad) {
-    if(!Value(theX, theVal))
-      return Standard_False;
-    
-    if(!Gradient(theX, theGrad))
-      return Standard_False;
-    
-    return Standard_True;
-  }
-  //
-  virtual Standard_Boolean Values(const math_Vector& theX,
-                                  Standard_Real& theVal,
-                                  math_Vector& theGrad,
-                                  math_Matrix& theHessian) {
-    if(!Value(theX, theVal))
-      return Standard_False;
-    
-    if(!Gradient(theX, theGrad))
-      return Standard_False;
-    
-    theHessian(1,1) = theGrad(1);
-    
-    return Standard_True;
-  }
-  //
- private:
-
-  Standard_Boolean CheckParameter(const Standard_Real theParam)
-  {
-    if (theParam < my3DCurve->FirstParameter() ||
-        theParam > my3DCurve->LastParameter()  ||
-        theParam < my2DCurve->FirstParameter() ||
-        theParam > my2DCurve->LastParameter() )
-    {
-      return Standard_False;
-    }
-    return Standard_True;
-  }
-
-  Handle(Geom_Curve) my3DCurve;
-  Handle(Geom2d_Curve) my2DCurve;
-  Handle(Geom_Surface) mySurf;
-};
-//=======================================================================
 //
 //=======================================================================
 // Function : CorrectTolerances
@@ -1165,124 +1043,6 @@ void UpdateEdges(const TopoDS_Face& aF)
     }
   }
 }
-//=======================================================================
-// Function : MinComputing
-// purpose : 
-//=======================================================================
-static Standard_Boolean MinComputing( BOPTools_CheckCurveOnSurface& theFunction,
-                                      const Standard_Real theFirst,
-                                      const Standard_Real theLast,
-                                      const Standard_Real theEpsilon, //1.0e-3
-                                      Standard_Real & theBestValue,
-                                      Standard_Real & theBestParameter)
-{
-  //Standard_Real aPrevValue = theBestValue;
-  const Standard_Real aStepMin = 1.0e-2;
-  math_Vector aFirstV(1, 1), aLastV(1, 1), anOutputParam(1, 1);
-  aFirstV(1) = theFirst;
-  aLastV(1) = theLast;
-
-  math_GlobOptMin aFinder(&theFunction, aFirstV, aLastV);
-  aFinder.SetTol(aStepMin, theEpsilon);
-  aFinder.Perform();
-
-  const Standard_Integer aNbExtr = aFinder.NbExtrema();
-  for(Standard_Integer i = 1; i <= aNbExtr; i++)
-  {
-    Standard_Real aValue = 0.0;
-    aFinder.Points(i, anOutputParam);
-    theFunction.Value(anOutputParam, aValue);
-    
-    if(aValue < theBestValue)
-    {
-      theBestValue = aValue;
-      theBestParameter = anOutputParam(1);
-    }
-  }
-
-  return Standard_True;
-}
-
-//=======================================================================
-// Function : ComputeTolerance
-// purpose : 
-//=======================================================================
-Standard_Boolean BOPTools_AlgoTools::ComputeTolerance
-  (const Handle(Geom_Curve)& theCurve3D,
-   const Handle(Geom2d_Curve)& theCurve2D,
-   const Handle(Geom_Surface)& theSurf,
-   Standard_Real& theMaxDist,
-   Standard_Real& theMaxPar)
-{
-  if (theCurve3D.IsNull() ||
-      theCurve2D.IsNull() ||
-      theSurf.IsNull()) {
-    return Standard_False;
-  }
-
-  const Standard_Real anEpsilonRange = 1.0e-3, aMinDelta = 1.0e-5;
-
-  //
-  try {
-    Standard_Real aFirst = theCurve3D->FirstParameter(),
-                  aLast = theCurve3D->LastParameter();
-
-    BOPTools_CheckCurveOnSurface aFunc(theCurve3D, theCurve2D, theSurf);
-    //
-    math_Vector anOutputParam(1, 1);
-    anOutputParam(1) = theMaxPar = aFirst;
-    //
-    theMaxDist = 0.;
-    MinComputing(aFunc, aFirst, aLast, anEpsilonRange, theMaxDist, theMaxPar);
-
-    Standard_Integer aNbIteration = 100;
-    Standard_Boolean aStatus = Standard_True;
-    while((aNbIteration-- >= 0) && aStatus)
-    {
-      Standard_Real aValue = theMaxDist, aParam = theMaxPar;
-      Standard_Real aBP = theMaxPar - aMinDelta;
-      MinComputing(aFunc, aFirst, aBP, anEpsilonRange, theMaxDist, theMaxPar);
-
-      if(theMaxDist < aValue)
-      {
-        aLast = aBP;
-        aStatus = Standard_True;
-      }
-      else 
-      {
-        theMaxDist = aValue;
-        theMaxPar = aParam;
-        aStatus = Standard_False;
-      }
-
-      if(!aStatus)
-      {
-        aBP = theMaxPar + aMinDelta;
-        MinComputing(aFunc, aBP, aLast, 1.0e-3, theMaxDist, theMaxPar);
-
-        if(theMaxDist < aValue)
-        {
-          aFirst = aBP;
-          aStatus = Standard_True;
-        }
-        else
-        {
-          theMaxDist = aValue;
-          theMaxPar = aParam;
-          aStatus = Standard_False;
-        }
-      }
-    }
-
-    theMaxDist = sqrt(Abs(theMaxDist));
-  }
-  catch (Standard_Failure) {
-    return Standard_False;
-  }
-  //
-  return Standard_True;
-}
-
 //=======================================================================
 // Function : ComputeTolerance
 // purpose : 
@@ -1291,111 +1051,18 @@ Standard_Boolean BOPTools_AlgoTools::ComputeTolerance
   (const TopoDS_Face& theFace,
    const TopoDS_Edge& theEdge,
    Standard_Real& theMaxDist,
-   Standard_Real& theParameter)
+   Standard_Real& theMaxPar)
 {
-  Standard_Boolean bRet;
-  Standard_Real aT, aD, aFirst, aLast;
-  TopLoc_Location aLocC, aLocS;
+  BRepLib_CheckCurveOnSurface aCS;
   //
-  theMaxDist = 0.;
-  theParameter = 0.;
-  bRet = Standard_False;
-  //
-  const Handle(BRep_TEdge)& aTE = *((Handle(BRep_TEdge)*)&theEdge.TShape());
-  //The edge is considered to be same range and not degenerated
-  if ((!aTE->SameRange() && aTE->SameParameter()) ||
-      aTE->Degenerated()) {
-    return bRet;
+  aCS.Init(theEdge, theFace);
+  aCS.Perform();
+  if (!aCS.IsDone()) {
+    return Standard_False;
   }
   //
-  Handle(Geom_Curve) aC = Handle(Geom_Curve)::
-    DownCast(BRep_Tool::Curve(theEdge, aLocC, aFirst, aLast)->Copy());
-  aC = new Geom_TrimmedCurve(aC, aFirst, aLast);
-  aC->Transform(aLocC.Transformation());
+  theMaxDist = aCS.MaxDistance();
+  theMaxPar  = aCS.MaxParameter();
   //
-  const Handle(Geom_Surface)& aSurfF = BRep_Tool::Surface(theFace, aLocS);
-  const Handle(Geom_Surface)& aSurf = Handle(Geom_Surface)::
-    DownCast(aSurfF->Copy()->Transformed(aLocS.Transformation()));
-  //
-  Standard_Boolean isPCurveFound = Standard_False;
-  BRep_ListIteratorOfListOfCurveRepresentation itcr(aTE->Curves());
-  for (; itcr.More(); itcr.Next()) {
-    const Handle(BRep_CurveRepresentation)& cr = itcr.Value();
-    if (!(cr->IsCurveOnSurface(aSurfF, aLocS.Predivided(theEdge.Location())))) {
-      continue;
-    }
-    isPCurveFound = Standard_True;
-    //
-    Handle(Geom2d_Curve) aC2d = Handle(Geom2d_Curve)::
-      DownCast(cr->PCurve()->Copy());
-    aC2d = new Geom2d_TrimmedCurve(aC2d, aFirst, aLast);
-      //
-    if(BOPTools_AlgoTools::ComputeTolerance
-       (aC, aC2d, aSurf, aD, aT)) {
-      bRet = Standard_True;
-      if (aD > theMaxDist) {
-        theMaxDist = aD;
-        theParameter = aT;
-      }
-    }
-    //
-    if (cr->IsCurveOnClosedSurface()) {
-      Handle(Geom2d_Curve) aC2d = Handle(Geom2d_Curve)::
-        DownCast(cr->PCurve2()->Copy());
-      aC2d = new Geom2d_TrimmedCurve(aC2d, aFirst, aLast);
-      //
-      if(BOPTools_AlgoTools::ComputeTolerance
-         (aC, aC2d, aSurf, aD, aT)) {
-        bRet = Standard_True;
-        if (aD > theMaxDist) {
-          theMaxDist = aD;
-          theParameter = aT;
-        }
-      }
-    }
-  }
-  //
-  if (isPCurveFound) {
-    return bRet;
-  }
-  //
-  Handle(Geom_Plane) aPlane;
-  Handle(Standard_Type) dtyp = aSurf->DynamicType();
-  //
-  if (dtyp == STANDARD_TYPE(Geom_RectangularTrimmedSurface)) {
-    aPlane = Handle(Geom_Plane)::
-      DownCast(Handle(Geom_RectangularTrimmedSurface)::
-               DownCast(aSurf)->BasisSurface()->Copy());
-  }
-  else {
-    aPlane = Handle(Geom_Plane)::DownCast(aSurf->Copy());
-  }
-  //
-  if (aPlane.IsNull()) { // not a plane
-    return bRet;
-  }
-  //
-  aPlane = Handle(Geom_Plane)::DownCast(aPlane);//
-  //
-  Handle(GeomAdaptor_HSurface) GAHS = new GeomAdaptor_HSurface(aPlane);
-  Handle(Geom_Curve) ProjOnPlane = 
-    GeomProjLib::ProjectOnPlane (new Geom_TrimmedCurve(aC, aFirst, aLast), 
-                                 aPlane, aPlane->Position().Direction(), 
-                                 Standard_True);
-  Handle(GeomAdaptor_HCurve) aHCurve = new GeomAdaptor_HCurve(ProjOnPlane);
-  //
-  ProjLib_ProjectedCurve proj(GAHS,aHCurve);
-  Handle(Geom2d_Curve) aC2d = Geom2dAdaptor::MakeCurve(proj);
-  aC2d = new Geom2d_TrimmedCurve(aC2d, aFirst, aLast);
-  //
-  if(BOPTools_AlgoTools::ComputeTolerance
-     (aC, aC2d, aPlane, aD, aT)) {
-    bRet = Standard_True;
-    if (aD > theMaxDist) {
-      theMaxDist = aD;
-      theParameter = aT;
-    }
-  }
-  //
-  return bRet;
+  return Standard_True;
 }
diff --git a/src/BRepLib/BRepLib.cdl b/src/BRepLib/BRepLib.cdl
index d950f3eba7..cc7935ea9e 100644
--- a/src/BRepLib/BRepLib.cdl
+++ b/src/BRepLib/BRepLib.cdl
@@ -133,6 +133,10 @@ is
     class FuseEdges;
 	---Purpose: 
 
+    class CheckCurveOnSurface;
+    ---Purpose:
+    -- Computes the max distance between edge
+    -- and its 2d representation on the face.
 
     --
     --   Default precison methods.
diff --git a/src/BRepLib/BRepLib_CheckCurveOnSurface.cdl b/src/BRepLib/BRepLib_CheckCurveOnSurface.cdl
new file mode 100644
index 0000000000..cf893a41a5
--- /dev/null
+++ b/src/BRepLib/BRepLib_CheckCurveOnSurface.cdl
@@ -0,0 +1,168 @@
+-- Created by: Eugeny MALTCHIKOV
+-- Copyright (c) 2014 OPEN CASCADE SAS
+--
+-- This file is part of Open CASCADE Technology software library.
+--
+-- This library is free software; you can redistribute it and/or modify it under
+-- the terms of the GNU Lesser General Public License version 2.1 as published
+-- by the Free Software Foundation, with special exception defined in the file
+-- OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
+-- distribution for complete text of the license and disclaimer of any warranty.
+--
+-- Alternatively, this file may be used under the terms of Open CASCADE
+-- commercial license or contractual agreement.
+
+class CheckCurveOnSurface from BRepLib
+
+  ---Purpose:
+  -- Computes the max distance between edge and its
+  -- 2d representation on the face.
+  --
+  -- The algorithm can be initialized in the following ways:
+  -- 1. Input args are Edge and Face;
+  -- 2. Input args are 3D curve, 2d curve, Surface and
+  --    parametric range of the curve (first and last values).
+
+uses
+
+    Edge    from TopoDS,
+    Face    from TopoDS,
+    Curve   from Geom,
+    Curve   from Geom2d,
+    Surface from Geom
+
+is
+
+    Create
+    returns CheckCurveOnSurface from BRepLib;
+    ---Purpose:
+    -- Empty contructor
+
+    Create(
+        theEdge : Edge from TopoDS;
+        theFace : Face from TopoDS);
+    ---Purpose:
+    -- Contructor 
+
+    Create(
+        theCurve   : Curve   from Geom;
+        thePCurve  : Curve   from Geom2d;
+        theSurface : Surface from Geom;
+        theFirst   : Real from Standard;
+        theLast    : Real from Standard);
+    ---Purpose:
+    -- Contructor
+
+    Init(me:out;
+        theEdge : Edge from TopoDS;
+        theFace : Face from TopoDS);
+    ---Purpose:
+    -- Sets the data for the algorithm
+
+    Init(me:out;
+        theCurve   : Curve   from Geom;
+        thePCurve  : Curve   from Geom2d;
+        theSurface : Surface from Geom;
+        theFirst   : Real from Standard;
+        theLast    : Real from Standard);
+    ---Purpose:
+    -- Sets the data for the algorithm
+
+    Curve(me)
+    returns Curve from Geom;
+    ---C++: inline
+    ---C++: return const &
+    ---Purpose:
+    -- Returns my3DCurve
+
+    PCurve(me)
+    returns Curve from Geom2d;
+    ---C++: inline
+    ---C++: return const &
+    ---Purpose:
+    -- Returns my2DCurve
+
+    PCurve2(me)
+    returns Curve from Geom2d;
+    ---C++: inline
+    ---C++: return const &
+    ---Purpose:
+    -- Returns my2DCurve
+
+    Surface(me)
+    returns Surface from Geom;
+    ---C++: inline
+    ---C++: return const &
+    ---Purpose:
+    -- Returns mySurface
+
+    Range(me:out;
+        theFirst : out Real from Standard;
+        theLast  : out Real from Standard);
+    ---C++: inline
+    ---Purpose:
+    -- Returns the range
+
+    -- computations
+    --
+    Perform(me:out);
+    ---Purpose:
+    -- Performs the calculation
+
+    CheckData(me:out)
+    is protected;
+    ---Purpose:
+    -- Checks the data
+
+    Compute(me:out;
+        thePCurve : Curve from Geom2d)
+    is protected;
+    ---Purpose:
+    -- Computes the max distance for the 3d curve <myCurve>
+    -- and 2d curve <thePCurve>
+
+    -- results
+    --
+    IsDone(me)
+    returns Boolean from Standard;
+    ---C++: inline
+    ---Purpose:
+    -- Returns true if the max distance has been found
+
+    ErrorStatus(me)
+    returns Integer from Standard;
+    ---C++: inline
+    ---Purpose:
+    -- Returns error status
+    -- The possible values are:
+    -- 0 - OK;
+    -- 1 - null curve or surface or 2d curve;
+    -- 2 - invalid parametric range;
+    -- 3 - error in calculations.
+
+    MaxDistance(me)
+    returns Real from Standard;
+    ---C++: inline
+    ---Purpose:
+    -- Returns max distance
+
+    MaxParameter(me)
+    returns Real from Standard;
+    ---C++: inline
+    ---Purpose:
+    -- Returns parameter in which the distance is maximal
+
+fields
+    -- source data
+    myCurve   : Curve   from Geom;
+    myPCurve  : Curve   from Geom2d;
+    myPCurve2 : Curve   from Geom2d;
+    mySurface : Surface from Geom;
+    myFirst   : Real    from Standard;
+    myLast    : Real    from Standard;
+    -- result
+    myErrorStatus  : Integer from Standard;
+    myMaxDistance  : Real    from Standard;
+    myMaxParameter : Real    from Standard;
+
+end CheckCurveOnSurface;
diff --git a/src/BRepLib/BRepLib_CheckCurveOnSurface.cxx b/src/BRepLib/BRepLib_CheckCurveOnSurface.cxx
new file mode 100644
index 0000000000..815ad6fcb8
--- /dev/null
+++ b/src/BRepLib/BRepLib_CheckCurveOnSurface.cxx
@@ -0,0 +1,456 @@
+// Created by: Eugeny MALTCHIKOV
+// Copyright (c) 2014 OPEN CASCADE SAS
+//
+// This file is part of Open CASCADE Technology software library.
+//
+// This library is free software; you can redistribute it and/or modify it under
+// the terms of the GNU Lesser General Public License version 2.1 as published
+// by the Free Software Foundation, with special exception defined in the file
+// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
+// distribution for complete text of the license and disclaimer of any warranty.
+//
+// Alternatively, this file may be used under the terms of Open CASCADE
+// commercial license or contractual agreement.
+
+#include <BRepLib_CheckCurveOnSurface.ixx>
+
+#include <math_GlobOptMin.hxx>
+#include <math_MultipleVarFunctionWithHessian.hxx>
+#include <math_Matrix.hxx>
+
+#include <Geom_Plane.hxx>
+#include <Geom_RectangularTrimmedSurface.hxx>
+#include <Geom_TrimmedCurve.hxx>
+
+#include <Geom2dAdaptor.hxx>
+
+#include <GeomAdaptor_HSurface.hxx>
+#include <GeomAdaptor_HCurve.hxx>
+
+#include <GeomProjLib.hxx>
+
+#include <BRep_Tool.hxx>
+#include <BRep_TEdge.hxx>
+#include <BRep_CurveRepresentation.hxx>
+#include <BRep_ListIteratorOfListOfCurveRepresentation.hxx>
+
+#include <TopLoc_Location.hxx>
+
+#include <ProjLib_ProjectedCurve.hxx>
+
+
+//=======================================================================
+//class   : BRepLib_CheckCurveOnSurface_GlobOptFunc
+//purpose : provides necessary methods to be used in math_GlobOptMin
+//=======================================================================
+class BRepLib_CheckCurveOnSurface_GlobOptFunc :
+  public math_MultipleVarFunctionWithHessian
+{
+ public:
+  BRepLib_CheckCurveOnSurface_GlobOptFunc
+    (BRepLib_CheckCurveOnSurface_GlobOptFunc&);
+  BRepLib_CheckCurveOnSurface_GlobOptFunc
+    (const Handle(Geom_Curve)& theC3D,
+     const Handle(Geom2d_Curve)& theC2D,
+     const Handle(Geom_Surface)& theSurf,
+     const Standard_Real theFirst,
+     const Standard_Real theLast)
+    :
+      myCurve(theC3D),
+      myPCurve(theC2D),
+      mySurf(theSurf),
+      myFirst(theFirst),
+      myLast(theLast)
+  {
+  }
+  //
+  virtual Standard_Integer NbVariables() const {
+    return 1;
+  }
+  //
+  virtual Standard_Boolean Value(const math_Vector& theX,
+                                 Standard_Real& theFVal) {
+    try {
+      const Standard_Real aPar = theX(1);
+      if (!CheckParameter(aPar))
+        return Standard_False;
+      gp_Pnt aP1, aP2;
+      gp_Pnt2d aP2d;
+      //
+      myCurve->D0(aPar, aP1);
+      myPCurve->D0(aPar, aP2d);
+      mySurf->D0(aP2d.X(), aP2d.Y(), aP2);
+      //
+      theFVal = -1.0*aP1.SquareDistance(aP2);
+    }
+    catch(Standard_Failure) {
+      return Standard_False;
+    }
+    //
+    return Standard_True;
+  }
+  //
+  virtual Standard_Integer GetStateNumber() {
+    return 0;
+  }
+  //
+  virtual Standard_Boolean Gradient(const math_Vector& theX,
+                                    math_Vector& theGrad) {
+    try {
+      const Standard_Real aPar = theX(1);
+      if (!CheckParameter(aPar)) {
+        return Standard_False;
+      }
+      //
+      gp_Pnt aP1, aP2;
+      gp_Vec aDC3D, aDSU, aDSV;
+      gp_Pnt2d aP2d;
+      gp_Vec2d aDC2D;
+      //
+      myCurve->D1(aPar, aP1, aDC3D);
+      myPCurve->D1(aPar, aP2d, aDC2D);
+      mySurf->D1(aP2d.X(), aP2d.Y(), aP2, aDSU, aDSV);
+      //
+      aP1.SetXYZ(aP1.XYZ() - aP2.XYZ());
+      aP2.SetXYZ(aDC3D.XYZ() - aDC2D.X()*aDSU.XYZ() - aDC2D.Y()*aDSV.XYZ());
+      //
+      theGrad(1) = -2.0*aP1.XYZ().Dot(aP2.XYZ());
+    }
+    catch(Standard_Failure) {
+      return Standard_False;
+    }
+    //
+    return Standard_True;
+  }
+  //   
+  virtual Standard_Boolean Values(const math_Vector& theX,
+                                  Standard_Real& theVal,
+                                  math_Vector& theGrad) {
+    if (!Value(theX, theVal)) {
+      return Standard_False;
+    }
+    //
+    if (!Gradient(theX, theGrad)) {
+      return Standard_False;
+    }
+    //
+    return Standard_True;
+  }
+  //
+  virtual Standard_Boolean Values(const math_Vector& theX,
+                                  Standard_Real& theVal,
+                                  math_Vector& theGrad,
+                                  math_Matrix& theHessian) {
+    if (!Value(theX, theVal)) {
+      return Standard_False;
+    }
+    //    
+    if (!Gradient(theX, theGrad)) {
+      return Standard_False;
+    }
+    //
+    theHessian(1,1) = theGrad(1);
+    //
+    return Standard_True;
+  }
+  //
+ private:
+
+  Standard_Boolean CheckParameter(const Standard_Real theParam) {
+    return ((myFirst <= theParam) && (theParam <= myLast));
+  }
+
+  Handle(Geom_Curve) myCurve;
+  Handle(Geom2d_Curve) myPCurve;
+  Handle(Geom_Surface) mySurf;
+  Standard_Real myFirst;
+  Standard_Real myLast;
+};
+
+static 
+  void MinComputing(BRepLib_CheckCurveOnSurface_GlobOptFunc& theFunction,
+                    const Standard_Real theFirst,
+                    const Standard_Real theLast,
+                    const Standard_Real theEpsilon,
+                    Standard_Real& theBestValue,
+                    Standard_Real& theBestParameter);
+
+
+//=======================================================================
+//function : BRepLib_CheckCurveOnSurface
+//purpose  : 
+//=======================================================================
+BRepLib_CheckCurveOnSurface::BRepLib_CheckCurveOnSurface()
+:
+  myFirst(0.),
+  myLast(0.),
+  myErrorStatus(0),
+  myMaxDistance(0.),
+  myMaxParameter(0.)
+{
+}
+
+//=======================================================================
+//function : BRepLib_CheckCurveOnSurface
+//purpose  : 
+//=======================================================================
+BRepLib_CheckCurveOnSurface::BRepLib_CheckCurveOnSurface
+  (const TopoDS_Edge& theEdge,
+   const TopoDS_Face& theFace)
+:
+  myErrorStatus(0),
+  myMaxDistance(0.),
+  myMaxParameter(0.)
+{
+  Init(theEdge, theFace);
+}
+
+//=======================================================================
+//function : BRepLib_CheckCurveOnSurface
+//purpose  : 
+//=======================================================================
+BRepLib_CheckCurveOnSurface::BRepLib_CheckCurveOnSurface
+  (const Handle(Geom_Curve)& the3DCurve,
+   const Handle(Geom2d_Curve)& the2DCurve,
+   const Handle(Geom_Surface)& theSurface,
+   const Standard_Real theFirst,
+   const Standard_Real theLast)
+:
+  myErrorStatus(0),
+  myMaxDistance(0.),
+  myMaxParameter(0.)
+{
+  Init(the3DCurve, the2DCurve, theSurface, theFirst, theLast);
+}
+
+//=======================================================================
+//function : Init
+//purpose  : 
+//=======================================================================
+void BRepLib_CheckCurveOnSurface::Init
+  (const TopoDS_Edge& theEdge,
+   const TopoDS_Face& theFace)
+{
+  if (theEdge.IsNull() || theFace.IsNull()) {
+    return;
+  }
+  //
+  if (BRep_Tool::Degenerated(theEdge) ||
+      !BRep_Tool::IsGeometric(theEdge)) {
+    return;
+  }
+  //
+  Standard_Boolean isPCurveFound;
+  TopLoc_Location aLocE, aLocF, aLocC2D;
+  //
+  // 3D curve initialization
+  myCurve = Handle(Geom_Curve)::
+    DownCast(BRep_Tool::Curve(theEdge, aLocE, myFirst, myLast)->Copy());
+  myCurve->Transform(aLocE.Transformation());
+  //
+  // Surface initialization
+  const Handle(Geom_Surface)& aS = BRep_Tool::Surface(theFace, aLocF);
+  mySurface = Handle(Geom_Surface)::
+    DownCast(aS->Copy()->Transformed(aLocF.Transformation()));
+  //
+  // 2D curves initialization 
+  isPCurveFound = Standard_False;
+  aLocC2D = aLocF.Predivided(aLocE);
+  const Handle(BRep_TEdge)& aTE = *((Handle(BRep_TEdge)*)&theEdge.TShape());
+  BRep_ListIteratorOfListOfCurveRepresentation itcr(aTE->Curves());
+  //
+  for (; itcr.More(); itcr.Next()) {
+    const Handle(BRep_CurveRepresentation)& cr = itcr.Value();
+    if (cr->IsCurveOnSurface(aS, aLocC2D)) {
+      isPCurveFound = Standard_True;
+      myPCurve = cr->PCurve();
+      //
+      if (cr->IsCurveOnClosedSurface()) {
+        myPCurve2 = cr->PCurve2();
+      }
+      break;
+    }
+  }
+  //
+  if (isPCurveFound) {
+    return;
+  }
+  //
+  Handle(Geom_Plane) aPlane;
+  Handle(Standard_Type) dtyp = mySurface->DynamicType();
+  //
+  if (dtyp == STANDARD_TYPE(Geom_RectangularTrimmedSurface)) {
+    aPlane = Handle(Geom_Plane)::
+      DownCast(Handle(Geom_RectangularTrimmedSurface)::
+               DownCast(mySurface)->BasisSurface()->Copy());
+  }
+  else {
+    aPlane = Handle(Geom_Plane)::DownCast(mySurface->Copy());
+  }
+  //
+  if (aPlane.IsNull()) { // not a plane
+    return;
+  }
+  //
+  aPlane = Handle(Geom_Plane)::DownCast(aPlane);
+  //
+  Handle(GeomAdaptor_HSurface) aGAHS = new GeomAdaptor_HSurface(aPlane);
+  Handle(Geom_Curve) aProjOnPlane = 
+    GeomProjLib::ProjectOnPlane (new Geom_TrimmedCurve(myCurve, myFirst, myLast), 
+                                 aPlane, aPlane->Position().Direction(), 
+                                 Standard_True);
+  Handle(GeomAdaptor_HCurve) aHCurve = new GeomAdaptor_HCurve(aProjOnPlane);
+  //
+  ProjLib_ProjectedCurve aProj(aGAHS, aHCurve);
+  myPCurve = Geom2dAdaptor::MakeCurve(aProj);
+}
+
+//=======================================================================
+//function : Init
+//purpose  : 
+//=======================================================================
+void BRepLib_CheckCurveOnSurface::Init
+  (const Handle(Geom_Curve)& the3DCurve,
+   const Handle(Geom2d_Curve)& the2DCurve,
+   const Handle(Geom_Surface)& theSurface,
+   const Standard_Real theFirst,
+   const Standard_Real theLast)
+{
+  myCurve = the3DCurve;
+  myPCurve = the2DCurve;
+  mySurface = theSurface;
+  myFirst = theFirst;
+  myLast  = theLast;
+}
+
+//=======================================================================
+//function : Perform
+//purpose  : 
+//=======================================================================
+void BRepLib_CheckCurveOnSurface::Perform()
+{
+  try {
+    //
+    // 1. Check data
+    CheckData();
+    if (myErrorStatus) {
+      return;
+    }
+    //
+    // 2. Compute the max distance
+    Compute(myPCurve);
+    //
+    if (!myPCurve2.IsNull()) {
+      // compute max distance for myPCurve2
+      // (for the second curve on closed surface)
+      Compute(myPCurve2);
+    }
+  }
+  catch (Standard_Failure) {
+    myErrorStatus = 3;
+  }
+}
+
+//=======================================================================
+//function : Compute
+//purpose  : 
+//=======================================================================
+void BRepLib_CheckCurveOnSurface::Compute
+  (const Handle(Geom2d_Curve)& thePCurve)
+{
+  Standard_Integer aNbIt, aStatus;
+  Standard_Real anEpsilonRange, aMinDelta;
+  Standard_Real aFirst, aLast;
+  Standard_Real aValue, aParam, aBP;
+  Standard_Real theMaxDist, theMaxPar;
+  //
+  anEpsilonRange = 1.e-3;
+  aMinDelta = 1.e-5;
+  aFirst = myFirst;
+  aLast  = myLast;
+  //
+  BRepLib_CheckCurveOnSurface_GlobOptFunc aFunc
+    (myCurve, thePCurve, mySurface, myFirst, myLast);
+  //
+  math_Vector anOutputParam(1, 1);
+  anOutputParam(1) = aFirst;
+  theMaxDist = 0.;
+  theMaxPar = aFirst;
+  aNbIt = 100;
+  aStatus = Standard_True;
+  //
+  MinComputing(aFunc, aFirst, aLast, anEpsilonRange, theMaxDist, theMaxPar);
+  //
+  while((aNbIt-- >= 0) && aStatus) {
+    aValue = theMaxDist;
+    aParam = theMaxPar;
+    aBP = theMaxPar - aMinDelta;
+    MinComputing(aFunc, aFirst, aBP, anEpsilonRange, theMaxDist, theMaxPar);
+    //
+    if(theMaxDist < aValue) {
+      aLast = aBP;
+      aStatus = Standard_True;
+    }
+    else {
+      theMaxDist = aValue;
+      theMaxPar = aParam;
+      aStatus = Standard_False;
+    }
+    //
+    if(!aStatus) {
+      aBP = theMaxPar + aMinDelta;
+      MinComputing(aFunc, aBP, aLast, 1.0e-3, theMaxDist, theMaxPar);
+      //
+      if(theMaxDist < aValue) {
+        aFirst = aBP;
+        aStatus = Standard_True;
+      }
+      else {
+        theMaxDist = aValue;
+        theMaxPar = aParam;
+        aStatus = Standard_False;
+      }
+    }
+  }
+  //
+  theMaxDist = sqrt(Abs(theMaxDist));
+  if (theMaxDist > myMaxDistance) {
+    myMaxDistance = theMaxDist;
+    myMaxParameter = theMaxPar;
+  }  
+}
+
+//=======================================================================
+// Function : MinComputing
+// purpose : 
+//=======================================================================
+void MinComputing
+  (BRepLib_CheckCurveOnSurface_GlobOptFunc& theFunction,
+   const Standard_Real theFirst,
+   const Standard_Real theLast,
+   const Standard_Real theEpsilon, //1.0e-3
+   Standard_Real& theBestValue,
+   Standard_Real& theBestParameter)
+{
+  const Standard_Real aStepMin = 1.0e-2;
+  math_Vector aFirstV(1, 1), aLastV(1, 1), anOutputParam(1, 1);
+  aFirstV(1) = theFirst;
+  aLastV(1) = theLast;
+  //
+  math_GlobOptMin aFinder(&theFunction, aFirstV, aLastV);
+  aFinder.SetTol(aStepMin, theEpsilon);
+  aFinder.Perform();
+  //
+  const Standard_Integer aNbExtr = aFinder.NbExtrema();
+  for(Standard_Integer i = 1; i <= aNbExtr; i++)
+  {
+    Standard_Real aValue = 0.0;
+    aFinder.Points(i, anOutputParam);
+    theFunction.Value(anOutputParam, aValue);
+    //
+    if(aValue < theBestValue) {
+      theBestValue = aValue;
+      theBestParameter = anOutputParam(1);
+    }
+  }
+}
diff --git a/src/BRepLib/BRepLib_CheckCurveOnSurface.lxx b/src/BRepLib/BRepLib_CheckCurveOnSurface.lxx
new file mode 100644
index 0000000000..c59a29e2a1
--- /dev/null
+++ b/src/BRepLib/BRepLib_CheckCurveOnSurface.lxx
@@ -0,0 +1,109 @@
+// Created by: Eugeny MALTCHIKOV
+// Copyright (c) 2014 OPEN CASCADE SAS
+//
+// This file is part of Open CASCADE Technology software library.
+//
+// This library is free software; you can redistribute it and/or modify it under
+// the terms of the GNU Lesser General Public License version 2.1 as published
+// by the Free Software Foundation, with special exception defined in the file
+// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
+// distribution for complete text of the license and disclaimer of any warranty.
+//
+// Alternatively, this file may be used under the terms of Open CASCADE
+// commercial license or contractual agreement.
+
+//=======================================================================
+//function : Curve
+//purpose  : 
+//=======================================================================
+inline const Handle(Geom_Curve)& BRepLib_CheckCurveOnSurface::Curve() const
+{
+  return myCurve;
+}
+
+//=======================================================================
+//function : PCurve
+//purpose  : 
+//=======================================================================
+inline const Handle(Geom2d_Curve)& BRepLib_CheckCurveOnSurface::PCurve() const
+{
+  return myPCurve;
+}
+
+//=======================================================================
+//function : PCurve2
+//purpose  : 
+//=======================================================================
+inline const Handle(Geom2d_Curve)& BRepLib_CheckCurveOnSurface::PCurve2() const
+{
+  return myPCurve2;
+}
+
+//=======================================================================
+//function : Surface
+//purpose  : 
+//=======================================================================
+inline const Handle(Geom_Surface)& BRepLib_CheckCurveOnSurface::Surface() const
+{
+  return mySurface;
+}
+
+//=======================================================================
+//function : Range
+//purpose  : 
+//=======================================================================
+inline void BRepLib_CheckCurveOnSurface::Range
+  (Standard_Real& theFirst,
+   Standard_Real& theLast)
+{
+  theFirst = myFirst;
+  theLast  = myLast;
+}
+
+//=======================================================================
+//function : CheckData
+//purpose  : 
+//=======================================================================
+inline void BRepLib_CheckCurveOnSurface::CheckData()
+{
+  if (myCurve.IsNull() ||
+      myPCurve.IsNull() ||
+      mySurface.IsNull()) {
+    myErrorStatus = 1;
+    return;
+  }
+  //
+  if ((myCurve->FirstParameter() > myFirst) ||
+      (myCurve->LastParameter() < myLast) ||
+      (myPCurve->FirstParameter() > myFirst) ||
+      (myPCurve->LastParameter() < myLast)) {
+    myErrorStatus = 2;
+  }
+}
+
+//=======================================================================
+//function : IsDone
+//purpose  : 
+//=======================================================================
+inline Standard_Boolean BRepLib_CheckCurveOnSurface::IsDone() const
+{
+  return (myErrorStatus == 0);
+}
+
+//=======================================================================
+//function : MaxDistance
+//purpose  : 
+//=======================================================================
+inline Standard_Real BRepLib_CheckCurveOnSurface::MaxDistance() const
+{
+  return myMaxDistance;
+}
+
+//=======================================================================
+//function : MaxParameter
+//purpose  : 
+//=======================================================================
+inline Standard_Real BRepLib_CheckCurveOnSurface::MaxParameter() const
+{
+  return myMaxParameter;
+}
diff --git a/src/IntTools/IntTools_FaceFace.cdl b/src/IntTools/IntTools_FaceFace.cdl
index ef9573c3a3..2fa6a074a4 100644
--- a/src/IntTools/IntTools_FaceFace.cdl
+++ b/src/IntTools/IntTools_FaceFace.cdl
@@ -133,6 +133,10 @@ is
     ComputeTolReached3d(me:out) 
      is  protected; 
  
+    ComputeTolerance(me:out)  
+    returns Real from Standard
+    is  protected; 
+ 
     SetContext(me:out; 
        aContext : Context from IntTools); 
     ---Purpose:
diff --git a/src/IntTools/IntTools_FaceFace.cxx b/src/IntTools/IntTools_FaceFace.cxx
index 2f658e4e0a..a06269587c 100644
--- a/src/IntTools/IntTools_FaceFace.cxx
+++ b/src/IntTools/IntTools_FaceFace.cxx
@@ -281,34 +281,27 @@ static
   Standard_Integer IndexType(const GeomAbs_SurfaceType aType);
 
 //
-static
-  Standard_Real MaxSquareDistance (const Standard_Real aT,
-                                   const Handle(Geom_Curve)& aC3D,
-                                   const Handle(Geom2d_Curve)& aC2D1,
-                                   const Handle(Geom2d_Curve)& aC2D2,
-                                   const Handle(GeomAdaptor_HSurface) myHS1,
-                                   const Handle(GeomAdaptor_HSurface) myHS2,
-                                   const TopoDS_Face& aF1,
-                                   const TopoDS_Face& aF2,
-                                   const Handle(IntTools_Context)& aCtx);
-
 static
   Standard_Boolean CheckPCurve(const Handle(Geom2d_Curve)& aPC, 
                                const TopoDS_Face& aFace);
 
 //
 static
-  Standard_Real FindMaxSquareDistance (const Standard_Real aA,
-                                       const Standard_Real aB,
-                                       const Standard_Real aEps,
-                                       const Handle(Geom_Curve)& aC3D,
-                                       const Handle(Geom2d_Curve)& aC2D1,
-                                       const Handle(Geom2d_Curve)& aC2D2,
-                                       const Handle(GeomAdaptor_HSurface)& myHS1,
-                                       const Handle(GeomAdaptor_HSurface)& myHS2,
-                                       const TopoDS_Face& aF1,
-                                       const TopoDS_Face& aF2,
-                                       const Handle(IntTools_Context)& aCtx);
+  Standard_Real MaxDistance(const Handle(Geom_Curve)& theC,
+                            const Standard_Real aT,
+                            GeomAPI_ProjectPointOnSurf& theProjPS);
+static
+  Standard_Real FindMaxDistance(const Handle(Geom_Curve)& theC,
+                                const Standard_Real theFirst,
+                                const Standard_Real theLast,
+                                GeomAPI_ProjectPointOnSurf& theProjPS,
+                                const Standard_Real theEps);
+static
+  Standard_Real FindMaxDistance(const Handle(Geom_Curve)& theCurve,
+                                const Standard_Real theFirst,
+                                const Standard_Real theLast,
+                                const TopoDS_Face& theFace,
+                                const Handle(IntTools_Context)& theContext);
 
 //=======================================================================
 //function : 
@@ -822,13 +815,66 @@ static Standard_Boolean isTreatAnalityc(const TopoDS_Face& theF1,
   }
 }
 
+//=======================================================================
+//function : ComputeTolerance
+//purpose  : 
+//=======================================================================
+Standard_Real IntTools_FaceFace::ComputeTolerance()
+{
+  Standard_Integer i, j, aNbLin;
+  Standard_Real aFirst, aLast, aD, aDMax, aT, aDelta;
+  Handle(Geom_Surface) aS1, aS2;
+  //
+  aDMax = 0;
+  aDelta = Precision::PConfusion();
+  aNbLin = mySeqOfCurve.Length();
+  //
+  aS1 = myHS1->ChangeSurface().Surface();
+  aS2 = myHS2->ChangeSurface().Surface();
+  //
+  for (i = 1; i <= aNbLin; ++i) {
+    const IntTools_Curve& aIC = mySeqOfCurve(i);
+    const Handle(Geom_Curve)& aC3D = aIC.Curve();
+    if (aC3D.IsNull()) {
+      continue;
+    }
+    //
+    aFirst = aC3D->FirstParameter();
+    aLast  = aC3D->LastParameter();
+    //
+    const Handle(Geom2d_Curve)& aC2D1 = aIC.FirstCurve2d();
+    const Handle(Geom2d_Curve)& aC2D2 = aIC.SecondCurve2d();
+    //
+    for (j = 0; j < 2; ++j) {
+      const Handle(Geom2d_Curve)& aC2D = !j ? aC2D1 : aC2D2;
+      const Handle(Geom_Surface)& aS = !j ? aS1 : aS2;
+      //
+      if (!aC2D.IsNull()) {
+        if (IntTools_Tools::ComputeTolerance
+            (aC3D, aC2D, aS, aFirst, aLast, aD, aT)) {
+          if (aD > aDMax) {
+            aDMax = aD;
+          }
+        }
+      }
+      //
+      const TopoDS_Face& aF = !i ? myFace1 : myFace2;
+      aD = FindMaxDistance(aC3D, aFirst, aLast, aF, myContext);
+      if (aD > aDMax) {
+        aDMax = aD;
+      }
+    }
+  }
+  //
+  return aDMax;
+}
+
 //=======================================================================
 //function :ComputeTolReached3d 
 //purpose  : 
 //=======================================================================
   void IntTools_FaceFace::ComputeTolReached3d()
 {
-  Standard_Boolean bCase1;
   Standard_Integer aNbLin, i;
   GeomAbs_SurfaceType aType1, aType2;
   //
@@ -840,9 +886,6 @@ static Standard_Boolean isTreatAnalityc(const TopoDS_Face& theF1,
   aType1=myHS1->Surface().GetType();
   aType2=myHS2->Surface().GetType();
   //
-  bCase1=((aType1==GeomAbs_Plane && aType2==GeomAbs_SurfaceOfExtrusion) ||
-          (aType2==GeomAbs_Plane && aType1==GeomAbs_SurfaceOfExtrusion));
-  //
   if (aType1==GeomAbs_Cylinder && aType2==GeomAbs_Cylinder) {
     if (aNbLin==2){ 
       Handle(IntPatch_Line) aIL1, aIL2;
@@ -871,45 +914,12 @@ static Standard_Boolean isTreatAnalityc(const TopoDS_Face& theF1,
     }
     //ZZ
     if (aNbLin) {// Check the distances
-      Standard_Integer  aNbP, j ;
-      Standard_Real aT1, aT2, dT, aD2, aD2Max, aEps, aT11, aT12;
+      Standard_Real aDMax;
       //
-      aD2Max=0.;
-      aNbP=10;
-      aNbLin=mySeqOfCurve.Length();
-      //
-      for (i=1; i<=aNbLin; ++i) {
-        const IntTools_Curve& aIC=mySeqOfCurve(i);
-        const Handle(Geom_Curve)& aC3D=aIC.Curve();
-        const Handle(Geom2d_Curve)& aC2D1=aIC.FirstCurve2d();
-        const Handle(Geom2d_Curve)& aC2D2=aIC.SecondCurve2d();
-        //
-        if (aC3D.IsNull()) {
-          continue;
-        }
-        const Handle(Geom_BSplineCurve)& aBC=
-          Handle(Geom_BSplineCurve)::DownCast(aC3D);
-        if (aBC.IsNull()) {
-          continue;
-        }
-        //
-        aT1=aBC->FirstParameter();
-        aT2=aBC->LastParameter();
-        //
-        aEps=0.01*(aT2-aT1);
-        dT=(aT2-aT1)/aNbP;
-        for (j=1; j<aNbP; ++j) {
-          aT11=aT1+j*dT;
-          aT12=aT11+dT;
-          aD2=FindMaxSquareDistance(aT11, aT12, aEps, aC3D, aC2D1, aC2D2,
-                                    myHS1, myHS2, myFace1, myFace2, myContext);
-          if (aD2>aD2Max) {
-            aD2Max=aD2;
-          }
-        }
-      }//for (i=1; i<=aNbLin; ++i) {
-      //
-      myTolReached3d=sqrt(aD2Max);
+      aDMax = ComputeTolerance();
+      if (aDMax > 0.) {
+        myTolReached3d = aDMax;
+      }
     }// if (aNbLin) 
   }// if (aType1==GeomAbs_Cylinder && aType2==GeomAbs_Cylinder) {
   //
@@ -1023,151 +1033,24 @@ static Standard_Boolean isTreatAnalityc(const TopoDS_Face& theF1,
   }// if ((aType1==GeomAbs_Plane && aType2==GeomAbs_Torus) ||
   //
   else if ((aType1==GeomAbs_SurfaceOfRevolution && aType2==GeomAbs_Cylinder) ||
-           (aType2==GeomAbs_SurfaceOfRevolution && aType1==GeomAbs_Cylinder)) {
-    Standard_Integer j, aNbP;
-    Standard_Real aT, aT1, aT2, dT, aD2max, aD2;
+           (aType2==GeomAbs_SurfaceOfRevolution && aType1==GeomAbs_Cylinder) ||
+           (aType1==GeomAbs_Plane && aType2==GeomAbs_Sphere) ||
+           (aType2==GeomAbs_Plane && aType1==GeomAbs_Sphere) ||
+           (aType1==GeomAbs_Plane && aType2==GeomAbs_SurfaceOfExtrusion) ||
+           (aType2==GeomAbs_Plane && aType1==GeomAbs_SurfaceOfExtrusion) ||
+           (aType1==GeomAbs_Plane && aType2==GeomAbs_BSplineSurface) ||
+           (aType2==GeomAbs_Plane && aType1==GeomAbs_BSplineSurface) ||
+           !myApprox) {
     //
-    aNbLin=mySeqOfCurve.Length();
-    aD2max=0.;
-    aNbP=11;
+    Standard_Real aDMax;
     //
-    for (i=1; i<=aNbLin; ++i) {
-      const IntTools_Curve& aIC=mySeqOfCurve(i);
-      const Handle(Geom_Curve)& aC3D=aIC.Curve();
-      const Handle(Geom2d_Curve)& aC2D1=aIC.FirstCurve2d();
-      const Handle(Geom2d_Curve)& aC2D2=aIC.SecondCurve2d();
-      //
-      if (aC3D.IsNull()) {
-        continue;
-      }
-      const Handle(Geom_BSplineCurve)& aBC=
-        Handle(Geom_BSplineCurve)::DownCast(aC3D);
-      if (aBC.IsNull()) {
-        return;
-      }
-      //
-      aT1=aBC->FirstParameter();
-      aT2=aBC->LastParameter();
-      //
-      dT=(aT2-aT1)/(aNbP-1);
-      for (j=0; j<aNbP; ++j) {
-        aT=aT1+j*dT;
-        if (j==aNbP-1) {
-          aT=aT2;
-        }
-        //
-        aD2=MaxSquareDistance(aT, aC3D, aC2D1, aC2D2,
-                              myHS1, myHS2, myFace1, myFace2, myContext);
-        if (aD2>aD2max) {
-          aD2max=aD2;
-        }
-      }//for (j=0; j<aNbP; ++j) {
-     
-    }//for (i=1; i<=aNbLin; ++i) {
-    //
-    aD2=myTolReached3d*myTolReached3d;
-    if (aD2max > aD2) {
-      myTolReached3d=sqrt(aD2max);
+    aDMax = ComputeTolerance();
+    if (aDMax > myTolReached3d) {
+      myTolReached3d = aDMax;
     }
-  }//if((aType1==GeomAbs_SurfaceOfRevolution ...
-  else if ((aType1==GeomAbs_Plane && aType2==GeomAbs_Sphere) ||
-           (aType2==GeomAbs_Plane && aType1==GeomAbs_Sphere)) {
-    Standard_Integer  j, aNbP;
-    Standard_Real aT1, aT2, dT, aD2max, aD2, aEps, aT11, aT12;
-    //
-    aNbLin=mySeqOfCurve.Length();
-    aD2max=0.;
-    aNbP=10;
-    //
-    for (i=1; i<=aNbLin; ++i) {
-      const IntTools_Curve& aIC=mySeqOfCurve(i);
-      const Handle(Geom_Curve)& aC3D=aIC.Curve();
-      const Handle(Geom2d_Curve)& aC2D1=aIC.FirstCurve2d();
-      const Handle(Geom2d_Curve)& aC2D2=aIC.SecondCurve2d();
-      //
-      const Handle(Geom2d_BSplineCurve)& aBC2D1=
-        Handle(Geom2d_BSplineCurve)::DownCast(aC2D1);
-      const Handle(Geom2d_BSplineCurve)& aBC2D2=
-        Handle(Geom2d_BSplineCurve)::DownCast(aC2D2);
-      //
-      if (aBC2D1.IsNull() && aBC2D2.IsNull()) {
-        return;
-      }
-      //
-      if (!aBC2D1.IsNull()) {
-        aT1=aBC2D1->FirstParameter();
-        aT2=aBC2D1->LastParameter();
-      }
-      else {
-        aT1=aBC2D2->FirstParameter();
-        aT2=aBC2D2->LastParameter();
-      }
-      //
-      aEps=0.01*(aT2-aT1);
-      dT=(aT2-aT1)/aNbP;
-      for (j=0; j<aNbP; ++j) {
-        aT11=aT1+j*dT;
-        aT12=aT11+dT;
-        if (j==aNbP-1) {
-          aT12=aT2;
-        }
-        //
-        aD2=FindMaxSquareDistance(aT11, aT12, aEps, aC3D, aC2D1, aC2D2,
-                                  myHS1, myHS2, myFace1, myFace2, myContext);
-        if (aD2>aD2max) {
-          aD2max=aD2;
-        }
-      }//for (j=0; j<aNbP; ++j) {
-     
-    }//for (i=1; i<=aNbLin; ++i) {
-    //
-    aD2=myTolReached3d*myTolReached3d;
-    if (aD2max > aD2) {
-      myTolReached3d=sqrt(aD2max);
-    }
-  }//else if ((aType1==GeomAbs_Plane && aType2==GeomAbs_Sphere) ...
-  else if (!myApprox || bCase1) {
-  //else if (!myApprox) {
-    Standard_Integer aNbP, j;
-    Standard_Real aT1, aT2, dT, aD2, aD2Max, aEps, aT11, aT12;
-    //
-    aD2Max=0.;
-    aNbLin=mySeqOfCurve.Length();
-    //
-    for (i=1; i<=aNbLin; ++i) {
-      const IntTools_Curve& aIC=mySeqOfCurve(i);
-      const Handle(Geom_Curve)& aC3D=aIC.Curve();
-      const Handle(Geom2d_Curve)& aC2D1=aIC.FirstCurve2d();
-      const Handle(Geom2d_Curve)& aC2D2=aIC.SecondCurve2d();
-      //
-      if (aC3D.IsNull()) {
-        continue;
-}
-      const Handle(Geom_BSplineCurve)& aBC=
-        Handle(Geom_BSplineCurve)::DownCast(aC3D);
-      if (aBC.IsNull()) {
-        continue;
-      }
-      //
-      aT1=aBC->FirstParameter();
-      aT2=aBC->LastParameter();
-      //
-      aEps=0.0001*(aT2-aT1);
-      aNbP=11;
-      dT=(aT2-aT1)/aNbP;
-      for (j=1; j<aNbP-1; ++j) {
-        aT11=aT1+j*dT;
-        aT12=aT11+dT;
-        aD2=FindMaxSquareDistance(aT11, aT12, aEps, aC3D, aC2D1, aC2D2,
-                                  myHS1, myHS2, myFace1, myFace2, myContext);
-        if (aD2>aD2Max) {
-          aD2Max=aD2;
-        }
-      }
-    }//for (i=1; i<=aNbLin; ++i) {
-    myTolReached3d=sqrt(aD2Max);
   }
 }
+
 //=======================================================================
 //function : MakeCurve
 //purpose  : 
@@ -4925,118 +4808,112 @@ void RefineVector(gp_Vec2d& aV2D)
   aV2D.SetCoord(aC[0], aC[1]);
 }
 //=======================================================================
-//function : FindMaxSquareDistance
-//purpose  : 
+// Function : FindMaxDistance
+// purpose : 
 //=======================================================================
-Standard_Real FindMaxSquareDistance (const Standard_Real aT1,
-                                     const Standard_Real aT2,
-                                     const Standard_Real aEps,
-                                     const Handle(Geom_Curve)& aC3D,
-                                     const Handle(Geom2d_Curve)& aC2D1,
-                                     const Handle(Geom2d_Curve)& aC2D2,
-                                     const Handle(GeomAdaptor_HSurface)& myHS1,
-                                     const Handle(GeomAdaptor_HSurface)& myHS2,
-                                     const TopoDS_Face& myFace1,
-                                     const TopoDS_Face& myFace2,
-                                     const Handle(IntTools_Context)& myContext)
+Standard_Real FindMaxDistance(const Handle(Geom_Curve)& theCurve,
+                              const Standard_Real theFirst,
+                              const Standard_Real theLast,
+                              const TopoDS_Face& theFace,
+                              const Handle(IntTools_Context)& theContext)
 {
-  Standard_Real aA, aB, aCf, aX1, aX2, aF1, aF2, aX, aF;
+  Standard_Integer aNbS;
+  Standard_Real aT1, aT2, aDt, aD, aDMax, anEps;
   //
-  aCf=1.6180339887498948482045868343656;// =0.5*(1.+sqrt(5.));
-  aA=aT1;
-  aB=aT2;
-  aX1=aB-(aB-aA)/aCf;  
-  aF1=MaxSquareDistance(aX1, 
-                        aC3D, aC2D1, aC2D2, myHS1, myHS2, myFace1, myFace2, myContext);
-  aX2=aA+(aB-aA)/aCf;
-  aF2=MaxSquareDistance(aX2, 
-                        aC3D, aC2D1, aC2D2, myHS1, myHS2, myFace1, myFace2, myContext);
+  aNbS = 11;
+  aDt = (theLast - theFirst) / aNbS;
+  aDMax = 0.;
+  anEps = 1.e-4 * aDt;
   //
-  for(;;) {
+  GeomAPI_ProjectPointOnSurf& aProjPS = theContext->ProjPS(theFace);
+  aT2 = theFirst;
+  for (;;) {
+    aT1 = aT2;
+    aT2 += aDt;
     //
-    if (fabs(aA-aB)<aEps) {
-      aX=0.5*(aA+aB);
-      aF=MaxSquareDistance(aX, 
-                        aC3D, aC2D1, aC2D2, myHS1, myHS2, myFace1, myFace2, myContext);
+    if (aT2 > theLast) {
       break;
     }
-    if (aF1<aF2){
-      aA=aX1;
-      aX1=aX2;
-      aF1=aF2;
-      aX2=aA+(aB-aA)/aCf;
-      aF2=MaxSquareDistance(aX2, 
-                            aC3D, aC2D1, aC2D2, myHS1, myHS2, myFace1, myFace2, myContext);
-      
-    }
-    else {
-      aB=aX2;
-      aX2=aX1;
-      aF2=aF1;
-      aX1=aB-(aB-aA)/aCf; 
-      aF1=MaxSquareDistance(aX1, 
-                            aC3D, aC2D1, aC2D2, myHS1, myHS2, myFace1, myFace2, myContext);
+    //
+    aD = FindMaxDistance(theCurve, aT1, aT2, aProjPS, anEps);
+    if (aD > aDMax) {
+      aDMax = aD;
     }
   }
+  //
+  return aDMax;
+}
+//=======================================================================
+// Function : FindMaxDistance
+// purpose : 
+//=======================================================================
+Standard_Real FindMaxDistance(const Handle(Geom_Curve)& theC,
+                              const Standard_Real theFirst,
+                              const Standard_Real theLast,
+                              GeomAPI_ProjectPointOnSurf& theProjPS,
+                              const Standard_Real theEps)
+{
+  Standard_Real aA, aB, aCf, aX, aX1, aX2, aF1, aF2, aF;
+  //
+  aCf = 0.61803398874989484820458683436564;//(sqrt(5.)-1)/2.;
+  aA = theFirst;
+  aB = theLast;
+  //
+  aX1 = aB - aCf * (aB - aA);
+  aF1 = MaxDistance(theC, aX1, theProjPS);
+  aX2 = aA + aCf * (aB - aA);
+  aF2 = MaxDistance(theC, aX2, theProjPS);
+  //
+  for (;;) {
+    if ((aB - aA) < theEps) {
+      break;
+    }
+    //
+    if (aF1 > aF2) {
+      aB = aX2;
+      aX2 = aX1;
+      aF2 = aF1;
+      aX1 = aB - aCf * (aB - aA); 
+      aF1 = MaxDistance(theC, aX1, theProjPS);
+    }
+    else {
+      aA = aX1;
+      aX1 = aX2;
+      aF1 = aF2;
+      aX2 = aA + aCf * (aB - aA);
+      aF2 = MaxDistance(theC, aX2, theProjPS);
+    }
+  }
+  //
+  aX = 0.5 * (aA + aB);
+  aF = MaxDistance(theC, aX, theProjPS);
+  //
+  if (aF1 > aF) {
+    aF = aF1;
+  }
+  //
+  if (aF2 > aF) {
+    aF = aF2;
+  }
+  //
   return aF;
 }
 //=======================================================================
-//function : MaxSquareDistance
-//purpose  : 
+// Function : MaxDistance
+// purpose : 
 //=======================================================================
-Standard_Real MaxSquareDistance (const Standard_Real aT,
-                                 const Handle(Geom_Curve)& aC3D,
-                                 const Handle(Geom2d_Curve)& aC2D1,
-                                 const Handle(Geom2d_Curve)& aC2D2,
-                                 const Handle(GeomAdaptor_HSurface) myHS1,
-                                 const Handle(GeomAdaptor_HSurface) myHS2,
-                                 const TopoDS_Face& aF1,
-                                 const TopoDS_Face& aF2,
-                                 const Handle(IntTools_Context)& aCtx)
+Standard_Real MaxDistance(const Handle(Geom_Curve)& theC,
+                          const Standard_Real aT,
+                          GeomAPI_ProjectPointOnSurf& theProjPS)
 {
-  Standard_Boolean bIsDone;
-  Standard_Integer i;
-  Standard_Real aU, aV, aD2Max, aD2;
-  gp_Pnt2d aP2D;
-  gp_Pnt aP, aPS;
+  Standard_Real aD;
+  gp_Pnt aP;
   //
-  aD2Max=0.;
+  theC->D0(aT, aP);
+  theProjPS.Perform(aP);
+  aD = theProjPS.NbPoints() ? theProjPS.LowerDistance() : 0.;
   //
-  aC3D->D0(aT, aP);
-  if (aC3D.IsNull()) {
-    return aD2Max;
-  }
-  //
-  for (i=0; i<2; ++i) {
-    const Handle(GeomAdaptor_HSurface)& aGHS=(!i) ? myHS1 : myHS2;
-    const TopoDS_Face &aF=(!i) ? aF1 : aF2;
-    const Handle(Geom2d_Curve)& aC2D=(!i) ? aC2D1 : aC2D2;
-    //
-    if (!aC2D.IsNull()) {
-      aC2D->D0(aT, aP2D);
-      aP2D.Coord(aU, aV);
-      aGHS->D0(aU, aV, aPS);
-      aD2=aP.SquareDistance(aPS);
-      if (aD2>aD2Max) {
-        aD2Max=aD2;
-      }
-    }
-    //
-    GeomAPI_ProjectPointOnSurf& aProjector=aCtx->ProjPS(aF);
-    //
-    aProjector.Perform(aP);
-    bIsDone=aProjector.IsDone();
-    if (bIsDone) {
-      aProjector.LowerDistanceParameters(aU, aV);
-      aGHS->D0(aU, aV, aPS);
-      aD2=aP.SquareDistance(aPS);
-      if (aD2>aD2Max) {
-        aD2Max=aD2;
-      }
-    }
-  }
-  //
-  return aD2Max;
+  return aD;
 }
 
 //=======================================================================
diff --git a/src/IntTools/IntTools_Tools.cdl b/src/IntTools/IntTools_Tools.cdl
index 1c5bcf7390..bdd6941495 100644
--- a/src/IntTools/IntTools_Tools.cdl
+++ b/src/IntTools/IntTools_Tools.cdl
@@ -35,6 +35,8 @@ uses
     Range             from IntTools, 
     SequenceOfCurves  from IntTools,
     Curve   from  Geom,
+    Curve   from  Geom2d,
+    Surface from  Geom,
     State   from  TopAbs,
     Box from Bnd 
     
@@ -237,4 +239,17 @@ is
         theTmax  :out Real from Standard) 
     returns Integer from Standard;  
  
+    ComputeTolerance(myclass;
+        theCurve3D : Curve from Geom;
+        theCurve2D : Curve from Geom2d;
+        theSurf    : Surface from Geom; 
+        theFirst   : Real from Standard;
+        theLast    : Real from Standard;
+        theMaxDist : out Real from Standard; 
+        theMaxPar  : out Real from Standard)
+    returns Boolean from Standard; 
+    ---Purpose:  
+    -- Computes the max distance between points  
+    -- taken from 3D and 2D curves by the same parameter
+ 
 end Tools;
diff --git a/src/IntTools/IntTools_Tools.cxx b/src/IntTools/IntTools_Tools.cxx
index f9c005a544..7045298a71 100644
--- a/src/IntTools/IntTools_Tools.cxx
+++ b/src/IntTools/IntTools_Tools.cxx
@@ -516,6 +516,7 @@ static
 #include <Geom_Parabola.hxx>
 #include <gp_Parab.hxx>
 #include <BndLib_Add3dCurve.hxx>
+#include <BRepLib_CheckCurveOnSurface.hxx>
 //=======================================================================
 //function : ParabolaTolerance
 //purpose  : 
@@ -777,3 +778,30 @@ Standard_Integer IntTools_Tools::SegPln(const gp_Lin& theLin,
   iRet=0; // intersection point
   return iRet;
 }
+
+//=======================================================================
+// Function : ComputeTolerance
+// purpose : 
+//=======================================================================
+Standard_Boolean IntTools_Tools::ComputeTolerance
+  (const Handle(Geom_Curve)& theCurve3D,
+   const Handle(Geom2d_Curve)& theCurve2D,
+   const Handle(Geom_Surface)& theSurf,
+   const Standard_Real theFirst,
+   const Standard_Real theLast,
+   Standard_Real& theMaxDist,
+   Standard_Real& theMaxPar)
+{
+  BRepLib_CheckCurveOnSurface aCS;
+  //
+  aCS.Init(theCurve3D, theCurve2D, theSurf, theFirst, theLast);
+  aCS.Perform();
+  if (!aCS.IsDone()) {
+    return Standard_False;
+  }
+  //
+  theMaxDist = aCS.MaxDistance();
+  theMaxPar  = aCS.MaxParameter();
+  //
+  return Standard_True;
+}
diff --git a/tests/boolean/bsection/N2 b/tests/boolean/bsection/N2
index 95096b00a5..ae758156b8 100644
--- a/tests/boolean/bsection/N2
+++ b/tests/boolean/bsection/N2
@@ -3,4 +3,4 @@ plane p 0 0 0 1 0 0
 mkface f p
 bsection result a f
 
-set length 18981.4
+set length 20674.3
diff --git a/tests/bugs/modalg_2/bug22967 b/tests/bugs/modalg_2/bug22967
index 08115ced95..56b4924896 100755
--- a/tests/bugs/modalg_2/bug22967
+++ b/tests/bugs/modalg_2/bug22967
@@ -1,3 +1,4 @@
+puts "TODO OCC25597 ALL: OCC22967: Faulty"
 puts "============"
 puts "OCC22967"
 puts "============"
diff --git a/tests/bugs/modalg_2/bug23218 b/tests/bugs/modalg_2/bug23218
index 7039562cca..ed7b465c30 100755
--- a/tests/bugs/modalg_2/bug23218
+++ b/tests/bugs/modalg_2/bug23218
@@ -1,3 +1,4 @@
+puts "TODO OCC25597 ALL: OCC23218: Faulty"
 puts "============"
 puts "OCC23218"
 puts "============"
diff --git a/tests/bugs/modalg_5/bug24558 b/tests/bugs/modalg_5/bug24558
index abd263363d..ae0a5584f0 100644
--- a/tests/bugs/modalg_5/bug24558
+++ b/tests/bugs/modalg_5/bug24558
@@ -21,14 +21,14 @@ bbuild result
 
 set square 134338
 
-set nb_v_good 109
-set nb_e_good 189
-set nb_w_good 95
-set nb_f_good 88
-set nb_sh_good 13
+set nb_v_good 108
+set nb_e_good 187
+set nb_w_good 94
+set nb_f_good 87
+set nb_sh_good 11
 set nb_sol_good 5
 set nb_compsol_good 0
 set nb_compound_good 1
-set nb_shape_good 500
+set nb_shape_good 493
 
 set 2dviewer 1
diff --git a/tests/bugs/modalg_5/bug24915 b/tests/bugs/modalg_5/bug24915
index adba90f4e3..62b4e5fb8f 100755
--- a/tests/bugs/modalg_5/bug24915
+++ b/tests/bugs/modalg_5/bug24915
@@ -1,3 +1,4 @@
+puts "TODO OCC25597 ALL: Error: Tolerance is too big!"
 puts "========="
 puts "CR24915"
 puts "========="
diff --git a/tests/bugs/modalg_5/bug25292_32 b/tests/bugs/modalg_5/bug25292_32
index 5dad932ed7..3080aa7d1a 100644
--- a/tests/bugs/modalg_5/bug25292_32
+++ b/tests/bugs/modalg_5/bug25292_32
@@ -1,3 +1,4 @@
+puts "TODO OCC25597 ALL: Error: Tolerance is too big!"
 puts "================"
 puts "OCC25292"
 puts "================"
diff --git a/tests/bugs/modalg_5/bug25597 b/tests/bugs/modalg_5/bug25597
new file mode 100755
index 0000000000..2ae3286319
--- /dev/null
+++ b/tests/bugs/modalg_5/bug25597
@@ -0,0 +1,25 @@
+puts "============"
+puts "OCC25597"
+puts "============"
+puts ""
+######################################################
+# Invalid curve on surface in the result of General Fuse operation
+######################################################
+
+restore [locate_data_file bug25597_c1ext.brep] b1
+restore [locate_data_file bug25597_c2ext.brep] b2
+
+bclearobjects
+bcleartools
+baddobjects b1
+baddtools b2
+bfillds -s
+bbuild r
+
+set info [bopargcheck r]
+
+if { [regexp "to be valid for BOP" ${info}] == 1 } {
+    puts "OK : Created curve is correct"
+} else {
+    puts "Error : Created curve is not correct"
+}