0024650: Wrong intersection curves obtained for a surface of revolution and a plane.

Test case for issue CR24650

src/GeomInt/GeomInt_IntSS_1.cxx

@@ -35,6 +35,7 @@
 #include <IntPatch_Line.hxx>
 #include <IntPatch_WLine.hxx>
 #include <IntPatch_GLine.hxx>
+#include <IntPatch_RLine.hxx>
 #include <IntPatch_ALineToWLine.hxx>
 #include <IntPatch_IType.hxx>
 #include <NCollection_IncAllocator.hxx>
@@ -48,6 +49,7 @@
 #include <GeomAbs_SurfaceType.hxx>
 #include <GeomAbs_CurveType.hxx>
 
+#include <GeomAdaptor.hxx>
 #include <Geom_Curve.hxx>
 #include <Geom_Line.hxx>
 #include <Geom_Parabola.hxx>
@@ -57,6 +59,8 @@
 #include <Geom_Ellipse.hxx>
 #include <Geom_BSplineCurve.hxx>
 
+#include <Geom2dAdaptor.hxx>
+#include <Adaptor2d_HCurve2d.hxx>
 #include <Geom2d_Curve.hxx>
 #include <Geom2d_BSplineCurve.hxx>
 #include <Geom2d_TrimmedCurve.hxx>
@@ -64,7 +68,7 @@
 #include <GeomLib_CheckBSplineCurve.hxx>
 #include <GeomLib_Check2dBSplineCurve.hxx>
 #include <GeomProjLib.hxx>
-
+#include <Approx_CurveOnSurface.hxx>
 
 #include <ElSLib.hxx>
 
@@ -141,6 +145,16 @@ static
   void GetQuadric(const Handle(GeomAdaptor_HSurface)& HS1,
 		  IntSurf_Quadric& quad1);
 
+static
+    void TreatRLine(const Handle(IntPatch_RLine)& theRL, 
+                    const Handle(GeomAdaptor_HSurface)& theHS1, 
+                    const Handle(GeomAdaptor_HSurface)& theHS2, 
+                    const Standard_Boolean theApprox1, 
+                    const Standard_Boolean theApprox2, 
+                    Handle(Geom_Curve)& theC3d,
+                    Handle(Geom2d_Curve)& theC2d1, 
+                    Handle(Geom2d_Curve)& theC2d2, 
+                    Standard_Real& theTolReached);
 
 //
 
@@ -308,13 +322,6 @@ Standard_Real ProjectPointOnSurf::LowerDistance() const
 
       if (!Precision::IsNegativeInfinite(fprm) && 
         !Precision::IsPositiveInfinite(lprm)) {
-	// /cto/900/F1
-	//if (typl == IntPatch_Lin) {
-	  //Standard_Real dPrm=1.3e-6;
-	  //fprm=fprm-dPrm;
-	  //lprm=lprm+dPrm;
-	//}
-	//
           Handle(Geom_TrimmedCurve) aCT3D=new Geom_TrimmedCurve(newc, fprm, lprm);
           sline.Append(aCT3D);
           //
@@ -533,7 +540,7 @@ Standard_Real ProjectPointOnSurf::LowerDistance() const
         }//  end of for (Standard_Integer j=0; j<=17; j++)
       }//  end of else { on regarde si on garde
     }// for (i=1; i<=myLConstruct.NbParts(); i++)
-  }// IntPatch_Circle: IntPatch_Ellipse:
+  }// IntPatch_Circle: IntPatch_Ellipse
   break;
 
                          //########################################  
@@ -1004,6 +1011,45 @@ Standard_Real ProjectPointOnSurf::LowerDistance() const
   break;
 
   case IntPatch_Restriction: 
+    {
+      GeomAbs_SurfaceType typS1 = myHS1->Surface().GetType();
+      GeomAbs_SurfaceType typS2 = myHS2->Surface().GetType();
+      Standard_Boolean isAnalS1 = Standard_False;
+      switch (typS1)
+      {
+      case GeomAbs_Plane:
+      case GeomAbs_Cylinder:
+      case GeomAbs_Sphere:
+      case GeomAbs_Cone: 
+      case GeomAbs_Torus: isAnalS1 = Standard_True; break;
+      default: break;
+      }
+
+      Standard_Integer isAnalS2 = Standard_False;
+      switch (typS2)
+      {
+      case GeomAbs_Plane:
+      case GeomAbs_Cylinder:
+      case GeomAbs_Sphere:
+      case GeomAbs_Cone: 
+      case GeomAbs_Torus: isAnalS2 = Standard_True; break;
+      default: break;
+      }
+
+      if(!isAnalS1 || !isAnalS2)
+      {
+        Handle(IntPatch_RLine) RL = 
+          Handle(IntPatch_RLine)::DownCast(L);
+        Handle(Geom_Curve) aC3d;
+        Handle(Geom2d_Curve) aC2d1, aC2d2;
+        Standard_Real aTolReached;
+        TreatRLine(RL, myHS1, myHS2, myApprox1, myApprox2, aC3d,
+          aC2d1, aC2d2, aTolReached);
+        sline.Append(aC3d);
+        slineS1.Append(aC2d1);
+        slineS2.Append(aC2d2);
+      }
+    }
     break;
 
   }
@@ -1012,6 +1058,74 @@ Standard_Real ProjectPointOnSurf::LowerDistance() const
 //=======================================================================
 //function : AdjustUPeriodic
 //purpose  : 
+//=======================================================================
+void TreatRLine(const Handle(IntPatch_RLine)& theRL, 
+                const Handle(GeomAdaptor_HSurface)& theHS1, 
+                const Handle(GeomAdaptor_HSurface)& theHS2, 
+                const Standard_Boolean theApprox1, 
+                const Standard_Boolean theApprox2, 
+                Handle(Geom_Curve)& theC3d,
+                Handle(Geom2d_Curve)& theC2d1, 
+                Handle(Geom2d_Curve)& theC2d2, 
+                Standard_Real& theTolReached)
+{
+  Handle(GeomAdaptor_HSurface) aGAHS;
+  Handle(Adaptor2d_HCurve2d) anAHC2d;
+  Standard_Real tf, tl;
+  gp_Lin2d aL;
+  // It is assumed that 2d curve is 2d line (rectangular surface domain)
+  if(theRL->IsArcOnS1())
+  {
+    aGAHS = theHS1;
+    anAHC2d = theRL->ArcOnS1();
+    theRL->ParamOnS1(tf, tl);
+    theC2d1 = Geom2dAdaptor::MakeCurve(theRL->ArcOnS1()->Curve2d());
+  }
+  else if (theRL->IsArcOnS2())
+  {
+    aGAHS = theHS2;
+    anAHC2d = theRL->ArcOnS2();
+    theRL->ParamOnS2(tf, tl);
+    theC2d2 = Geom2dAdaptor::MakeCurve(theRL->ArcOnS2()->Curve2d());
+  }
+  else
+  {
+    return;
+  }
+  //
+  //To provide sameparameter it is necessary to get 3d curve as
+  //approximation of curve on surface.
+  Standard_Integer aMaxDeg = 8;
+  Standard_Integer aMaxSeg = 1000;
+  Approx_CurveOnSurface anApp(anAHC2d, aGAHS, tf, tl, Precision::Confusion(),
+                              GeomAbs_C1, aMaxDeg, aMaxSeg, 
+                              Standard_True, Standard_False);
+  if(anApp.HasResult())
+  {
+    theC3d = anApp.Curve3d();
+    theTolReached = anApp.MaxError3d();
+    Standard_Real aTol = Precision::Confusion();
+    if(theRL->IsArcOnS1() && theApprox2)
+    {
+      Handle(Geom_Surface) aS = GeomAdaptor::MakeSurface(theHS1->Surface());
+      BuildPCurves (tf, tl, aTol, 
+                    aS, theC3d, theC2d2);
+    }
+    if(theRL->IsArcOnS2() && theApprox1)
+    {
+      Handle(Geom_Surface) aS = GeomAdaptor::MakeSurface(theHS2->Surface());
+      BuildPCurves (tf, tl, aTol, 
+                    aS, theC3d, theC2d1);
+    }
+    theTolReached = Max(theTolReached, aTol);
+  }
+  //
+
+}
+//
+//=======================================================================
+//function : AdjustUPeriodic
+//purpose  : 
 //=======================================================================
  void AdjustUPeriodic (const Handle(Geom_Surface)& aS, Handle(Geom2d_Curve)& aC2D)
 {

src/IntPatch/IntPatch_ImpPrmIntersection.cxx

@@ -22,6 +22,7 @@
 #include <IntSurf_PntOn2S.hxx>
 #include <IntSurf_LineOn2S.hxx>
 #include <IntSurf.hxx>
+#include <IntSurf_InteriorPoint.hxx>
 
 #include <Adaptor2d_HCurve2d.hxx>
 #include <IntSurf_PathPoint.hxx>
@@ -529,7 +530,49 @@ void IntPatch_ImpPrmIntersection::Perform (const Handle(Adaptor3d_HSurface)& Sur
     }
   }
   //
+  Standard_Boolean SearchIns = Standard_True;
+  if(Quad.TypeQuadric() == GeomAbs_Plane && solrst.NbSegments() > 0)
+  {
+    //For such kind of cases it is possible that whole surface is on one side of plane,
+    //plane only touches surface and does not cross it,
+    //so no inner points exist.
+    SearchIns = Standard_False;
+    Handle(Adaptor3d_TopolTool) T;
+    if(reversed)
+    {
+      T = D1;
+    }
+    else
+    {
+      T = D2;
+    }
+    Standard_Integer aNbSamples = 0;
+    aNbSamples = T->NbSamples();
+    gp_Pnt2d s2d;
+    gp_Pnt s3d;
+    Standard_Real aValf[1], aUVap[2];
+    math_Vector Valf(aValf,1,1), UVap(aUVap,1,2);
+    T->SamplePoint(1,s2d, s3d);
+    UVap(1)=s2d.X(); 
+    UVap(2)=s2d.Y();
+    Func.Value(UVap,Valf);
+    Standard_Real rvalf = Sign(1.,Valf(1));
+    for(i = 2; i <= aNbSamples; ++i)
+    {
+      D1->SamplePoint(i,s2d, s3d);
+      UVap(1)=s2d.X(); 
+      UVap(2)=s2d.Y();
+      Func.Value(UVap,Valf);
+      if(rvalf * Valf(1) < 0.)
+      {
+        SearchIns = Standard_True;
+        break;
+      }   
+    }
+  }
   // Recherche des points interieurs
+  NbPointIns = 0;
+  if(SearchIns) {
     if (!reversed) {
       if (myIsStartPnt)
         solins.Perform(Func,Surf2,myUStart,myVStart);
@@ -542,11 +585,12 @@ void IntPatch_ImpPrmIntersection::Perform (const Handle(Adaptor3d_HSurface)& Sur
       else
         solins.Perform(Func,Surf1,D1,TolTang);
     }
-  //
     NbPointIns = solins.NbPoints();
     for (i=1; i <= NbPointIns; i++) {
       seqpins.Append(solins.Value(i));
     }
+  }
+  //
   NbPointDep=seqpdep.Length();
   //
   if (NbPointDep || NbPointIns) {
@@ -722,7 +766,8 @@ void IntPatch_ImpPrmIntersection::Perform (const Handle(Adaptor3d_HSurface)& Sur
         wline = new IntPatch_WLine(thelin,Standard_False,trans1,trans2);
 
         if (   iwline->HasFirstPoint() 
-          && iwline->IsTangentAtBegining() == Standard_False) {
+          && iwline->IsTangentAtBegining() == Standard_False) 
+        {
           indfirst = iwline->FirstPointIndex();
           PPoint = seqpdep(indfirst);
           tgline = PPoint.Direction3d();
@@ -784,7 +829,8 @@ void IntPatch_ImpPrmIntersection::Perform (const Handle(Adaptor3d_HSurface)& Sur
             }
           }
         }
-        else if (iwline->IsTangentAtBegining()) {
+        else if (iwline->IsTangentAtBegining()) 
+        {
           gp_Pnt psol = thelin->Value(1).Value();
           thelin->Value(1).ParametersOnS1(U1,V1);
           thelin->Value(1).ParametersOnS2(U2,V2);
@@ -794,7 +840,8 @@ void IntPatch_ImpPrmIntersection::Perform (const Handle(Adaptor3d_HSurface)& Sur
           wline->AddVertex(ptdeb);
           wline->SetFirstPoint(wline->NbVertex());
         }
-        else { 
+        else 
+        { 
           gp_Pnt psol = thelin->Value(1).Value();
           thelin->Value(1).ParametersOnS1(U1,V1);
           thelin->Value(1).ParametersOnS2(U2,V2);
@@ -807,7 +854,8 @@ void IntPatch_ImpPrmIntersection::Perform (const Handle(Adaptor3d_HSurface)& Sur
 
 
         if (   iwline->HasLastPoint() 
-          && iwline->IsTangentAtEnd() == Standard_False) {
+          && iwline->IsTangentAtEnd() == Standard_False) 
+        {
           indlast = iwline->LastPointIndex();
           PPoint = seqpdep(indlast);
           tgline = PPoint.Direction3d().Reversed();
@@ -871,7 +919,8 @@ void IntPatch_ImpPrmIntersection::Perform (const Handle(Adaptor3d_HSurface)& Sur
             }
           }
         }
-        else if (iwline->IsTangentAtEnd()) {
+        else if (iwline->IsTangentAtEnd()) 
+        {
           gp_Pnt psol = thelin->Value(Nbpts).Value();
           thelin->Value(Nbpts).ParametersOnS1(U1,V1);
           thelin->Value(Nbpts).ParametersOnS2(U2,V2);
@@ -881,7 +930,8 @@ void IntPatch_ImpPrmIntersection::Perform (const Handle(Adaptor3d_HSurface)& Sur
           wline->AddVertex(ptfin);
           wline->SetLastPoint(wline->NbVertex());
         }
-        else { 
+        else 
+        { 
           gp_Pnt psol = thelin->Value(Nbpts).Value();
           thelin->Value(Nbpts).ParametersOnS1(U1,V1);
           thelin->Value(Nbpts).ParametersOnS2(U2,V2);
@@ -983,6 +1033,38 @@ void IntPatch_ImpPrmIntersection::Perform (const Handle(Adaptor3d_HSurface)& Sur
   if (NbSegm) {
     for(i=1; i<=NbSegm; i++) {
       thesegm = solrst.Segment(i);  
+      //Check if segment is degenerated
+      if(thesegm.HasFirstPoint() && thesegm.HasLastPoint())
+      {
+        Standard_Real tol2 = Precision::Confusion();
+        tol2 *= tol2;
+        const gp_Pnt& aPf = thesegm.FirstPoint().Value();
+        const gp_Pnt& aPl = thesegm.LastPoint().Value();
+        if(aPf.SquareDistance(aPl) <= tol2)
+        {
+          //segment can be degenerated - check inner point
+          paramf = thesegm.FirstPoint().Parameter();
+          paraml = thesegm.LastPoint().Parameter();
+          gp_Pnt2d _p2d = 
+            thesegm.Curve()->Value(.57735 * paramf + 0.42265 * paraml);
+          gp_Pnt aPm;
+          if(reversed)
+          {
+            Surf1->D0(_p2d.X(), _p2d.Y(), aPm);
+          }
+          else
+          {
+            Surf2->D0(_p2d.X(), _p2d.Y(), aPm);
+          }
+          if(aPm.SquareDistance(aPf) <= tol2)
+          {
+            //Degenerated
+            continue;
+          }
+        }
+      }
+
+
       //----------------------------------------------------------------------      
       // on cree une ligne d intersection contenant uniquement le segment.
       // VOIR POUR LA TRANSITION DE LA LIGNE

tests/bugs/modalg_2/bug497_5

@@ -1,4 +1,3 @@
-
 puts "========="
 puts " OCC497 "
 puts "(case 5)"
@@ -19,5 +18,5 @@ if [catch {bfuse result a_1 a_2 } catch_result] {
 } else {
     puts "OCC497 : function FUSE works without hangs up "
 }
-set square 3233.42
+set square 3280.73
 set 2dviewer 0

tests/bugs/modalg_5/bug24650

@@ -0,0 +1,34 @@
+puts "==========="
+puts "OCC24650"
+puts "==========="
+puts ""
+##############################################################
+# Wrong intersection curves obtained for a surface of revolution and a plane.
+##############################################################
+
+restore [locate_data_file bug24650_fz1365.brep] b1
+restore [locate_data_file bug24650_fz2495.brep] b2
+
+mksurface sb1 b1
+mksurface sb2 b2
+trimu sb1tu sb1 0.242 0.291
+
+set N [intersect i sb1tu sb2]
+
+if { [llength ${N}] != 1} {
+   puts "Error : Wrong intersection curves"
+} else {
+   set info1 [dump i]
+   set type_of_curve "BSplineCurve"
+   if { [regexp ${type_of_curve} ${info1}] != 1 } {
+      puts "Error : Wrong type of intersection curve"
+   } else {
+      puts "OK : Good type of intersection curve"
+   }
+}
+
+smallview
+donly i
+fit
+set only_screen_axo 1
+

tests/bugs/moddata_2/bug254

@@ -1,4 +1,3 @@
-
 puts "================"
 puts "OCC254"
 puts "================"
@@ -13,5 +12,5 @@ checkshape shape3
 common result shape1 shape3
 checkshape result
 
-set length 1231.07
+set length 1194.29
 set 2dviewer 0