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//======================================================================= 31static 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 NORMALIZED 76// (if theIsTo3DTgCompute == TRUE then new computed vector will 77// always have magnitude 1.0). 78//======================================================================= 79static 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 normalized (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 normalized. 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 normalized. 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 NORMALIZED. 290//======================================================================= 291static 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//-------------------------------------------------------------------------------- 330ApproxInt_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 SetUseSolver(Standard_True); 340} 341//-------------------------------------------------------------------------------- 342ApproxInt_ImpPrmSvSurfaces::ApproxInt_ImpPrmSvSurfaces( const ThePSurface& PSurf 343 ,const TheISurface& ISurf): 344 MyIsTangent(Standard_False), 345 MyHasBeenComputed(Standard_False), 346 MyIsTangentbis(Standard_False), 347 MyHasBeenComputedbis(Standard_False), 348 MyImplicitFirst(Standard_False), 349 MyZerImpFunc(PSurf,ISurf) 350{ 351 SetUseSolver(Standard_True); 352} 353//-------------------------------------------------------------------------------- 354void ApproxInt_ImpPrmSvSurfaces::Pnt(const Standard_Real u1, 355 const Standard_Real v1, 356 const Standard_Real u2, 357 const Standard_Real v2, 358 gp_Pnt& P) { 359 gp_Pnt aP; 360 gp_Vec aT; 361 gp_Vec2d aTS1,aTS2; 362 Standard_Real tu1=u1; 363 Standard_Real tu2=u2; 364 Standard_Real tv1=v1; 365 Standard_Real tv2=v2; 366 this->Compute(tu1,tv1,tu2,tv2,aP,aT,aTS1,aTS2); 367 P=MyPnt; 368} 369//-------------------------------------------------------------------------------- 370Standard_Boolean ApproxInt_ImpPrmSvSurfaces::Tangency(const Standard_Real u1, 371 const Standard_Real v1, 372 const Standard_Real u2, 373 const Standard_Real v2, 374 gp_Vec& T) { 375 gp_Pnt aP; 376 gp_Vec aT; 377 gp_Vec2d aTS1,aTS2; 378 Standard_Real tu1=u1; 379 Standard_Real tu2=u2; 380 Standard_Real tv1=v1; 381 Standard_Real tv2=v2; 382 Standard_Boolean t=this->Compute(tu1,tv1,tu2,tv2,aP,aT,aTS1,aTS2); 383 T=MyTg; 384 return(t); 385} 386//-------------------------------------------------------------------------------- 387Standard_Boolean ApproxInt_ImpPrmSvSurfaces::TangencyOnSurf1(const Standard_Real u1, 388 const Standard_Real v1, 389 const Standard_Real u2, 390 const Standard_Real v2, 391 gp_Vec2d& T) { 392 gp_Pnt aP; 393 gp_Vec aT; 394 gp_Vec2d aTS1,aTS2; 395 Standard_Real tu1=u1; 396 Standard_Real tu2=u2; 397 Standard_Real tv1=v1; 398 Standard_Real tv2=v2; 399 Standard_Boolean t=this->Compute(tu1,tv1,tu2,tv2,aP,aT,aTS1,aTS2); 400 T=MyTguv1; 401 return(t); 402} 403//-------------------------------------------------------------------------------- 404Standard_Boolean ApproxInt_ImpPrmSvSurfaces::TangencyOnSurf2(const Standard_Real u1, 405 const Standard_Real v1, 406 const Standard_Real u2, 407 const Standard_Real v2, 408 gp_Vec2d& T) { 409 gp_Pnt aP; 410 gp_Vec aT; 411 gp_Vec2d aTS1,aTS2; 412 Standard_Real tu1=u1; 413 Standard_Real tu2=u2; 414 Standard_Real tv1=v1; 415 Standard_Real tv2=v2; 416 Standard_Boolean t=this->Compute(tu1,tv1,tu2,tv2,aP,aT,aTS1,aTS2); 417 T=MyTguv2; 418 return(t); 419} 420 421//======================================================================= 422//function : Compute 423//purpose : Computes point on curve, 3D and 2D-tangents of a curve and 424// parameters on the surfaces. 425//======================================================================= 426Standard_Boolean ApproxInt_ImpPrmSvSurfaces::Compute( Standard_Real& u1, 427 Standard_Real& v1, 428 Standard_Real& u2, 429 Standard_Real& v2, 430 gp_Pnt& P, 431 gp_Vec& Tg, 432 gp_Vec2d& Tguv1, 433 gp_Vec2d& Tguv2) 434{ 435 const IntSurf_Quadric& aQSurf = MyZerImpFunc.ISurface(); 436 const ThePSurface& aPSurf = MyZerImpFunc.PSurface(); 437 gp_Vec2d& aQuadTg = MyImplicitFirst ? Tguv1 : Tguv2; 438 gp_Vec2d& aPrmTg = MyImplicitFirst ? Tguv2 : Tguv1; 439 440 //for square 441 const Standard_Real aNullValue = Precision::Approximation()* 442 Precision::Approximation(), 443 anAngTol = Precision::Angular(); 444 445 Standard_Real tu1=u1; 446 Standard_Real tu2=u2; 447 Standard_Real tv1=v1; 448 Standard_Real tv2=v2; 449 450 if(MyHasBeenComputed) { 451 if( (MyParOnS1.X()==u1)&&(MyParOnS1.Y()==v1) 452 &&(MyParOnS2.X()==u2)&&(MyParOnS2.Y()==v2)) { 453 return(MyIsTangent); 454 } 455 else if(MyHasBeenComputedbis == Standard_False) { 456 MyTgbis = MyTg; 457 MyTguv1bis = MyTguv1; 458 MyTguv2bis = MyTguv2; 459 MyPntbis = MyPnt; 460 MyParOnS1bis = MyParOnS1; 461 MyParOnS2bis = MyParOnS2; 462 MyIsTangentbis = MyIsTangent; 463 MyHasBeenComputedbis = MyHasBeenComputed; 464 } 465 } 466 467 if(MyHasBeenComputedbis) { 468 if( (MyParOnS1bis.X()==u1)&&(MyParOnS1bis.Y()==v1) 469 &&(MyParOnS2bis.X()==u2)&&(MyParOnS2bis.Y()==v2)) { 470 471 gp_Vec TV(MyTg); 472 gp_Vec2d TV1(MyTguv1); 473 gp_Vec2d TV2(MyTguv2); 474 gp_Pnt TP(MyPnt); 475 gp_Pnt2d TP1(MyParOnS1); 476 gp_Pnt2d TP2(MyParOnS2); 477 Standard_Boolean TB=MyIsTangent; 478 479 MyTg = MyTgbis; 480 MyTguv1 = MyTguv1bis; 481 MyTguv2 = MyTguv2bis; 482 MyPnt = MyPntbis; 483 MyParOnS1 = MyParOnS1bis; 484 MyParOnS2 = MyParOnS2bis; 485 MyIsTangent = MyIsTangentbis; 486 487 MyTgbis = TV; 488 MyTguv1bis = TV1; 489 MyTguv2bis = TV2; 490 MyPntbis = TP; 491 MyParOnS1bis = TP1; 492 MyParOnS2bis = TP2; 493 MyIsTangentbis = TB; 494 495 return(MyIsTangent); 496 } 497 } 498 499 math_Vector X(1,2); 500 math_Vector BornInf(1,2), BornSup(1,2), Tolerance(1,2); 501 //--- ThePSurfaceTool::GetResolution(aPSurf,Tolerance(1),Tolerance(2)); 502 Tolerance(1) = 1.0e-8; Tolerance(2) = 1.0e-8; 503 Standard_Real binfu,bsupu,binfv,bsupv; 504 binfu = ThePSurfaceTool::FirstUParameter(aPSurf); 505 binfv = ThePSurfaceTool::FirstVParameter(aPSurf); 506 bsupu = ThePSurfaceTool::LastUParameter(aPSurf); 507 bsupv = ThePSurfaceTool::LastVParameter(aPSurf); 508 BornInf(1) = binfu; BornSup(1) = bsupu; 509 BornInf(2) = binfv; BornSup(2) = bsupv; 510 Standard_Real TranslationU = 0., TranslationV = 0.; 511 512 if (!FillInitialVectorOfSolution(u1, v1, u2, v2, 513 binfu, bsupu, binfv, bsupv, 514 X, 515 TranslationU, TranslationV)) 516 { 517 MyIsTangent=MyIsTangentbis=Standard_False; 518 MyHasBeenComputed = MyHasBeenComputedbis = Standard_False; 519 return(Standard_False); 520 } 521 522 Standard_Boolean aRsnldIsDone = Standard_False; 523 Standard_Real PourTesterU = X(1); 524 Standard_Real PourTesterV = X(2); 525 if (GetUseSolver()) 526 { 527 math_FunctionSetRoot Rsnld(MyZerImpFunc); 528 Rsnld.SetTolerance(Tolerance); 529 Rsnld.Perform(MyZerImpFunc, X, BornInf, BornSup); 530 aRsnldIsDone = Rsnld.IsDone(); 531 if (aRsnldIsDone) 532 Rsnld.Root(X); 533 } 534 if(aRsnldIsDone || !GetUseSolver()) 535 { 536 MyHasBeenComputed = Standard_True; 537 538 Standard_Real DistAvantApresU = Abs(PourTesterU-X(1)); 539 Standard_Real DistAvantApresV = Abs(PourTesterV-X(2)); 540 541 MyPnt = P = ThePSurfaceTool::Value(aPSurf, X(1), X(2)); 542 543 if( (DistAvantApresV <= 0.001 ) && 544 (DistAvantApresU <= 0.001 )) 545 { 546 gp_Vec aD1uPrm,aD1vPrm; 547 gp_Vec aD1uQuad,aD1vQuad; 548 549 if(MyImplicitFirst) 550 { 551 u2 = X(1)-TranslationU; 552 v2 = X(2)-TranslationV; 553 554 if(aQSurf.TypeQuadric() != GeomAbs_Plane) 555 { 556 while(u1-tu1>M_PI) u1-=M_PI+M_PI; 557 while(tu1-u1>M_PI) u1+=M_PI+M_PI; 558 } 559 560 MyParOnS1.SetCoord(tu1,tv1); 561 MyParOnS2.SetCoord(tu2,tv2); 562 563 gp_Pnt aP2; 564 565 ThePSurfaceTool::D1(aPSurf, X(1), X(2), P, aD1uPrm, aD1vPrm); 566 aQSurf.D1(u1,v1, aP2, aD1uQuad, aD1vQuad); 567 568 //Middle-point of P-P2 segment 569 P.BaryCenter(1.0, aP2, 1.0); 570 } 571 else 572 { 573 u1 = X(1)-TranslationU; 574 v1 = X(2)-TranslationV; 575 //aQSurf.Parameters(P, u2, v2); 576 if(aQSurf.TypeQuadric() != GeomAbs_Plane) 577 { 578 while(u2-tu2>M_PI) u2-=M_PI+M_PI; 579 while(tu2-u2>M_PI) u2+=M_PI+M_PI; 580 } 581 582 MyParOnS1.SetCoord(tu1,tv1); 583 MyParOnS2.SetCoord(tu2,tu2); 584 585 gp_Pnt aP2; 586 ThePSurfaceTool::D1(aPSurf, X(1), X(2), P, aD1uPrm, aD1vPrm); 587 588 aQSurf.D1(u2, v2, aP2, aD1uQuad, aD1vQuad); 589 590 //Middle-point of P-P2 segment 591 P.BaryCenter(1.0, aP2, 1.0); 592 } 593 594 MyPnt = P; 595 596 //Normals to the surfaces 597 gp_Vec aNormalPrm(aD1uPrm.Crossed(aD1vPrm)), 598 aNormalImp(aQSurf.Normale(MyPnt)); 599 600 const Standard_Real aSQMagnPrm = aNormalPrm.SquareMagnitude(), 601 aSQMagnImp = aNormalImp.SquareMagnitude(); 602 603 Standard_Boolean isPrmSingular = Standard_False, 604 isImpSingular = Standard_False; 605 606 if(IsSingular(aD1uPrm, aD1vPrm, aNullValue, anAngTol)) 607 { 608 isPrmSingular = Standard_True; 609 if(!SingularProcessing(aD1uPrm, aD1vPrm, Standard_True, 610 aNullValue, anAngTol, Tg, aPrmTg)) 611 { 612 MyIsTangent = Standard_False; 613 MyHasBeenComputed = MyHasBeenComputedbis = Standard_False; 614 return Standard_False; 615 } 616 617 MyTg = Tg; 618 } 619 else 620 { 621 aNormalPrm.Divide(sqrt(aSQMagnPrm)); 622 } 623 624 //Analogically for implicit surface 625 if(aSQMagnImp < aNullValue) 626 { 627 isImpSingular = Standard_True; 628 629 if(!SingularProcessing(aD1uQuad, aD1vQuad, !isPrmSingular, 630 aNullValue, anAngTol, Tg, aQuadTg)) 631 { 632 MyIsTangent = Standard_False; 633 MyHasBeenComputed = MyHasBeenComputedbis = Standard_False; 634 return Standard_False; 635 } 636 637 MyTg = Tg; 638 } 639 else 640 { 641 aNormalImp.Divide(sqrt(aSQMagnImp)); 642 } 643 644 if(isImpSingular && isPrmSingular) 645 { 646 //All is OK. All abnormal cases were processed above. 647 648 MyTguv1 = Tguv1; 649 MyTguv2 = Tguv2; 650 651 MyIsTangent=Standard_True; 652 return MyIsTangent; 653 } 654 else if(!(isImpSingular || isPrmSingular)) 655 { 656 //Processing pure non-singular case 657 //(3D- and 2D-tangents are still not defined) 658 659 //Ask to pay attention to the fact that here 660 //aNormalImp and aNormalPrm are normalized. 661 //Therefore, @ \left \| \vec{Tg} \right \| = 0.0 @ 662 //if and only if (aNormalImp || aNormalPrm). 663 Tg = aNormalImp.Crossed(aNormalPrm); 664 } 665 666 const Standard_Real aSQMagnTg = Tg.SquareMagnitude(); 667 668 if(aSQMagnTg < aNullValue) 669 { 670 MyIsTangent = Standard_False; 671 MyHasBeenComputed = MyHasBeenComputedbis = Standard_False; 672 return Standard_False; 673 } 674 675 //Normalize Tg vector 676 Tg.Divide(sqrt(aSQMagnTg)); 677 MyTg = Tg; 678 679 if(!isPrmSingular) 680 { 681 //If isPrmSingular==TRUE then aPrmTg has already been computed. 682 683 if(!NonSingularProcessing(aD1uPrm, aD1vPrm, Tg, aNullValue, anAngTol, aPrmTg)) 684 { 685 MyIsTangent = Standard_False; 686 MyHasBeenComputed = MyHasBeenComputedbis = Standard_False; 687 return Standard_False; 688 } 689 } 690 691 if(!isImpSingular) 692 { 693 //If isImpSingular==TRUE then aQuadTg has already been computed. 694 695 if(!NonSingularProcessing(aD1uQuad, aD1vQuad, Tg, aNullValue, anAngTol, aQuadTg)) 696 { 697 MyIsTangent = Standard_False; 698 MyHasBeenComputed = MyHasBeenComputedbis = Standard_False; 699 return Standard_False; 700 } 701 } 702 703 MyTguv1 = Tguv1; 704 MyTguv2 = Tguv2; 705 706 MyIsTangent=Standard_True; 707 708#ifdef OCCT_DEBUG 709 //cout << "+++++++++++++++++ ApproxInt_ImpPrmSvSurfaces::Compute(...) ++++++++++" << endl; 710 //printf( "P2d_1(%+10.20f, %+10.20f); P2d_2(%+10.20f, %+10.20f)\n" 711 // "P(%+10.20f, %+10.20f, %+10.20f);\n" 712 // "Tg = {%+10.20f, %+10.20f, %+10.20f};\n" 713 // "Tguv1 = {%+10.20f, %+10.20f};\n" 714 // "Tguv2 = {%+10.20f, %+10.20f}\n", 715 // u1, v1, u2, v2, 716 // P.X(), P.Y(), P.Z(), 717 // Tg.X(), Tg.Y(), Tg.Z(), 718 // Tguv1.X(), Tguv1.Y(), Tguv2.X(), Tguv2.Y()); 719 //cout << "-----------------------------------------------------------------------" << endl; 720#endif 721 722 return Standard_True; 723 } 724 else { 725 //-- cout<<" ApproxInt_ImpImpSvSurfaces.gxx : Distance apres recadrage Trop Grande "<<endl; 726 727 MyIsTangent=Standard_False; 728 MyHasBeenComputed = MyHasBeenComputedbis = Standard_False; 729 return Standard_False; 730 } 731 } 732 else { 733 MyIsTangent = Standard_False; 734 MyHasBeenComputed = MyHasBeenComputedbis = Standard_False; 735 return Standard_False; 736 } 737} 738 739//======================================================================= 740//function : SeekPoint 741//purpose : Computes point on curve and 742// parameters on the surfaces. 743//======================================================================= 744Standard_Boolean ApproxInt_ImpPrmSvSurfaces::SeekPoint(const Standard_Real u1, 745 const Standard_Real v1, 746 const Standard_Real u2, 747 const Standard_Real v2, 748 IntSurf_PntOn2S& Point) { 749 gp_Pnt aP; 750 gp_Vec aT; 751 gp_Vec2d aTS1,aTS2; 752 const IntSurf_Quadric& aQSurf = MyZerImpFunc.ISurface(); 753 const ThePSurface& aPSurf = MyZerImpFunc.PSurface(); 754 755 math_Vector X(1,2); 756 math_Vector BornInf(1,2), BornSup(1,2), Tolerance(1,2); 757 //--- ThePSurfaceTool::GetResolution(aPSurf,Tolerance(1),Tolerance(2)); 758 Tolerance(1) = 1.0e-8; Tolerance(2) = 1.0e-8; 759 Standard_Real binfu,bsupu,binfv,bsupv; 760 binfu = ThePSurfaceTool::FirstUParameter(aPSurf); 761 binfv = ThePSurfaceTool::FirstVParameter(aPSurf); 762 bsupu = ThePSurfaceTool::LastUParameter(aPSurf); 763 bsupv = ThePSurfaceTool::LastVParameter(aPSurf); 764 BornInf(1) = binfu; BornSup(1) = bsupu; 765 BornInf(2) = binfv; BornSup(2) = bsupv; 766 Standard_Real TranslationU = 0., TranslationV = 0.; 767 768 if (!FillInitialVectorOfSolution(u1, v1, u2, v2, 769 binfu, bsupu, binfv, bsupv, 770 X, 771 TranslationU, TranslationV)) 772 return Standard_False; 773 774 Standard_Real NewU1, NewV1, NewU2, NewV2; 775 776 math_FunctionSetRoot Rsnld(MyZerImpFunc); 777 Rsnld.SetTolerance(Tolerance); 778 Rsnld.Perform(MyZerImpFunc,X,BornInf,BornSup); 779 if(Rsnld.IsDone()) { 780 MyHasBeenComputed = Standard_True; 781 Rsnld.Root(X); 782 783 MyPnt = ThePSurfaceTool::Value(aPSurf, X(1), X(2)); 784 785 if(MyImplicitFirst) 786 { 787 NewU2 = X(1)-TranslationU; 788 NewV2 = X(2)-TranslationV; 789 790 aQSurf.Parameters(MyPnt, NewU1, NewV1); 791 //adjust U 792 if (aQSurf.TypeQuadric() != GeomAbs_Plane) 793 { 794 Standard_Real sign = (NewU1 > u1)? -1 : 1; 795 while (Abs(u1 - NewU1) > M_PI) 796 NewU1 += sign*(M_PI+M_PI); 797 } 798 } 799 else 800 { 801 NewU1 = X(1)-TranslationU; 802 NewV1 = X(2)-TranslationV; 803 804 aQSurf.Parameters(MyPnt, NewU2, NewV2); 805 //adjust U 806 if (aQSurf.TypeQuadric() != GeomAbs_Plane) 807 { 808 Standard_Real sign = (NewU2 > u2)? -1 : 1; 809 while (Abs(u2 - NewU2) > M_PI) 810 NewU2 += sign*(M_PI+M_PI); 811 } 812 } 813 } 814 else 815 return Standard_False; 816 817 Point.SetValue(MyPnt, NewU1, NewV1, NewU2, NewV2); 818 return Standard_True; 819} 820//-------------------------------------------------------------------------------- 821 822Standard_Boolean 823ApproxInt_ImpPrmSvSurfaces::FillInitialVectorOfSolution(const Standard_Real u1, 824 const Standard_Real v1, 825 const Standard_Real u2, 826 const Standard_Real v2, 827 const Standard_Real binfu, 828 const Standard_Real bsupu, 829 const Standard_Real binfv, 830 const Standard_Real bsupv, 831 math_Vector& X, 832 Standard_Real& TranslationU, 833 Standard_Real& TranslationV) 834{ 835 const ThePSurface& aPSurf = MyZerImpFunc.PSurface(); 836 837 math_Vector F(1,1); 838 839 TranslationU = 0.0; 840 TranslationV = 0.0; 841 842 if(MyImplicitFirst) { 843 if(u2<binfu-0.0000000001) { 844 if(ThePSurfaceTool::IsUPeriodic(aPSurf)) { 845 Standard_Real d = ThePSurfaceTool::UPeriod(aPSurf); 846 do { TranslationU+=d; } while(u2+TranslationU < binfu); 847 } 848 else 849 return(Standard_False); 850 } 851 else if(u2>bsupu+0.0000000001) { 852 if(ThePSurfaceTool::IsUPeriodic(aPSurf)) { 853 Standard_Real d = ThePSurfaceTool::UPeriod(aPSurf); 854 do { TranslationU-=d; } while(u2+TranslationU > bsupu); 855 } 856 else 857 return(Standard_False); 858 } 859 if(v2<binfv-0.0000000001) { 860 if(ThePSurfaceTool::IsVPeriodic(aPSurf)) { 861 Standard_Real d = ThePSurfaceTool::VPeriod(aPSurf); 862 do { TranslationV+=d; } while(v2+TranslationV < binfv); 863 } 864 else 865 return(Standard_False); 866 } 867 else if(v2>bsupv+0.0000000001) { 868 if(ThePSurfaceTool::IsVPeriodic(aPSurf)) { 869 Standard_Real d = ThePSurfaceTool::VPeriod(aPSurf); 870 do { TranslationV-=d; } while(v2+TranslationV > bsupv); 871 } 872 else 873 return(Standard_False); 874 } 875 X(1) = u2+TranslationU; 876 X(2) = v2+TranslationV; 877 } 878 else { 879 if(u1<binfu-0.0000000001) { 880 if(ThePSurfaceTool::IsUPeriodic(aPSurf)) { 881 Standard_Real d = ThePSurfaceTool::UPeriod(aPSurf); 882 do { TranslationU+=d; } while(u1+TranslationU < binfu); 883 } 884 else 885 return(Standard_False); 886 } 887 else if(u1>bsupu+0.0000000001) { 888 if(ThePSurfaceTool::IsUPeriodic(aPSurf)) { 889 Standard_Real d = ThePSurfaceTool::UPeriod(aPSurf); 890 do { TranslationU-=d; } while(u1+TranslationU > bsupu); 891 } 892 else 893 return(Standard_False); 894 } 895 if(v1<binfv-0.0000000001) { 896 if(ThePSurfaceTool::IsVPeriodic(aPSurf)) { 897 Standard_Real d = ThePSurfaceTool::VPeriod(aPSurf); 898 do { TranslationV+=d; } while(v1+TranslationV < binfv); 899 } 900 else 901 return(Standard_False); 902 } 903 else if(v1>bsupv+0.0000000001) { 904 if(ThePSurfaceTool::IsVPeriodic(aPSurf)) { 905 Standard_Real d = ThePSurfaceTool::VPeriod(aPSurf); 906 do { TranslationV-=d; } while(v1+TranslationV > bsupv); 907 } 908 else 909 return(Standard_False); 910 } 911 X(1) = u1+TranslationU; 912 X(2) = v1+TranslationV; 913 } 914 915 //---------------------------------------------------- 916 //Make a small step from boundaries in order to avoid 917 //finding "outboundaried" solution (Rsnld -> NotDone). 918 if (GetUseSolver()) 919 { 920 Standard_Real du = Max(Precision::Confusion(), ThePSurfaceTool::UResolution(aPSurf, Precision::Confusion())); 921 Standard_Real dv = Max(Precision::Confusion(), ThePSurfaceTool::VResolution(aPSurf, Precision::Confusion())); 922 if (X(1) - 0.0000000001 <= binfu) X(1) = X(1) + du; 923 if (X(1) + 0.0000000001 >= bsupu) X(1) = X(1) - du; 924 if (X(2) - 0.0000000001 <= binfv) X(2) = X(2) + dv; 925 if (X(2) + 0.0000000001 >= bsupv) X(2) = X(2) - dv; 926 } 927 928 return Standard_True; 929} 930