diff --git a/src/BOPAlgo/BOPAlgo_PaveFiller_5.cxx b/src/BOPAlgo/BOPAlgo_PaveFiller_5.cxx
index 66aec9f8e8..785b53d5c7 100644
--- a/src/BOPAlgo/BOPAlgo_PaveFiller_5.cxx
+++ b/src/BOPAlgo/BOPAlgo_PaveFiller_5.cxx
@@ -332,13 +332,14 @@ void BOPAlgo_PaveFiller::PerformEF()
             }
             //
             const gp_Pnt& aPnew = BRep_Tool::Pnt(aVnew);
-            if (!myContext->IsValidPointForFace(aPnew, 
-                                                aF, 
-                                                aTolE+aTolF)) {
+            Standard_Real aTolV = BRep_Tool::Tolerance(aVnew);
+            aTolV = Max(aTolV, Max(aTolE, aTolF));
+            //
+            if (!myContext->IsPointInFace(aPnew, aF, aTolV)) {
               continue;
             }
             //
-            aBB.UpdateVertex(aVnew, aTolE);
+            aBB.UpdateVertex(aVnew, aTolV);
             //
             aMIEFC.Add(nF);
             // 1
diff --git a/src/BOPTools/BOPTools_AlgoTools_1.cxx b/src/BOPTools/BOPTools_AlgoTools_1.cxx
index 6fa2144fef..ea06e2d336 100644
--- a/src/BOPTools/BOPTools_AlgoTools_1.cxx
+++ b/src/BOPTools/BOPTools_AlgoTools_1.cxx
@@ -617,7 +617,7 @@ Standard_Real IntersectCurves2d(const gp_Pnt& aPV,
   gp_Pnt2d aP2D;
   //
   aDist = 0.;
-  aTol2d = Precision::Confusion();
+  aTol2d = 1.e-10;//Precision::Confusion();
   //
   const Handle(Geom2d_Curve)& aC2D1=
     BRep_Tool::CurveOnSurface(aE1, aF, aT11, aT12);
@@ -637,13 +637,13 @@ Standard_Real IntersectCurves2d(const gp_Pnt& aPV,
     }
     aNbPnt = aInter.NbPoints();
     if (aNbPnt) {
-      aDist = Precision::Infinite();
+      aDist = -Precision::Infinite();
       for (j = 1; j <= aNbPnt; ++j) {
         aP2D = aInter.Point(j).Value();
         aS->D0(aP2D.X(), aP2D.Y(), aP);
         aD=aPV.SquareDistance(aP);
-        if (aD < aDist) {
-          aDist = aD;
+        if (aD > aDist) {
+          aDist = 1.01 * aD;
         }
       }
     }
diff --git a/src/IntTools/IntTools_Context.cdl b/src/IntTools/IntTools_Context.cdl
index 48d4ca0c1e..7bbcd9b953 100644
--- a/src/IntTools/IntTools_Context.cdl
+++ b/src/IntTools/IntTools_Context.cdl
@@ -188,7 +188,17 @@ is
     --- inside the boundaries of the face aF,
     --- otherwise returns false
     ---
-
+ 
+    IsPointInFace(me:mutable;    
+       aP3D :  Pnt   from  gp; 
+       aF   :  Face  from TopoDS; 
+       aTol :  Real from Standard) 
+    returns Boolean from Standard; 
+    ---Purpose:
+    --- Returns true if the point aP2D is
+    --- inside the boundaries of the face aF,
+    --- otherwise returns false
+    ---
     
     IsPointInOnFace(me:mutable;     
        aF   :  Face   from  TopoDS;
diff --git a/src/IntTools/IntTools_Context.cxx b/src/IntTools/IntTools_Context.cxx
index d83b49c522..71432b39a1 100644
--- a/src/IntTools/IntTools_Context.cxx
+++ b/src/IntTools/IntTools_Context.cxx
@@ -553,6 +553,35 @@ Standard_Boolean IntTools_Context::IsPointInFace
   return Standard_True;
 }
 //=======================================================================
+//function : IsPointInFace
+//purpose  : 
+//=======================================================================
+Standard_Boolean IntTools_Context::IsPointInFace
+  (const gp_Pnt& aP,
+   const TopoDS_Face& aF,
+   const Standard_Real aTol) 
+{
+  Standard_Boolean bIn;
+  Standard_Real aDist;
+  //
+  GeomAPI_ProjectPointOnSurf& aProjector=ProjPS(aF);
+  aProjector.Perform(aP);
+  //
+  bIn = aProjector.IsDone();
+  if (bIn) {
+    aDist = aProjector.LowerDistance();
+    if (aDist < aTol) {
+      Standard_Real U, V;
+      //
+      aProjector.LowerDistanceParameters(U, V);
+      gp_Pnt2d aP2D(U, V);
+      bIn = IsPointInFace(aF, aP2D);
+    }
+  }
+  //
+  return bIn;
+}
+//=======================================================================
 //function : IsPointInOnFace
 //purpose  : 
 //=======================================================================
diff --git a/src/IntTools/IntTools_EdgeEdge.cxx b/src/IntTools/IntTools_EdgeEdge.cxx
index f7ce4a6944..580a271736 100644
--- a/src/IntTools/IntTools_EdgeEdge.cxx
+++ b/src/IntTools/IntTools_EdgeEdge.cxx
@@ -48,11 +48,11 @@ static
   Standard_Real PointBoxDistance(const Bnd_Box& aB,
                                  const gp_Pnt& aP);
 static 
-  void SplitRangeOnSegments(const Standard_Real aT1, 
-                            const Standard_Real aT2,
-                            const Standard_Real theResolution,
-                            const Standard_Integer theNbSeg,
-                            IntTools_SequenceOfRanges& theSegments);
+  Standard_Integer SplitRangeOnSegments(const Standard_Real aT1, 
+                                        const Standard_Real aT2,
+                                        const Standard_Real theResolution,
+                                        const Standard_Integer theNbSeg,
+                                        IntTools_SequenceOfRanges& theSegments);
 static
  Standard_Integer DistPC(const Standard_Real aT1, 
                          const Handle(Geom_Curve)& theC1,
@@ -268,11 +268,9 @@ void IntTools_EdgeEdge::FindSolutions(IntTools_SequenceOfRanges& theRanges1,
   aNb1 = IsClosed(myGeom1, aT11, aT12, myTol1, myRes1) ? 2 : 1;
   aNb2 = 2;
   //
-  SplitRangeOnSegments(aT11, aT12, myRes1, aNb1, aSegments1);
-  SplitRangeOnSegments(aT21, aT22, myRes2, aNb2, aSegments2);
+  aNb1 = SplitRangeOnSegments(aT11, aT12, myRes1, aNb1, aSegments1);
+  aNb2 = SplitRangeOnSegments(aT21, aT22, myRes2, aNb2, aSegments2);
   //
-  aNb1 = aSegments1.Length();
-  aNb2 = aSegments2.Length();
   for (i = 1; i <= aNb1; ++i) {
     const IntTools_Range& aR1 = aSegments1(i);
     for (j = 1; j <= aNb2; ++j) {
@@ -433,8 +431,7 @@ void IntTools_EdgeEdge::FindSolutions(const IntTools_Range& theR1,
   IntTools_Range aR2(aT21, aT22);
   BndBuildBox(myCurve2, aT21, aT22, myTol2, aB2);
   //
-  SplitRangeOnSegments(aT11, aT12, myRes1, 3, aSegments1);
-  aNb1 = aSegments1.Length();
+  aNb1 = SplitRangeOnSegments(aT11, aT12, myRes1, 3, aSegments1);
   for (i = 1; i <= aNb1; ++i) {
     const IntTools_Range& aR1 = aSegments1(i);
     FindSolutions(aR1, aR2, aB2, theRanges1, theRanges2);
@@ -685,37 +682,64 @@ void IntTools_EdgeEdge::FindBestSolution(const Standard_Real aT11,
                                          Standard_Real& aT2)
 {
   Standard_Integer i, aNbS, iErr;
-  Standard_Real aDMin, aD, aCrit, aRes1;
-  Standard_Real aT1x, aT2x, aT1p, aT2p;
-  GeomAPI_ProjectPointOnCurve aProj;
-  IntTools_SequenceOfRanges aSeg1;
+  Standard_Real aDMin, aD, aRes1, aSolCriteria, aTouchCriteria;
+  Standard_Real aT1A, aT1B, aT1Min, aT2Min;
+  Standard_Real aT1Im, aT2Im, aT1Touch;
+  GeomAPI_ProjectPointOnCurve aProjPC;
+  IntTools_SequenceOfRanges aRanges;
+  Standard_Boolean bTouch;
   //
-  aT1 = (aT11 + aT12) * .5;
-  aT2 = (aT21 + aT22) * .5;
-  //
-  aDMin = 100.;
-  aD = 100.;
-  aCrit = 0.;//1.e-16;
+  aDMin = Precision::Infinite();
+  aSolCriteria   = 5.e-16;
+  aTouchCriteria = 5.e-13;
+  bTouch = Standard_False;
+  aT1Touch = aT11;
   //
   aRes1 = Resolution(myCurve1.Curve().Curve(), 
                      myCurve1.GetType(), myResCoeff1, myTol);
   aNbS = 10;
-  SplitRangeOnSegments(aT11, aT12, 3*aRes1, aNbS, aSeg1);
-  aNbS = aSeg1.Length();
+  aNbS = SplitRangeOnSegments(aT11, aT12, 3*aRes1, aNbS, aRanges);
+  //
+  aProjPC.Init(myGeom2, aT21, aT22);
+  //
+  aT1 = (aT11 + aT12) * 0.5;
+  iErr = DistPC(aT1, myGeom1, aSolCriteria, aProjPC, aD, aT2, -1);
+  if (iErr == 1) {
+    aT2 = (aT21 + aT22) * 0.5;
+  }
+  //
+  aT1Im = aT1;
+  aT2Im = aT2;
   //
-  aProj.Init(myGeom2, aT21, aT22);
   for (i = 1; i <= aNbS; ++i) {
-    const IntTools_Range& aR1 = aSeg1(i);
-    aR1.Range(aT1x, aT2x);
+    const IntTools_Range& aR1 = aRanges(i);
+    aR1.Range(aT1A, aT1B);
     //
-    iErr = FindDistPC(aT1x, aT2x, myGeom1, aCrit, myPTol1,
-                      aProj, aD, aT1p, aT2p, Standard_False);
-    if (iErr != 1 && aD < aDMin) {
-      aT1 = aT1p;
-      aT2 = aT2p;
-      aDMin = aD;
-      if (aDMin <= aCrit) {
-        break;
+    aD = myTol;
+    iErr = FindDistPC(aT1A, aT1B, myGeom1, aSolCriteria, myPTol1,
+                      aProjPC, aD, aT1Min, aT2Min, Standard_False);
+    if (iErr != 1) {
+      if (aD < aDMin) {
+        aT1 = aT1Min;
+        aT2 = aT2Min;
+        aDMin = aD;
+      }
+      //
+      if (aD < aTouchCriteria) {
+        if (bTouch) {
+          aT1A = (aT1Touch + aT1Min) * 0.5;
+          iErr = DistPC(aT1A, myGeom1, aTouchCriteria, 
+                        aProjPC, aD, aT2Min, -1);
+          if (aD > aTouchCriteria) {
+            aT1 = aT1Im;
+            aT2 = aT2Im;
+            break;
+          }
+        }
+        else {
+          aT1Touch = aT1Min;
+          bTouch = Standard_True;
+        }
       }
     }
   }
@@ -920,13 +944,14 @@ Standard_Boolean IntTools_EdgeEdge::IsIntersection(const Standard_Real aT11,
     //
     if (((anAngle1 < anAngleCriteria) || ((M_PI - anAngle1) < anAngleCriteria)) ||
         ((anAngle2 < anAngleCriteria) || ((M_PI - anAngle2) < anAngleCriteria))) {
-      GeomAPI_ProjectPointOnCurve aProj;
+      GeomAPI_ProjectPointOnCurve aProjPC;
       Standard_Integer iErr;
-      Standard_Real aD, aT1p, aT2p;
+      Standard_Real aD, aT1Min, aT2Min;
       //
-      aD = 100.;
-      aProj.Init(myGeom2, aT21, aT22);
-      iErr = FindDistPC(aT11, aT12, myGeom1, myTol, myRes1, aProj, aD, aT1p, aT2p, Standard_False);
+      aD = Precision::Infinite();
+      aProjPC.Init(myGeom2, aT21, aT22);
+      iErr = FindDistPC(aT11, aT12, myGeom1, myTol, myRes1, 
+                        aProjPC, aD, aT1Min, aT2Min, Standard_False);
       bRet = (iErr == 2);
     }
   }
@@ -947,7 +972,7 @@ Standard_Integer IntTools_EdgeEdge::CheckCoincidence(const Standard_Real aT11,
   Standard_Integer iErr, aNb, aNb1, i;
   Standard_Real aT1A, aT1B, aT1max, aT2max, aDmax;
   GeomAPI_ProjectPointOnCurve aProjPC;
-  IntTools_SequenceOfRanges aSeg1;
+  IntTools_SequenceOfRanges aRanges;
   //
   iErr  = 0;
   aDmax = -1.;
@@ -955,10 +980,9 @@ Standard_Integer IntTools_EdgeEdge::CheckCoincidence(const Standard_Real aT11,
   //
   // 1. Express evaluation
   aNb = 10; // Number of intervals on the curve #1
-  SplitRangeOnSegments(aT11, aT12, theCurveRes1, aNb, aSeg1);
-  aNb1 = aSeg1.Length();
+  aNb1 = SplitRangeOnSegments(aT11, aT12, theCurveRes1, aNb, aRanges);
   for (i = 1; i < aNb1; ++i) {
-    const IntTools_Range& aR1 = aSeg1(i);
+    const IntTools_Range& aR1 = aRanges(i);
     aR1.Range(aT1A, aT1B);
     //
     iErr = DistPC(aT1B, myGeom1, theCriteria, aProjPC, aDmax, aT2max);
@@ -967,7 +991,7 @@ Standard_Integer IntTools_EdgeEdge::CheckCoincidence(const Standard_Real aT11,
     }
   }
   //
-  // if the ranges in aSeg1 are less than theCurveRes1,
+  // if the ranges in aRanges are less than theCurveRes1,
   // there is no need to do step 2 (deep evaluation)
   if (aNb1 < aNb) {
     return iErr;
@@ -975,7 +999,7 @@ Standard_Integer IntTools_EdgeEdge::CheckCoincidence(const Standard_Real aT11,
   //
   // 2. Deep evaluation
   for (i = 2; i < aNb1; ++i) {
-    const IntTools_Range& aR1 = aSeg1(i);
+    const IntTools_Range& aR1 = aRanges(i);
     aR1.Range(aT1A, aT1B);
     //
     iErr = FindDistPC(aT1A, aT1B, myGeom1, theCriteria, theCurveRes1, 
@@ -1017,12 +1041,14 @@ Standard_Integer FindDistPC(const Standard_Real aT1A,
   aB = aT1B;
   //
   // check bounds
-  iErr = DistPC(aA, theC1, theCriteria, theProjPC, aYP, aT2P, aDmax, aT1max, aT2max, iC);
+  iErr = DistPC(aA, theC1, theCriteria, theProjPC, 
+                aYP, aT2P, aDmax, aT1max, aT2max, iC);
   if (iErr == 2) {
     return iErr;
   }
   //
-  iErr = DistPC(aB, theC1, theCriteria, theProjPC, aYL, aT2L, aDmax, aT1max, aT2max, iC);
+  iErr = DistPC(aB, theC1, theCriteria, theProjPC, 
+                aYL, aT2L, aDmax, aT1max, aT2max, iC);
   if (iErr == 2) {
     return iErr;
   }
@@ -1030,12 +1056,14 @@ Standard_Integer FindDistPC(const Standard_Real aT1A,
   aXP = aA + (aB - aA)*aGS;
   aXL = aB - (aB - aA)*aGS;
   //
-  iErr = DistPC(aXP, theC1, theCriteria, theProjPC, aYP, aT2P, aDmax, aT1max, aT2max, iC);
+  iErr = DistPC(aXP, theC1, theCriteria, theProjPC, 
+                aYP, aT2P, aDmax, aT1max, aT2max, iC);
   if (iErr) {
     return iErr;
   }
   //
-  iErr = DistPC(aXL, theC1, theCriteria, theProjPC, aYL, aT2L, aDmax, aT1max, aT2max, iC);
+  iErr = DistPC(aXL, theC1, theCriteria, theProjPC, 
+                aYL, aT2L, aDmax, aT1max, aT2max, iC);
   if (iErr) {
     return iErr;
   }
@@ -1046,20 +1074,25 @@ Standard_Integer FindDistPC(const Standard_Real aT1A,
       aXL = aXP;
       aYL = aYP;
       aXP = aA + (aB - aA)*aGS;
-      iErr = DistPC(aXP, theC1, theCriteria, theProjPC, aYP, aT2P, aDmax, aT1max, aT2max, iC);
-      if (iErr) {
-        return iErr;
-      }
+      iErr = DistPC(aXP, theC1, theCriteria, theProjPC, 
+                    aYP, aT2P, aDmax, aT1max, aT2max, iC);
     }
     else {
       aB = aXP;
       aXP = aXL;
       aYP = aYL;
       aXL = aB - (aB - aA)*aGS;
-      iErr = DistPC(aXL, theC1, theCriteria, theProjPC, aYL, aT2L, aDmax, aT1max, aT2max, iC);
-      if (iErr) {
-        return iErr;
+      iErr = DistPC(aXL, theC1, theCriteria, theProjPC, 
+                    aYL, aT2L, aDmax, aT1max, aT2max, iC);
+    }
+    //
+    if (iErr) {
+      if ((iErr == 2) && !bMaxDist) {
+        aXP = (aA + aB) * 0.5;
+        DistPC(aXP, theC1, theCriteria, theProjPC, 
+               aYP, aT2P, aDmax, aT1max, aT2max, iC);
       }
+      return iErr;
     }
     //
     if ((aB - aA) < theEps) {
@@ -1087,7 +1120,7 @@ Standard_Integer DistPC(const Standard_Real aT1,
   Standard_Integer iErr;
   //
   iErr = DistPC(aT1, theC1, theCriteria, theProjPC, aD, aT2, iC);
-  if (iErr) {
+  if (iErr == 1) {
     return iErr;
   }
   //
@@ -1137,16 +1170,16 @@ Standard_Integer DistPC(const Standard_Real aT1,
 //function : SplitRangeOnSegments
 //purpose  : 
 //=======================================================================
-void SplitRangeOnSegments(const Standard_Real aT1, 
-                          const Standard_Real aT2,
-                          const Standard_Real theResolution,
-                          const Standard_Integer theNbSeg,
-                          IntTools_SequenceOfRanges& theSegments)
+Standard_Integer SplitRangeOnSegments(const Standard_Real aT1, 
+                                      const Standard_Real aT2,
+                                      const Standard_Real theResolution,
+                                      const Standard_Integer theNbSeg,
+                                      IntTools_SequenceOfRanges& theSegments)
 {
   Standard_Real aDiff = aT2 - aT1;
   if (aDiff < theResolution || theNbSeg == 1) {
     theSegments.Append(IntTools_Range(aT1, aT2));
-    return;
+    return 1;
   }
   //
   Standard_Real aDt, aT1x, aT2x, aSeg;
@@ -1172,6 +1205,8 @@ void SplitRangeOnSegments(const Standard_Real aT1,
   //
   IntTools_Range aR(aT1x, aT2);
   theSegments.Append(aR);
+  //
+  return aNbSegments;
 }
 
 //=======================================================================
diff --git a/src/IntTools/IntTools_EdgeFace.cxx b/src/IntTools/IntTools_EdgeFace.cxx
index 85bc50a721..399e8b42bc 100644
--- a/src/IntTools/IntTools_EdgeFace.cxx
+++ b/src/IntTools/IntTools_EdgeFace.cxx
@@ -70,14 +70,15 @@
 
 static
   Standard_Boolean IsCoplanar (const BRepAdaptor_Curve&  ,
-          const BRepAdaptor_Surface& );
+                               const BRepAdaptor_Surface& );
 static
   Standard_Boolean IsRadius (const BRepAdaptor_Curve& aCurve ,
-        const BRepAdaptor_Surface& aSurface);
+                             const BRepAdaptor_Surface& aSurface,
+                             const Standard_Real aCriteria);
 static
   Standard_Integer AdaptiveDiscret (const Standard_Integer iDiscret,
-        const BRepAdaptor_Curve& aCurve ,
-        const BRepAdaptor_Surface& aSurface);
+                                    const BRepAdaptor_Curve& aCurve ,
+                                    const BRepAdaptor_Surface& aSurface);
 
 //=======================================================================
 //function : IntTools_EdgeFace::IntTools_EdgeFace
@@ -1309,14 +1310,14 @@ void IntTools_EdgeFace::Perform()
           if (bIsTouch) {
             aCP.SetType(TopAbs_VERTEX);
             aCP.SetVertexParameter1(aTx);
-            aCP.SetRange1 (aTx, aTx);
+            //aCP.SetRange1 (aTx, aTx);
           }
         }
-        if (aType==TopAbs_VERTEX) {
+        else if (aType==TopAbs_VERTEX) {
           bIsTouch=CheckTouchVertex (aCP, aTx);
           if (bIsTouch) {
             aCP.SetVertexParameter1(aTx);
-            aCP.SetRange1 (aTx, aTx);
+            //aCP.SetRange1 (aTx, aTx);
           }
         }
       }
@@ -1327,7 +1328,7 @@ void IntTools_EdgeFace::Perform()
     if (aCType==GeomAbs_Circle && aSType==GeomAbs_Plane) {
       Standard_Boolean bIsCoplanar, bIsRadius;
       bIsCoplanar=IsCoplanar(myC, myS);
-      bIsRadius=IsRadius(myC, myS);
+      bIsRadius=IsRadius(myC, myS, myCriteria);
       if (!bIsCoplanar && !bIsRadius) {
         for (i=1; i<=aNb; i++) {
           IntTools_CommonPrt& aCP=mySeqOfCommonPrts(i);
@@ -1337,7 +1338,14 @@ void IntTools_EdgeFace::Perform()
             if (bIsTouch) {
               aCP.SetType(TopAbs_VERTEX);
               aCP.SetVertexParameter1(aTx);
-              aCP.SetRange1 (aTx, aTx);
+              //aCP.SetRange1 (aTx, aTx);
+            }
+          }
+          else if (aType==TopAbs_VERTEX) {
+            bIsTouch=CheckTouchVertex (aCP, aTx);
+            if (bIsTouch) {
+              aCP.SetVertexParameter1(aTx);
+              //aCP.SetRange1 (aTx, aTx);
             }
           }
         }
@@ -1483,8 +1491,9 @@ Standard_Boolean IsCoplanar (const BRepAdaptor_Curve& aCurve ,
 //function :  IsRadius
 //purpose  : 
 //=======================================================================
-Standard_Boolean IsRadius (const BRepAdaptor_Curve& aCurve ,
-                           const BRepAdaptor_Surface& aSurface)
+Standard_Boolean IsRadius (const BRepAdaptor_Curve& aCurve,
+                           const BRepAdaptor_Surface& aSurface,
+                           const Standard_Real aCriteria)
 {
   Standard_Boolean bFlag=Standard_False;
 
@@ -1500,7 +1509,7 @@ Standard_Boolean IsRadius (const BRepAdaptor_Curve& aCurve ,
     Standard_Real aR=aCirc.Radius();
     gp_Pln aPln=aSurface.Plane();
     Standard_Real aD=aPln.Distance(aCenter);
-    if (fabs (aD-aR) < 1.e-7) {
+    if (fabs (aD-aR) < aCriteria) {
       return !bFlag;
     }
   }
diff --git a/tests/bugs/modalg_5/bug24286 b/tests/bugs/modalg_5/bug24286
index e8981aed04..aa631d3aec 100644
--- a/tests/bugs/modalg_5/bug24286
+++ b/tests/bugs/modalg_5/bug24286
@@ -9,6 +9,8 @@ puts ""
 restore [locate_data_file bug24286_pipeFiss.brep] b1
 restore [locate_data_file bug24286_shellFiss.brep] b2
 
+bfuzzyvalue 5.e-6
+
 bclearobjects
 bcleartools
 baddobjects b1 b2
diff --git a/tests/bugs/modalg_5/bug24646_1 b/tests/bugs/modalg_5/bug24646_1
new file mode 100644
index 0000000000..18c060a419
--- /dev/null
+++ b/tests/bugs/modalg_5/bug24646_1
@@ -0,0 +1,33 @@
+puts "============"
+puts "OCC24646"
+puts "============"
+puts ""
+######################################################
+# Wrong result done by Boolean Operation algorithm
+######################################################
+
+restore [locate_data_file bug24646_b1.brep] b1
+restore [locate_data_file bug24646_b2.brep] b2
+
+bop b1 b2
+bopfuse result
+
+bopcheck result
+
+set square 8.76332e+06
+
+set nbshapes_expected "
+Number of shapes in shape
+ VERTEX    : 88
+ EDGE      : 180
+ WIRE      : 86
+ FACE      : 86
+ SHELL     : 1
+ SOLID     : 1
+ COMPSOLID : 0
+ COMPOUND  : 1
+ SHAPE     : 443
+"
+checknbshapes result ${nbshapes_expected} 1 "Result done by Boolean Operation algorithm"
+
+set 2dviewer 1
diff --git a/tests/bugs/modalg_5/bug24646_2 b/tests/bugs/modalg_5/bug24646_2
new file mode 100644
index 0000000000..aba94ebca4
--- /dev/null
+++ b/tests/bugs/modalg_5/bug24646_2
@@ -0,0 +1,31 @@
+puts "============"
+puts "OCC24646"
+puts "============"
+puts ""
+######################################################
+# Wrong result done by Boolean Operation algorithm
+######################################################
+
+restore [locate_data_file bug24646_b3.brep] b1
+restore [locate_data_file bug24646_b4.brep] b2
+
+bop b1 b2
+bopfuse result
+
+set square 3.05204e+06
+
+set nbshapes_expected "
+Number of shapes in shape
+ VERTEX    : 22
+ EDGE      : 36
+ WIRE      : 17
+ FACE      : 14
+ SHELL     : 1
+ SOLID     : 1
+ COMPSOLID : 0
+ COMPOUND  : 1
+ SHAPE     : 92
+"
+checknbshapes result ${nbshapes_expected} 1 "Result done by Boolean Operation algorithm"
+
+set 2dviewer 1
diff --git a/tests/bugs/modalg_5/bug24646_3 b/tests/bugs/modalg_5/bug24646_3
new file mode 100644
index 0000000000..6350a5bd4e
--- /dev/null
+++ b/tests/bugs/modalg_5/bug24646_3
@@ -0,0 +1,39 @@
+puts "============"
+puts "OCC24646"
+puts "============"
+puts ""
+######################################################
+# Wrong result done by Boolean Operation algorithm
+######################################################
+
+restore [locate_data_file bug24646_b1.brep] b1
+restore [locate_data_file bug24646_b2.brep] b2
+
+nurbsconvert b1 b1
+nurbsconvert b2 b2
+
+bclearobjects
+bcleartools
+baddobjects b1
+baddtools b2
+
+bfillds
+bbuild result
+
+set square 8.82625e+06
+
+set nbshapes_expected "
+Number of shapes in shape
+ VERTEX    : 88
+ EDGE      : 180
+ WIRE      : 88
+ FACE      : 88
+ SHELL     : 2
+ SOLID     : 2
+ COMPSOLID : 0
+ COMPOUND  : 1
+ SHAPE     : 449
+"
+checknbshapes result ${nbshapes_expected} 1 "Result done by Boolean Operation algorithm"
+
+set 2dviewer 1
diff --git a/tests/bugs/modalg_5/bug24646_4 b/tests/bugs/modalg_5/bug24646_4
new file mode 100644
index 0000000000..fdfd963ea7
--- /dev/null
+++ b/tests/bugs/modalg_5/bug24646_4
@@ -0,0 +1,39 @@
+puts "============"
+puts "OCC24646"
+puts "============"
+puts ""
+######################################################
+# Wrong result done by Boolean Operation algorithm
+######################################################
+
+restore [locate_data_file bug24646_b3.brep] b1
+restore [locate_data_file bug24646_b4.brep] b2
+
+nurbsconvert b1 b1
+nurbsconvert b2 b2
+
+bclearobjects
+bcleartools
+baddobjects b1
+baddtools b2
+
+bfillds
+bbuild result
+
+set square 3.1449e+06
+
+set nbshapes_expected "
+Number of shapes in shape
+ VERTEX    : 22
+ EDGE      : 38
+ WIRE      : 20
+ FACE      : 18
+ SHELL     : 3
+ SOLID     : 3
+ COMPSOLID : 0
+ COMPOUND  : 1
+ SHAPE     : 105
+"
+checknbshapes result ${nbshapes_expected} 1 "Result done by Boolean Operation algorithm"
+
+set 2dviewer 1
diff --git a/tests/bugs/modalg_5/bug25319_1 b/tests/bugs/modalg_5/bug25319_1
index 1ee8ada52a..893fda7ca9 100644
--- a/tests/bugs/modalg_5/bug25319_1
+++ b/tests/bugs/modalg_5/bug25319_1
@@ -1,5 +1,3 @@
-puts "TODO OCC25319 ALL: Faulty shapes in variables faulty_1 to faulty_"
-
 puts "================"
 puts "OCC25319"
 puts "================"
diff --git a/tests/bugs/modalg_5/bug25319_2 b/tests/bugs/modalg_5/bug25319_2
index 794cbcb5b6..601d0f5435 100644
--- a/tests/bugs/modalg_5/bug25319_2
+++ b/tests/bugs/modalg_5/bug25319_2
@@ -1,5 +1,3 @@
-puts "TODO OCC25319 ALL: Faulty shapes in variables faulty_1 to faulty_"
-
 puts "================"
 puts "OCC25319"
 puts "================"