diff --git a/src/BOPTest/BOPTest_BOPCommands.cxx b/src/BOPTest/BOPTest_BOPCommands.cxx
index aa269f17fc..39aba2fab1 100644
--- a/src/BOPTest/BOPTest_BOPCommands.cxx
+++ b/src/BOPTest/BOPTest_BOPCommands.cxx
@@ -35,6 +35,7 @@
 
 #include <IntTools_FaceFace.hxx>
 #include <IntTools_Curve.hxx>
+#include <IntTools_PntOn2Faces.hxx>
 
 #include <BOPCol_ListOfShape.hxx>
 
@@ -641,9 +642,9 @@ Standard_Integer bopcurves (Draw_Interpretor& di,
   const TopoDS_Face& aF2=*(TopoDS_Face*)(&S2);
   //
   Standard_Boolean aToApproxC3d, aToApproxC2dOnS1, aToApproxC2dOnS2, anIsDone;
-  Standard_Integer i, aNbCurves;
+  Standard_Integer i, aNbCurves, aNbPoints;
   Standard_Real anAppTol, aTolR;
-  TCollection_AsciiString aNm("c_");
+  TCollection_AsciiString aNm("c_"), aNp("p_");
   //
   anAppTol = 0.0000001;
   aToApproxC3d = Standard_True;
@@ -684,71 +685,92 @@ Standard_Integer bopcurves (Draw_Interpretor& di,
   //
   aFF.PrepareLines3D(Standard_False);
   const IntTools_SequenceOfCurves& aSCs=aFF.Lines();
-
+  const IntTools_SequenceOfPntOn2Faces& aSPs = aFF.Points();
+  //
+  aNbCurves = aSCs.Length();
+  aNbPoints = aSPs.Length();
+  if (!aNbCurves && !aNbPoints) {
+    di << " has no 3d curves\n";
+    di << " has no 3d points\n";
+    return 0;
+  }
   //
   aTolR=aFF.TolReached3d();
   di << "Tolerance Reached=" << aTolR << "\n";
-
-  aNbCurves=aSCs.Length();
-  if (!aNbCurves) {
-    di << " has no 3d curve\n";
-    return 0;
-  }
-  else
-  {
+  //
+  // curves
+  if (aNbCurves) {
     di << aNbCurves << " curve(s) found.\n";
-  }
-
-  for (i=1; i<=aNbCurves; i++) {
-    const IntTools_Curve& anIC=aSCs(i);
-
-    Handle (Geom_Curve)  aC3D = anIC.Curve();
-
-    if (aC3D.IsNull()) {
-      di << " has Null 3d curve# " << i << "\n";
-      continue;
-    }
-
-    TCollection_AsciiString anIndx(i), aNmx;
-    aNmx = aNm + anIndx;
-
-    Standard_CString nameC = aNmx.ToCString();
-
-    DrawTrSurf::Set(nameC, aC3D);
-    di << nameC << " ";
     //
-    Handle(Geom2d_Curve) aPC1 = anIC.FirstCurve2d();
-    Handle(Geom2d_Curve) aPC2 = anIC.SecondCurve2d();
-    //
-    if (!aPC1.IsNull() || !aPC2.IsNull()) {
-      di << "(";
-      //
-      if (!aPC1.IsNull()) {
-        TCollection_AsciiString pc1N("c2d1_"), pc1Nx;
-        pc1Nx = pc1N + anIndx;
-        Standard_CString nameC2d1 = pc1Nx.ToCString();
-        //
-        DrawTrSurf::Set(nameC2d1, aPC1);
-        di << nameC2d1;
+    for (i=1; i<=aNbCurves; i++) {
+      const IntTools_Curve& anIC=aSCs(i);
+
+      Handle (Geom_Curve)  aC3D = anIC.Curve();
+
+      if (aC3D.IsNull()) {
+        di << " has Null 3d curve# " << i << "\n";
+        continue;
       }
+
+      TCollection_AsciiString anIndx(i), aNmx;
+      aNmx = aNm + anIndx;
+
+      Standard_CString nameC = aNmx.ToCString();
+
+      DrawTrSurf::Set(nameC, aC3D);
+      di << nameC << " ";
       //
-      if (!aPC2.IsNull()) {
-        TCollection_AsciiString pc2N("c2d2_"), pc2Nx;
-        pc2Nx = pc2N + anIndx;
-        Standard_CString nameC2d2 = pc2Nx.ToCString();
-        //
-        DrawTrSurf::Set(nameC2d2, aPC2);
+      Handle(Geom2d_Curve) aPC1 = anIC.FirstCurve2d();
+      Handle(Geom2d_Curve) aPC2 = anIC.SecondCurve2d();
+      //
+      if (!aPC1.IsNull() || !aPC2.IsNull()) {
+        di << "(";
         //
         if (!aPC1.IsNull()) {
-          di << ", ";
+          TCollection_AsciiString pc1N("c2d1_"), pc1Nx;
+          pc1Nx = pc1N + anIndx;
+          Standard_CString nameC2d1 = pc1Nx.ToCString();
+          //
+          DrawTrSurf::Set(nameC2d1, aPC1);
+          di << nameC2d1;
         }
-        di << nameC2d2;
+        //
+        if (!aPC2.IsNull()) {
+          TCollection_AsciiString pc2N("c2d2_"), pc2Nx;
+          pc2Nx = pc2N + anIndx;
+          Standard_CString nameC2d2 = pc2Nx.ToCString();
+          //
+          DrawTrSurf::Set(nameC2d2, aPC2);
+          //
+          if (!aPC1.IsNull()) {
+            di << ", ";
+          }
+          di << nameC2d2;
+        }
+        di << ") ";
       }
-      di << ") ";
     }
+    di << "\n";
+  }
+  //
+  // points
+  if (aNbPoints) {
+    di << aNbPoints << " point(s) found.\n";
+    //
+    for (i = 1; i <= aNbPoints; i++) {
+      const IntTools_PntOn2Faces& aPi = aSPs(i);
+      const gp_Pnt& aP = aPi.P1().Pnt();
+      //
+      TCollection_AsciiString anIndx(i), aNmx;
+      aNmx = aNp + anIndx;
+      Standard_CString nameP = aNmx.ToCString();
+      //
+      DrawTrSurf::Set(nameP, aP);
+      di << nameP << " ";
+    }
+    di << "\n";
   }
   //
-  di << "\n";
   return 0;
 }
 //=======================================================================
diff --git a/src/IntTools/IntTools_FaceFace.cxx b/src/IntTools/IntTools_FaceFace.cxx
index 8880e66498..8fe4b73a95 100644
--- a/src/IntTools/IntTools_FaceFace.cxx
+++ b/src/IntTools/IntTools_FaceFace.cxx
@@ -766,6 +766,7 @@ void IntTools_FaceFace::Perform(const TopoDS_Face& aF1,
     }
 
     // Points
+    Standard_Boolean bValid2D1, bValid2D2;
     Standard_Real U1,V1,U2,V2;
     IntTools_PntOnFace aPntOnF1, aPntOnF2;
     IntTools_PntOn2Faces aPntOn2Faces;
@@ -776,6 +777,19 @@ void IntTools_FaceFace::Perform(const TopoDS_Face& aF1,
       const IntSurf_PntOn2S& aISPnt=myIntersector.Point(i).PntOn2S();
       const gp_Pnt& aPnt=aISPnt.Value();
       aISPnt.Parameters(U1,V1,U2,V2);
+      //
+      // check the validity of the intersection point for the faces
+      bValid2D1 = myContext->IsPointInOnFace(myFace1, gp_Pnt2d(U1, V1));
+      if (!bValid2D1) {
+        continue;
+      }
+      //
+      bValid2D2 = myContext->IsPointInOnFace(myFace2, gp_Pnt2d(U2, V2));
+      if (!bValid2D2) {
+        continue;
+      }
+      //
+      // add the intersection point
       aPntOnF1.Init(myFace1, aPnt, U1, V1);
       aPntOnF2.Init(myFace2, aPnt, U2, V2);
       //
diff --git a/tests/bugs/modalg_6/bug26316 b/tests/bugs/modalg_6/bug26316
new file mode 100755
index 0000000000..dea361cc0e
--- /dev/null
+++ b/tests/bugs/modalg_6/bug26316
@@ -0,0 +1,48 @@
+puts "============"
+puts "OCC26316"
+puts "============"
+puts ""
+###############################
+## Infinite modification of the shapes in General Fuse operation
+###############################
+
+restore [locate_data_file bug26316_shape1.brep] b1
+restore [locate_data_file bug26316_shape2.brep] b2
+
+bclearobjects
+bcleartools
+baddobjects b1
+baddtools b2
+bfillds
+bbuild r
+ 
+explode r so
+ 
+bclearobjects
+bcleartools
+baddobjects r_1
+baddtools r_2
+bfillds
+bbuild r1
+
+explode r1 so
+
+set info1 [compare r_1 r1_1]
+
+set info2 [compare r_2 r1_2]
+
+if { [regexp "same shapes" $info1] == 1  } {
+    puts "OK: r_1 and r1_1 shapes are equal"
+} else {
+    puts "Error : r_1 and r1_1 shapes are not equal"
+}
+
+if { [regexp "same shapes" $info2] == 1  } {
+    puts "OK: r_2 and r1_2 shapes are equal"
+} else {
+    puts "Error : r_2 and r1_2 shapes are not equal"
+}
+
+smallview
+fit
+xwd ${imagedir}/${test_image}.png