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