From 2da6c2bc1ec47f323da8a053ceb76ba1de0b636b Mon Sep 17 00:00:00 2001 From: abn Date: Wed, 9 Jan 2019 15:30:10 +0100 Subject: [PATCH] Intersec bug fix: correct polygons with flat corners, they crash residual computation ... also added test. --- .../InterpKernelGeo2DComposedEdge.cxx | 12 +++++++ .../InterpKernelGeo2DComposedEdge.hxx | 1 + .../InterpKernelGeo2DElementaryEdge.cxx | 8 +++++ .../InterpKernelGeo2DElementaryEdge.hxx | 1 + .../InterpKernelGeo2DQuadraticPolygon.cxx | 33 +++++++++++++++++-- .../InterpKernelGeo2DQuadraticPolygon.hxx | 1 + .../MEDCouplingUMesh_intersection.cxx | 5 +++ .../MEDCouplingIntersectTest.py | 31 +++++++++++++++++ 8 files changed, 89 insertions(+), 3 deletions(-) diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DComposedEdge.cxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DComposedEdge.cxx index 3d0f0a468..82c8f7b9c 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DComposedEdge.cxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DComposedEdge.cxx @@ -97,6 +97,18 @@ void ComposedEdge::pushBack(ComposedEdge *elem) _sub_edges.insert(_sub_edges.end(),elemsOfElem->begin(),elemsOfElem->end()); } +/*! + * Warning! This is highly inefficient ... + */ +void ComposedEdge::erase(int index) +{ + // Not the most efficient thing to do, but rarely called ... + auto it = _sub_edges.begin(); + for (int i=0; i < index; i++, it++); + delete (*it); + _sub_edges.erase(it++, it); +} + ElementaryEdge *ComposedEdge::operator[](int i) const { std::list::const_iterator iter=_sub_edges.begin(); diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DComposedEdge.hxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DComposedEdge.hxx index 4d04aa7b9..c484b0020 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DComposedEdge.hxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DComposedEdge.hxx @@ -98,6 +98,7 @@ namespace INTERP_KERNEL INTERPKERNEL_EXPORT void pushBack(Edge *edge, bool direction=true); INTERPKERNEL_EXPORT void pushBack(ElementaryEdge *elem); INTERPKERNEL_EXPORT void pushBack(ComposedEdge *elem); + INTERPKERNEL_EXPORT void erase(int index); INTERPKERNEL_EXPORT int size() const { return (int)_sub_edges.size(); } INTERPKERNEL_EXPORT ElementaryEdge *operator[](int i) const; INTERPKERNEL_EXPORT Node *getEndNode() const; diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DElementaryEdge.cxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DElementaryEdge.cxx index 0f62a67c0..3f99fcde6 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DElementaryEdge.cxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DElementaryEdge.cxx @@ -72,6 +72,14 @@ void ElementaryEdge::getAllNodes(std::set& output) const output.insert(_ptr->getEndNode()); } +bool ElementaryEdge::hasSameExtremities(const ElementaryEdge& other) const +{ + std::set s1, s2; + getAllNodes(s1); + other.getAllNodes(s2); + return (s1 == s2); +} + void ElementaryEdge::getBarycenter(double *bary, double& weigh) const { _ptr->getBarycenter(bary); diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DElementaryEdge.hxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DElementaryEdge.hxx index 18f6efafd..f585b9aa1 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DElementaryEdge.hxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DElementaryEdge.hxx @@ -50,6 +50,7 @@ namespace INTERP_KERNEL INTERPKERNEL_EXPORT void applySimilarity(double xBary, double yBary, double dimChar) { _ptr->applySimilarity(xBary,yBary,dimChar); } INTERPKERNEL_EXPORT void unApplySimilarity(double xBary, double yBary, double dimChar) { _ptr->unApplySimilarity(xBary,yBary,dimChar); } INTERPKERNEL_EXPORT void getAllNodes(std::set& output) const; + INTERPKERNEL_EXPORT bool hasSameExtremities(const ElementaryEdge& other) const; INTERPKERNEL_EXPORT void getBarycenter(double *bary, double& weigh) const; INTERPKERNEL_EXPORT ElementaryEdge *clone() const; INTERPKERNEL_EXPORT void initLocations() const; diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.cxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.cxx index 26c7a7223..b4353458d 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.cxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.cxx @@ -463,7 +463,7 @@ void QuadraticPolygon::appendSubEdgeFromCrudeDataArray(Edge *baseEdge, std::size * This method builds from descending conn of a quadratic polygon stored in crude mode (MEDCoupling). Descending conn is in FORTRAN relative mode in order to give the * orientation of edge. */ -void QuadraticPolygon::buildFromCrudeDataArray2(const std::map& mapp, bool isQuad, const int *nodalBg, const double *coords, const int *descBg, const int *descEnd, const std::vector >& intersectEdges, +void QuadraticPolygon::buildFromCrudeDataArray2(const std::map& mapp, bool isQuad, const int *nodalBg, const double *coords, const int *descBg, const int *descEnd, const std::vector >& intersectEdges2, const INTERP_KERNEL::QuadraticPolygon& pol1, const int *descBg1, const int *descEnd1, const std::vector >& intersectEdges1, const std::vector< std::vector >& colinear1, std::map >& alreadyExistingIn2) @@ -518,7 +518,7 @@ void QuadraticPolygon::buildFromCrudeDataArray2(const std::map& subEdge=intersectEdges[edgeId]; + const std::vector& subEdge=intersectEdges2[edgeId]; std::size_t nbOfSubEdges=subEdge.size()/2; for(std::size_t j=0;j 2) + for(it2ii.first();!it2ii.finished();it2ii.next()) + { + ElementaryEdge * cur = it2ii.current(); + if (prevEdge && prevEdge->hasSameExtremities(*cur)) + { + // Delete the two 'identical' edges: + it2ii.previousLoop(); it2ii.previousLoop(); + erase(--kk); erase(kk); + prevEdge = it2ii.current(); + } + else + { + kk++; + prevEdge = cur; + } + } +} + /*! * Warning This method is \b NOT const. 'this' and 'other' are modified after call of this method. * 'other' is a QuadraticPolygon of \b non closed edges. diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.hxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.hxx index a5afe97e5..ee286fd39 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.hxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.hxx @@ -93,6 +93,7 @@ namespace INTERP_KERNEL INTERPKERNEL_EXPORT bool haveIAChanceToBeCompletedBy(const QuadraticPolygon& pol1NotSplitted, const QuadraticPolygon& pol2Splitted, bool& direction) const; INTERPKERNEL_EXPORT static void ComputeResidual(const QuadraticPolygon& pol1, const std::set& notUsedInPol1, const std::set& edgesInPol2OnBoundary, const std::map& mapp, int offset, int idThis, std::vector& addCoordsQuadratic, std::vector& conn, std::vector& connI, std::vector& nb1, std::vector& nb2); + INTERPKERNEL_EXPORT void cleanDegeneratedConsecutiveEdges(); protected: std::list zipConsecutiveInSegments() const; void dumpInXfigFile(std::ostream& stream, int resolution, const Bounds& box) const; diff --git a/src/MEDCoupling/MEDCouplingUMesh_intersection.cxx b/src/MEDCoupling/MEDCouplingUMesh_intersection.cxx index 39e7cb4ee..f61742e0b 100644 --- a/src/MEDCoupling/MEDCouplingUMesh_intersection.cxx +++ b/src/MEDCoupling/MEDCouplingUMesh_intersection.cxx @@ -1360,6 +1360,11 @@ void MEDCouplingUMesh::BuildIntersecting2DCellsFromEdges(double eps, const MEDCo pol2s[ii].buildFromCrudeDataArray2(mappRev,cm2.isQuadratic(),conn2+connI2[*it2]+1,coo2,desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2, pol1,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,colinear2, /* output */ edgesIn2ForShare); } + // The cleaning below must be done after the full construction of all pol2s to correctly deal with shared edges: + for (auto &p: pol2s) + p.cleanDegeneratedConsecutiveEdges(); + edgesIn2ForShare.clear(); // removing temptation to use it further since it might now contain invalid edges. + /// ii=0; // Now rebuild intersected cells from all this: for(std::vector::const_iterator it2=candidates2.begin();it2!=candidates2.end();it2++,ii++) diff --git a/src/MEDCoupling_Swig/MEDCouplingIntersectTest.py b/src/MEDCoupling_Swig/MEDCouplingIntersectTest.py index 7694d7870..f19780274 100644 --- a/src/MEDCoupling_Swig/MEDCouplingIntersectTest.py +++ b/src/MEDCoupling_Swig/MEDCouplingIntersectTest.py @@ -548,6 +548,37 @@ class MEDCouplingIntersectTest(unittest.TestCase): self.assertEqual(map2.getValues(), [0]) pass + def testIntersect2DMeshesTmp9(self): + """ Tricky case: two triangular shapes intersecting, but not perfectly, at their tips. Several issues fixed: + - Bug fix: seg seg intersector epsilon is to be taken absolutely for colinearity test (even for very small vectors + we don't want to have false positive on colinearity. So go back to a comparison with an angle.) + - when intersecting nodes are merged, they were not properly added on pol2. + - bug fix in compute residual: the stop condition is really on pol1Zip only. + - correcting polygons with flat corners, they were crashing residual computation + """ + eps = 1.0e-6 # This is the key parameter. DO NOT CHANGE IT. + back = MEDCouplingUMesh('crh8_rse3', 2) + coo = DataArrayDouble([(-31.313754538446631,-32.512813836330515),(-31.531462871779969,-32.135731941766032),(-31.422608705113298,-32.324272889048274),(-31.690836433011114,-32.295105502997181),(-31.621640616088342,-32.204927758688783),(-31.502295485728872,-32.403959669663848)]) + back.setCoords(coo) + c = DataArrayInt([32, 0, 3, 1, 5, 4, 2]) + cI = DataArrayInt([0, 7]) + back.setConnectivity(c, cI) + + tool = MEDCouplingUMesh('TA-536193G_expl_20181022_merged', 2) + coo = DataArrayDouble([(-29.918137808525149,-26.883223901634544),(-32.919909136264039,-26.939612990540404),(-27.866900000000001,-28.016680435212603),(-31.313800000000001,-32.512799999999999),(-27.866900000000001,-28.933918793630923)]) + tool.setCoords(coo) + c = DataArrayInt([5, 1, 0, 3, 5, 0, 2, 3, 5, 4, 3, 2]) + cI = DataArrayInt([0, 4, 8, 12]) + tool.setConnectivity(c, cI) + + inter, res2Back, res2Tool = MEDCouplingUMesh.Intersect2DMeshes(back, tool, eps) + + self.assertEqual(inter.getNodalConnectivity().getValues(), [5, 14, 13, 11, 12, 5, 13, 15, 11, 32, 12, 3, 1, 14, 16, 17, 18, 19, 5, 15, 0, 11]) + self.assertEqual(inter.getNodalConnectivityIndex().getValues(), [0, 5, 9, 18, 22]) + 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.]) -- 2.39.2