src/Extrema/Extrema_ExtPExtS.cxx

   1 // Created on: 1999-09-16
   2 // Created by: Edward AGAPOV
   3 // Copyright (c) 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 #include <Standard_NotImplemented.hxx>
  23 #include <Standard_OutOfRange.hxx>
  24 #include <StdFail_NotDone.hxx>
  25 #include <Adaptor3d_HCurve.hxx>
  26 #include <ElCLib.hxx>
  27 #include <Extrema_ExtPElC.hxx>
  28 #include <Extrema_ExtPExtS.hxx>
  29 #include <Extrema_POnCurv.hxx>
  30 #include <Extrema_POnSurf.hxx>
  31 #include <Precision.hxx>
  32 #include <gp.hxx>
  33 #include <gp_Ax2.hxx>
  34 #include <gp_Dir.hxx>
  35 #include <gp_Lin.hxx>
  36 #include <gp_Pln.hxx>
  37 #include <gp_Pnt.hxx>
  38 #include <gp_Vec.hxx>
  39 #include <math_FunctionSetRoot.hxx>
  40 #include <math_Vector.hxx>
  41 #include <Adaptor3d_SurfaceOfLinearExtrusion.hxx>
  42
  43 static gp_Ax2 GetPosition (const Handle(Adaptor3d_HCurve)& C);
  44
  45 static void PerformExtPElC (Extrema_ExtPElC& E,
  46                             const gp_Pnt& P,
  47                             const Handle(Adaptor3d_HCurve)& C,
  48                             const Standard_Real Tol);
  49
  50 static Standard_Boolean
  51   IsCaseAnalyticallyComputable (const GeomAbs_CurveType& theType,
  52                                 const gp_Ax2& theCurvePos,
  53                                 const gp_Dir& theSurfaceDirection) ;
  54
  55 static gp_Pnt GetValue(const Standard_Real U,
  56                        const Handle(Adaptor3d_HCurve)& C);
  57 //=======================================================================
  58 //function : Project
  59 //purpose  : Returns the projection of a point <Point> on a plane
  60 //           <ThePlane>  along a direction <TheDir>.
  61 //=======================================================================
  62 static gp_Pnt ProjectPnt(const gp_Ax2& ThePlane,
  63                          const gp_Dir& TheDir,
  64                          const gp_Pnt& Point)
  65 {
  66   gp_Vec PO(Point,ThePlane.Location());
  67   Standard_Real Alpha = PO * gp_Vec(ThePlane.Direction());
  68   Alpha /= TheDir * ThePlane.Direction();
  69   gp_Pnt P;
  70   P.SetXYZ(Point.XYZ() + Alpha * TheDir.XYZ());
  71   return P;
  72 }
  73
  74 //=======================================================================
  75 //function : IsOriginalPnt
  76 //purpose  :
  77 //=======================================================================
  78
  79 static Standard_Boolean IsOriginalPnt (const gp_Pnt& P,
  80                                        const Extrema_POnSurf* Points,
  81                                        const Standard_Integer NbPoints)
  82 {
  83   for (Standard_Integer i=1; i<=NbPoints; i++) {
  84     if (Points[i-1].Value().IsEqual(P, Precision::Confusion())) {
  85       return Standard_False;
  86     }
  87   }
  88   return Standard_True;
  89 }
  90
  91 //=======================================================================
  92 //function : MakePreciser
  93 //purpose  :
  94 //=======================================================================
  95
  96 void Extrema_ExtPExtS::MakePreciser (Standard_Real& U,
  97                                      const gp_Pnt& P,
  98                                      const Standard_Boolean isMin,
  99                                      const gp_Ax2& OrtogSection) const
 100 {
 101   if (U>myusup) {
 102     U = myusup;
 103   } else if (U<myuinf) {
 104     U = myuinf;
 105   } else {
 106
 107     Standard_Real step = (myusup - myuinf) / 30, D2e, D2next ,D2prev;
 108     gp_Pnt
 109       Pe = ProjectPnt (OrtogSection, myDirection, GetValue(U,myC)),
 110       Pprev = ProjectPnt (OrtogSection, myDirection, GetValue(U-step, myC)),
 111       Pnext = ProjectPnt (OrtogSection, myDirection, GetValue(U+step, myC));
 112     D2e = P.SquareDistance(Pe),
 113     D2next = P.SquareDistance(Pnext),
 114     D2prev = P.SquareDistance(Pprev);
 115     Standard_Boolean notFound;
 116     if (isMin)
 117       notFound =  (D2e > D2prev || D2e > D2next);
 118     else
 119       notFound = (D2e < D2prev || D2e < D2next);
 120
 121     if (notFound && (D2e < D2next && isMin)) {
 122       step = -step;
 123       D2next = D2prev;
 124       Pnext = Pprev;
 125     }
 126     while (notFound) {
 127       U = U + step;
 128       if (U > myusup) {
 129               U = myusup;
 130               break;
 131       }
 132       if (U < myuinf) {
 133               U = myuinf;
 134               break;
 135       }
 136       D2e = D2next;
 137       Pe = Pnext;
 138       Pnext = ProjectPnt (OrtogSection, myDirection, GetValue(U+step, myC));
 139       D2next = P.SquareDistance(Pnext);
 140       if (isMin)
 141               notFound = D2e > D2next;
 142       else
 143               notFound = D2e < D2next;
 144     }
 145   }
 146 }
 147 //=============================================================================
 148
 149 Extrema_ExtPExtS::Extrema_ExtPExtS ()
 150 {
 151   myDone = Standard_False;
 152 }
 153
 154 //=============================================================================
 155
 156 Extrema_ExtPExtS::Extrema_ExtPExtS (const gp_Pnt&       P,
 157                                     const Adaptor3d_SurfaceOfLinearExtrusion&  S,
 158                                     const Standard_Real    Umin,
 159                                     const Standard_Real    Usup,
 160                                     const Standard_Real    Vmin,
 161                                     const Standard_Real    Vsup,
 162                                     const Standard_Real    TolU,
 163                                     const Standard_Real    TolV)
 164 {
 165   Initialize (S,
 166               Umin, Usup, Vmin, Vsup,
 167               TolU, TolV);
 168   Perform(P);
 169 }
 170 //=============================================================================
 171
 172 Extrema_ExtPExtS::Extrema_ExtPExtS (const gp_Pnt&       P,
 173                                     const Adaptor3d_SurfaceOfLinearExtrusion&  S,
 174                                     const Standard_Real    TolU,
 175                                     const Standard_Real    TolV)
 176 {
 177   Initialize (S,
 178               S.FirstUParameter(), S.LastUParameter(),
 179               S.FirstVParameter(), S.LastVParameter(),
 180               TolU, TolV);
 181   Perform(P);
 182 }
 183 //=======================================================================
 184 //function : Initialize
 185 //purpose  :
 186 //=======================================================================
 187
 188 void Extrema_ExtPExtS::Initialize(const Adaptor3d_SurfaceOfLinearExtrusion& S,
 189                                   const Standard_Real Uinf,
 190                                   const Standard_Real Usup,
 191                                   const Standard_Real Vinf,
 192                                   const Standard_Real Vsup,
 193                                   const Standard_Real TolU,
 194                                   const Standard_Real TolV)
 195 {
 196   myuinf=Uinf;
 197   myusup=Usup;
 198   mytolu=TolU;
 199
 200   myvinf=Vinf;
 201   myvsup=Vsup;
 202   mytolv=TolV;
 203
 204   Handle(Adaptor3d_HCurve) anACurve = S.BasisCurve();
 205
 206   myF.Initialize(S);
 207   myC = anACurve;
 208   myS = (Adaptor3d_SurfacePtr)&S;
 209   myPosition = GetPosition(myC);
 210   myDirection = S.Direction();
 211   myIsAnalyticallyComputable = //Standard_False;
 212     IsCaseAnalyticallyComputable (myC->GetType(),myPosition,myDirection);
 213
 214   if (!myIsAnalyticallyComputable)
 215
 216     myExtPS.Initialize(S, 32, 32,
 217                        Uinf, Usup, Vinf, Vsup,
 218                        TolU, TolV);
 219 }
 220
 221
 222 //=======================================================================
 223 //function : Perform
 224 //purpose  :
 225 //=======================================================================
 226
 227 void Extrema_ExtPExtS::Perform (const gp_Pnt& P)
 228 {
 229   const Standard_Integer NbExtMax = 4; //dimension of arrays
 230                                        //myPoint[] and mySqDist[]
 231                                        //For "analytical" case
 232   myDone = Standard_False;
 233   myNbExt = 0;
 234
 235   if (!myIsAnalyticallyComputable) {
 236     myExtPS.Perform(P);
 237     myDone = myExtPS.IsDone();
 238 //  modified by NIZHNY-EAP Wed Nov 17 12:59:08 1999 ___BEGIN___
 239     myNbExt = myExtPS.NbExt();
 240 //  modified by NIZHNY-EAP Wed Nov 17 12:59:09 1999 ___END___
 241     return;
 242   }
 243
 244   gp_Pnt Pe, Pp = ProjectPnt(myPosition,myDirection,P);
 245   Extrema_ExtPElC anExt;
 246   PerformExtPElC(anExt, Pp, myC, mytolu);
 247   if (!anExt.IsDone()) return;
 248
 249   gp_Ax2 anOrtogSection (P, myDirection);
 250   Standard_Real U,V;
 251   Standard_Boolean
 252     isMin,
 253     isSimpleCase =
 254       myDirection.IsParallel(myPosition.Direction(),Precision::Angular());
 255   Standard_Integer i, aNbExt = anExt.NbExt();
 256   math_Vector UV(1,2), Tol(1,2), UVinf(1,2), UVsup(1,2);
 257   Tol(1) = mytolu;   Tol(2) = mytolv;
 258   UVinf(1) = myuinf; UVinf(2) = myvinf;
 259   UVsup(1) = myusup; UVsup(2) = myvsup;
 260
 261
 262   for (i=1; i<=aNbExt; i++) {
 263     Extrema_POnCurv POC=anExt.Point(i);
 264     U = POC.Parameter();
 265     //// modified by jgv, 23.12.2008 for OCC17194 ////
 266     if (myC->IsPeriodic())
 267     {
 268       Standard_Real U2 = U;
 269       ElCLib::AdjustPeriodic(myuinf, myuinf + 2.*M_PI, Precision::PConfusion(), U, U2);
 270     }
 271     //////////////////////////////////////////////////
 272     gp_Pnt E = POC.Value();
 273     Pe = ProjectPnt(anOrtogSection, myDirection, E);
 274
 275     if (isSimpleCase) {
 276       V = gp_Vec(E,Pe) * gp_Vec(myDirection);
 277       // modified by NIZHNY-MKK  Thu Sep 18 14:46:14 2003.BEGIN
 278       //       myPoint[++myNbExt] = Extrema_POnSurf(U, V, Pe);
 279       //       myValue[myNbExt] = anExt.Value(i);
 280       myPoint[myNbExt] = Extrema_POnSurf(U, V, Pe);
 281       mySqDist[myNbExt] = anExt.SquareDistance(i);
 282       myNbExt++;
 283       if(myNbExt == NbExtMax)
 284       {
 285         break;
 286       }
 287       // modified by NIZHNY-MKK  Thu Sep 18 14:46:18 2003.END
 288     }
 289     else {
 290       myF.SetPoint(P);
 291       isMin = anExt.IsMin(i);//( Pp.Distance(GetValue(U+10,myC)) > anExt.Value(i) );
 292
 293       MakePreciser(U, P, isMin, anOrtogSection);
 294       E = GetValue(U, myC);
 295       Pe = ProjectPnt (anOrtogSection, myDirection, E),
 296       V = gp_Vec(E,Pe) * gp_Vec(myDirection);
 297       UV(1) = U;         UV(2) = V;
 298       math_FunctionSetRoot aFSR (myF,UV,Tol,UVinf,UVsup);
 299 //      for (Standard_Integer k=1 ; k <= myF.NbExt();
 300       Standard_Integer k;
 301       for ( k=1 ; k <= myF.NbExt(); k++) {
 302               if (IsOriginalPnt(myF.Point(k).Value(), myPoint, myNbExt)) {
 303                 // modified by NIZHNY-MKK  Thu Sep 18 14:46:41 2003.BEGIN
 304                 //        myPoint[++myNbExt] = myF.Point(k);
 305                 //        myValue[myNbExt] = myF.Value(k);
 306                 myPoint[myNbExt] = myF.Point(k);
 307                 mySqDist[myNbExt] = myF.SquareDistance(k);
 308                 myNbExt++;
 309           if(myNbExt == NbExtMax)
 310           {
 311             break;
 312           }
 313                 // modified by NIZHNY-MKK  Thu Sep 18 14:46:43 2003.END
 314               }
 315       }
 316       if(myNbExt == NbExtMax)
 317       {
 318         break;
 319       }
 320       // try symmetric point
 321       myF.SetPoint(P); //To clear previous solutions
 322       U *= -1;
 323       MakePreciser(U, P, isMin, anOrtogSection);
 324       E = GetValue(U, myC);
 325       Pe = ProjectPnt (anOrtogSection, myDirection, E),
 326       V = gp_Vec(E,Pe) * gp_Vec(myDirection);
 327       UV(1) = U; UV(2) = V;
 328
 329       aFSR.Perform (myF,UV,UVinf,UVsup);
 330
 331       for (k=1 ; k <= myF.NbExt(); k++) {
 332         if(myF.SquareDistance(k) > Precision::Confusion()*Precision::Confusion())
 333         {
 334           //Additional checking solution: FSR sometimes is wrong
 335           //when starting point is far from solution.
 336           Standard_Real dist = Sqrt(myF.SquareDistance(k));
 337           math_Vector Vals(1, 2);
 338           const Extrema_POnSurf& PonS=myF.Point(k);
 339           Standard_Real u, v;
 340           PonS.Parameter(u, v);
 341           UV(1) = u;
 342           UV(2) = v;
 343           myF.Value(UV, Vals);
 344           gp_Vec du, dv;
 345           myS->D1(u, v, Pe, du, dv);
 346           Standard_Real mdu = du.Magnitude();
 347           Standard_Real mdv = dv.Magnitude();
 348           u = Abs(Vals(1));
 349           v = Abs(Vals(2));
 350           if(mdu > Precision::PConfusion())
 351           {
 352             if(u / dist / mdu > Precision::PConfusion())
 353             {
 354               continue;
 355             }
 356           }
 357           if(mdv > Precision::PConfusion())
 358           {
 359             if(v / dist / mdv > Precision::PConfusion())
 360             {
 361               continue;
 362             }
 363           }
 364
 365         }
 366               if (IsOriginalPnt(myF.Point(k).Value(), myPoint, myNbExt)) {
 367                 // modified by NIZHNY-MKK  Thu Sep 18 14:46:59 2003.BEGIN
 368                 //        myPoint[++myNbExt] = myF.Point(k);
 369                 //        myValue[myNbExt] = myF.Value(k);
 370                 myPoint[myNbExt] = myF.Point(k);
 371                 mySqDist[myNbExt] = myF.SquareDistance(k);
 372                 myNbExt++;
 373           if(myNbExt == NbExtMax)
 374           {
 375             break;
 376           }
 377                 // modified by NIZHNY-MKK  Thu Sep 18 14:47:04 2003.END
 378               }
 379       }
 380       if(myNbExt == NbExtMax)
 381       {
 382         break;
 383       }
 384     }
 385   }
 386   myDone = Standard_True;
 387   return;
 388 }
 389
 390 //=============================================================================
 391
 392 Standard_Boolean Extrema_ExtPExtS::IsDone () const { return myDone; }
 393 //=============================================================================
 394
 395 Standard_Integer Extrema_ExtPExtS::NbExt () const
 396 {
 397   if (!IsDone()) { StdFail_NotDone::Raise(); }
 398   if (myIsAnalyticallyComputable)
 399     return myNbExt;
 400   else
 401     return myExtPS.NbExt();
 402 }
 403 //=============================================================================
 404
 405 Standard_Real Extrema_ExtPExtS::SquareDistance (const Standard_Integer N) const
 406 {
 407   if (!IsDone()) { StdFail_NotDone::Raise(); }
 408   if ((N < 1) || (N > myNbExt)) { Standard_OutOfRange::Raise(); }
 409   if (myIsAnalyticallyComputable)
 410     // modified by NIZHNY-MKK  Thu Sep 18 14:48:39 2003.BEGIN
 411     //     return myValue[N];
 412     return mySqDist[N-1];
 413   // modified by NIZHNY-MKK  Thu Sep 18 14:48:42 2003.END
 414   else
 415     return myExtPS.SquareDistance(N);
 416 }
 417 //=============================================================================
 418
 419 const Extrema_POnSurf& Extrema_ExtPExtS::Point (const Standard_Integer N) const
 420 {
 421   if (!IsDone()) { StdFail_NotDone::Raise(); }
 422   if ((N < 1) || (N > myNbExt)) { Standard_OutOfRange::Raise(); }
 423   if (myIsAnalyticallyComputable) {
 424     // modified by NIZHNY-MKK  Thu Sep 18 14:47:40 2003.BEGIN
 425     //     return myPoint[N];
 426     return myPoint[N-1];
 427   }
 428   // modified by NIZHNY-MKK  Thu Sep 18 14:47:43 2003.END
 429   else
 430     return myExtPS.Point(N);
 431 }
 432 //=============================================================================
 433
 434
 435 static gp_Ax2 GetPosition (const Handle(Adaptor3d_HCurve)& C)
 436 {
 437   switch (C->GetType()) {
 438   case GeomAbs_Line: {
 439     gp_Lin L = C->Line();
 440     gp_Pln Pln = gp_Pln(L.Location(), L.Direction());
 441     //:abv 30.05.02: OCC  - use constructor instead of Set...s() to avoid exception
 442     gp_Ax2 Pos ( Pln.Location(), Pln.Position().Direction(), Pln.Position().XDirection() );
 443 //     Pos.SetAxis(Pln.XAxis());
 444 //     Pos.SetXDirection(Pln.YAxis().Direction());
 445 //     Pos.SetYDirection(Pln.Position().Direction());
 446     return Pos;
 447   }
 448   case GeomAbs_Circle:
 449     return  C->Circle().Position();
 450   case GeomAbs_Ellipse:
 451     return  C->Ellipse().Position();
 452   case GeomAbs_Hyperbola:
 453     return  C->Hyperbola().Position();
 454   case GeomAbs_Parabola:
 455     return  C->Parabola().Position();
 456   default:
 457     return gp_Ax2 ();
 458   }
 459 }
 460 //=============================================================================
 461
 462 static void PerformExtPElC (Extrema_ExtPElC& E,
 463                             const gp_Pnt& P,
 464                             const Handle(Adaptor3d_HCurve)& C,
 465                             const Standard_Real Tol)
 466 {
 467   switch (C->GetType()) {
 468   case GeomAbs_Hyperbola:
 469     E.Perform(P, C->Hyperbola(), Tol, -Precision::Infinite(),Precision::Infinite());
 470     return;
 471   case GeomAbs_Line:
 472     E.Perform(P, C->Line(), Tol, -Precision::Infinite(),Precision::Infinite());
 473     return;
 474   case GeomAbs_Circle:
 475     E.Perform(P, C->Circle(), Tol, 0.0, 2.0 * M_PI);
 476     return;
 477   case GeomAbs_Ellipse:
 478     E.Perform(P, C->Ellipse(), Tol, 0.0, 2.0 * M_PI);
 479     return;
 480   case GeomAbs_Parabola:
 481     E.Perform(P, C->Parabola(), Tol, -Precision::Infinite(),Precision::Infinite());
 482     return;
 483   default:
 484     return;
 485   }
 486 }
 487
 488 //=======================================================================
 489 //function : IsCaseAnalyticallyComputable
 490 //purpose  :
 491 //=======================================================================
 492
 493 static Standard_Boolean IsCaseAnalyticallyComputable
 494   (const GeomAbs_CurveType& theType,
 495    const gp_Ax2& theCurvePos,
 496    const gp_Dir& theSurfaceDirection)
 497 {
 498   // check type
 499   switch (theType) {
 500   case GeomAbs_Line:
 501   case GeomAbs_Circle:
 502   case GeomAbs_Ellipse:
 503   case GeomAbs_Hyperbola:
 504   case GeomAbs_Parabola:
 505     break;
 506   default:
 507     return  Standard_False;
 508   }
 509   // check if it is a plane
 510   if (Abs(theCurvePos.Direction() * theSurfaceDirection) <= gp::Resolution())
 511     return Standard_False;
 512   else
 513     return Standard_True;
 514 }
 515 //=======================================================================
 516 //function : GetValue
 517 //purpose  :
 518 //=======================================================================
 519
 520 static gp_Pnt GetValue(const Standard_Real U,
 521                        const Handle(Adaptor3d_HCurve)& C)
 522 {
 523   switch (C->GetType()) {
 524   case GeomAbs_Line:
 525     return ElCLib::Value(U, C->Line());
 526   case GeomAbs_Circle:
 527     return ElCLib::Value(U, C->Circle());
 528   case GeomAbs_Ellipse:
 529     return ElCLib::Value(U, C->Ellipse());
 530   case GeomAbs_Hyperbola:
 531     return ElCLib::Value(U, C->Hyperbola());
 532   case GeomAbs_Parabola:
 533     return ElCLib::Value(U, C->Parabola());
 534   default:
 535     return gp_Pnt();
 536   }
 537 }
 538 //=======================================================================
 539 //function : GetU
 540 //purpose  :
 541 //=======================================================================
 542 //#ifdef DEB
 543 //static Standard_Real GetU(const gp_Vec& vec,
 544 //                        const gp_Pnt& P,
 545 //                        const Handle(Adaptor3d_HCurve)& C)
 546 //{
 547 //  switch (C->GetType()) {
 548 //  case GeomAbs_Line:
 549 //    return ElCLib::Parameter(C->Line().Translated(vec), P);
 550 //  case GeomAbs_Circle:
 551 //    return ElCLib::Parameter(C->Circle().Translated(vec), P);
 552 //  case GeomAbs_Ellipse:
 553 //    return ElCLib::Parameter(C->Ellipse().Translated(vec), P);
 554 //  case GeomAbs_Hyperbola:
 555 //    return ElCLib::Parameter(C->Hyperbola().Translated(vec), P);
 556 //  case GeomAbs_Parabola:
 557 //    return ElCLib::Parameter(C->Parabola().Translated(vec), P);
 558 //  default:
 559 //    return 0;
 560 //  }
 561 //}
 562 //#endif