[Intersect2D] Bug fix: residual cell construction

+ was buggy when mesh2 fully included in mesh1 and points
of mesh2 on edges of mesh1

diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.cxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.cxx
index 658521b05c6d7d0f1eac89a2d6acef68e4abc3fd..38dbea803b6b5dd542d0f8027819e12140ff472c 100644 (file)
--- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.cxx
+++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.cxx
@@ -1304,24 +1304,33 @@ void QuadraticPolygon::ComputeResidual(const QuadraticPolygon& pol1, const std::
       // retPolsUnderConstruction initially empty -> see if(!pol1Zip.empty()) below ...
       for(std::list<QuadraticPolygon *>::iterator itConstr=retPolsUnderContruction.begin();itConstr!=retPolsUnderContruction.end();)
         {
-          if((*itConstr)->getStartNode()==(*itConstr)->getEndNode())  // reconstruction of a cell is finished
-            {
-              itConstr++;
-              continue;
-            }
-          Node *curN=(*itConstr)->getEndNode();
-          bool smthHappened=false;
+          Node *startN = (*itConstr)->getStartNode();
+          Node *curN = (*itConstr)->getEndNode();
+          if(startN == curN)  // reconstruction of a cell is finished
+            {  itConstr++; continue;  }
+
+          bool smthHappened=false, doneEarly=false;
           // Complete a partially reconstructed polygon with boundary edges by matching nodes:
           for(std::list<Edge *>::iterator it2=edgesInPol2OnBoundaryL.begin();it2!=edgesInPol2OnBoundaryL.end();)
             {
-              if(curN==(*it2)->getStartNode())
-                { (*it2)->incrRef(); (*itConstr)->pushBack(new ElementaryEdge(*it2,true)); curN=(*it2)->getEndNode(); smthHappened=true; it2=edgesInPol2OnBoundaryL.erase(it2); }
-              else if(curN==(*it2)->getEndNode())
-                { (*it2)->incrRef(); (*itConstr)->pushBack(new ElementaryEdge(*it2,false)); curN=(*it2)->getStartNode(); smthHappened=true; it2=edgesInPol2OnBoundaryL.erase(it2); }
+              if(curN==(*it2)->getEndNode())  // only end node should be considered if orientation is correct for input meshes
+                                              // in the funny case of cells exactly included (see test case testIntersect2DMeshesTmp13) this is mandatory to take edges from pol2 in the right order.
+                {
+                  (*it2)->incrRef();
+                  (*itConstr)->pushBack(new ElementaryEdge(*it2,false));
+                  curN=(*it2)->getStartNode();
+                  smthHappened=true;
+                  it2=edgesInPol2OnBoundaryL.erase(it2);
+                }
               else
                 it2++;
+              if (curN == startN) // we might end here
+                { doneEarly = true;  break;  }
             }
-          if(smthHappened)
+          if (doneEarly)
+            {  itConstr++; continue;  }
+
+          if(smthHappened) // Now continue the construction by finding the next bit in pol1Zip. Not too sure what are the cases where the boolean if False ...
             {
               for(std::list<QuadraticPolygon *>::iterator itZip=pol1Zip.begin();itZip!=pol1Zip.end();)
                 {

diff --git a/src/MEDCoupling_Swig/MEDCouplingIntersectTest.py b/src/MEDCoupling_Swig/MEDCouplingIntersectTest.py
index 8ecb6e48d05bf12f005873e0d06f0da554d496ec..3887e6f723c1f5796f98481bad8b4539eecce6d5 100644 (file)
--- a/src/MEDCoupling_Swig/MEDCouplingIntersectTest.py
+++ b/src/MEDCoupling_Swig/MEDCouplingIntersectTest.py
@@ -697,6 +697,43 @@ class MEDCouplingIntersectTest(unittest.TestCase):
         self.assertEqual(res2Tool.getValues(), [0, -1, -1])
         pass
 
+    def testIntersect2DMeshesTmp13(self):
+        """ Bug fix: when part of mesh2 is fully included in mesh1 and some points of mesh2 are touching edges of mesh1, reconstruction of
+        residual cells was not properly done.
+        """
+        eps = 1.0e-6
+        # First case: only two points of mesh2 touching edges of mesh1
+        coo1 = DataArrayDouble([(0,0) , (1,1),  (1,0)])
+        mesh1 = MEDCouplingUMesh("mesh1", 2)
+        mesh1.setCoords(coo1)
+        c = DataArrayInt([NORM_TRI3, 0,1,2])
+        cI = DataArrayInt([0,4])
+        mesh1.setConnectivity(c, cI)
+
+        coo2 = DataArrayDouble([(0.5,0) , (0.5,0.5),  (0.7,0.25)])
+        mesh2 = MEDCouplingUMesh("mesh1", 2)
+        mesh2.setCoords(coo2)
+        c = DataArrayInt([NORM_TRI3, 0,1,2])
+        cI = DataArrayInt([0,4])
+        mesh2.setConnectivity(c, cI)
+
+        meshinter, res2Back, res2Tool = MEDCouplingUMesh.Intersect2DMeshes(mesh1,mesh2,eps)
+        self.assertEqual(meshinter.getNodalConnectivity().getValues(), [5, 3, 4, 5, 5, 4, 1, 2, 3, 5, 5, 3, 0, 4])
+        self.assertEqual(meshinter.getNodalConnectivityIndex().getValues(), [0, 4, 10, 14])
+        self.assertEqual(res2Back.getValues(), [0, 0, 0])
+        self.assertEqual(res2Tool.getValues(), [0, -1, -1])
+
+        # Second case: all three points of mesh2 touching edges of mesh1
+        coo2[2, 0] = 1.0
+        mesh2.setConnectivity(c, cI)
+
+        meshinter, res2Back, res2Tool = MEDCouplingUMesh.Intersect2DMeshes(mesh1,mesh2,eps)
+        self.assertEqual(meshinter.getNodalConnectivity().getValues(), [5, 3, 4, 5, 5, 4, 1, 5, 5, 5, 2, 3, 5, 3, 0, 4])
+        self.assertEqual(meshinter.getNodalConnectivityIndex().getValues(), [0, 4, 8, 12, 16])
+        self.assertEqual(res2Back.getValues(), [0,0,0,0])
+        self.assertEqual(res2Tool.getValues(), [0,-1,-1,-1])
+        pass
+
     def testSwig2Intersect2DMeshWith1DLine1(self):
         """A basic test with no colinearity between m1 and m2."""
         i=MEDCouplingIMesh("mesh",2,[5,5],[0.,0.],[1.,1.])
@@ -1238,12 +1275,12 @@ class MEDCouplingIntersectTest(unittest.TestCase):
         mesh.setConnectivity(c, cI)
 
         tool = MEDCouplingUMesh('tool', 1)
-        coo = DataArrayDouble([(0.0, 182.78400),  (182.78400, 182.78400)]) 
+        coo = DataArrayDouble([(0.0, 182.78400),  (182.78400, 182.78400)])
         tool.setCoords(coo)
         c = DataArrayInt([NORM_SEG2,0,1])
         cI = DataArrayInt([0, 3])
         tool.setConnectivity(c, cI)
-        
+
         res2D, res1D, resToSelf, mapLeftRight = MEDCouplingUMesh.Intersect2DMeshWith1DLine(mesh, tool, eps)
         self.assertEqual(res1D.getNodalConnectivity().getValues(), [1, 6, 1, 1, 1, 8, 1, 8, 9, 1, 9, 7])
         self.assertEqual(res1D.getNodalConnectivityIndex().getValues(),[0, 3, 6, 9, 12])
@@ -1253,7 +1290,7 @@ class MEDCouplingIntersectTest(unittest.TestCase):
         self.assertEqual(resToSelf.getValues(), [0, 1, 1])
         self.assertEqual(mapLeftRight.getValues(), [-1, -1, 1, 0, 1, 2, -1, -1])
 
-            
+
 
     def testSwig2Conformize2D1(self):
         eps = 1.0e-8