Back to index

salome-med  6.5.0
InterpKernelGeo2DEdgeArcCircle.cxx
Go to the documentation of this file.
00001 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D
00002 //
00003 // This library is free software; you can redistribute it and/or
00004 // modify it under the terms of the GNU Lesser General Public
00005 // License as published by the Free Software Foundation; either
00006 // version 2.1 of the License.
00007 //
00008 // This library is distributed in the hope that it will be useful,
00009 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00010 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00011 // Lesser General Public License for more details.
00012 //
00013 // You should have received a copy of the GNU Lesser General Public
00014 // License along with this library; if not, write to the Free Software
00015 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
00016 //
00017 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
00018 //
00019 
00020 #include "InterpKernelGeo2DEdgeArcCircle.hxx"
00021 #include "InterpKernelGeo2DEdgeLin.hxx"
00022 #include "InterpKernelException.hxx"
00023 #include "InterpKernelGeo2DNode.hxx"
00024 #include "NormalizedUnstructuredMesh.hxx"
00025 
00026 #include <sstream>
00027 #include <algorithm>
00028 
00029 using namespace INTERP_KERNEL;
00030 
00031 ArcCArcCIntersector::ArcCArcCIntersector(const EdgeArcCircle& e1, const EdgeArcCircle& e2):SameTypeEdgeIntersector(e1,e2),_dist(0.)
00032 {
00033 }
00034 
00035 bool ArcCArcCIntersector::haveTheySameDirection() const
00036 {
00037   return (getE1().getAngle()>0. &&  getE2().getAngle()>0.) || (getE1().getAngle()<0. &&  getE2().getAngle()<0.);
00038 }
00039 
00043 void ArcCArcCIntersector::getPlacements(Node *start, Node *end, TypeOfLocInEdge& whereStart, TypeOfLocInEdge& whereEnd, MergePoints& commonNode) const
00044 {
00045   bool obvious1,obvious2;
00046   obviousCaseForCurvAbscisse(start,whereStart,commonNode,obvious1);
00047   obviousCaseForCurvAbscisse(end,whereEnd,commonNode,obvious2);
00048   if(obvious1 && obvious2)
00049     return ;
00050   double angleInRadStart=getAngle(start);
00051   double angleInRadEnd=getAngle(end);
00052   if(obvious1 || obvious2)
00053     {
00054       if(obvious1)
00055         {
00056           if(EdgeArcCircle::IsIn2Pi(getE1().getAngle0(),getE1().getAngle(),angleInRadEnd))
00057             whereEnd=INSIDE;
00058           else
00059             whereEnd=OUT_AFTER;
00060           return ;
00061         }
00062       else
00063         {
00064           if(EdgeArcCircle::IsIn2Pi(getE1().getAngle0(),getE1().getAngle(),angleInRadStart))
00065             whereStart=INSIDE;
00066           else
00067             whereStart=OUT_BEFORE;
00068           return ;
00069         }
00070     }
00071   if(EdgeArcCircle::IsIn2Pi(getE1().getAngle0(),getE1().getAngle(),angleInRadStart))
00072     {
00073       whereStart=INSIDE;
00074       if(EdgeArcCircle::IsIn2Pi(getE1().getAngle0(),getE1().getAngle(),angleInRadEnd))
00075         whereEnd=INSIDE;
00076       else
00077         whereEnd=OUT_AFTER;
00078     }
00079   else
00080     {//we are out in start.
00081       if(EdgeArcCircle::IsIn2Pi(getE1().getAngle0(),getE1().getAngle(),angleInRadEnd))
00082         {
00083           whereStart=OUT_BEFORE;
00084           whereEnd=INSIDE;
00085         }
00086       else
00087         {
00088           if(EdgeArcCircle::IsIn2Pi(getE2().getAngle0(),getE2().getAngle(),getE1().getAngle0()))
00089             {//_e2 contains stictly _e1
00090               whereStart=OUT_BEFORE;
00091               whereEnd=OUT_AFTER;
00092             }
00093           else
00094             {//_e2 is outside from _e1
00095               whereStart=OUT_BEFORE;
00096               whereEnd=OUT_BEFORE;
00097             }
00098         }
00099     }
00100 }
00101 
00105 double ArcCArcCIntersector::getAngle(Node *node) const
00106 {
00107   return EdgeArcCircle::GetAbsoluteAngleOfNormalizedVect(((*node)[0]-getE1().getCenter()[0])/getE1().getRadius(),((*node)[1]-getE1().getCenter()[1])/getE1().getRadius());
00108 }
00109 
00110 bool ArcCArcCIntersector::areArcsOverlapped(const EdgeArcCircle& a1, const EdgeArcCircle& a2)
00111 {
00112   double centerL[2],radiusL,angle0L,angleL;
00113   double centerB[2],radiusB;
00114   double lgth1=fabs(a1.getAngle()*a1.getRadius());
00115   double lgth2=fabs(a2.getAngle()*a2.getRadius());
00116   if(lgth1<lgth2)
00117     {//a1 is the little one ('L') and a2 the big one ('B')
00118       a1.getCenter(centerL); radiusL=a1.getRadius(); angle0L=a1.getAngle0(); angleL=a1.getAngle();
00119       a2.getCenter(centerB); radiusB=a2.getRadius();
00120     }
00121   else
00122     {
00123       a2.getCenter(centerL); radiusL=a2.getRadius(); angle0L=a2.getAngle0(); angleL=a2.getAngle();
00124       a1.getCenter(centerB); radiusB=a1.getRadius();
00125     }
00126   // dividing from the begining by radiusB^2 to keep precision
00127   double tmp=Node::distanceBtw2PtSq(centerL,centerB);
00128   double cst=tmp/(radiusB*radiusB);
00129   cst+=radiusL*radiusL/(radiusB*radiusB);
00130   if(!Node::areDoubleEqualsWP(cst,1.,2.))
00131     return false;
00132   //
00133   Bounds *merge=a1.getBounds().nearlyAmIIntersectingWith(a2.getBounds());
00134   merge->getInterceptedArc(centerL,radiusL,angle0L,angleL);
00135   delete merge;
00136   //
00137   tmp=sqrt(tmp);
00138   if(Node::areDoubleEqualsWP(tmp,0.,1/(10*std::max(radiusL,radiusB))))
00139     {
00140       if(Node::areDoubleEquals(radiusL,radiusB))
00141         return true;
00142       else
00143         return false;
00144     }
00145   double phi=EdgeArcCircle::GetAbsoluteAngleOfNormalizedVect((centerL[0]-centerB[0])/tmp,(centerL[1]-centerB[1])/tmp);
00146   double cst2=2*radiusL*tmp/(radiusB*radiusB);
00147   double cmpContainer[4];
00148   int sizeOfCmpContainer=2;
00149   cmpContainer[0]=cst+cst2*cos(phi-angle0L);
00150   cmpContainer[1]=cst+cst2*cos(phi-angle0L+angleL);
00151   double a=EdgeArcCircle::NormalizeAngle(phi-angle0L);
00152   if(EdgeArcCircle::IsIn2Pi(angle0L,angleL,a))
00153     cmpContainer[sizeOfCmpContainer++]=cst+cst2;
00154   a=EdgeArcCircle::NormalizeAngle(phi-angle0L+M_PI);
00155   if(EdgeArcCircle::IsIn2Pi(angle0L,angleL,a))
00156     cmpContainer[sizeOfCmpContainer++]=cst-cst2;
00157   a=*std::max_element(cmpContainer,cmpContainer+sizeOfCmpContainer);
00158   return Node::areDoubleEqualsWP(a,1.,2.);
00159 }
00160 
00161 void ArcCArcCIntersector::areOverlappedOrOnlyColinears(const Bounds *whereToFind, bool& obviousNoIntersection, bool& areOverlapped)
00162 {
00163   _dist=Node::distanceBtw2Pt(getE1().getCenter(),getE2().getCenter());
00164   double radius1=getE1().getRadius(); double radius2=getE2().getRadius();
00165   if(_dist>radius1+radius2+QUADRATIC_PLANAR::_precision || _dist+std::min(radius1,radius2)+QUADRATIC_PLANAR::_precision<std::max(radius1,radius2))
00166     {
00167       obviousNoIntersection=true;
00168       areOverlapped=false;
00169       return ;
00170     }
00171   if(areArcsOverlapped(getE1(),getE2()))//(Node::areDoubleEquals(_dist,0.) && Node::areDoubleEquals(radius1,radius2))
00172     {
00173       obviousNoIntersection=false;
00174       areOverlapped=true;
00175     }
00176   else
00177     {
00178       obviousNoIntersection=false;
00179       areOverlapped=false;
00180     }
00181 }
00182 
00183 std::list< IntersectElement > ArcCArcCIntersector::getIntersectionsCharacteristicVal() const
00184 {
00185   std::list< IntersectElement > ret;
00186   const double *center1=getE1().getCenter();
00187   const double *center2=getE2().getCenter();
00188   double radius1=getE1().getRadius(); double radius2=getE2().getRadius();
00189   double d1_1=(_dist*_dist-radius2*radius2+radius1*radius1)/(2.*_dist);
00190   double u[2];//u is normalized vector from center1 to center2.
00191   u[0]=(center2[0]-center1[0])/_dist; u[1]=(center2[1]-center1[1])/_dist;
00192   double d1_1y=EdgeArcCircle::SafeSqrt(radius1*radius1-d1_1*d1_1);
00193   double angleE1=EdgeArcCircle::NormalizeAngle(getE1().getAngle0()+getE1().getAngle());
00194   double angleE2=EdgeArcCircle::NormalizeAngle(getE2().getAngle0()+getE2().getAngle());
00195   if(!Node::areDoubleEquals(d1_1y,0))
00196     {
00197       //2 intersections
00198       double v1[2],v2[2];
00199       v1[0]=u[0]*d1_1-u[1]*d1_1y; v1[1]=u[1]*d1_1+u[0]*d1_1y;
00200       v2[0]=u[0]*d1_1+u[1]*d1_1y; v2[1]=u[1]*d1_1-u[0]*d1_1y;
00201       Node *node1=new Node(center1[0]+v1[0],center1[1]+v1[1]); node1->declareOn();
00202       Node *node2=new Node(center1[0]+v2[0],center1[1]+v2[1]); node2->declareOn();
00203       double angle1_1=EdgeArcCircle::GetAbsoluteAngleOfNormalizedVect(v1[0]/radius1,v1[1]/radius1);
00204       double angle2_1=EdgeArcCircle::GetAbsoluteAngleOfNormalizedVect(v2[0]/radius1,v2[1]/radius1);
00205       double v3[2],v4[2];
00206       v3[0]=center1[0]-center2[0]+v1[0]; v3[1]=center1[1]-center2[1]+v1[1];
00207       v4[0]=center1[0]-center2[0]+v2[0]; v4[1]=center1[1]-center2[1]+v2[1];
00208       double angle1_2=EdgeArcCircle::GetAbsoluteAngleOfNormalizedVect(v3[0]/radius2,v3[1]/radius2);
00209       double angle2_2=EdgeArcCircle::GetAbsoluteAngleOfNormalizedVect(v4[0]/radius2,v4[1]/radius2);
00210       //
00211       bool e1_1S=Node::areDoubleEqualsWP(angle1_1,getE1().getAngle0(),radius1);
00212       bool e1_1E=Node::areDoubleEqualsWP(angle1_1,angleE1,radius1);
00213       bool e1_2S=Node::areDoubleEqualsWP(angle1_2,getE2().getAngle0(),radius1);
00214       bool e1_2E=Node::areDoubleEqualsWP(angle1_2,angleE2,radius1);
00215       //
00216       bool e2_1S=Node::areDoubleEqualsWP(angle2_1,getE1().getAngle0(),radius2);
00217       bool e2_1E=Node::areDoubleEqualsWP(angle2_1,angleE1,radius2);
00218       bool e2_2S=Node::areDoubleEqualsWP(angle2_2,getE2().getAngle0(),radius2);
00219       bool e2_2E=Node::areDoubleEqualsWP(angle2_2,angleE2,radius2);
00220       ret.push_back(IntersectElement(angle1_1,angle1_2,e1_1S,e1_1E,e1_2S,e1_2E,node1,_e1,_e2,keepOrder()));
00221       ret.push_back(IntersectElement(angle2_1,angle2_2,e2_1S,e2_1E,e2_2S,e2_2E,node2,_e1,_e2,keepOrder()));
00222     }
00223   else
00224     {
00225       //tangent intersection
00226       double v1[2],v2[2];
00227       v1[0]=d1_1*u[0]; v1[1]=d1_1*u[1];
00228       v2[0]=center1[0]-center2[0]+v1[0]; v2[1]=center1[1]-center2[1]+v1[1];
00229       double angle0_1=EdgeArcCircle::GetAbsoluteAngleOfNormalizedVect(v1[0]/radius1,v1[1]/radius1);
00230       double angle0_2=EdgeArcCircle::GetAbsoluteAngleOfNormalizedVect(v2[0]/radius2,v2[1]/radius2);
00231       bool e0_1S=Node::areDoubleEqualsWP(angle0_1,getE1().getAngle0(),radius1);
00232       bool e0_1E=Node::areDoubleEqualsWP(angle0_1,angleE1,radius1);
00233       bool e0_2S=Node::areDoubleEqualsWP(angle0_2,getE2().getAngle0(),radius2);
00234       bool e0_2E=Node::areDoubleEqualsWP(angle0_2,angleE2,radius2);
00235       Node *node=new Node(center1[0]+d1_1*u[0],center1[1]+d1_1*u[1]); node->declareOnTangent();
00236       ret.push_back(IntersectElement(angle0_1,angle0_2,e0_1S,e0_1E,e0_2S,e0_2E,node,_e1,_e2,keepOrder()));
00237     }
00238   return ret;
00239 }
00240 /*double angle0_2;
00241   double signDeltaAngle2;
00242   double d1_2;
00243   if(u[1]<0.)
00244   angle0_1=-angle0_1;
00245   if(d1_1>=0.)
00246   {
00247   if(_dist>radius1)
00248   {
00249   angle0_2=angle0_1+M_PI;
00250   signDeltaAngle2=-1.;
00251   }
00252   else
00253   {
00254   angle0_2=angle0_1;
00255   signDeltaAngle2=1.;
00256   }
00257   }
00258   else
00259   {
00260   angle0_1+=M_PI;
00261   angle0_2=angle0_1;
00262   signDeltaAngle2=1.;
00263   }
00264   angle0_1=NormalizeAngle(angle0_1);
00265   angle0_2=NormalizeAngle(angle0_2);
00266   double angleE1=NormalizeAngle(getE1().getAngle0()+getE1().getAngle());
00267   double angleE2=NormalizeAngle(getE2().getAngle0()+getE2().getAngle());
00268   if(!(Node::areDoubleEquals(d1_1,radius1) || Node::areDoubleEquals(d1_1,-radius1)) )
00269   {
00270   //2 intersections   
00271   double deltaAngle1=EdgeArcCircle::SafeAcos(fabs(d1_1)/radius1); //owns to 0;Pi/2 by construction
00272   double deltaAngle2=EdgeArcCircle::SafeAcos(fabs(d1_2)/radius2); //owns to 0;Pi/2 by construction
00273   double angle1_1=NormalizeAngle(angle0_1+deltaAngle1);// Intersection 1 seen for _e1
00274   double angle2_1=NormalizeAngle(angle0_1-deltaAngle1);// Intersection 2 seen for _e1
00275   double angle1_2=NormalizeAngle(angle0_2+signDeltaAngle2*deltaAngle2);// Intersection 1 seen for _e2
00276   double angle2_2=NormalizeAngle(angle0_2-signDeltaAngle2*deltaAngle2);// Intersection 2 seen for _e2
00277   //
00278   bool e1_1S=Node::areDoubleEqualsWP(angle1_1,getE1().getAngle0(),radius1);
00279   bool e1_1E=Node::areDoubleEqualsWP(angle1_1,angleE1,radius1);
00280   bool e1_2S=Node::areDoubleEqualsWP(angle1_2,getE2().getAngle0(),radius1);
00281   bool e1_2E=Node::areDoubleEqualsWP(angle1_2,angleE2,radius1);
00282   //
00283   bool e2_1S=Node::areDoubleEqualsWP(angle2_1,getE1().getAngle0(),radius2);
00284   bool e2_1E=Node::areDoubleEqualsWP(angle2_1,angleE1,radius2);
00285   bool e2_2S=Node::areDoubleEqualsWP(angle2_2,getE2().getAngle0(),radius2);
00286   bool e2_2E=Node::areDoubleEqualsWP(angle2_2,angleE2,radius2);
00287   Node *node1=new Node(center1[0]+radius1*cos(angle1_1),center1[0]+radius1*sin(angle1_1)); node1->declareOn();
00288   Node *node2=new Node(center1[0]+radius1*cos(angle2_1),center1[0]+radius1*sin(angle2_1)); node2->declareOn();
00289   ret.push_back(IntersectElement(angle1_1,angle1_2,e1_1S,e1_1E,e1_2S,e1_2E,node1,_e1,_e2,keepOrder()));
00290   ret.push_back(IntersectElement(angle2_1,angle2_2,e2_1S,e2_1E,e2_2S,e2_2E,node2,_e1,_e2,keepOrder()));
00291   }
00292   else
00293   //tangent intersection
00294   {
00295   bool e0_1S=Node::areDoubleEqualsWP(angle0_1,getE1().getAngle0(),radius1);
00296   bool e0_1E=Node::areDoubleEqualsWP(angle0_1,angleE1,radius1);
00297   bool e0_2S=Node::areDoubleEqualsWP(angle0_2,getE2().getAngle0(),radius2);
00298   bool e0_2E=Node::areDoubleEqualsWP(angle0_2,angleE2,radius2);
00299   Node *node=new Node(center1[0]+radius1*cos(angle0_1),center1[0]+radius1*sin(angle0_1)); node->declareOnTangent();
00300   ret.push_back(IntersectElement(angle0_1,angle0_2,e0_1S,e0_1E,e0_2S,e0_2E,node,_e1,_e2,keepOrder()));
00301   }
00302   return ret;*/
00303 
00304 ArcCSegIntersector::ArcCSegIntersector(const EdgeArcCircle& e1, const EdgeLin& e2, bool reverse):CrossTypeEdgeIntersector(e1,e2,reverse)
00305 {
00306 }
00307 
00308 void ArcCSegIntersector::areOverlappedOrOnlyColinears(const Bounds *whereToFind, bool& obviousNoIntersection, bool& areOverlapped)
00309 {
00310   areOverlapped=false;//No overlapping by contruction
00311   const double *center=getE1().getCenter();
00312   _dx=(*(_e2.getEndNode()))[0]-(*(_e2.getStartNode()))[0];
00313   _dy=(*(_e2.getEndNode()))[1]-(*(_e2.getStartNode()))[1];
00314   _drSq=_dx*_dx+_dy*_dy;
00315   _cross=
00316     ((*(_e2.getStartNode()))[0]-center[0])*((*(_e2.getEndNode()))[1]-center[1])-
00317     ((*(_e2.getStartNode()))[1]-center[1])*((*(_e2.getEndNode()))[0]-center[0]);
00318   _determinant=getE1().getRadius()*getE1().getRadius()/_drSq-_cross*_cross/(_drSq*_drSq);
00319   if(_determinant>-2*QUADRATIC_PLANAR::_precision)//QUADRATIC_PLANAR::_precision*QUADRATIC_PLANAR::_precision*_drSq*_drSq/(2.*_dx*_dx))
00320     obviousNoIntersection=false;
00321   else
00322     obviousNoIntersection=true;   
00323 }
00324 
00325 void ArcCSegIntersector::getPlacements(Node *start, Node *end, TypeOfLocInEdge& whereStart, TypeOfLocInEdge& whereEnd, MergePoints& commonNode) const
00326 {
00327   throw Exception("Internal error. Should never been called : no overlapping possible between arc of circle and a segment.");
00328 }
00329 
00330 std::list< IntersectElement > ArcCSegIntersector::getIntersectionsCharacteristicVal() const
00331 {
00332   std::list< IntersectElement > ret;
00333   const double *center=getE1().getCenter();
00334   if(!(fabs(_determinant)<(2.*QUADRATIC_PLANAR::_precision)))//QUADRATIC_PLANAR::_precision*QUADRATIC_PLANAR::_precision*_drSq*_drSq/(2.*_dx*_dx))
00335     {
00336       double determinant=EdgeArcCircle::SafeSqrt(_determinant);
00337       double x1=(_cross*_dy/_drSq+Node::sign(_dy)*_dx*determinant)+center[0];
00338       double y1=(-_cross*_dx/_drSq+fabs(_dy)*determinant)+center[1];
00339       Node *intersect1=new Node(x1,y1); intersect1->declareOn();
00340       bool i1_1S=_e1.getStartNode()->isEqual(*intersect1);
00341       bool i1_1E=_e1.getEndNode()->isEqual(*intersect1);
00342       bool i1_2S=_e2.getStartNode()->isEqual(*intersect1);
00343       bool i1_2E=_e2.getEndNode()->isEqual(*intersect1);
00344       ret.push_back(IntersectElement(getE1().getCharactValue(*intersect1),getE2().getCharactValue(*intersect1),i1_1S,i1_1E,i1_2S,i1_2E,intersect1,_e1,_e2,keepOrder()));
00345       //
00346       double x2=(_cross*_dy/_drSq-Node::sign(_dy)*_dx*determinant)+center[0];
00347       double y2=(-_cross*_dx/_drSq-fabs(_dy)*determinant)+center[1];
00348       Node *intersect2=new Node(x2,y2); intersect2->declareOn();
00349       bool i2_1S=_e1.getStartNode()->isEqual(*intersect2);
00350       bool i2_1E=_e1.getEndNode()->isEqual(*intersect2);
00351       bool i2_2S=_e2.getStartNode()->isEqual(*intersect2);
00352       bool i2_2E=_e2.getEndNode()->isEqual(*intersect2);
00353       ret.push_back(IntersectElement(getE1().getCharactValue(*intersect2),getE2().getCharactValue(*intersect2),i2_1S,i2_1E,i2_2S,i2_2E,intersect2,_e1,_e2,keepOrder()));
00354     }
00355   else//tangent intersection
00356     {
00357       double x=(_cross*_dy)/_drSq+center[0];
00358       double y=(-_cross*_dx)/_drSq+center[1];
00359       Node *intersect3=new Node(x,y); intersect3->declareOnTangent();
00360       bool i_1S=_e1.getStartNode()->isEqual(*intersect3);
00361       bool i_1E=_e1.getEndNode()->isEqual(*intersect3);
00362       bool i_2S=_e2.getStartNode()->isEqual(*intersect3);
00363       bool i_2E=_e2.getEndNode()->isEqual(*intersect3);
00364       ret.push_back(IntersectElement(_e1.getCharactValue(*intersect3),_e2.getCharactValue(*intersect3),i_1S,i_1E,i_2S,i_2E,intersect3,_e1,_e2,keepOrder()));
00365     }
00366   return ret;
00367 }
00368 
00369 EdgeArcCircle::EdgeArcCircle(std::istream& lineInXfig)
00370 {
00371   const unsigned NB_OF_SKIP_FIELDS=15;
00372   std::string tmpS;
00373   for(unsigned i=0;i<NB_OF_SKIP_FIELDS;i++)
00374     lineInXfig >> tmpS;
00375   _start=new Node(lineInXfig);
00376   Node *middle=new Node(lineInXfig);
00377   _end=new Node(lineInXfig);
00378   GetArcOfCirclePassingThru(*_start,*middle,*_end,_center,_radius,_angle,_angle0);
00379   middle->decrRef();
00380   updateBounds();
00381 }
00382 
00383 EdgeArcCircle::EdgeArcCircle(Node *start, Node *middle, Node *end, bool direction):Edge(start,end, direction)
00384 {
00385   GetArcOfCirclePassingThru(*_start,*middle,*_end,_center,_radius,_angle,_angle0);
00386   updateBounds();
00387 }
00388 
00389 EdgeArcCircle::EdgeArcCircle(double sX, double sY, double mX, double mY, double eX, double eY):Edge(sX,sY,eX,eY)
00390 {
00391   double middle[2]; middle[0]=mX; middle[1]=mY;
00392   GetArcOfCirclePassingThru(*_start,middle,*_end,_center,_radius,_angle,_angle0);
00393   updateBounds();
00394 }
00395 
00400 EdgeArcCircle::EdgeArcCircle(Node *start, Node *end, const double *center, double radius, double angle0, double deltaAngle, bool direction):Edge(start,end,direction),_angle(deltaAngle),
00401                                                                                                                                             _angle0(angle0),_radius(radius)
00402 {
00403   _center[0]=center[0];
00404   _center[1]=center[1];
00405   updateBounds();
00406 }
00407 
00408 void EdgeArcCircle::changeMiddle(Node *newMiddle)
00409 {
00410   GetArcOfCirclePassingThru(*_start,*newMiddle,*_end,_center,_radius,_angle,_angle0);
00411   updateBounds();
00412 }
00413 
00414 Edge *EdgeArcCircle::buildEdgeLyingOnMe(Node *start, Node *end, bool direction) const
00415 {
00416   double sx=((*start)[0]-_center[0])/_radius;
00417   double sy=((*start)[1]-_center[1])/_radius;
00418   double ex=((*end)[0]-_center[0])/_radius;
00419   double ey=((*end)[1]-_center[1])/_radius;
00420   double angle0=GetAbsoluteAngleOfNormalizedVect(direction?sx:ex,direction?sy:ey);
00421   double deltaAngle=GetAbsoluteAngleOfNormalizedVect(sx*ex+sy*ey,sx*ey-sy*ex);
00422   if(deltaAngle>0. && _angle<0.)
00423     deltaAngle-=2.*M_PI;
00424   else if(deltaAngle<0. && _angle>0.)
00425     deltaAngle+=2.*M_PI;
00426   deltaAngle=direction?deltaAngle:-deltaAngle;
00427   return new EdgeArcCircle(start,end,_center,_radius,angle0,deltaAngle,direction);
00428 }
00429 
00430 void EdgeArcCircle::applySimilarity(double xBary, double yBary, double dimChar)
00431 {
00432   Edge::applySimilarity(xBary,yBary,dimChar);
00433   _radius/=dimChar;
00434   _center[0]=(_center[0]-xBary)/dimChar;
00435   _center[1]=(_center[1]-yBary)/dimChar;
00436 }
00437 
00438 void EdgeArcCircle::unApplySimilarity(double xBary, double yBary, double dimChar)
00439 {
00440   Edge::unApplySimilarity(xBary,yBary,dimChar);
00441   _radius*=dimChar;
00442   _center[0]=_center[0]*dimChar+xBary;
00443   _center[1]=_center[1]*dimChar+yBary;
00444 }
00445 
00451 void EdgeArcCircle::tesselate(const int *conn, int offset, double eps, std::vector<int>& newConn, std::vector<double>& addCoo) const
00452 {
00453   newConn.push_back(INTERP_KERNEL::NORM_POLYL);
00454   int nbOfSubDiv=fabs(_angle)/eps;
00455   if(nbOfSubDiv<=2)
00456     {
00457       newConn.push_back(conn[0]); newConn.push_back(conn[2]); newConn.push_back(conn[1]);
00458       return ;
00459     }
00460   double signOfAngle=_angle>0.?1.:-1.;
00461   int offset2=offset+((int)addCoo.size())/2;
00462   newConn.push_back(conn[0]);
00463   for(int i=1;i<nbOfSubDiv;i++,offset2++)
00464     {
00465       double angle=_angle0+i*eps*signOfAngle;
00466       newConn.push_back(offset2);
00467       addCoo.push_back(_center[0]+_radius*cos(angle)); addCoo.push_back(_center[1]+_radius*sin(angle));
00468     }
00469   newConn.push_back(conn[1]);
00470 }
00471 
00472 EdgeArcCircle *EdgeArcCircle::BuildFromNodes(Node *start, Node *middle, Node *end)
00473 {
00474   EdgeLin *e1,*e2;
00475   e1=new EdgeLin(start,middle);
00476   e2=new EdgeLin(middle,end);
00477   SegSegIntersector inters(*e1,*e2);
00478   bool colinearity=inters.areColinears();
00479   delete e1; delete e2;
00480   if(colinearity)
00481     {
00482       start->decrRef(); middle->decrRef(); end->decrRef();
00483       return 0;
00484     }
00485   else
00486     {
00487       EdgeArcCircle *ret=new EdgeArcCircle(start,middle,end);
00488       start->decrRef(); middle->decrRef(); end->decrRef();
00489       return ret;
00490     }
00491 }
00492 
00497 double EdgeArcCircle::GetAbsoluteAngle(const double *vect, double& normVect)
00498 {
00499   normVect=Node::norm(vect);
00500   return GetAbsoluteAngleOfNormalizedVect(vect[0]/normVect,vect[1]/normVect);
00501 }
00502 
00510 double EdgeArcCircle::GetAbsoluteAngleOfNormalizedVect(double ux, double uy)
00511 {
00512   //When arc is lower than 0.707 Using Asin 
00513   if(fabs(ux)<0.707)
00514     {
00515       double ret=SafeAcos(ux);
00516       if(uy>0.)
00517         return ret;
00518       ret=-ret;
00519       return ret;
00520     }
00521   else
00522     {
00523       double ret=SafeAsin(uy);
00524       if(ux>0.)
00525         return ret;
00526       if(ret>0.)
00527         return M_PI-ret;
00528       else
00529         return -M_PI-ret;
00530     }
00531 }
00532 
00533 void EdgeArcCircle::GetArcOfCirclePassingThru(const double *start, const double *middle, const double *end, 
00534                                               double *center, double& radius, double& angleInRad, double& angleInRad0)
00535 {
00536   double delta=(middle[0]-start[0])*(end[1]-middle[1])-(end[0]-middle[0])*(middle[1]-start[1]);
00537   double b1=(middle[1]*middle[1]+middle[0]*middle[0]-start[0]*start[0]-start[1]*start[1])/2;
00538   double b2=(end[1]*end[1]+end[0]*end[0]-middle[0]*middle[0]-middle[1]*middle[1])/2;
00539   center[0]=((end[1]-middle[1])*b1+(start[1]-middle[1])*b2)/delta;
00540   center[1]=((middle[0]-end[0])*b1+(middle[0]-start[0])*b2)/delta;
00541   radius=SafeSqrt((start[0]-center[0])*(start[0]-center[0])+(start[1]-center[1])*(start[1]-center[1]));
00542   angleInRad0=GetAbsoluteAngleOfNormalizedVect((start[0]-center[0])/radius,(start[1]-center[1])/radius);
00543   double angleInRadM=GetAbsoluteAngleOfNormalizedVect((middle[0]-center[0])/radius,(middle[1]-center[1])/radius);
00544   angleInRad=GetAbsoluteAngleOfNormalizedVect(((start[0]-center[0])*(end[0]-center[0])+(start[1]-center[1])*(end[1]-center[1]))/(radius*radius),
00545                                               ((start[0]-center[0])*(end[1]-center[1])-(start[1]-center[1])*(end[0]-center[0]))/(radius*radius));
00546   if(IsAngleNotIn(angleInRad0,angleInRad,angleInRadM))
00547     angleInRad=angleInRad<0?2*M_PI+angleInRad:angleInRad-2*M_PI;
00548 }
00549 
00550 void EdgeArcCircle::dumpInXfigFile(std::ostream& stream, bool direction, int resolution, const Bounds& box) const
00551 {
00552   stream << "5 1 0 1 ";
00553   fillXfigStreamForLoc(stream);
00554   stream << " 7 50 -1 -1 0.000 0 ";
00555   if( (direction && (-_angle)>=0) || (!direction && (-_angle)<0))
00556     stream << '0';//'0'
00557   else
00558     stream << '1';//'1'
00559   stream << " 1 0 ";
00560   stream << box.fitXForXFigD(_center[0],resolution) << " " << box.fitYForXFigD(_center[1],resolution) << " ";
00561   direction?_start->dumpInXfigFile(stream,resolution,box):_end->dumpInXfigFile(stream,resolution,box);
00562   Node *middle=buildRepresentantOfMySelf();
00563   middle->dumpInXfigFile(stream,resolution,box);
00564   middle->decrRef();
00565   direction?_end->dumpInXfigFile(stream,resolution,box):_start->dumpInXfigFile(stream,resolution,box);
00566   stream << std::endl << "1 1 2.00 120.00 180.00" << std::endl;
00567 }
00568 
00569 void EdgeArcCircle::update(Node *m)
00570 {
00571   GetArcOfCirclePassingThru(*_start,*m,*_end,_center,_radius,_angle,_angle0);
00572   updateBounds();
00573 }
00574 
00581 double EdgeArcCircle::getAreaOfZone() const
00582 {
00583   return -_radius*_radius*(sin(_angle)-_angle)/2.+((*_start)[0]-(*_end)[0])*((*_start)[1]+(*_end)[1])/2.;
00584 }
00585 
00586 double EdgeArcCircle::getCurveLength() const
00587 {
00588   return fabs(_angle*_radius);
00589 }
00590 
00591 void EdgeArcCircle::getBarycenter(double *bary) const
00592 {
00593   bary[0]=_center[0]+_radius*cos(_angle0+_angle/2.);
00594   bary[1]=_center[1]+_radius*sin(_angle0+_angle/2.);
00595 }
00596 
00615 void EdgeArcCircle::getBarycenterOfZone(double *bary) const
00616 {
00617   double x0=_center[0];
00618   double y0=_center[1];
00619   double angle1=_angle0+_angle;
00620   double tmp1=sin(angle1);
00621   double tmp0=sin(_angle0);
00622   double tmp2=_radius*_radius*_radius;
00623   double tmp3=cos(angle1);
00624   double tmp4=cos(_angle0);
00625   bary[0]=_radius*x0*y0*(tmp4-tmp3)+_radius*_radius*(y0*(cos(2*_angle0)-cos(2*angle1))/4.+
00626                                                      x0*(_angle/2.+(sin(2.*_angle0)-sin(2.*angle1))/4.))
00627     +tmp2*(tmp1*tmp1*tmp1-tmp0*tmp0*tmp0)/3.;
00628   bary[1]=y0*y0*_radius*(tmp4-tmp3)/2.+_radius*_radius*y0*(_angle/2.+(sin(2.*_angle0)-sin(2.*angle1))/4.)
00629     +tmp2*(tmp4-tmp3+(tmp3*tmp3*tmp3-tmp4*tmp4*tmp4)/3.)/2.;
00630 }
00631 
00635 bool EdgeArcCircle::isIn(double characterVal) const
00636 {
00637   return IsIn2Pi(_angle0,_angle,characterVal);
00638 }
00639 
00640 Node *EdgeArcCircle::buildRepresentantOfMySelf() const
00641 {
00642   return new Node(_center[0]+_radius*cos(_angle0+_angle/2.),_center[1]+_radius*sin(_angle0+_angle/2.));
00643 }
00644 
00649 bool EdgeArcCircle::isLower(double val1, double val2) const
00650 {
00651   double myDelta1=val1-_angle0;
00652   double myDelta2=val2-_angle0;
00653   if(_angle>0.)
00654     {
00655       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.
00656       myDelta2=myDelta2>-(_radius*QUADRATIC_PLANAR::_precision)?myDelta2:myDelta2+2.*M_PI;
00657       return myDelta1<myDelta2;
00658     }
00659   else
00660     {
00661       myDelta1=myDelta1<(_radius*QUADRATIC_PLANAR::_precision)?myDelta1:myDelta1-2.*M_PI;
00662       myDelta2=myDelta2<(_radius*QUADRATIC_PLANAR::_precision)?myDelta2:myDelta2-2.*M_PI;
00663       return myDelta2<myDelta1;
00664     }
00665 }
00666 
00670 double EdgeArcCircle::getCharactValue(const Node& node) const
00671 {
00672   double dx=(node[0]-_center[0])/_radius;
00673   double dy=(node[1]-_center[1])/_radius;
00674   return GetAbsoluteAngleOfNormalizedVect(dx,dy);
00675 }
00676 
00677 double EdgeArcCircle::getCharactValueBtw0And1(const Node& node) const
00678 {
00679   double dx=(node[0]-_center[0])/_radius;
00680   double dy=(node[1]-_center[1])/_radius;
00681   double angle=GetAbsoluteAngleOfNormalizedVect(dx,dy);
00682   //
00683   double myDelta=angle-_angle0;
00684   if(_angle>0.)
00685     myDelta=myDelta>=0.?myDelta:myDelta+2.*M_PI;
00686   else
00687     myDelta=myDelta<=0.?myDelta:myDelta-2.*M_PI;
00688   return myDelta/_angle;
00689 }
00690 
00691 double EdgeArcCircle::getDistanceToPoint(const double *pt) const
00692 {
00693   double angle=Node::computeAngle(_center,pt);
00694   if(IsIn2Pi(_angle0,_angle,angle))
00695     return fabs(Node::distanceBtw2Pt(_center,pt)-_radius);
00696   else
00697     {
00698       double dist1=Node::distanceBtw2Pt(*_start,pt);
00699       double dist2=Node::distanceBtw2Pt(*_end,pt);
00700       return std::min(dist1,dist2);
00701     }
00702 }
00703 
00704 bool EdgeArcCircle::isNodeLyingOn(const double *coordOfNode) const
00705 {
00706   double dist=Node::distanceBtw2Pt(_center,coordOfNode);
00707   if(Node::areDoubleEquals(dist,_radius))
00708     {
00709       double angle=Node::computeAngle(_center,coordOfNode);
00710       return IsIn2Pi(_angle0,_angle,angle);
00711     }
00712   else
00713     return false;
00714 }
00715 
00720 bool EdgeArcCircle::IsIn2Pi(double start, double delta, double angleIn)
00721 {
00722   double myDelta=angleIn-start;
00723   if(delta>0.)
00724     {
00725       myDelta=myDelta>=0.?myDelta:myDelta+2.*M_PI;
00726       return myDelta>0. && myDelta<delta;
00727     }
00728   else
00729     {
00730       myDelta=myDelta<=0.?myDelta:myDelta-2.*M_PI;
00731       return myDelta<0. && myDelta>delta;
00732     }
00733 }
00734 
00738 bool EdgeArcCircle::IsAngleNotIn(double start, double delta, double angleIn)
00739 {
00740   double tmp=start;
00741   if(tmp<0.)
00742     tmp+=2*M_PI;
00743   double tmp2=angleIn;
00744   if(tmp2<0.)
00745     tmp2+=2*M_PI;
00746   if(tmp+delta>=2.*M_PI)
00747     return (tmp2<tmp) && (tmp2>tmp+delta-2*M_PI);
00748   else if(tmp+delta>=0.)
00749     return (tmp2<std::min(tmp,tmp+delta) || tmp2>std::max(tmp,tmp+delta));
00750   else
00751     return (tmp2>tmp) && (tmp2<(tmp+delta+2.*M_PI));
00752 }
00753 
00754 void EdgeArcCircle::updateBounds()
00755 {
00756   _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]));
00757   if(IsIn2Pi(_angle0,_angle,M_PI/2))
00758     _bounds[3]=_center[1]+_radius;
00759   if(IsIn2Pi(_angle0,_angle,-M_PI/2))
00760     _bounds[2]=_center[1]-_radius;
00761   if(IsIn2Pi(_angle0,_angle,0.))
00762     _bounds[1]=_center[0]+_radius;
00763   if(IsIn2Pi(_angle0,_angle,M_PI))
00764     _bounds[0]=_center[0]-_radius;
00765 }
00766 
00767 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,
00768                                       std::vector<int>& edgesThis, std::vector<double>& addCoo, std::map<INTERP_KERNEL::Node *,int> mapAddCoo) const
00769 {
00770   int tmp[2];
00771   _start->fillGlobalInfoAbs(mapThis,mapOther,offset1,offset2,fact,baryX,baryY,addCoo,mapAddCoo,tmp);
00772   _end->fillGlobalInfoAbs(mapThis,mapOther,offset1,offset2,fact,baryX,baryY,addCoo,mapAddCoo,tmp+1);
00773   if(direction)
00774     {
00775       edgesThis.push_back(tmp[0]);
00776       edgesThis.push_back(tmp[1]);
00777     }
00778   else
00779     {
00780       edgesThis.push_back(tmp[1]);
00781       edgesThis.push_back(tmp[0]);
00782     }
00783 }
00784 
00785 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,
00786                                        std::vector<int>& edgesOther, std::vector<double>& addCoo, std::map<INTERP_KERNEL::Node *,int>& mapAddCoo) const
00787 {
00788   _start->fillGlobalInfoAbs2(mapThis,mapOther,offset1,offset2,fact,baryX,baryY,addCoo,mapAddCoo,edgesOther);
00789   _end->fillGlobalInfoAbs2(mapThis,mapOther,offset1,offset2,fact,baryX,baryY,addCoo,mapAddCoo,edgesOther);
00790 }