From 305e7a33835941b70e0c8037351c92476ed4a08f Mon Sep 17 00:00:00 2001 From: abn Date: Mon, 9 Apr 2018 16:16:16 +0200 Subject: [PATCH] Bug fix: Intersect2DMeshWith1DLine now correctly handling closed lines ... --- .../Geometric2D/InterpKernelGeo2DEdge.hxx | 8 ++-- .../InterpKernelGeo2DQuadraticPolygon.cxx | 5 ++- .../MEDCouplingUMesh_intersection.cxx | 42 ++++++++++++++++--- .../MEDCouplingIntersectTest.py | 27 ++++++++++++ 4 files changed, 70 insertions(+), 12 deletions(-) diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdge.hxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdge.hxx index 43eed2a87..c847530f2 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdge.hxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdge.hxx @@ -124,10 +124,10 @@ namespace INTERP_KERNEL Node *getNodeAndReleaseIt() { Node *tmp=_node; _node=0; return tmp; } ~IntersectElement(); private: - bool _1S; - bool _1E; - bool _2S; - bool _2E; + bool _1S; // true if starting point of edge 1 is located exactly on edge 2 (not nearby) + bool _1E; // true if ending point of edge 1 is located exactly on edge 2 (not nearby) + bool _2S; // true if starting point of edge 2 is located exactly on edge 1 (not nearby) + bool _2E; // true if ending point of edge 2 is located exactly on edge 1 (not nearby) double _chararct_val_for_e1; double _chararct_val_for_e2; Node *_node; diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.cxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.cxx index 3cb86c501..ef169e3d6 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.cxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.cxx @@ -293,6 +293,7 @@ void QuadraticPolygon::splitAbs(QuadraticPolygon& other, { double xBaryBB, yBaryBB; double fact=normalizeExt(&other, xBaryBB, yBaryBB); + // IteratorOnComposedEdge it1(this),it3(&other); MergePoints merge; @@ -300,13 +301,13 @@ void QuadraticPolygon::splitAbs(QuadraticPolygon& other, ComposedEdge *c2=new ComposedEdge; int i=0; std::map mapAddCoo; - for(it3.first();!it3.finished();it3.next(),i++)//iteration over 'other' _sub_edges + for(it3.first();!it3.finished();it3.next(),i++)//iteration over 'other->_sub_edges' { QuadraticPolygon otherTmp; ElementaryEdge* curE3=it3.current(); otherTmp.pushBack(new ElementaryEdge(curE3->getPtr(),curE3->getDirection())); curE3->getPtr()->incrRef(); IteratorOnComposedEdge it2(&otherTmp); - for(it2.first();!it2.finished();it2.next())//iteration on subedges of 'other->_sub_edge' + for(it2.first();!it2.finished();it2.next())//iteration on subedges of 'otherTmp->_sub_edge' { ElementaryEdge* curE2=it2.current(); if(!curE2->isThereStartPoint()) diff --git a/src/MEDCoupling/MEDCouplingUMesh_intersection.cxx b/src/MEDCoupling/MEDCouplingUMesh_intersection.cxx index 9aaae278b..6f2918148 100644 --- a/src/MEDCoupling/MEDCouplingUMesh_intersection.cxx +++ b/src/MEDCoupling/MEDCouplingUMesh_intersection.cxx @@ -423,7 +423,7 @@ bool IsColinearOfACellOf(const std::vector< std::vector >& intersectEdge1, * This method has 4 inputs : * - a mesh 'm1' with meshDim==1 and a SpaceDim==2 * - a mesh 'm2' with meshDim==1 and a SpaceDim==2 - * - subDiv of size 'm2->getNumberOfCells()' that lists for each seg cell in 'm' the splitting node ids randomly sorted. + * - subDiv of size 'm2->getNumberOfCells()' that lists for each seg cell in 'm2' the splitting node ids randomly sorted. * The aim of this method is to sort the splitting nodes, if any, and to put them in 'intersectEdge' output parameter based on edges of mesh 'm2' * Nodes end up lying consecutively on a cutted edge. * \param m1 is expected to be a mesh of meshDimension equal to 1 and spaceDim equal to 2. No check of that is performed by this method. @@ -998,9 +998,9 @@ CellInfo& VectorOfCellInfo::get(int pos) * * Algorithm : \a splitMesh1D is cut into contiguous parts. Each contiguous parts will build incrementally the output 2D cells. * - * \param [in] allEdges a list of pairs (beginNode, endNode). Linked with \a allEdgesPtr to get the equation of edge. + * \param [in] allEdges a list of pairs (beginNode, endNode). Represents all edges (already cut) in the single 2D cell being handled here. Linked with \a allEdgesPtr to get the equation of edge. */ -MEDCouplingUMesh *BuildMesh2DCutInternal(double eps, const MEDCouplingUMesh *splitMesh1D, const std::vector& allEdges, const std::vector< MCAuto >& allEdgesPtr, int offset, +MEDCouplingUMesh *BuildMesh2DCutInternal(double eps, MEDCouplingUMesh *splitMesh1D, const std::vector& allEdges, const std::vector< MCAuto >& allEdgesPtr, int offset, MCAuto& idsLeftRight) { int nbCellsInSplitMesh1D(splitMesh1D->getNumberOfCells()); @@ -1020,6 +1020,32 @@ MEDCouplingUMesh *BuildMesh2DCutInternal(double eps, const MEDCouplingUMesh *spl idsLeftRight=DataArrayInt::New(); idsLeftRight->alloc(nbCellsInSplitMesh1D*2); idsLeftRight->fillWithValue(-2); idsLeftRight->rearrange(2); int *idsLeftRightPtr(idsLeftRight->getPointer()); VectorOfCellInfo pool(edge1Bis,edge1BisPtr); + + // Compute contiguous parts of splitMesh1D. We can not make the full assumption that segments are consecutive in the connectivity + // (even if the user correctly called orderConsecutiveCells1D()). Indeed the tool might be a closed line whose junction point is in + // splitMesh1D. There can be only one such a point, and if this happens this is necessarily at the start + // of the connectivity. + MCAuto renumb(DataArrayInt::New()); + renumb->alloc(nbCellsInSplitMesh1D,1); + const int * renumbP(renumb->begin()); + + int i, first=cSplitPtr[1]; + // Follow 1D line backward as long as it is connected: + for (i=nbCellsInSplitMesh1D-1; cSplitPtr[ciSplitPtr[i]+2] == first; i--) + first=cSplitPtr[ciSplitPtr[i]+1]; + if (i < nbCellsInSplitMesh1D-1) + { + // Build circular permutation to shift consecutive edges together + renumb->iota(i+1); + renumb->applyModulus(nbCellsInSplitMesh1D); + splitMesh1D->renumberCells(renumbP, false); + cSplitPtr = splitMesh1D->getNodalConnectivity()->begin(); + ciSplitPtr = splitMesh1D->getNodalConnectivityIndex()->begin(); + } + else + renumb->iota(); + // + for(int iStart=0;iStart partOfSplitMesh1D(static_cast(splitMesh1D->buildPartOfMySelfSlice(iStart,iEnd,1,true))); int pos(pool.getPositionOf(eps,partOfSplitMesh1D)); // @@ -1052,11 +1078,15 @@ MEDCouplingUMesh *BuildMesh2DCutInternal(double eps, const MEDCouplingUMesh *spl iStart=iEnd; } for(int mm=0;mm >& intersectEdge1, int offset, MCAuto& idsLeftRight) { diff --git a/src/MEDCoupling_Swig/MEDCouplingIntersectTest.py b/src/MEDCoupling_Swig/MEDCouplingIntersectTest.py index 36a1e5573..7c2eae6f7 100644 --- a/src/MEDCoupling_Swig/MEDCouplingIntersectTest.py +++ b/src/MEDCoupling_Swig/MEDCouplingIntersectTest.py @@ -892,6 +892,33 @@ class MEDCouplingIntersectTest(unittest.TestCase): self.assertTrue([0,1,2], c.getValues()) self.assertEqual([2,1], d.getValues()) + def testSwig2Intersect2DMeshWith1DLine18(self): + """ Rare case: a *closed* line used as a tool, with the closing point inside a 2D cell ... """ + tool = MEDCouplingUMesh('circle', 1) + coo = DataArrayDouble([(39.35,0),(27.8247,27.8247),(2.40949e-15,39.35),(-27.8247,27.8247),(-39.35,4.81899e-15),(-27.8247,-27.8247),(-7.22848e-15,-39.35),(27.8247,-27.8247),(39.35,7.39805e-15)]) + tool.setCoords(coo) + c = DataArrayInt([2, 3, 5, 8, 2, 5, 3, 4]) + cI = DataArrayInt([0, 4, 8]) + tool.setConnectivity(c, cI) + + meh = MEDCouplingUMesh('meh', 2) + coo = DataArrayDouble([(-26.4275,36.6199),(-23.5868,31.6996),(-34.1861,41.0993),(-30.3383,25.0214),(-40.1861,30.707),(-35.2622,27.8642),(-37.1861,35.9032),(-30.3068,38.8596),(-25.0071,34.1598),(-26.9625,28.3605),(-25.7138,32.5128),(-27.354,36.4726),(-36.9138,32.5128),(-27.354,28.553),(-26.8908,36.5462),(-28.8461,26.7872)]) + meh.setCoords(coo) + c = DataArrayInt([32, 0, 1, 3, 13, 11, 8, 9, 15, 10, 14, + 32, 3, 4, 2, 0, 11, 13, 5, 6, 7, 14, 12, 15]) + cI = DataArrayInt([0, 11, 24]) + meh.setConnectivity(c, cI) + + res2D, res1D, m1, m2 = MEDCouplingUMesh.Intersect2DMeshWith1DLine(meh, tool, 1e-12) + self.assertEqual(4, res2D.getNumberOfCells()) + self.assertEqual(res2D.getNodalConnectivity().getValues(),[32, 13, 11, 0, 1, 25, 19, 26, 33, 34, 35, 36, 37, 38, 39, 32, 3, 26, 19, 25, 40, 41, 42, 43, + 32, 4, 2, 0, 11, 13, 26, 27, 44, 45, 46, 47, 48, 49, 50, 32, 3, 27, 26, 51, 52, 53]) + self.assertEqual(res2D.getNodalConnectivityIndex().getValues(),[0, 15, 24, 39, 46]) + self.assertEqual(res1D.getNodalConnectivity().getValues(),[2, 19, 25, 28, 2, 25, 21, 29, 2, 21, 27, 30, 2, 27, 26, 31, 2, 26, 19, 32]) + self.assertEqual(res1D.getNodalConnectivityIndex().getValues(),[0, 4, 8, 12, 16, 20]) + self.assertEqual(m1.getValues(), [0,0,1,1]) + self.assertEqual(m2.getValues(), [0,1, -1,-1, -1,-1, 2,3, 0,1]) + def testSwig2Conformize2D1(self): eps = 1.0e-8 coo = [0.,-0.5,0.,0.,0.5,0.,0.5,-0.5,0.25, -- 2.39.2