src/ApproxInt/ApproxInt_ImpPrmSvSurfaces.gxx

   1 // Created on: 1993-03-17
   2 // Created by: Laurent BUCHARD
   3 // Copyright (c) 1993-1999 Matra Datavision
   4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
   5 //
   6 // This file is part of Open CASCADE Technology software library.
   7 //
   8 // This library is free software; you can redistribute it and/or modify it under
   9 // the terms of the GNU Lesser General Public License version 2.1 as published
  10 // by the Free Software Foundation, with special exception defined in the file
  11 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
  12 // distribution for complete text of the license and disclaimer of any warranty.
  13 //
  14 // Alternatively, this file may be used under the terms of Open CASCADE
  15 // commercial license or contractual agreement.
  16
  17 #include <IntSurf_PntOn2S.hxx>
  18 #include <TColStd_Array1OfReal.hxx>
  19 #include <math_FunctionSetRoot.hxx>
  20 #include <StdFail_NotDone.hxx>
  21 #include <GeomAbs_SurfaceType.hxx>
  22 #include <Precision.hxx>
  23
  24 //=======================================================================
  25 //function : IsSingular
  26 //purpose  :    Returns TRUE if vectors theDU || theDV or if at least one
  27 //            of them has null-magnitude.
  28 //              theSqLinTol is square of linear tolerance.
  29 //              theAngTol is angular tolerance.
  30 //=======================================================================
  31 static Standard_Boolean IsSingular( const gp_Vec& theDU,
  32                                     const gp_Vec& theDV,
  33                                     const Standard_Real theSqLinTol,
  34                                     const Standard_Real theAngTol)
  35 {
  36   gp_Vec aDU(theDU), aDV(theDV);
  37
  38   const Standard_Real aSqMagnDU = aDU.SquareMagnitude(),
  39                       aSqMagnDV = aDV.SquareMagnitude();
  40
  41   if(aSqMagnDU < theSqLinTol)
  42     return Standard_True;
  43
  44   aDU.Divide(sqrt(aSqMagnDU));
  45
  46   if(aSqMagnDV < theSqLinTol)
  47     return Standard_True;
  48
  49   aDV.Divide(sqrt(aSqMagnDV));
  50
  51   //Here aDU and aDV vectors have magnitude 1.0.
  52
  53   if(aDU.Crossed(aDV).SquareMagnitude() < theAngTol*theAngTol)
  54     return Standard_True;
  55
  56   return Standard_False;
  57 }
  58
  59 //=======================================================================
  60 //function : SingularProcessing
  61 //purpose  :    Computes 2D-representation (in UV-coordinates) of
  62 //            theTg3D vector on the surface in case when
  63 //            theDU.Crossed(theDV).Magnitude() == 0.0. Stores result in
  64 //            theTg2D variable.
  65 //              theDU and theDV are vectors of 1st derivative
  66 //            (with respect to U and V variables correspondingly).
  67 //              If theIsTo3DTgCompute == TRUE then theTg3D has not been
  68 //            defined yet (it should be computed).
  69 //              theLinTol is SQUARE of the tolerance.
  70 //
  71 //Algorithm:
  72 //              Condition
  73 //                  Tg=theDU*theTg2D.X()+theDV*theTg2D.Y()
  74 //            has to be satisfied strictly.
  75 //              More over, vector Tg has to be NORMALYZED
  76 //            (if theIsTo3DTgCompute == TRUE then new computed vector will
  77 //            always have magnitude 1.0).
  78 //=======================================================================
  79 static Standard_Boolean SingularProcessing( const gp_Vec& theDU,
  80                                             const gp_Vec& theDV,
  81                                             const Standard_Boolean theIsTo3DTgCompute,
  82                                             const Standard_Real theLinTol,
  83                                             const Standard_Real theAngTol,
  84                                             gp_Vec& theTg3D,
  85                                             gp_Vec2d& theTg2D)
  86 {
  87   //Attention: @ \sin theAngTol \approx theAngTol @ (for cross-product)
  88
  89   //Really, vector theTg3D has to be normalyzed (if theIsTo3DTgCompute == FALSE).
  90   const Standard_Real aSQTan = theTg3D.SquareMagnitude();
  91
  92   const Standard_Real aSqMagnDU = theDU.SquareMagnitude(),
  93                       aSqMagnDV = theDV.SquareMagnitude();
  94
  95   //There are some reasons of singularity
  96
  97   //1.
  98   if((aSqMagnDU < theLinTol) && (aSqMagnDV < theLinTol))
  99   {
 100     //For future, this case can be processed as same as in case of
 101     //osculating surfaces (expanding in Taylor series). Here,
 102     //we return only.
 103
 104     return Standard_False;
 105   }
 106
 107   //2.
 108   if(aSqMagnDU < theLinTol)
 109   {
 110     //In this case, theTg3D vector will be parallel with theDV.
 111     //Its true direction shall be precised later (the algorithm is
 112     //based on array of Walking-points).
 113
 114     if(theIsTo3DTgCompute)
 115     {
 116       //theTg3D will be normalyzed. Its magnitude is
 117       const Standard_Real aTgMagn = 1.0;
 118
 119       const Standard_Real aNorm = sqrt(aSqMagnDV);
 120       theTg3D = theDV.Divided(aNorm);
 121       theTg2D.SetCoord(0.0, aTgMagn/aNorm);
 122     }
 123     else
 124     {
 125       //theTg3D is already defined.
 126       //Here we check only, if this tangent is parallel to theDV.
 127
 128       if(theDV.Crossed(theTg3D).SquareMagnitude() <
 129                     theAngTol*theAngTol*aSqMagnDV*aSQTan)
 130       {
 131         //theTg3D is parallel to theDV
 132
 133         //Use sign "+" if theTg3D and theDV are codirectional
 134         //and sign "-" if opposite
 135         const Standard_Real aDP = theTg3D.Dot(theDV);
 136         theTg2D.SetCoord(0.0, Sign(sqrt(aSQTan/aSqMagnDV), aDP));
 137       }
 138       else
 139       {
 140         //theTg3D is not parallel to theDV
 141         //It is abnormal
 142
 143         return Standard_False;
 144       }
 145     }
 146
 147     return Standard_True;
 148   }
 149
 150   //3.
 151   if(aSqMagnDV < theLinTol)
 152   {
 153     //In this case, theTg3D vector will be parallel with theDU.
 154     //Its true direction shall be precised later (the algorithm is
 155     //based on array of Walking-points).
 156
 157     if(theIsTo3DTgCompute)
 158     {
 159       //theTg3D will be normalyzed. Its magnitude is
 160       const Standard_Real aTgMagn = 1.0;
 161
 162       const Standard_Real aNorm = sqrt(aSqMagnDU);
 163       theTg3D = theDU.Divided(aNorm);
 164       theTg2D.SetCoord(aTgMagn/aNorm, 0.0);
 165     }
 166     else
 167     {
 168       //theTg3D is already defined.
 169       //Here we check only, if this tangent is parallel to theDU.
 170
 171       if(theDU.Crossed(theTg3D).SquareMagnitude() <
 172                     theAngTol*theAngTol*aSqMagnDU*aSQTan)
 173       {
 174         //theTg3D is parallel to theDU
 175
 176         //Use sign "+" if theTg3D and theDU are codirectional
 177         //and sign "-" if opposite
 178         const Standard_Real aDP = theTg3D.Dot(theDU);
 179         theTg2D.SetCoord(Sign(sqrt(aSQTan/aSqMagnDU), aDP), 0.0);
 180       }
 181       else
 182       {
 183         //theTg3D is not parallel to theDU
 184         //It is abnormal
 185
 186         return Standard_False;
 187       }
 188     }
 189
 190     return Standard_True;
 191   }
 192
 193   //4. If aSqMagnDU > 0.0 && aSqMagnDV > 0.0 but theDV || theDU.
 194
 195   const Standard_Real aLenU = sqrt(aSqMagnDU),
 196                       aLenV = sqrt(aSqMagnDV);
 197
 198   //aLenSum > 0.0 definitely
 199   const Standard_Real aLenSum = aLenU + aLenV;
 200
 201   if(theDV.Dot(theDU) > 0.0)
 202   {
 203     //Vectors theDV and theDU are codirectional.
 204
 205     if(theIsTo3DTgCompute)
 206     {
 207       theTg2D.SetCoord(1.0/aLenSum, 1.0/aLenSum);
 208       theTg3D = theDU*theTg2D.X() + theDV*theTg2D.Y();
 209     }
 210     else
 211     {
 212       //theTg3D is already defined.
 213       //Here we check only, if this tangent is parallel to theDU
 214       //(and theDV together).
 215
 216       if(theDU.Crossed(theTg3D).SquareMagnitude() <
 217                     theAngTol*theAngTol*aSqMagnDU*aSQTan)
 218       {
 219         //theTg3D is parallel to theDU
 220
 221         const Standard_Real aDP = theTg3D.Dot(theDU);
 222         const Standard_Real aLenTg = Sign(sqrt(aSQTan), aDP);
 223         theTg2D.SetCoord(aLenTg/aLenSum, aLenTg/aLenSum);
 224       }
 225       else
 226       {
 227         //theTg3D is not parallel to theDU
 228         //It is abnormal
 229
 230         return Standard_False;
 231       }
 232     }
 233   }
 234   else
 235   {
 236     //Vectors theDV and theDU are opposite.
 237
 238     if(theIsTo3DTgCompute)
 239     {
 240       //Here we chose theDU as direction of theTg3D.
 241       //True direction shall be precised later (the algorithm is
 242       //based on array of Walking-points).
 243
 244       theTg2D.SetCoord(1.0/aLenSum, -1.0/aLenSum);
 245       theTg3D = theDU*theTg2D.X() + theDV*theTg2D.Y();
 246     }
 247     else
 248     {
 249       //theTg3D is already defined.
 250       //Here we check only, if this tangent is parallel to theDU
 251       //(and theDV together).
 252
 253       if(theDU.Crossed(theTg3D).SquareMagnitude() <
 254                     theAngTol*theAngTol*aSqMagnDU*aSQTan)
 255       {
 256         //theTg3D is parallel to theDU
 257
 258         const Standard_Real aDP = theTg3D.Dot(theDU);
 259         const Standard_Real aLenTg = Sign(sqrt(aSQTan), aDP);
 260         theTg2D.SetCoord(aLenTg/aLenSum, -aLenTg/aLenSum);
 261       }
 262       else
 263       {
 264         //theTg3D is not parallel to theDU
 265         //It is abnormal
 266
 267         return Standard_False;
 268       }
 269     }
 270   }
 271
 272   return Standard_True;
 273 }
 274
 275 //=======================================================================
 276 //function : NonSingularProcessing
 277 //purpose  :    Computes 2D-representation (in UV-coordinates) of
 278 //            theTg3D vector on the surface in case when
 279 //            theDU.Crossed(theDV).Magnitude() > 0.0. Stores result in
 280 //            theTg2D variable.
 281 //              theDU and theDV are vectors of 1st derivative
 282 //            (with respect to U and V variables correspondingly).
 283 //              theLinTol is SQUARE of the tolerance.
 284 //
 285 //Algorithm:
 286 //              Condition
 287 //                  Tg=theDU*theTg2D.X()+theDV*theTg2D.Y()
 288 //            has to be satisfied strictly.
 289 //              More over, vector Tg has always to be NORMALYZED.
 290 //=======================================================================
 291 static Standard_Boolean NonSingularProcessing(const gp_Vec& theDU,
 292                                               const gp_Vec& theDV,
 293                                               const gp_Vec& theTg3D,
 294                                               const Standard_Real theLinTol,
 295                                               const Standard_Real theAngTol,
 296                                               gp_Vec2d& theTg2D)
 297 {
 298   const gp_Vec aNormal = theDU.Crossed(theDV);
 299   const Standard_Real aSQMagn = aNormal.SquareMagnitude();
 300
 301   if(IsSingular(theDU, theDV, theLinTol, theAngTol))
 302   {
 303     gp_Vec aTg(theTg3D);
 304     return
 305       SingularProcessing(theDU, theDV, Standard_False,
 306                           theLinTol, theAngTol, aTg, theTg2D);
 307   }
 308
 309   //If @\vec{T}=\vec{A}*U+\vec{B}*V@ then
 310
 311   //  \left\{\begin{matrix}
 312   //  \vec{A} \times \vec{T} = (\vec{A} \times \vec{B})*V
 313   //  \vec{B} \times \vec{T} = (\vec{B} \times \vec{A})*U
 314   //  \end{matrix}\right.
 315
 316   //From here, values of U and V can be found very easily
 317   //(if @\left \| \vec{A} \times \vec{B} \right \| > 0.0 @,
 318   //else it is singular case).
 319
 320   const gp_Vec aTgU(theTg3D.Crossed(theDU)), aTgV(theTg3D.Crossed(theDV));
 321   const Standard_Real aDeltaU = aTgV.SquareMagnitude()/aSQMagn;
 322   const Standard_Real aDeltaV = aTgU.SquareMagnitude()/aSQMagn;
 323
 324   theTg2D.SetCoord(Sign(sqrt(aDeltaU), aTgV.Dot(aNormal)), -Sign(sqrt(aDeltaV), aTgU.Dot(aNormal)));
 325
 326   return Standard_True;
 327 }
 328
 329 //--------------------------------------------------------------------------------
 330 ApproxInt_ImpPrmSvSurfaces::ApproxInt_ImpPrmSvSurfaces( const TheISurface& ISurf
 331                                                        ,const ThePSurface& PSurf):
 332        MyIsTangent(Standard_False),
 333        MyHasBeenComputed(Standard_False),
 334        MyIsTangentbis(Standard_False),
 335        MyHasBeenComputedbis(Standard_False),
 336        MyImplicitFirst(Standard_True),
 337        MyZerImpFunc(PSurf,ISurf)
 338 {
 339 }
 340 //--------------------------------------------------------------------------------
 341 ApproxInt_ImpPrmSvSurfaces::ApproxInt_ImpPrmSvSurfaces( const ThePSurface& PSurf
 342                                                        ,const TheISurface& ISurf):
 343        MyIsTangent(Standard_False),
 344        MyHasBeenComputed(Standard_False),
 345        MyIsTangentbis(Standard_False),
 346        MyHasBeenComputedbis(Standard_False),
 347        MyImplicitFirst(Standard_False),
 348        MyZerImpFunc(PSurf,ISurf)
 349 {
 350 }
 351 //--------------------------------------------------------------------------------
 352 void ApproxInt_ImpPrmSvSurfaces::Pnt(const Standard_Real u1,
 353                                      const Standard_Real v1,
 354                                      const Standard_Real u2,
 355                                      const Standard_Real v2,
 356                                      gp_Pnt& P) {
 357   gp_Pnt aP;
 358   gp_Vec aT;
 359   gp_Vec2d aTS1,aTS2;
 360   Standard_Real tu1=u1;
 361   Standard_Real tu2=u2;
 362   Standard_Real tv1=v1;
 363   Standard_Real tv2=v2;
 364   this->Compute(tu1,tv1,tu2,tv2,aP,aT,aTS1,aTS2);
 365   P=MyPnt;
 366 }
 367 //--------------------------------------------------------------------------------
 368 Standard_Boolean ApproxInt_ImpPrmSvSurfaces::Tangency(const Standard_Real u1,
 369                                                       const Standard_Real v1,
 370                                                       const Standard_Real u2,
 371                                                       const Standard_Real v2,
 372                                                       gp_Vec& T) {
 373   gp_Pnt aP;
 374   gp_Vec aT;
 375   gp_Vec2d aTS1,aTS2;
 376   Standard_Real tu1=u1;
 377   Standard_Real tu2=u2;
 378   Standard_Real tv1=v1;
 379   Standard_Real tv2=v2;
 380   Standard_Boolean t=this->Compute(tu1,tv1,tu2,tv2,aP,aT,aTS1,aTS2);
 381   T=MyTg;
 382   return(t);
 383 }
 384 //--------------------------------------------------------------------------------
 385 Standard_Boolean ApproxInt_ImpPrmSvSurfaces::TangencyOnSurf1(const Standard_Real u1,
 386                                                              const Standard_Real v1,
 387                                                              const Standard_Real u2,
 388                                                              const Standard_Real v2,
 389                                                              gp_Vec2d& T) {
 390   gp_Pnt aP;
 391   gp_Vec aT;
 392   gp_Vec2d aTS1,aTS2;
 393   Standard_Real tu1=u1;
 394   Standard_Real tu2=u2;
 395   Standard_Real tv1=v1;
 396   Standard_Real tv2=v2;
 397   Standard_Boolean t=this->Compute(tu1,tv1,tu2,tv2,aP,aT,aTS1,aTS2);
 398   T=MyTguv1;
 399   return(t);
 400 }
 401 //--------------------------------------------------------------------------------
 402 Standard_Boolean ApproxInt_ImpPrmSvSurfaces::TangencyOnSurf2(const Standard_Real u1,
 403                                                              const Standard_Real v1,
 404                                                              const Standard_Real u2,
 405                                                              const Standard_Real v2,
 406                                                              gp_Vec2d& T) {
 407   gp_Pnt aP;
 408   gp_Vec aT;
 409   gp_Vec2d aTS1,aTS2;
 410   Standard_Real tu1=u1;
 411   Standard_Real tu2=u2;
 412   Standard_Real tv1=v1;
 413   Standard_Real tv2=v2;
 414   Standard_Boolean t=this->Compute(tu1,tv1,tu2,tv2,aP,aT,aTS1,aTS2);
 415   T=MyTguv2;
 416   return(t);
 417 }
 418
 419 //=======================================================================
 420 //function : Compute
 421 //purpose  :    Computes point on curve, 3D and 2D-tangents of a curve and
 422 //            parameters on the surfaces.
 423 //=======================================================================
 424 Standard_Boolean ApproxInt_ImpPrmSvSurfaces::Compute( Standard_Real& u1,
 425                                                       Standard_Real& v1,
 426                                                       Standard_Real& u2,
 427                                                       Standard_Real& v2,
 428                                                       gp_Pnt& P,
 429                                                       gp_Vec& Tg,
 430                                                       gp_Vec2d& Tguv1,
 431                                                       gp_Vec2d& Tguv2)
 432 {
 433   const IntSurf_Quadric&  aQSurf = MyZerImpFunc.ISurface();
 434   const ThePSurface&      aPSurf = MyZerImpFunc.PSurface();
 435   gp_Vec2d& aQuadTg = MyImplicitFirst ? Tguv1 : Tguv2;
 436   gp_Vec2d& aPrmTg = MyImplicitFirst ? Tguv2 : Tguv1;
 437
 438   //for square
 439   const Standard_Real aNullValue =  Precision::Approximation()*
 440                                     Precision::Approximation(),
 441                       anAngTol = Precision::Angular();
 442
 443   Standard_Real tu1=u1;
 444   Standard_Real tu2=u2;
 445   Standard_Real tv1=v1;
 446   Standard_Real tv2=v2;
 447
 448   if(MyHasBeenComputed) {
 449     if(  (MyParOnS1.X()==u1)&&(MyParOnS1.Y()==v1)
 450        &&(MyParOnS2.X()==u2)&&(MyParOnS2.Y()==v2)) {
 451       return(MyIsTangent);
 452     }
 453     else if(MyHasBeenComputedbis == Standard_False) {
 454       MyTgbis         = MyTg;
 455       MyTguv1bis      = MyTguv1;
 456       MyTguv2bis      = MyTguv2;
 457       MyPntbis        = MyPnt;
 458       MyParOnS1bis    = MyParOnS1;
 459       MyParOnS2bis    = MyParOnS2;
 460       MyIsTangentbis  = MyIsTangent;
 461       MyHasBeenComputedbis = MyHasBeenComputed;
 462     }
 463   }
 464
 465   if(MyHasBeenComputedbis) {
 466     if(  (MyParOnS1bis.X()==u1)&&(MyParOnS1bis.Y()==v1)
 467        &&(MyParOnS2bis.X()==u2)&&(MyParOnS2bis.Y()==v2)) {
 468
 469       gp_Vec            TV(MyTg);
 470       gp_Vec2d          TV1(MyTguv1);
 471       gp_Vec2d          TV2(MyTguv2);
 472       gp_Pnt            TP(MyPnt);
 473       gp_Pnt2d          TP1(MyParOnS1);
 474       gp_Pnt2d          TP2(MyParOnS2);
 475       Standard_Boolean  TB=MyIsTangent;
 476
 477       MyTg        = MyTgbis;
 478       MyTguv1     = MyTguv1bis;
 479       MyTguv2     = MyTguv2bis;
 480       MyPnt       = MyPntbis;
 481       MyParOnS1   = MyParOnS1bis;
 482       MyParOnS2   = MyParOnS2bis;
 483       MyIsTangent = MyIsTangentbis;
 484
 485       MyTgbis         = TV;
 486       MyTguv1bis      = TV1;
 487       MyTguv2bis      = TV2;
 488       MyPntbis        = TP;
 489       MyParOnS1bis    = TP1;
 490       MyParOnS2bis    = TP2;
 491       MyIsTangentbis  = TB;
 492
 493       return(MyIsTangent);
 494     }
 495   }
 496
 497   math_Vector X(1,2);
 498   math_Vector BornInf(1,2), BornSup(1,2), Tolerance(1,2);
 499   //--- ThePSurfaceTool::GetResolution(aPSurf,Tolerance(1),Tolerance(2));
 500   Tolerance(1) = 1.0e-8; Tolerance(2) = 1.0e-8;
 501   Standard_Real binfu,bsupu,binfv,bsupv;
 502   binfu = ThePSurfaceTool::FirstUParameter(aPSurf);
 503   binfv = ThePSurfaceTool::FirstVParameter(aPSurf);
 504   bsupu = ThePSurfaceTool::LastUParameter(aPSurf);
 505   bsupv = ThePSurfaceTool::LastVParameter(aPSurf);
 506   BornInf(1) = binfu; BornSup(1) = bsupu;
 507   BornInf(2) = binfv; BornSup(2) = bsupv;
 508   Standard_Real TranslationU = 0., TranslationV = 0.;
 509
 510   if (!FillInitialVectorOfSolution(u1, v1, u2, v2,
 511                                    binfu, bsupu, binfv, bsupv,
 512                                    X,
 513                                    TranslationU, TranslationV))
 514   {
 515     MyIsTangent=MyIsTangentbis=Standard_False;
 516     MyHasBeenComputed = MyHasBeenComputedbis = Standard_False;
 517     return(Standard_False);
 518   }
 519
 520   Standard_Real PourTesterU = X(1);
 521   Standard_Real PourTesterV = X(2);
 522
 523   math_FunctionSetRoot  Rsnld(MyZerImpFunc);
 524   Rsnld.SetTolerance(Tolerance);
 525   Rsnld.Perform(MyZerImpFunc,X,BornInf,BornSup);
 526   if(Rsnld.IsDone()) {
 527     MyHasBeenComputed = Standard_True;
 528     Rsnld.Root(X);
 529
 530     Standard_Real DistAvantApresU = Abs(PourTesterU-X(1));
 531     Standard_Real DistAvantApresV = Abs(PourTesterV-X(2));
 532
 533     MyPnt = P = ThePSurfaceTool::Value(aPSurf, X(1), X(2));
 534
 535     if( (DistAvantApresV <= 0.001 ) &&
 536         (DistAvantApresU <= 0.001 ))
 537     {
 538       gp_Vec aD1uPrm,aD1vPrm;
 539       gp_Vec aD1uQuad,aD1vQuad;
 540
 541       if(MyImplicitFirst)
 542       {
 543         u2 = X(1)-TranslationU;
 544         v2 = X(2)-TranslationV;
 545
 546         if(aQSurf.TypeQuadric() != GeomAbs_Plane)
 547         {
 548           while(u1-tu1>M_PI)  u1-=M_PI+M_PI;
 549           while(tu1-u1>M_PI)  u1+=M_PI+M_PI;
 550         }
 551
 552         MyParOnS1.SetCoord(tu1,tv1);
 553         MyParOnS2.SetCoord(tu2,tv2);
 554
 555         gp_Pnt aP2;
 556
 557         ThePSurfaceTool::D1(aPSurf, X(1), X(2), P, aD1uPrm, aD1vPrm);
 558         aQSurf.D1(u1,v1, aP2, aD1uQuad, aD1vQuad);
 559
 560         //Middle-point of P-P2 segment
 561         P.BaryCenter(1.0, aP2, 1.0);
 562       }
 563       else
 564       {
 565         u1 = X(1)-TranslationU;
 566         v1 = X(2)-TranslationV;
 567         //aQSurf.Parameters(P, u2, v2);
 568         if(aQSurf.TypeQuadric() != GeomAbs_Plane)
 569         {
 570           while(u2-tu2>M_PI)  u2-=M_PI+M_PI;
 571           while(tu2-u2>M_PI)  u2+=M_PI+M_PI;
 572         }
 573
 574         MyParOnS1.SetCoord(tu1,tv1);
 575         MyParOnS2.SetCoord(tu2,tu2);
 576
 577         gp_Pnt aP2;
 578         ThePSurfaceTool::D1(aPSurf, X(1), X(2), P, aD1uPrm, aD1vPrm);
 579
 580         aQSurf.D1(u2, v2, aP2, aD1uQuad, aD1vQuad);
 581
 582         //Middle-point of P-P2 segment
 583         P.BaryCenter(1.0, aP2, 1.0);
 584       }
 585
 586       MyPnt = P;
 587
 588       //Normals to the surfaces
 589       gp_Vec  aNormalPrm(aD1uPrm.Crossed(aD1vPrm)),
 590               aNormalImp(aQSurf.Normale(MyPnt));
 591
 592       const Standard_Real aSQMagnPrm = aNormalPrm.SquareMagnitude(),
 593                           aSQMagnImp = aNormalImp.SquareMagnitude();
 594
 595       Standard_Boolean  isPrmSingular = Standard_False,
 596                         isImpSingular = Standard_False;
 597
 598       if(IsSingular(aD1uPrm, aD1vPrm, aNullValue, anAngTol))
 599       {
 600         isPrmSingular = Standard_True;
 601         if(!SingularProcessing(aD1uPrm, aD1vPrm, Standard_True,
 602                                 aNullValue, anAngTol, Tg, aPrmTg))
 603         {
 604           MyIsTangent = Standard_False;
 605           MyHasBeenComputed = MyHasBeenComputedbis = Standard_False;
 606           return Standard_False;
 607         }
 608
 609         MyTg = Tg;
 610       }
 611       else
 612       {
 613         aNormalPrm.Divide(sqrt(aSQMagnPrm));
 614       }
 615
 616       //Analogicaly for implicit surface
 617       if(aSQMagnImp < aNullValue)
 618       {
 619         isImpSingular = Standard_True;
 620
 621         if(!SingularProcessing(aD1uQuad, aD1vQuad, !isPrmSingular,
 622                                   aNullValue, anAngTol, Tg, aQuadTg))
 623         {
 624           MyIsTangent = Standard_False;
 625           MyHasBeenComputed = MyHasBeenComputedbis = Standard_False;
 626           return Standard_False;
 627         }
 628
 629         MyTg = Tg;
 630       }
 631       else
 632       {
 633         aNormalImp.Divide(sqrt(aSQMagnImp));
 634       }
 635
 636       if(isImpSingular && isPrmSingular)
 637       {
 638         //All is OK. All abnormal cases were processed above.
 639
 640         MyTguv1 = Tguv1;
 641         MyTguv2 = Tguv2;
 642
 643         MyIsTangent=Standard_True;
 644         return MyIsTangent;
 645       }
 646       else if(!(isImpSingular || isPrmSingular))
 647       {
 648         //Processing pure non-singular case
 649         //(3D- and 2D-tangents are still not defined)
 650
 651         //Ask to pay attention to the fact that here
 652         //aNormalImp and aNormalPrm are normalyzed.
 653         //Therefore, @ \left \| \vec{Tg} \right \| = 0.0 @
 654         //if and only if (aNormalImp || aNormalPrm).
 655         Tg = aNormalImp.Crossed(aNormalPrm);
 656       }
 657
 658       const Standard_Real aSQMagnTg = Tg.SquareMagnitude();
 659
 660       if(aSQMagnTg < aNullValue)
 661       {
 662         MyIsTangent = Standard_False;
 663         MyHasBeenComputed = MyHasBeenComputedbis = Standard_False;
 664         return Standard_False;
 665       }
 666
 667       //Normalyze Tg vector
 668       Tg.Divide(sqrt(aSQMagnTg));
 669       MyTg = Tg;
 670
 671       if(!isPrmSingular)
 672       {
 673         //If isPrmSingular==TRUE then aPrmTg has already been computed.
 674
 675         if(!NonSingularProcessing(aD1uPrm, aD1vPrm, Tg, aNullValue, anAngTol, aPrmTg))
 676         {
 677           MyIsTangent = Standard_False;
 678           MyHasBeenComputed = MyHasBeenComputedbis = Standard_False;
 679           return Standard_False;
 680         }
 681       }
 682
 683       if(!isImpSingular)
 684       {
 685         //If isImpSingular==TRUE then aQuadTg has already been computed.
 686
 687         if(!NonSingularProcessing(aD1uQuad, aD1vQuad, Tg, aNullValue, anAngTol, aQuadTg))
 688         {
 689           MyIsTangent = Standard_False;
 690           MyHasBeenComputed = MyHasBeenComputedbis = Standard_False;
 691           return Standard_False;
 692         }
 693       }
 694
 695       MyTguv1 = Tguv1;
 696       MyTguv2 = Tguv2;
 697
 698       MyIsTangent=Standard_True;
 699
 700 #ifdef OCCT_DEBUG
 701       //cout << "+++++++++++++++++  ApproxInt_ImpPrmSvSurfaces::Compute(...)  ++++++++++" << endl;
 702       //printf( "P2d_1(%+10.20f, %+10.20f); P2d_2(%+10.20f, %+10.20f)\n"
 703       //        "P(%+10.20f, %+10.20f, %+10.20f);\n"
 704       //        "Tg = {%+10.20f, %+10.20f, %+10.20f};\n"
 705       //        "Tguv1 = {%+10.20f, %+10.20f};\n"
 706       //        "Tguv2 = {%+10.20f, %+10.20f}\n",
 707       //        u1, v1, u2, v2,
 708       //        P.X(), P.Y(), P.Z(),
 709       //        Tg.X(), Tg.Y(), Tg.Z(),
 710       //        Tguv1.X(), Tguv1.Y(), Tguv2.X(), Tguv2.Y());
 711       //cout << "-----------------------------------------------------------------------" << endl;
 712 #endif
 713
 714       return Standard_True;
 715     }
 716     else {
 717       //-- cout<<" ApproxInt_ImpImpSvSurfaces.gxx : Distance apres recadrage Trop Grande "<<endl;
 718
 719       MyIsTangent=Standard_False;
 720       MyHasBeenComputed = MyHasBeenComputedbis = Standard_False;
 721       return Standard_False;
 722     }
 723   }
 724   else {
 725     MyIsTangent = Standard_False;
 726     MyHasBeenComputed = MyHasBeenComputedbis = Standard_False;
 727     return Standard_False;
 728   }
 729 }
 730
 731 //=======================================================================
 732 //function : SeekPoint
 733 //purpose  :    Computes point on curve and
 734 //            parameters on the surfaces.
 735 //=======================================================================
 736 Standard_Boolean ApproxInt_ImpPrmSvSurfaces::SeekPoint(const Standard_Real u1,
 737                                                        const Standard_Real v1,
 738                                                        const Standard_Real u2,
 739                                                        const Standard_Real v2,
 740                                                        IntSurf_PntOn2S& Point) {
 741   gp_Pnt aP;
 742   gp_Vec aT;
 743   gp_Vec2d aTS1,aTS2;
 744   const IntSurf_Quadric&  aQSurf = MyZerImpFunc.ISurface();
 745   const ThePSurface&      aPSurf = MyZerImpFunc.PSurface();
 746
 747   math_Vector X(1,2);
 748   math_Vector BornInf(1,2), BornSup(1,2), Tolerance(1,2);
 749   //--- ThePSurfaceTool::GetResolution(aPSurf,Tolerance(1),Tolerance(2));
 750   Tolerance(1) = 1.0e-8; Tolerance(2) = 1.0e-8;
 751   Standard_Real binfu,bsupu,binfv,bsupv;
 752   binfu = ThePSurfaceTool::FirstUParameter(aPSurf);
 753   binfv = ThePSurfaceTool::FirstVParameter(aPSurf);
 754   bsupu = ThePSurfaceTool::LastUParameter(aPSurf);
 755   bsupv = ThePSurfaceTool::LastVParameter(aPSurf);
 756   BornInf(1) = binfu; BornSup(1) = bsupu;
 757   BornInf(2) = binfv; BornSup(2) = bsupv;
 758   Standard_Real TranslationU = 0., TranslationV = 0.;
 759
 760   if (!FillInitialVectorOfSolution(u1, v1, u2, v2,
 761                                    binfu, bsupu, binfv, bsupv,
 762                                    X,
 763                                    TranslationU, TranslationV))
 764     return Standard_False;
 765
 766   Standard_Real NewU1, NewV1, NewU2, NewV2;
 767
 768   math_FunctionSetRoot  Rsnld(MyZerImpFunc);
 769   Rsnld.SetTolerance(Tolerance);
 770   Rsnld.Perform(MyZerImpFunc,X,BornInf,BornSup);
 771   if(Rsnld.IsDone()) {
 772     MyHasBeenComputed = Standard_True;
 773     Rsnld.Root(X);
 774
 775     MyPnt = ThePSurfaceTool::Value(aPSurf, X(1), X(2));
 776
 777     if(MyImplicitFirst)
 778     {
 779       NewU2 = X(1)-TranslationU;
 780       NewV2 = X(2)-TranslationV;
 781
 782       aQSurf.Parameters(MyPnt, NewU1, NewV1);
 783       //adjust U
 784       if (aQSurf.TypeQuadric() != GeomAbs_Plane)
 785       {
 786         Standard_Real sign = (NewU1 > u1)? -1 : 1;
 787         while (Abs(u1 - NewU1) > M_PI)
 788           NewU1 += sign*(M_PI+M_PI);
 789       }
 790     }
 791     else
 792     {
 793       NewU1 = X(1)-TranslationU;
 794       NewV1 = X(2)-TranslationV;
 795
 796       aQSurf.Parameters(MyPnt, NewU2, NewV2);
 797       //adjust U
 798       if (aQSurf.TypeQuadric() != GeomAbs_Plane)
 799       {
 800         Standard_Real sign = (NewU2 > u2)? -1 : 1;
 801         while (Abs(u2 - NewU2) > M_PI)
 802           NewU2 += sign*(M_PI+M_PI);
 803       }
 804     }
 805   }
 806   else
 807     return Standard_False;
 808
 809   Point.SetValue(MyPnt, NewU1, NewV1, NewU2, NewV2);
 810   return Standard_True;
 811 }
 812 //--------------------------------------------------------------------------------
 813
 814 Standard_Boolean
 815 ApproxInt_ImpPrmSvSurfaces::FillInitialVectorOfSolution(const Standard_Real u1,
 816                                                         const Standard_Real v1,
 817                                                         const Standard_Real u2,
 818                                                         const Standard_Real v2,
 819                                                         const Standard_Real binfu,
 820                                                         const Standard_Real bsupu,
 821                                                         const Standard_Real binfv,
 822                                                         const Standard_Real bsupv,
 823                                                         math_Vector& X,
 824                                                         Standard_Real& TranslationU,
 825                                                         Standard_Real& TranslationV)
 826 {
 827   const ThePSurface&      aPSurf = MyZerImpFunc.PSurface();
 828
 829   math_Vector F(1,1);
 830
 831   TranslationU = 0.0;
 832   TranslationV = 0.0;
 833
 834   if(MyImplicitFirst) {
 835     if(u2<binfu-0.0000000001) {
 836       if(ThePSurfaceTool::IsUPeriodic(aPSurf)) {
 837         Standard_Real d = ThePSurfaceTool::UPeriod(aPSurf);
 838         do {  TranslationU+=d; } while(u2+TranslationU < binfu);
 839       }
 840       else
 841         return(Standard_False);
 842     }
 843     else if(u2>bsupu+0.0000000001) {
 844       if(ThePSurfaceTool::IsUPeriodic(aPSurf)) {
 845         Standard_Real d = ThePSurfaceTool::UPeriod(aPSurf);
 846         do { TranslationU-=d; } while(u2+TranslationU > bsupu);
 847       }
 848       else
 849         return(Standard_False);
 850     }
 851     if(v2<binfv-0.0000000001) {
 852       if(ThePSurfaceTool::IsVPeriodic(aPSurf)) {
 853         Standard_Real d = ThePSurfaceTool::VPeriod(aPSurf);
 854         do { TranslationV+=d; } while(v2+TranslationV < binfv);
 855       }
 856       else
 857         return(Standard_False);
 858     }
 859     else if(v2>bsupv+0.0000000001) {
 860       if(ThePSurfaceTool::IsVPeriodic(aPSurf)) {
 861         Standard_Real d = ThePSurfaceTool::VPeriod(aPSurf);
 862         do { TranslationV-=d; } while(v2+TranslationV > bsupv);
 863       }
 864       else
 865         return(Standard_False);
 866     }
 867     X(1) = u2+TranslationU;
 868     X(2) = v2+TranslationV;
 869   }
 870   else {
 871     if(u1<binfu-0.0000000001) {
 872       if(ThePSurfaceTool::IsUPeriodic(aPSurf)) {
 873         Standard_Real d = ThePSurfaceTool::UPeriod(aPSurf);
 874         do {  TranslationU+=d;  } while(u1+TranslationU < binfu);
 875       }
 876       else
 877         return(Standard_False);
 878     }
 879     else if(u1>bsupu+0.0000000001) {
 880       if(ThePSurfaceTool::IsUPeriodic(aPSurf)) {
 881         Standard_Real d = ThePSurfaceTool::UPeriod(aPSurf);
 882         do { TranslationU-=d; } while(u1+TranslationU > bsupu);
 883       }
 884       else
 885         return(Standard_False);
 886     }
 887     if(v1<binfv-0.0000000001) {
 888       if(ThePSurfaceTool::IsVPeriodic(aPSurf)) {
 889         Standard_Real d = ThePSurfaceTool::VPeriod(aPSurf);
 890         do { TranslationV+=d; } while(v1+TranslationV < binfv);
 891       }
 892       else
 893         return(Standard_False);
 894     }
 895     else if(v1>bsupv+0.0000000001) {
 896       if(ThePSurfaceTool::IsVPeriodic(aPSurf)) {
 897         Standard_Real d = ThePSurfaceTool::VPeriod(aPSurf);
 898         do { TranslationV-=d; } while(v1+TranslationV > bsupv);
 899       }
 900       else
 901         return(Standard_False);
 902     }
 903     X(1) = u1+TranslationU;
 904     X(2) = v1+TranslationV;
 905   }
 906
 907   //----------------------------------------------------
 908   //Make a small step from boundaries in order to avoid
 909   //finding "outboundaried" solution (Rsnld -> NotDone).
 910   if(X(1)-0.0000000001 <= binfu) X(1)=X(1)+0.0000001;
 911   if(X(1)+0.0000000001 >= bsupu) X(1)=X(1)-0.0000001;
 912   if(X(2)-0.0000000001 <= binfv) X(2)=X(2)+0.0000001;
 913   if(X(2)+0.0000000001 >= bsupv) X(2)=X(2)-0.0000001;
 914
 915   return Standard_True;
 916 }