From 61efdf9ade5f0018da96eccdb3fcb05b43afe040 Mon Sep 17 00:00:00 2001 From: abn Date: Wed, 19 Dec 2018 16:11:30 +0100 Subject: [PATCH] To exchange with FB --- .../Geometric2D/InterpKernelGeo2DEdge.cxx | 25 ++++++ .../Geometric2D/InterpKernelGeo2DEdge.hxx | 9 ++- .../InterpKernelGeo2DEdgeArcCircle.cxx | 25 ------ .../InterpKernelGeo2DEdgeArcCircle.hxx | 5 +- .../Geometric2D/InterpKernelGeo2DEdgeLin.cxx | 81 +++++++++---------- .../Geometric2D/InterpKernelGeo2DEdgeLin.hxx | 4 - .../InterpKernelGeo2DQuadraticPolygon.cxx | 9 ++- .../MEDCouplingUMesh_intersection.cxx | 3 + 8 files changed, 79 insertions(+), 82 deletions(-) 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..57dde706f 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeLin.cxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeLin.cxx @@ -116,14 +116,22 @@ std::list< IntersectElement > SegSegIntersector::getIntersectionsCharacteristicV */ bool SegSegIntersector::areColinears() const { - Bounds b; - b.prepareForAggregation(); - b.aggregate(_e1.getBounds()); - b.aggregate(_e2.getBounds()); - double determinant=_matrix[0]*_matrix[3]-_matrix[1]*_matrix[2]; - double dimChar=b.getCaracteristicDim(); +// Bounds b; +// b.prepareForAggregation(); +// b.aggregate(_e1.getBounds()); +// b.aggregate(_e2.getBounds()); +// double dimChar=b.getCaracteristicDim(); - return fabs(determinant)< 2.*dimChar*QuadraticPlanarPrecision::getPrecision(); // same criteria as in areOverlappedOrOnlyColinears, see comment below + double determinant=_matrix[0]*_matrix[3]-_matrix[1]*_matrix[2]; + Bounds b1, b2; + b1.prepareForAggregation(); + b2.prepareForAggregation(); + b1.aggregate(_e1.getBounds()); + b2.aggregate(_e2.getBounds()); + double dimCharE1(b1.getCaracteristicDim()) ,dimCharE2(b2.getCaracteristicDim()); + +// return fabs(determinant)< 2.*dimChar*QuadraticPlanarPrecision::getPrecision(); // same criteria as in areOverlappedOrOnlyColinears, see comment below + return fabs(determinant)< 2.*dimCharE1*dimCharE2*QuadraticPlanarPrecision::getPrecision(); // same criteria as in areOverlappedOrOnlyColinears, see comment below } /*! @@ -137,25 +145,39 @@ bool SegSegIntersector::areColinears() const void SegSegIntersector::areOverlappedOrOnlyColinears(const Bounds *whereToFind, bool& obviousNoIntersection, bool& areOverlapped) { double determinant=_matrix[0]*_matrix[3]-_matrix[1]*_matrix[2]; - Bounds b; - b.prepareForAggregation(); - b.aggregate(_e1.getBounds()); - b.aggregate(_e2.getBounds()); - double dimChar=b.getCaracteristicDim(); + Bounds b1, b2; + b1.prepareForAggregation(); + b2.prepareForAggregation(); + b1.aggregate(_e1.getBounds()); + b2.aggregate(_e2.getBounds()); + double dimCharE1(b1.getCaracteristicDim()) ,dimCharE2(b2.getCaracteristicDim()); // Same criteria as in areColinears(), see doc. // [ABN] the 2 is not really justified, but the initial tests from Tony were written so closely to precision that I can't bother to change all of them ... - if(fabs(determinant)>2.*dimChar*QuadraticPlanarPrecision::getPrecision()) + if(fabs(determinant)>2.*dimCharE1*dimCharE2*QuadraticPlanarPrecision::getPrecision()) { obviousNoIntersection=false; areOverlapped=false; _matrix[0]/=determinant; _matrix[1]/=determinant; _matrix[2]/=determinant; _matrix[3]/=determinant; } else // colinear vectors { - 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 - + // Compute vectors joining tips of e1 and e2 + double xS=(*(_e1.getStartNode()))[0]-(*(_e2.getStartNode()))[0]; + double yS=(*(_e1.getStartNode()))[1]-(*(_e2.getStartNode()))[1]; + double xE=(*(_e1.getEndNode()))[0]-(*(_e2.getEndNode()))[0]; + double yE=(*(_e1.getEndNode()))[1]-(*(_e2.getEndNode()))[1]; + double maxDimS(std::max(fabs(xS),fabs(yS))), maxDimE(std::max(fabs(xE), fabs(yE))); + bool isS = (maxDimS > maxDimE), isE1 = (dimCharE1 >= dimCharE2); + double x = isS ? xS : xE; + double y = isS ? yS : yE; + unsigned shift = isE1 ? 0 : 2; + // test colinearity of the greatest tip-joining vector and greatest vector among {e1, e2} + areOverlapped = fabs(x*_matrix[1+shift]-y*_matrix[0+shift]) < dimCharE1*dimCharE2*QuadraticPlanarPrecision::getPrecision(); + if (areOverlapped) + { + std::cout << "ARE OVER\n"; + std::cout << maxDimS << " " << maxDimE << " " << dimCharE1 << " " << dimCharE2 << "\n"; + } // 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 +352,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 a27582d8f..6e8e5345b 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.cxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.cxx @@ -304,6 +304,7 @@ void QuadraticPolygon::splitAbs(QuadraticPolygon& other, 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. + std::cout << " with tool seg #" << otherEdgeIds[i] << "\n"; QuadraticPolygon otherTmp; ElementaryEdge* curOther=itOther.current(); otherTmp.pushBack(new ElementaryEdge(curOther->getPtr(),curOther->getDirection())); curOther->getPtr()->incrRef(); @@ -312,7 +313,7 @@ void QuadraticPolygon::splitAbs(QuadraticPolygon& other, { ElementaryEdge* curOtherTmp=itOtherTmp.current(); if(!curOtherTmp->isThereStartPoint()) - itThis.first(); + itThis.first(); // rewind iterator on 'this' else itThis=curOtherTmp->getIterator(); for(;!itThis.finished();) @@ -329,15 +330,17 @@ void QuadraticPolygon::splitAbs(QuadraticPolygon& other, { 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); - //Substitution of simple edge by sub-edges. 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(); // - itThis.assignMySelfToAllElems(cOther);//To avoid that others + itThis.assignMySelfToAllElems(cOther); SoftDelete(cThis); SoftDelete(cOther); cThis=new ComposedEdge; diff --git a/src/MEDCoupling/MEDCouplingUMesh_intersection.cxx b/src/MEDCoupling/MEDCouplingUMesh_intersection.cxx index 570cb1090..97f6711db 100644 --- a/src/MEDCoupling/MEDCouplingUMesh_intersection.cxx +++ b/src/MEDCoupling/MEDCouplingUMesh_intersection.cxx @@ -1236,6 +1236,7 @@ void MEDCouplingUMesh::Intersect1DMeshes(const MEDCouplingUMesh *m1Desc, const M int offset2(offset1+m2Desc->getNumberOfNodes()); for(int i=0;i candidates2; // edges of mesh2 candidate for intersection myTree.getIntersectingElems(bbox1+i*2*SPACEDIM,candidates2); if(!candidates2.empty()) // candidates2 holds edges from the second mesh potentially intersecting current edge i in mesh1 @@ -1631,6 +1632,8 @@ MEDCouplingUMesh *MEDCouplingUMesh::Intersect2DMeshes(const MEDCouplingUMesh *m1 revDesc1->decrRef(); revDescIndx1->decrRef(); revDesc2->decrRef(); revDescIndx2->decrRef(); MCAuto dd1(desc1),dd2(descIndx1),dd3(desc2),dd4(descIndx2); MCAuto dd5(m1Desc),dd6(m2Desc); +// subDiv2[4].push_back(12);subDiv2[4].push_back(14); +// subDiv2[5][0] = 12; subDiv2[5].push_back(15); // Step 2: re-order newly created nodes according to the ordering found in m2 std::vector< std::vector > intersectEdge2; -- 2.39.2