ArcCSegIntersector fix: tangent cases were not properly handled.

(also removed a duplicated test case)

diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.cxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.cxx
index ba2d601e7a731fb6d5fda37770c49f63e428190e..429cd0cb698226bd493c816c2df38d99b3b4198a 100644 (file)
--- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.cxx
+++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.cxx
@@ -26,6 +26,7 @@
 
 #include <sstream>
 #include <algorithm>
+#include <limits>
 
 using namespace INTERP_KERNEL;
 
@@ -324,21 +325,41 @@ ArcCSegIntersector::ArcCSegIntersector(const EdgeArcCircle& e1, const EdgeLin& e
 {
 }
 
+/**
+  See http://mathworld.wolfram.com/Circle-LineIntersection.html
+  _cross is 'D', the computation is done with the translation to put back the circle at the origin.s
+*/
 void ArcCSegIntersector::areOverlappedOrOnlyColinears(const Bounds *whereToFind, bool& obviousNoIntersection, bool& areOverlapped)
 {
   areOverlapped=false;//No overlapping by construction
   const double *center=getE1().getCenter();
+  const double R = getE1().getRadius();
   _dx=(*(_e2.getEndNode()))[0]-(*(_e2.getStartNode()))[0];
   _dy=(*(_e2.getEndNode()))[1]-(*(_e2.getStartNode()))[1];
   _drSq=_dx*_dx+_dy*_dy;
   _cross=
       ((*(_e2.getStartNode()))[0]-center[0])*((*(_e2.getEndNode()))[1]-center[1])-
       ((*(_e2.getStartNode()))[1]-center[1])*((*(_e2.getEndNode()))[0]-center[0]);
-  _determinant = (getE1().getRadius()*getE1().getRadius()-_cross*_cross/_drSq) / _drSq;
-  if(_determinant > -2*QuadraticPlanarPrecision::getPrecision())//QuadraticPlanarPrecision::getPrecision()*QuadraticPlanarPrecision::getPrecision()*_drSq*_drSq/(2.*_dx*_dx))
+
+  // We need to compute d = R*R-_cross*_cross/_drSq
+  // In terms of numerical precision, this can trigger 'catastrophic cancellation' and is hence better expressed as:
+  double _dr = sqrt(_drSq);
+  double diff = (R-_cross/_dr), add=(R+_cross/_dr);
+  // Ah ah: we will be taking a square root later. If we want the user to be able to use an epsilon finer than 1.0e-8, then we need
+  // to prevent ourselves going below machine precision (typ. 1.0e-16 for double).
+  const double eps_machine = std::numeric_limits<double>::epsilon();
+  diff = fabs(diff) < 10*eps_machine ? 0.0 : diff;
+  add  = fabs(add)  < 10*eps_machine ? 0.0 : add;
+  double d = add*diff;
+  // Compute deltaRoot_div_dr := sqrt(delta)/dr, where delta has the meaning of Wolfram.
+  // Then 2*deltaRoot_div_dr is the distance between the two intersection points of the line with the circle. This is what we compare to eps.
+  // We compute it in such a way that it can be used in tests too (a very negative value means we're far apart from intersection)
+  _deltaRoot_div_dr = Node::sign(d)*sqrt(fabs(d));
+
+  if( 2*_deltaRoot_div_dr > -QuadraticPlanarPrecision::getPrecision())
     obviousNoIntersection=false;
   else
-    obviousNoIntersection=true;   
+    obviousNoIntersection=true;
 }
 
 /*!
@@ -358,9 +379,9 @@ std::list< IntersectElement > ArcCSegIntersector::getIntersectionsCharacteristic
 {
   std::list< IntersectElement > ret;
   const double *center=getE1().getCenter();
-  if(!(fabs(_determinant)<(2.*QuadraticPlanarPrecision::getPrecision())))//QuadraticPlanarPrecision::getPrecision()*QuadraticPlanarPrecision::getPrecision()*_drSq*_drSq/(2.*_dx*_dx))
+  if(!(2*fabs(_deltaRoot_div_dr) < QuadraticPlanarPrecision::getPrecision())) // see comments in areOverlappedOrOnlyColinears()
     {
-      double determinant=EdgeArcCircle::SafeSqrt(_determinant);
+      double determinant=fabs(_deltaRoot_div_dr)/sqrt(_drSq);
       double x1=(_cross*_dy/_drSq+Node::sign(_dy)*_dx*determinant)+center[0];
       double y1=(-_cross*_dx/_drSq+fabs(_dy)*determinant)+center[1];
       Node *intersect1=new Node(x1,y1); intersect1->declareOn();

diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.hxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.hxx
index 436920c60ce91dd1ce0ecea1a4f149e7756adb35..fc32223e20ee2176e8034d486741bff6a02d586b 100644 (file)
--- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.hxx
+++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.hxx
@@ -67,7 +67,7 @@ namespace INTERP_KERNEL
     double _dy;           //!< Y extent of the segment
     double _drSq;         //!< Square of the norm of the seg
     double _cross;        //!< See areOverlappedOrOnlyColinears()
-    double _determinant;  //!< See areOverlappedOrOnlyColinears()
+    double _deltaRoot_div_dr;    //!< See areOverlappedOrOnlyColinears()
   };
 
   class INTERPKERNEL_EXPORT EdgeArcCircle : public Edge

diff --git a/src/MEDCoupling/MEDCouplingUMesh_intersection.cxx b/src/MEDCoupling/MEDCouplingUMesh_intersection.cxx
index 28652a6bbdf68891a18187e479a414b1f6e7fb58..52cd1f1d05ffa0950be2e02afbc7f89b88ac282e 100644 (file)
--- a/src/MEDCoupling/MEDCouplingUMesh_intersection.cxx
+++ b/src/MEDCoupling/MEDCouplingUMesh_intersection.cxx
@@ -900,7 +900,9 @@ int VectorOfCellInfo::getPositionOf(double eps, const MEDCouplingUMesh *mesh) co
   return zeMesh->getCellContainingPoint(barys->begin(),eps);
 }
 
-void VectorOfCellInfo::setMeshAt(std::size_t pos, const MCAuto<MEDCouplingUMesh>& mesh, int istart, int iend, const MCAuto<MEDCouplingUMesh>& mesh1DInCase, const std::vector< std::vector<int> >& edges, const std::vector< std::vector< MCAuto<INTERP_KERNEL::Edge> > >& edgePtrs)
+void VectorOfCellInfo::setMeshAt(std::size_t pos, const MCAuto<MEDCouplingUMesh>& mesh, int istart, int iend,
+                                 const MCAuto<MEDCouplingUMesh>& mesh1DInCase, const std::vector< std::vector<int> >& edges,
+                                 const std::vector< std::vector< MCAuto<INTERP_KERNEL::Edge> > >& edgePtrs)
 {
   get(pos);//to check pos
   bool isFast(pos==0 && _pool.size()==1);
@@ -967,7 +969,7 @@ int VectorOfCellInfo::getZePosOfEdgeGivenItsGlobalId(int pos) const
 
 void VectorOfCellInfo::updateEdgeInfo(int pos, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newLeft, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newRight)
 {
-  get(pos);//to check;
+  get(pos);//to perform the sanity check;
   if(_edge_info.empty())
     return ;
   std::size_t sz(_edge_info.size()-1);

diff --git a/src/MEDCoupling_Swig/MEDCouplingBasicsTest3.py b/src/MEDCoupling_Swig/MEDCouplingBasicsTest3.py
index 89abd42f1e90ba848c1a75374e83edf569306f48..80f5f6bb0b3382240dc2d78d50013cbb20b8dba9 100644 (file)
--- a/src/MEDCoupling_Swig/MEDCouplingBasicsTest3.py
+++ b/src/MEDCoupling_Swig/MEDCouplingBasicsTest3.py
@@ -130,7 +130,7 @@ class MEDCouplingBasicsTest3(unittest.TestCase):
         da.iota(7)
         da.rearrange(5)
         da.setInfoOnComponent(0,"X [m]") ; da.setInfoOnComponent(1,"Y [km]") ; da.setInfoOnComponent(2,"Y [m]")
-        da.setInfoOnComponent(3,"Z [W]") ; da.setInfoOnComponent(4,"ZZ [km]") ; 
+        da.setInfoOnComponent(3,"Z [W]") ; da.setInfoOnComponent(4,"ZZ [km]") ;
         da[:,2]=3
         self.assertEqual([7, 8, 3, 10, 11, 12, 13, 3, 15, 16, 17, 18, 3, 20, 21, 22, 23, 3, 25, 26],da.getValues())
         da.rearrange(1) ; da.iota(7) ; da.rearrange(5)
@@ -186,7 +186,7 @@ class MEDCouplingBasicsTest3(unittest.TestCase):
         da.iota(7)
         da.rearrange(5)
         da.setInfoOnComponent(0,"X [m]") ; da.setInfoOnComponent(1,"Y [km]") ; da.setInfoOnComponent(2,"Y [m]")
-        da.setInfoOnComponent(3,"Z [W]") ; da.setInfoOnComponent(4,"ZZ [km]") ; 
+        da.setInfoOnComponent(3,"Z [W]") ; da.setInfoOnComponent(4,"ZZ [km]") ;
         da[:,2]=3.
         self.assertEqual([7., 8., 3., 10., 11., 12., 13., 3., 15., 16., 17., 18., 3., 20., 21., 22., 23., 3., 25., 26.],da.getValues())
         da.rearrange(1) ; da.iota(7) ; da.rearrange(5)
@@ -1926,7 +1926,7 @@ class MEDCouplingBasicsTest3(unittest.TestCase):
         self.assertTrue(f1.isEqual(f2,1e-12,1e-12));
         #
         pass
-    
+
     def testGetDistributionOfTypes1(self):
         m=MEDCouplingDataForTest.build2DTargetMesh_1();
         tab1=[2,0,1,3,4]
@@ -2076,7 +2076,7 @@ class MEDCouplingBasicsTest3(unittest.TestCase):
         m3.setName(m.getName());
         self.assertTrue(m.isEqual(m3,1e-12));
         pass
-    
+
     def testChangeUnderlyingMeshWithCMesh1(self):
         mesh=MEDCouplingCMesh.New();
         coordsX=DataArrayDouble.New();
@@ -2154,7 +2154,7 @@ class MEDCouplingBasicsTest3(unittest.TestCase):
         self.assertEqual(0,c.getNbOfElems());
         self.assertEqual(1,cI.getNbOfElems());
         self.assertEqual([0],cI.getValues())
-        
+
         array12=[0.]*(6*5)
         da.setValues(array12,6,5) #bad NumberOfComponents
         self.assertRaises(InterpKernelException, da.findCommonTuples, 1e-2);
@@ -2390,7 +2390,7 @@ class MEDCouplingBasicsTest3(unittest.TestCase):
             self.assertAlmostEqual(expected5[i],m3.getCoords().getIJ(0,i),12);
             pass
         pass
-    
+
     def testBuildPartOfMySelfSafe1(self):
         mesh=MEDCouplingDataForTest.build2DTargetMesh_1()
         self.assertRaises(InterpKernelException,mesh.buildPartOfMySelf,[0,-1,4,2],True)
@@ -2509,68 +2509,12 @@ class MEDCouplingBasicsTest3(unittest.TestCase):
 
         myCoords = DataArrayDouble.New(mcoords, 3, 2)
         m1.setCoords(myCoords)
-        
+
         m2 = m1.deepCopy()
         m2.tessellate2D(0.1)
         # If the following raises, the test will fail automatically:
         m2.checkConsistency(0.0) # eps param not used
 
-    def testIntersect2DMeshesTmp4(self):
-        m1Coords=[0.,0.,1.,0.,1.5,0.,0.,1.,0.,1.5,-1.,0.,-1.5,0.,0.,-1,0.,-1.5,0.5,0.,1.25,0.,0.70710678118654757,0.70710678118654757,1.0606601717798214,1.0606601717798214,0.,0.5,0.,1.25,-0.70710678118654757,0.70710678118654757,-1.0606601717798214,1.0606601717798214,-0.5,0.,-1.25,0.,-0.70710678118654757,-0.70710678118654757,-1.0606601717798214,-1.0606601717798214,0.,-0.5,0.,-1.25,0.70710678118654757,-0.70710678118654757,1.0606601717798214,-1.0606601717798214];
-        m1Conn=[0,3,1,13,11,9, 3,4,2,1,14,12,10,11, 5,3,0,15,13,17, 6,4,3,5,16,14,15,18, 5,0,7,17,21,19, 6,5,7,8,18,19,22,20, 0,1,7,9,23,21, 1,2,8,7,10,24,22,23];
-        m1=MEDCouplingUMesh.New();
-        m1.setMeshDimension(2);
-        m1.allocateCells(8);
-        m1.insertNextCell(NORM_TRI6,6,m1Conn[0:6]);
-        m1.insertNextCell(NORM_QUAD8,8,m1Conn[6:14]);
-        m1.insertNextCell(NORM_TRI6,6,m1Conn[14:20]);
-        m1.insertNextCell(NORM_QUAD8,8,m1Conn[20:28]);
-        m1.insertNextCell(NORM_TRI6,6,m1Conn[28:34]);
-        m1.insertNextCell(NORM_QUAD8,8,m1Conn[34:42]);
-        m1.insertNextCell(NORM_TRI6,6,m1Conn[42:48]);
-        m1.insertNextCell(NORM_QUAD8,8,m1Conn[48:56]);
-        m1.finishInsertingCells();
-        myCoords1=DataArrayDouble.New();
-        myCoords1.setValues(m1Coords,25,2);
-        m1.setCoords(myCoords1);
-        #
-        m2Coords=[0.,0.,1.1,0.,1.1,1.,0.,1.,1.7,0.,1.7,1.,-1.1,1.,-1.1,0.,-1.7,0.,-1.7,1.,-1.7,-1,-1.1,-1.,0.,-1.,1.1,-1,1.7,-1.]
-        m2Conn=[0,3,2,1, 1,2,5,4, 7,6,3,0, 8,9,6,7, 7,0,12,11, 8,7,11,10, 0,1,13,12, 1,4,14,13]
-        m2=MEDCouplingUMesh.New();
-        m2.setMeshDimension(2);
-        m2.allocateCells(8);
-        for i in range(8):
-            m2.insertNextCell(NORM_QUAD4,4,m2Conn[4*i:4*(i+1)])
-            pass
-        m2.finishInsertingCells();
-        myCoords2=DataArrayDouble.New();
-        myCoords2.setValues(m2Coords,15,2);
-        m2.setCoords(myCoords2);
-        #
-        m3,d1,d2=MEDCouplingUMesh.Intersect2DMeshes(m2,m1,1e-10)
-        m3.unPolyze()
-        #
-        expected1=[0,0,1,1,2,2,3,3,4,4,5,5,6,6,7,7]
-        expected2=[0,1,1,-1,2,3,3,-1,4,5,5,-1,6,7,7,-1]
-        self.assertEqual(16,d1.getNumberOfTuples());
-        self.assertEqual(16,d2.getNumberOfTuples());
-        self.assertEqual(16,m3.getNumberOfCells());
-        self.assertEqual(104,m3.getNumberOfNodes());
-        self.assertEqual(2,m3.getSpaceDimension());
-        self.assertEqual(expected1,d1.getValues());
-        self.assertEqual(expected2,d2.getValues());
-        expected3=[6,16,15,18,44,45,46,8,18,2,1,16,47,48,49,50,8,17,1,2,40,51,52,53,54,8,40,5,4,17,55,56,57,58,6,18,15,20,59,60,61,8,20,7,6,18,62,63,64,65,8,41,6,7,21,66,67,68,69,8,21,8,9,41,70,71,72,73,6,20,15,22,74,75,76,8,22,11,7,20,77,78,79,80,8,21,7,11,42,81,82,83,84,8,42,10,8,21,85,86,87,88,6,22,15,16,89,90,91,8,16,1,13,22,92,93,94,95,8,43,13,1,17,96,97,98,99,8,17,4,14,43,100,101,102,103]
-        expected4=[0,7,16,25,34,41,50,59,68,75,84,93,102,109,118,127,136]
-        expected5=[0.,0.,1.1, 0.,1.1,1.,0.,1.,1.7,0.,1.7,1.,-1.1,1.,-1.1,0.,-1.7,0.,-1.7,1.,-1.7,-1.,-1.1,-1.,0.,-1.,1.1,-1.,1.7,-1.,0.,0.,1.,0.,1.5,0.,0.,1.,0.,1.5,-1.,0.,-1.5,0.,0.,-1.,0.,-1.5,0.5,0.,1.25,0.,0.7071067811865476,0.7071067811865476,1.0606601717798214,1.0606601717798214,0.,0.5,0.,1.25,-0.7071067811865476,0.7071067811865476,-1.0606601717798214,1.0606601717798214,-0.5,0.,-1.25,0.,-0.7071067811865476,-0.7071067811865476,-1.0606601717798214,-1.0606601717798214,0.,-0.5,0.,-1.25,0.7071067811865476,-0.7071067811865476,1.0606601717798214,-1.0606601717798214,1.1180339887498951,1.,-1.1180339887498951,1.,-1.1180339887498951,-1.,1.1180339887498951,-1.,0.5,0.,0.,0.5,0.7071067811865477,0.7071067811865476,0.55,1.,1.1,0.5,1.05,0.,0.7071067811865477,0.7071067811865475,1.3,0.,1.1,0.5,1.1090169943749475,1.,1.4012585384440737,0.535233134659635,1.4090169943749475,1.,1.7,0.5,1.6,0.,1.4012585384440737,0.535233134659635,0.,0.5,-0.5,0.,-0.7071067811865477,0.7071067811865476,-1.05,0.,-1.1,0.5,-0.55,1.,-0.7071067811865478,0.7071067811865475,-1.1090169943749475,1.,-1.1,0.5,-1.3,0.,-1.4012585384440737,0.5352331346596344,-1.6,0.,-1.7,0.5,-1.4090169943749475,1.,-1.4012585384440737,0.5352331346596344,-0.5,0.,0.,-0.5,-0.7071067811865475,-0.7071067811865477,-0.55,-1.,-1.1,-0.5,-1.05,0.,-0.7071067811865475,-0.7071067811865477,-1.3,0.,-1.1,-0.5,-1.1090169943749475,-1.,-1.4012585384440734,-0.5352331346596354,-1.4090169943749475,-1.,-1.7,-0.5,-1.6,0.,-1.4012585384440732,-0.5352331346596354,0.,-0.5,0.5,0.,0.7071067811865475,-0.7071067811865477,1.05,0.,1.1,-0.5,0.55,-1.,0.7071067811865475,-0.7071067811865477,1.1090169943749475,-1.,1.1,-0.5,1.3,0.,1.4012585384440737,-0.535233134659635,1.6,0.,1.7,-0.5,1.4090169943749475,-1.,1.4012585384440737,-0.535233134659635]
-        self.assertEqual(136,m3.getNodalConnectivity().getNumberOfTuples());
-        self.assertEqual(17,m3.getNodalConnectivityIndex().getNumberOfTuples());
-        self.assertEqual(expected3,m3.getNodalConnectivity().getValues());
-        self.assertEqual(expected4,m3.getNodalConnectivityIndex().getValues());
-        for i in range(208):
-            self.assertAlmostEqual(expected5[i],m3.getCoords().getIJ(0,i),12);
-            pass
-        pass
-
     def testGetCellIdsCrossingPlane1(self):
         mesh3D,mesh2D=MEDCouplingDataForTest.build3DExtrudedUMesh_1();
         vec=[-0.07,1.,0.07]
@@ -2994,7 +2938,7 @@ class MEDCouplingBasicsTest3(unittest.TestCase):
             self.assertEqual(expected3[i],m1c.getNodeIdsOfCell(i))
             pass
         pass
-          
+
     pass
 
 if __name__ == '__main__':

diff --git a/src/MEDCoupling_Swig/MEDCouplingIntersectTest.py b/src/MEDCoupling_Swig/MEDCouplingIntersectTest.py
index ca079bc596ff292525a061936e5cb8378209cf53..2bb2f88511b2a16757b2e6ed308a67ca097b9b1d 100644 (file)
--- a/src/MEDCoupling_Swig/MEDCouplingIntersectTest.py
+++ b/src/MEDCoupling_Swig/MEDCouplingIntersectTest.py
@@ -945,6 +945,34 @@ class MEDCouplingIntersectTest(unittest.TestCase):
         self.assertEqual(m1.getValues(), [0,0,1,1])
         self.assertEqual(m2.getValues(), [0,1,   -1,-1,  -1,-1,   2,3,  0,1])
 
+    def testSwig2Intersect2DMeshWith1DLine19(self):
+        """ Intersection arc of circle / segment was not properly detecting tangent cases """
+        eps=1.0e-5  # was working at 1.0e-8, but should also really work with 1.0e-5
+        mesh = MEDCouplingUMesh('layer_1', 2)
+        coo = DataArrayDouble([(55.4,3.7239),(61.4,7.188),(61.4,13.943),(49.55,7.1014),
+                                  (61.4,10.5655),(58.4,5.45595),(52.475,5.41265),(55.475,10.5222),
+                                  (56.9,9.34),(56.3343,7.97431),(56.9,7.74),(57.4657,7.97431),(59.4328,7.58116),
+                                  (55.8672,5.84911),    (0.,0.)])
+        mesh.setCoords(coo)
+        c = DataArrayInt([32, 0, 3, 2, 1, 11, 9,     6, 7, 4, 12, 8, 13])
+        cI = DataArrayInt([0, 13])
+        mesh.setConnectivity(c, cI)
+        tool = MEDCouplingUMesh('segment', 1)
+        coo = DataArrayDouble([(-166.611,-119.951),(269.611,131.902)])
+        tool.setCoords(coo)
+        c = DataArrayInt([1, 0, 1])
+        cI = DataArrayInt([0, 3])
+        tool.setConnectivity(c, cI)
+
+        res2D, res1D, m1, m2 = MEDCouplingUMesh.Intersect2DMeshWith1DLine(mesh, tool, eps)
+
+        self.assertEqual(res2D.getNodalConnectivity().getValues(),[32, 19, 17, 3, 2, 18, 20, 33, 34, 35, 36, 37, 38, 32, 1, 11, 20, 18, 39, 40, 41, 42, 32, 9, 0, 17, 19, 29, 30, 31, 32])
+        self.assertEqual(res2D.getNodalConnectivityIndex().getValues(),[0, 13, 22, 31])
+        self.assertEqual(res1D.getNodalConnectivity().getValues(),[1, 15, 17, 1, 17, 19, 1, 19, 20, 1, 20, 18, 1, 18, 16])
+        self.assertEqual(res1D.getNodalConnectivityIndex().getValues(),[0, 3, 6, 9, 12, 15])
+        self.assertEqual(m1.getValues(), [0, 0, 0])
+        self.assertEqual(m2.getValues(), [-1, -1, 0, 2, -1, -1, 0, 1, -1, -1])
+
     def testSwig2Conformize2D1(self):
         eps = 1.0e-8
         coo = [0.,-0.5,0.,0.,0.5,0.,0.5,-0.5,0.25,