From 086b2cb05ad483dd98e0dfc32a8f744e7bf1be37 Mon Sep 17 00:00:00 2001 From: abn Date: Mon, 9 Apr 2018 17:24:59 +0200 Subject: [PATCH] Simplifying GetAbsoluteAngle* functions. atan2 is robust. --- .../Geometric2D/InterpKernelGeo2DBounds.cxx | 5 +- .../InterpKernelGeo2DEdgeArcCircle.cxx | 68 ++++++++----------- .../InterpKernelGeo2DEdgeArcCircle.hxx | 3 +- .../Geometric2D/InterpKernelGeo2DNode.cxx | 3 +- src/MEDCoupling/MEDCouplingUMesh.cxx | 3 +- 5 files changed, 34 insertions(+), 48 deletions(-) diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DBounds.cxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DBounds.cxx index 4f2915ed4..9a3d39e72 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DBounds.cxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DBounds.cxx @@ -116,16 +116,15 @@ void Bounds::getInterceptedArc(const double *center, double radius, double& intr w1[0]=v1[0]; w1[1]=_y_min-center[1]; w2[0]=v2[0]; w2[1]=_y_max-center[1]; double delta1=EdgeArcCircle::SafeAsin(v1[0]*v2[1]-v1[1]*v2[0]); double delta2=EdgeArcCircle::SafeAsin(w1[0]*w2[1]-w1[1]*w2[0]); - double tmp; if(fabs(delta1)>fabs(delta2)) { intrcptArcDelta=delta1; - intrcptArcAngle0=EdgeArcCircle::GetAbsoluteAngle(v1,tmp); + intrcptArcAngle0=EdgeArcCircle::GetAbsoluteAngleOfVector(v1[0],v1[1]); } else { intrcptArcDelta=delta2; - intrcptArcAngle0=EdgeArcCircle::GetAbsoluteAngle(w1,tmp); + intrcptArcAngle0=EdgeArcCircle::GetAbsoluteAngleOfVector(w1[0],w1[1]); } } } diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.cxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.cxx index 0ca8f1f3d..b39a630c5 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.cxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.cxx @@ -113,7 +113,7 @@ void ArcCArcCIntersector::getPlacements(Node *start, Node *end, TypeOfLocInEdge& */ double ArcCArcCIntersector::getAngle(Node *node) const { - return EdgeArcCircle::GetAbsoluteAngleOfNormalizedVect(((*node)[0]-getE1().getCenter()[0])/getE1().getRadius(),((*node)[1]-getE1().getCenter()[1])/getE1().getRadius()); + return EdgeArcCircle::GetAbsoluteAngleOfVector((*node)[0]-getE1().getCenter()[0],(*node)[1]-getE1().getCenter()[1]); } bool ArcCArcCIntersector::internalAreColinears(const EdgeArcCircle& a1, const EdgeArcCircle& a2, double& distBetweenCenters, double& cst, @@ -154,7 +154,7 @@ bool ArcCArcCIntersector::areArcsOverlapped(const EdgeArcCircle& a1, const EdgeA tmp=sqrt(tmp); 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 phi=EdgeArcCircle::GetAbsoluteAngleOfVector(centerL[0]-centerB[0],centerL[1]-centerB[1]); double cst2=2*radiusL*tmp/(radiusB*radiusB); double cmpContainer[4]; int sizeOfCmpContainer=2; @@ -218,13 +218,13 @@ std::list< IntersectElement > ArcCArcCIntersector::getIntersectionsCharacteristi 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(); Node *node2=new Node(center1[0]+v2[0],center1[1]+v2[1]); node2->declareOn(); - double angle1_1=EdgeArcCircle::GetAbsoluteAngleOfNormalizedVect(v1[0]/radius1,v1[1]/radius1); - double angle2_1=EdgeArcCircle::GetAbsoluteAngleOfNormalizedVect(v2[0]/radius1,v2[1]/radius1); + double angle1_1=EdgeArcCircle::GetAbsoluteAngleOfVector(v1[0],v1[1]); + double angle2_1=EdgeArcCircle::GetAbsoluteAngleOfVector(v2[0],v2[1]); double v3[2],v4[2]; v3[0]=center1[0]-center2[0]+v1[0]; v3[1]=center1[1]-center2[1]+v1[1]; 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); + double angle1_2=EdgeArcCircle::GetAbsoluteAngleOfVector(v3[0],v3[1]); + double angle2_2=EdgeArcCircle::GetAbsoluteAngleOfVector(v4[0],v4[1]); // 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); @@ -245,8 +245,8 @@ std::list< IntersectElement > ArcCArcCIntersector::getIntersectionsCharacteristi double v1[2],v2[2]; v1[0]=d1_1*u[0]; v1[1]=d1_1*u[1]; 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); + double angle0_1=EdgeArcCircle::GetAbsoluteAngleOfVector(v1[0],v1[1]); + double angle0_2=EdgeArcCircle::GetAbsoluteAngleOfVector(v2[0],v2[1]); 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); @@ -440,12 +440,12 @@ void EdgeArcCircle::changeMiddle(Node *newMiddle) Edge *EdgeArcCircle::buildEdgeLyingOnMe(Node *start, Node *end, bool direction) const { - double sx=((*start)[0]-_center[0])/_radius; - double sy=((*start)[1]-_center[1])/_radius; - double ex=((*end)[0]-_center[0])/_radius; - double ey=((*end)[1]-_center[1])/_radius; - double angle0=GetAbsoluteAngleOfNormalizedVect(direction?sx:ex,direction?sy:ey); - double deltaAngle=GetAbsoluteAngleOfNormalizedVect(sx*ex+sy*ey,sx*ey-sy*ex); + double sx=(*start)[0]-_center[0]; + double sy=(*start)[1]-_center[1]; + double ex=(*end)[0]-_center[0]; + double ey=(*end)[1]-_center[1]; + double angle0=GetAbsoluteAngleOfVector(direction?sx:ex,direction?sy:ey); + double deltaAngle=GetAbsoluteAngleOfVector(sx*ex+sy*ey,sx*ey-sy*ex); if(deltaAngle>0. && _angle<0.) deltaAngle-=2.*M_PI; else if(deltaAngle<0. && _angle>0.) @@ -517,21 +517,11 @@ EdgeArcCircle *EdgeArcCircle::BuildFromNodes(Node *start, Node *middle, Node *en } } -/*! - * Given an \b NON normalized vector 'vect', returns its norm 'normVect' and its - * angle in ]-Pi,Pi] relative to Ox axe. - */ -double EdgeArcCircle::GetAbsoluteAngle(const double *vect, double& normVect) -{ - normVect=Node::norm(vect); - return GetAbsoluteAngleOfNormalizedVect(vect[0]/normVect,vect[1]/normVect); -} - /*! * Given a \b normalized vector defined by (ux,uy) returns its angle in ]-Pi;Pi]. * Actually in the current implementation, the vector does not even need to be normalized ... */ -double EdgeArcCircle::GetAbsoluteAngleOfNormalizedVect(double ux, double uy) +double EdgeArcCircle::GetAbsoluteAngleOfVector(double ux, double uy) { return atan2(uy, ux); } @@ -545,10 +535,10 @@ void EdgeArcCircle::GetArcOfCirclePassingThru(const double *start, const double center[0]=((end[1]-middle[1])*b1+(start[1]-middle[1])*b2)/delta; center[1]=((middle[0]-end[0])*b1+(middle[0]-start[0])*b2)/delta; radius=SafeSqrt((start[0]-center[0])*(start[0]-center[0])+(start[1]-center[1])*(start[1]-center[1])); - 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)); + angleInRad0=GetAbsoluteAngleOfVector(start[0]-center[0],start[1]-center[1]); + double angleInRadM=GetAbsoluteAngleOfVector(middle[0]-center[0],middle[1]-center[1]); + angleInRad=GetAbsoluteAngleOfVector((start[0]-center[0])*(end[0]-center[0])+(start[1]-center[1])*(end[1]-center[1]), + (start[0]-center[0])*(end[1]-center[1])-(start[1]-center[1])*(end[0]-center[0])); if(IsAngleNotIn(angleInRad0,angleInRad,angleInRadM)) angleInRad=angleInRad<0?2*M_PI+angleInRad:angleInRad-2*M_PI; } @@ -643,8 +633,8 @@ void EdgeArcCircle::getBarycenterOfZone(double *bary) const */ 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 angle1(GetAbsoluteAngleOfNormalizedVect(dx1,dy1)),angle2(GetAbsoluteAngleOfNormalizedVect(dx2,dy2)); + double dx1(p1[0]-_center[0]),dy1(p1[1]-_center[1]),dx2(p2[0]-_center[0]),dy2(p2[1]-_center[1]); + double angle1(GetAbsoluteAngleOfVector(dx1,dy1)),angle2(GetAbsoluteAngleOfVector(dx2,dy2)); // double myDelta1(angle1-_angle0),myDelta2(angle2-_angle0); if(_angle>0.) @@ -665,8 +655,8 @@ void EdgeArcCircle::getMiddleOfPoints(const double *p1, const double *p2, double */ 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)); + double dx1(p1[0]-_center[0]),dy1(p1[1]-_center[1]),dx2(p2[0]-_center[0]),dy2(p2[1]-_center[1]); + double angle1(GetAbsoluteAngleOfVector(dx1,dy1)),angle2(GetAbsoluteAngleOfVector(dx2,dy2)); if (angle1 <= 0.0) angle1 += 2.*M_PI; @@ -724,16 +714,16 @@ bool EdgeArcCircle::isLower(double val1, double val2) const */ double EdgeArcCircle::getCharactValue(const Node& node) const { - double dx=(node[0]-_center[0])/_radius; - double dy=(node[1]-_center[1])/_radius; - return GetAbsoluteAngleOfNormalizedVect(dx,dy); + double dx=node[0]-_center[0]; + double dy=node[1]-_center[1]; + return GetAbsoluteAngleOfVector(dx,dy); } double EdgeArcCircle::getCharactValueBtw0And1(const Node& node) const { - double dx=(node[0]-_center[0])/_radius; - double dy=(node[1]-_center[1])/_radius; - double angle=GetAbsoluteAngleOfNormalizedVect(dx,dy); + double dx=node[0]-_center[0]; + double dy=node[1]-_center[1]; + double angle=GetAbsoluteAngleOfVector(dx,dy); // double myDelta=angle-_angle0; if(_angle>0.) diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.hxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.hxx index 19559f9bc..45fc84259 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.hxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.hxx @@ -104,8 +104,7 @@ namespace INTERP_KERNEL double getAngle() const { return _angle; } void tesselate(const int *conn, int offset, double eps, std::vector& newConn, std::vector& addCoo) const; static EdgeArcCircle *BuildFromNodes(Node *start, Node *middle, Node *end); - static double GetAbsoluteAngle(const double *vect, double& normVect); - static double GetAbsoluteAngleOfNormalizedVect(double ux, double uy); + static double GetAbsoluteAngleOfVector(double ux, double uy); static void GetArcOfCirclePassingThru(const double *start, const double *middle, const double *end, double *center, double& radius, double& angleInRad, double& angleInRad0); //! To avoid in aggressive optimizations nan. diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DNode.cxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DNode.cxx index 56c0c894c..80d225f97 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DNode.cxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DNode.cxx @@ -116,8 +116,7 @@ double Node::computeAngle(const double *pt1, const double *pt2) { double x=pt2[0]-pt1[0]; double y=pt2[1]-pt1[1]; - double norm=sqrt(x*x+y*y); - return EdgeArcCircle::GetAbsoluteAngleOfNormalizedVect(x/norm,y/norm); + return EdgeArcCircle::GetAbsoluteAngleOfVector(x,y); } /*! diff --git a/src/MEDCoupling/MEDCouplingUMesh.cxx b/src/MEDCoupling/MEDCouplingUMesh.cxx index 1bb16f76b..c30481178 100644 --- a/src/MEDCoupling/MEDCouplingUMesh.cxx +++ b/src/MEDCoupling/MEDCouplingUMesh.cxx @@ -7644,7 +7644,6 @@ bool MEDCouplingUMesh::BuildConvexEnvelopOf2DCellJarvis(const double *coords, co } std::vector tmpOut; tmpOut.reserve(sz); tmpOut.push_back(startNode); refX=1e300; - double tmp1; double tmp2[2]; double angle0=-M_PI/2; // @@ -7661,7 +7660,7 @@ bool MEDCouplingUMesh::BuildConvexEnvelopOf2DCellJarvis(const double *coords, co if(*node!=tmpOut.back() && *node!=prevNode) { tmp2[0]=coords[2*(*node)]-coords[2*tmpOut.back()]; tmp2[1]=coords[2*(*node)+1]-coords[2*tmpOut.back()+1]; - double angleM=INTERP_KERNEL::EdgeArcCircle::GetAbsoluteAngle(tmp2,tmp1); + double angleM=INTERP_KERNEL::EdgeArcCircle::GetAbsoluteAngleOfVector(tmp2[0], tmp2[1]); double res; if(angleM<=angle0) res=angle0-angleM; -- 2.39.2