-// Copyright (C) 2007-2014 CEA/DEN, EDF R&D
+// Copyright (C) 2007-2016 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
a2.getCenter(centerL); radiusL=a2.getRadius();
a1.getCenter(centerB); radiusB=a1.getRadius();
}
- // dividing from the begining by radiusB^2 to keep precision
+ // dividing from the beginning by radiusB^2 to keep precision
distBetweenCenters=Node::distanceBtw2PtSq(centerL,centerB);
cst=distBetweenCenters/(radiusB*radiusB);
cst+=radiusL*radiusL/(radiusB*radiusB);
{
_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;
((*(_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))
+ if(_determinant>-2*QuadraticPlanarPrecision::getPrecision())//QuadraticPlanarPrecision::getPrecision()*QuadraticPlanarPrecision::getPrecision()*_drSq*_drSq/(2.*_dx*_dx))
obviousNoIntersection=false;
else
obviousNoIntersection=true;
{
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))
+ if(!(fabs(_determinant)<(2.*QuadraticPlanarPrecision::getPrecision())))//QuadraticPlanarPrecision::getPrecision()*QuadraticPlanarPrecision::getPrecision()*_drSq*_drSq/(2.*_dx*_dx))
{
double determinant=EdgeArcCircle::SafeSqrt(_determinant);
double x1=(_cross*_dy/_drSq+Node::sign(_dy)*_dx*determinant)+center[0];
/*!
* '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
{
*/
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,
+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;
}
}
