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