src/IntAna/IntAna_Curve.cxx

   1 // Created on: 1992-06-30
   2 // Created by: Laurent BUCHARD
   3 // Copyright (c) 1992-1999 Matra Datavision
   4 // Copyright (c) 1999-2012 OPEN CASCADE SAS
   5 //
   6 // The content of this file is subject to the Open CASCADE Technology Public
   7 // License Version 6.5 (the "License"). You may not use the content of this file
   8 // except in compliance with the License. Please obtain a copy of the License
   9 // at http://www.opencascade.org and read it completely before using this file.
  10 //
  11 // The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
  12 // main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
  13 //
  14 // The Original Code and all software distributed under the License is
  15 // distributed on an "AS IS" basis, without warranty of any kind, and the
  16 // Initial Developer hereby disclaims all such warranties, including without
  17 // limitation, any warranties of merchantability, fitness for a particular
  18 // purpose or non-infringement. Please see the License for the specific terms
  19 // and conditions governing the rights and limitations under the License.
  20
  21
  22 #ifndef DEB
  23 #define No_Standard_RangeError
  24 #define No_Standard_OutOfRange
  25 #endif
  26
  27 //----------------------------------------------------------------------
  28 //-- Differents constructeurs sont proposes qui correspondent aux
  29 //-- polynomes en Z :
  30 //--    A(Sin(Theta),Cos(Theta)) Z**2
  31 //--  + B(Sin(Theta),Cos(Theta)) Z
  32 //--  + C(Sin(Theta),Cos(Theta))
  33 //--
  34 //-- Une Courbe est definie sur un domaine
  35 //--
  36 //-- Value retourne le point de parametre U(Theta),V(Theta)
  37 //--       ou V est la solution du polynome A V**2 + B V + C
  38 //--       (Selon les cas, on prend V+ ou V-)
  39 //--
  40 //-- D1u   calcule le vecteur tangent a la courbe
  41 //--       et retourne le booleen Standard_False si ce calcul ne peut
  42 //--       pas etre mene a bien.
  43 //----------------------------------------------------------------------
  44
  45 #include <Precision.hxx>
  46
  47 #include <IntAna_Curve.ixx>
  48
  49 #include <Standard_DomainError.hxx>
  50 #include <math_DirectPolynomialRoots.hxx>
  51 #include <ElSLib.hxx>
  52 #include <gp_XYZ.hxx>
  53
  54 //=======================================================================
  55 //function : IntAna_Curve
  56 //purpose  :
  57 //=======================================================================
  58   IntAna_Curve::IntAna_Curve()
  59 {
  60   typequadric=GeomAbs_OtherSurface;
  61   firstbounded=Standard_False;
  62   lastbounded=Standard_False;
  63 }
  64 //=======================================================================
  65 //function : SetConeQuadValues
  66 //purpose  : Description de l intersection Cone Quadrique
  67 //=======================================================================
  68   void IntAna_Curve::SetConeQuadValues(const gp_Cone& Cone,
  69                                        const Standard_Real Qxx,
  70                                        const Standard_Real Qyy,
  71                                        const Standard_Real Qzz,
  72                                        const Standard_Real Qxy,
  73                                        const Standard_Real Qxz,
  74                                        const Standard_Real Qyz,
  75                                        const Standard_Real Qx,
  76                                        const Standard_Real Qy,
  77                                        const Standard_Real Qz,
  78                                        const Standard_Real Q1,
  79                                        const Standard_Real TOL,
  80                                        const Standard_Real DomInf,
  81                                        const Standard_Real DomSup,
  82                                        const Standard_Boolean twocurves,
  83                                        const Standard_Boolean takezpositive)
  84 {
  85
  86   Ax3        = Cone.Position();
  87   RCyl       = Cone.RefRadius();
  88
  89   Angle      = Cone.SemiAngle();
  90   Standard_Real UnSurTgAngle = 1.0/(Tan(Cone.SemiAngle()));
  91
  92   typequadric= GeomAbs_Cone;
  93
  94   TwoCurves     = twocurves;         //-- deux  Z pour un meme parametre
  95   TakeZPositive = takezpositive;     //-- Prendre sur la courbe le Z Positif
  96                                      //--   ( -B + Sqrt()) et non (-B - Sqrt())
  97
  98
  99   Z0Cte      = Q1;                   //-- Attention On a    Z?Cos Cos(t)
 100   Z0Sin      = 0.0;                  //-- et Non          2 Z?Cos Cos(t) !!!
 101   Z0Cos      = 0.0;                  //-- Ce pour tous les Parametres
 102   Z0CosCos   = 0.0;                  //--  ie pas de Coefficient 2
 103   Z0SinSin   = 0.0;                  //--     devant les termes CS C S
 104   Z0CosSin   = 0.0;
 105
 106   Z1Cte      = 2.0*(UnSurTgAngle)*Qz;
 107   Z1Sin      = Qy+Qy;
 108   Z1Cos      = Qx+Qx;
 109   Z1CosCos   = 0.0;
 110   Z1SinSin   = 0.0;
 111   Z1CosSin   = 0.0;
 112
 113   Z2Cte      = Qzz * UnSurTgAngle*UnSurTgAngle;
 114   Z2Sin      = (UnSurTgAngle+UnSurTgAngle)*Qyz;
 115   Z2Cos      = (UnSurTgAngle+UnSurTgAngle)*Qxz;
 116   Z2CosCos   = Qxx;
 117   Z2SinSin   = Qyy;
 118   Z2CosSin   = Qxy+Qxy;
 119
 120   Tolerance  = TOL;
 121   DomainInf  = DomInf;
 122   DomainSup  = DomSup;
 123
 124   RestrictedInf = RestrictedSup = Standard_True;   //-- Le Domaine est Borne
 125   firstbounded = lastbounded = Standard_False;
 126 }
 127
 128 //=======================================================================
 129 //function : SetCylinderQuadValues
 130 //purpose  : Description de l intersection Cylindre Quadrique
 131 //=======================================================================
 132   void IntAna_Curve::SetCylinderQuadValues(const gp_Cylinder& Cyl,
 133                                            const Standard_Real Qxx,
 134                                            const Standard_Real Qyy,
 135                                            const Standard_Real Qzz,
 136                                            const Standard_Real Qxy,
 137                                            const Standard_Real Qxz,
 138                                            const Standard_Real Qyz,
 139                                            const Standard_Real Qx,
 140                                            const Standard_Real Qy,
 141                                            const Standard_Real Qz,
 142                                            const Standard_Real Q1,
 143                                            const Standard_Real TOL,
 144                                            const Standard_Real DomInf,
 145                                            const Standard_Real DomSup,
 146                                            const Standard_Boolean twocurves,
 147                                            const Standard_Boolean takezpositive)
 148 {
 149
 150   Ax3        = Cyl.Position();
 151   RCyl       = Cyl.Radius();
 152   typequadric= GeomAbs_Cylinder;
 153
 154   TwoCurves     = twocurves;         //-- deux  Z pour un meme parametre
 155   TakeZPositive = takezpositive;     //-- Prendre sur la courbe le Z Positif
 156   Standard_Real RCylmul2 = RCyl+RCyl;         //--   ( -B + Sqrt())
 157
 158   Z0Cte      = Q1;
 159   Z0Sin      = RCylmul2*Qy;
 160   Z0Cos      = RCylmul2*Qx;
 161   Z0CosCos   = Qxx*RCyl*RCyl;
 162   Z0SinSin   = Qyy*RCyl*RCyl;
 163   Z0CosSin   = RCylmul2*RCyl*Qxy;
 164
 165   Z1Cte      = Qz+Qz;
 166   Z1Sin      = RCylmul2*Qyz;
 167   Z1Cos      = RCylmul2*Qxz;
 168   Z1CosCos   = 0.0;
 169   Z1SinSin   = 0.0;
 170   Z1CosSin   = 0.0;
 171
 172   Z2Cte      = Qzz;
 173   Z2Sin      = 0.0;
 174   Z2Cos      = 0.0;
 175   Z2CosCos   = 0.0;
 176   Z2SinSin   = 0.0;
 177   Z2CosSin   = 0.0;
 178
 179   Tolerance  = TOL;
 180   DomainInf  = DomInf;
 181   DomainSup  = DomSup;
 182
 183   RestrictedInf = RestrictedSup = Standard_True;
 184   firstbounded = lastbounded = Standard_False;
 185 }
 186
 187 //=======================================================================
 188 //function : IsOpen
 189 //purpose  :
 190 //=======================================================================
 191   Standard_Boolean IntAna_Curve::IsOpen() const
 192 {
 193   return(RestrictedInf && RestrictedSup);
 194 }
 195
 196 //=======================================================================
 197 //function : Domain
 198 //purpose  :
 199 //=======================================================================
 200   void IntAna_Curve::Domain(Standard_Real& DInf,
 201                             Standard_Real& DSup) const
 202 {
 203   if(RestrictedInf && RestrictedSup) {
 204     DInf=DomainInf;
 205     DSup=DomainSup;
 206     if(TwoCurves) {
 207       DSup+=DSup-DInf;
 208     }
 209   }
 210   else {
 211     Standard_DomainError::Raise("IntAna_Curve::Domain");
 212   }
 213 }
 214 //=======================================================================
 215 //function : IsConstant
 216 //purpose  :
 217 //=======================================================================
 218   Standard_Boolean IntAna_Curve::IsConstant() const
 219 {
 220   //-- ???  Pas facile de decider a la seule vue des Param.
 221   return(Standard_False);
 222 }
 223
 224 //=======================================================================
 225 //function : IsFirstOpen
 226 //purpose  :
 227 //=======================================================================
 228   Standard_Boolean IntAna_Curve::IsFirstOpen() const
 229 {
 230   return(firstbounded);
 231 }
 232
 233 //=======================================================================
 234 //function : IsLastOpen
 235 //purpose  :
 236 //=======================================================================
 237   Standard_Boolean IntAna_Curve::IsLastOpen() const
 238 {
 239   return(lastbounded);
 240 }
 241 //=======================================================================
 242 //function : SetIsFirstOpen
 243 //purpose  :
 244 //=======================================================================
 245   void IntAna_Curve::SetIsFirstOpen(const Standard_Boolean Flag)
 246 {
 247   firstbounded = Flag;
 248 }
 249
 250 //=======================================================================
 251 //function : SetIsLastOpen
 252 //purpose  :
 253 //=======================================================================
 254   void IntAna_Curve::SetIsLastOpen(const Standard_Boolean Flag)
 255 {
 256   lastbounded = Flag;
 257 }
 258
 259 //=======================================================================
 260 //function : InternalUVValue
 261 //purpose  :
 262 //=======================================================================
 263   void IntAna_Curve::InternalUVValue(const Standard_Real theta,
 264                                      Standard_Real& Param1,
 265                                      Standard_Real& Param2,
 266                                      Standard_Real& A,
 267                                      Standard_Real& B,
 268                                      Standard_Real& C,
 269                                      Standard_Real& cost,
 270                                      Standard_Real& sint,
 271                                      Standard_Real& SigneSqrtDis) const
 272 {
 273   Standard_Real Theta=theta;
 274   Standard_Boolean SecondSolution=Standard_False;
 275
 276   if((Theta<DomainInf)  ||
 277      ((Theta>DomainSup) && (!TwoCurves)) ||
 278      (Theta>(DomainSup+DomainSup-DomainInf+0.00000000000001))) {
 279     Standard_DomainError::Raise("IntAna_Curve::Domain");
 280   }
 281
 282   if(Theta>DomainSup) {
 283     Theta=DomainSup+DomainSup-Theta;
 284     SecondSolution=Standard_True;
 285   }
 286
 287   Param1=Theta;
 288
 289   if(!TwoCurves) {
 290     SecondSolution=TakeZPositive;
 291   }
 292   //
 293   cost = Cos(Theta);
 294   sint = Sin(Theta);
 295   Standard_Real costsint = cost*sint;
 296
 297   A=Z2Cte+sint*(Z2Sin+sint*Z2SinSin)+cost*(Z2Cos+cost*Z2CosCos)
 298     +Z2CosSin*costsint;
 299
 300   B=Z1Cte+sint*(Z1Sin+sint*Z1SinSin)+cost*(Z1Cos+cost*Z1CosCos)
 301     +Z1CosSin*costsint;
 302
 303   C=Z0Cte+sint*(Z0Sin+sint*Z0SinSin)+cost*(Z0Cos+cost*Z0CosCos)
 304     +Z0CosSin*costsint;
 305
 306
 307   Standard_Real Discriminant = B*B-4.0*A*C;
 308
 309   if(Abs(A)<=0.000000001) {
 310     //-- cout<<" IntAna_Curve:: Internal UV Value : A="<<A<<"  -> Abs(A)="<<Abs(A)<<endl;
 311     if(Abs(B)<=0.0000000001) {
 312       //-- cout<<" Probleme :  Pas de solutions "<<endl;
 313       Param2=0.0;
 314     }
 315     else {
 316       //modified by NIZNHY-PKV Fri Dec  2 16:02:46 2005f
 317       Param2 = -C/B;
 318       /*
 319       if(!SecondSolution) {
 320         //-- Cas Param2 = (-B+Sqrt(Discriminant))/(A+A);
 321         //--            = (-B+Sqrt(B**2 - Eps)) / 2A
 322         //--            = -C / B
 323         Param2 = -C/B;
 324       }
 325       else {
 326         //-- Cas Param2 = (-B-Sqrt(Discriminant))/(A+A);
 327         //--            = (-B-Sqrt(B**2 - Eps)) / 2A
 328         if(A) {
 329           Param2 = -B/A;
 330         }
 331         else {
 332           Param2 = -B*10000000.0;
 333         }
 334       }
 335       */
 336       //modified by NIZNHY-PKV Fri Dec  2 16:02:54 2005t
 337     }
 338   }
 339   else {
 340     if(Discriminant<=0.0000000001 ||
 341        Abs(Discriminant/(4*A))<=0.0000000001) Discriminant=0.0;
 342     SigneSqrtDis = (SecondSolution)? Sqrt(Discriminant)
 343       : -Sqrt(Discriminant);
 344     Param2=(-B+SigneSqrtDis)/(A+A);
 345   }
 346 }
 347
 348 //=======================================================================
 349 //function : Value
 350 //purpose  :
 351 //=======================================================================
 352   gp_Pnt IntAna_Curve::Value(const Standard_Real theta)
 353 {
 354   Standard_Real A, B, C, U, V, sint, cost, SigneSqrtDis;
 355   //
 356   A=0.0;  B=0.0;   C=0.0;
 357   U=0.0;  V=0.0;
 358   sint=0.0;  cost=0.0;
 359   SigneSqrtDis=0.0;
 360   InternalUVValue(theta,U,V,A,B,C,cost,sint,SigneSqrtDis);
 361   //-- checked the parameter U and Raises Domain Error if Error
 362   return(InternalValue(U,V));
 363 }
 364 //=======================================================================
 365 //function : D1u
 366 //purpose  :
 367 //=======================================================================
 368   Standard_Boolean IntAna_Curve::D1u(const Standard_Real theta,
 369                                      gp_Pnt& Pt,
 370                                      gp_Vec& Vec)
 371 {
 372   //-- Pour detecter le cas ou le calcul est impossible
 373   Standard_Real A, B, C, U, V, sint, cost, SigneSqrtDis;
 374   A=0.0;  B=0.0;   C=0.0;
 375   U=0.0;  V=0.0;
 376   sint=0.0;  cost=0.0;
 377   //
 378   InternalUVValue(theta,U,V,A,B,C,cost,sint,SigneSqrtDis);
 379   //
 380   Pt = Value(theta);
 381   if(Abs(SigneSqrtDis)<0.0000000001  || Abs(A)<0.0000001) return(Standard_False);
 382
 383
 384   //-- Approximation de la derivee (mieux que le calcul mathematique!)
 385   Standard_Real dtheta = (DomainSup-DomainInf)*0.000001;
 386   Standard_Real theta2 = theta+dtheta;
 387   if((theta2<DomainInf)  ||  ((theta2>DomainSup) && (!TwoCurves))
 388      ||  (theta2>(DomainSup+DomainSup-DomainInf+0.00000000000001))) {
 389     dtheta = -dtheta;
 390     theta2 = theta+dtheta;
 391   }
 392   gp_Pnt P2 = Value(theta2);
 393   dtheta = 1.0/dtheta;
 394   Vec.SetCoord((P2.X()-Pt.X())*dtheta,
 395                (P2.Y()-Pt.Y())*dtheta,
 396                (P2.Z()-Pt.Z())*dtheta);
 397
 398   return(Standard_True);
 399 }
 400 //=======================================================================
 401 //function : FindParameter
 402 //purpose  : Para est en sortie le parametre sur la courbe
 403 //=======================================================================
 404   Standard_Boolean IntAna_Curve::FindParameter (const gp_Pnt& P,
 405                                                 Standard_Real& Para) const
 406 {
 407   Standard_Real theta,z, aTolPrecision=0.0001;
 408   Standard_Real PIpPI = M_PI + M_PI;
 409   //
 410   switch (typequadric) {
 411
 412   case GeomAbs_Cylinder:
 413     {
 414       ElSLib::CylinderParameters(Ax3,RCyl,P,theta,z);
 415     }
 416     break;
 417
 418   case GeomAbs_Cone :
 419     {
 420       ElSLib::ConeParameters(Ax3,RCyl,Angle,P,theta,z);
 421     }
 422     break;
 423
 424   default:
 425     break;
 426   }
 427   //
 428   Standard_Real epsAng = 1.e-8;
 429   Standard_Real tmin = DomainInf;
 430   Standard_Real tmax = DomainSup;
 431   Standard_Real U,V,A,B,C,sint,cost,SigneSqrtDis;
 432   Standard_Real z1,z2;
 433
 434   A=0.0;  B=0.0;   C=0.0;
 435   U=0.0;  V=0.0;
 436   sint=0.0;  cost=0.0;
 437   SigneSqrtDis=0.0;
 438   //U=V=A=B=C=sint=cost=SigneSqrtDis=0.0;
 439   //
 440   if (!firstbounded && tmin > theta && (tmin-theta) <= epsAng) {
 441     theta = tmin;
 442   }
 443   else if (!lastbounded && theta > tmax && (theta-tmax) <= epsAng) {
 444     theta = tmax;
 445   }
 446   //
 447   if (theta < tmin ) {
 448     theta = theta + PIpPI;
 449   }
 450   else if (theta > tmax) {
 451     theta = theta - PIpPI;
 452   }
 453   if (theta < tmin || theta > tmax) {
 454     if(theta>tmax) {
 455       InternalUVValue(tmax,U,V,A,B,C,cost,sint,SigneSqrtDis);
 456       gp_Pnt PMax(InternalValue(U,V));
 457       if(PMax.Distance(P) < aTolPrecision) {
 458         Para = tmax;
 459         return(Standard_True);
 460       }
 461     }
 462     if(theta<tmin) {
 463       InternalUVValue(tmin,U,V,A,B,C,cost,sint,SigneSqrtDis);
 464       gp_Pnt PMin(InternalValue(U,V));
 465       if(PMin.Distance(P) < aTolPrecision) {
 466         Para = tmin;
 467         return(Standard_True);
 468       }
 469     }
 470     //-- lbr le 14 Fev 96 : On teste malgre tout si le point n est pas le
 471     //-- point de debut ou de fin
 472     //-- cout<<"False 1 "<<endl;
 473     //    theta = tmin;                               le 25 Nov 96
 474   }
 475
 476   if (TwoCurves) {
 477     if(theta > tmax)
 478       theta =  tmax;
 479     if(theta < tmin)
 480       theta = tmin;
 481     InternalUVValue(theta,U,z1,A,B,C,cost,sint,SigneSqrtDis);
 482     A = B = C = sint = cost = SigneSqrtDis = 0.0;
 483     InternalUVValue(tmax+tmax - theta,U,z2,A,B,C,cost,sint,SigneSqrtDis);
 484
 485     if (Abs(z-z1) <= Abs(z-z2)) {
 486       Para = theta;
 487     }
 488     else {
 489       Para = tmax+tmax - theta;
 490     }
 491   }
 492   else {
 493     Para = theta;
 494   }
 495
 496   if((Para<DomainInf)  ||  ((Para>DomainSup) && (!TwoCurves))
 497      ||  (Para>(DomainSup+DomainSup-DomainInf+0.00000000000001))) {
 498     return(Standard_False);
 499   }
 500
 501   InternalUVValue(Para,U,V,A,B,C,cost,sint,SigneSqrtDis);
 502   gp_Pnt PPara = InternalValue(U,V);
 503   Standard_Real Dist = PPara.Distance(P);
 504   if(Dist > aTolPrecision) {
 505     //-- Il y a eu un probleme
 506     //-- On teste si le point est un point double
 507     InternalUVValue(tmin,U,V,A,B,C,cost,sint,SigneSqrtDis);
 508     PPara = InternalValue(U,V);
 509     Dist = PPara.Distance(P);
 510     if(Dist <= aTolPrecision) {
 511       Para = tmin;
 512       return(Standard_True);
 513     }
 514
 515     InternalUVValue(tmax,U,V,A,B,C,cost,sint,SigneSqrtDis);
 516     PPara = InternalValue(U,V);
 517     Dist = PPara.Distance(P);
 518     if(Dist <= aTolPrecision) {
 519       Para = tmax;
 520       return(Standard_True);
 521     }
 522     if (TwoCurves) {
 523       Standard_Real Theta = DomainSup+DomainSup-DomainInf;
 524       InternalUVValue(Theta,U,V,A,B,C,cost,sint,SigneSqrtDis);
 525       PPara = InternalValue(U,V);
 526       Dist = PPara.Distance(P);
 527       if(Dist <= aTolPrecision) {
 528         Para = Theta;
 529         return(Standard_True);
 530       }
 531     }
 532     return(Standard_False);
 533   }
 534   return(Standard_True);
 535 }
 536 //=======================================================================
 537 //function : InternalValue
 538 //purpose  :
 539 //=======================================================================
 540   gp_Pnt IntAna_Curve::InternalValue(const Standard_Real U,
 541                                      const Standard_Real _V) const
 542 {
 543   //-- cout<<" ["<<U<<","<<V<<"]";
 544   Standard_Real V = _V;
 545   if(V > 100000.0 )   {   V= 100000.0; }
 546   if(V < -100000.0 )  {   V=-100000.0; }
 547
 548   switch(typequadric) {
 549   case GeomAbs_Cone:
 550     {
 551       //------------------------------------------------
 552       //-- Parametrage : X = V * Cos(U)              ---
 553       //--               Y = V * Sin(U)              ---
 554       //--               Z = (V-RCyl) / Tan(SemiAngle)--
 555       //------------------------------------------------
 556       //-- Angle Vaut Cone.SemiAngle()
 557       return(ElSLib::ConeValue(U,(V-RCyl)/Sin(Angle),Ax3,RCyl,Angle));
 558     }
 559     break;
 560
 561   case GeomAbs_Cylinder:
 562     return(ElSLib::CylinderValue(U,V,Ax3,RCyl));
 563   case GeomAbs_Sphere:
 564     return(ElSLib::SphereValue(U,V,Ax3,RCyl));
 565   default:
 566     return(gp_Pnt(0.0,0.0,0.0));
 567   }
 568   return(gp_Pnt(0.0,0.0,0.0));
 569 }
 570
 571 //
 572 //=======================================================================
 573 //function : SetDomain
 574 //purpose  :
 575 //=======================================================================
 576   void IntAna_Curve::SetDomain(const Standard_Real DInf,
 577                                const Standard_Real DSup)
 578 {
 579   if(DInf>=DSup) {
 580     Standard_DomainError::Raise("IntAna_Curve::Domain");
 581   }
 582   //
 583   DomainInf=DInf;
 584   DomainSup=DSup;
 585 }