From 922c4272c03cf8618ffa3964f0411c10922cf12f Mon Sep 17 00:00:00 2001 From: abn Date: Fri, 15 Feb 2019 14:08:51 +0100 Subject: [PATCH] Intersector: Earlier detection of merged nodes. Use a BBTreePoints to detect merged nodes. Allows to more efficiently use identifyEarlyIntersection() and avoid degenerated situations again. --- .../InterpKernelGeo2DEdgeArcCircle.hxx | 2 +- .../Geometric2D/InterpKernelGeo2DNode.cxx | 12 ++-- .../MEDCouplingUMesh_intersection.cxx | 59 ++++++++++++++++++- 3 files changed, 65 insertions(+), 8 deletions(-) diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.hxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.hxx index 481fb8c6f..ea2373dc0 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.hxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.hxx @@ -44,7 +44,7 @@ namespace INTERP_KERNEL const EdgeArcCircle& getE1() const { return (const EdgeArcCircle&)_e1; } const EdgeArcCircle& getE2() const { return (const EdgeArcCircle&)_e2; } private: - double _dist; + double _dist; // distance between the two arc centers }; /** diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DNode.cxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DNode.cxx index 4d0ce0802..c48d783f9 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DNode.cxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DNode.cxx @@ -151,16 +151,16 @@ void Node::unApplySimilarity(double xBary, double yBary, double dimChar) void Node::fillGlobalInfoAbs(const std::map& mapThis, const std::map& mapOther, int offset1, int offset2, double fact, double baryX, double baryY, std::vector& addCoo, std::map& mapAddCoo, int *nodeId) const { - std::map::const_iterator it=mapThis.find(const_cast(this)); - if(it!=mapThis.end()) + std::map::const_iterator it=mapOther.find(const_cast(this)); + if(it!=mapOther.end()) // order matters, try in mapOther first. { - *nodeId=(*it).second; + *nodeId=(*it).second+offset1; return; } - it=mapOther.find(const_cast(this)); - if(it!=mapOther.end()) + it=mapThis.find(const_cast(this)); + if(it!=mapThis.end()) { - *nodeId=(*it).second+offset1; + *nodeId=(*it).second; return; } it=mapAddCoo.find(const_cast(this)); diff --git a/src/MEDCoupling/MEDCouplingUMesh_intersection.cxx b/src/MEDCoupling/MEDCouplingUMesh_intersection.cxx index d6986299b..aa20071e1 100644 --- a/src/MEDCoupling/MEDCouplingUMesh_intersection.cxx +++ b/src/MEDCoupling/MEDCouplingUMesh_intersection.cxx @@ -271,6 +271,57 @@ namespace MEDCoupling return ret; } + INTERP_KERNEL::QuadraticPolygon *MEDCouplingUMeshBuildQPFromMeshWithTree(const MEDCouplingUMesh *mDesc, const std::vector& candidates, + std::map& mapp, + const BBTreePts<2,int> & nodeTree, + const std::map& mapRev) + { + mapp.clear(); + std::map mapp2; // the last var is a flag specifying if node is an extreme node of the seg (LINEAR) or only a middle for SEG3 (QUADRATIC_ONLY). + const double *coo=mDesc->getCoords()->getConstPointer(); + const int *c=mDesc->getNodalConnectivity()->getConstPointer(); + const int *cI=mDesc->getNodalConnectivityIndex()->getConstPointer(); + std::set s; + for(std::vector::const_iterator it=candidates.begin();it!=candidates.end();it++) + s.insert(c+cI[*it]+1,c+cI[(*it)+1]); + for(std::set::const_iterator it2=s.begin();it2!=s.end();it2++) + { + INTERP_KERNEL::Node *n; + // Look for a potential node to merge + std::vector candNode; + nodeTree.getElementsAroundPoint(coo+2*(*it2), candNode); + if (candNode.size() > 2) + throw INTERP_KERNEL::Exception("MEDCouplingUMesh::MEDCouplingUMeshBuildQPFromMeshWithTree(): some nodes are not properly merged (within eps) in input mesh!"); + bool node_created=false; + if (candNode.size()) + { + auto itt=mapRev.find(candNode[0]); + if (itt != mapRev.end()) // we might hit a node which is in the coords array but not used in the connectivity in which case it won't be in the revMap + { + node_created=true; + n = (*itt).second; + n->incrRef(); + } + } + if(!node_created) + n = new INTERP_KERNEL::Node(coo[2*(*it2)],coo[2*(*it2)+1]); + mapp2[*it2]=INTERP_KERNEL::NodeWithUsage(n,INTERP_KERNEL::USAGE_UNKNOWN); + } + INTERP_KERNEL::QuadraticPolygon *ret=new INTERP_KERNEL::QuadraticPolygon; + for(std::vector::const_iterator it=candidates.begin();it!=candidates.end();it++) + { + INTERP_KERNEL::NormalizedCellType typ=(INTERP_KERNEL::NormalizedCellType)c[cI[*it]]; + ret->pushBack(MEDCouplingUMeshBuildQPFromEdge(typ,mapp2,c+cI[*it]+1)); // this call will set quad points to false in the map + } + for(std::map::const_iterator it2=mapp2.begin();it2!=mapp2.end();it2++) + { + if((*it2).second.second == INTERP_KERNEL::USAGE_LINEAR) + mapp[(*it2).second.first]=(*it2).first; + ((*it2).second.first)->decrRef(); + } + return ret; + } + INTERP_KERNEL::Node *MEDCouplingUMeshBuildQPNode(int nodeId, const double *coo1, int offset1, const double *coo2, int offset2, const std::vector& addCoo) { if(nodeId>=offset2) @@ -1235,6 +1286,7 @@ void MEDCouplingUMesh::Intersect1DMeshes(const MEDCouplingUMesh *m1Desc, const M colinear2.resize(nDescCell2); subDiv2.resize(nDescCell2); BBTree myTree(bbox2,0,0,m2Desc->getNumberOfCells(),-eps); + BBTreePts treeNodes2(m2Desc->getCoords()->begin(),0,0,m2Desc->getCoords()->getNumberOfTuples(),eps); std::vector candidates1(1); int offset1(m1Desc->getNumberOfNodes()); @@ -1246,10 +1298,15 @@ void MEDCouplingUMesh::Intersect1DMeshes(const MEDCouplingUMesh *m1Desc, const M if(!candidates2.empty()) // candidates2 holds edges from the second mesh potentially intersecting current edge i in mesh1 { std::map map1,map2; + std::map revMap2; // pol2 is not necessarily a closed polygon: just a set of (quadratic) edges (same as candidates2) in the Geometric DS format INTERP_KERNEL::QuadraticPolygon *pol2=MEDCouplingUMeshBuildQPFromMesh(m2Desc,candidates2,map2); + // Build revMap2 + for (auto& kv : map2) + revMap2[kv.second] = kv.first; candidates1[0]=i; - INTERP_KERNEL::QuadraticPolygon *pol1=MEDCouplingUMeshBuildQPFromMesh(m1Desc,candidates1,map1); + // In the construction of pol1 we might reuse nodes from pol2, that we have identified as to be merged. + INTERP_KERNEL::QuadraticPolygon *pol1=MEDCouplingUMeshBuildQPFromMeshWithTree(m1Desc,candidates1,map1,treeNodes2, revMap2); // This following part is to avoid that some removed nodes (for example due to a merge between pol1 and pol2) are replaced by a newly created one // This trick guarantees that Node * are discriminant (i.e. form a unique identifier) std::set nodes; -- 2.39.2