From bf608d33613f52a5209c15c286ab8a7238cbf8d7 Mon Sep 17 00:00:00 2001 From: abn Date: Thu, 29 Apr 2021 09:26:21 +0200 Subject: [PATCH] [Intersect2D] Bug fix: residual cell construction + was buggy when mesh2 fully included in mesh1 and points of mesh2 on edges of mesh1 --- .../InterpKernelGeo2DQuadraticPolygon.cxx | 33 ++++++++------ .../MEDCouplingIntersectTest.py | 43 +++++++++++++++++-- 2 files changed, 61 insertions(+), 15 deletions(-) diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.cxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.cxx index 658521b05..38dbea803 100644 --- 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::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::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::iterator itZip=pol1Zip.begin();itZip!=pol1Zip.end();) { diff --git a/src/MEDCoupling_Swig/MEDCouplingIntersectTest.py b/src/MEDCoupling_Swig/MEDCouplingIntersectTest.py index 8ecb6e48d..3887e6f72 100644 --- 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 -- 2.39.2