Simplifying GetAbsoluteAngle* functions. atan2 is robust.

diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DBounds.cxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DBounds.cxx
index 4f2915ed4b495c47bee1715b8123847c620a1677..9a3d39e728ad8ff437928dc516d550a1eb6a2a2a 100644 (file)
--- 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 0ca8f1f3dac180f03ac2c1300dbc0c983d23ac8a..b39a630c59b3e27e3bd5a680b8dd04653e18aa1f 100644 (file)
--- 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 19559f9bc08f1bb7f4f7afb8b921b889797d5d19..45fc84259b5fa425daa625ca8153485308a6bc94 100644 (file)
--- 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<int>& newConn, std::vector<double>& 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 56c0c894c0148d246e3fa058cf328c35093fe875..80d225f970c4fa11d6b324615e1fb0e996852e70 100644 (file)
--- 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 1bb16f76b14af6e8962e49291568bd6e50192e7d..c30481178e5cc8f0b861d40d010249beabb17353 100644 (file)
--- a/src/MEDCoupling/MEDCouplingUMesh.cxx
+++ b/src/MEDCoupling/MEDCouplingUMesh.cxx
@@ -7644,7 +7644,6 @@ bool MEDCouplingUMesh::BuildConvexEnvelopOf2DCellJarvis(const double *coords, co
             }
           std::vector<int> 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;