Bug fix: points located on the arc where wrongly detected by arc/arc intersector.

diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.cxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.cxx
index 0e2ec6ce9799212b2dc3d674b7b3dcfcc9c305c5..0ca8f1f3dac180f03ac2c1300dbc0c983d23ac8a 100644 (file)
--- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.cxx
+++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.cxx
@@ -135,7 +135,7 @@ bool ArcCArcCIntersector::internalAreColinears(const EdgeArcCircle& a1, const Ed
   distBetweenCenters=Node::distanceBtw2PtSq(centerL,centerB);
   cst=distBetweenCenters/(radiusB*radiusB);
   cst+=radiusL*radiusL/(radiusB*radiusB);
-  return Node::areDoubleEqualsWP(cst,1.,2.);
+  return Node::areDoubleEqualsWPRight(cst,1.,2.);
 }
 
 bool ArcCArcCIntersector::areArcsOverlapped(const EdgeArcCircle& a1, const EdgeArcCircle& a2)
@@ -152,7 +152,7 @@ bool ArcCArcCIntersector::areArcsOverlapped(const EdgeArcCircle& a1, const EdgeA
   delete merge;
   //
   tmp=sqrt(tmp);
-  if(Node::areDoubleEqualsWP(tmp,0.,1/(10*std::max(radiusL,radiusB))))
+  if(Node::areDoubleEqualsWPLeft(tmp,0.,10*std::max(radiusL,radiusB)))
     return Node::areDoubleEquals(radiusL,radiusB);
   double phi=EdgeArcCircle::GetAbsoluteAngleOfNormalizedVect((centerL[0]-centerB[0])/tmp,(centerL[1]-centerB[1])/tmp);
   double cst2=2*radiusL*tmp/(radiusB*radiusB);
@@ -167,7 +167,7 @@ bool ArcCArcCIntersector::areArcsOverlapped(const EdgeArcCircle& a1, const EdgeA
   if(EdgeArcCircle::IsIn2Pi(angle0L,angleL,a))
     cmpContainer[sizeOfCmpContainer++]=cst-cst2;
   a=*std::max_element(cmpContainer,cmpContainer+sizeOfCmpContainer);
-  return Node::areDoubleEqualsWP(a,1.,2.);
+  return Node::areDoubleEqualsWPRight(a,1.,2.);
 }
 
 void ArcCArcCIntersector::areOverlappedOrOnlyColinears(const Bounds *whereToFind, bool& obviousNoIntersection, bool& areOverlapped)
@@ -192,22 +192,28 @@ void ArcCArcCIntersector::areOverlappedOrOnlyColinears(const Bounds *whereToFind
     }
 }
 
+/**
+ Heart of the algorithm for arc/arc intersection.
+ See http://mathworld.wolfram.com/Circle-CircleIntersection.html
+ The computation is done in the coordinate system where Ox is the line between the 2 circle centers.
+*/
 std::list< IntersectElement > ArcCArcCIntersector::getIntersectionsCharacteristicVal() const
 {
   std::list< IntersectElement > ret;
   const double *center1=getE1().getCenter();
   const double *center2=getE2().getCenter();
   double radius1=getE1().getRadius(); double radius2=getE2().getRadius();
-  double d1_1=(_dist*_dist-radius2*radius2+radius1*radius1)/(2.*_dist);
+  double d1_1=(_dist*_dist-radius2*radius2+radius1*radius1)/(2.*_dist);  // computation of 'x' on wolfram
   double u[2];//u is normalized vector from center1 to center2.
   u[0]=(center2[0]-center1[0])/_dist; u[1]=(center2[1]-center1[1])/_dist;
-  double d1_1y=EdgeArcCircle::SafeSqrt(radius1*radius1-d1_1*d1_1);
+  double d1_1y=EdgeArcCircle::SafeSqrt(radius1*radius1-d1_1*d1_1);      // computation of 'y' on wolfram
   double angleE1=EdgeArcCircle::NormalizeAngle(getE1().getAngle0()+getE1().getAngle());
   double angleE2=EdgeArcCircle::NormalizeAngle(getE2().getAngle0()+getE2().getAngle());
   if(!Node::areDoubleEquals(d1_1y,0))
     {
       //2 intersections
       double v1[2],v2[2];
+      // coming back to our coordinate system:
       v1[0]=u[0]*d1_1-u[1]*d1_1y; v1[1]=u[1]*d1_1+u[0]*d1_1y;
       v2[0]=u[0]*d1_1+u[1]*d1_1y; v2[1]=u[1]*d1_1-u[0]*d1_1y;
       Node *node1=new Node(center1[0]+v1[0],center1[1]+v1[1]); node1->declareOn();
@@ -219,16 +225,17 @@ std::list< IntersectElement > ArcCArcCIntersector::getIntersectionsCharacteristi
       v4[0]=center1[0]-center2[0]+v2[0]; v4[1]=center1[1]-center2[1]+v2[1];
       double angle1_2=EdgeArcCircle::GetAbsoluteAngleOfNormalizedVect(v3[0]/radius2,v3[1]/radius2);
       double angle2_2=EdgeArcCircle::GetAbsoluteAngleOfNormalizedVect(v4[0]/radius2,v4[1]/radius2);
+      // Check whether intersection points are exactly ON the other arc or not
+      //   -> the curvilinear distance (=radius*angle) must below eps
+      bool e1_1S=Node::areDoubleEqualsWPLeft(angle1_1,getE1().getAngle0(),radius1);
+      bool e1_1E=Node::areDoubleEqualsWPLeft(angle1_1,angleE1,radius1);
+      bool e1_2S=Node::areDoubleEqualsWPLeft(angle1_2,getE2().getAngle0(),radius1);
+      bool e1_2E=Node::areDoubleEqualsWPLeft(angle1_2,angleE2,radius1);
       //
-      bool e1_1S=Node::areDoubleEqualsWP(angle1_1,getE1().getAngle0(),radius1);
-      bool e1_1E=Node::areDoubleEqualsWP(angle1_1,angleE1,radius1);
-      bool e1_2S=Node::areDoubleEqualsWP(angle1_2,getE2().getAngle0(),radius1);
-      bool e1_2E=Node::areDoubleEqualsWP(angle1_2,angleE2,radius1);
-      //
-      bool e2_1S=Node::areDoubleEqualsWP(angle2_1,getE1().getAngle0(),radius2);
-      bool e2_1E=Node::areDoubleEqualsWP(angle2_1,angleE1,radius2);
-      bool e2_2S=Node::areDoubleEqualsWP(angle2_2,getE2().getAngle0(),radius2);
-      bool e2_2E=Node::areDoubleEqualsWP(angle2_2,angleE2,radius2);
+      bool e2_1S=Node::areDoubleEqualsWPLeft(angle2_1,getE1().getAngle0(),radius2);
+      bool e2_1E=Node::areDoubleEqualsWPLeft(angle2_1,angleE1,radius2);
+      bool e2_2S=Node::areDoubleEqualsWPLeft(angle2_2,getE2().getAngle0(),radius2);
+      bool e2_2E=Node::areDoubleEqualsWPLeft(angle2_2,angleE2,radius2);
       ret.push_back(IntersectElement(angle1_1,angle1_2,e1_1S,e1_1E,e1_2S,e1_2E,node1,_e1,_e2,keepOrder()));
       ret.push_back(IntersectElement(angle2_1,angle2_2,e2_1S,e2_1E,e2_2S,e2_2E,node2,_e1,_e2,keepOrder()));
     }
@@ -240,10 +247,10 @@ std::list< IntersectElement > ArcCArcCIntersector::getIntersectionsCharacteristi
       v2[0]=center1[0]-center2[0]+v1[0]; v2[1]=center1[1]-center2[1]+v1[1];
       double angle0_1=EdgeArcCircle::GetAbsoluteAngleOfNormalizedVect(v1[0]/radius1,v1[1]/radius1);
       double angle0_2=EdgeArcCircle::GetAbsoluteAngleOfNormalizedVect(v2[0]/radius2,v2[1]/radius2);
-      bool e0_1S=Node::areDoubleEqualsWP(angle0_1,getE1().getAngle0(),radius1);
-      bool e0_1E=Node::areDoubleEqualsWP(angle0_1,angleE1,radius1);
-      bool e0_2S=Node::areDoubleEqualsWP(angle0_2,getE2().getAngle0(),radius2);
-      bool e0_2E=Node::areDoubleEqualsWP(angle0_2,angleE2,radius2);
+      bool e0_1S=Node::areDoubleEqualsWPLeft(angle0_1,getE1().getAngle0(),radius1);
+      bool e0_1E=Node::areDoubleEqualsWPLeft(angle0_1,angleE1,radius1);
+      bool e0_2S=Node::areDoubleEqualsWPLeft(angle0_2,getE2().getAngle0(),radius2);
+      bool e0_2E=Node::areDoubleEqualsWPLeft(angle0_2,angleE2,radius2);
       Node *node=new Node(center1[0]+d1_1*u[0],center1[1]+d1_1*u[1]); node->declareOnTangent();
       ret.push_back(IntersectElement(angle0_1,angle0_2,e0_1S,e0_1E,e0_2S,e0_2E,node,_e1,_e2,keepOrder()));
     }
@@ -522,10 +529,7 @@ double EdgeArcCircle::GetAbsoluteAngle(const double *vect, double& normVect)
 
 /*!
  * Given a \b normalized vector defined by (ux,uy) returns its angle in ]-Pi;Pi].
- * So before using this method ux*ux+uy*uy should as much as possible close to 1.
- * This methods is quite time consuming in order to keep as much as possible precision.
- * It is NOT ALWAYS possible to do that only in one call of acos. Sometimes call to asin is necessary
- * due to imperfection of acos near 0. and Pi (cos x ~ 1-x*x/2.)
+ * Actually in the current implementation, the vector does not even need to be normalized ...
  */
 double EdgeArcCircle::GetAbsoluteAngleOfNormalizedVect(double ux, double uy)
 {

diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.hxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.hxx
index 31f5fafc130bfaf1314026ce581ce886a7893907..19559f9bc08f1bb7f4f7afb8b921b889797d5d19 100644 (file)
--- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.hxx
+++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.hxx
@@ -126,11 +126,13 @@ namespace INTERP_KERNEL
     void fillGlobalInfoAbs2(const std::map<INTERP_KERNEL::Node *,int>& mapThis, const std::map<INTERP_KERNEL::Node *,int>& mapOther, int offset1, int offset2, double fact, double baryX, double baryY,
                             std::vector<int>& edgesOther, std::vector<double>& addCoo, std::map<INTERP_KERNEL::Node *,int>& mapAddCoo) const;
   protected:
-    //!Value between -2Pi and 2Pi
-    double _angle;
-    //!Value between -Pi and Pi
+    //! Absolute angle where the arc starts. Value between -Pi and Pi
     double _angle0;
+    //! Angular span of the arc. Value between -2Pi and 2Pi
+    double _angle;
+    //! Radius of the arc of circle
     double _radius;
+    //! Center of the arc of circle
     double _center[2];
   };
 }

diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DNode.hxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DNode.hxx
index 970b0f3280d154afa099e068490efb9daa237d72..cbe1e0471924a136bd3eb9c4ceaba2b49e632e66 100644 (file)
--- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DNode.hxx
+++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DNode.hxx
@@ -90,7 +90,9 @@ namespace INTERP_KERNEL
     static double norm(const double *vect) { return sqrt(vect[0]*vect[0]+vect[1]*vect[1]); }
     static bool areDoubleEquals(double a, double b) { return fabs(a-b) < QuadraticPlanarPrecision::getPrecision(); }
     //! idem areDoubleEquals except that precision of comparison is modified.
-    static bool areDoubleEqualsWP(double a, double b, double k) { return fabs(a-b) < k*QuadraticPlanarPrecision::getPrecision(); }
+    static bool areDoubleEqualsWPLeft(double a, double b, double k) { return k*fabs(a-b) < QuadraticPlanarPrecision::getPrecision(); }
+    //! idem areDoubleEquals except that precision of comparison is modified.
+    static bool areDoubleEqualsWPRight(double a, double b, double k) { return fabs(a-b) < k*QuadraticPlanarPrecision::getPrecision(); }
     static double distanceBtw2Pt(const double *a, const double *b) { return sqrt((a[0]-b[0])*(a[0]-b[0])+(a[1]-b[1])*(a[1]-b[1])); }
     static double distanceBtw2PtSq(const double *a, const double *b) { return (a[0]-b[0])*(a[0]-b[0])+(a[1]-b[1])*(a[1]-b[1]); }
     //

diff --git a/src/INTERP_KERNELTest/QuadraticPlanarInterpTest.hxx b/src/INTERP_KERNELTest/QuadraticPlanarInterpTest.hxx
index b5cb15293f30edee2105b4fac16d6a8e71a930ad..87b0ba653e27a42b3c817bf8b3c4bb6a76a0b56d 100644 (file)
--- a/src/INTERP_KERNELTest/QuadraticPlanarInterpTest.hxx
+++ b/src/INTERP_KERNELTest/QuadraticPlanarInterpTest.hxx
@@ -110,6 +110,8 @@ namespace INTERP_TEST
     CPPUNIT_TEST( checkIsInOrOut );
     CPPUNIT_TEST( checkGetMiddleOfPoints );
     CPPUNIT_TEST( checkGetMiddleOfPointsOriented );
+    CPPUNIT_TEST( checkArcArcIntersection1 );
+
     CPPUNIT_TEST_SUITE_END();
   public:  
     void setUp();
@@ -202,6 +204,7 @@ namespace INTERP_TEST
     void checkIsInOrOut();
     void checkGetMiddleOfPoints();
     void checkGetMiddleOfPointsOriented();
+    void checkArcArcIntersection1();
 
   private:
     INTERP_KERNEL::QuadraticPolygon *buildQuadraticPolygonCoarseInfo(const double *coords, const int *conn, int lgth);

diff --git a/src/INTERP_KERNELTest/QuadraticPlanarInterpTest4.cxx b/src/INTERP_KERNELTest/QuadraticPlanarInterpTest4.cxx
index 019eeb84e0cddbeb41c7e337fac8a67368244354..67f0a110496177ff4818ddd26e2a7e48c94cd888 100644 (file)
--- a/src/INTERP_KERNELTest/QuadraticPlanarInterpTest4.cxx
+++ b/src/INTERP_KERNELTest/QuadraticPlanarInterpTest4.cxx
@@ -1682,4 +1682,37 @@ void QuadraticPlanarInterpTest::checkMakePartitionAbs1()
   pol1.dumpInXfigFileWithOther(pol2,"tony.fig");
 }
 
+/**
+ * Arc/arc intersection was buggy: a point was detected OFF when used in a linear segment, but
+ * detected ON when used in an (almost flat) arc of circle.
+ */
+void QuadraticPlanarInterpTest::checkArcArcIntersection1()
+{
+  double eps=1.0e-8;
+  INTERP_KERNEL::QuadraticPlanarPrecision::setPrecision(eps);
+  INTERP_KERNEL::QuadraticPlanarArcDetectionPrecision::setArcDetectionPrecision(eps);
+
+  Node *n0=new Node(6.37533,38.8928);            Node *n3=new Node(6.29194,39.2789);
+  Node *n1=new Node(6.13158,38.8308);            Node *n4=new Node(6.31919,38.7607);
+  Node *n2=new Node(6.25346,38.8618);            Node *n5=new Node(6.38778,39.0241);
+
+  Node *n6=new Node(6.2534549999999998, 38.861800000000002);  // to have a linear edge e1
+
+  //EdgeArcCircle *e1=new EdgeArcCircle(n0, n2, n6, true);  // to have a linear edge e1
+  EdgeArcCircle *e1=new EdgeArcCircle(n0, n2, n1, true);
+  EdgeArcCircle *e2=new EdgeArcCircle(n3, n5, n4, true);
+
+  MergePoints merge;
+  QuadraticPolygon c1,c2;
+  e1->intersectWith(e2,merge,c1,c2);
+  CPPUNIT_ASSERT_EQUAL(2,c1.size());
+  CPPUNIT_ASSERT_EQUAL(2,c2.size());
+  //clean-up
+  n0->decrRef(); n1->decrRef(); n2->decrRef(); n3->decrRef(); n4->decrRef(); n5->decrRef(); n6->decrRef();
+  e1->decrRef(); e2->decrRef();
 }
+
+
+}
+
+