-// Copyright (C) 2007-2013 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
// License as published by the Free Software Foundation; either
-// version 2.1 of the License.
+// version 2.1 of the License, or (at your option) any later version.
//
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
#include <sstream>
#include <algorithm>
+#include <limits>
using namespace INTERP_KERNEL;
return (getE1().getAngle()>0. && getE2().getAngle()>0.) || (getE1().getAngle()<0. && getE2().getAngle()<0.);
}
+bool ArcCArcCIntersector::areColinears() const
+{
+ double radiusL,radiusB;
+ double centerL[2],centerB[2];
+ double tmp,cst;
+ return internalAreColinears(getE1(),getE2(),tmp,cst,radiusL,centerL,radiusB,centerB);
+}
+
/*!
* Precondition 'start' and 'end' are on the same curve than this.
*/
return EdgeArcCircle::GetAbsoluteAngleOfNormalizedVect(((*node)[0]-getE1().getCenter()[0])/getE1().getRadius(),((*node)[1]-getE1().getCenter()[1])/getE1().getRadius());
}
-bool ArcCArcCIntersector::areArcsOverlapped(const EdgeArcCircle& a1, const EdgeArcCircle& a2)
+bool ArcCArcCIntersector::internalAreColinears(const EdgeArcCircle& a1, const EdgeArcCircle& a2, double& distBetweenCenters, double& cst,
+ double& radiusL, double centerL[2], double& radiusB, double centerB[2])
{
- double centerL[2],radiusL,angle0L,angleL;
- double centerB[2],radiusB;
double lgth1=fabs(a1.getAngle()*a1.getRadius());
double lgth2=fabs(a2.getAngle()*a2.getRadius());
if(lgth1<lgth2)
{//a1 is the little one ('L') and a2 the big one ('B')
- a1.getCenter(centerL); radiusL=a1.getRadius(); angle0L=a1.getAngle0(); angleL=a1.getAngle();
+ a1.getCenter(centerL); radiusL=a1.getRadius();
a2.getCenter(centerB); radiusB=a2.getRadius();
}
else
{
- a2.getCenter(centerL); radiusL=a2.getRadius(); angle0L=a2.getAngle0(); angleL=a2.getAngle();
+ a2.getCenter(centerL); radiusL=a2.getRadius();
a1.getCenter(centerB); radiusB=a1.getRadius();
}
- // dividing from the begining by radiusB^2 to keep precision
- double tmp=Node::distanceBtw2PtSq(centerL,centerB);
- double cst=tmp/(radiusB*radiusB);
+ // dividing from the beginning by radiusB^2 to keep precision
+ distBetweenCenters=Node::distanceBtw2PtSq(centerL,centerB);
+ cst=distBetweenCenters/(radiusB*radiusB);
cst+=radiusL*radiusL/(radiusB*radiusB);
- if(!Node::areDoubleEqualsWP(cst,1.,2.))
+ return Node::areDoubleEqualsWPRight(cst,1.,2.);
+}
+
+bool ArcCArcCIntersector::areArcsOverlapped(const EdgeArcCircle& a1, const EdgeArcCircle& a2)
+{
+ double radiusL,radiusB;
+ double centerL[2],centerB[2];
+ double tmp(0.),cst(0.);
+ if(!internalAreColinears(a1,a2,tmp,cst,radiusL,centerL,radiusB,centerB))
return false;
//
+ double angle0L,angleL;
Bounds *merge=a1.getBounds().nearlyAmIIntersectingWith(a2.getBounds());
merge->getInterceptedArc(centerL,radiusL,angle0L,angleL);
delete merge;
//
tmp=sqrt(tmp);
- if(Node::areDoubleEqualsWP(tmp,0.,1/(10*std::max(radiusL,radiusB))))
- {
- if(Node::areDoubleEquals(radiusL,radiusB))
- return true;
- else
- return false;
- }
+ if(Node::areDoubleEqualsWPLeft(tmp,0.,10*std::max(radiusL,radiusB)))
+ return Node::areDoubleEquals(radiusL,radiusB);
double phi=EdgeArcCircle::GetAbsoluteAngleOfNormalizedVect((centerL[0]-centerB[0])/tmp,(centerL[1]-centerB[1])/tmp);
double cst2=2*radiusL*tmp/(radiusB*radiusB);
double cmpContainer[4];
if(EdgeArcCircle::IsIn2Pi(angle0L,angleL,a))
cmpContainer[sizeOfCmpContainer++]=cst-cst2;
a=*std::max_element(cmpContainer,cmpContainer+sizeOfCmpContainer);
- return Node::areDoubleEqualsWP(a,1.,2.);
+ 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();
- if(_dist>radius1+radius2+QUADRATIC_PLANAR::_precision || _dist+std::min(radius1,radius2)+QUADRATIC_PLANAR::_precision<std::max(radius1,radius2))
+ if(_dist>radius1+radius2+QuadraticPlanarPrecision::getPrecision() || _dist+std::min(radius1,radius2)+QuadraticPlanarPrecision::getPrecision()<std::max(radius1,radius2))
{
obviousNoIntersection=true;
areOverlapped=false;
}
}
+/**
+ Heart of the algorithm for arc/arc intersection.
+ See http://mathworld.wolfram.com/Circle-CircleIntersection.html
+ The computation is done in the coordinate system where Ox is the line between the 2 circle centers.
+*/
std::list< IntersectElement > ArcCArcCIntersector::getIntersectionsCharacteristicVal() const
{
std::list< IntersectElement > ret;
const double *center1=getE1().getCenter();
const double *center2=getE2().getCenter();
double radius1=getE1().getRadius(); double radius2=getE2().getRadius();
- double d1_1=(_dist*_dist-radius2*radius2+radius1*radius1)/(2.*_dist);
+ double d1_1=(_dist*_dist-radius2*radius2+radius1*radius1)/(2.*_dist); // computation of 'x' on wolfram
double u[2];//u is normalized vector from center1 to center2.
u[0]=(center2[0]-center1[0])/_dist; u[1]=(center2[1]-center1[1])/_dist;
- double d1_1y=EdgeArcCircle::SafeSqrt(radius1*radius1-d1_1*d1_1);
+ double d1_1y=EdgeArcCircle::SafeSqrt(radius1*radius1-d1_1*d1_1); // computation of 'y' on wolfram
double angleE1=EdgeArcCircle::NormalizeAngle(getE1().getAngle0()+getE1().getAngle());
double angleE2=EdgeArcCircle::NormalizeAngle(getE2().getAngle0()+getE2().getAngle());
if(!Node::areDoubleEquals(d1_1y,0))
{
//2 intersections
double v1[2],v2[2];
+ // coming back to our coordinate system:
v1[0]=u[0]*d1_1-u[1]*d1_1y; v1[1]=u[1]*d1_1+u[0]*d1_1y;
v2[0]=u[0]*d1_1+u[1]*d1_1y; v2[1]=u[1]*d1_1-u[0]*d1_1y;
Node *node1=new Node(center1[0]+v1[0],center1[1]+v1[1]); node1->declareOn();
v4[0]=center1[0]-center2[0]+v2[0]; v4[1]=center1[1]-center2[1]+v2[1];
double angle1_2=EdgeArcCircle::GetAbsoluteAngleOfNormalizedVect(v3[0]/radius2,v3[1]/radius2);
double angle2_2=EdgeArcCircle::GetAbsoluteAngleOfNormalizedVect(v4[0]/radius2,v4[1]/radius2);
+ // Check whether intersection points are exactly ON the other arc or not
+ // -> the curvilinear distance (=radius*angle) must below eps
+ bool e1_1S=Node::areDoubleEqualsWPLeft(angle1_1,getE1().getAngle0(),radius1);
+ bool e1_1E=Node::areDoubleEqualsWPLeft(angle1_1,angleE1,radius1);
+ bool e1_2S=Node::areDoubleEqualsWPLeft(angle1_2,getE2().getAngle0(),radius1);
+ bool e1_2E=Node::areDoubleEqualsWPLeft(angle1_2,angleE2,radius1);
//
- bool e1_1S=Node::areDoubleEqualsWP(angle1_1,getE1().getAngle0(),radius1);
- bool e1_1E=Node::areDoubleEqualsWP(angle1_1,angleE1,radius1);
- bool e1_2S=Node::areDoubleEqualsWP(angle1_2,getE2().getAngle0(),radius1);
- bool e1_2E=Node::areDoubleEqualsWP(angle1_2,angleE2,radius1);
- //
- bool e2_1S=Node::areDoubleEqualsWP(angle2_1,getE1().getAngle0(),radius2);
- bool e2_1E=Node::areDoubleEqualsWP(angle2_1,angleE1,radius2);
- bool e2_2S=Node::areDoubleEqualsWP(angle2_2,getE2().getAngle0(),radius2);
- bool e2_2E=Node::areDoubleEqualsWP(angle2_2,angleE2,radius2);
+ bool e2_1S=Node::areDoubleEqualsWPLeft(angle2_1,getE1().getAngle0(),radius2);
+ bool e2_1E=Node::areDoubleEqualsWPLeft(angle2_1,angleE1,radius2);
+ bool e2_2S=Node::areDoubleEqualsWPLeft(angle2_2,getE2().getAngle0(),radius2);
+ bool e2_2E=Node::areDoubleEqualsWPLeft(angle2_2,angleE2,radius2);
ret.push_back(IntersectElement(angle1_1,angle1_2,e1_1S,e1_1E,e1_2S,e1_2E,node1,_e1,_e2,keepOrder()));
ret.push_back(IntersectElement(angle2_1,angle2_2,e2_1S,e2_1E,e2_2S,e2_2E,node2,_e1,_e2,keepOrder()));
}
v2[0]=center1[0]-center2[0]+v1[0]; v2[1]=center1[1]-center2[1]+v1[1];
double angle0_1=EdgeArcCircle::GetAbsoluteAngleOfNormalizedVect(v1[0]/radius1,v1[1]/radius1);
double angle0_2=EdgeArcCircle::GetAbsoluteAngleOfNormalizedVect(v2[0]/radius2,v2[1]/radius2);
- bool e0_1S=Node::areDoubleEqualsWP(angle0_1,getE1().getAngle0(),radius1);
- bool e0_1E=Node::areDoubleEqualsWP(angle0_1,angleE1,radius1);
- bool e0_2S=Node::areDoubleEqualsWP(angle0_2,getE2().getAngle0(),radius2);
- bool e0_2E=Node::areDoubleEqualsWP(angle0_2,angleE2,radius2);
+ bool e0_1S=Node::areDoubleEqualsWPLeft(angle0_1,getE1().getAngle0(),radius1);
+ bool e0_1E=Node::areDoubleEqualsWPLeft(angle0_1,angleE1,radius1);
+ bool e0_2S=Node::areDoubleEqualsWPLeft(angle0_2,getE2().getAngle0(),radius2);
+ bool e0_2E=Node::areDoubleEqualsWPLeft(angle0_2,angleE2,radius2);
Node *node=new Node(center1[0]+d1_1*u[0],center1[1]+d1_1*u[1]); node->declareOnTangent();
ret.push_back(IntersectElement(angle0_1,angle0_2,e0_1S,e0_1E,e0_2S,e0_2E,node,_e1,_e2,keepOrder()));
}
}
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];
_drSq=_dx*_dx+_dy*_dy;
_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*QUADRATIC_PLANAR::_precision)//QUADRATIC_PLANAR::_precision*QUADRATIC_PLANAR::_precision*_drSq*_drSq/(2.*_dx*_dx))
+ ((*(_e2.getStartNode()))[0]-center[0])*((*(_e2.getEndNode()))[1]-center[1])-
+ ((*(_e2.getStartNode()))[1]-center[1])*((*(_e2.getEndNode()))[0]-center[0]);
+}
+
+/**
+ 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<double>::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;
+}
+
+/*!
+ * By construction, no chance that an arc of circle and line to be colinear.
+ */
+bool ArcCSegIntersector::areColinears() const
+{
+ return false;
}
void ArcCSegIntersector::getPlacements(Node *start, Node *end, TypeOfLocInEdge& whereStart, TypeOfLocInEdge& whereEnd, MergePoints& commonNode) const
{
std::list< IntersectElement > ret;
const double *center=getE1().getCenter();
- if(!(fabs(_determinant)<(2.*QUADRATIC_PLANAR::_precision)))//QUADRATIC_PLANAR::_precision*QUADRATIC_PLANAR::_precision*_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();
* @param deltaAngle in ]-2.*Pi;2.*Pi[
*/
EdgeArcCircle::EdgeArcCircle(Node *start, Node *end, const double *center, double radius, double angle0, double deltaAngle, bool direction):Edge(start,end,direction),_angle(deltaAngle),
- _angle0(angle0),_radius(radius)
+ _angle0(angle0),_radius(radius)
{
_center[0]=center[0];
_center[1]=center[1];
/*!
* 'eps' is expected to be > 0.
* 'conn' is of size 3. conn[0] is start id, conn[1] is end id and conn[2] is middle id.
- * 'offset' is typically the number of nodes already existing in global 2D curve mesh. Additionnal coords 'addCoo' ids will be put after the already existing.
+ * 'offset' is typically the number of nodes already existing in global 2D curve mesh. Additional coords 'addCoo' ids will be put after the already existing.
*/
void EdgeArcCircle::tesselate(const int *conn, int offset, double eps, std::vector<int>& newConn, std::vector<double>& addCoo) const
{
/*!
* Given a \b normalized vector defined by (ux,uy) returns its angle in ]-Pi;Pi].
- * So before using this method ux*ux+uy*uy should as much as possible close to 1.
- * This methods is quite time consuming in order to keep as much as possible precision.
- * It is NOT ALWAYS possible to do that only in one call of acos. Sometimes call to asin is necessary
- * due to imperfection of acos near 0. and Pi (cos x ~ 1-x*x/2.)
+ * Actually in the current implementation, the vector does not even need to be normalized ...
*/
double EdgeArcCircle::GetAbsoluteAngleOfNormalizedVect(double ux, double uy)
{
- //When arc is lower than 0.707 Using Asin
- if(fabs(ux)<0.707)
- {
- double ret=SafeAcos(ux);
- if(uy>0.)
- return ret;
- ret=-ret;
- return ret;
- }
- else
- {
- double ret=SafeAsin(uy);
- if(ux>0.)
- return ret;
- if(ret>0.)
- return M_PI-ret;
- else
- return -M_PI-ret;
- }
+ return atan2(uy, ux);
}
void EdgeArcCircle::GetArcOfCirclePassingThru(const double *start, const double *middle, const double *end,
angleInRad0=GetAbsoluteAngleOfNormalizedVect((start[0]-center[0])/radius,(start[1]-center[1])/radius);
double angleInRadM=GetAbsoluteAngleOfNormalizedVect((middle[0]-center[0])/radius,(middle[1]-center[1])/radius);
angleInRad=GetAbsoluteAngleOfNormalizedVect(((start[0]-center[0])*(end[0]-center[0])+(start[1]-center[1])*(end[1]-center[1]))/(radius*radius),
- ((start[0]-center[0])*(end[1]-center[1])-(start[1]-center[1])*(end[0]-center[0]))/(radius*radius));
+ ((start[0]-center[0])*(end[1]-center[1])-(start[1]-center[1])*(end[0]-center[0]))/(radius*radius));
if(IsAngleNotIn(angleInRad0,angleInRad,angleInRadM))
angleInRad=angleInRad<0?2*M_PI+angleInRad:angleInRad-2*M_PI;
}
double tmp3=cos(angle1);
double tmp4=cos(_angle0);
bary[0]=_radius*x0*y0*(tmp4-tmp3)+_radius*_radius*(y0*(cos(2*_angle0)-cos(2*angle1))/4.+
- x0*(_angle/2.+(sin(2.*_angle0)-sin(2.*angle1))/4.))
- +tmp2*(tmp1*tmp1*tmp1-tmp0*tmp0*tmp0)/3.;
+ x0*(_angle/2.+(sin(2.*_angle0)-sin(2.*angle1))/4.))
+ +tmp2*(tmp1*tmp1*tmp1-tmp0*tmp0*tmp0)/3.;
bary[1]=y0*y0*_radius*(tmp4-tmp3)/2.+_radius*_radius*y0*(_angle/2.+(sin(2.*_angle0)-sin(2.*angle1))/4.)
- +tmp2*(tmp4-tmp3+(tmp3*tmp3*tmp3-tmp4*tmp4*tmp4)/3.)/2.;
+ +tmp2*(tmp4-tmp3+(tmp3*tmp3*tmp3-tmp4*tmp4*tmp4)/3.)/2.;
}
+/**
+ * Compute the "middle" of two points on the arc of circle.
+ * The order (p1,p2) or (p2,p1) doesn't matter. p1 and p2 have to be localized on the edge defined by this.
+ * \param[out] mid the point located half-way between p1 and p2 on the arc defined by this.
+ * \sa getMiddleOfPointsOriented() a generalisation working also when p1 and p2 are not on the arc.
+ */
void EdgeArcCircle::getMiddleOfPoints(const double *p1, const double *p2, double *mid) const
{
double dx1((p1[0]-_center[0])/_radius),dy1((p1[1]-_center[1])/_radius),dx2((p2[0]-_center[0])/_radius),dy2((p2[1]-_center[1])/_radius);
//
double myDelta1(angle1-_angle0),myDelta2(angle2-_angle0);
if(_angle>0.)
- { myDelta1=myDelta1>=0.?myDelta1:myDelta1+2.*M_PI; myDelta2=myDelta2>=0.?myDelta2:myDelta2+2.*M_PI; }
+ { myDelta1=myDelta1>-QuadraticPlanarPrecision::getPrecision()?myDelta1:myDelta1+2.*M_PI; myDelta2=myDelta2>-QuadraticPlanarPrecision::getPrecision()?myDelta2:myDelta2+2.*M_PI; }
else
- { myDelta1=myDelta1<=0.?myDelta1:myDelta1-2.*M_PI; myDelta2=myDelta2<=0.?myDelta2:myDelta2-2.*M_PI; }
+ { myDelta1=myDelta1<QuadraticPlanarPrecision::getPrecision()?myDelta1:myDelta1-2.*M_PI; myDelta2=myDelta2<QuadraticPlanarPrecision::getPrecision()?myDelta2:myDelta2-2.*M_PI; }
////
mid[0]=_center[0]+_radius*cos(_angle0+(myDelta1+myDelta2)/2.);
mid[1]=_center[1]+_radius*sin(_angle0+(myDelta1+myDelta2)/2.);
}
+/**
+ * Compute the "middle" of two points on the arc of circle.
+ * Walk on the circle from p1 to p2 using the rotation direction indicated by this->_angle (i.e. by the orientation of the arc).
+ * This function is sensitive to the ordering of p1 and p2.
+ * \param[out] mid the point located half-way between p1 and p2
+ * \sa getMiddleOfPoints() to be used when the order of p1 and p2 is not relevant.
+ */
+void EdgeArcCircle::getMiddleOfPointsOriented(const double *p1, const double *p2, double *mid) const
+{
+ double dx1((p1[0]-_center[0])/_radius),dy1((p1[1]-_center[1])/_radius),dx2((p2[0]-_center[0])/_radius),dy2((p2[1]-_center[1])/_radius);
+ double angle1(GetAbsoluteAngleOfNormalizedVect(dx1,dy1)),angle2(GetAbsoluteAngleOfNormalizedVect(dx2,dy2));
+
+ if (angle1 <= 0.0)
+ angle1 += 2.*M_PI;
+ if (angle2 <= 0.0)
+ angle2 += 2.*M_PI;
+
+ double avg;
+ if((_angle>0. && angle1 <= angle2) || (_angle<=0. && angle1 >= angle2))
+ avg = (angle1+angle2)/2.;
+ else
+ avg = (angle1+angle2)/2. - M_PI;
+
+ mid[0]=_center[0]+_radius*cos(avg);
+ mid[1]=_center[1]+_radius*sin(avg);
+}
+
+
/*!
* Characteristic value used is angle in ]_Pi;Pi[ from axe 0x.
*/
double myDelta2=val2-_angle0;
if(_angle>0.)
{
- myDelta1=myDelta1>-(_radius*QUADRATIC_PLANAR::_precision)?myDelta1:myDelta1+2.*M_PI;//in some cases val1 or val2 are so close to angle0 that myDelta is close to 0. but negative.
- myDelta2=myDelta2>-(_radius*QUADRATIC_PLANAR::_precision)?myDelta2:myDelta2+2.*M_PI;
+ myDelta1=myDelta1>-(_radius*QuadraticPlanarPrecision::getPrecision())?myDelta1:myDelta1+2.*M_PI;//in some cases val1 or val2 are so close to angle0 that myDelta is close to 0. but negative.
+ myDelta2=myDelta2>-(_radius*QuadraticPlanarPrecision::getPrecision())?myDelta2:myDelta2+2.*M_PI;
return myDelta1<myDelta2;
}
else
{
- myDelta1=myDelta1<(_radius*QUADRATIC_PLANAR::_precision)?myDelta1:myDelta1-2.*M_PI;
- myDelta2=myDelta2<(_radius*QUADRATIC_PLANAR::_precision)?myDelta2:myDelta2-2.*M_PI;
+ myDelta1=myDelta1<(_radius*QuadraticPlanarPrecision::getPrecision())?myDelta1:myDelta1-2.*M_PI;
+ myDelta2=myDelta2<(_radius*QuadraticPlanarPrecision::getPrecision())?myDelta2:myDelta2-2.*M_PI;
return myDelta2<myDelta1;
}
}
if(IsIn2Pi(_angle0,_angle,M_PI))
_bounds[0]=_center[0]-_radius;
}
-
-void EdgeArcCircle::fillGlobalInfoAbs(bool direction, const std::map<INTERP_KERNEL::Node *,int>& mapThis, const std::map<INTERP_KERNEL::Node *,int>& mapOther, int offset1, int offset2, double fact, double baryX, double baryY,
- std::vector<int>& edgesThis, std::vector<double>& addCoo, std::map<INTERP_KERNEL::Node *,int> 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<INTERP_KERNEL::Node *,int>& mapThis, const std::map<INTERP_KERNEL::Node *,int>& mapOther, int offset1, int offset2, double fact, double baryX, double baryY,
- std::vector<int>& edgesOther, std::vector<double>& addCoo, std::map<INTERP_KERNEL::Node *,int>& 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);
-}
