0028394: Not precise extrema solution of line and circle lying in the same plane

If the line is in the circle-plane completely (or parallel to the circle-plane) then extremas and intersections in 2D-space are looked for. These case are pure analytical and solutions will be found precisely.

src/Extrema/Extrema_ExtElC.cxx

@@ -15,20 +15,24 @@
 
 #include <ElCLib.hxx>
 #include <Extrema_ExtElC.hxx>
+#include <Extrema_ExtElC2d.hxx>
 #include <Extrema_ExtPElC.hxx>
 #include <Extrema_POnCurv.hxx>
 #include <gp_Ax1.hxx>
 #include <gp_Ax2.hxx>
 #include <gp_Ax3.hxx>
 #include <gp_Circ.hxx>
+#include <gp_Circ2d.hxx>
 #include <gp_Dir.hxx>
 #include <gp_Elips.hxx>
 #include <gp_Hypr.hxx>
 #include <gp_Lin.hxx>
+#include <gp_Lin2d.hxx>
 #include <gp_Parab.hxx>
 #include <gp_Pln.hxx>
 #include <gp_Pnt.hxx>
 #include <IntAna_QuadQuadGeo.hxx>
+#include <IntAna2d_AnaIntersection.hxx>
 #include <math_DirectPolynomialRoots.hxx>
 #include <math_TrigonometricFunctionRoots.hxx>
 #include <Precision.hxx>
@@ -323,6 +327,90 @@ Extrema_ExtElC::Extrema_ExtElC (const gp_Lin& theC1,
   myNbExt = 1;
   myDone = Standard_True;
 }
+
+//=======================================================================
+//function : PlanarLineCircleExtrema
+//purpose  : 
+//=======================================================================
+Standard_Boolean Extrema_ExtElC::PlanarLineCircleExtrema(const gp_Lin& theLin,
+                                                         const gp_Circ& theCirc)
+{
+  const gp_Dir &aDirC = theCirc.Axis().Direction(),
+               &aDirL = theLin.Direction();
+
+  if (Abs(aDirC.Dot(aDirL)) > Precision::Angular())
+    return Standard_False;
+  
+  //The line is in the circle-plane completely
+  //(or parallel to the circle-plane).
+  //Therefore, we are looking for extremas and
+  //intersections in 2D-space.
+
+  const gp_XYZ &aCLoc = theCirc.Location().XYZ();
+  const gp_XYZ &aDCx = theCirc.Position().XDirection().XYZ(),
+               &aDCy = theCirc.Position().YDirection().XYZ();
+
+  const gp_XYZ &aLLoc = theLin.Location().XYZ();
+  const gp_XYZ &aLDir = theLin.Direction().XYZ();
+
+  const gp_XYZ aVecCL(aLLoc - aCLoc);
+
+  //Center of 2D-circle
+  const gp_Pnt2d aPC(0.0, 0.0);
+
+  gp_Ax22d aCircAxis(aPC, gp_Dir2d(1.0, 0.0), gp_Dir2d(0.0, 1.0));
+  gp_Circ2d aCirc2d(aCircAxis, theCirc.Radius());
+
+  gp_Pnt2d aPL(aVecCL.Dot(aDCx), aVecCL.Dot(aDCy));
+  gp_Dir2d aDL(aLDir.Dot(aDCx), aLDir.Dot(aDCy));
+  gp_Lin2d aLin2d(aPL, aDL);
+
+  // Extremas
+  Extrema_ExtElC2d anExt2d(aLin2d, aCirc2d, Precision::Confusion());
+  //Intersections
+  IntAna2d_AnaIntersection anInters(aLin2d, aCirc2d);
+
+  myDone = anExt2d.IsDone() || anInters.IsDone();
+
+  if (!myDone)
+    return Standard_True;
+
+  const Standard_Integer aNbExtr = anExt2d.NbExt();
+  const Standard_Integer aNbSol = anInters.NbPoints();
+
+  const Standard_Integer aNbSum = aNbExtr + aNbSol;
+  for (Standard_Integer anExtrID = 1; anExtrID <= aNbSum; anExtrID++)
+  {
+    const Standard_Integer aDelta = anExtrID - aNbExtr;
+
+    Standard_Real aLinPar = 0.0, aCircPar = 0.0;
+
+    if (aDelta < 1)
+    {
+      Extrema_POnCurv2d aPLin2d, aPCirc2d;
+      anExt2d.Points(anExtrID, aPLin2d, aPCirc2d);
+      aLinPar = aPLin2d.Parameter();
+      aCircPar = aPCirc2d.Parameter();
+    }
+    else
+    {
+      aLinPar = anInters.Point(aDelta).ParamOnFirst();
+      aCircPar = anInters.Point(aDelta).ParamOnSecond();
+    }
+
+    const gp_Pnt aPOnL(ElCLib::LineValue(aLinPar, theLin.Position())),
+                 aPOnC(ElCLib::CircleValue(aCircPar,
+                                            theCirc.Position(), theCirc.Radius()));
+
+    mySqDist[myNbExt] = aPOnL.SquareDistance(aPOnC);
+    myPoint[myNbExt][0].SetValues(aLinPar, aPOnL);
+    myPoint[myNbExt][1].SetValues(aCircPar, aPOnC);
+    myNbExt++;
+  }
+
+  return Standard_True;
+}
+
 //=======================================================================
 //function : Extrema_ExtElC
 //purpose  : 
@@ -367,6 +455,11 @@ Extrema_ExtElC::Extrema_ExtElC (const gp_Lin& C1,
   myDone = Standard_False;
   myNbExt = 0;
 
+  if (PlanarLineCircleExtrema(C1, C2))
+  {
+    return;
+  }
+
   // Calculate T1 in the reference of the circle ...
   D = C1.Direction();
   D1 = D;

src/Extrema/Extrema_ExtElC.hxx

@@ -130,6 +130,9 @@ public:
 
 protected:
 
+  //! Computes extrema in case when considered line and circle are in one plane
+  Standard_EXPORT Standard_Boolean PlanarLineCircleExtrema(const gp_Lin& C1,
+                                                           const gp_Circ& C2);
 
 
 

tests/bugs/modalg_6/bug28394_1

@@ -0,0 +1,15 @@
+puts "========"
+puts "OCC28394"
+puts "========"
+puts ""
+##############################################
+# Not precise extrema solution of line and circle lying in the same plane
+##############################################
+
+restore [locate_data_file bug28394_edges.brep] e
+explode e
+set anInfo [distmini d e_1 e_2]
+
+if {[dval d_val] > 1.0e-7} {
+  puts "Error: Extrema cannot find minimal distance"
+}

tests/bugs/modalg_6/bug28394_2

@@ -0,0 +1,42 @@
+puts "========"
+puts "OCC28394"
+puts "========"
+puts ""
+##############################################
+# Not precise extrema solution of line and circle lying in the same plane
+##############################################
+
+set GoodNbExtremas 4
+
+circle c1 5 5 10 0 1 1 20
+mkedge e1 c1
+bmirror e1 e1 5 5 10 1 0 0
+mkcurve c1 e1
+
+cvalue c1 0.63 x1 y1 z1
+cvalue c1 5.47 x2 y2 z2
+
+dset dlx x2-x1
+dset dly y2-y1
+dset dlz z2-z1
+line l1 x1 y1+5 z1+5 dlx dly dlz
+
+set extrema_length [ llength [ extrema c1 l1 1 ] ]
+
+# Amount check
+if {${extrema_length} != [ expr 5*${GoodNbExtremas}] } {
+  puts "Error: Number of extremas computed is wrong"
+}
+
+for {set i 1} {${i} <= 4} {incr i} {
+  regexp {Axis   :([-0-9.+eE]+), ([-0-9.+eE]+), ([-0-9.+eE]+)} [dump ext_${i}] full dx dy dz
+  
+  cvalue c1 prm_1_${i} x y z dcx dcy dcz
+  
+  set DPL [ dval ${dx}*dlx+${dy}*dly+${dz}*dlz ]
+  set DPC [ dval ${dx}*dcx+${dy}*dcy+${dz}*dcz ]
+  
+  if { (abs(${DPL}) > 1.0e-12) || (abs(${DPC}) > 1.0e-12) } {
+    puts "Error: extrema ext_${i} was computed wrong (is not the normal to the curves)"
+  }
+}