src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.cxx

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