X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;f=src%2FINTERP_KERNEL%2FGeometric2D%2FInterpKernelGeo2DEdgeArcCircle.cxx;h=f73b9e154295082f7cbe64569e2786a5132d962e;hb=9727e779d56acece98be02cdccd0f99cc5ef0fa2;hp=0ca8f1f3dac180f03ac2c1300dbc0c983d23ac8a;hpb=22180620a030d1e31592593fb396e3e2700986a7;p=tools%2Fmedcoupling.git diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.cxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.cxx index 0ca8f1f3d..f73b9e154 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.cxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.cxx @@ -1,4 +1,4 @@ -// Copyright (C) 2007-2016 CEA/DEN, EDF R&D +// Copyright (C) 2007-2019 CEA/DEN, EDF R&D // // This library is free software; you can redistribute it and/or // modify it under the terms of the GNU Lesser General Public @@ -26,6 +26,7 @@ #include #include +#include using namespace INTERP_KERNEL; @@ -170,7 +171,7 @@ bool ArcCArcCIntersector::areArcsOverlapped(const EdgeArcCircle& a1, const EdgeA return Node::areDoubleEqualsWPRight(a,1.,2.); } -void ArcCArcCIntersector::areOverlappedOrOnlyColinears(const Bounds *whereToFind, bool& obviousNoIntersection, bool& areOverlapped) +void ArcCArcCIntersector::areOverlappedOrOnlyColinears(bool& obviousNoIntersection, bool& areOverlapped) { _dist=Node::distanceBtw2Pt(getE1().getCenter(),getE2().getCenter()); double radius1=getE1().getRadius(); double radius2=getE2().getRadius(); @@ -320,13 +321,11 @@ std::list< IntersectElement > ArcCArcCIntersector::getIntersectionsCharacteristi } return ret;*/ -ArcCSegIntersector::ArcCSegIntersector(const EdgeArcCircle& e1, const EdgeLin& e2, bool reverse):CrossTypeEdgeIntersector(e1,e2,reverse) +ArcCSegIntersector::ArcCSegIntersector(const EdgeArcCircle& e1, const EdgeLin& e2, bool reverse): + CrossTypeEdgeIntersector(e1,e2,reverse), + _deltaRoot_div_dr(0.), + _i1S2E(false),_i1E2E(false) { -} - -void ArcCSegIntersector::areOverlappedOrOnlyColinears(const Bounds *whereToFind, bool& obviousNoIntersection, bool& areOverlapped) -{ - areOverlapped=false;//No overlapping by construction const double *center=getE1().getCenter(); _dx=(*(_e2.getEndNode()))[0]-(*(_e2.getStartNode()))[0]; _dy=(*(_e2.getEndNode()))[1]-(*(_e2.getStartNode()))[1]; @@ -334,11 +333,41 @@ void ArcCSegIntersector::areOverlappedOrOnlyColinears(const Bounds *whereToFind, _cross= ((*(_e2.getStartNode()))[0]-center[0])*((*(_e2.getEndNode()))[1]-center[1])- ((*(_e2.getStartNode()))[1]-center[1])*((*(_e2.getEndNode()))[0]-center[0]); - _determinant=getE1().getRadius()*getE1().getRadius()/_drSq-_cross*_cross/(_drSq*_drSq); - if(_determinant>-2*QuadraticPlanarPrecision::getPrecision())//QuadraticPlanarPrecision::getPrecision()*QuadraticPlanarPrecision::getPrecision()*_drSq*_drSq/(2.*_dx*_dx)) +} + +/** + See http://mathworld.wolfram.com/Circle-LineIntersection.html + _cross is 'D', the computation is done with the translation to put back the circle at the origin +*/ +void ArcCSegIntersector::areOverlappedOrOnlyColinears(bool& obviousNoIntersection, bool& areOverlapped) +{ + areOverlapped=false;//No overlapping by construction + + // Similar optimisation than SegSegIntersector::areOverlappedOrOnlyColinears() + bool dnu1, dnu2; + identifyEarlyIntersection(dnu1, dnu2, _i1S2E, _i1E2E); + + const double R = getE1().getRadius(); + + // We need to compute d = R*R-_cross*_cross/_drSq + // In terms of numerical precision, this can trigger 'catastrophic cancellation' and is hence better expressed as: + double _dr = sqrt(_drSq); + double diff = (R-_cross/_dr), add=(R+_cross/_dr); + // Ah ah: we will be taking a square root later. If we want the user to be able to use an epsilon finer than 1.0e-8, then we need + // to prevent ourselves going below machine precision (typ. 1.0e-16 for double). + const double eps_machine = std::numeric_limits::epsilon(); + diff = fabs(diff/R) < eps_machine ? 0.0 : diff; + add = fabs(add/R) < eps_machine ? 0.0 : add; + double d = add*diff; + // Compute deltaRoot_div_dr := sqrt(delta)/dr, where delta has the meaning of Wolfram. + // Then 2*deltaRoot_div_dr is the distance between the two intersection points of the line with the circle. This is what we compare to eps. + // We compute it in such a way that it can be used in boolean tests too (a very negative value means we're far apart from intersection) + _deltaRoot_div_dr = Node::sign(d)*sqrt(fabs(d)); + + if( 2*_deltaRoot_div_dr > -QuadraticPlanarPrecision::getPrecision()) obviousNoIntersection=false; else - obviousNoIntersection=true; + obviousNoIntersection=true; } /*! @@ -358,29 +387,71 @@ std::list< IntersectElement > ArcCSegIntersector::getIntersectionsCharacteristic { std::list< IntersectElement > ret; const double *center=getE1().getCenter(); - if(!(fabs(_determinant)<(2.*QuadraticPlanarPrecision::getPrecision())))//QuadraticPlanarPrecision::getPrecision()*QuadraticPlanarPrecision::getPrecision()*_drSq*_drSq/(2.*_dx*_dx)) - { - double determinant=EdgeArcCircle::SafeSqrt(_determinant); + if(!(2*fabs(_deltaRoot_div_dr) < QuadraticPlanarPrecision::getPrecision())) // see comments in areOverlappedOrOnlyColinears() + { // Two intersection nodes + // -> if a common node found, there is a chance that this is the only one (i.e. second intersection point is outside e1 and e2) + if(_earlyInter) + { + // Check tangent vector of the arc circle at the common node with the linear segment. + // There we can tell if the arc of circle is 'moving away' from the seg, or if it might intersect it twice + const Node &n(*_earlyInter->getNodeOnly()); + const double * center(getE1().getCenter()); + + double tang[2] = {-(n[1]-center[1]), n[0]-center[0]}; // (-y, x) is the tangent vector in the trigo direction with (x,y) = (center->node) + bool invSeg = _i1S2E || _i1E2E; + double linEdge[2] = {invSeg ? (-_dx) : _dx, invSeg ? (-_dy) : _dy}; + if(tang[1]*linEdge[0]-tang[0]*linEdge[1] < 0) + { + ret.push_back(*_earlyInter); + return ret; + } + } + + double determinant=fabs(_deltaRoot_div_dr)/sqrt(_drSq); double x1=(_cross*_dy/_drSq+Node::sign(_dy)*_dx*determinant)+center[0]; double y1=(-_cross*_dx/_drSq+fabs(_dy)*determinant)+center[1]; Node *intersect1=new Node(x1,y1); intersect1->declareOn(); - bool i1_1S=_e1.getStartNode()->isEqual(*intersect1); - bool i1_1E=_e1.getEndNode()->isEqual(*intersect1); - bool i1_2S=_e2.getStartNode()->isEqual(*intersect1); - bool i1_2E=_e2.getEndNode()->isEqual(*intersect1); - ret.push_back(IntersectElement(getE1().getCharactValue(*intersect1),getE2().getCharactValue(*intersect1),i1_1S,i1_1E,i1_2S,i1_2E,intersect1,_e1,_e2,keepOrder())); - // double x2=(_cross*_dy/_drSq-Node::sign(_dy)*_dx*determinant)+center[0]; double y2=(-_cross*_dx/_drSq-fabs(_dy)*determinant)+center[1]; Node *intersect2=new Node(x2,y2); intersect2->declareOn(); - bool i2_1S=_e1.getStartNode()->isEqual(*intersect2); - bool i2_1E=_e1.getEndNode()->isEqual(*intersect2); - bool i2_2S=_e2.getStartNode()->isEqual(*intersect2); - bool i2_2E=_e2.getEndNode()->isEqual(*intersect2); - ret.push_back(IntersectElement(getE1().getCharactValue(*intersect2),getE2().getCharactValue(*intersect2),i2_1S,i2_1E,i2_2S,i2_2E,intersect2,_e1,_e2,keepOrder())); + + bool isN1(false), isN2(false); + if (_earlyInter) + { + // Which node do we actually already found? Assume this is the closest ... + const Node &iN = *(_earlyInter->getNodeOnly()); + const Node &n1(*intersect1), &n2(*intersect2); + double d1 = std::max(fabs(iN[0]-n1[0]), fabs(iN[1]-n1[1])); + double d2 = std::max(fabs(iN[0]-n2[0]), fabs(iN[1]-n2[1])); + isN1 = d1 < d2; isN2 = !isN1; + if (isN1) intersect1->decrRef(); + if (isN2) intersect2->decrRef(); + ret.push_back(*_earlyInter); + } + if (!isN1) + { + bool i1_1S=_e1.getStartNode()->isEqual(*intersect1); + bool i1_1E=_e1.getEndNode()->isEqual(*intersect1); + bool i1_2S=_e2.getStartNode()->isEqual(*intersect1); + bool i1_2E=_e2.getEndNode()->isEqual(*intersect1); + ret.push_back(IntersectElement(getE1().getCharactValue(*intersect1),getE2().getCharactValue(*intersect1),i1_1S,i1_1E,i1_2S,i1_2E,intersect1,_e1,_e2,keepOrder())); + } + if(!isN2) + { + bool i2_1S=_e1.getStartNode()->isEqual(*intersect2); + bool i2_1E=_e1.getEndNode()->isEqual(*intersect2); + bool i2_2S=_e2.getStartNode()->isEqual(*intersect2); + bool i2_2E=_e2.getEndNode()->isEqual(*intersect2); + ret.push_back(IntersectElement(getE1().getCharactValue(*intersect2),getE2().getCharactValue(*intersect2),i2_1S,i2_1E,i2_2S,i2_2E,intersect2,_e1,_e2,keepOrder())); + } } else//tangent intersection { + if (_earlyInter) + { + ret.push_back(*_earlyInter); + return ret; + } double x=(_cross*_dy)/_drSq+center[0]; double y=(-_cross*_dx)/_drSq+center[1]; Node *intersect3=new Node(x,y); intersect3->declareOnTangent(); @@ -818,28 +889,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); -}