From 90b7c9c73054965b8085ab0ecd2646a3d59f475a Mon Sep 17 00:00:00 2001 From: Anida KHIZAR Date: Mon, 19 Oct 2020 16:24:22 +0200 Subject: [PATCH] Intersect2DMeshWith1DLine: bug fix (collinear edges not always detected) + with test case + and doc --- .../MEDCouplingUMesh_intersection.cxx | 67 +++++++++++++++++-- .../MEDCouplingIntersectTest.py | 28 ++++++++ 2 files changed, 89 insertions(+), 6 deletions(-) diff --git a/src/MEDCoupling/MEDCouplingUMesh_intersection.cxx b/src/MEDCoupling/MEDCouplingUMesh_intersection.cxx index 3fce3e6f9..b480bec30 100644 --- a/src/MEDCoupling/MEDCouplingUMesh_intersection.cxx +++ b/src/MEDCoupling/MEDCouplingUMesh_intersection.cxx @@ -556,7 +556,15 @@ void MEDCouplingUMesh::BuildIntersectEdges(const MEDCouplingUMesh *m1, const MED } } -MEDCouplingUMesh *BuildMesh1DCutFrom(const MEDCouplingUMesh *mesh1D, const std::vector< std::vector >& intersectEdge2, const DataArrayDouble *coords1, const std::vector& addCoo, const std::map& mergedNodes, const std::vector< std::vector >& colinear2, const std::vector< std::vector >& intersectEdge1, +/* + * Build the final 1D mesh resulting from the newly created points after intersection with the segments of the descending 2D mesh. + * @param[out] idsInRetColinear IDs of edges in the result (ret) that are colinears to one of the segment of the descending 2D mesh. Indexing scheme + * is the one of the ret 1D mesh. + * @param[out] idsInMesh1DForIdsInRetColinear same IDs as above in the descending 2D mesh + */ +MEDCouplingUMesh *BuildMesh1DCutFrom(const MEDCouplingUMesh *mesh1D, const std::vector< std::vector >& intersectEdge2, + const DataArrayDouble *coords1, const std::vector& addCoo, const std::map& mergedNodes, + const std::vector< std::vector >& colinear2, const std::vector< std::vector >& intersectEdge1, MCAuto& idsInRetColinear, MCAuto& idsInMesh1DForIdsInRetColinear) { idsInRetColinear=DataArrayIdType::New(); idsInRetColinear->alloc(0,1); @@ -581,7 +589,8 @@ MEDCouplingUMesh *BuildMesh1DCutFrom(const MEDCouplingUMesh *mesh1D, const std:: mcIdType nbSubEdge=ToIdType(subEdges.size()/2); for(mcIdType j=0;j n1(MEDCouplingUMeshBuildQPNode(subEdges[2*j],coords1->begin(),offset1,coo2Ptr,offset2,addCoo)),n2(MEDCouplingUMeshBuildQPNode(subEdges[2*j+1],coords1->begin(),offset1,coo2Ptr,offset2,addCoo)); + MCAuto n1(MEDCouplingUMeshBuildQPNode(subEdges[2*j],coords1->begin(),offset1,coo2Ptr,offset2,addCoo)), + n2(MEDCouplingUMeshBuildQPNode(subEdges[2*j+1],coords1->begin(),offset1,coo2Ptr,offset2,addCoo)); MCAuto e2(e->buildEdgeLyingOnMe(n1,n2)); INTERP_KERNEL::Edge *e2Ptr(e2); std::map::const_iterator itm; @@ -619,6 +628,7 @@ MEDCouplingUMesh *BuildMesh1DCutFrom(const MEDCouplingUMesh *mesh1D, const std:: } e->decrRef(); } + MCAuto ret(MEDCouplingUMesh::New(mesh1D->getName(),1)); ret->setConnectivity(cOut,ciOut,true); MCAuto arr3(DataArrayDouble::New()); @@ -1821,13 +1831,56 @@ void MEDCouplingUMesh::Intersect2DMeshWith1DLine(const MEDCouplingUMesh *mesh2D, MCAuto idsInRet1Colinear,idsInDescMesh2DForIdsInRetColinear; MCAuto ret2(DataArrayIdType::New()); ret2->alloc(0,1); MCAuto ret1(BuildMesh1DCutFrom(mesh1D,intersectEdge2,mesh2D->getCoords(),addCoo,mergedNodes,colinear2,intersectEdge1, - idsInRet1Colinear,idsInDescMesh2DForIdsInRetColinear)); + idsInRet1Colinear,idsInDescMesh2DForIdsInRetColinear)); + + // ### Colinearity fix : + // if a node in ret1 has been merged with an already existing node in mesh2D, + // we might end up with edges in ret1 that are colinear with some edges in mesh2D + // even if none of the edges of the two original meshes were colinear. + // this procedure detects such edges and adds them in idsInRet1Colinear/idsInDescMesh2DForIdsInRetColinear + // a- find edges in ret1 that are in contact with 2 cells + MCAuto centerOfMassRet1(ret1->computeCellCenterOfMass()); + MCAuto cells,cellsIndex; + mesh2D->getCellsContainingPoints(centerOfMassRet1->begin(),centerOfMassRet1->getNumberOfTuples(),eps,cells,cellsIndex); + MCAuto cellsIndex2(DataArrayIdType::New()); cellsIndex2->alloc(0,1); + if (cellsIndex->getNumberOfTuples() > 1) + cellsIndex2 = cellsIndex->deltaShiftIndex(); + MCAuto idsInRet1With2Contacts(cellsIndex2->findIdsEqual(2)); + + MCAuto realIdsInDesc2D(desc1->deepCopy()); + realIdsInDesc2D->abs(); realIdsInDesc2D->applyLin(1,-1); + const mcIdType *cRet1(ret1->getNodalConnectivity()->begin()),*ciRet1(ret1->getNodalConnectivityIndex()->begin()); + for(const mcIdType *it=idsInRet1With2Contacts->begin();it!=idsInRet1With2Contacts->end();it++) + { + // b- find the edge that the 2 cells in m1Desc have in common: + // this is the edge which is colinear with the one in ret1 + const mcIdType* cellId1 = cells->begin() + cellsIndex->begin()[*it]; + const mcIdType* cellId2 = cells->begin() + cellsIndex->begin()[*it]+1; + + std::set s1(realIdsInDesc2D->begin()+dd2->begin()[*cellId1], realIdsInDesc2D->begin()+dd2->begin()[*cellId1+1]); + std::set s2(realIdsInDesc2D->begin()+dd2->begin()[*cellId2], realIdsInDesc2D->begin()+dd2->begin()[*cellId2+1]); + mcIdType commonEdgeId; + std::set_intersection(s1.begin(),s1.end(),s2.begin(),s2.end(), &commonEdgeId); + + // c- find correct orientation for commonEdgeId + const mcIdType* firstNodeInColinearEdgeRet1 = cRet1 + ciRet1[*it]+1; + const mcIdType* secondNodeInColinearEdgeRet1 = cRet1 + ciRet1[*it]+2; + std::vector v; v.push_back(commonEdgeId); + if(IsColinearOfACellOf(intersectEdge1, v, *firstNodeInColinearEdgeRet1,*secondNodeInColinearEdgeRet1,commonEdgeId)) + { + idsInRet1Colinear->pushBackSilent(*it); + idsInDescMesh2DForIdsInRetColinear->pushBackSilent(commonEdgeId); + } + } + // ### End colinearity fix + MCAuto ret3(DataArrayIdType::New()); ret3->alloc(ret1->getNumberOfCells()*2,1); ret3->fillWithValue(std::numeric_limits::max()); ret3->rearrange(2); MCAuto idsInRet1NotColinear(idsInRet1Colinear->buildComplement(ret1->getNumberOfCells())); // deal with cells in mesh2D that are not cut but only some of their edges are MCAuto idsInDesc2DToBeRefined(idsInDescMesh2DForIdsInRetColinear->deepCopy()); idsInDesc2DToBeRefined->abs(); idsInDesc2DToBeRefined->applyLin(1,-1); idsInDesc2DToBeRefined=idsInDesc2DToBeRefined->buildUnique(); + MCAuto out0s;//ids in mesh2D that are impacted by the fact that some edges of \a mesh1D are part of the edges of those cells if(!idsInDesc2DToBeRefined->empty()) { @@ -1840,9 +1893,11 @@ void MEDCouplingUMesh::Intersect2DMeshWith1DLine(const MEDCouplingUMesh *mesh2D, } // MCAuto ret1NonCol(static_cast(ret1->buildPartOfMySelf(idsInRet1NotColinear->begin(),idsInRet1NotColinear->end()))); - MCAuto baryRet1(ret1NonCol->computeCellCenterOfMass()); - MCAuto elts,eltsIndex; - mesh2D->getCellsContainingPoints(baryRet1->begin(),baryRet1->getNumberOfTuples(),eps,elts,eltsIndex); + MCAuto baryRet1(centerOfMassRet1->selectByTupleId(idsInRet1NotColinear->begin(), idsInRet1NotColinear->end())); + DataArrayIdType *out(0),*outi(0); + DataArrayIdType::ExtractFromIndexedArrays(idsInRet1NotColinear->begin(),idsInRet1NotColinear->end(),cells,cellsIndex,out,outi); + MCAuto elts(out),eltsIndex(outi); + MCAuto eltsIndex2(DataArrayIdType::New()); eltsIndex2->alloc(0,1); if (eltsIndex->getNumberOfTuples() > 1) eltsIndex2 = eltsIndex->deltaShiftIndex(); diff --git a/src/MEDCoupling_Swig/MEDCouplingIntersectTest.py b/src/MEDCoupling_Swig/MEDCouplingIntersectTest.py index 2033bd948..6064d93de 100644 --- a/src/MEDCoupling_Swig/MEDCouplingIntersectTest.py +++ b/src/MEDCoupling_Swig/MEDCouplingIntersectTest.py @@ -1227,6 +1227,34 @@ class MEDCouplingIntersectTest(unittest.TestCase): self.assertEqual(mapLeftRight.getValues(), [-1, -1, 0, 4, -1, -1, 1, 4, -1, -1, 1, 3, -1, -1, 2, 3, -1, -1]) pass + def testSwig2Intersect2DMeshWith1DLine21(self): + """ A line intersecting a cell very close to one of its node (collinearity not detected) """ + eps=1.0e-5 # was working at 1.0e-8, but should also really work with 1.0e-5 + mesh = MEDCouplingUMesh('mesh', 2) + coo = DataArrayDouble([(110.65324,180.56968),(112.01128,182.78580),(113.36932,185.00192),(118.27200,181.90669),(118.27200,178.79852),(118.27200,175.67380)]) + mesh.setCoords(coo) + c = DataArrayInt([NORM_QUAD4, 0, 1, 4, 5, NORM_QUAD4, 1, 2, 3, 4]) + cI = DataArrayInt([0, 5, 10]) + mesh.setConnectivity(c, cI) + + tool = MEDCouplingUMesh('tool', 1) + 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]) + self.assertEqual(res2D.getNodalConnectivity().getValues(), [5, 0, 1, 8, 4, 5, 5, 1, 2, 9, 8, 5, 3, 4, 8, 9]) + self.assertEqual(res2D.getNodalConnectivityIndex().getValues(),[0, 6, 11, 16]) + + 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 coo = [0.,-0.5,0.,0.,0.5,0.,0.5,-0.5,0.25, -- 2.39.2