From ed26f73bb07b412b350ed1e42ec92bbb4f08a0e0 Mon Sep 17 00:00:00 2001 From: abn Date: Wed, 19 Dec 2018 12:23:44 +0100 Subject: [PATCH] Intersection: renaming some variables and refactor to make the algo easier to read. --- .../InterpKernelGeo2DComposedEdge.cxx | 23 ++ .../InterpKernelGeo2DComposedEdge.hxx | 1 + .../Geometric2D/InterpKernelGeo2DEdge.cxx | 25 +++ .../Geometric2D/InterpKernelGeo2DEdge.hxx | 9 +- .../InterpKernelGeo2DEdgeArcCircle.cxx | 25 --- .../InterpKernelGeo2DEdgeArcCircle.hxx | 5 +- .../Geometric2D/InterpKernelGeo2DEdgeLin.cxx | 26 --- .../Geometric2D/InterpKernelGeo2DEdgeLin.hxx | 4 - .../InterpKernelGeo2DQuadraticPolygon.cxx | 207 ++++++++++-------- .../InterpKernelGeo2DQuadraticPolygon.hxx | 7 +- .../MEDCouplingUMesh_intersection.cxx | 5 +- 11 files changed, 177 insertions(+), 160 deletions(-) diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DComposedEdge.cxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DComposedEdge.cxx index f445f5e44..55eda23f6 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DComposedEdge.cxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DComposedEdge.cxx @@ -310,6 +310,29 @@ void ComposedEdge::dumpInXfigFile(std::ostream& stream, int resolution, const Bo (*iter)->dumpInXfigFile(stream,resolution,box); } +void ComposedEdge::dumpToCout(const std::map& mapp) const +{ + int i=0; + for(std::list::const_iterator iter=_sub_edges.begin();iter!=_sub_edges.end();iter++, i++) + { + const ElementaryEdge& e(*(*iter)); + auto sI(mapp.find(e.getStartNode())), eI(mapp.find(e.getEndNode())); + int start = (sI == mapp.end() ? -1 : sI->second), end = (eI == mapp.end() ? -1 : eI->second); + std::string locs; + switch (e.getLoc()) + { + case FULL_IN_1: locs="FULL_IN_1"; break; + case FULL_ON_1: locs="FULL_ON_1"; break; + case FULL_OUT_1: locs="FULL_OUT_1"; break; + case FULL_UNKNOWN: locs="FULL_UNKNOWN"; break; + default: locs="oh my God! This is so wrong."; + } + std::cout << "Edge [" << i << "] : ("<< std::hex << &e << std::dec << ") -> (" << start << ", " << end << ")\t" << locs << std::endl; + } + std::cout << std::endl; +} + + Node *ComposedEdge::getEndNode() const { return _sub_edges.back()->getEndNode(); diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DComposedEdge.hxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DComposedEdge.hxx index a475986e1..4d04aa7b9 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DComposedEdge.hxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DComposedEdge.hxx @@ -105,6 +105,7 @@ namespace INTERP_KERNEL INTERPKERNEL_EXPORT bool changeEndNodeWith(Node *node) const; INTERPKERNEL_EXPORT bool changeStartNodeWith(Node *node) const; INTERPKERNEL_EXPORT void dumpInXfigFile(std::ostream& stream, int resolution, const Bounds& box) const; + INTERPKERNEL_EXPORT void dumpToCout(const std::map& mapp) const; INTERPKERNEL_EXPORT bool isInOrOut(Node *nodeToTest) const; INTERPKERNEL_EXPORT bool isInOrOut2(Node *nodeToTest) const; INTERPKERNEL_EXPORT bool getDirection() const; diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdge.cxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdge.cxx index 15687f762..24051cbdb 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdge.cxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdge.cxx @@ -1034,3 +1034,28 @@ void Edge::sortIdsAbs(const std::vector& addNodes, const edgesThis.push_back(tmpp2[i+1]); } } + +void Edge::fillGlobalInfoAbs(bool direction, const std::map& mapThis, const std::map& mapOther, int offset1, int offset2, double fact, double baryX, double baryY, + std::vector& edgesThis, std::vector& addCoo, std::map mapAddCoo) const +{ + int tmp[2]; + _start->fillGlobalInfoAbs(mapThis,mapOther,offset1,offset2,fact,baryX,baryY,addCoo,mapAddCoo,tmp); + _end->fillGlobalInfoAbs(mapThis,mapOther,offset1,offset2,fact,baryX,baryY,addCoo,mapAddCoo,tmp+1); + if(direction) + { + edgesThis.push_back(tmp[0]); + edgesThis.push_back(tmp[1]); + } + else + { + edgesThis.push_back(tmp[1]); + edgesThis.push_back(tmp[0]); + } +} + +void Edge::fillGlobalInfoAbs2(const std::map& mapThis, const std::map& mapOther, int offset1, int offset2, double fact, double baryX, double baryY, + std::vector& edgesOther, std::vector& addCoo, std::map& mapAddCoo) const +{ + _start->fillGlobalInfoAbs2(mapThis,mapOther,offset1,offset2,fact,baryX,baryY,addCoo,mapAddCoo,edgesOther); + _end->fillGlobalInfoAbs2(mapThis,mapOther,offset1,offset2,fact,baryX,baryY,addCoo,mapAddCoo,edgesOther); +} diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdge.hxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdge.hxx index c847530f2..e5501ad1a 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdge.hxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdge.hxx @@ -280,11 +280,12 @@ namespace INTERP_KERNEL public: bool sortSubNodesAbs(const double *coo, std::vector& subNodes); void sortIdsAbs(const std::vector& addNodes, const std::map& mapp1, const std::map& mapp2, std::vector& edgesThis); - virtual void fillGlobalInfoAbs(bool direction, const std::map& mapThis, const std::map& mapOther, int offset1, int offset2, double fact, double baryX, double baryY, - std::vector& edgesThis, std::vector& addCoo, std::map mapAddCoo) const = 0; - virtual void fillGlobalInfoAbs2(const std::map& mapThis, const std::map& mapOther, int offset1, int offset2, double fact, double baryX, double baryY, - std::vector& edgesOther, std::vector& addCoo, std::map& mapAddCoo) const = 0; virtual Edge *buildEdgeLyingOnMe(Node *start, Node *end, bool direction=true) const = 0; + void fillGlobalInfoAbs(bool direction, const std::map& mapThis, const std::map& mapOther, int offset1, int offset2, double fact, double baryX, double baryY, + std::vector& edgesThis, std::vector& addCoo, std::map mapAddCoo) const; + void fillGlobalInfoAbs2(const std::map& mapThis, const std::map& mapOther, int offset1, int offset2, double fact, double baryX, double baryY, + std::vector& edgesOther, std::vector& addCoo, std::map& mapAddCoo) const; + protected: Edge():_cnt(1),_loc(FULL_UNKNOWN),_start(0),_end(0) { } virtual ~Edge(); diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.cxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.cxx index a0b3b9d23..39c86ad43 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.cxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.cxx @@ -839,28 +839,3 @@ void EdgeArcCircle::updateBounds() if(IsIn2Pi(_angle0,_angle,M_PI)) _bounds[0]=_center[0]-_radius; } - -void EdgeArcCircle::fillGlobalInfoAbs(bool direction, const std::map& mapThis, const std::map& mapOther, int offset1, int offset2, double fact, double baryX, double baryY, - std::vector& edgesThis, std::vector& addCoo, std::map mapAddCoo) const -{ - int tmp[2]; - _start->fillGlobalInfoAbs(mapThis,mapOther,offset1,offset2,fact,baryX,baryY,addCoo,mapAddCoo,tmp); - _end->fillGlobalInfoAbs(mapThis,mapOther,offset1,offset2,fact,baryX,baryY,addCoo,mapAddCoo,tmp+1); - if(direction) - { - edgesThis.push_back(tmp[0]); - edgesThis.push_back(tmp[1]); - } - else - { - edgesThis.push_back(tmp[1]); - edgesThis.push_back(tmp[0]); - } -} - -void EdgeArcCircle::fillGlobalInfoAbs2(const std::map& mapThis, const std::map& mapOther, int offset1, int offset2, double fact, double baryX, double baryY, - std::vector& edgesOther, std::vector& addCoo, std::map& mapAddCoo) const -{ - _start->fillGlobalInfoAbs2(mapThis,mapOther,offset1,offset2,fact,baryX,baryY,addCoo,mapAddCoo,edgesOther); - _end->fillGlobalInfoAbs2(mapThis,mapOther,offset1,offset2,fact,baryX,baryY,addCoo,mapAddCoo,edgesOther); -} diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.hxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.hxx index fc32223e2..ee5e3bb42 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.hxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.hxx @@ -124,10 +124,7 @@ namespace INTERP_KERNEL protected: void updateBounds(); Edge *buildEdgeLyingOnMe(Node *start, Node *end, bool direction=true) const; - void fillGlobalInfoAbs(bool direction, const std::map& mapThis, const std::map& mapOther, int offset1, int offset2, double fact, double baryX, double baryY, - std::vector& edgesThis, std::vector& addCoo, std::map mapAddCoo) const; - void fillGlobalInfoAbs2(const std::map& mapThis, const std::map& mapOther, int offset1, int offset2, double fact, double baryX, double baryY, - std::vector& edgesOther, std::vector& addCoo, std::map& mapAddCoo) const; + protected: //! Absolute angle where the arc starts. Value between -Pi and Pi double _angle0; diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeLin.cxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeLin.cxx index 724593474..39b95ba01 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeLin.cxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeLin.cxx @@ -155,7 +155,6 @@ void SegSegIntersector::areOverlappedOrOnlyColinears(const Bounds *whereToFind, double x=(*(_e1.getStartNode()))[0]-(*(_e2.getStartNode()))[0]; double y=(*(_e1.getStartNode()))[1]-(*(_e2.getStartNode()))[1]; // (x,y) is the vector between the two start points of e1 and e2 areOverlapped = fabs(-_matrix[0]*y+_matrix[1]*x) < dimChar*QuadraticPlanarPrecision::getPrecision(); // test colinearity of (x,y) with e1 - // explanation: if areOverlapped is true, we don't know yet if there will be an intersection (see meaning of areOverlapped in method doxy above) // if areOverlapped is false, we have two colinear vectors, not lying on the same line, so we're sure there is no intersec obviousNoIntersection = !areOverlapped; @@ -330,28 +329,3 @@ double EdgeLin::getCharactValueEng(const double *node) const double car1_1y=node[1]-(*(_start))[1]; double car1_2y=(*(_end))[1]-(*(_start))[1]; return (car1_1x*car1_2x+car1_1y*car1_2y)/(car1_2x*car1_2x+car1_2y*car1_2y); } - -void EdgeLin::fillGlobalInfoAbs(bool direction, const std::map& mapThis, const std::map& mapOther, int offset1, int offset2, double fact, double baryX, double baryY, - std::vector& edgesThis, std::vector& addCoo, std::map mapAddCoo) const -{ - int tmp[2]; - _start->fillGlobalInfoAbs(mapThis,mapOther,offset1,offset2,fact,baryX,baryY,addCoo,mapAddCoo,tmp); - _end->fillGlobalInfoAbs(mapThis,mapOther,offset1,offset2,fact,baryX,baryY,addCoo,mapAddCoo,tmp+1); - if(direction) - { - edgesThis.push_back(tmp[0]); - edgesThis.push_back(tmp[1]); - } - else - { - edgesThis.push_back(tmp[1]); - edgesThis.push_back(tmp[0]); - } -} - -void EdgeLin::fillGlobalInfoAbs2(const std::map& mapThis, const std::map& mapOther, int offset1, int offset2, double fact, double baryX, double baryY, - std::vector& edgesOther, std::vector& addCoo, std::map& mapAddCoo) const -{ - _start->fillGlobalInfoAbs2(mapThis,mapOther,offset1,offset2,fact,baryX,baryY,addCoo,mapAddCoo,edgesOther); - _end->fillGlobalInfoAbs2(mapThis,mapOther,offset1,offset2,fact,baryX,baryY,addCoo,mapAddCoo,edgesOther); -} diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeLin.hxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeLin.hxx index e721a0126..f76ba58f9 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeLin.hxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeLin.hxx @@ -77,10 +77,6 @@ namespace INTERP_KERNEL EdgeLin() { } void updateBounds(); Edge *buildEdgeLyingOnMe(Node *start, Node *end, bool direction) const; - void fillGlobalInfoAbs(bool direction, const std::map& mapThis, const std::map& mapOther, int offset1, int offset2, double fact, double baryX, double baryY, - std::vector& edgesThis, std::vector& addCoo, std::map mapAddCoo) const; - void fillGlobalInfoAbs2(const std::map& mapThis, const std::map& mapOther, int offset1, int offset2, double fact, double baryX, double baryY, - std::vector& edgesOther, std::vector& addCoo, std::map& mapAddCoo) const; }; } diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.cxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.cxx index bbd3d0097..9648a3ddc 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.cxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.cxx @@ -295,69 +295,79 @@ void QuadraticPolygon::splitAbs(QuadraticPolygon& other, double fact=normalizeExt(&other, xBaryBB, yBaryBB); // - IteratorOnComposedEdge it1(this),it3(&other); + IteratorOnComposedEdge itThis(this),itOther(&other); // other is (part of) the tool mesh MergePoints merge; - ComposedEdge *c1=new ComposedEdge; - ComposedEdge *c2=new ComposedEdge; + ComposedEdge *cThis=new ComposedEdge; + ComposedEdge *cOther=new ComposedEdge; int i=0; std::map mapAddCoo; - for(it3.first();!it3.finished();it3.next(),i++)//iteration over 'other->_sub_edges' + for(itOther.first();!itOther.finished();itOther.next(),i++) { + // For each edge of 'other', proceed with intersections: the edge might split into sub-edges, 'otherTmp' will hold the final split result. + // In the process of going through all edges of 'other', 'this' (which contains initially only one edge) + // is sub-divised into several edges : each of them has to be tested when intersecting the next candidate stored in 'other'. 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 'otherTmp->_sub_edge' + ElementaryEdge* curOther=itOther.current(); + otherTmp.pushBack(new ElementaryEdge(curOther->getPtr(),curOther->getDirection())); curOther->getPtr()->incrRef(); + IteratorOnComposedEdge itOtherTmp(&otherTmp); + for(itOtherTmp.first();!itOtherTmp.finished();itOtherTmp.next()) { - ElementaryEdge* curE2=it2.current(); - if(!curE2->isThereStartPoint()) - it1.first(); + ElementaryEdge* curOtherTmp=itOtherTmp.current(); + if(!curOtherTmp->isThereStartPoint()) + itThis.first(); // reset iterator on 'this' else - it1=curE2->getIterator(); - for(;!it1.finished();)//iteration over 'this' _sub_edges + itThis=curOtherTmp->getIterator(); + for(;!itThis.finished();) { - ElementaryEdge* curE1=it1.current(); + ElementaryEdge* curThis=itThis.current(); merge.clear(); // - std::map::const_iterator thisStart(mapThis.find(curE1->getStartNode())),thisEnd(mapThis.find(curE1->getEndNode())),otherStart(mapOther.find(curE2->getStartNode())),otherEnd(mapOther.find(curE2->getEndNode())); - int thisStart2(thisStart==mapThis.end()?-1:(*thisStart).second),thisEnd2(thisEnd==mapThis.end()?-1:(*thisEnd).second),otherStart2(otherStart==mapOther.end()?-1:(*otherStart).second+offset1),otherEnd2(otherEnd==mapOther.end()?-1:(*otherEnd).second+offset1); + std::map::const_iterator thisStart(mapThis.find(curThis->getStartNode())),thisEnd(mapThis.find(curThis->getEndNode())), + otherStart(mapOther.find(curOtherTmp->getStartNode())),otherEnd(mapOther.find(curOtherTmp->getEndNode())); + int thisStart2(thisStart==mapThis.end()?-1:(*thisStart).second), thisEnd2(thisEnd==mapThis.end()?-1:(*thisEnd).second), + otherStart2(otherStart==mapOther.end()?-1:(*otherStart).second+offset1),otherEnd2(otherEnd==mapOther.end()?-1:(*otherEnd).second+offset1); // - if(curE1->getPtr()->intersectWith(curE2->getPtr(),merge,*c1,*c2)) + if(curThis->getPtr()->intersectWith(curOtherTmp->getPtr(),merge,*cThis,*cOther)) { - if(!curE1->getDirection()) c1->reverse(); - if(!curE2->getDirection()) c2->reverse(); - UpdateNeighbours(merge,it1,it2,c1,c2); - //Substitution of simple edge by sub-edges. - delete curE1; // <-- destroying simple edge coming from pol1 - delete curE2; // <-- destroying simple edge coming from pol2 - it1.insertElemEdges(c1,true);// <-- 2nd param is true to go next. - it2.insertElemEdges(c2,false);// <-- 2nd param is false to avoid to go next. - curE2=it2.current(); + if(!curThis->getDirection()) cThis->reverse(); + if(!curOtherTmp->getDirection()) cOther->reverse(); + // Substitution of a single simple edge by two sub-edges resulting from the intersection + // First modify the edges currently pointed by itThis and itOtherTmp so that the newly created node + // becomes the end of the previous sub-edge and the beginning of the next one. + UpdateNeighbours(merge,itThis,itOtherTmp,cThis,cOther); + delete curThis; // <-- destroying simple edge coming from pol1 + delete curOtherTmp; // <-- destroying simple edge coming from pol2 + // Then insert second part of the intersection. + itThis.insertElemEdges(cThis,true); // <-- 2nd param is true to go next. + itOtherTmp.insertElemEdges(cOther,false); // <-- 2nd param is false to avoid to go next. + curOtherTmp=itOtherTmp.current(); // - it1.assignMySelfToAllElems(c2);//To avoid that others - SoftDelete(c1); - SoftDelete(c2); - c1=new ComposedEdge; - c2=new ComposedEdge; + itThis.assignMySelfToAllElems(cOther); + SoftDelete(cThis); + SoftDelete(cOther); + cThis=new ComposedEdge; + cOther=new ComposedEdge; } else { - UpdateNeighbours(merge,it1,it2,curE1,curE2); - it1.next(); + UpdateNeighbours(merge,itThis,itOtherTmp,curThis,curOtherTmp); + itThis.next(); } merge.updateMergedNodes(thisStart2,thisEnd2,otherStart2,otherEnd2,mergedNodes); } } + // If one sub-edge of otherTmp is "ON" an edge of this, then we have colinearity (all edges in otherTmp are //) if(otherTmp.presenceOfOn()) edgesInOtherColinearWithThis[otherEdgeIds[i]].push_back(cellIdThis); - if(otherTmp._sub_edges.size()>1) + // Converting back to integer connectivity: + if(otherTmp._sub_edges.size()>1) // only if a new point has been added (i.e. an actual intersection was done) { for(std::list::const_iterator it=otherTmp._sub_edges.begin();it!=otherTmp._sub_edges.end();it++) (*it)->fillGlobalInfoAbs2(mapThis,mapOther,offset1,offset2,/**/fact,xBaryBB,yBaryBB,/**/subDivOther[otherEdgeIds[i]],addCoo,mapAddCoo); } } - Delete(c1); - Delete(c2); + Delete(cThis); + Delete(cOther); // for(std::list::const_iterator it=_sub_edges.begin();it!=_sub_edges.end();it++) (*it)->fillGlobalInfoAbs(mapThis,mapOther,offset1,offset2,/**/fact,xBaryBB,yBaryBB,/**/edgesThis,addCoo,mapAddCoo); @@ -659,13 +669,16 @@ void QuadraticPolygon::appendCrudeData(const std::map * @param [in,out] edgesThis, parameter that keep informed the caller about the edges in this not shared by the result of intersection of \a this with \a other * @param [in,out] edgesBoundaryOther, parameter that stores all edges in result of intersection that are not */ -void QuadraticPolygon::buildPartitionsAbs(QuadraticPolygon& other, std::set& edgesThis, std::set& edgesBoundaryOther, const std::map& mapp, int idThis, int idOther, int offset, std::vector& addCoordsQuadratic, std::vector& conn, std::vector& connI, std::vector& nbThis, std::vector& nbOther) +void QuadraticPolygon::buildPartitionsAbs(QuadraticPolygon& other, std::set& edgesThis, std::set& edgesBoundaryOther, + const std::map& mapp, int idThis, int idOther, int offset, + std::vector& addCoordsQuadratic, std::vector& conn, std::vector& connI, + std::vector& nbThis, std::vector& nbOther) { double xBaryBB, yBaryBB; double fact=normalizeExt(&other, xBaryBB, yBaryBB); //Locate \a this relative to \a other (edges of \a this, aka \a pol1 are marked as IN or OUT) other.performLocatingOperationSlow(*this); // without any assumption - std::vector res=buildIntersectionPolygons(other,*this); + std::vector res=buildIntersectionPolygons(*this,other); for(std::vector::iterator it=res.begin();it!=res.end();it++) { (*it)->appendCrudeData(mapp,xBaryBB,yBaryBB,fact,offset,addCoordsQuadratic,conn,connI); @@ -889,6 +902,7 @@ void QuadraticPolygon::intersectForPoint(const QuadraticPolygon& other, std::vec * \b WARNING this method is const and other is const too. \b BUT location of Edges in 'this' and 'other' are nevertheless modified. * This is possible because loc attribute in Edge class is mutable. * This implies that if 'this' or/and 'other' are reused for intersect* method initLocations has to be called on each of this/them. + * This method is currently not used by any high level functionality. */ std::vector QuadraticPolygon::intersectMySelfWith(const QuadraticPolygon& other) const { @@ -897,7 +911,7 @@ std::vector QuadraticPolygon::intersectMySelfWith(const Quad SplitPolygonsEachOther(cpyOfThis,cpyOfOther,nbOfSplits); //At this point cpyOfThis and cpyOfOther have been splited at maximum edge so that in/out can been done. performLocatingOperation(cpyOfOther); - return other.buildIntersectionPolygons(cpyOfThis,cpyOfOther); + return other.buildIntersectionPolygons(cpyOfOther, cpyOfThis); } /*! @@ -976,32 +990,34 @@ void QuadraticPolygon::performLocatingOperationSlow(QuadraticPolygon& pol2) cons /*! * Given 2 polygons \a pol1 and \a pol2 (localized) the resulting polygons are returned. * - * this : pol2 simplified. + * this : pol1 simplified. * @param [in] pol1 pol1 split. * @param [in] pol2 pol2 split. */ std::vector QuadraticPolygon::buildIntersectionPolygons(const QuadraticPolygon& pol1, const QuadraticPolygon& pol2) const { std::vector ret; - std::list pol2Zip=pol2.zipConsecutiveInSegments(); - if(!pol2Zip.empty()) - ClosePolygons(pol2Zip,pol1,*this,ret); + // Extract from pol1, and pol1 only, all consecutive edges. + // pol1Zip contains concatenated pieces of pol1 which are part of the resulting intersecting cell being built. + std::list pol1Zip=pol1.zipConsecutiveInSegments(); + if(!pol1Zip.empty()) + ClosePolygons(pol1Zip,*this,pol2,ret); else - {//borders of pol2 do not cross pol1,and pol2 borders are outside of pol1. That is to say, either pol2 and pol1 - //do not overlap or pol1 is fully inside pol2. So in the first case no intersection, in the other case - //the intersection is pol1. - ElementaryEdge *e1FromPol1=pol1[0]; + {//borders of pol1 do not cross pol2,and pol1 borders are outside of pol2. That is to say, either pol1 and pol2 + //do not overlap or pol2 is fully inside pol1. So in the first case no intersection, in the other case + //the intersection is pol2. + ElementaryEdge *e1FromPol2=pol2[0]; TypeOfEdgeLocInPolygon loc=FULL_ON_1; - loc=e1FromPol1->locateFullyMySelf(*this,loc); + loc=e1FromPol2->locateFullyMySelf(*this,loc); if(loc==FULL_IN_1) - ret.push_back(new QuadraticPolygon(pol1)); + ret.push_back(new QuadraticPolygon(pol2)); } return ret; } /*! * Returns parts of potentially non closed-polygons. Each returned polygons are not mergeable. - * this : pol2 split and locallized. + * this : pol1 split and localized. */ std::list QuadraticPolygon::zipConsecutiveInSegments() const { @@ -1036,72 +1052,77 @@ std::list QuadraticPolygon::zipConsecutiveInSegments() const } /*! - * @param [in] pol2zip is a list of set of edges (=an opened polygon) coming from split polygon 2. - * @param [in] pol1 is split pol1. - * @param [in] pol2 should be considered as pol2Simplified. + * @param [in] pol1zip is a list of set of edges (=an opened polygon) coming from split polygon 1. + * @param [in] pol1 should be considered as pol1Simplified. + * @param [in] pol2 is split pol2. * @param [out] results the resulting \b CLOSED polygons. */ -void QuadraticPolygon::ClosePolygons(std::list& pol2Zip, const QuadraticPolygon& pol1, const QuadraticPolygon& pol2, +void QuadraticPolygon::ClosePolygons(std::list& pol1Zip, const QuadraticPolygon& pol1, const QuadraticPolygon& pol2, std::vector& results) { - bool directionKnownInPol1=false; - bool directionInPol1; - for(std::list::iterator iter=pol2Zip.begin();iter!=pol2Zip.end();) + bool directionKnownInPol2=false; + bool directionInPol2; + for(std::list::iterator iter=pol1Zip.begin();iter!=pol1Zip.end();) { + // Build incrementally the full closed cells from the consecutive line parts already built in pol1Zip. + // At the end of the process the item currently iterated has been totally completed (start_node=end_node) + // This process can produce several cells. if((*iter)->completed()) { results.push_back(*iter); - directionKnownInPol1=false; - iter=pol2Zip.erase(iter); + directionKnownInPol2=false; + iter=pol1Zip.erase(iter); continue; } - if(!directionKnownInPol1) + if(!directionKnownInPol2) { - if(!(*iter)->haveIAChanceToBeCompletedBy(pol1,pol2,directionInPol1)) - { delete *iter; iter=pol2Zip.erase(iter); continue; } + if(!(*iter)->haveIAChanceToBeCompletedBy(pol1,pol2,directionInPol2)) + { delete *iter; iter=pol1Zip.erase(iter); continue; } else - directionKnownInPol1=true; + directionKnownInPol2=true; } std::list::iterator iter2=iter; iter2++; - std::list::iterator iter3=(*iter)->fillAsMuchAsPossibleWith(pol1,iter2,pol2Zip.end(),directionInPol1); - if(iter3!=pol2Zip.end()) + // Fill as much as possible the current iterate (=a part of pol1) with consecutive pieces from pol2: + std::list::iterator iter3=(*iter)->fillAsMuchAsPossibleWith(pol2,iter2,pol1Zip.end(),directionInPol2); + // and now add a full connected piece from pol1Zip: + if(iter3!=pol1Zip.end()) { (*iter)->pushBack(*iter3); SoftDelete(*iter3); - pol2Zip.erase(iter3); + pol1Zip.erase(iter3); } } } /*! - * 'this' is expected to be set of edges (not closed) of pol2 split. + * 'this' is expected to be set of edges (not closed) of pol1 split. */ -bool QuadraticPolygon::haveIAChanceToBeCompletedBy(const QuadraticPolygon& pol1Splitted,const QuadraticPolygon& pol2NotSplitted, bool& direction) +bool QuadraticPolygon::haveIAChanceToBeCompletedBy(const QuadraticPolygon& pol1NotSplitted, const QuadraticPolygon& pol2Splitted, bool& direction) const { - IteratorOnComposedEdge it(const_cast(&pol1Splitted)); + IteratorOnComposedEdge it2(const_cast(&pol2Splitted)); bool found=false; Node *n=getEndNode(); - ElementaryEdge *cur=it.current(); - for(it.first();!it.finished() && !found;) + ElementaryEdge *cur=it2.current(); + for(it2.first();!it2.finished() && !found;) { - cur=it.current(); + cur=it2.current(); found=(cur->getStartNode()==n); if(!found) - it.next(); + it2.next(); } if(!found) throw Exception("Internal error: polygons incompatible with each others. Should never happen!"); - //Ok we found correspondence between this and pol1. Searching for right direction to close polygon. + //Ok we found correspondence between this and pol2. Searching for right direction to close polygon. ElementaryEdge *e=_sub_edges.back(); if(e->getLoc()==FULL_ON_1) { if(e->getPtr()==cur->getPtr()) { direction=false; - it.previousLoop(); - cur=it.current(); + it2.previousLoop(); + cur=it2.current(); Node *repr=cur->getPtr()->buildRepresentantOfMySelf(); - bool ret=pol2NotSplitted.isInOrOut(repr); + bool ret=pol1NotSplitted.isInOrOut(repr); repr->decrRef(); return ret; } @@ -1109,49 +1130,49 @@ bool QuadraticPolygon::haveIAChanceToBeCompletedBy(const QuadraticPolygon& pol1S { direction=true; Node *repr=cur->getPtr()->buildRepresentantOfMySelf(); - bool ret=pol2NotSplitted.isInOrOut(repr); + bool ret=pol1NotSplitted.isInOrOut(repr); repr->decrRef(); return ret; } } else - direction=cur->locateFullyMySelfAbsolute(pol2NotSplitted)==FULL_IN_1; + direction=cur->locateFullyMySelfAbsolute(pol1NotSplitted)==FULL_IN_1; return true; } /*! - * This method fills as much as possible \a this (a sub-part of pol2 split) with edges of \a pol1Splitted. + * This method fills as much as possible \a this (a sub-part of pol1 split) with edges of \a pol2Splitted. */ -std::list::iterator QuadraticPolygon::fillAsMuchAsPossibleWith(const QuadraticPolygon& pol1Splitted, +std::list::iterator QuadraticPolygon::fillAsMuchAsPossibleWith(const QuadraticPolygon& pol2Splitted, std::list::iterator iStart, std::list::iterator iEnd, bool direction) { - IteratorOnComposedEdge it(const_cast(&pol1Splitted)); + IteratorOnComposedEdge it1(const_cast(&pol2Splitted)); bool found=false; Node *n=getEndNode(); - ElementaryEdge *cur; - for(it.first();!it.finished() && !found;) + ElementaryEdge *cur1; + for(it1.first();!it1.finished() && !found;) { - cur=it.current(); - found=(cur->getStartNode()==n); + cur1=it1.current(); + found=(cur1->getStartNode()==n); if(!found) - it.next(); + it1.next(); } if(!direction) - it.previousLoop(); + it1.previousLoop(); Node *nodeToTest; - int szMax(pol1Splitted.size()+1),ii(0);// here a protection against aggressive users of IntersectMeshes of invalid input meshes + int szMax(pol2Splitted.size()+1),ii(0); // protection against aggressive users of IntersectMeshes using invalid input meshes ... std::list::iterator ret; do - { - cur=it.current(); - ElementaryEdge *tmp=cur->clone(); + { // Stack (consecutive) edges of pol1 into the result (no need to care about ordering, edges from pol1 are already consecutive) + cur1=it1.current(); + ElementaryEdge *tmp=cur1->clone(); if(!direction) tmp->reverse(); pushBack(tmp); nodeToTest=tmp->getEndNode(); - direction?it.nextLoop():it.previousLoop(); + direction?it1.nextLoop():it1.previousLoop(); ret=CheckInList(nodeToTest,iStart,iEnd); if(completed()) return iEnd; diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.hxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.hxx index f44f95b98..a5afe97e5 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.hxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.hxx @@ -90,17 +90,18 @@ namespace INTERP_KERNEL INTERPKERNEL_EXPORT void performLocatingOperationSlow(QuadraticPolygon& pol2) const; INTERPKERNEL_EXPORT static void SplitPolygonsEachOther(QuadraticPolygon& pol1, QuadraticPolygon& pol2, int& nbOfSplits); INTERPKERNEL_EXPORT std::vector buildIntersectionPolygons(const QuadraticPolygon& pol1, const QuadraticPolygon& pol2) const; - INTERPKERNEL_EXPORT bool haveIAChanceToBeCompletedBy(const QuadraticPolygon& pol1Splitted, const QuadraticPolygon& pol2NotSplitted, bool& direction); + 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); protected: std::list zipConsecutiveInSegments() const; void dumpInXfigFile(std::ostream& stream, int resolution, const Bounds& box) const; - static void ClosePolygons(std::list& pol2Zip, const QuadraticPolygon& pol1, const QuadraticPolygon& pol2, std::vector& results); + static void ClosePolygons(std::list& pol1Zip, const QuadraticPolygon& pol1, const QuadraticPolygon& pol2, + std::vector& results); template static void UpdateNeighbours(const MergePoints& merger, IteratorOnComposedEdge it1, IteratorOnComposedEdge it2, const EDGES *e1, const EDGES *e2); - std::list::iterator fillAsMuchAsPossibleWith(const QuadraticPolygon& pol1Splitted, + std::list::iterator fillAsMuchAsPossibleWith(const QuadraticPolygon& pol2Splitted, std::list::iterator iStart, std::list::iterator iEnd, bool direction); diff --git a/src/MEDCoupling/MEDCouplingUMesh_intersection.cxx b/src/MEDCoupling/MEDCouplingUMesh_intersection.cxx index 570cb1090..39e7cb4ee 100644 --- a/src/MEDCoupling/MEDCouplingUMesh_intersection.cxx +++ b/src/MEDCoupling/MEDCouplingUMesh_intersection.cxx @@ -1335,7 +1335,7 @@ void MEDCouplingUMesh::BuildIntersecting2DCellsFromEdges(double eps, const MEDCo const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(typ); // Populate mapp and mappRev with nodes from the current cell (i) from mesh1 - this also builds the Node* objects: MEDCouplingUMeshBuildQPFromMesh3(coo1,offset1,coo2,offset2,addCoords,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,/* output */mapp,mappRev); - // pol1 is the full cell from mesh2, in QP format, with all the additional intersecting nodes. + // pol1 is the full cell from mesh1, in QP format, with all the additional intersecting nodes. pol1.buildFromCrudeDataArray(mappRev,cm.isQuadratic(),conn1+connI1[i]+1,coo1, desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1); // @@ -1348,6 +1348,8 @@ void MEDCouplingUMesh::BuildIntersecting2DCellsFromEdges(double eps, const MEDCo std::map > edgesIn2ForShare; // common edges std::vector pol2s(candidates2.size()); int ii=0; + // Build, for each intersecting cell candidate from mesh2, the corresponding QP. + // Again all the additional intersecting nodes are there. for(std::vector::const_iterator it2=candidates2.begin();it2!=candidates2.end();it2++,ii++) { INTERP_KERNEL::NormalizedCellType typ2=(INTERP_KERNEL::NormalizedCellType)conn2[connI2[*it2]]; @@ -1359,6 +1361,7 @@ void MEDCouplingUMesh::BuildIntersecting2DCellsFromEdges(double eps, const MEDCo pol1,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,colinear2, /* output */ edgesIn2ForShare); } ii=0; + // Now rebuild intersected cells from all this: for(std::vector::const_iterator it2=candidates2.begin();it2!=candidates2.end();it2++,ii++) { INTERP_KERNEL::ComposedEdge::InitLocationsWithOther(pol1,pol2s[ii]); -- 2.39.2