1 // Copyright (C) 2007-2024 CEA, EDF
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.
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.
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
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
19 // Author : Anthony Geay (CEA/DEN)
21 #include "InterpKernelGeo2DEdgeArcCircle.hxx"
22 #include "InterpKernelGeo2DEdgeLin.hxx"
23 #include "InterpKernelException.hxx"
24 #include "InterpKernelGeo2DNode.hxx"
25 #include "NormalizedUnstructuredMesh.hxx"
31 using namespace INTERP_KERNEL;
33 ArcCArcCIntersector::ArcCArcCIntersector(const EdgeArcCircle& e1, const EdgeArcCircle& e2):SameTypeEdgeIntersector(e1,e2),_dist(0.)
37 bool ArcCArcCIntersector::haveTheySameDirection() const
39 return (getE1().getAngle()>0. && getE2().getAngle()>0.) || (getE1().getAngle()<0. && getE2().getAngle()<0.);
42 bool ArcCArcCIntersector::areColinears() const
44 double radiusL,radiusB;
45 double centerL[2],centerB[2];
47 return internalAreColinears(getE1(),getE2(),tmp,cst,radiusL,centerL,radiusB,centerB);
51 * Precondition 'start' and 'end' are on the same curve than this.
53 void ArcCArcCIntersector::getPlacements(Node *start, Node *end, TypeOfLocInEdge& whereStart, TypeOfLocInEdge& whereEnd, MergePoints& commonNode) const
55 bool obvious1,obvious2;
56 obviousCaseForCurvAbscisse(start,whereStart,commonNode,obvious1);
57 obviousCaseForCurvAbscisse(end,whereEnd,commonNode,obvious2);
58 if(obvious1 && obvious2)
60 double angleInRadStart=getAngle(start);
61 double angleInRadEnd=getAngle(end);
62 if(obvious1 || obvious2)
66 if(EdgeArcCircle::IsIn2Pi(getE1().getAngle0(),getE1().getAngle(),angleInRadEnd))
74 if(EdgeArcCircle::IsIn2Pi(getE1().getAngle0(),getE1().getAngle(),angleInRadStart))
77 whereStart=OUT_BEFORE;
81 if(EdgeArcCircle::IsIn2Pi(getE1().getAngle0(),getE1().getAngle(),angleInRadStart))
84 if(EdgeArcCircle::IsIn2Pi(getE1().getAngle0(),getE1().getAngle(),angleInRadEnd))
90 {//we are out in start.
91 if(EdgeArcCircle::IsIn2Pi(getE1().getAngle0(),getE1().getAngle(),angleInRadEnd))
93 whereStart=OUT_BEFORE;
98 if(EdgeArcCircle::IsIn2Pi(getE2().getAngle0(),getE2().getAngle(),getE1().getAngle0()))
99 {//_e2 contains stictly _e1
100 whereStart=OUT_BEFORE;
104 {//_e2 is outside from _e1
105 whereStart=OUT_BEFORE;
113 * Return angle between ]-Pi;Pi[
115 double ArcCArcCIntersector::getAngle(Node *node) const
117 return EdgeArcCircle::GetAbsoluteAngleOfNormalizedVect(((*node)[0]-getE1().getCenter()[0])/getE1().getRadius(),((*node)[1]-getE1().getCenter()[1])/getE1().getRadius());
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])
123 double lgth1=fabs(a1.getAngle()*a1.getRadius());
124 double lgth2=fabs(a2.getAngle()*a2.getRadius());
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();
132 a2.getCenter(centerL); radiusL=a2.getRadius();
133 a1.getCenter(centerB); radiusB=a1.getRadius();
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.);
142 bool ArcCArcCIntersector::areArcsOverlapped(const EdgeArcCircle& a1, const EdgeArcCircle& a2)
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))
150 double angle0L,angleL;
151 Bounds *merge=a1.getBounds().nearlyAmIIntersectingWith(a2.getBounds());
152 merge->getInterceptedArc(centerL,radiusL,angle0L,angleL);
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.);
174 void ArcCArcCIntersector::areOverlappedOrOnlyColinears(bool& obviousNoIntersection, bool& areOverlapped)
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))
180 obviousNoIntersection=true;
184 if(areArcsOverlapped(getE1(),getE2()))//(Node::areDoubleEquals(_dist,0.) && Node::areDoubleEquals(radius1,radius2))
186 obviousNoIntersection=false;
191 obviousNoIntersection=false;
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.
201 std::list< IntersectElement > ArcCArcCIntersector::getIntersectionsCharacteristicVal() const
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))
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);
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);
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()));
245 //tangent intersection
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()));
261 double signDeltaAngle2;
269 angle0_2=angle0_1+M_PI;
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)) )
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
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);
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()));
313 //tangent intersection
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()));
324 ArcCSegIntersector::ArcCSegIntersector(const EdgeArcCircle& e1, const EdgeLin& e2, bool reverse):
325 CrossTypeEdgeIntersector(e1,e2,reverse),
326 _deltaRoot_div_dr(0.),
327 _i1S2E(false),_i1E2E(false)
329 const double *center=getE1().getCenter();
330 _dx=(*(_e2.getEndNode()))[0]-(*(_e2.getStartNode()))[0];
331 _dy=(*(_e2.getEndNode()))[1]-(*(_e2.getStartNode()))[1];
332 _drSq=_dx*_dx+_dy*_dy;
334 ((*(_e2.getStartNode()))[0]-center[0])*((*(_e2.getEndNode()))[1]-center[1])-
335 ((*(_e2.getStartNode()))[1]-center[1])*((*(_e2.getEndNode()))[0]-center[0]);
339 See http://mathworld.wolfram.com/Circle-LineIntersection.html
340 _cross is 'D', the computation is done with the translation to put back the circle at the origin
342 void ArcCSegIntersector::areOverlappedOrOnlyColinears(bool& obviousNoIntersection, bool& areOverlapped)
344 areOverlapped=false;//No overlapping by construction
346 // Similar optimisation than SegSegIntersector::areOverlappedOrOnlyColinears()
348 identifyEarlyIntersection(dnu1, dnu2, _i1S2E, _i1E2E);
350 const double R = getE1().getRadius();
352 // We need to compute d = R*R-_cross*_cross/_drSq
353 // In terms of numerical precision, this can trigger 'catastrophic cancellation' and is hence better expressed as:
354 double _dr = sqrt(_drSq);
355 double diff = (R-_cross/_dr), add=(R+_cross/_dr);
356 // 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
357 // to prevent ourselves going below machine precision (typ. 1.0e-16 for double).
358 const double eps_machine = std::numeric_limits<double>::epsilon();
359 diff = fabs(diff/R) < eps_machine ? 0.0 : diff;
360 add = fabs(add/R) < eps_machine ? 0.0 : add;
362 // Compute deltaRoot_div_dr := sqrt(delta)/dr, where delta has the meaning of Wolfram.
363 // 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.
364 // We compute it in such a way that it can be used in boolean tests too (a very negative value means we're far apart from intersection)
365 _deltaRoot_div_dr = Node::sign(d)*sqrt(fabs(d));
367 if( 2*_deltaRoot_div_dr > -QuadraticPlanarPrecision::getPrecision())
368 obviousNoIntersection=false;
370 obviousNoIntersection=true;
374 * By construction, no chance that an arc of circle and line to be colinear.
376 bool ArcCSegIntersector::areColinears() const
381 void ArcCSegIntersector::getPlacements(Node *start, Node *end, TypeOfLocInEdge& whereStart, TypeOfLocInEdge& whereEnd, MergePoints& commonNode) const
383 throw Exception("Internal error. Should never been called : no overlapping possible between arc of circle and a segment.");
386 std::list< IntersectElement > ArcCSegIntersector::getIntersectionsCharacteristicVal() const
388 std::list< IntersectElement > ret;
389 const double *center=getE1().getCenter();
390 if(!(2*fabs(_deltaRoot_div_dr) < QuadraticPlanarPrecision::getPrecision())) // see comments in areOverlappedOrOnlyColinears()
391 { // Two intersection nodes
392 // -> if a common node found, there is a chance that this is the only one (i.e. second intersection point is outside e1 and e2)
395 // Check tangent vector of the arc circle at the common node with the linear segment.
396 // There we can tell if the arc of circle is 'moving away' from the seg, or if it might intersect it twice
397 const Node &n(*_earlyInter->getNodeOnly());
399 double tang[2] = {-(n[1]-center[1]), n[0]-center[0]}; // (-y, x) is the tangent vector in the trigo direction with (x,y) = (center->node)
400 bool invSeg = _i1S2E || _i1E2E;
401 double linEdge[2] = {invSeg ? (-_dx) : _dx, invSeg ? (-_dy) : _dy};
402 if(tang[1]*linEdge[0]-tang[0]*linEdge[1] < 0)
404 ret.push_back(*_earlyInter);
409 double determinant=fabs(_deltaRoot_div_dr)/sqrt(_drSq);
410 double x1=(_cross*_dy/_drSq+Node::sign(_dy)*_dx*determinant)+center[0];
411 double y1=(-_cross*_dx/_drSq+fabs(_dy)*determinant)+center[1];
412 Node *intersect1=new Node(x1,y1); intersect1->declareOn();
413 double x2=(_cross*_dy/_drSq-Node::sign(_dy)*_dx*determinant)+center[0];
414 double y2=(-_cross*_dx/_drSq-fabs(_dy)*determinant)+center[1];
415 Node *intersect2=new Node(x2,y2); intersect2->declareOn();
417 bool isN1(false), isN2(false);
420 // Which node do we actually already found? Assume this is the closest ...
421 const Node &iN = *(_earlyInter->getNodeOnly());
422 const Node &n1(*intersect1), &n2(*intersect2);
423 double d1 = std::max(fabs(iN[0]-n1[0]), fabs(iN[1]-n1[1]));
424 double d2 = std::max(fabs(iN[0]-n2[0]), fabs(iN[1]-n2[1]));
425 isN1 = d1 < d2; isN2 = !isN1;
426 if (isN1) intersect1->decrRef();
427 if (isN2) intersect2->decrRef();
428 ret.push_back(*_earlyInter);
432 bool i1_1S=_e1.getStartNode()->isEqual(*intersect1);
433 bool i1_1E=_e1.getEndNode()->isEqual(*intersect1);
434 bool i1_2S=_e2.getStartNode()->isEqual(*intersect1);
435 bool i1_2E=_e2.getEndNode()->isEqual(*intersect1);
436 ret.push_back(IntersectElement(getE1().getCharactValue(*intersect1),getE2().getCharactValue(*intersect1),i1_1S,i1_1E,i1_2S,i1_2E,intersect1,_e1,_e2,keepOrder()));
440 bool i2_1S=_e1.getStartNode()->isEqual(*intersect2);
441 bool i2_1E=_e1.getEndNode()->isEqual(*intersect2);
442 bool i2_2S=_e2.getStartNode()->isEqual(*intersect2);
443 bool i2_2E=_e2.getEndNode()->isEqual(*intersect2);
444 ret.push_back(IntersectElement(getE1().getCharactValue(*intersect2),getE2().getCharactValue(*intersect2),i2_1S,i2_1E,i2_2S,i2_2E,intersect2,_e1,_e2,keepOrder()));
447 else//tangent intersection
451 ret.push_back(*_earlyInter);
454 double x=(_cross*_dy)/_drSq+center[0];
455 double y=(-_cross*_dx)/_drSq+center[1];
456 Node *intersect3=new Node(x,y); intersect3->declareOnTangent();
457 bool i_1S=_e1.getStartNode()->isEqual(*intersect3);
458 bool i_1E=_e1.getEndNode()->isEqual(*intersect3);
459 bool i_2S=_e2.getStartNode()->isEqual(*intersect3);
460 bool i_2E=_e2.getEndNode()->isEqual(*intersect3);
461 ret.push_back(IntersectElement(_e1.getCharactValue(*intersect3),_e2.getCharactValue(*intersect3),i_1S,i_1E,i_2S,i_2E,intersect3,_e1,_e2,keepOrder()));
466 EdgeArcCircle::EdgeArcCircle(std::istream& lineInXfig)
468 const unsigned NB_OF_SKIP_FIELDS=15;
470 for(unsigned i=0;i<NB_OF_SKIP_FIELDS;i++)
472 _start=new Node(lineInXfig);
473 Node *middle=new Node(lineInXfig);
474 _end=new Node(lineInXfig);
475 GetArcOfCirclePassingThru(*_start,*middle,*_end,_center,_radius,_angle,_angle0);
480 EdgeArcCircle::EdgeArcCircle(Node *start, Node *middle, Node *end, bool direction):Edge(start,end, direction)
482 GetArcOfCirclePassingThru(*_start,*middle,*_end,_center,_radius,_angle,_angle0);
486 EdgeArcCircle::EdgeArcCircle(double sX, double sY, double mX, double mY, double eX, double eY):Edge(sX,sY,eX,eY)
488 double middle[2]; middle[0]=mX; middle[1]=mY;
489 GetArcOfCirclePassingThru(*_start,middle,*_end,_center,_radius,_angle,_angle0);
494 * @param angle0 in ]-Pi;Pi[
495 * @param deltaAngle in ]-2.*Pi;2.*Pi[
497 EdgeArcCircle::EdgeArcCircle(Node *start, Node *end, const double *center, double radius, double angle0, double deltaAngle, bool direction):Edge(start,end,direction),_angle0(angle0),_angle(deltaAngle),
500 _center[0]=center[0];
501 _center[1]=center[1];
505 void EdgeArcCircle::changeMiddle(Node *newMiddle)
507 GetArcOfCirclePassingThru(*_start,*newMiddle,*_end,_center,_radius,_angle,_angle0);
511 Edge *EdgeArcCircle::buildEdgeLyingOnMe(Node *start, Node *end, bool direction) const
513 double sx=((*start)[0]-_center[0])/_radius;
514 double sy=((*start)[1]-_center[1])/_radius;
515 double ex=((*end)[0]-_center[0])/_radius;
516 double ey=((*end)[1]-_center[1])/_radius;
517 double angle0=GetAbsoluteAngleOfNormalizedVect(direction?sx:ex,direction?sy:ey);
518 double deltaAngle=GetAbsoluteAngleOfNormalizedVect(sx*ex+sy*ey,sx*ey-sy*ex);
519 if(deltaAngle>0. && _angle<0.)
521 else if(deltaAngle<0. && _angle>0.)
523 deltaAngle=direction?deltaAngle:-deltaAngle;
524 return new EdgeArcCircle(start,end,_center,_radius,angle0,deltaAngle,direction);
527 void EdgeArcCircle::applySimilarity(double xBary, double yBary, double dimChar)
529 Edge::applySimilarity(xBary,yBary,dimChar);
531 _center[0]=(_center[0]-xBary)/dimChar;
532 _center[1]=(_center[1]-yBary)/dimChar;
535 void EdgeArcCircle::unApplySimilarity(double xBary, double yBary, double dimChar)
537 Edge::unApplySimilarity(xBary,yBary,dimChar);
539 _center[0]=_center[0]*dimChar+xBary;
540 _center[1]=_center[1]*dimChar+yBary;
544 * 'eps' is expected to be > 0.
545 * 'conn' is of size 3. conn[0] is start id, conn[1] is end id and conn[2] is middle id.
546 * '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.
548 void EdgeArcCircle::tesselate(const mcIdType *conn, mcIdType offset, double eps, std::vector<mcIdType>& newConn, std::vector<double>& addCoo) const
550 newConn.push_back(INTERP_KERNEL::NORM_POLYL);
551 int nbOfSubDiv=(int)(fabs(_angle)/eps);
554 newConn.push_back(conn[0]); newConn.push_back(conn[2]); newConn.push_back(conn[1]);
557 double signOfAngle=_angle>0.?1.:-1.;
558 mcIdType offset2=offset+(ToIdType(addCoo.size()))/2;
559 newConn.push_back(conn[0]);
560 for(int i=1;i<nbOfSubDiv;i++,offset2++)
562 double angle=_angle0+i*eps*signOfAngle;
563 newConn.push_back(offset2);
564 addCoo.push_back(_center[0]+_radius*cos(angle)); addCoo.push_back(_center[1]+_radius*sin(angle));
566 newConn.push_back(conn[1]);
569 EdgeArcCircle *EdgeArcCircle::BuildFromNodes(Node *start, Node *middle, Node *end)
572 e1=new EdgeLin(start,middle);
573 e2=new EdgeLin(middle,end);
574 SegSegIntersector inters(*e1,*e2);
575 bool colinearity=inters.areColinears();
576 delete e1; delete e2;
579 start->decrRef(); middle->decrRef(); end->decrRef();
584 EdgeArcCircle *ret=new EdgeArcCircle(start,middle,end);
585 start->decrRef(); middle->decrRef(); end->decrRef();
591 * Given an \b NON normalized vector 'vect', returns its norm 'normVect' and its
592 * angle in ]-Pi,Pi] relative to Ox axe.
594 double EdgeArcCircle::GetAbsoluteAngle(const double *vect, double& normVect)
596 normVect=Node::norm(vect);
597 return GetAbsoluteAngleOfNormalizedVect(vect[0]/normVect,vect[1]/normVect);
601 * Given a \b normalized vector defined by (ux,uy) returns its angle in ]-Pi;Pi].
602 * Actually in the current implementation, the vector does not even need to be normalized ...
604 double EdgeArcCircle::GetAbsoluteAngleOfNormalizedVect(double ux, double uy)
606 return atan2(uy, ux);
609 void EdgeArcCircle::GetArcOfCirclePassingThru(const double *start, const double *middle, const double *end,
610 double *center, double& radius, double& angleInRad, double& angleInRad0)
612 double delta=(middle[0]-start[0])*(end[1]-middle[1])-(end[0]-middle[0])*(middle[1]-start[1]);
613 double b1=(middle[1]*middle[1]+middle[0]*middle[0]-start[0]*start[0]-start[1]*start[1])/2;
614 double b2=(end[1]*end[1]+end[0]*end[0]-middle[0]*middle[0]-middle[1]*middle[1])/2;
615 center[0]=((end[1]-middle[1])*b1+(start[1]-middle[1])*b2)/delta;
616 center[1]=((middle[0]-end[0])*b1+(middle[0]-start[0])*b2)/delta;
617 radius=SafeSqrt((start[0]-center[0])*(start[0]-center[0])+(start[1]-center[1])*(start[1]-center[1]));
618 angleInRad0=GetAbsoluteAngleOfNormalizedVect((start[0]-center[0])/radius,(start[1]-center[1])/radius);
619 double angleInRadM=GetAbsoluteAngleOfNormalizedVect((middle[0]-center[0])/radius,(middle[1]-center[1])/radius);
620 angleInRad=GetAbsoluteAngleOfNormalizedVect(((start[0]-center[0])*(end[0]-center[0])+(start[1]-center[1])*(end[1]-center[1]))/(radius*radius),
621 ((start[0]-center[0])*(end[1]-center[1])-(start[1]-center[1])*(end[0]-center[0]))/(radius*radius));
622 if(IsAngleNotIn(angleInRad0,angleInRad,angleInRadM))
623 angleInRad=angleInRad<0?2*M_PI+angleInRad:angleInRad-2*M_PI;
626 void EdgeArcCircle::dumpInXfigFile(std::ostream& stream, bool direction, int resolution, const Bounds& box) const
628 stream << "5 1 0 1 ";
629 fillXfigStreamForLoc(stream);
630 stream << " 7 50 -1 -1 0.000 0 ";
631 if( (direction && (-_angle)>=0) || (!direction && (-_angle)<0))
636 stream << box.fitXForXFigD(_center[0],resolution) << " " << box.fitYForXFigD(_center[1],resolution) << " ";
637 direction?_start->dumpInXfigFile(stream,resolution,box):_end->dumpInXfigFile(stream,resolution,box);
638 Node *middle=buildRepresentantOfMySelf();
639 middle->dumpInXfigFile(stream,resolution,box);
641 direction?_end->dumpInXfigFile(stream,resolution,box):_start->dumpInXfigFile(stream,resolution,box);
642 stream << std::endl << "1 1 2.00 120.00 180.00" << std::endl;
645 void EdgeArcCircle::update(Node *m)
647 GetArcOfCirclePassingThru(*_start,*m,*_end,_center,_radius,_angle,_angle0);
652 * This methods computes :
654 * \int_{Current Edge} -ydx
657 double EdgeArcCircle::getAreaOfZone() const
659 return -_radius*_radius*(sin(_angle)-_angle)/2.+((*_start)[0]-(*_end)[0])*((*_start)[1]+(*_end)[1])/2.;
662 double EdgeArcCircle::getCurveLength() const
664 return fabs(_angle*_radius);
667 void EdgeArcCircle::getBarycenter(double *bary) const
669 bary[0]=_center[0]+_radius*cos(_angle0+_angle/2.);
670 bary[1]=_center[1]+_radius*sin(_angle0+_angle/2.);
675 * bary[0]=\int_{Current Edge} -yxdx
678 * bary[1]=\int_{Current Edge} -\frac{y^{2}}{2}dx
680 * To compute these 2 expressions in this class we have :
682 * x=x_{0}+Radius \cdot cos(\theta)
685 * y=y_{0}+Radius \cdot sin(\theta)
688 * dx=-Radius \cdot sin(\theta) \cdot d\theta
691 void EdgeArcCircle::getBarycenterOfZone(double *bary) const
693 double x0=_center[0];
694 double y0=_center[1];
695 double angle1=_angle0+_angle;
696 double tmp1=sin(angle1);
697 double tmp0=sin(_angle0);
698 double tmp2=_radius*_radius*_radius;
699 double tmp3=cos(angle1);
700 double tmp4=cos(_angle0);
701 bary[0]=_radius*x0*y0*(tmp4-tmp3)+_radius*_radius*(y0*(cos(2*_angle0)-cos(2*angle1))/4.+
702 x0*(_angle/2.+(sin(2.*_angle0)-sin(2.*angle1))/4.))
703 +tmp2*(tmp1*tmp1*tmp1-tmp0*tmp0*tmp0)/3.;
704 bary[1]=y0*y0*_radius*(tmp4-tmp3)/2.+_radius*_radius*y0*(_angle/2.+(sin(2.*_angle0)-sin(2.*angle1))/4.)
705 +tmp2*(tmp4-tmp3+(tmp3*tmp3*tmp3-tmp4*tmp4*tmp4)/3.)/2.;
709 * Compute the "middle" of two points on the arc of circle.
710 * The order (p1,p2) or (p2,p1) doesn't matter. p1 and p2 have to be localized on the edge defined by this.
711 * \param[out] mid the point located half-way between p1 and p2 on the arc defined by this.
712 * \sa getMiddleOfPointsOriented() a generalisation working also when p1 and p2 are not on the arc.
714 void EdgeArcCircle::getMiddleOfPoints(const double *p1, const double *p2, double *mid) const
716 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);
717 double angle1(GetAbsoluteAngleOfNormalizedVect(dx1,dy1)),angle2(GetAbsoluteAngleOfNormalizedVect(dx2,dy2));
719 double myDelta1(angle1-_angle0),myDelta2(angle2-_angle0);
721 { myDelta1=myDelta1>-QuadraticPlanarPrecision::getPrecision()?myDelta1:myDelta1+2.*M_PI; myDelta2=myDelta2>-QuadraticPlanarPrecision::getPrecision()?myDelta2:myDelta2+2.*M_PI; }
723 { myDelta1=myDelta1<QuadraticPlanarPrecision::getPrecision()?myDelta1:myDelta1-2.*M_PI; myDelta2=myDelta2<QuadraticPlanarPrecision::getPrecision()?myDelta2:myDelta2-2.*M_PI; }
725 mid[0]=_center[0]+_radius*cos(_angle0+(myDelta1+myDelta2)/2.);
726 mid[1]=_center[1]+_radius*sin(_angle0+(myDelta1+myDelta2)/2.);
730 * Compute the "middle" of two points on the arc of circle.
731 * Walk on the circle from p1 to p2 using the rotation direction indicated by this->_angle (i.e. by the orientation of the arc).
732 * This function is sensitive to the ordering of p1 and p2.
733 * \param[out] mid the point located half-way between p1 and p2
734 * \sa getMiddleOfPoints() to be used when the order of p1 and p2 is not relevant.
736 void EdgeArcCircle::getMiddleOfPointsOriented(const double *p1, const double *p2, double *mid) const
738 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);
739 double angle1(GetAbsoluteAngleOfNormalizedVect(dx1,dy1)),angle2(GetAbsoluteAngleOfNormalizedVect(dx2,dy2));
747 if((_angle>0. && angle1 <= angle2) || (_angle<=0. && angle1 >= angle2))
748 avg = (angle1+angle2)/2.;
750 avg = (angle1+angle2)/2. - M_PI;
752 mid[0]=_center[0]+_radius*cos(avg);
753 mid[1]=_center[1]+_radius*sin(avg);
758 * Characteristic value used is angle in ]_Pi;Pi[ from axe 0x.
760 bool EdgeArcCircle::isIn(double characterVal) const
762 return IsIn2Pi(_angle0,_angle,characterVal);
765 Node *EdgeArcCircle::buildRepresentantOfMySelf() const
767 return new Node(_center[0]+_radius*cos(_angle0+_angle/2.),_center[1]+_radius*sin(_angle0+_angle/2.));
771 * Characteristic value used is angle in ]_Pi;Pi[ from axe 0x.
772 * 'val1' and 'val2' have been detected previously as owning to this.
774 bool EdgeArcCircle::isLower(double val1, double val2) const
776 double myDelta1=val1-_angle0;
777 double myDelta2=val2-_angle0;
780 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.
781 myDelta2=myDelta2>-(_radius*QuadraticPlanarPrecision::getPrecision())?myDelta2:myDelta2+2.*M_PI;
782 return myDelta1<myDelta2;
786 myDelta1=myDelta1<(_radius*QuadraticPlanarPrecision::getPrecision())?myDelta1:myDelta1-2.*M_PI;
787 myDelta2=myDelta2<(_radius*QuadraticPlanarPrecision::getPrecision())?myDelta2:myDelta2-2.*M_PI;
788 return myDelta2<myDelta1;
793 * For Arc circle the caract value is angle with Ox between -Pi and Pi.
795 double EdgeArcCircle::getCharactValue(const Node& node) const
797 double dx=(node[0]-_center[0])/_radius;
798 double dy=(node[1]-_center[1])/_radius;
799 return GetAbsoluteAngleOfNormalizedVect(dx,dy);
802 double EdgeArcCircle::getCharactValueBtw0And1(const Node& node) const
804 double dx=(node[0]-_center[0])/_radius;
805 double dy=(node[1]-_center[1])/_radius;
806 double angle=GetAbsoluteAngleOfNormalizedVect(dx,dy);
808 double myDelta=angle-_angle0;
810 myDelta=myDelta>=0.?myDelta:myDelta+2.*M_PI;
812 myDelta=myDelta<=0.?myDelta:myDelta-2.*M_PI;
813 return myDelta/_angle;
816 double EdgeArcCircle::getDistanceToPoint(const double *pt) const
818 double angle=Node::computeAngle(_center,pt);
819 if(IsIn2Pi(_angle0,_angle,angle))
820 return fabs(Node::distanceBtw2Pt(_center,pt)-_radius);
823 double dist1=Node::distanceBtw2Pt(*_start,pt);
824 double dist2=Node::distanceBtw2Pt(*_end,pt);
825 return std::min(dist1,dist2);
829 bool EdgeArcCircle::isNodeLyingOn(const double *coordOfNode) const
831 double dist=Node::distanceBtw2Pt(_center,coordOfNode);
832 if(Node::areDoubleEquals(dist,_radius))
834 double angle=Node::computeAngle(_center,coordOfNode);
835 return IsIn2Pi(_angle0,_angle,angle);
842 * Idem IsAngleNotIn except that here 'start' in ]-Pi;Pi[ and delta in ]-2*Pi;2Pi[.
843 * @param angleIn in ]-Pi;Pi[.
845 bool EdgeArcCircle::IsIn2Pi(double start, double delta, double angleIn)
847 double myDelta=angleIn-start;
850 myDelta=myDelta>=0.?myDelta:myDelta+2.*M_PI;
851 return myDelta>0. && myDelta<delta;
855 myDelta=myDelta<=0.?myDelta:myDelta-2.*M_PI;
856 return myDelta<0. && myDelta>delta;
861 * 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'.
863 bool EdgeArcCircle::IsAngleNotIn(double start, double delta, double angleIn)
871 if(tmp+delta>=2.*M_PI)
872 return (tmp2<tmp) && (tmp2>tmp+delta-2*M_PI);
873 else if(tmp+delta>=0.)
874 return (tmp2<std::min(tmp,tmp+delta) || tmp2>std::max(tmp,tmp+delta));
876 return (tmp2>tmp) && (tmp2<(tmp+delta+2.*M_PI));
879 void EdgeArcCircle::updateBounds()
881 _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]));
882 if(IsIn2Pi(_angle0,_angle,M_PI/2))
883 _bounds[3]=_center[1]+_radius;
884 if(IsIn2Pi(_angle0,_angle,-M_PI/2))
885 _bounds[2]=_center[1]-_radius;
886 if(IsIn2Pi(_angle0,_angle,0.))
887 _bounds[1]=_center[0]+_radius;
888 if(IsIn2Pi(_angle0,_angle,M_PI))
889 _bounds[0]=_center[0]-_radius;