diff --git a/src/IntCurveSurface/IntCurveSurface_Inter.gxx b/src/IntCurveSurface/IntCurveSurface_Inter.gxx
index e875d94308..3038684199 100644
--- a/src/IntCurveSurface/IntCurveSurface_Inter.gxx
+++ b/src/IntCurveSurface/IntCurveSurface_Inter.gxx
@@ -306,204 +306,65 @@ void IntCurveSurface_Inter::DoNewBounds(
 //purpose  : Decompose la surface si besoin est 
 //=======================================================================
 void IntCurveSurface_Inter::Perform(const TheCurve&   curve,
-				    const TheSurface& surface) {
+                                    const TheSurface& surface)
+{
   ResetFields();
   done = Standard_True;
   Standard_Integer NbUOnS = TheSurfaceTool::NbUIntervals(surface,GeomAbs_C2);
   Standard_Integer NbVOnS = TheSurfaceTool::NbVIntervals(surface,GeomAbs_C2);
   Standard_Real U0,U1,V0,V1; 
-  
-  if(NbUOnS > 1) { 
+
+  if(NbUOnS > 1)
+  {
     TColStd_Array1OfReal TabU(1,NbUOnS+1);
     TheSurfaceTool::UIntervals(surface,TabU,GeomAbs_C2);
-    for(Standard_Integer iu = 1;iu <= NbUOnS; iu++) {
+    for(Standard_Integer iu = 1;iu <= NbUOnS; iu++)
+    {
       U0 = TabU.Value(iu);
       U1 = TabU.Value(iu+1);
-      if(NbVOnS > 1) { 
-	TColStd_Array1OfReal TabV(1,NbVOnS+1);
-	TheSurfaceTool::VIntervals(surface,TabV,GeomAbs_C2);
-	for(Standard_Integer iv = 1;iv <= NbVOnS; iv++) {
-	  V0 = TabV.Value(iv);
-	  V1 = TabV.Value(iv+1);
-	  Perform(curve,surface,U0,V0,U1,V1);
-	}
+      if(NbVOnS > 1)
+      {
+        TColStd_Array1OfReal TabV(1,NbVOnS+1);
+        TheSurfaceTool::VIntervals(surface,TabV,GeomAbs_C2);
+        for(Standard_Integer iv = 1;iv <= NbVOnS; iv++)
+        {
+          // More than one interval on U and V param space.
+          V0 = TabV.Value(iv);
+          V1 = TabV.Value(iv+1);
+          Perform(curve,surface,U0,V0,U1,V1);
+        }
       }
-      else { //------ composite en U     non composite en V
-	V0 = TheSurfaceTool::FirstVParameter(surface);
-	V1 = TheSurfaceTool::LastVParameter(surface);
-	Perform(curve,surface,U0,V0,U1,V1);
+      else
+      {
+        // More than one interval only on U param space.
+        V0 = TheSurfaceTool::FirstVParameter(surface);
+        V1 = TheSurfaceTool::LastVParameter(surface);
+        Perform(curve,surface,U0,V0,U1,V1);
       }
     }
   }
-  else if(NbVOnS > 1) { //---------- non composite en U  composite en V
+  else if(NbVOnS > 1)
+  {
+    // More than one interval only on V param space.
     U0 = TheSurfaceTool::FirstUParameter(surface);
     U1 = TheSurfaceTool::LastUParameter(surface);   
     TColStd_Array1OfReal TabV(1,NbVOnS+1);
     TheSurfaceTool::VIntervals(surface,TabV,GeomAbs_C2);
-    for(Standard_Integer iv = 1;iv <= NbVOnS; iv++) {
+    for(Standard_Integer iv = 1;iv <= NbVOnS; iv++)
+    {
       V0 = TabV.Value(iv);
       V1 = TabV.Value(iv+1);
       Perform(curve,surface,U0,V0,U1,V1);
     }
-  } 
-  else {
+  }
+  else
+  {
+    // One interval on U and V param space.
     V0 = TheSurfaceTool::FirstVParameter(surface);
     V1 = TheSurfaceTool::LastVParameter(surface);
     U0 = TheSurfaceTool::FirstUParameter(surface);
-    U1 = TheSurfaceTool::LastUParameter(surface); 
+    U1 = TheSurfaceTool::LastUParameter(surface);
 
-
-#if 0
-    //-- jgv patch (from)
-    //Computing of local bounds
-    Standard_Real LocalU0 = U0, LocalU1 = U1, LocalV0 = V0, LocalV1 = V1;
-    Standard_Real Umin = RealLast(), Vmin = RealLast(), Umax = RealFirst(), Vmax = RealFirst();
-    Bnd_Box CurveBox;
-    Standard_Integer i, j, k;
-    //Making GeomAdaptor_Curve
-    Standard_Real f = TheCurveTool::FirstParameter(curve);
-    Standard_Real l = TheCurveTool::LastParameter(curve);
-    GeomAbs_CurveType CurveType = TheCurveTool::GetType(curve);
-    Handle(Geom_Curve) theCurve;
-    switch (CurveType)
-      {
-      case GeomAbs_Line:
-	theCurve = new Geom_Line( TheCurveTool::Line(curve) );
-	break;
-      case GeomAbs_Circle:
-	theCurve = new Geom_Circle( TheCurveTool::Circle(curve) );
-	break;
-      case GeomAbs_Ellipse:
-	theCurve = new Geom_Ellipse( TheCurveTool::Ellipse(curve) );
-	break;
-      case GeomAbs_Hyperbola:
-	theCurve = new Geom_Hyperbola( TheCurveTool::Hyperbola(curve) );
-	break;
-      case GeomAbs_Parabola:
-	theCurve = new Geom_Parabola( TheCurveTool::Parabola(curve) );
-	break;
-      case GeomAbs_BezierCurve:
-	theCurve = TheCurveTool::Bezier(curve);
-	break;
-      case GeomAbs_BSplineCurve:
-	theCurve = TheCurveTool::BSpline(curve);
-	break;
-      }
-    if (!theCurve.IsNull())
-      {
-	GeomAdaptor_Curve GACurve( theCurve, f, l );
-	BndLib_Add3dCurve::Add( GACurve, Precision::Confusion(), CurveBox );
-      }
-    else
-      {
-	Standard_Integer nbp = TheCurveTool::NbSamples(curve,f,l);
-	Standard_Real delta = (f-l)/nbp;
-	for (i = 0; i <= nbp; i++)
-	  {
-	    gp_Pnt aPnt = TheCurveTool::Value(curve, f + i*delta);
-	    CurveBox.Add(aPnt);
-	  }
-      }
-    Standard_Real X[2], Y[2], Z[2];
-    CurveBox.Get( X[0], Y[0], Z[0], X[1], Y[1], Z[1] );
-    Standard_Real Ures = TheSurfaceTool::UResolution(surface, Precision::Confusion());
-    Standard_Real Vres = TheSurfaceTool::VResolution(surface, Precision::Confusion());
-    //Making GeomAdaptor_Surface
-    GeomAbs_SurfaceType SurfaceType = TheSurfaceTool::GetType(surface);
-    Handle(Geom_Surface) theSurface;
-    switch (SurfaceType)
-      {
-      case GeomAbs_Plane:
-	theSurface = new Geom_Plane( TheSurfaceTool::Plane(surface) );
-	break;
-      case GeomAbs_Cylinder:
-	theSurface = new Geom_CylindricalSurface( TheSurfaceTool::Cylinder(surface) );
-	break;
-      case GeomAbs_Cone:
-	theSurface = new Geom_ConicalSurface( TheSurfaceTool::Cone(surface) );
-	break;
-      case GeomAbs_Torus:
-	theSurface = new Geom_ToroidalSurface( TheSurfaceTool::Torus(surface) );
-	break;
-      case GeomAbs_Sphere:
-	theSurface = new Geom_SphericalSurface( TheSurfaceTool::Sphere(surface) );
-	break;
-      case GeomAbs_BezierSurface:
-	theSurface = TheSurfaceTool::Bezier(surface);
-	break;
-      case GeomAbs_BSplineSurface:
-	theSurface = TheSurfaceTool::BSpline(surface);
-	break;
-      case GeomAbs_SurfaceOfRevolution:
-	{
-	  gp_Ax1 Axis = TheSurfaceTool::AxeOfRevolution(surface);
-	  Handle(Adaptor3d_HCurve) AdBC = TheSurfaceTool::BasisCurve(surface);
-	  Handle(Geom_Curve) BC = GetCurve(AdBC);
-	  if (!BC.IsNull())
-	    theSurface = new Geom_SurfaceOfRevolution( BC, Axis );
-	  break;
-	}
-      case GeomAbs_SurfaceOfExtrusion:
-	{
-	  gp_Dir Direction = TheSurfaceTool::Direction(surface);
-	  Handle(Adaptor3d_HCurve) AdBC = TheSurfaceTool::BasisCurve(surface);
-	  Handle(Geom_Curve) BC = GetCurve(AdBC);
-	  if (!BC.IsNull())
-	    theSurface = new Geom_SurfaceOfLinearExtrusion( BC, Direction );
-	  break;
-	}
-      case GeomAbs_OffsetSurface:
-	{
-	  Standard_Real OffsetValue = TheSurfaceTool::OffsetValue(surface);
-	  Handle(Adaptor3d_HSurface) AdBS = TheSurfaceTool::BasisSurface(surface);
-	  Handle(Geom_Surface) BS = GetSurface(AdBS);
-	  if (!BS.IsNull())
-	    theSurface = new Geom_OffsetSurface( BS, OffsetValue );
-	  break;
-	}
-      }
-    if (!theSurface.IsNull())
-      {
-	GeomAdaptor_Surface GASurface( theSurface );
-	Extrema_ExtPS Projector;
-	Projector.Initialize( GASurface, U0, U1, V0, V1, Ures, Vres );
-	for (i = 0; i <= 1; i++)
-	  for (j = 0; j <= 1; j++)
-	    for (k = 0; k <= 1; k++)
-	      {
-		gp_Pnt aPoint( X[i], Y[j], Z[k] );
-		Projector.Perform( aPoint );
-		if (Projector.IsDone() && Projector.NbExt() > 0)
-		  {
-		    Standard_Real mindist = RealLast();
-		    Standard_Integer minind, ind;
-		    for (ind = 1; ind <= Projector.NbExt(); ind++)
-		      {
-			Standard_Real dist = Projector.Value(ind);
-			if (dist < mindist)
-			  {
-			    mindist = dist;
-			    minind  = ind;
-			  }
-		      }
-		    Extrema_POnSurf pons = Projector.Point(minind);
-		    Standard_Real U, V;
-		    pons.Parameter( U, V );
-		    if (U < Umin) Umin = U;
-		    if (U > Umax) Umax = U;
-		    if (V < Vmin) Vmin = V;
-		    if (V > Vmax) Vmax = V;
-		  }
-	      }
-	Umin -= Ures; Umax += Ures; Vmin -= Vres; Vmax += Vres;
-	if (Umin > U0 && Umin <= U1) LocalU0 = Umin;
-	if (Umax < U1 && Umax >= U0) LocalU1 = Umax;
-	if (Vmin > V0 && Vmin <= V1) LocalV0 = Vmin;
-	if (Vmax < V1 && Vmax >= V0) LocalV1 = Vmax;
-	U0 = LocalU0; U1 = LocalU1; V0 = LocalV0; V1 = LocalV1;
-      }
-    //-- jgv patch (to)
-#endif
     Perform(curve,surface,U0,V0,U1,V1);
   }
 }
@@ -512,36 +373,49 @@ void IntCurveSurface_Inter::Perform(const TheCurve&   curve,
 //purpose  : 
 //=======================================================================
 void IntCurveSurface_Inter::Perform(const TheCurve&   curve,
-				    const TheSurface& surface,
-				    const Standard_Real U1,const Standard_Real V1,
-				    const Standard_Real U2,const Standard_Real V2) { 
+                                    const TheSurface& surface,
+                                    const Standard_Real U1,const Standard_Real V1,
+                                    const Standard_Real U2,const Standard_Real V2)
+{
+  // Protection from double type overflow.
+  // This may happen inside square magnitude computation based on normal,
+  // which was computed on bound parameteres (bug26525).
+  Standard_Real UU1 = U1, UU2 = U2, VV1 = V1, VV2 = V2;
+  if (U1 < -1.0e50)
+    UU1 = -1.0e50;
+  if (U2 >  1.0e50)
+    UU2 = 1.0e50;
+  if (V1 < -1.0e50)
+    VV1 = -1.0e50;
+  if (V2 >  1.0e50)
+    VV2 = 1.0e50;
 
   GeomAbs_CurveType CurveType = TheCurveTool::GetType(curve);
   
   switch(CurveType) { 
   case GeomAbs_Line:  
     { 
-      PerformConicSurf(TheCurveTool::Line(curve),curve,surface,U1,V1,U2,V2);
+      PerformConicSurf(TheCurveTool::Line(curve),curve,surface,UU1,VV1,UU2,VV2);
       break;
     }
   case GeomAbs_Circle: 
     { 
-      PerformConicSurf(TheCurveTool::Circle(curve),curve,surface,U1,V1,U2,V2);
+      PerformConicSurf(TheCurveTool::Circle(curve),curve,surface,UU1,VV1,UU2,VV2);
       break;
     }
   case GeomAbs_Ellipse: 
     { 
-      PerformConicSurf(TheCurveTool::Ellipse(curve),curve,surface,U1,V1,U2,V2);
+      PerformConicSurf(TheCurveTool::Ellipse(curve),curve,surface,UU1,VV1,UU2,VV2);
       break;
     }    
   case GeomAbs_Parabola:  
     {
-      PerformConicSurf(TheCurveTool::Parabola(curve),curve,surface,U1,V1,U2,V2);
+      PerformConicSurf(TheCurveTool::Parabola(curve),curve,surface,UU1,VV1,UU2,VV2);
       break;
     }
   case GeomAbs_Hyperbola: 
     {
-      PerformConicSurf(TheCurveTool::Hyperbola(curve),curve,surface,U1,V1,U2,V2);
+      PerformConicSurf(TheCurveTool::Hyperbola(curve),curve,surface,UU1,VV1,UU2,VV2);
       break;
     }
   default:    
@@ -568,7 +442,7 @@ void IntCurveSurface_Inter::Perform(const TheCurve&   curve,
 
 // 	    IntCurveSurface_ThePolygon  polygon(curve,u1,u2,TheCurveTool::NbSamples(curve,u1,u2));
 	    IntCurveSurface_ThePolygon  polygon(curve, aPars->Array1());
-	    InternalPerform(curve,polygon,surface,U1,V1,U2,V2);
+	    InternalPerform(curve,polygon,surface,UU1,VV1,UU2,VV2);
 	  }
 	}
 	else {
@@ -583,7 +457,7 @@ void IntCurveSurface_Inter::Perform(const TheCurve&   curve,
 
 // 	  IntCurveSurface_ThePolygon       polygon(curve,TheCurveTool::NbSamples(curve,u1,u2));
 	  IntCurveSurface_ThePolygon  polygon(curve, aPars->Array1());
-	  InternalPerform(curve,polygon,surface,U1,V1,U2,V2);
+	  InternalPerform(curve,polygon,surface,UU1,VV1,UU2,VV2);
 	}
       }
       else {   //-- la surface est une quadrique 
diff --git a/src/QABugs/QABugs_19.cxx b/src/QABugs/QABugs_19.cxx
index cead6a144f..9105b2e01d 100644
--- a/src/QABugs/QABugs_19.cxx
+++ b/src/QABugs/QABugs_19.cxx
@@ -4133,6 +4133,90 @@ static Standard_Integer OCC26313(Draw_Interpretor& di,Standard_Integer n,const c
   return 0;
 }
 
+//=======================================================================
+//function : OCC26525
+//purpose  : check number of intersection points
+//=======================================================================
+#include <BRepAdaptor_Curve.hxx>
+#include <BRepAdaptor_Surface.hxx>
+#include <BRepAdaptor_HCurve.hxx>
+#include <BRepAdaptor_HSurface.hxx>
+#include <IntCurveSurface_HInter.hxx>
+Standard_Integer OCC26525 (Draw_Interpretor& di, 
+                           Standard_Integer n, 
+                           const char** a)
+{
+  TopoDS_Shape aS1, aS2;
+  TopoDS_Edge aE; 
+  TopoDS_Face aF;
+
+  if (n<4)
+  {
+    di << " use OCC26525 r edge face \n";
+    return 1;
+  }
+
+  aS1 = DBRep::Get(a[2]);
+  aS2 = DBRep::Get(a[3]);
+
+  if (aS1.IsNull() || aS2.IsNull())
+  {
+    di << " Null shapes are not allowed \n";
+    return 0;
+  }
+  if (aS1.ShapeType()!=TopAbs_EDGE)
+  {
+    di << " Shape" << a[2] << " should be of type EDGE\n";
+    return 0;
+  }
+  if (aS2.ShapeType()!=TopAbs_FACE)
+  {
+    di << " Shape" << a[3] << " should be of type FACE\n";
+    return 0;
+  }
+
+  aE=TopoDS::Edge(aS1);
+  aF=TopoDS::Face(aS2);
+
+  char buf[128];
+  Standard_Boolean bIsDone;
+  Standard_Integer i, aNbPoints;
+  Standard_Real aU, aV, aT;
+  gp_Pnt aP;
+  BRepAdaptor_Curve aBAC;
+  BRepAdaptor_Surface aBAS;
+  IntCurveSurface_TransitionOnCurve aTC;
+  IntCurveSurface_HInter aHInter;
+
+  aBAC.Initialize(aE);
+  aBAS.Initialize(aF);
+
+  Handle(BRepAdaptor_HCurve) aHBAC=new BRepAdaptor_HCurve(aBAC);
+  Handle(BRepAdaptor_HSurface) aHBAS = new BRepAdaptor_HSurface(aBAS);
+
+  aHInter.Perform(aHBAC, aHBAS);
+  bIsDone=aHInter.IsDone();
+  if (!bIsDone)
+  {
+    di << " intersection is not done\n";
+    return 0;
+  }
+
+  aNbPoints=aHInter.NbPoints();
+  sprintf (buf, " Number of intersection points found: %d", aNbPoints);
+  di <<  buf << "\n";
+  for (i=1; i<=aNbPoints; ++i)
+  {
+    const IntCurveSurface_IntersectionPoint& aIP=aHInter.Point(i);
+    aIP.Values(aP, aU, aV, aT, aTC);
+    //
+    sprintf (buf, "point %s_%d %lg %lg %lg  ", a[1], i, aP.X(), aP.Y(), aP.Z());
+    di << buf << "\n";
+  }
+
+  return 0;
+}
+
 void QABugs::Commands_19(Draw_Interpretor& theCommands) {
   const char *group = "QABugs";
 
@@ -4222,6 +4306,7 @@ void QABugs::Commands_19(Draw_Interpretor& theCommands) {
                    __FILE__, OCC26462, group);
 
   theCommands.Add ("OCC26313", "OCC26313 result shape", __FILE__, OCC26313, group);
+  theCommands.Add ("OCC26525", "OCC26525 result edge face ", __FILE__, OCC26525, group);
 
   return;
 }
diff --git a/tests/bugs/modalg_6/bug26525_1 b/tests/bugs/modalg_6/bug26525_1
new file mode 100644
index 0000000000..796d609424
--- /dev/null
+++ b/tests/bugs/modalg_6/bug26525_1
@@ -0,0 +1,52 @@
+puts "================"
+puts "OCC26525"
+puts "================"
+puts ""
+#######################################################################
+# Wrong result obtained by curve / surface intersection algorithm.
+#######################################################################
+
+pload QAcommands
+
+restore [locate_data_file bug26525_a.brep] b1
+restore [locate_data_file bug26525_b.brep] b2
+
+mksurface sb1 b1
+trimv sb1t sb1 -30000 30000
+mkface b1t sb1t
+
+explode b2 e
+
+# Case 1. The curve is from the edge b2_3 
+# 1.1 The face b1 is based on untrimmed surface
+
+set log [OCC26525 p b2_3 b1]
+
+regexp {([-0-9]+)} $log full Number
+
+if { $Number == 2} {
+   set tol_abs 0.0001
+   set tol_rel 0.01
+
+   regexp {point p_1 +([-0-9.+eE]+) +([-0-9.+eE]+) +([-0-9.+eE]+)} $log full x1 y1 z1
+
+   set expected_x1 48.4205
+   set expected_y1 -22.5336
+   set expected_z1 82.7431
+
+   checkreal "x1" ${x1} ${expected_x1} ${tol_abs} ${tol_rel}
+   checkreal "y1" ${y1} ${expected_y1} ${tol_abs} ${tol_rel}
+   checkreal "z1" ${z1} ${expected_z1} ${tol_abs} ${tol_rel}
+
+   regexp {point p_2 +([-0-9.+eE]+) +([-0-9.+eE]+) +([-0-9.+eE]+)} $log full x2 y2 z2
+
+   set expected_x2 32.5621
+   set expected_y2 -5.89907
+   set expected_z2 82.7431
+
+   checkreal "x2" ${x2} ${expected_x2} ${tol_abs} ${tol_rel}
+   checkreal "y2" ${y2} ${expected_y2} ${tol_abs} ${tol_rel}
+   checkreal "z2" ${z2} ${expected_z2} ${tol_abs} ${tol_rel}
+} else {
+   puts "Error: Bad Number of intersection points"
+}
diff --git a/tests/bugs/modalg_6/bug26525_2 b/tests/bugs/modalg_6/bug26525_2
new file mode 100644
index 0000000000..2e0c01808e
--- /dev/null
+++ b/tests/bugs/modalg_6/bug26525_2
@@ -0,0 +1,52 @@
+puts "================"
+puts "OCC26525"
+puts "================"
+puts ""
+#######################################################################
+# Wrong result obtained by curve / surface intersection algorithm.
+#######################################################################
+
+pload QAcommands
+
+restore [locate_data_file bug26525_a.brep] b1
+restore [locate_data_file bug26525_b.brep] b2
+
+mksurface sb1 b1
+trimv sb1t sb1 -30000 30000
+mkface b1t sb1t
+
+explode b2 e
+
+# Case 1. The curve is from the edge b2_3 
+# 1.2 The face b1 is based on trimmed surface
+
+set log [OCC26525 pt b2_3 b1t]
+
+regexp {([-0-9]+)} $log full Number
+
+if { $Number == 2} {
+   set tol_abs 0.0001
+   set tol_rel 0.01
+
+   regexp {point pt_1 +([-0-9.+eE]+) +([-0-9.+eE]+) +([-0-9.+eE]+)} $log full x1 y1 z1
+
+   set expected_x1 48.4205
+   set expected_y1 -22.5336
+   set expected_z1 82.7431
+
+   checkreal "x1" ${x1} ${expected_x1} ${tol_abs} ${tol_rel}
+   checkreal "y1" ${y1} ${expected_y1} ${tol_abs} ${tol_rel}
+   checkreal "z1" ${z1} ${expected_z1} ${tol_abs} ${tol_rel}
+
+   regexp {point pt_2 +([-0-9.+eE]+) +([-0-9.+eE]+) +([-0-9.+eE]+)} $log full x2 y2 z2
+
+   set expected_x2 32.5621
+   set expected_y2 -5.89907
+   set expected_z2 82.7431
+
+   checkreal "x2" ${x2} ${expected_x2} ${tol_abs} ${tol_rel}
+   checkreal "y2" ${y2} ${expected_y2} ${tol_abs} ${tol_rel}
+   checkreal "z2" ${z2} ${expected_z2} ${tol_abs} ${tol_rel}
+} else {
+   puts "Error: Bad Number of intersection points"
+}
diff --git a/tests/bugs/modalg_6/bug26525_3 b/tests/bugs/modalg_6/bug26525_3
new file mode 100644
index 0000000000..bbc864ec31
--- /dev/null
+++ b/tests/bugs/modalg_6/bug26525_3
@@ -0,0 +1,52 @@
+puts "================"
+puts "OCC26525"
+puts "================"
+puts ""
+#######################################################################
+# Wrong result obtained by curve / surface intersection algorithm.
+#######################################################################
+
+pload QAcommands
+
+restore [locate_data_file bug26525_a.brep] b1
+restore [locate_data_file bug26525_b.brep] b2
+
+mksurface sb1 b1
+trimv sb1t sb1 -30000 30000
+mkface b1t sb1t
+
+explode b2 e
+
+# Case 2. The curve is from the edge: b2_1
+# 2.1 The face b1 is based on untrimmed surface
+
+set log [OCC26525 p b2_1 b1]
+
+regexp {([-0-9]+)} $log full Number
+
+if { $Number == 2} {
+   set tol_abs 0.0001
+   set tol_rel 0.01
+
+   regexp {point p_1 +([-0-9.+eE]+) +([-0-9.+eE]+) +([-0-9.+eE]+)} $log full x1 y1 z1
+
+   set expected_x1 39.0504
+   set expected_y1 -12.8696
+   set expected_z1 82.6099
+
+   checkreal "x1" ${x1} ${expected_x1} ${tol_abs} ${tol_rel}
+   checkreal "y1" ${y1} ${expected_y1} ${tol_abs} ${tol_rel}
+   checkreal "z1" ${z1} ${expected_z1} ${tol_abs} ${tol_rel}
+
+   regexp {point p_2 +([-0-9.+eE]+) +([-0-9.+eE]+) +([-0-9.+eE]+)} $log full x2 y2 z2
+
+   set expected_x2 43.1172
+   set expected_y2 -17.16
+   set expected_z2 82.6048
+
+   checkreal "x2" ${x2} ${expected_x2} ${tol_abs} ${tol_rel}
+   checkreal "y2" ${y2} ${expected_y2} ${tol_abs} ${tol_rel}
+   checkreal "z2" ${z2} ${expected_z2} ${tol_abs} ${tol_rel}
+} else {
+   puts "Error: Bad Number of intersection points"
+}
diff --git a/tests/bugs/modalg_6/bug26525_4 b/tests/bugs/modalg_6/bug26525_4
new file mode 100644
index 0000000000..49a03fd711
--- /dev/null
+++ b/tests/bugs/modalg_6/bug26525_4
@@ -0,0 +1,53 @@
+puts "================"
+puts "OCC26525"
+puts "================"
+puts ""
+#######################################################################
+# Wrong result obtained by curve / surface intersection algorithm.
+#######################################################################
+
+pload QAcommands
+
+restore [locate_data_file bug26525_a.brep] b1
+restore [locate_data_file bug26525_b.brep] b2
+
+mksurface sb1 b1
+trimv sb1t sb1 -30000 30000
+mkface b1t sb1t
+
+explode b2 e
+
+# Case 2. The curve is from the edge: b2_1
+#
+# 2.2 The face b1 is based on trimmed surface
+
+set log [OCC26525 pt b2_1 b1t]
+
+regexp {([-0-9]+)} $log full Number
+
+if { $Number == 2} {
+   set tol_abs 0.0001
+   set tol_rel 0.01
+
+   regexp {point pt_1 +([-0-9.+eE]+) +([-0-9.+eE]+) +([-0-9.+eE]+)} $log full x1 y1 z1
+
+   set expected_x1 39.0504
+   set expected_y1 -12.8696
+   set expected_z1 82.6099
+
+   checkreal "x1" ${x1} ${expected_x1} ${tol_abs} ${tol_rel}
+   checkreal "y1" ${y1} ${expected_y1} ${tol_abs} ${tol_rel}
+   checkreal "z1" ${z1} ${expected_z1} ${tol_abs} ${tol_rel}
+
+   regexp {point pt_2 +([-0-9.+eE]+) +([-0-9.+eE]+) +([-0-9.+eE]+)} $log full x2 y2 z2
+
+   set expected_x2 43.1172
+   set expected_y2 -17.16
+   set expected_z2 82.6048
+
+   checkreal "x2" ${x2} ${expected_x2} ${tol_abs} ${tol_rel}
+   checkreal "y2" ${y2} ${expected_y2} ${tol_abs} ${tol_rel}
+   checkreal "z2" ${z2} ${expected_z2} ${tol_abs} ${tol_rel}
+} else {
+   puts "Error: Bad Number of intersection points"
+}