src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.cxx

   1 // Copyright (C) 2007-2015  CEA/DEN, EDF R&D
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License, or (at your option) any later version.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19 // Author : Anthony Geay (CEA/DEN)
  20
  21 #include "InterpKernelGeo2DEdgeArcCircle.hxx"
  22 #include "InterpKernelGeo2DEdgeLin.hxx"
  23 #include "InterpKernelException.hxx"
  24 #include "InterpKernelGeo2DNode.hxx"
  25 #include "NormalizedUnstructuredMesh.hxx"
  26
  27 #include <sstream>
  28 #include <algorithm>
  29
  30 using namespace INTERP_KERNEL;
  31
  32 ArcCArcCIntersector::ArcCArcCIntersector(const EdgeArcCircle& e1, const EdgeArcCircle& e2):SameTypeEdgeIntersector(e1,e2),_dist(0.)
  33 {
  34 }
  35
  36 bool ArcCArcCIntersector::haveTheySameDirection() const
  37 {
  38   return (getE1().getAngle()>0. &&  getE2().getAngle()>0.) || (getE1().getAngle()<0. &&  getE2().getAngle()<0.);
  39 }
  40
  41 bool ArcCArcCIntersector::areColinears() const
  42 {
  43   double radiusL,radiusB;
  44   double centerL[2],centerB[2];
  45   double tmp,cst;
  46   return internalAreColinears(getE1(),getE2(),tmp,cst,radiusL,centerL,radiusB,centerB);
  47 }
  48
  49 /*!
  50  * Precondition 'start' and 'end' are on the same curve than this.
  51  */
  52 void ArcCArcCIntersector::getPlacements(Node *start, Node *end, TypeOfLocInEdge& whereStart, TypeOfLocInEdge& whereEnd, MergePoints& commonNode) const
  53 {
  54   bool obvious1,obvious2;
  55   obviousCaseForCurvAbscisse(start,whereStart,commonNode,obvious1);
  56   obviousCaseForCurvAbscisse(end,whereEnd,commonNode,obvious2);
  57   if(obvious1 && obvious2)
  58     return ;
  59   double angleInRadStart=getAngle(start);
  60   double angleInRadEnd=getAngle(end);
  61   if(obvious1 || obvious2)
  62     {
  63       if(obvious1)
  64         {
  65           if(EdgeArcCircle::IsIn2Pi(getE1().getAngle0(),getE1().getAngle(),angleInRadEnd))
  66             whereEnd=INSIDE;
  67           else
  68             whereEnd=OUT_AFTER;
  69           return ;
  70         }
  71       else
  72         {
  73           if(EdgeArcCircle::IsIn2Pi(getE1().getAngle0(),getE1().getAngle(),angleInRadStart))
  74             whereStart=INSIDE;
  75           else
  76             whereStart=OUT_BEFORE;
  77           return ;
  78         }
  79     }
  80   if(EdgeArcCircle::IsIn2Pi(getE1().getAngle0(),getE1().getAngle(),angleInRadStart))
  81     {
  82       whereStart=INSIDE;
  83       if(EdgeArcCircle::IsIn2Pi(getE1().getAngle0(),getE1().getAngle(),angleInRadEnd))
  84         whereEnd=INSIDE;
  85       else
  86         whereEnd=OUT_AFTER;
  87     }
  88   else
  89     {//we are out in start.
  90       if(EdgeArcCircle::IsIn2Pi(getE1().getAngle0(),getE1().getAngle(),angleInRadEnd))
  91         {
  92           whereStart=OUT_BEFORE;
  93           whereEnd=INSIDE;
  94         }
  95       else
  96         {
  97           if(EdgeArcCircle::IsIn2Pi(getE2().getAngle0(),getE2().getAngle(),getE1().getAngle0()))
  98             {//_e2 contains stictly _e1
  99               whereStart=OUT_BEFORE;
 100               whereEnd=OUT_AFTER;
 101             }
 102           else
 103             {//_e2 is outside from _e1
 104               whereStart=OUT_BEFORE;
 105               whereEnd=OUT_BEFORE;
 106             }
 107         }
 108     }
 109 }
 110
 111 /*!
 112  * Return angle between ]-Pi;Pi[
 113  */
 114 double ArcCArcCIntersector::getAngle(Node *node) const
 115 {
 116   return EdgeArcCircle::GetAbsoluteAngleOfNormalizedVect(((*node)[0]-getE1().getCenter()[0])/getE1().getRadius(),((*node)[1]-getE1().getCenter()[1])/getE1().getRadius());
 117 }
 118
 119 bool ArcCArcCIntersector::internalAreColinears(const EdgeArcCircle& a1, const EdgeArcCircle& a2, double& distBetweenCenters, double& cst,
 120                                                double& radiusL, double centerL[2], double& radiusB, double centerB[2])
 121 {
 122   double lgth1=fabs(a1.getAngle()*a1.getRadius());
 123   double lgth2=fabs(a2.getAngle()*a2.getRadius());
 124   if(lgth1<lgth2)
 125     {//a1 is the little one ('L') and a2 the big one ('B')
 126       a1.getCenter(centerL); radiusL=a1.getRadius();
 127       a2.getCenter(centerB); radiusB=a2.getRadius();
 128     }
 129   else
 130     {
 131       a2.getCenter(centerL); radiusL=a2.getRadius();
 132       a1.getCenter(centerB); radiusB=a1.getRadius();
 133     }
 134   // dividing from the begining by radiusB^2 to keep precision
 135   distBetweenCenters=Node::distanceBtw2PtSq(centerL,centerB);
 136   cst=distBetweenCenters/(radiusB*radiusB);
 137   cst+=radiusL*radiusL/(radiusB*radiusB);
 138   return Node::areDoubleEqualsWP(cst,1.,2.);
 139 }
 140
 141 bool ArcCArcCIntersector::areArcsOverlapped(const EdgeArcCircle& a1, const EdgeArcCircle& a2)
 142 {
 143   double radiusL,radiusB;
 144   double centerL[2],centerB[2];
 145   double tmp(0.),cst(0.);
 146   if(!internalAreColinears(a1,a2,tmp,cst,radiusL,centerL,radiusB,centerB))
 147     return false;
 148   //
 149   double angle0L,angleL;
 150   Bounds *merge=a1.getBounds().nearlyAmIIntersectingWith(a2.getBounds());
 151   merge->getInterceptedArc(centerL,radiusL,angle0L,angleL);
 152   delete merge;
 153   //
 154   tmp=sqrt(tmp);
 155   if(Node::areDoubleEqualsWP(tmp,0.,1/(10*std::max(radiusL,radiusB))))
 156     return Node::areDoubleEquals(radiusL,radiusB);
 157   double phi=EdgeArcCircle::GetAbsoluteAngleOfNormalizedVect((centerL[0]-centerB[0])/tmp,(centerL[1]-centerB[1])/tmp);
 158   double cst2=2*radiusL*tmp/(radiusB*radiusB);
 159   double cmpContainer[4];
 160   int sizeOfCmpContainer=2;
 161   cmpContainer[0]=cst+cst2*cos(phi-angle0L);
 162   cmpContainer[1]=cst+cst2*cos(phi-angle0L+angleL);
 163   double a=EdgeArcCircle::NormalizeAngle(phi-angle0L);
 164   if(EdgeArcCircle::IsIn2Pi(angle0L,angleL,a))
 165     cmpContainer[sizeOfCmpContainer++]=cst+cst2;
 166   a=EdgeArcCircle::NormalizeAngle(phi-angle0L+M_PI);
 167   if(EdgeArcCircle::IsIn2Pi(angle0L,angleL,a))
 168     cmpContainer[sizeOfCmpContainer++]=cst-cst2;
 169   a=*std::max_element(cmpContainer,cmpContainer+sizeOfCmpContainer);
 170   return Node::areDoubleEqualsWP(a,1.,2.);
 171 }
 172
 173 void ArcCArcCIntersector::areOverlappedOrOnlyColinears(const Bounds *whereToFind, bool& obviousNoIntersection, bool& areOverlapped)
 174 {
 175   _dist=Node::distanceBtw2Pt(getE1().getCenter(),getE2().getCenter());
 176   double radius1=getE1().getRadius(); double radius2=getE2().getRadius();
 177   if(_dist>radius1+radius2+QUADRATIC_PLANAR::_precision || _dist+std::min(radius1,radius2)+QUADRATIC_PLANAR::_precision<std::max(radius1,radius2))
 178     {
 179       obviousNoIntersection=true;
 180       areOverlapped=false;
 181       return ;
 182     }
 183   if(areArcsOverlapped(getE1(),getE2()))//(Node::areDoubleEquals(_dist,0.) && Node::areDoubleEquals(radius1,radius2))
 184     {
 185       obviousNoIntersection=false;
 186       areOverlapped=true;
 187     }
 188   else
 189     {
 190       obviousNoIntersection=false;
 191       areOverlapped=false;
 192     }
 193 }
 194
 195 std::list< IntersectElement > ArcCArcCIntersector::getIntersectionsCharacteristicVal() const
 196 {
 197   std::list< IntersectElement > ret;
 198   const double *center1=getE1().getCenter();
 199   const double *center2=getE2().getCenter();
 200   double radius1=getE1().getRadius(); double radius2=getE2().getRadius();
 201   double d1_1=(_dist*_dist-radius2*radius2+radius1*radius1)/(2.*_dist);
 202   double u[2];//u is normalized vector from center1 to center2.
 203   u[0]=(center2[0]-center1[0])/_dist; u[1]=(center2[1]-center1[1])/_dist;
 204   double d1_1y=EdgeArcCircle::SafeSqrt(radius1*radius1-d1_1*d1_1);
 205   double angleE1=EdgeArcCircle::NormalizeAngle(getE1().getAngle0()+getE1().getAngle());
 206   double angleE2=EdgeArcCircle::NormalizeAngle(getE2().getAngle0()+getE2().getAngle());
 207   if(!Node::areDoubleEquals(d1_1y,0))
 208     {
 209       //2 intersections
 210       double v1[2],v2[2];
 211       v1[0]=u[0]*d1_1-u[1]*d1_1y; v1[1]=u[1]*d1_1+u[0]*d1_1y;
 212       v2[0]=u[0]*d1_1+u[1]*d1_1y; v2[1]=u[1]*d1_1-u[0]*d1_1y;
 213       Node *node1=new Node(center1[0]+v1[0],center1[1]+v1[1]); node1->declareOn();
 214       Node *node2=new Node(center1[0]+v2[0],center1[1]+v2[1]); node2->declareOn();
 215       double angle1_1=EdgeArcCircle::GetAbsoluteAngleOfNormalizedVect(v1[0]/radius1,v1[1]/radius1);
 216       double angle2_1=EdgeArcCircle::GetAbsoluteAngleOfNormalizedVect(v2[0]/radius1,v2[1]/radius1);
 217       double v3[2],v4[2];
 218       v3[0]=center1[0]-center2[0]+v1[0]; v3[1]=center1[1]-center2[1]+v1[1];
 219       v4[0]=center1[0]-center2[0]+v2[0]; v4[1]=center1[1]-center2[1]+v2[1];
 220       double angle1_2=EdgeArcCircle::GetAbsoluteAngleOfNormalizedVect(v3[0]/radius2,v3[1]/radius2);
 221       double angle2_2=EdgeArcCircle::GetAbsoluteAngleOfNormalizedVect(v4[0]/radius2,v4[1]/radius2);
 222       //
 223       bool e1_1S=Node::areDoubleEqualsWP(angle1_1,getE1().getAngle0(),radius1);
 224       bool e1_1E=Node::areDoubleEqualsWP(angle1_1,angleE1,radius1);
 225       bool e1_2S=Node::areDoubleEqualsWP(angle1_2,getE2().getAngle0(),radius1);
 226       bool e1_2E=Node::areDoubleEqualsWP(angle1_2,angleE2,radius1);
 227       //
 228       bool e2_1S=Node::areDoubleEqualsWP(angle2_1,getE1().getAngle0(),radius2);
 229       bool e2_1E=Node::areDoubleEqualsWP(angle2_1,angleE1,radius2);
 230       bool e2_2S=Node::areDoubleEqualsWP(angle2_2,getE2().getAngle0(),radius2);
 231       bool e2_2E=Node::areDoubleEqualsWP(angle2_2,angleE2,radius2);
 232       ret.push_back(IntersectElement(angle1_1,angle1_2,e1_1S,e1_1E,e1_2S,e1_2E,node1,_e1,_e2,keepOrder()));
 233       ret.push_back(IntersectElement(angle2_1,angle2_2,e2_1S,e2_1E,e2_2S,e2_2E,node2,_e1,_e2,keepOrder()));
 234     }
 235   else
 236     {
 237       //tangent intersection
 238       double v1[2],v2[2];
 239       v1[0]=d1_1*u[0]; v1[1]=d1_1*u[1];
 240       v2[0]=center1[0]-center2[0]+v1[0]; v2[1]=center1[1]-center2[1]+v1[1];
 241       double angle0_1=EdgeArcCircle::GetAbsoluteAngleOfNormalizedVect(v1[0]/radius1,v1[1]/radius1);
 242       double angle0_2=EdgeArcCircle::GetAbsoluteAngleOfNormalizedVect(v2[0]/radius2,v2[1]/radius2);
 243       bool e0_1S=Node::areDoubleEqualsWP(angle0_1,getE1().getAngle0(),radius1);
 244       bool e0_1E=Node::areDoubleEqualsWP(angle0_1,angleE1,radius1);
 245       bool e0_2S=Node::areDoubleEqualsWP(angle0_2,getE2().getAngle0(),radius2);
 246       bool e0_2E=Node::areDoubleEqualsWP(angle0_2,angleE2,radius2);
 247       Node *node=new Node(center1[0]+d1_1*u[0],center1[1]+d1_1*u[1]); node->declareOnTangent();
 248       ret.push_back(IntersectElement(angle0_1,angle0_2,e0_1S,e0_1E,e0_2S,e0_2E,node,_e1,_e2,keepOrder()));
 249     }
 250   return ret;
 251 }
 252 /*double angle0_2;
 253   double signDeltaAngle2;
 254   double d1_2;
 255   if(u[1]<0.)
 256   angle0_1=-angle0_1;
 257   if(d1_1>=0.)
 258   {
 259   if(_dist>radius1)
 260   {
 261   angle0_2=angle0_1+M_PI;
 262   signDeltaAngle2=-1.;
 263   }
 264   else
 265   {
 266   angle0_2=angle0_1;
 267   signDeltaAngle2=1.;
 268   }
 269   }
 270   else
 271   {
 272   angle0_1+=M_PI;
 273   angle0_2=angle0_1;
 274   signDeltaAngle2=1.;
 275   }
 276   angle0_1=NormalizeAngle(angle0_1);
 277   angle0_2=NormalizeAngle(angle0_2);
 278   double angleE1=NormalizeAngle(getE1().getAngle0()+getE1().getAngle());
 279   double angleE2=NormalizeAngle(getE2().getAngle0()+getE2().getAngle());
 280   if(!(Node::areDoubleEquals(d1_1,radius1) || Node::areDoubleEquals(d1_1,-radius1)) )
 281   {
 282   //2 intersections
 283   double deltaAngle1=EdgeArcCircle::SafeAcos(fabs(d1_1)/radius1); //owns to 0;Pi/2 by construction
 284   double deltaAngle2=EdgeArcCircle::SafeAcos(fabs(d1_2)/radius2); //owns to 0;Pi/2 by construction
 285   double angle1_1=NormalizeAngle(angle0_1+deltaAngle1);// Intersection 1 seen for _e1
 286   double angle2_1=NormalizeAngle(angle0_1-deltaAngle1);// Intersection 2 seen for _e1
 287   double angle1_2=NormalizeAngle(angle0_2+signDeltaAngle2*deltaAngle2);// Intersection 1 seen for _e2
 288   double angle2_2=NormalizeAngle(angle0_2-signDeltaAngle2*deltaAngle2);// Intersection 2 seen for _e2
 289   //
 290   bool e1_1S=Node::areDoubleEqualsWP(angle1_1,getE1().getAngle0(),radius1);
 291   bool e1_1E=Node::areDoubleEqualsWP(angle1_1,angleE1,radius1);
 292   bool e1_2S=Node::areDoubleEqualsWP(angle1_2,getE2().getAngle0(),radius1);
 293   bool e1_2E=Node::areDoubleEqualsWP(angle1_2,angleE2,radius1);
 294   //
 295   bool e2_1S=Node::areDoubleEqualsWP(angle2_1,getE1().getAngle0(),radius2);
 296   bool e2_1E=Node::areDoubleEqualsWP(angle2_1,angleE1,radius2);
 297   bool e2_2S=Node::areDoubleEqualsWP(angle2_2,getE2().getAngle0(),radius2);
 298   bool e2_2E=Node::areDoubleEqualsWP(angle2_2,angleE2,radius2);
 299   Node *node1=new Node(center1[0]+radius1*cos(angle1_1),center1[0]+radius1*sin(angle1_1)); node1->declareOn();
 300   Node *node2=new Node(center1[0]+radius1*cos(angle2_1),center1[0]+radius1*sin(angle2_1)); node2->declareOn();
 301   ret.push_back(IntersectElement(angle1_1,angle1_2,e1_1S,e1_1E,e1_2S,e1_2E,node1,_e1,_e2,keepOrder()));
 302   ret.push_back(IntersectElement(angle2_1,angle2_2,e2_1S,e2_1E,e2_2S,e2_2E,node2,_e1,_e2,keepOrder()));
 303   }
 304   else
 305   //tangent intersection
 306   {
 307   bool e0_1S=Node::areDoubleEqualsWP(angle0_1,getE1().getAngle0(),radius1);
 308   bool e0_1E=Node::areDoubleEqualsWP(angle0_1,angleE1,radius1);
 309   bool e0_2S=Node::areDoubleEqualsWP(angle0_2,getE2().getAngle0(),radius2);
 310   bool e0_2E=Node::areDoubleEqualsWP(angle0_2,angleE2,radius2);
 311   Node *node=new Node(center1[0]+radius1*cos(angle0_1),center1[0]+radius1*sin(angle0_1)); node->declareOnTangent();
 312   ret.push_back(IntersectElement(angle0_1,angle0_2,e0_1S,e0_1E,e0_2S,e0_2E,node,_e1,_e2,keepOrder()));
 313   }
 314   return ret;*/
 315
 316 ArcCSegIntersector::ArcCSegIntersector(const EdgeArcCircle& e1, const EdgeLin& e2, bool reverse):CrossTypeEdgeIntersector(e1,e2,reverse)
 317 {
 318 }
 319
 320 void ArcCSegIntersector::areOverlappedOrOnlyColinears(const Bounds *whereToFind, bool& obviousNoIntersection, bool& areOverlapped)
 321 {
 322   areOverlapped=false;//No overlapping by construction
 323   const double *center=getE1().getCenter();
 324   _dx=(*(_e2.getEndNode()))[0]-(*(_e2.getStartNode()))[0];
 325   _dy=(*(_e2.getEndNode()))[1]-(*(_e2.getStartNode()))[1];
 326   _drSq=_dx*_dx+_dy*_dy;
 327   _cross=
 328       ((*(_e2.getStartNode()))[0]-center[0])*((*(_e2.getEndNode()))[1]-center[1])-
 329       ((*(_e2.getStartNode()))[1]-center[1])*((*(_e2.getEndNode()))[0]-center[0]);
 330   _determinant=getE1().getRadius()*getE1().getRadius()/_drSq-_cross*_cross/(_drSq*_drSq);
 331   if(_determinant>-2*QUADRATIC_PLANAR::_precision)//QUADRATIC_PLANAR::_precision*QUADRATIC_PLANAR::_precision*_drSq*_drSq/(2.*_dx*_dx))
 332     obviousNoIntersection=false;
 333   else
 334     obviousNoIntersection=true;
 335 }
 336
 337 /*!
 338  * By construction, no chance that an arc of circle and line to be colinear.
 339  */
 340 bool ArcCSegIntersector::areColinears() const
 341 {
 342   return false;
 343 }
 344
 345 void ArcCSegIntersector::getPlacements(Node *start, Node *end, TypeOfLocInEdge& whereStart, TypeOfLocInEdge& whereEnd, MergePoints& commonNode) const
 346 {
 347   throw Exception("Internal error. Should never been called : no overlapping possible between arc of circle and a segment.");
 348 }
 349
 350 std::list< IntersectElement > ArcCSegIntersector::getIntersectionsCharacteristicVal() const
 351 {
 352   std::list< IntersectElement > ret;
 353   const double *center=getE1().getCenter();
 354   if(!(fabs(_determinant)<(2.*QUADRATIC_PLANAR::_precision)))//QUADRATIC_PLANAR::_precision*QUADRATIC_PLANAR::_precision*_drSq*_drSq/(2.*_dx*_dx))
 355     {
 356       double determinant=EdgeArcCircle::SafeSqrt(_determinant);
 357       double x1=(_cross*_dy/_drSq+Node::sign(_dy)*_dx*determinant)+center[0];
 358       double y1=(-_cross*_dx/_drSq+fabs(_dy)*determinant)+center[1];
 359       Node *intersect1=new Node(x1,y1); intersect1->declareOn();
 360       bool i1_1S=_e1.getStartNode()->isEqual(*intersect1);
 361       bool i1_1E=_e1.getEndNode()->isEqual(*intersect1);
 362       bool i1_2S=_e2.getStartNode()->isEqual(*intersect1);
 363       bool i1_2E=_e2.getEndNode()->isEqual(*intersect1);
 364       ret.push_back(IntersectElement(getE1().getCharactValue(*intersect1),getE2().getCharactValue(*intersect1),i1_1S,i1_1E,i1_2S,i1_2E,intersect1,_e1,_e2,keepOrder()));
 365       //
 366       double x2=(_cross*_dy/_drSq-Node::sign(_dy)*_dx*determinant)+center[0];
 367       double y2=(-_cross*_dx/_drSq-fabs(_dy)*determinant)+center[1];
 368       Node *intersect2=new Node(x2,y2); intersect2->declareOn();
 369       bool i2_1S=_e1.getStartNode()->isEqual(*intersect2);
 370       bool i2_1E=_e1.getEndNode()->isEqual(*intersect2);
 371       bool i2_2S=_e2.getStartNode()->isEqual(*intersect2);
 372       bool i2_2E=_e2.getEndNode()->isEqual(*intersect2);
 373       ret.push_back(IntersectElement(getE1().getCharactValue(*intersect2),getE2().getCharactValue(*intersect2),i2_1S,i2_1E,i2_2S,i2_2E,intersect2,_e1,_e2,keepOrder()));
 374     }
 375   else//tangent intersection
 376     {
 377       double x=(_cross*_dy)/_drSq+center[0];
 378       double y=(-_cross*_dx)/_drSq+center[1];
 379       Node *intersect3=new Node(x,y); intersect3->declareOnTangent();
 380       bool i_1S=_e1.getStartNode()->isEqual(*intersect3);
 381       bool i_1E=_e1.getEndNode()->isEqual(*intersect3);
 382       bool i_2S=_e2.getStartNode()->isEqual(*intersect3);
 383       bool i_2E=_e2.getEndNode()->isEqual(*intersect3);
 384       ret.push_back(IntersectElement(_e1.getCharactValue(*intersect3),_e2.getCharactValue(*intersect3),i_1S,i_1E,i_2S,i_2E,intersect3,_e1,_e2,keepOrder()));
 385     }
 386   return ret;
 387 }
 388
 389 EdgeArcCircle::EdgeArcCircle(std::istream& lineInXfig)
 390 {
 391   const unsigned NB_OF_SKIP_FIELDS=15;
 392   std::string tmpS;
 393   for(unsigned i=0;i<NB_OF_SKIP_FIELDS;i++)
 394     lineInXfig >> tmpS;
 395   _start=new Node(lineInXfig);
 396   Node *middle=new Node(lineInXfig);
 397   _end=new Node(lineInXfig);
 398   GetArcOfCirclePassingThru(*_start,*middle,*_end,_center,_radius,_angle,_angle0);
 399   middle->decrRef();
 400   updateBounds();
 401 }
 402
 403 EdgeArcCircle::EdgeArcCircle(Node *start, Node *middle, Node *end, bool direction):Edge(start,end, direction)
 404 {
 405   GetArcOfCirclePassingThru(*_start,*middle,*_end,_center,_radius,_angle,_angle0);
 406   updateBounds();
 407 }
 408
 409 EdgeArcCircle::EdgeArcCircle(double sX, double sY, double mX, double mY, double eX, double eY):Edge(sX,sY,eX,eY)
 410 {
 411   double middle[2]; middle[0]=mX; middle[1]=mY;
 412   GetArcOfCirclePassingThru(*_start,middle,*_end,_center,_radius,_angle,_angle0);
 413   updateBounds();
 414 }
 415
 416 /*!
 417  * @param angle0 in ]-Pi;Pi[
 418  * @param deltaAngle in ]-2.*Pi;2.*Pi[
 419  */
 420 EdgeArcCircle::EdgeArcCircle(Node *start, Node *end, const double *center, double radius, double angle0, double deltaAngle, bool direction):Edge(start,end,direction),_angle(deltaAngle),
 421     _angle0(angle0),_radius(radius)
 422 {
 423   _center[0]=center[0];
 424   _center[1]=center[1];
 425   updateBounds();
 426 }
 427
 428 void EdgeArcCircle::changeMiddle(Node *newMiddle)
 429 {
 430   GetArcOfCirclePassingThru(*_start,*newMiddle,*_end,_center,_radius,_angle,_angle0);
 431   updateBounds();
 432 }
 433
 434 Edge *EdgeArcCircle::buildEdgeLyingOnMe(Node *start, Node *end, bool direction) const
 435 {
 436   double sx=((*start)[0]-_center[0])/_radius;
 437   double sy=((*start)[1]-_center[1])/_radius;
 438   double ex=((*end)[0]-_center[0])/_radius;
 439   double ey=((*end)[1]-_center[1])/_radius;
 440   double angle0=GetAbsoluteAngleOfNormalizedVect(direction?sx:ex,direction?sy:ey);
 441   double deltaAngle=GetAbsoluteAngleOfNormalizedVect(sx*ex+sy*ey,sx*ey-sy*ex);
 442   if(deltaAngle>0. && _angle<0.)
 443     deltaAngle-=2.*M_PI;
 444   else if(deltaAngle<0. && _angle>0.)
 445     deltaAngle+=2.*M_PI;
 446   deltaAngle=direction?deltaAngle:-deltaAngle;
 447   return new EdgeArcCircle(start,end,_center,_radius,angle0,deltaAngle,direction);
 448 }
 449
 450 void EdgeArcCircle::applySimilarity(double xBary, double yBary, double dimChar)
 451 {
 452   Edge::applySimilarity(xBary,yBary,dimChar);
 453   _radius/=dimChar;
 454   _center[0]=(_center[0]-xBary)/dimChar;
 455   _center[1]=(_center[1]-yBary)/dimChar;
 456 }
 457
 458 void EdgeArcCircle::unApplySimilarity(double xBary, double yBary, double dimChar)
 459 {
 460   Edge::unApplySimilarity(xBary,yBary,dimChar);
 461   _radius*=dimChar;
 462   _center[0]=_center[0]*dimChar+xBary;
 463   _center[1]=_center[1]*dimChar+yBary;
 464 }
 465
 466 /*!
 467  * 'eps' is expected to be > 0.
 468  * 'conn' is of size 3. conn[0] is start id, conn[1] is end id and conn[2] is middle id.
 469  * '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.
 470  */
 471 void EdgeArcCircle::tesselate(const int *conn, int offset, double eps, std::vector<int>& newConn, std::vector<double>& addCoo) const
 472 {
 473   newConn.push_back(INTERP_KERNEL::NORM_POLYL);
 474   int nbOfSubDiv=(int)(fabs(_angle)/eps);
 475   if(nbOfSubDiv<=2)
 476     {
 477       newConn.push_back(conn[0]); newConn.push_back(conn[2]); newConn.push_back(conn[1]);
 478       return ;
 479     }
 480   double signOfAngle=_angle>0.?1.:-1.;
 481   int offset2=offset+((int)addCoo.size())/2;
 482   newConn.push_back(conn[0]);
 483   for(int i=1;i<nbOfSubDiv;i++,offset2++)
 484     {
 485       double angle=_angle0+i*eps*signOfAngle;
 486       newConn.push_back(offset2);
 487       addCoo.push_back(_center[0]+_radius*cos(angle)); addCoo.push_back(_center[1]+_radius*sin(angle));
 488     }
 489   newConn.push_back(conn[1]);
 490 }
 491
 492 EdgeArcCircle *EdgeArcCircle::BuildFromNodes(Node *start, Node *middle, Node *end)
 493 {
 494   EdgeLin *e1,*e2;
 495   e1=new EdgeLin(start,middle);
 496   e2=new EdgeLin(middle,end);
 497   SegSegIntersector inters(*e1,*e2);
 498   bool colinearity=inters.areColinears();
 499   delete e1; delete e2;
 500   if(colinearity)
 501     {
 502       start->decrRef(); middle->decrRef(); end->decrRef();
 503       return 0;
 504     }
 505   else
 506     {
 507       EdgeArcCircle *ret=new EdgeArcCircle(start,middle,end);
 508       start->decrRef(); middle->decrRef(); end->decrRef();
 509       return ret;
 510     }
 511 }
 512
 513 /*!
 514  * Given an \b NON normalized vector 'vect', returns its norm 'normVect' and its
 515  * angle in ]-Pi,Pi] relative to Ox axe.
 516  */
 517 double EdgeArcCircle::GetAbsoluteAngle(const double *vect, double& normVect)
 518 {
 519   normVect=Node::norm(vect);
 520   return GetAbsoluteAngleOfNormalizedVect(vect[0]/normVect,vect[1]/normVect);
 521 }
 522
 523 /*!
 524  * Given a \b normalized vector defined by (ux,uy) returns its angle in ]-Pi;Pi].
 525  * So before using this method ux*ux+uy*uy should as much as possible close to 1.
 526  * This methods is quite time consuming in order to keep as much as possible precision.
 527  * It is NOT ALWAYS possible to do that only in one call of acos. Sometimes call to asin is necessary
 528  * due to imperfection of acos near 0. and Pi (cos x ~ 1-x*x/2.)
 529  */
 530 double EdgeArcCircle::GetAbsoluteAngleOfNormalizedVect(double ux, double uy)
 531 {
 532   return atan2(uy, ux);
 533 }
 534
 535 void EdgeArcCircle::GetArcOfCirclePassingThru(const double *start, const double *middle, const double *end,
 536                                               double *center, double& radius, double& angleInRad, double& angleInRad0)
 537 {
 538   double delta=(middle[0]-start[0])*(end[1]-middle[1])-(end[0]-middle[0])*(middle[1]-start[1]);
 539   double b1=(middle[1]*middle[1]+middle[0]*middle[0]-start[0]*start[0]-start[1]*start[1])/2;
 540   double b2=(end[1]*end[1]+end[0]*end[0]-middle[0]*middle[0]-middle[1]*middle[1])/2;
 541   center[0]=((end[1]-middle[1])*b1+(start[1]-middle[1])*b2)/delta;
 542   center[1]=((middle[0]-end[0])*b1+(middle[0]-start[0])*b2)/delta;
 543   radius=SafeSqrt((start[0]-center[0])*(start[0]-center[0])+(start[1]-center[1])*(start[1]-center[1]));
 544   angleInRad0=GetAbsoluteAngleOfNormalizedVect((start[0]-center[0])/radius,(start[1]-center[1])/radius);
 545   double angleInRadM=GetAbsoluteAngleOfNormalizedVect((middle[0]-center[0])/radius,(middle[1]-center[1])/radius);
 546   angleInRad=GetAbsoluteAngleOfNormalizedVect(((start[0]-center[0])*(end[0]-center[0])+(start[1]-center[1])*(end[1]-center[1]))/(radius*radius),
 547       ((start[0]-center[0])*(end[1]-center[1])-(start[1]-center[1])*(end[0]-center[0]))/(radius*radius));
 548   if(IsAngleNotIn(angleInRad0,angleInRad,angleInRadM))
 549     angleInRad=angleInRad<0?2*M_PI+angleInRad:angleInRad-2*M_PI;
 550 }
 551
 552 void EdgeArcCircle::dumpInXfigFile(std::ostream& stream, bool direction, int resolution, const Bounds& box) const
 553 {
 554   stream << "5 1 0 1 ";
 555   fillXfigStreamForLoc(stream);
 556   stream << " 7 50 -1 -1 0.000 0 ";
 557   if( (direction && (-_angle)>=0) || (!direction && (-_angle)<0))
 558     stream << '0';//'0'
 559   else
 560     stream << '1';//'1'
 561   stream << " 1 0 ";
 562   stream << box.fitXForXFigD(_center[0],resolution) << " " << box.fitYForXFigD(_center[1],resolution) << " ";
 563   direction?_start->dumpInXfigFile(stream,resolution,box):_end->dumpInXfigFile(stream,resolution,box);
 564   Node *middle=buildRepresentantOfMySelf();
 565   middle->dumpInXfigFile(stream,resolution,box);
 566   middle->decrRef();
 567   direction?_end->dumpInXfigFile(stream,resolution,box):_start->dumpInXfigFile(stream,resolution,box);
 568   stream << std::endl << "1 1 2.00 120.00 180.00" << std::endl;
 569 }
 570
 571 void EdgeArcCircle::update(Node *m)
 572 {
 573   GetArcOfCirclePassingThru(*_start,*m,*_end,_center,_radius,_angle,_angle0);
 574   updateBounds();
 575 }
 576
 577 /*!
 578  * This methods computes :
 579  * \f[
 580  * \int_{Current Edge} -ydx
 581  * \f]
 582  */
 583 double EdgeArcCircle::getAreaOfZone() const
 584 {
 585   return -_radius*_radius*(sin(_angle)-_angle)/2.+((*_start)[0]-(*_end)[0])*((*_start)[1]+(*_end)[1])/2.;
 586 }
 587
 588 double EdgeArcCircle::getCurveLength() const
 589 {
 590   return fabs(_angle*_radius);
 591 }
 592
 593 void EdgeArcCircle::getBarycenter(double *bary) const
 594 {
 595   bary[0]=_center[0]+_radius*cos(_angle0+_angle/2.);
 596   bary[1]=_center[1]+_radius*sin(_angle0+_angle/2.);
 597 }
 598
 599 /*!
 600  * \f[
 601  * bary[0]=\int_{Current Edge} -yxdx
 602  * \f]
 603  * \f[
 604  * bary[1]=\int_{Current Edge} -\frac{y^{2}}{2}dx
 605  * \f]
 606  * To compute these 2 expressions in this class we have :
 607  * \f[
 608  * x=x_{0}+Radius \cdot cos(\theta)
 609  * \f]
 610  * \f[
 611  * y=y_{0}+Radius \cdot sin(\theta)
 612  * \f]
 613  * \f[
 614  * dx=-Radius \cdot sin(\theta) \cdot d\theta
 615  * \f]
 616  */
 617 void EdgeArcCircle::getBarycenterOfZone(double *bary) const
 618 {
 619   double x0=_center[0];
 620   double y0=_center[1];
 621   double angle1=_angle0+_angle;
 622   double tmp1=sin(angle1);
 623   double tmp0=sin(_angle0);
 624   double tmp2=_radius*_radius*_radius;
 625   double tmp3=cos(angle1);
 626   double tmp4=cos(_angle0);
 627   bary[0]=_radius*x0*y0*(tmp4-tmp3)+_radius*_radius*(y0*(cos(2*_angle0)-cos(2*angle1))/4.+
 628       x0*(_angle/2.+(sin(2.*_angle0)-sin(2.*angle1))/4.))
 629       +tmp2*(tmp1*tmp1*tmp1-tmp0*tmp0*tmp0)/3.;
 630   bary[1]=y0*y0*_radius*(tmp4-tmp3)/2.+_radius*_radius*y0*(_angle/2.+(sin(2.*_angle0)-sin(2.*angle1))/4.)
 631         +tmp2*(tmp4-tmp3+(tmp3*tmp3*tmp3-tmp4*tmp4*tmp4)/3.)/2.;
 632 }
 633
 634 /**
 635  * Compute the "middle" of two points on the arc of circle.
 636  * The order (p1,p2) or (p2,p1) doesn't matter. p1 and p2 have to be localized on the edge defined by this.
 637  * \param[out] mid the point located half-way between p1 and p2 on the arc defined by this.
 638  * \sa getMiddleOfPointsOriented() a generalisation working also when p1 and p2 are not on the arc.
 639  */
 640 void EdgeArcCircle::getMiddleOfPoints(const double *p1, const double *p2, double *mid) const
 641 {
 642   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);
 643   double angle1(GetAbsoluteAngleOfNormalizedVect(dx1,dy1)),angle2(GetAbsoluteAngleOfNormalizedVect(dx2,dy2));
 644   //
 645   double myDelta1(angle1-_angle0),myDelta2(angle2-_angle0);
 646   if(_angle>0.)
 647     { myDelta1=myDelta1>-QUADRATIC_PLANAR::_precision?myDelta1:myDelta1+2.*M_PI; myDelta2=myDelta2>-QUADRATIC_PLANAR::_precision?myDelta2:myDelta2+2.*M_PI; }
 648   else
 649     { myDelta1=myDelta1<QUADRATIC_PLANAR::_precision?myDelta1:myDelta1-2.*M_PI; myDelta2=myDelta2<QUADRATIC_PLANAR::_precision?myDelta2:myDelta2-2.*M_PI; }
 650   ////
 651   mid[0]=_center[0]+_radius*cos(_angle0+(myDelta1+myDelta2)/2.);
 652   mid[1]=_center[1]+_radius*sin(_angle0+(myDelta1+myDelta2)/2.);
 653 }
 654
 655 /**
 656  * Compute the "middle" of two points on the arc of circle.
 657  * Walk on the circle from p1 to p2 using the rotation direction indicated by this->_angle (i.e. by the orientation of the arc).
 658  * This function is sensitive to the ordering of p1 and p2.
 659  * \param[out] mid the point located half-way between p1 and p2
 660  * \sa getMiddleOfPoints() to be used when the order of p1 and p2 is not relevant.
 661  */
 662 void EdgeArcCircle::getMiddleOfPointsOriented(const double *p1, const double *p2, double *mid) const
 663 {
 664   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);
 665   double angle1(GetAbsoluteAngleOfNormalizedVect(dx1,dy1)),angle2(GetAbsoluteAngleOfNormalizedVect(dx2,dy2));
 666
 667   if (angle1 <= 0.0)
 668     angle1 += 2.*M_PI;
 669   if (angle2 <= 0.0)
 670     angle2 += 2.*M_PI;
 671
 672   double avg;
 673   if((_angle>0. && angle1 <= angle2) || (_angle<=0. && angle1 >= angle2))
 674     avg = (angle1+angle2)/2.;
 675   else
 676     avg = (angle1+angle2)/2. - M_PI;
 677
 678   mid[0]=_center[0]+_radius*cos(avg);
 679   mid[1]=_center[1]+_radius*sin(avg);
 680 }
 681
 682
 683 /*!
 684  * Characteristic value used is angle in ]_Pi;Pi[ from axe 0x.
 685  */
 686 bool EdgeArcCircle::isIn(double characterVal) const
 687 {
 688   return IsIn2Pi(_angle0,_angle,characterVal);
 689 }
 690
 691 Node *EdgeArcCircle::buildRepresentantOfMySelf() const
 692 {
 693   return new Node(_center[0]+_radius*cos(_angle0+_angle/2.),_center[1]+_radius*sin(_angle0+_angle/2.));
 694 }
 695
 696 /*!
 697  * Characteristic value used is angle in ]_Pi;Pi[ from axe 0x.
 698  * 'val1' and 'val2' have been detected previously as owning to this.
 699  */
 700 bool EdgeArcCircle::isLower(double val1, double val2) const
 701 {
 702   double myDelta1=val1-_angle0;
 703   double myDelta2=val2-_angle0;
 704   if(_angle>0.)
 705     {
 706       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.
 707       myDelta2=myDelta2>-(_radius*QUADRATIC_PLANAR::_precision)?myDelta2:myDelta2+2.*M_PI;
 708       return myDelta1<myDelta2;
 709     }
 710   else
 711     {
 712       myDelta1=myDelta1<(_radius*QUADRATIC_PLANAR::_precision)?myDelta1:myDelta1-2.*M_PI;
 713       myDelta2=myDelta2<(_radius*QUADRATIC_PLANAR::_precision)?myDelta2:myDelta2-2.*M_PI;
 714       return myDelta2<myDelta1;
 715     }
 716 }
 717
 718 /*!
 719  * For Arc circle the caract value is angle with Ox between -Pi and Pi.
 720  */
 721 double EdgeArcCircle::getCharactValue(const Node& node) const
 722 {
 723   double dx=(node[0]-_center[0])/_radius;
 724   double dy=(node[1]-_center[1])/_radius;
 725   return GetAbsoluteAngleOfNormalizedVect(dx,dy);
 726 }
 727
 728 double EdgeArcCircle::getCharactValueBtw0And1(const Node& node) const
 729 {
 730   double dx=(node[0]-_center[0])/_radius;
 731   double dy=(node[1]-_center[1])/_radius;
 732   double angle=GetAbsoluteAngleOfNormalizedVect(dx,dy);
 733   //
 734   double myDelta=angle-_angle0;
 735   if(_angle>0.)
 736     myDelta=myDelta>=0.?myDelta:myDelta+2.*M_PI;
 737   else
 738     myDelta=myDelta<=0.?myDelta:myDelta-2.*M_PI;
 739   return myDelta/_angle;
 740 }
 741
 742 double EdgeArcCircle::getDistanceToPoint(const double *pt) const
 743 {
 744   double angle=Node::computeAngle(_center,pt);
 745   if(IsIn2Pi(_angle0,_angle,angle))
 746     return fabs(Node::distanceBtw2Pt(_center,pt)-_radius);
 747   else
 748     {
 749       double dist1=Node::distanceBtw2Pt(*_start,pt);
 750       double dist2=Node::distanceBtw2Pt(*_end,pt);
 751       return std::min(dist1,dist2);
 752     }
 753 }
 754
 755 bool EdgeArcCircle::isNodeLyingOn(const double *coordOfNode) const
 756 {
 757   double dist=Node::distanceBtw2Pt(_center,coordOfNode);
 758   if(Node::areDoubleEquals(dist,_radius))
 759     {
 760       double angle=Node::computeAngle(_center,coordOfNode);
 761       return IsIn2Pi(_angle0,_angle,angle);
 762     }
 763   else
 764     return false;
 765 }
 766
 767 /*!
 768  * Idem IsAngleNotIn except that here 'start' in ]-Pi;Pi[ and delta in ]-2*Pi;2Pi[.
 769  * @param angleIn in ]-Pi;Pi[.
 770  */
 771 bool EdgeArcCircle::IsIn2Pi(double start, double delta, double angleIn)
 772 {
 773   double myDelta=angleIn-start;
 774   if(delta>0.)
 775     {
 776       myDelta=myDelta>=0.?myDelta:myDelta+2.*M_PI;
 777       return myDelta>0. && myDelta<delta;
 778     }
 779   else
 780     {
 781       myDelta=myDelta<=0.?myDelta:myDelta-2.*M_PI;
 782       return myDelta<0. && myDelta>delta;
 783     }
 784 }
 785
 786 /*!
 787  * Given the arc 'a' defined by 'start' angle and a 'delta' [-Pi;Pi] states for the angle 'angleIn' [-Pi;Pi] if it owns or not 'a'.
 788  */
 789 bool EdgeArcCircle::IsAngleNotIn(double start, double delta, double angleIn)
 790 {
 791   double tmp=start;
 792   if(tmp<0.)
 793     tmp+=2*M_PI;
 794   double tmp2=angleIn;
 795   if(tmp2<0.)
 796     tmp2+=2*M_PI;
 797   if(tmp+delta>=2.*M_PI)
 798     return (tmp2<tmp) && (tmp2>tmp+delta-2*M_PI);
 799   else if(tmp+delta>=0.)
 800     return (tmp2<std::min(tmp,tmp+delta) || tmp2>std::max(tmp,tmp+delta));
 801   else
 802     return (tmp2>tmp) && (tmp2<(tmp+delta+2.*M_PI));
 803 }
 804
 805 void EdgeArcCircle::updateBounds()
 806 {
 807   _bounds.setValues(std::min((*_start)[0],(*_end)[0]),std::max((*_start)[0],(*_end)[0]),std::min((*_start)[1],(*_end)[1]),std::max((*_start)[1],(*_end)[1]));
 808   if(IsIn2Pi(_angle0,_angle,M_PI/2))
 809     _bounds[3]=_center[1]+_radius;
 810   if(IsIn2Pi(_angle0,_angle,-M_PI/2))
 811     _bounds[2]=_center[1]-_radius;
 812   if(IsIn2Pi(_angle0,_angle,0.))
 813     _bounds[1]=_center[0]+_radius;
 814   if(IsIn2Pi(_angle0,_angle,M_PI))
 815     _bounds[0]=_center[0]-_radius;
 816 }
 817
 818 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,
 819                                       std::vector<int>& edgesThis, std::vector<double>& addCoo, std::map<INTERP_KERNEL::Node *,int> mapAddCoo) const
 820 {
 821   int tmp[2];
 822   _start->fillGlobalInfoAbs(mapThis,mapOther,offset1,offset2,fact,baryX,baryY,addCoo,mapAddCoo,tmp);
 823   _end->fillGlobalInfoAbs(mapThis,mapOther,offset1,offset2,fact,baryX,baryY,addCoo,mapAddCoo,tmp+1);
 824   if(direction)
 825     {
 826       edgesThis.push_back(tmp[0]);
 827       edgesThis.push_back(tmp[1]);
 828     }
 829   else
 830     {
 831       edgesThis.push_back(tmp[1]);
 832       edgesThis.push_back(tmp[0]);
 833     }
 834 }
 835
 836 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,
 837                                        std::vector<int>& edgesOther, std::vector<double>& addCoo, std::map<INTERP_KERNEL::Node *,int>& mapAddCoo) const
 838 {
 839   _start->fillGlobalInfoAbs2(mapThis,mapOther,offset1,offset2,fact,baryX,baryY,addCoo,mapAddCoo,edgesOther);
 840   _end->fillGlobalInfoAbs2(mapThis,mapOther,offset1,offset2,fact,baryX,baryY,addCoo,mapAddCoo,edgesOther);
 841 }