1 // Copyright (c) 1995-1999 Matra Datavision
2 // Copyright (c) 1999-2014 OPEN CASCADE SAS
4 // This file is part of Open CASCADE Technology software library.
6 // This library is free software; you can redistribute it and/or modify it under
7 // the terms of the GNU Lesser General Public License version 2.1 as published
8 // by the Free Software Foundation, with special exception defined in the file
9 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
10 // distribution for complete text of the license and disclaimer of any warranty.
12 // Alternatively, this file may be used under the terms of Open CASCADE
13 // commercial license or contractual agreement.
16 #include <Adaptor3d_HSurface.hxx>
17 #include <Adaptor3d_HSurfaceTool.hxx>
18 #include <Extrema_GenLocateExtPS.hxx>
19 #include <Geom_Surface.hxx>
22 #include <gp_Pnt2d.hxx>
23 #include <IntImp_ComputeTangence.hxx>
24 #include <IntSurf_LineOn2S.hxx>
25 #include <IntSurf_PntOn2S.hxx>
26 #include <IntWalk_PWalking.hxx>
27 #include <IntWalk_StatusDeflection.hxx>
28 #include <math_FunctionSetRoot.hxx>
29 #include <Precision.hxx>
30 #include <Standard_Failure.hxx>
31 #include <Standard_OutOfRange.hxx>
32 #include <StdFail_NotDone.hxx>
33 #include <TColgp_Array1OfPnt.hxx>
34 #include <TColStd_Array1OfReal.hxx>
36 //==================================================================================
37 // function : ComputePasInit
38 // purpose : estimate of max step : To avoid abrupt changes during change of isos
39 //==================================================================================
40 void IntWalk_PWalking::ComputePasInit(const Standard_Real theDeltaU1,
41 const Standard_Real theDeltaV1,
42 const Standard_Real theDeltaU2,
43 const Standard_Real theDeltaV2)
45 const Standard_Real aRangePart = 0.01;
46 const Standard_Real Increment = 2.0*pasMax;
47 const Handle(Adaptor3d_HSurface)&
48 Caro1 = myIntersectionOn2S.Function().AuxillarSurface1();
49 const Handle(Adaptor3d_HSurface)&
50 Caro2 = myIntersectionOn2S.Function().AuxillarSurface2();
52 const Standard_Real aDeltaU1=Abs(UM1-Um1);
53 const Standard_Real aDeltaV1=Abs(VM1-Vm1);
54 const Standard_Real aDeltaU2=Abs(UM2-Um2);
55 const Standard_Real aDeltaV2=Abs(VM2-Vm2);
57 //-- limit the reduction of uv box estimate to 0.01 natural box
58 //-- theDeltaU1 : On box of Inter
59 //-- aDeltaU1 : On parametric space
60 if(!Precision::IsInfinite(aDeltaU1))
61 pasuv[0]=Max(Increment*Max(theDeltaU1, aRangePart*aDeltaU1), pasuv[0]);
63 pasuv[0]=Max(Increment*theDeltaU1, pasuv[0]);
65 if(!Precision::IsInfinite(aDeltaV1))
66 pasuv[1]=Max(Increment*Max(theDeltaV1, aRangePart*aDeltaV1), pasuv[1]);
68 pasuv[1]=Max(Increment*theDeltaV1, pasuv[1]);
70 if(!Precision::IsInfinite(aDeltaU2))
71 pasuv[2]=Max(Increment*Max(theDeltaU2, aRangePart*aDeltaU2), pasuv[2]);
73 pasuv[2]=Max(Increment*theDeltaU2, pasuv[2]);
75 if(!Precision::IsInfinite(aDeltaV2))
76 pasuv[3]=Max(Increment*Max(theDeltaV2, aRangePart*aDeltaV2), pasuv[3]);
78 pasuv[3]=Max(Increment*theDeltaV2, pasuv[3]);
80 const Standard_Real ResoU1tol = Adaptor3d_HSurfaceTool::UResolution(Caro1, tolconf);
81 const Standard_Real ResoV1tol = Adaptor3d_HSurfaceTool::VResolution(Caro1, tolconf);
82 const Standard_Real ResoU2tol = Adaptor3d_HSurfaceTool::UResolution(Caro2, tolconf);
83 const Standard_Real ResoV2tol = Adaptor3d_HSurfaceTool::VResolution(Caro2, tolconf);
85 myStepMin[0] = Max(myStepMin[0], 2.0*ResoU1tol);
86 myStepMin[1] = Max(myStepMin[1], 2.0*ResoV1tol);
87 myStepMin[2] = Max(myStepMin[2], 2.0*ResoU2tol);
88 myStepMin[3] = Max(myStepMin[3], 2.0*ResoV2tol);
90 for(Standard_Integer i = 0; i < 4; i++)
92 pasuv[i]=Max(myStepMin[i], pasuv[i]);
96 //=======================================================================
97 //function : IsParallel
98 //purpose : Checks if theLine is parallel of some boundary of given
99 // surface (it is determined by theCheckSurf1 flag).
100 // Parallelism assumes small oscillations (swing is less or
101 // equal than theToler).
102 // Small lines (if first and last parameters in the Surface
103 // are almost equal) are classified as parallel (as same as
104 // any point can be considered as parallel of any line).
105 //=======================================================================
106 static void IsParallel(const Handle(IntSurf_LineOn2S)& theLine,
107 const Standard_Boolean theCheckSurf1,
108 const Standard_Real theToler,
109 Standard_Boolean& theIsUparallel,
110 Standard_Boolean& theIsVparallel)
112 const Standard_Integer aNbPointsMAX = 23;
114 theIsUparallel = theIsVparallel = Standard_True;
116 Standard_Integer aNbPoints = theLine->NbPoints();
117 if(aNbPoints > aNbPointsMAX)
119 aNbPoints = aNbPointsMAX;
121 else if(aNbPoints < 3)
123 //Here we cannot estimate parallelism.
124 //Do all same as for small lines
128 Standard_Real aStep = IntToReal(theLine->NbPoints()) / aNbPoints;
129 Standard_Real aNPoint = 1.0;
131 Standard_Real aUmin = RealLast(), aUmax = RealFirst(), aVmin = RealLast(), aVmax = RealFirst();
132 for(Standard_Integer aNum = 1; aNum <= aNbPoints; aNum++, aNPoint += aStep)
134 if(aNPoint > aNbPoints)
141 theLine->Value(RealToInt(aNPoint)).ParametersOnS1(u, v);
143 theLine->Value(RealToInt(aNPoint)).ParametersOnS2(u, v);
158 theIsVparallel = ((aUmax - aUmin) < theToler);
159 theIsUparallel = ((aVmax - aVmin) < theToler);
162 //=======================================================================
163 //function : Checking
164 //purpose : Check, if given point is in surface's boundaries.
165 // If "yes" then theFactTol = 0.0, else theFactTol is
166 // equal maximal deviation.
167 //=======================================================================
168 static Standard_Boolean Checking( const Handle(Adaptor3d_HSurface)& theASurf1,
169 const Handle(Adaptor3d_HSurface)& theASurf2,
170 Standard_Real& theU1,
171 Standard_Real& theV1,
172 Standard_Real& theU2,
173 Standard_Real& theV2,
174 Standard_Real& theFactTol)
176 const Standard_Real aTol = Precision::PConfusion();
177 const Standard_Real aU1bFirst = theASurf1->FirstUParameter();
178 const Standard_Real aU1bLast = theASurf1->LastUParameter();
179 const Standard_Real aU2bFirst = theASurf2->FirstUParameter();
180 const Standard_Real aU2bLast = theASurf2->LastUParameter();
181 const Standard_Real aV1bFirst = theASurf1->FirstVParameter();
182 const Standard_Real aV1bLast = theASurf1->LastVParameter();
183 const Standard_Real aV2bFirst = theASurf2->FirstVParameter();
184 const Standard_Real aV2bLast = theASurf2->LastVParameter();
186 Standard_Boolean isOnOrIn = Standard_True;
189 Standard_Real aDelta = aU1bFirst - theU1;
193 theFactTol = Max(theFactTol, aDelta);
194 isOnOrIn = Standard_False;
197 aDelta = theU1 - aU1bLast;
201 theFactTol = Max(theFactTol, aDelta);
202 isOnOrIn = Standard_False;
205 aDelta = aV1bFirst - theV1;
209 theFactTol = Max(theFactTol, aDelta);
210 isOnOrIn = Standard_False;
213 aDelta = theV1 - aV1bLast;
217 theFactTol = Max(theFactTol, aDelta);
218 isOnOrIn = Standard_False;
221 aDelta = aU2bFirst - theU2;
225 theFactTol = Max(theFactTol, aDelta);
226 isOnOrIn = Standard_False;
229 aDelta = theU2 - aU2bLast;
233 theFactTol = Max(theFactTol, aDelta);
234 isOnOrIn = Standard_False;
237 aDelta = aV2bFirst - theV2;
241 theFactTol = Max(theFactTol, aDelta);
242 isOnOrIn = Standard_False;
245 aDelta = theV2 - aV2bLast;
249 theFactTol = Max(theFactTol, aDelta);
250 isOnOrIn = Standard_False;
256 //==================================================================================
257 // function : IntWalk_PWalking::IntWalk_PWalking
259 //==================================================================================
260 IntWalk_PWalking::IntWalk_PWalking(const Handle(Adaptor3d_HSurface)& Caro1,
261 const Handle(Adaptor3d_HSurface)& Caro2,
262 const Standard_Real TolTangency,
263 const Standard_Real Epsilon,
264 const Standard_Real Deflection,
265 const Standard_Real Increment )
269 close(Standard_False),
272 myTolTang(TolTangency),
274 myIntersectionOn2S(Caro1,Caro2,TolTangency),
275 STATIC_BLOCAGE_SUR_PAS_TROP_GRAND(0),
276 STATIC_PRECEDENT_INFLEXION(0)
278 Standard_Real KELARG=20.;
280 pasMax=Increment*0.2; //-- June 25 99 after problems with precision
281 Um1 = Adaptor3d_HSurfaceTool::FirstUParameter(Caro1);
282 Vm1 = Adaptor3d_HSurfaceTool::FirstVParameter(Caro1);
283 UM1 = Adaptor3d_HSurfaceTool::LastUParameter(Caro1);
284 VM1 = Adaptor3d_HSurfaceTool::LastVParameter(Caro1);
286 Um2 = Adaptor3d_HSurfaceTool::FirstUParameter(Caro2);
287 Vm2 = Adaptor3d_HSurfaceTool::FirstVParameter(Caro2);
288 UM2 = Adaptor3d_HSurfaceTool::LastUParameter(Caro2);
289 VM2 = Adaptor3d_HSurfaceTool::LastVParameter(Caro2);
291 ResoU1 = Adaptor3d_HSurfaceTool::UResolution(Caro1,Precision::Confusion());
292 ResoV1 = Adaptor3d_HSurfaceTool::VResolution(Caro1,Precision::Confusion());
294 ResoU2 = Adaptor3d_HSurfaceTool::UResolution(Caro2,Precision::Confusion());
295 ResoV2 = Adaptor3d_HSurfaceTool::VResolution(Caro2,Precision::Confusion());
297 Standard_Real NEWRESO;
298 Standard_Real MAXVAL;
299 Standard_Real MAXVAL2;
301 MAXVAL = Abs(Um1); MAXVAL2 = Abs(UM1);
302 if(MAXVAL2 > MAXVAL) MAXVAL = MAXVAL2;
303 NEWRESO = ResoU1 * MAXVAL ;
304 if(NEWRESO > ResoU1 &&NEWRESO<10) { ResoU1 = NEWRESO; }
307 MAXVAL = Abs(Um2); MAXVAL2 = Abs(UM2);
308 if(MAXVAL2 > MAXVAL) MAXVAL = MAXVAL2;
309 NEWRESO = ResoU2 * MAXVAL ;
310 if(NEWRESO > ResoU2 && NEWRESO<10) { ResoU2 = NEWRESO; }
313 MAXVAL = Abs(Vm1); MAXVAL2 = Abs(VM1);
314 if(MAXVAL2 > MAXVAL) MAXVAL = MAXVAL2;
315 NEWRESO = ResoV1 * MAXVAL ;
316 if(NEWRESO > ResoV1 && NEWRESO<10) { ResoV1 = NEWRESO; }
319 MAXVAL = Abs(Vm2); MAXVAL2 = Abs(VM2);
320 if(MAXVAL2 > MAXVAL) MAXVAL = MAXVAL2;
321 NEWRESO = ResoV2 * MAXVAL ;
322 if(NEWRESO > ResoV2 && NEWRESO<10) { ResoV2 = NEWRESO; }
324 pasuv[0]=pasMax*Abs(UM1-Um1);
325 pasuv[1]=pasMax*Abs(VM1-Vm1);
326 pasuv[2]=pasMax*Abs(UM2-Um2);
327 pasuv[3]=pasMax*Abs(VM2-Vm2);
329 if(ResoU1>0.0001*pasuv[0]) ResoU1=0.00001*pasuv[0];
330 if(ResoV1>0.0001*pasuv[1]) ResoV1=0.00001*pasuv[1];
331 if(ResoU2>0.0001*pasuv[2]) ResoU2=0.00001*pasuv[2];
332 if(ResoV2>0.0001*pasuv[3]) ResoV2=0.00001*pasuv[3];
335 if(Adaptor3d_HSurfaceTool::IsUPeriodic(Caro1)==Standard_False) {
336 //UM1+=KELARG*pasuv[0]; Um1-=KELARG*pasuv[0];
339 Standard_Real t = UM1-Um1;
340 if(t<Adaptor3d_HSurfaceTool::UPeriod(Caro1)) {
341 t=0.5*(Adaptor3d_HSurfaceTool::UPeriod(Caro1)-t);
342 t=(t>KELARG*pasuv[0])? KELARG*pasuv[0] : t;
347 if(Adaptor3d_HSurfaceTool::IsVPeriodic(Caro1)==Standard_False) {
348 //VM1+=KELARG*pasuv[1]; Vm1-=KELARG*pasuv[1];
351 Standard_Real t = VM1-Vm1;
352 if(t<Adaptor3d_HSurfaceTool::VPeriod(Caro1)) {
353 t=0.5*(Adaptor3d_HSurfaceTool::VPeriod(Caro1)-t);
354 t=(t>KELARG*pasuv[1])? KELARG*pasuv[1] : t;
359 if(Adaptor3d_HSurfaceTool::IsUPeriodic(Caro2)==Standard_False) {
360 //UM2+=KELARG*pasuv[2]; Um2-=KELARG*pasuv[2];
363 Standard_Real t = UM2-Um2;
364 if(t<Adaptor3d_HSurfaceTool::UPeriod(Caro2)) {
365 t=0.5*(Adaptor3d_HSurfaceTool::UPeriod(Caro2)-t);
366 t=(t>KELARG*pasuv[2])? KELARG*pasuv[2] : t;
371 if(Adaptor3d_HSurfaceTool::IsVPeriodic(Caro2)==Standard_False) {
372 //VM2+=KELARG*pasuv[3]; Vm2-=KELARG*pasuv[3];
375 Standard_Real t = VM2-Vm2;
376 if(t<Adaptor3d_HSurfaceTool::VPeriod(Caro2)) {
377 t=0.5*(Adaptor3d_HSurfaceTool::VPeriod(Caro2)-t);
378 t=(t>KELARG*pasuv[3])? KELARG*pasuv[3] : t;
383 myStepMin[0] = 100.0*ResoU1;
384 myStepMin[1] = 100.0*ResoV1;
385 myStepMin[2] = 100.0*ResoU2;
386 myStepMin[3] = 100.0*ResoV2;
388 //-- ComputePasInit(pasuv,Um1,UM1,Vm1,VM1,Um2,UM2,Vm2,VM2,Caro1,Caro2);
390 for (Standard_Integer i = 0; i<=3;i++) {
393 pasInit[i] = pasSav[i] = pasuv[i];
398 //==================================================================================
399 // function : IntWalk_PWalking
401 //==================================================================================
402 IntWalk_PWalking::IntWalk_PWalking(const Handle(Adaptor3d_HSurface)& Caro1,
403 const Handle(Adaptor3d_HSurface)& Caro2,
404 const Standard_Real TolTangency,
405 const Standard_Real Epsilon,
406 const Standard_Real Deflection,
407 const Standard_Real Increment,
408 const Standard_Real U1,
409 const Standard_Real V1,
410 const Standard_Real U2,
411 const Standard_Real V2)
415 close(Standard_False),
418 myTolTang(TolTangency),
420 myIntersectionOn2S(Caro1,Caro2,TolTangency),
421 STATIC_BLOCAGE_SUR_PAS_TROP_GRAND(0),
422 STATIC_PRECEDENT_INFLEXION(0)
424 Standard_Real KELARG=20.;
426 pasMax=Increment*0.2; //-- June 25 99 after problems with precision
428 Um1 = Adaptor3d_HSurfaceTool::FirstUParameter(Caro1);
429 Vm1 = Adaptor3d_HSurfaceTool::FirstVParameter(Caro1);
430 UM1 = Adaptor3d_HSurfaceTool::LastUParameter(Caro1);
431 VM1 = Adaptor3d_HSurfaceTool::LastVParameter(Caro1);
433 Um2 = Adaptor3d_HSurfaceTool::FirstUParameter(Caro2);
434 Vm2 = Adaptor3d_HSurfaceTool::FirstVParameter(Caro2);
435 UM2 = Adaptor3d_HSurfaceTool::LastUParameter(Caro2);
436 VM2 = Adaptor3d_HSurfaceTool::LastVParameter(Caro2);
438 ResoU1 = Adaptor3d_HSurfaceTool::UResolution(Caro1,Precision::Confusion());
439 ResoV1 = Adaptor3d_HSurfaceTool::VResolution(Caro1,Precision::Confusion());
441 ResoU2 = Adaptor3d_HSurfaceTool::UResolution(Caro2,Precision::Confusion());
442 ResoV2 = Adaptor3d_HSurfaceTool::VResolution(Caro2,Precision::Confusion());
444 Standard_Real NEWRESO, MAXVAL, MAXVAL2;
448 if(MAXVAL2 > MAXVAL) {
451 NEWRESO = ResoU1 * MAXVAL ;
452 if(NEWRESO > ResoU1) {
458 if(MAXVAL2 > MAXVAL){
461 NEWRESO = ResoU2 * MAXVAL ;
462 if(NEWRESO > ResoU2) {
468 if(MAXVAL2 > MAXVAL) {
471 NEWRESO = ResoV1 * MAXVAL ;
472 if(NEWRESO > ResoV1) {
478 if(MAXVAL2 > MAXVAL){
481 NEWRESO = ResoV2 * MAXVAL ;
482 if(NEWRESO > ResoV2) {
486 pasuv[0]=pasMax*Abs(UM1-Um1);
487 pasuv[1]=pasMax*Abs(VM1-Vm1);
488 pasuv[2]=pasMax*Abs(UM2-Um2);
489 pasuv[3]=pasMax*Abs(VM2-Vm2);
491 if(Adaptor3d_HSurfaceTool::IsUPeriodic(Caro1)==Standard_False) {
492 UM1+=KELARG*pasuv[0];
493 Um1-=KELARG*pasuv[0];
496 Standard_Real t = UM1-Um1;
497 if(t<Adaptor3d_HSurfaceTool::UPeriod(Caro1)) {
498 t=0.5*(Adaptor3d_HSurfaceTool::UPeriod(Caro1)-t);
499 t=(t>KELARG*pasuv[0])? KELARG*pasuv[0] : t;
505 if(Adaptor3d_HSurfaceTool::IsVPeriodic(Caro1)==Standard_False) {
506 VM1+=KELARG*pasuv[1];
507 Vm1-=KELARG*pasuv[1];
510 Standard_Real t = VM1-Vm1;
511 if(t<Adaptor3d_HSurfaceTool::VPeriod(Caro1)) {
512 t=0.5*(Adaptor3d_HSurfaceTool::VPeriod(Caro1)-t);
513 t=(t>KELARG*pasuv[1])? KELARG*pasuv[1] : t;
518 if(Adaptor3d_HSurfaceTool::IsUPeriodic(Caro2)==Standard_False) {
519 UM2+=KELARG*pasuv[2];
520 Um2-=KELARG*pasuv[2];
523 Standard_Real t = UM2-Um2;
524 if(t<Adaptor3d_HSurfaceTool::UPeriod(Caro2)) {
525 t=0.5*(Adaptor3d_HSurfaceTool::UPeriod(Caro2)-t);
526 t=(t>KELARG*pasuv[2])? KELARG*pasuv[2] : t;
532 if(Adaptor3d_HSurfaceTool::IsVPeriodic(Caro2)==Standard_False) {
533 VM2+=KELARG*pasuv[3];
534 Vm2-=KELARG*pasuv[3];
537 Standard_Real t = VM2-Vm2;
538 if(t<Adaptor3d_HSurfaceTool::VPeriod(Caro2)) {
539 t=0.5*(Adaptor3d_HSurfaceTool::VPeriod(Caro2)-t);
540 t=(t>KELARG*pasuv[3])? KELARG*pasuv[3] : t;
545 //-- ComputePasInit(pasuv,Um1,UM1,Vm1,VM1,Um2,UM2,Vm2,VM2,Caro1,Caro2);
547 for (Standard_Integer i = 0; i<=3;i++) {
548 pasInit[i] = pasSav[i] = pasuv[i];
551 if(ResoU1>0.0001*pasuv[0]) ResoU1=0.00001*pasuv[0];
552 if(ResoV1>0.0001*pasuv[1]) ResoV1=0.00001*pasuv[1];
553 if(ResoU2>0.0001*pasuv[2]) ResoU2=0.00001*pasuv[2];
554 if(ResoV2>0.0001*pasuv[3]) ResoV2=0.00001*pasuv[3];
556 myStepMin[0] = 100.0*ResoU1;
557 myStepMin[1] = 100.0*ResoV1;
558 myStepMin[2] = 100.0*ResoU2;
559 myStepMin[3] = 100.0*ResoV2;
562 TColStd_Array1OfReal Par(1,4);
570 //==================================================================================
571 // function : PerformFirstPoint
573 //==================================================================================
574 Standard_Boolean IntWalk_PWalking::PerformFirstPoint (const TColStd_Array1OfReal& ParDep,
575 IntSurf_PntOn2S& FirstPoint)
578 close = Standard_False;
581 TColStd_Array1OfReal Param(1,4);
583 for (i=1; i<=4; ++i) {
584 Param(i) = ParDep(i);
586 //-- calculate the first solution point
587 math_FunctionSetRoot Rsnld(myIntersectionOn2S.Function());
589 myIntersectionOn2S.Perform(Param,Rsnld);
590 if (!myIntersectionOn2S.IsDone()) {
591 return Standard_False;
594 if (myIntersectionOn2S.IsEmpty()) {
595 return Standard_False;
598 FirstPoint = myIntersectionOn2S.Point();
599 return Standard_True;
601 //==================================================================================
602 // function : Perform
604 //==================================================================================
605 void IntWalk_PWalking::Perform(const TColStd_Array1OfReal& ParDep)
607 Perform(ParDep,Um1,Vm1,Um2,Vm2,UM1,VM1,UM2,VM2);
610 //=======================================================================
611 //function : SQDistPointSurface
612 //purpose : Returns square distance between thePnt and theSurf.
613 // (theU0, theV0) is initial point for extrema
614 //=======================================================================
615 static Standard_Real SQDistPointSurface(const gp_Pnt &thePnt,
616 const Adaptor3d_Surface& theSurf,
617 const Standard_Real theU0,
618 const Standard_Real theV0)
620 const Extrema_GenLocateExtPS aExtPS(thePnt, theSurf, theU0, theV0,
621 Precision::PConfusion(), Precision::PConfusion());
625 return aExtPS.SquareDistance();
628 //==================================================================================
629 // function : IsTangentExtCheck
630 // purpose : Additional check if the surfaces are tangent.
631 // Checks if any point in one surface lie in another surface
632 // (with given tolerance)
633 //==================================================================================
634 static Standard_Boolean IsTangentExtCheck(const Handle(Adaptor3d_HSurface)& theSurf1,
635 const Handle(Adaptor3d_HSurface)& theSurf2,
636 const Standard_Real theU10,
637 const Standard_Real theV10,
638 const Standard_Real theU20,
639 const Standard_Real theV20,
640 const Standard_Real theToler,
641 const Standard_Real theArrStep[])
645 gp_Vec aDu1, aDv1, aDu2, aDv2;
646 theSurf1->D1(theU10, theV10, aPt, aDu1, aDv1);
647 theSurf2->D1(theU20, theV20, aPt, aDu2, aDv2);
649 const gp_Vec aN1(aDu1.Crossed(aDv1)),
650 aN2(aDu2.Crossed(aDv2));
651 const Standard_Real aDP = aN1.Dot(aN2),
652 aSQ1 = aN1.SquareMagnitude(),
653 aSQ2 = aN2.SquareMagnitude();
655 if((aSQ1 < RealSmall()) || (aSQ2 < RealSmall()))
656 return Standard_True; //Tangent
658 if(aDP*aDP < 0.9998*aSQ1*aSQ2)
659 {//cos(ang N1<->N2) < 0.9999
660 return Standard_False; //Not tangent
664 //For two faces (2^2 = 4)
665 const Standard_Real aSQToler = 4.0*theToler*theToler;
666 const Standard_Integer aNbItems = 4;
667 const Standard_Real aParUS1[aNbItems] = { theU10 + theArrStep[0],
668 theU10 - theArrStep[0],
670 const Standard_Real aParVS1[aNbItems] = { theV10, theV10,
671 theV10 + theArrStep[1],
672 theV10 - theArrStep[1]};
673 const Standard_Real aParUS2[aNbItems] = { theU20 + theArrStep[2],
674 theU20 - theArrStep[2],
676 const Standard_Real aParVS2[aNbItems] = { theV20, theV20,
677 theV20 + theArrStep[3],
678 theV20 - theArrStep[3]};
680 for(Standard_Integer i = 0; i < aNbItems; i++)
682 gp_Pnt aP(theSurf1->Value(aParUS1[i], aParVS1[i]));
683 const Standard_Real aSqDist = SQDistPointSurface(aP, theSurf2->Surface(), theU20, theV20);
684 if(aSqDist > aSQToler)
685 return Standard_False;
688 for(Standard_Integer i = 0; i < aNbItems; i++)
690 gp_Pnt aP(theSurf2->Value(aParUS2[i], aParVS2[i]));
691 const Standard_Real aSqDist = SQDistPointSurface(aP, theSurf1->Surface(), theU10, theV10);
692 if(aSqDist > aSQToler)
693 return Standard_False;
696 return Standard_True;
699 //==================================================================================
700 // function : Perform
702 //==================================================================================
703 void IntWalk_PWalking::Perform(const TColStd_Array1OfReal& ParDep,
704 const Standard_Real u1min,
705 const Standard_Real v1min,
706 const Standard_Real u2min,
707 const Standard_Real v2min,
708 const Standard_Real u1max,
709 const Standard_Real v1max,
710 const Standard_Real u2max,
711 const Standard_Real v2max)
713 const Standard_Real aSQDistMax = 1.0e-14;
716 Standard_Integer NbPasOKConseq=0;
717 TColStd_Array1OfReal Param(1,4);
718 IntImp_ConstIsoparametric ChoixIso;
721 done = Standard_False;
724 const Handle(Adaptor3d_HSurface)& Caro1 =myIntersectionOn2S.Function().AuxillarSurface1();
725 const Handle(Adaptor3d_HSurface)& Caro2 =myIntersectionOn2S.Function().AuxillarSurface2();
727 const Standard_Real UFirst1 = Adaptor3d_HSurfaceTool::FirstUParameter(Caro1);
728 const Standard_Real VFirst1 = Adaptor3d_HSurfaceTool::FirstVParameter(Caro1);
729 const Standard_Real ULast1 = Adaptor3d_HSurfaceTool::LastUParameter (Caro1);
730 const Standard_Real VLast1 = Adaptor3d_HSurfaceTool::LastVParameter (Caro1);
732 const Standard_Real UFirst2 = Adaptor3d_HSurfaceTool::FirstUParameter(Caro2);
733 const Standard_Real VFirst2 = Adaptor3d_HSurfaceTool::FirstVParameter(Caro2);
734 const Standard_Real ULast2 = Adaptor3d_HSurfaceTool::LastUParameter (Caro2);
735 const Standard_Real VLast2 = Adaptor3d_HSurfaceTool::LastVParameter (Caro2);
737 ComputePasInit(u1max - u1min,v1max - v1min,u2max - u2min,v2max - v2min);
739 for (Standard_Integer i=0; i<4; ++i)
746 pasInit[i] = pasSav[i] = pasuv[i];
749 line = new IntSurf_LineOn2S ();
751 for (Standard_Integer i=1; i<=4; ++i)
755 //-- reproduce steps uv connected to surfaces Caro1 and Caro2
756 //-- pasuv[] and pasSav[] are modified during the marching
757 for(Standard_Integer i = 0; i < 4; ++i)
759 pasSav[i] = pasuv[i] = pasInit[i];
762 //-- calculate the first solution point
763 math_FunctionSetRoot Rsnld(myIntersectionOn2S.Function());
765 ChoixIso = myIntersectionOn2S.Perform(Param,Rsnld);
766 if (!myIntersectionOn2S.IsDone())
772 if (myIntersectionOn2S.IsEmpty())
777 if(myIntersectionOn2S.IsTangent())
782 Standard_Boolean Arrive, DejaReparti;
783 const Standard_Integer RejectIndexMAX = 250000;
784 Standard_Integer IncKey, RejectIndex;
787 DejaReparti = Standard_False;
791 previousPoint = myIntersectionOn2S.Point();
792 previoustg = Standard_False;
793 previousd = myIntersectionOn2S.Direction();
794 previousd1 = myIntersectionOn2S.DirectionOnS1();
795 previousd2 = myIntersectionOn2S.DirectionOnS2();
798 firstd1 = previousd1;
799 firstd2 = previousd2;
800 tgfirst = tglast = Standard_False;
801 choixIsoSav = ChoixIso;
802 //------------------------------------------------------------
803 //-- Test if the first point of marching corresponds
804 //-- to a point on borders.
805 //-- In this case, DejaReparti is initialized as True
807 pf = previousPoint.Value();
808 Standard_Boolean bTestFirstPoint = Standard_True;
810 previousPoint.Parameters(Param(1),Param(2),Param(3),Param(4));
812 if(IsTangentExtCheck(Caro1, Caro2, Param(1), Param(2), Param(3), Param(4), myTolTang, pasuv))
815 AddAPoint(line,previousPoint);
817 IntWalk_StatusDeflection Status = IntWalk_OK;
818 Standard_Boolean NoTestDeflection = Standard_False;
819 Standard_Real SvParam[4], f;
820 Standard_Integer LevelOfEmptyInmyIntersectionOn2S=0;
821 Standard_Integer LevelOfPointConfondu = 0;
822 Standard_Integer LevelOfIterWithoutAppend = -1;
825 const Standard_Real aTol[4] = { Epsilon(u1max - u1min),
826 Epsilon(v1max - v1min),
827 Epsilon(u2max - u2min),
828 Epsilon(v2max - v2min)};
829 Arrive = Standard_False;
832 LevelOfIterWithoutAppend++;
833 if(LevelOfIterWithoutAppend>20)
835 Arrive = Standard_True;
839 RepartirOuDiviser(DejaReparti,ChoixIso,Arrive);
840 LevelOfIterWithoutAppend = 0;
846 case IntImp_UIsoparametricOnCaro1: f = Abs(previousd1.X()); break;
847 case IntImp_VIsoparametricOnCaro1: f = Abs(previousd1.Y()); break;
848 case IntImp_UIsoparametricOnCaro2: f = Abs(previousd2.X()); break;
849 case IntImp_VIsoparametricOnCaro2: f = Abs(previousd2.Y()); break;
857 previousPoint.Parameters(Param(1),Param(2),Param(3),Param(4));
860 Standard_Real aIncKey, aEps, dP1, dP2, dP3, dP4;
862 dP1 = sensCheminement * pasuv[0] * previousd1.X() /f;
863 dP2 = sensCheminement * pasuv[1] * previousd1.Y() /f;
864 dP3 = sensCheminement * pasuv[2] * previousd2.X() /f;
865 dP4 = sensCheminement * pasuv[3] * previousd2.Y() /f;
867 aIncKey=5.*(Standard_Real)IncKey;
869 if(ChoixIso == IntImp_UIsoparametricOnCaro1 && Abs(dP1) < aEps)
874 if(ChoixIso == IntImp_VIsoparametricOnCaro1 && Abs(dP2) < aEps)
879 if(ChoixIso == IntImp_UIsoparametricOnCaro2 && Abs(dP3) < aEps)
884 if(ChoixIso == IntImp_VIsoparametricOnCaro2 && Abs(dP4) < aEps)
894 //==========================
900 Standard_Integer aTryNumber = 0;
901 Standard_Real isBadPoint = Standard_False;
902 IntImp_ConstIsoparametric aBestIso = ChoixIso;
905 isBadPoint = Standard_False;
907 ChoixIso= myIntersectionOn2S.Perform(Param, Rsnld, aBestIso);
909 if (myIntersectionOn2S.IsDone() && !myIntersectionOn2S.IsEmpty())
911 //If we go along any surface boundary then it is possible
912 //to find "outboundaried" point.
913 //Nevertheless, if this deflection is quite small, we will be
914 //able to adjust this point to the boundary.
916 Standard_Real aNewPnt[4], anAbsParamDist[4];
917 myIntersectionOn2S.Point().Parameters(aNewPnt[0], aNewPnt[1], aNewPnt[2], aNewPnt[3]);
918 const Standard_Real aParMin[4] = {u1min, v1min, u2min, v2min};
919 const Standard_Real aParMax[4] = {u1max, v1max, u2max, v2max};
921 for(Standard_Integer i = 0; i < 4; i++)
923 if(Abs(aNewPnt[i] - aParMin[i]) < aTol[i])
924 aNewPnt[i] = aParMin[i];
925 else if(Abs(aNewPnt[i] - aParMax[i]) < aTol[i])
926 aNewPnt[i] = aParMax[i];
929 if (aNewPnt[0] < u1min || aNewPnt[0] > u1max ||
930 aNewPnt[1] < v1min || aNewPnt[1] > v1max ||
931 aNewPnt[2] < u2min || aNewPnt[2] > u2max ||
932 aNewPnt[3] < v2min || aNewPnt[3] > v2max)
934 break; // Out of borders, handle this later.
937 myIntersectionOn2S.ChangePoint().SetValue(aNewPnt[0],
942 anAbsParamDist[0] = Abs(Param(1) - dP1 - aNewPnt[0]);
943 anAbsParamDist[1] = Abs(Param(2) - dP2 - aNewPnt[1]);
944 anAbsParamDist[2] = Abs(Param(3) - dP3 - aNewPnt[2]);
945 anAbsParamDist[3] = Abs(Param(4) - dP4 - aNewPnt[3]);
946 if (anAbsParamDist[0] < ResoU1 &&
947 anAbsParamDist[1] < ResoV1 &&
948 anAbsParamDist[2] < ResoU2 &&
949 anAbsParamDist[3] < ResoV2 &&
950 Status != IntWalk_PasTropGrand)
952 isBadPoint = Standard_True;
953 aBestIso = IntImp_ConstIsoparametric((aBestIso + 1) % 4);
956 } while (isBadPoint && ++aTryNumber <= 4);
958 if (!myIntersectionOn2S.IsDone())
960 //end of line, division
961 Arrive = Standard_False;
966 RepartirOuDiviser(DejaReparti, ChoixIso, Arrive);
970 //== Calculation of exact point from Param(.) is possible
971 if (myIntersectionOn2S.IsEmpty())
973 Standard_Real u1,v1,u2,v2;
974 previousPoint.Parameters(u1,v1,u2,v2);
976 Arrive = Standard_False;
977 if(u1<UFirst1 || u1>ULast1)
979 Arrive=Standard_True;
982 if(u2<UFirst2 || u2>ULast2)
984 Arrive=Standard_True;
987 if(v1<VFirst1 || v1>VLast1)
989 Arrive=Standard_True;
992 if(v2<VFirst2 || v2>VLast2)
994 Arrive=Standard_True;
997 RepartirOuDiviser(DejaReparti,ChoixIso,Arrive);
998 LevelOfEmptyInmyIntersectionOn2S++;
1000 if(LevelOfEmptyInmyIntersectionOn2S>10)
1010 //============================================================
1011 //== A point has been found : T E S T D E F L E C T I O N
1012 //============================================================
1013 if(NoTestDeflection)
1015 NoTestDeflection = Standard_False;
1019 if(--LevelOfEmptyInmyIntersectionOn2S<=0)
1021 LevelOfEmptyInmyIntersectionOn2S=0;
1022 if(LevelOfIterWithoutAppend < 10)
1024 Status = TestDeflection(ChoixIso);
1036 //============================================================
1037 //== T r a i t e m e n t s u r S t a t u s ==
1038 //============================================================
1039 if(LevelOfPointConfondu > 5)
1041 Status = IntWalk_ArretSurPoint;
1042 LevelOfPointConfondu = 0;
1045 if(Status==IntWalk_OK)
1048 if(NbPasOKConseq >= 5)
1051 Standard_Boolean pastroppetit;
1056 pastroppetit=Standard_True;
1058 if(pasuv[0]<pasInit[0])
1060 t = (pasInit[0]-pasuv[0])*0.25;
1061 if(t>0.1*pasInit[0])
1067 pastroppetit=Standard_False;
1070 if(pasuv[1]<pasInit[1])
1072 t = (pasInit[1]-pasuv[1])*0.25;
1073 if(t>0.1*pasInit[1]) {
1078 pastroppetit=Standard_False;
1081 if(pasuv[2]<pasInit[2])
1083 t = (pasInit[2]-pasuv[2])*0.25;
1084 if(t>0.1*pasInit[2])
1090 pastroppetit=Standard_False;
1093 if(pasuv[3]<pasInit[3])
1095 t = (pasInit[3]-pasuv[3])*0.25;
1096 if(t>0.1*pasInit[3]) {
1100 pastroppetit=Standard_False;
1114 pastroppetit=Standard_False;
1118 while(pastroppetit);
1120 }//Status==IntWalk_OK
1127 case IntWalk_ArretSurPointPrecedent:
1129 Arrive = Standard_False;
1130 RepartirOuDiviser(DejaReparti, ChoixIso, Arrive);
1133 case IntWalk_PasTropGrand:
1135 Param(1)=SvParam[0];
1136 Param(2)=SvParam[1];
1137 Param(3)=SvParam[2];
1138 Param(4)=SvParam[3];
1140 if(LevelOfIterWithoutAppend > 5)
1142 for (Standard_Integer i = 0; i < 4; i++)
1144 if (pasSav[i] > pasInit[i])
1147 const Standard_Real aDelta = (pasInit[i]-pasSav[i])*0.25;
1149 if(aDelta > Epsilon(pasInit[i]))
1151 pasInit[i] -= aDelta;
1152 LevelOfIterWithoutAppend=0;
1159 case IntWalk_PointConfondu:
1161 LevelOfPointConfondu++;
1163 if(LevelOfPointConfondu>5)
1165 Standard_Boolean pastroppetit;
1169 pastroppetit=Standard_True;
1171 for(Standard_Integer i = 0; i < 4; i++)
1173 if(pasuv[i]<pasInit[i])
1175 pasuv[i]+=(pasInit[i]-pasuv[i])*0.25;
1176 pastroppetit=Standard_False;
1192 pastroppetit=Standard_False;
1196 while(pastroppetit);
1201 case IntWalk_StepTooSmall:
1203 Standard_Boolean hasStepBeenIncreased = Standard_False;
1205 for(Standard_Integer i = 0; i < 4; i++)
1207 const Standard_Real aNewStep = Min(1.5*pasuv[i], pasInit[i]);
1208 if(aNewStep > pasuv[i])
1210 pasuv[i] = aNewStep;
1211 hasStepBeenIncreased = Standard_True;
1215 if(hasStepBeenIncreased)
1217 Param(1)=SvParam[0];
1218 Param(2)=SvParam[1];
1219 Param(3)=SvParam[2];
1220 Param(4)=SvParam[3];
1222 LevelOfIterWithoutAppend = 0;
1228 case IntWalk_ArretSurPoint://006
1230 //=======================================================
1231 //== Stop Test t : Frame on Param(.) ==
1232 //=======================================================
1234 Arrive = TestArret(DejaReparti,Param,ChoixIso);
1235 // JMB 30th December 1999.
1236 // Some statement below should not be put in comment because they are useful.
1237 // See grid CTO 909 A1 which infinitely loops
1238 if(Arrive==Standard_False && Status==IntWalk_ArretSurPoint)
1240 Arrive=Standard_True;
1242 cout << "IntWalk_PWalking_1.gxx: Problems with intersection"<<endl;
1248 NbPasOKConseq = -10;
1253 //=====================================================
1254 //== Param(.) is in the limits ==
1255 //== and does not end a closed line ==
1256 //=====================================================
1257 //== Check on the current point of myInters
1258 Standard_Boolean pointisvalid = Standard_False;
1260 Standard_Real u1,v1,u2,v2;
1261 myIntersectionOn2S.Point().Parameters(u1,v1,u2,v2);
1264 if(u1 <= UM1 && u2 <= UM2 && v1 <= VM1 &&
1265 v2 <= VM2 && u1 >= Um1 && u2 >= Um2 &&
1266 v1 >= Vm1 && v2 >= Vm2)
1268 pointisvalid=Standard_True;
1275 previousPoint = myIntersectionOn2S.Point();
1276 previoustg = myIntersectionOn2S.IsTangent();
1280 previousd = myIntersectionOn2S.Direction();
1281 previousd1 = myIntersectionOn2S.DirectionOnS1();
1282 previousd2 = myIntersectionOn2S.DirectionOnS2();
1284 //=====================================================
1285 //== Check on the previous Point
1287 Standard_Real u1,v1,u2,v2;
1288 previousPoint.Parameters(u1,v1,u2,v2);
1289 if( u1 <= UM1 && u2 <= UM2 && v1 <= VM1 &&
1290 v2 <= VM2 && u1 >= Um1 && u2 >= Um2 &&
1291 v1 >= Vm1 && v2 >= Vm2)
1293 pl = previousPoint.Value();
1296 if(pf.SquareDistance(pl) < aSQDistMax)
1306 bTestFirstPoint = Standard_False;
1310 AddAPoint(line,previousPoint);
1313 if(RejectIndex >= RejectIndexMAX)
1315 Arrive = Standard_True;
1320 LevelOfIterWithoutAppend = 0;
1324 //====================================================
1326 if(Status == IntWalk_ArretSurPoint)
1328 RepartirOuDiviser(DejaReparti,ChoixIso,Arrive);
1332 if (line->NbPoints() == 2)
1334 pasSav[0] = pasuv[0];
1335 pasSav[1] = pasuv[1];
1336 pasSav[2] = pasuv[2];
1337 pasSav[3] = pasuv[3];
1345 //================= la ligne est fermee ===============
1346 AddAPoint(line,line->Value(1)); //ligne fermee
1347 LevelOfIterWithoutAppend=0;
1351 //====================================================
1352 //== Param was not in the limits (was reframed)
1353 //====================================================
1354 Standard_Boolean bPrevNotTangent = !previoustg || !myIntersectionOn2S.IsTangent();
1356 IntImp_ConstIsoparametric SauvChoixIso = ChoixIso;
1357 ChoixIso = myIntersectionOn2S.Perform(Param,Rsnld,ChoixIso);
1359 if(!myIntersectionOn2S.IsEmpty()) //002
1361 // mutially outpasses in the square or intersection in corner
1363 if(TestArret(Standard_True,Param,ChoixIso))
1365 NbPasOKConseq = -10;
1366 ChoixIso = myIntersectionOn2S.Perform(Param,Rsnld,ChoixIso);
1368 if(!myIntersectionOn2S.IsEmpty())
1370 previousPoint = myIntersectionOn2S.Point();
1371 previoustg = myIntersectionOn2S.IsTangent();
1375 previousd = myIntersectionOn2S.Direction();
1376 previousd1 = myIntersectionOn2S.DirectionOnS1();
1377 previousd2 = myIntersectionOn2S.DirectionOnS2();
1380 pl = previousPoint.Value();
1384 if(pf.SquareDistance(pl) < aSQDistMax)
1394 bTestFirstPoint = Standard_False;
1398 AddAPoint(line,previousPoint);
1401 if(RejectIndex >= RejectIndexMAX)
1403 Arrive = Standard_True;
1408 LevelOfIterWithoutAppend=0;
1409 RepartirOuDiviser(DejaReparti,ChoixIso,Arrive);
1413 //fail framing divides the step
1414 Arrive = Standard_False;
1415 RepartirOuDiviser(DejaReparti,ChoixIso,Arrive);
1416 NoTestDeflection = Standard_True;
1417 ChoixIso = SauvChoixIso;
1422 // save the last point
1423 // to revert to it if the current point is out of bounds
1425 IntSurf_PntOn2S previousPointSave = previousPoint;
1426 Standard_Boolean previoustgSave = previoustg;
1427 gp_Dir previousdSave = previousd;
1428 gp_Dir2d previousd1Save = previousd1;
1429 gp_Dir2d previousd2Save = previousd2;
1431 previousPoint = myIntersectionOn2S.Point();
1432 previoustg = myIntersectionOn2S.IsTangent();
1433 Arrive = Standard_False;
1437 previousd = myIntersectionOn2S.Direction();
1438 previousd1 = myIntersectionOn2S.DirectionOnS1();
1439 previousd2 = myIntersectionOn2S.DirectionOnS2();
1442 //========================================
1443 //== Check on PreviousPoint @@
1446 Standard_Real u1,v1,u2,v2;
1447 previousPoint.Parameters(u1,v1,u2,v2);
1449 //To save initial 2d points
1450 gp_Pnt2d ParamPntOnS1(Param(1), Param(2));
1451 gp_Pnt2d ParamPntOnS2(Param(3), Param(4));
1453 ///////////////////////////
1461 Standard_Boolean bFlag1, bFlag2;
1462 Standard_Real aTol2D=1.e-11;
1464 bFlag1=u1 >= Um1-aTol2D && v1 >= Vm1-aTol2D && u1 <= UM1+aTol2D && v1 <= VM1+aTol2D;
1465 bFlag2=u2 >= Um2-aTol2D && v2 >= Vm2-aTol2D && u2 <= UM2+aTol2D && v2 <= VM2+aTol2D;
1466 if (bFlag1 && bFlag2)
1468 if (line->NbPoints() > 1)
1470 IntSurf_PntOn2S prevprevPoint = line->Value(line->NbPoints()-1);
1471 Standard_Real ppU1, ppV1, ppU2, ppV2;
1472 prevprevPoint.Parameters(ppU1, ppV1, ppU2, ppV2);
1473 Standard_Real pU1, pV1, pU2, pV2;
1474 previousPointSave.Parameters(pU1, pV1, pU2, pV2);
1475 gp_Vec2d V1onS1(gp_Pnt2d(ppU1, ppV1), gp_Pnt2d(pU1, pV1));
1476 gp_Vec2d V2onS1(gp_Pnt2d(pU1, pV1), gp_Pnt2d(u1, v1));
1477 gp_Vec2d V1onS2(gp_Pnt2d(ppU2, ppV2), gp_Pnt2d(pU2, pV2));
1478 gp_Vec2d V2onS2(gp_Pnt2d(pU2, pV2), gp_Pnt2d(u2, v2));
1480 const Standard_Real aDot1 = V1onS1 * V2onS1;
1481 const Standard_Real aDot2 = V1onS2 * V2onS2;
1483 if ((aDot1 < 0.0) || (aDot2 < 0.0))
1485 Arrive = Standard_True;
1490 if(u1 <= UM1 && u2 <= UM2 && v1 <= VM1 &&
1491 v2 <= VM2 && u1 >= Um1 && u2 >= Um2 &&
1492 v1 >= Vm1 && v2 >= Vm2) {
1495 pl = previousPoint.Value();
1499 if(pf.SquareDistance(pl) < aSQDistMax)
1510 bTestFirstPoint = Standard_False;
1514 //To avoid walking around the same point
1515 //in the tangent zone near a border
1519 //There are three consecutive points:
1520 //previousPointSave -> ParamPnt -> curPnt.
1522 Standard_Real prevU1, prevV1, prevU2, prevV2;
1523 previousPointSave.Parameters(prevU1, prevV1, prevU2, prevV2);
1524 gp_Pnt2d prevPntOnS1(prevU1, prevV1), prevPntOnS2(prevU2, prevV2);
1525 gp_Pnt2d curPntOnS1(u1, v1), curPntOnS2(u2, v2);
1526 gp_Vec2d PrevToParamOnS1(prevPntOnS1, ParamPntOnS1);
1527 gp_Vec2d PrevToCurOnS1(prevPntOnS1, curPntOnS1);
1528 gp_Vec2d PrevToParamOnS2(prevPntOnS2, ParamPntOnS2);
1529 gp_Vec2d PrevToCurOnS2(prevPntOnS2, curPntOnS2);
1530 Standard_Real MaxAngle = 3*M_PI/4;
1531 Standard_Real anAngleS1 = 0.0, anAngleS2 = 0.0;
1532 const Standard_Real aSQMParS1 = PrevToParamOnS1.SquareMagnitude();
1533 const Standard_Real aSQMParS2 = PrevToParamOnS2.SquareMagnitude();
1534 const Standard_Real aSQMCurS1 = PrevToCurOnS1.SquareMagnitude();
1535 const Standard_Real aSQMCurS2 = PrevToCurOnS2.SquareMagnitude();
1537 if(aSQMCurS1 < gp::Resolution())
1539 //We came back to the one of previos point.
1540 //Therefore, we must break;
1544 else if(aSQMParS1 < gp::Resolution())
1546 //We are walking along tangent zone.
1547 //It should be continued.
1552 anAngleS1 = Abs(PrevToParamOnS1.Angle(PrevToCurOnS1));
1555 if(aSQMCurS2 < gp::Resolution())
1557 //We came back to the one of previos point.
1558 //Therefore, we must break;
1562 else if(aSQMParS2 < gp::Resolution())
1564 //We are walking along tangent zone.
1565 //It should be continued;
1570 anAngleS2 = Abs(PrevToParamOnS2.Angle(PrevToCurOnS2));
1573 if ((anAngleS1 > MaxAngle) && (anAngleS2 > MaxAngle))
1575 Arrive = Standard_True;
1580 //Check singularity.
1581 //I.e. check if we are walking along direction, which does not
1582 //result in comming to any point (i.e. derivative
1583 //3D-intersection curve along this direction is equal to 0).
1584 //A sphere with direction {dU=1, dV=0} from point
1585 //(U=0, V=M_PI/2) can be considered as example for
1586 //this case (we cannot find another 3D-point if we go thus).
1588 //Direction chosen along 1st and 2nd surface correspondingly
1589 const gp_Vec2d aDirS1(prevPntOnS1, curPntOnS1),
1590 aDirS2(prevPntOnS2, curPntOnS2);
1593 gp_Vec aDuS1, aDvS1, aDuS2, aDvS2;
1595 myIntersectionOn2S.Function().AuxillarSurface1()->
1596 D1(curPntOnS1.X(), curPntOnS1.Y(), aPtemp, aDuS1, aDvS1);
1597 myIntersectionOn2S.Function().AuxillarSurface2()->
1598 D1(curPntOnS2.X(), curPntOnS2.Y(), aPtemp, aDuS2, aDvS2);
1600 //Derivative WLine along (it is vector-function indeed)
1602 //(https://en.wikipedia.org/wiki/Directional_derivative#Variation_using_only_direction_of_vector).
1603 //F1 - on the 1st surface, F2 - on the 2nd surface.
1604 //x, y, z - coordinates of derivative vector.
1605 const Standard_Real aF1x = aDuS1.X()*aDirS1.X() +
1606 aDvS1.X()*aDirS1.Y();
1607 const Standard_Real aF1y = aDuS1.Y()*aDirS1.X() +
1608 aDvS1.Y()*aDirS1.Y();
1609 const Standard_Real aF1z = aDuS1.Z()*aDirS1.X() +
1610 aDvS1.Z()*aDirS1.Y();
1611 const Standard_Real aF2x = aDuS2.X()*aDirS2.X() +
1612 aDvS2.X()*aDirS2.Y();
1613 const Standard_Real aF2y = aDuS2.Y()*aDirS2.X() +
1614 aDvS2.Y()*aDirS2.Y();
1615 const Standard_Real aF2z = aDuS2.Z()*aDirS2.X() +
1616 aDvS2.Z()*aDirS2.Y();
1618 const Standard_Real aF1 = aF1x*aF1x + aF1y*aF1y + aF1z*aF1z;
1619 const Standard_Real aF2 = aF2x*aF2x + aF2y*aF2y + aF2z*aF2z;
1621 if((aF1 < gp::Resolution()) && (aF2 < gp::Resolution()))
1623 //All derivative are equal to 0. Therefore, there is
1624 //no point in going along direction chosen.
1625 Arrive = Standard_True;
1629 }//if (previoustg) cond.
1631 ////////////////////////////////////////
1632 AddAPoint(line,previousPoint);
1635 if(RejectIndex >= RejectIndexMAX)
1637 Arrive = Standard_True;
1643 LevelOfIterWithoutAppend=0;
1644 Arrive = Standard_True;
1648 // revert to the last correctly calculated point
1649 previousPoint = previousPointSave;
1650 previoustg = previoustgSave;
1651 previousd = previousdSave;
1652 previousd1 = previousd1Save;
1653 previousd2 = previousd2Save;
1658 Standard_Boolean wasExtended = Standard_False;
1660 if(Arrive && myIntersectionOn2S.IsTangent() && bPrevNotTangent)
1662 if(ExtendLineInCommonZone(SauvChoixIso, DejaReparti))
1664 wasExtended = Standard_True;
1665 Arrive = Standard_False;
1666 ChoixIso = SauvChoixIso;
1670 RepartirOuDiviser(DejaReparti,ChoixIso,Arrive);
1673 myIntersectionOn2S.IsDone() && !myIntersectionOn2S.IsEmpty() &&
1674 myIntersectionOn2S.IsTangent() && bPrevNotTangent &&
1677 if(ExtendLineInCommonZone(SauvChoixIso, DejaReparti))
1679 wasExtended = Standard_True;
1680 Arrive = Standard_False;
1681 ChoixIso = SauvChoixIso;
1684 }//else !TestArret() $
1685 }//$$ end successful framing on border (!myIntersectionOn2S.IsEmpty())
1688 //echec framing on border; division of step
1689 Arrive = Standard_False;
1690 NoTestDeflection = Standard_True;
1691 RepartirOuDiviser(DejaReparti,ChoixIso,Arrive);
1693 }//$$$ end framing on border (!close)
1694 }//004 fin TestArret return Arrive = True
1695 } // 006case IntWalk_ArretSurPoint: end Processing Status = OK or ArretSurPoint
1696 } //007 switch(Status)
1697 } //008 end processing point (TEST DEFLECTION)
1698 } //009 end processing line (else if myIntersectionOn2S.IsDone())
1699 } //010 end if first departure point allows marching while (!Arrive)
1701 done = Standard_True;
1703 // ===========================================================================================================
1704 // function: ExtendLineInCommonZone
1705 // purpose: Extends already computed line inside tangent zone in the direction given by theChoixIso.
1706 // Returns Standard_True if the line was extended through tangent zone and the last computed point
1707 // is outside the tangent zone (but it is not put into the line). Otherwise returns Standard_False.
1708 // ===========================================================================================================
1709 Standard_Boolean IntWalk_PWalking::ExtendLineInCommonZone(const IntImp_ConstIsoparametric theChoixIso,
1710 const Standard_Boolean theDirectionFlag)
1712 Standard_Boolean bOutOfTangentZone = Standard_False;
1713 Standard_Boolean bStop = !myIntersectionOn2S.IsTangent();
1714 Standard_Integer dIncKey = 1;
1715 TColStd_Array1OfReal Param(1,4);
1716 IntWalk_StatusDeflection Status = IntWalk_OK;
1717 Standard_Integer nbIterWithoutAppend = 0;
1718 Standard_Integer nbEqualPoints = 0;
1719 Standard_Integer parit = 0;
1720 Standard_Integer uvit = 0;
1721 IntSurf_SequenceOfPntOn2S aSeqOfNewPoint;
1724 nbIterWithoutAppend++;
1726 if((nbIterWithoutAppend > 20) || (nbEqualPoints > 20)) {
1728 cout<<"Infinite loop detected. Stop iterations (IntWalk_PWalking_1.gxx)" << endl;
1730 bStop = Standard_True;
1733 Standard_Real f = 0.;
1735 switch (theChoixIso)
1737 case IntImp_UIsoparametricOnCaro1: f = Abs(previousd1.X()); break;
1738 case IntImp_VIsoparametricOnCaro1: f = Abs(previousd1.Y()); break;
1739 case IntImp_UIsoparametricOnCaro2: f = Abs(previousd2.X()); break;
1740 case IntImp_VIsoparametricOnCaro2: f = Abs(previousd2.Y()); break;
1745 previousPoint.Parameters(Param(1),Param(2),Param(3),Param(4));
1747 Standard_Real dP1 = sensCheminement * pasuv[0] * previousd1.X() /f;
1748 Standard_Real dP2 = sensCheminement * pasuv[1] * previousd1.Y() /f;
1749 Standard_Real dP3 = sensCheminement * pasuv[2] * previousd2.X() /f;
1750 Standard_Real dP4 = sensCheminement * pasuv[3] * previousd2.Y() /f;
1752 if(theChoixIso == IntImp_UIsoparametricOnCaro1 && Abs(dP1) < 1.e-7) dP1 *= (5. * (Standard_Real)dIncKey);
1753 if(theChoixIso == IntImp_VIsoparametricOnCaro1 && Abs(dP2) < 1.e-7) dP2 *= (5. * (Standard_Real)dIncKey);
1754 if(theChoixIso == IntImp_UIsoparametricOnCaro2 && Abs(dP3) < 1.e-7) dP3 *= (5. * (Standard_Real)dIncKey);
1755 if(theChoixIso == IntImp_VIsoparametricOnCaro2 && Abs(dP4) < 1.e-7) dP4 *= (5. * (Standard_Real)dIncKey);
1761 Standard_Real SvParam[4];
1762 IntImp_ConstIsoparametric ChoixIso = theChoixIso;
1764 for(parit = 0; parit < 4; parit++) {
1765 SvParam[parit] = Param(parit+1);
1767 math_FunctionSetRoot Rsnld(myIntersectionOn2S.Function());
1768 ChoixIso = myIntersectionOn2S.Perform(Param,Rsnld, theChoixIso);
1770 if (!myIntersectionOn2S.IsDone()) {
1771 return bOutOfTangentZone;
1774 if (myIntersectionOn2S.IsEmpty()) {
1775 return bOutOfTangentZone;
1778 Status = TestDeflection(ChoixIso);
1780 if(Status == IntWalk_OK) {
1782 for(uvit = 0; uvit < 4; uvit++) {
1783 if(pasuv[uvit] < pasInit[uvit]) {
1784 pasuv[uvit] = pasInit[uvit];
1790 case IntWalk_ArretSurPointPrecedent:
1792 bStop = Standard_True;
1793 bOutOfTangentZone = !myIntersectionOn2S.IsTangent();
1796 case IntWalk_PasTropGrand:
1798 for(parit = 0; parit < 4; parit++) {
1799 Param(parit+1) = SvParam[parit];
1801 Standard_Boolean bDecrease = Standard_False;
1803 for(uvit = 0; uvit < 4; uvit++) {
1804 if(pasSav[uvit] < pasInit[uvit]) {
1805 pasInit[uvit] -= (pasInit[uvit] - pasSav[uvit]) * 0.1;
1806 bDecrease = Standard_True;
1810 if(bDecrease) nbIterWithoutAppend--;
1813 case IntWalk_PointConfondu:
1815 for(uvit = 0; uvit < 4; uvit++) {
1816 if(pasuv[uvit] < pasInit[uvit]) {
1817 pasuv[uvit] += (pasInit[uvit] - pasuv[uvit]) * 0.1;
1823 case IntWalk_ArretSurPoint:
1826 bStop = TestArret(theDirectionFlag, Param, ChoixIso);
1831 Standard_Real u11,v11,u12,v12;
1832 myIntersectionOn2S.Point().Parameters(u11,v11,u12,v12);
1833 Standard_Real u21,v21,u22,v22;
1834 previousPoint.Parameters(u21,v21,u22,v22);
1836 if(((fabs(u11-u21) < ResoU1) && (fabs(v11-v21) < ResoV1)) ||
1837 ((fabs(u12-u22) < ResoU2) && (fabs(v12-v22) < ResoV2))) {
1846 bStop = bStop || !myIntersectionOn2S.IsTangent();
1847 bOutOfTangentZone = !myIntersectionOn2S.IsTangent();
1850 Standard_Boolean pointisvalid = Standard_False;
1851 Standard_Real u1,v1,u2,v2;
1852 myIntersectionOn2S.Point().Parameters(u1,v1,u2,v2);
1854 if(u1 <= UM1 && u2 <= UM2 && v1 <= VM1 &&
1855 v2 <= VM2 && u1 >= Um1 && u2 >= Um2 &&
1856 v1 >= Vm1 && v2 >= Vm2)
1857 pointisvalid = Standard_True;
1860 previousPoint = myIntersectionOn2S.Point();
1861 previoustg = myIntersectionOn2S.IsTangent();
1864 previousd = myIntersectionOn2S.Direction();
1865 previousd1 = myIntersectionOn2S.DirectionOnS1();
1866 previousd2 = myIntersectionOn2S.DirectionOnS2();
1868 Standard_Boolean bAddPoint = Standard_True;
1870 if(line->NbPoints() >= 1) {
1871 gp_Pnt pf = line->Value(1).Value();
1872 gp_Pnt pl = previousPoint.Value();
1874 if(pf.Distance(pl) < Precision::Confusion()) {
1876 if(dIncKey == 5000) return bOutOfTangentZone;
1877 else bAddPoint = Standard_False;
1882 aSeqOfNewPoint.Append(previousPoint);
1883 nbIterWithoutAppend = 0;
1887 if (line->NbPoints() == 2) {
1888 for(uvit = 0; uvit < 4; uvit++) {
1889 pasSav[uvit] = pasuv[uvit];
1893 if ( !pointisvalid ) {
1894 // decrease step if out of bounds
1895 // otherwise the same calculations will be
1896 // repeated several times
1897 if ( ( u1 > UM1 ) || ( u1 < Um1 ) )
1900 if ( ( v1 > VM1 ) || ( v1 < Vm1 ) )
1903 if ( ( u2 > UM2 ) || ( u2 < Um2 ) )
1906 if ( ( v2 > VM2 ) || ( v2 < Vm2 ) )
1911 if(close && (line->NbPoints() >= 1)) {
1913 if(!bOutOfTangentZone) {
1914 aSeqOfNewPoint.Append(line->Value(1)); // line end
1916 nbIterWithoutAppend = 0;
1919 ChoixIso = myIntersectionOn2S.Perform(Param, Rsnld, theChoixIso);
1921 if(myIntersectionOn2S.IsEmpty()) {
1922 bStop = !myIntersectionOn2S.IsTangent();
1923 bOutOfTangentZone = !myIntersectionOn2S.IsTangent();
1926 Standard_Boolean bAddPoint = Standard_True;
1927 Standard_Boolean pointisvalid = Standard_False;
1929 previousPoint = myIntersectionOn2S.Point();
1930 Standard_Real u1,v1,u2,v2;
1931 previousPoint.Parameters(u1,v1,u2,v2);
1933 if(u1 <= UM1 && u2 <= UM2 && v1 <= VM1 &&
1934 v2 <= VM2 && u1 >= Um1 && u2 >= Um2 &&
1935 v1 >= Vm1 && v2 >= Vm2)
1936 pointisvalid = Standard_True;
1940 if(line->NbPoints() >= 1) {
1941 gp_Pnt pf = line->Value(1).Value();
1942 gp_Pnt pl = previousPoint.Value();
1944 if(pf.Distance(pl) < Precision::Confusion()) {
1946 if(dIncKey == 5000) return bOutOfTangentZone;
1947 else bAddPoint = Standard_False;
1951 if(bAddPoint && !bOutOfTangentZone) {
1952 aSeqOfNewPoint.Append(previousPoint);
1953 nbIterWithoutAppend = 0;
1968 Standard_Boolean bExtendLine = Standard_False;
1969 Standard_Real u1 = 0., v1 = 0., u2 = 0., v2 = 0.;
1971 Standard_Integer pit = 0;
1973 for(pit = 0; !bExtendLine && (pit < 2); pit++) {
1975 previousPoint.Parameters(u1,v1,u2,v2);
1977 if(aSeqOfNewPoint.Length() > 0)
1978 aSeqOfNewPoint.Value(aSeqOfNewPoint.Length()).Parameters(u1,v1,u2,v2);
1983 if(((u1 - Um1) < ResoU1) ||
1984 ((UM1 - u1) < ResoU1) ||
1985 ((u2 - Um2) < ResoU2) ||
1986 ((UM2 - u2) < ResoU2) ||
1987 ((v1 - Vm1) < ResoV1) ||
1988 ((VM1 - v1) < ResoV1) ||
1989 ((v2 - Vm2) < ResoV2) ||
1990 ((VM2 - v2) < ResoV2))
1991 bExtendLine = Standard_True;
1995 // if(Status == IntWalk_OK || Status == IntWalk_ArretSurPoint) {
1996 if(Status == IntWalk_OK) {
1997 bExtendLine = Standard_True;
1999 if(aSeqOfNewPoint.Length() > 1) {
2000 TColStd_Array1OfReal FirstParams(0, 3), LastParams(0, 3), Resolutions(0, 3);
2001 Resolutions(0) = ResoU1; Resolutions(1) = ResoV1; Resolutions(2) = ResoU2; Resolutions(3) = ResoV2;
2003 aSeqOfNewPoint(1).Parameters(FirstParams.ChangeValue(0), FirstParams.ChangeValue(1),
2004 FirstParams.ChangeValue(2), FirstParams.ChangeValue(3));
2005 aSeqOfNewPoint(aSeqOfNewPoint.Length()).Parameters(LastParams.ChangeValue(0),
2006 LastParams.ChangeValue(1),
2007 LastParams.ChangeValue(2),
2008 LastParams.ChangeValue(3));
2009 Standard_Integer indexofiso = 0;
2011 if(theChoixIso == IntImp_UIsoparametricOnCaro1) indexofiso = 0;
2012 if(theChoixIso == IntImp_VIsoparametricOnCaro1) indexofiso = 1;
2013 if(theChoixIso == IntImp_UIsoparametricOnCaro2) indexofiso = 2;
2014 if(theChoixIso == IntImp_VIsoparametricOnCaro2) indexofiso = 3;
2016 Standard_Integer afirstindex = (indexofiso < 2) ? 0 : 2;
2017 gp_Vec2d aTangentZoneDir(gp_Pnt2d(FirstParams.Value(afirstindex), FirstParams.Value(afirstindex + 1)),
2018 gp_Pnt2d(LastParams.Value(afirstindex), LastParams.Value(afirstindex + 1)));
2020 gp_Dir2d anIsoDir(0, 1);
2022 if((indexofiso == 1) || (indexofiso == 3))
2023 anIsoDir = gp_Dir2d(1, 0);
2025 if(aTangentZoneDir.SquareMagnitude() > gp::Resolution()) {
2026 Standard_Real piquota = M_PI*0.25;
2028 if(fabs(aTangentZoneDir.Angle(anIsoDir)) > piquota) {
2029 Standard_Integer ii = 1, nextii = 2;
2031 Standard_Real asqresol = gp::Resolution();
2032 asqresol *= asqresol;
2035 aSeqOfNewPoint(ii).Parameters(FirstParams.ChangeValue(0), FirstParams.ChangeValue(1),
2036 FirstParams.ChangeValue(2), FirstParams.ChangeValue(3));
2037 aSeqOfNewPoint(ii + 1).Parameters(LastParams.ChangeValue(0), LastParams.ChangeValue(1),
2038 LastParams.ChangeValue(2), LastParams.ChangeValue(3));
2039 d1 = gp_Vec2d(gp_Pnt2d(FirstParams.Value(afirstindex),
2040 FirstParams.Value(afirstindex + 1)),
2041 gp_Pnt2d(LastParams.Value(afirstindex),
2042 LastParams.Value(afirstindex + 1)));
2045 while((d1.SquareMagnitude() < asqresol) &&
2046 (ii < aSeqOfNewPoint.Length()));
2050 while(nextii < aSeqOfNewPoint.Length()) {
2052 gp_Vec2d nextd1(0, 0);
2053 Standard_Integer jj = nextii;
2056 aSeqOfNewPoint(jj).Parameters(FirstParams.ChangeValue(0), FirstParams.ChangeValue(1),
2057 FirstParams.ChangeValue(2), FirstParams.ChangeValue(3));
2058 aSeqOfNewPoint(jj + 1).Parameters(LastParams.ChangeValue(0), LastParams.ChangeValue(1),
2059 LastParams.ChangeValue(2), LastParams.ChangeValue(3));
2060 nextd1 = gp_Vec2d(gp_Pnt2d(FirstParams.Value(afirstindex),
2061 FirstParams.Value(afirstindex + 1)),
2062 gp_Pnt2d(LastParams.Value(afirstindex),
2063 LastParams.Value(afirstindex + 1)));
2067 while((nextd1.SquareMagnitude() < asqresol) &&
2068 (jj < aSeqOfNewPoint.Length()));
2071 if(fabs(d1.Angle(nextd1)) > piquota) {
2072 bExtendLine = Standard_False;
2078 // end if(fabs(aTangentZoneDir.Angle(anIsoDir)
2085 return Standard_False;
2087 Standard_Integer i = 0;
2089 for(i = 1; i <= aSeqOfNewPoint.Length(); i++) {
2090 AddAPoint(line, aSeqOfNewPoint.Value(i));
2093 return bOutOfTangentZone;
2096 //=======================================================================
2097 //function : DistanceMinimizeByGradient
2099 //=======================================================================
2100 Standard_Boolean IntWalk_PWalking::
2101 DistanceMinimizeByGradient( const Handle(Adaptor3d_HSurface)& theASurf1,
2102 const Handle(Adaptor3d_HSurface)& theASurf2,
2103 Standard_Real& theU1,
2104 Standard_Real& theV1,
2105 Standard_Real& theU2,
2106 Standard_Real& theV2,
2107 const Standard_Real theStep0U1V1,
2108 const Standard_Real theStep0U2V2)
2110 const Standard_Integer aNbIterMAX = 60;
2111 const Standard_Real aTol = 1.0e-14;
2112 Handle(Geom_Surface) aS1, aS2;
2114 if (theASurf1->GetType() != GeomAbs_BezierSurface &&
2115 theASurf1->GetType() != GeomAbs_BSplineSurface)
2116 return Standard_True;
2117 if (theASurf2->GetType() != GeomAbs_BezierSurface &&
2118 theASurf2->GetType() != GeomAbs_BSplineSurface)
2119 return Standard_True;
2121 Standard_Boolean aStatus = Standard_False;
2124 gp_Vec aD1u, aD1v, aD2U, aD2V;
2126 theASurf1->D1(theU1, theV1, aP1, aD1u, aD1v);
2127 theASurf2->D1(theU2, theV2, aP2, aD2U, aD2V);
2129 Standard_Real aSQDistPrev = aP1.SquareDistance(aP2);
2131 gp_Vec aP12(aP1, aP2);
2133 Standard_Real aGradFu(-aP12.Dot(aD1u));
2134 Standard_Real aGradFv(-aP12.Dot(aD1v));
2135 Standard_Real aGradFU( aP12.Dot(aD2U));
2136 Standard_Real aGradFV( aP12.Dot(aD2V));
2138 Standard_Real aSTEPuv = theStep0U1V1, aStepUV = theStep0U2V2;
2140 Standard_Boolean flRepeat = Standard_True;
2141 Standard_Integer aNbIter = aNbIterMAX;
2145 Standard_Real anAdd = aGradFu*aSTEPuv;
2146 Standard_Real aPARu = (anAdd >= 0.0)?
2147 (theU1 - Max(anAdd, Epsilon(theU1))) :
2148 (theU1 + Max(-anAdd, Epsilon(theU1)));
2149 anAdd = aGradFv*aSTEPuv;
2150 Standard_Real aPARv = (anAdd >= 0.0)?
2151 (theV1 - Max(anAdd, Epsilon(theV1))) :
2152 (theV1 + Max(-anAdd, Epsilon(theV1)));
2153 anAdd = aGradFU*aStepUV;
2154 Standard_Real aParU = (anAdd >= 0.0)?
2155 (theU2 - Max(anAdd, Epsilon(theU2))) :
2156 (theU2 + Max(-anAdd, Epsilon(theU2)));
2157 anAdd = aGradFV*aStepUV;
2158 Standard_Real aParV = (anAdd >= 0.0)?
2159 (theV2 - Max(anAdd, Epsilon(theV2))) :
2160 (theV2 + Max(-anAdd, Epsilon(theV2)));
2164 theASurf1->D1(aPARu, aPARv, aPt1, aD1u, aD1v);
2165 theASurf2->D1(aParU, aParV, aPt2, aD2U, aD2V);
2167 Standard_Real aSQDist = aPt1.SquareDistance(aPt2);
2169 if(aSQDist < aSQDistPrev)
2171 aSQDistPrev = aSQDist;
2177 aStatus = aSQDistPrev < aTol;
2185 flRepeat = Standard_False;
2189 theASurf1->D1(theU1, theV1, aPt1, aD1u, aD1v);
2190 theASurf2->D1(theU2, theV2, aPt2, aD2U, aD2V);
2192 gp_Vec aPt12(aPt1, aPt2);
2193 aGradFu = -aPt12.Dot(aD1u);
2194 aGradFv = -aPt12.Dot(aD1v);
2195 aGradFU = aPt12.Dot(aD2U);
2196 aGradFV = aPt12.Dot(aD2V);
2197 aSTEPuv = theStep0U1V1;
2198 aStepUV = theStep0U2V2;
2206 //=======================================================================
2207 //function : DistanceMinimizeByExtrema
2209 //=======================================================================
2210 Standard_Boolean IntWalk_PWalking::
2211 DistanceMinimizeByExtrema(const Handle(Adaptor3d_HSurface)& theASurf,
2212 const gp_Pnt& theP0,
2213 Standard_Real& theU0,
2214 Standard_Real& theV0,
2215 const Standard_Real theStep0U,
2216 const Standard_Real theStep0V)
2218 const Standard_Real aTol = 1.0e-14;
2220 gp_Vec aD1Su, aD1Sv, aD2Su, aD2Sv, aD2SuvTemp;
2221 Standard_Real aSQDistPrev = RealLast();
2222 Standard_Real aU = theU0, aV = theV0;
2224 Standard_Integer aNbIter = 10;
2227 theASurf->D2(aU, aV, aPS, aD1Su, aD1Sv, aD2Su, aD2Sv, aD2SuvTemp);
2229 gp_Vec aVec(theP0, aPS);
2231 Standard_Real aSQDist = aVec.SquareMagnitude();
2233 if(aSQDist >= aSQDistPrev)
2236 aSQDistPrev = aSQDist;
2241 if(aSQDistPrev < aTol)
2245 const Standard_Real aF1 = aD1Su.Dot(aVec), aF2 = aD1Sv.Dot(aVec);
2248 const Standard_Real aDf1u = aD2Su.Dot(aVec) + aD1Su.Dot(aD1Su),
2249 aDf1v = aD2Su.Dot(aD1Sv),
2251 aDf2v = aD2Sv.Dot(aVec) + aD1Sv.Dot(aD1Sv);
2253 const Standard_Real aDet = aDf1u*aDf2v - aDf1v*aDf2u;
2254 aU -= theStep0U*(aDf2v*aF1 - aDf1v*aF2)/aDet;
2255 aV += theStep0V*(aDf2u*aF1 - aDf1u*aF2)/aDet;
2259 return (aSQDistPrev < aTol);
2262 //=======================================================================
2263 //function : SeekPointOnBoundary
2265 //=======================================================================
2266 Standard_Boolean IntWalk_PWalking::
2267 SeekPointOnBoundary(const Handle(Adaptor3d_HSurface)& theASurf1,
2268 const Handle(Adaptor3d_HSurface)& theASurf2,
2269 const Standard_Real theU1,
2270 const Standard_Real theV1,
2271 const Standard_Real theU2,
2272 const Standard_Real theV2,
2273 const Standard_Boolean isTheFirst)
2275 const Standard_Real aTol = 1.0e-14;
2276 Standard_Boolean isOK = Standard_False;
2277 Standard_Real U1prec = theU1, V1prec = theV1, U2prec = theU2, V2prec = theV2;
2279 Standard_Boolean flFinish = Standard_False;
2281 Standard_Integer aNbIter = 20;
2284 flFinish = Standard_False;
2285 Standard_Boolean aStatus = Standard_False;
2290 aStatus = DistanceMinimizeByGradient(theASurf1, theASurf2, U1prec, V1prec, U2prec, V2prec);
2296 aStatus = DistanceMinimizeByExtrema(theASurf1, theASurf2->Value(U2prec, V2prec), U1prec, V1prec);
2302 aStatus = DistanceMinimizeByExtrema(theASurf2, theASurf1->Value(U1prec, V1prec), U2prec, V2prec);
2308 while(!aStatus && (aNbIter > 0));
2312 const Standard_Real aTolMax = 1.0e-8;
2313 Standard_Real aTolF = 0.0;
2315 Standard_Real u1 = U1prec, v1 = V1prec, u2 = U2prec, v2 = V2prec;
2317 flFinish = Checking(theASurf1, theASurf2, U1prec, V1prec, U2prec, V2prec, aTolF);
2319 if(aTolF <= aTolMax)
2321 gp_Pnt aP1 = theASurf1->Value(u1, v1),
2322 aP2 = theASurf2->Value(u2, v2);
2323 gp_Pnt aPInt(0.5*(aP1.XYZ() + aP2.XYZ()));
2325 const Standard_Real aSQDist1 = aPInt.SquareDistance(aP1),
2326 aSQDist2 = aPInt.SquareDistance(aP2);
2327 if((aSQDist1 < aTol) && (aSQDist2 < aTol))
2329 IntSurf_PntOn2S anIP;
2330 anIP.SetValue(aPInt, u1, v1, u2, v2);
2333 line->InsertBefore(1,anIP);
2337 isOK = Standard_True;
2353 //=======================================================================
2354 //function : PutToBoundary
2356 //=======================================================================
2357 Standard_Boolean IntWalk_PWalking::
2358 PutToBoundary(const Handle(Adaptor3d_HSurface)& theASurf1,
2359 const Handle(Adaptor3d_HSurface)& theASurf2)
2361 const Standard_Real aTolMin = Precision::Confusion();
2363 Standard_Boolean hasBeenAdded = Standard_False;
2365 const Standard_Real aU1bFirst = theASurf1->FirstUParameter();
2366 const Standard_Real aU1bLast = theASurf1->LastUParameter();
2367 const Standard_Real aU2bFirst = theASurf2->FirstUParameter();
2368 const Standard_Real aU2bLast = theASurf2->LastUParameter();
2369 const Standard_Real aV1bFirst = theASurf1->FirstVParameter();
2370 const Standard_Real aV1bLast = theASurf1->LastVParameter();
2371 const Standard_Real aV2bFirst = theASurf2->FirstVParameter();
2372 const Standard_Real aV2bLast = theASurf2->LastVParameter();
2374 Standard_Real aTol = 1.0;
2375 aTol = Min(aTol, aU1bLast - aU1bFirst);
2376 aTol = Min(aTol, aU2bLast - aU2bFirst);
2377 aTol = Min(aTol, aV1bLast - aV1bFirst);
2378 aTol = Min(aTol, aV2bLast - aV2bFirst)*1.0e-3;
2380 if(aTol <= 2.0*aTolMin)
2381 return hasBeenAdded;
2383 Standard_Boolean isNeedAdding = Standard_False;
2384 Standard_Boolean isU1parallel = Standard_False, isV1parallel = Standard_False;
2385 Standard_Boolean isU2parallel = Standard_False, isV2parallel = Standard_False;
2386 IsParallel(line, Standard_True, aTol, isU1parallel, isV1parallel);
2387 IsParallel(line, Standard_False, aTol, isU2parallel, isV2parallel);
2389 Standard_Real u1, v1, u2, v2;
2390 line->Value(1).Parameters(u1, v1, u2, v2);
2391 Standard_Real aDelta = 0.0;
2395 aDelta = u1 - aU1bFirst;
2396 if((aTolMin < aDelta) && (aDelta < aTol))
2398 u1 = aU1bFirst - aDelta;
2399 isNeedAdding = Standard_True;
2403 aDelta = aU1bLast - u1;
2404 if((aTolMin < aDelta) && (aDelta < aTol))
2406 u1 = aU1bLast + aDelta;
2407 isNeedAdding = Standard_True;
2414 aDelta = u2 - aU2bFirst;
2415 if((aTolMin < aDelta) && (aDelta < aTol))
2417 u2 = aU2bFirst - aDelta;
2418 isNeedAdding = Standard_True;
2422 aDelta = aU2bLast - u2;
2423 if((aTolMin < aDelta) && (aDelta < aTol))
2425 u2 = aU2bLast + aDelta;
2426 isNeedAdding = Standard_True;
2433 aDelta = v1 - aV1bFirst;
2434 if((aTolMin < aDelta) && (aDelta < aTol))
2436 v1 = aV1bFirst - aDelta;
2437 isNeedAdding = Standard_True;
2441 aDelta = aV1bLast - v1;
2442 if((aTolMin < aDelta) && (aDelta < aTol))
2444 v1 = aV1bLast + aDelta;
2445 isNeedAdding = Standard_True;
2452 aDelta = v2 - aV2bFirst;
2453 if((aTolMin < aDelta) && (aDelta < aTol))
2455 v2 = aV2bFirst - aDelta;
2456 isNeedAdding = Standard_True;
2460 aDelta = aV2bLast - v2;
2461 if((aTolMin < aDelta) && (aDelta < aTol))
2463 v2 = aV2bLast + aDelta;
2464 isNeedAdding = Standard_True;
2472 SeekPointOnBoundary(theASurf1, theASurf2, u1,
2473 v1, u2, v2, Standard_True);
2476 const Standard_Integer aNbPnts = line->NbPoints();
2477 isNeedAdding = Standard_False;
2478 line->Value(aNbPnts).Parameters(u1, v1, u2, v2);
2482 aDelta = u1 - aU1bFirst;
2483 if((aTolMin < aDelta) && (aDelta < aTol))
2485 u1 = aU1bFirst - aDelta;
2486 isNeedAdding = Standard_True;
2490 aDelta = aU1bLast - u1;
2491 if((aTolMin < aDelta) && (aDelta < aTol))
2493 u1 = aU1bLast + aDelta;
2494 isNeedAdding = Standard_True;
2501 aDelta = u2 - aU2bFirst;
2502 if((aTolMin < aDelta) && (aDelta < aTol))
2504 u2 = aU2bFirst - aDelta;
2505 isNeedAdding = Standard_True;
2509 aDelta = aU2bLast - u2;
2510 if((aTolMin < aDelta) && (aDelta < aTol))
2512 u2 = aU2bLast + aDelta;
2513 isNeedAdding = Standard_True;
2520 aDelta = v1 - aV1bFirst;
2521 if((aTolMin < aDelta) && (aDelta < aTol))
2523 v1 = aV1bFirst - aDelta;
2524 isNeedAdding = Standard_True;
2528 aDelta = aV1bLast - v1;
2529 if((aTolMin < aDelta) && (aDelta < aTol))
2531 v1 = aV1bLast + aDelta;
2532 isNeedAdding = Standard_True;
2539 aDelta = v2 - aV2bFirst;
2540 if((aTolMin < aDelta) && (aDelta < aTol))
2542 v2 = aV2bFirst - aDelta;
2543 isNeedAdding = Standard_True;
2547 aDelta = aV2bLast - v2;
2548 if((aTolMin < aDelta) && (aDelta < aTol))
2550 v2 = aV2bLast + aDelta;
2551 isNeedAdding = Standard_True;
2559 SeekPointOnBoundary(theASurf1, theASurf2, u1,
2560 v1, u2, v2, Standard_False);
2563 return hasBeenAdded;
2566 //=======================================================================
2567 //function : SeekAdditionalPoints
2569 //=======================================================================
2570 Standard_Boolean IntWalk_PWalking::
2571 SeekAdditionalPoints( const Handle(Adaptor3d_HSurface)& theASurf1,
2572 const Handle(Adaptor3d_HSurface)& theASurf2,
2573 const Standard_Integer theMinNbPoints)
2575 const Standard_Real aTol = 1.0e-14;
2576 Standard_Integer aNbPoints = line->NbPoints();
2577 if(aNbPoints > theMinNbPoints)
2578 return Standard_True;
2580 const Standard_Real aU1bFirst = theASurf1->FirstUParameter();
2581 const Standard_Real aU1bLast = theASurf1->LastUParameter();
2582 const Standard_Real aU2bFirst = theASurf2->FirstUParameter();
2583 const Standard_Real aU2bLast = theASurf2->LastUParameter();
2584 const Standard_Real aV1bFirst = theASurf1->FirstVParameter();
2585 const Standard_Real aV1bLast = theASurf1->LastVParameter();
2586 const Standard_Real aV2bFirst = theASurf2->FirstVParameter();
2587 const Standard_Real aV2bLast = theASurf2->LastVParameter();
2590 Standard_Boolean isPrecise = Standard_False;
2592 Standard_Real U1prec = 0.0, V1prec = 0.0, U2prec = 0.0, V2prec = 0.0;
2594 Standard_Integer aNbPointsPrev = 0;
2595 while(aNbPoints < theMinNbPoints && (aNbPoints != aNbPointsPrev))
2597 aNbPointsPrev = aNbPoints;
2598 for(Standard_Integer fp = 1, lp = 2; fp < aNbPoints; fp = lp + 1)
2600 Standard_Real U1f, V1f, U2f, V2f; //first point in 1st and 2nd surafaces
2601 Standard_Real U1l, V1l, U2l, V2l; //last point in 1st and 2nd surafaces
2604 line->Value(fp).Parameters(U1f, V1f, U2f, V2f);
2605 line->Value(lp).Parameters(U1l, V1l, U2l, V2l);
2607 U1prec = 0.5*(U1f+U1l);
2608 if(U1prec < aU1bFirst)
2610 if(U1prec > aU1bLast)
2613 V1prec = 0.5*(V1f+V1l);
2614 if(V1prec < aV1bFirst)
2616 if(V1prec > aV1bLast)
2619 U2prec = 0.5*(U2f+U2l);
2620 if(U2prec < aU2bFirst)
2622 if(U2prec > aU2bLast)
2625 V2prec = 0.5*(V2f+V2l);
2626 if(V2prec < aV2bFirst)
2628 if(V2prec > aV2bLast)
2631 Standard_Boolean aStatus = Standard_False;
2632 Standard_Integer aNbIter = 5;
2635 aStatus = DistanceMinimizeByGradient(theASurf1, theASurf2, U1prec, V1prec, U2prec, V2prec);
2641 aStatus = DistanceMinimizeByExtrema(theASurf1, theASurf2->Value(U2prec, V2prec), U1prec, V1prec);
2647 aStatus = DistanceMinimizeByExtrema(theASurf2, theASurf1->Value(U1prec, V1prec), U2prec, V2prec);
2653 while(!aStatus && (--aNbIter > 0));
2657 gp_Pnt aP1 = theASurf1->Value(U1prec, V1prec),
2658 aP2 = theASurf2->Value(U2prec, V2prec);
2659 gp_Pnt aPInt(0.5*(aP1.XYZ() + aP2.XYZ()));
2661 const Standard_Real aSQDist1 = aPInt.SquareDistance(aP1),
2662 aSQDist2 = aPInt.SquareDistance(aP2);
2664 if((aSQDist1 < aTol) && (aSQDist2 < aTol))
2666 IntSurf_PntOn2S anIP;
2667 anIP.SetValue(aPInt, U1prec, V1prec, U2prec, V2prec);
2668 line->InsertBefore(lp, anIP);
2670 isPrecise = Standard_True;
2672 if(++aNbPoints >= theMinNbPoints)
2686 void IntWalk_PWalking::
2687 RepartirOuDiviser(Standard_Boolean& DejaReparti,
2688 IntImp_ConstIsoparametric& ChoixIso,
2689 Standard_Boolean& Arrive)
2691 // at the neighborhood of a point, there is a fail of marching
2692 // it is required to divide the steps to try to continue
2693 // if the step is too small if we are on border
2694 // restart in another direction if it was not done, otherwise stop
2697 // Standard_Integer i;
2698 if (Arrive) { //restart in the other direction
2699 if (!DejaReparti ) {
2700 Arrive = Standard_False;
2701 DejaReparti = Standard_True;
2702 previousPoint = line->Value(1);
2703 previoustg = Standard_False;
2704 previousd1 = firstd1;
2705 previousd2 = firstd2;
2707 indextg = line->NbPoints();
2711 //-- printf("\nIntWalk_PWalking_2.gxx Reverse %3d\n",indextg);
2712 sensCheminement = -1;
2714 tglast = Standard_False;
2715 ChoixIso = choixIsoSav;
2722 Standard_Real u1,v1,u2,v2;
2723 Standard_Real U1,V1,U2,V2;
2724 Standard_Integer nn=line->NbPoints();
2726 line->Value(nn).Parameters(u1,v1,u2,v2);
2727 line->Value(nn-1).Parameters(U1,V1,U2,V2);
2728 pasuv[0]=Abs(u1-U1);
2729 pasuv[1]=Abs(v1-V1);
2730 pasuv[2]=Abs(u2-U2);
2731 pasuv[3]=Abs(v2-V2);
2738 if ( pasuv[0]*0.5 < ResoU1
2739 && pasuv[1]*0.5 < ResoV1
2740 && pasuv[2]*0.5 < ResoU2
2741 && pasuv[3]*0.5 < ResoV2
2744 tglast = Standard_True; // IS IT ENOUGH ????
2747 if (!DejaReparti) { //restart in the other direction
2748 DejaReparti = Standard_True;
2749 previousPoint = line->Value(1);
2750 previoustg = Standard_False;
2751 previousd1 = firstd1;
2752 previousd2 = firstd2;
2754 indextg = line->NbPoints();
2758 //-- printf("\nIntWalk_PWalking_2.gxx Reverse %3d\n",indextg);
2760 sensCheminement = -1;
2762 tglast = Standard_False;
2763 ChoixIso = choixIsoSav;
2771 Standard_Real u1,v1,u2,v2;
2772 Standard_Real U1,V1,U2,V2;
2773 Standard_Integer nn=line->NbPoints();
2775 line->Value(nn).Parameters(u1,v1,u2,v2);
2776 line->Value(nn-1).Parameters(U1,V1,U2,V2);
2777 pasuv[0]=Abs(u1-U1);
2778 pasuv[1]=Abs(v1-V1);
2779 pasuv[2]=Abs(u2-U2);
2780 pasuv[3]=Abs(v2-V2);
2784 else Arrive = Standard_True;
2796 //OCC431(apo): modified ->
2797 static const Standard_Real CosRef2D = Cos(M_PI/9.0), AngRef2D = M_PI/2.0;
2799 static const Standard_Real d = 7.0;
2802 IntWalk_StatusDeflection IntWalk_PWalking::TestDeflection(const IntImp_ConstIsoparametric choixIso)
2804 // test if vector is observed by calculating an increase of vector
2805 // or the previous point and its tangent, the new calculated point and its
2806 // tangent; it is possible to find a cube passing by the 2 points and having as a
2807 // derivative the tangents of the intersection
2808 // calculate the point with parameter 0.5 on cube=p1
2809 // calculate the medium point of 2 points of intersection=p2
2810 // if arrow/2<=||p1p2||<= arrow consider that the vector is observed
2811 // otherwise adjust the step depending on the ratio ||p1p2||/vector
2812 // and the previous step
2813 // test if in 2 tangent planes of surfaces there is no too great angle2d
2814 // grand : if yes divide the step
2815 // test if there is no change of side
2818 if(line->NbPoints() ==1 ) {
2819 STATIC_BLOCAGE_SUR_PAS_TROP_GRAND=STATIC_PRECEDENT_INFLEXION=0;
2822 IntWalk_StatusDeflection Status = IntWalk_OK;
2823 Standard_Real FlecheCourante , Ratio = 1.0;
2826 const Handle(Adaptor3d_HSurface)& Caro1 = myIntersectionOn2S.Function().AuxillarSurface1();
2827 const Handle(Adaptor3d_HSurface)& Caro2 = myIntersectionOn2S.Function().AuxillarSurface2();
2829 const IntSurf_PntOn2S& CurrentPoint = myIntersectionOn2S.Point();
2830 //==================================================================================
2831 //========= S t o p o n p o i n t ============
2832 //==================================================================================
2833 if (myIntersectionOn2S.IsTangent()) {
2834 return IntWalk_ArretSurPoint;
2837 const gp_Dir& TgCourante = myIntersectionOn2S.Direction();
2839 const Standard_Real aCosBetweenTangent = TgCourante.Dot(previousd);
2841 //==================================================================================
2842 //========= R i s k o f i n f l e x i o n p o i n t ============
2843 //==================================================================================
2844 if (aCosBetweenTangent < 0) {
2845 //------------------------------------------------------------
2846 //-- Risk of inflexion point : Divide the step by 2
2847 //-- Initialize STATIC_PRECEDENT_INFLEXION so that
2848 //-- at the next call to return Pas_OK if there is no
2849 //-- more risk of the point of inflexion
2850 //------------------------------------------------------------
2856 STATIC_PRECEDENT_INFLEXION+=3;
2857 if (pasuv[0] < ResoU1 && pasuv[1] <ResoV1 && pasuv[2] <ResoU2 && pasuv[3] < ResoV2)
2858 return IntWalk_ArretSurPointPrecedent;
2860 return IntWalk_PasTropGrand;
2863 if(STATIC_PRECEDENT_INFLEXION > 0) {
2864 STATIC_PRECEDENT_INFLEXION -- ;
2869 //==================================================================================
2870 //========= D e t e c t c o n f u s e d P o in t s ===========
2871 //==================================================================================
2873 const Standard_Real aSqDist = previousPoint.Value().
2874 SquareDistance(CurrentPoint.Value());
2877 if (aSqDist < tolconf*tolconf) {
2878 pasInit[0] = Max(pasInit[0], 5.0*ResoU1);
2879 pasInit[1] = Max(pasInit[1], 5.0*ResoV1);
2880 pasInit[2] = Max(pasInit[2], 5.0*ResoU2);
2881 pasInit[3] = Max(pasInit[3], 5.0*ResoV2);
2883 for(Standard_Integer i = 0; i < 4; i++)
2885 pasuv[i] = Max(pasuv[i], Min(1.5*pasuv[i], pasInit[i]));
2887 //Compute local resolution: for OCC26717
2888 if (Abs(pasuv[choixIso] - pasInit[choixIso]) <= Precision::Confusion())
2890 Standard_Real CurU, CurV;
2891 if (choixIso == IntImp_UIsoparametricOnCaro1 ||
2892 choixIso == IntImp_VIsoparametricOnCaro1)
2893 previousPoint.ParametersOnS1(CurU, CurV);
2895 previousPoint.ParametersOnS2(CurU, CurV);
2896 gp_Pnt CurPnt = (choixIso == IntImp_UIsoparametricOnCaro1 ||
2897 choixIso == IntImp_VIsoparametricOnCaro1)?
2898 Adaptor3d_HSurfaceTool::Value(Caro1, CurU, CurV) :
2899 Adaptor3d_HSurfaceTool::Value(Caro2, CurU, CurV);
2903 case IntImp_UIsoparametricOnCaro1:
2905 Adaptor3d_HSurfaceTool::Value(Caro1,
2906 CurU + sensCheminement*pasuv[0],
2909 case IntImp_VIsoparametricOnCaro1:
2911 Adaptor3d_HSurfaceTool::Value(Caro1,
2913 CurV + sensCheminement*pasuv[1]);
2915 case IntImp_UIsoparametricOnCaro2:
2917 Adaptor3d_HSurfaceTool::Value(Caro2,
2918 CurU + sensCheminement*pasuv[2],
2921 case IntImp_VIsoparametricOnCaro2:
2923 Adaptor3d_HSurfaceTool::Value(Caro2,
2925 CurV + sensCheminement*pasuv[3]);
2929 Standard_Real RefDist = CurPnt.Distance(OffsetPnt);
2930 Standard_Real LocalResol = 0.;
2931 if (RefDist > gp::Resolution())
2932 LocalResol = pasuv[choixIso] * tolconf / RefDist;
2933 if (pasuv[choixIso] < 2*LocalResol)
2934 pasuv[choixIso] = pasInit[choixIso] = 2*LocalResol;
2936 ////////////////////////////////////////
2937 Status = IntWalk_PointConfondu;
2940 //==================================================================================
2941 Standard_Real Up1,Vp1,Uc1,Vc1,Du1,Dv1,AbsDu1,AbsDu2,AbsDv1,AbsDv2;
2942 Standard_Real Up2,Vp2,Uc2,Vc2,Du2,Dv2;
2944 previousPoint.Parameters(Up1,Vp1,Up2,Vp2);
2945 CurrentPoint.Parameters(Uc1,Vc1,Uc2,Vc2);
2947 Du1 = Uc1 - Up1; Dv1 = Vc1 - Vp1;
2948 Du2 = Uc2 - Up2; Dv2 = Vc2 - Vp2;
2954 //=================================================================================
2955 //==== S t e p o f p r o g r e s s i o n (between previous and Current) =======
2956 //=================================================================================
2957 if ( AbsDu1 < ResoU1 && AbsDv1 < ResoV1
2958 && AbsDu2 < ResoU2 && AbsDv2 < ResoV2) {
2959 pasuv[0] = ResoU1; pasuv[1] = ResoV1; pasuv[2] = ResoU2; pasuv[3] = ResoV2;
2960 return(IntWalk_ArretSurPointPrecedent);
2962 //==================================================================================
2964 Standard_Real tolArea = 100.0;
2965 if (ResoU1 < Precision::PConfusion() ||
2966 ResoV1 < Precision::PConfusion() ||
2967 ResoU2 < Precision::PConfusion() ||
2968 ResoV2 < Precision::PConfusion() )
2969 tolArea = tolArea*2.0;
2971 Standard_Real Cosi1, CosRef1, Ang1, AngRef1, ResoUV1, Duv1, d1, tolCoeff1;
2972 Standard_Real Cosi2, CosRef2, Ang2, AngRef2, ResoUV2, Duv2, d2, tolCoeff2;
2973 Cosi1 = Du1*previousd1.X() + Dv1*previousd1.Y();
2974 Cosi2 = Du2*previousd2.X() + Dv2*previousd2.Y();
2975 Duv1 = Du1*Du1 + Dv1*Dv1;
2976 Duv2 = Du2*Du2 + Dv2*Dv2;
2977 ResoUV1 = ResoU1*ResoU1 + ResoV1*ResoV1;
2978 ResoUV2 = ResoU2*ResoU2 + ResoV2*ResoV2;
2980 //modified by NIZNHY-PKV Wed Nov 13 12:25:44 2002 f
2982 Standard_Real aMinDiv2=Precision::Confusion();
2983 aMinDiv2=aMinDiv2*aMinDiv2;
2986 if (Duv1>aMinDiv2) {
2987 d1 = Abs(ResoUV1/Duv1);
2988 d1 = Min(Sqrt(d1)*tolArea, d);
2990 //d1 = Abs(ResoUV1/Duv1);
2991 //d1 = Min(Sqrt(d1)*tolArea,d);
2992 //modified by NIZNHY-PKV Wed Nov 13 12:34:30 2002 t
2993 tolCoeff1 = Exp(d1);
2995 //modified by NIZNHY-PKV Wed Nov 13 12:34:43 2002 f
2997 if (Duv2>aMinDiv2) {
2998 d2 = Abs(ResoUV2/Duv2);
2999 d2 = Min(Sqrt(d2)*tolArea,d);
3001 //d2 = Abs(ResoUV2/Duv2);
3002 //d2 = Min(Sqrt(d2)*tolArea,d);
3003 //modified by NIZNHY-PKV Wed Nov 13 12:34:53 2002 t
3004 tolCoeff2 = Exp(d2);
3005 CosRef1 = CosRef2D/tolCoeff1;
3006 CosRef2 = CosRef2D/tolCoeff2;
3008 //==================================================================================
3009 //== The points are not confused : ==
3010 //== D e t e c t t h e S t o p a t p r e v i o u s p o i n t ==
3011 //== N o t T o o G r e a t (angle in space UV) ==
3012 //== C h a n g e o f s i d e ==
3013 //==================================================================================
3014 if (Status != IntWalk_PointConfondu) {
3015 if(Cosi1*Cosi1 < CosRef1*Duv1 || Cosi2*Cosi2 < CosRef2*Duv2) {
3016 pasuv[0]*=0.5; pasuv[1]*=0.5; pasuv[2]*=0.5; pasuv[3]*=0.5;
3017 if (pasuv[0]<ResoU1 && pasuv[1]<ResoV1 && pasuv[2]<ResoU2 && pasuv[3]<ResoV2) {
3018 return(IntWalk_ArretSurPointPrecedent);
3021 pasuv[0]*=0.5; pasuv[1]*=0.5; pasuv[2]*=0.5; pasuv[3]*=0.5;
3022 return(IntWalk_PasTropGrand);
3025 const gp_Dir2d& Tg2dcourante1 = myIntersectionOn2S.DirectionOnS1();
3026 const gp_Dir2d& Tg2dcourante2 = myIntersectionOn2S.DirectionOnS2();
3027 Cosi1 = Du1*Tg2dcourante1.X() + Dv1*Tg2dcourante1.Y();
3028 Cosi2 = Du2*Tg2dcourante2.X() + Dv2*Tg2dcourante2.Y();
3029 Ang1 = Abs(previousd1.Angle(Tg2dcourante1));
3030 Ang2 = Abs(previousd2.Angle(Tg2dcourante2));
3031 AngRef1 = AngRef2D*tolCoeff1;
3032 AngRef2 = AngRef2D*tolCoeff2;
3033 //-------------------------------------------------------
3034 //-- Test : Angle too great in space UV -----
3035 //-- Change of side -----
3036 //-------------------------------------------------------
3037 if(Cosi1*Cosi1 < CosRef1*Duv1 || Cosi2*Cosi2 < CosRef2*Duv2 || Ang1 > AngRef1 || Ang2 > AngRef2) {
3038 pasuv[0]*=0.5; pasuv[1]*=0.5; pasuv[2]*=0.5; pasuv[3]*=0.5;
3039 if (pasuv[0]<ResoU1 && pasuv[1]<ResoV1 && pasuv[2]<ResoU2 && pasuv[3]<ResoV2)
3040 return(IntWalk_ArretSurPoint);
3042 return(IntWalk_PasTropGrand);
3046 //==================================================================================
3047 //== D e t e c t i o n o f : Step Too Small
3049 //==================================================================================
3051 //---------------------------------------
3052 //-- Estimate of the vector --
3053 //---------------------------------------
3055 Sqrt(Abs((previousd.XYZ()-TgCourante.XYZ()).SquareModulus()*aSqDist))/8.;
3057 if ( FlecheCourante<= fleche*0.5) { //-- Current step too small
3058 if(FlecheCourante>1e-16) {
3059 Ratio = 0.5*(fleche/FlecheCourante);
3064 Standard_Real pasSu1 = pasuv[0];
3065 Standard_Real pasSv1 = pasuv[1];
3066 Standard_Real pasSu2 = pasuv[2];
3067 Standard_Real pasSv2 = pasuv[3];
3070 //-- a point at U+DeltaU is required, ....
3071 //-- return a point at U + Epsilon
3072 //-- Epsilon << DeltaU.
3074 if(pasuv[0]< AbsDu1) pasuv[0] = AbsDu1;
3075 if(pasuv[1]< AbsDv1) pasuv[1] = AbsDv1;
3076 if(pasuv[2]< AbsDu2) pasuv[2] = AbsDu2;
3077 if(pasuv[3]< AbsDv2) pasuv[3] = AbsDv2;
3079 if(pasuv[0]<ResoU1) pasuv[0]=ResoU1;
3080 if(pasuv[1]<ResoV1) pasuv[1]=ResoV1;
3081 if(pasuv[2]<ResoU2) pasuv[2]=ResoU2;
3082 if(pasuv[3]<ResoV2) pasuv[3]=ResoV2;
3083 //-- if(Ratio>10.0 ) { Ratio=10.0; }
3084 Standard_Real R1,R = pasInit[0]/pasuv[0];
3085 R1= pasInit[1]/pasuv[1]; if(R1<R) R=R1;
3086 R1= pasInit[2]/pasuv[2]; if(R1<R) R=R1;
3087 R1= pasInit[3]/pasuv[3]; if(R1<R) R=R1;
3088 if(Ratio > R) Ratio=R;
3089 pasuv[0] = Min(Ratio*pasuv[0],pasInit[0]);
3090 pasuv[1] = Min(Ratio*pasuv[1],pasInit[1]);
3091 pasuv[2] = Min(Ratio*pasuv[2],pasInit[2]);
3092 pasuv[3] = Min(Ratio*pasuv[3],pasInit[3]);
3093 if (pasuv[0] != pasSu1 || pasuv[2] != pasSu2||
3094 pasuv[1] != pasSv1 || pasuv[3] != pasSv2) {
3095 if(++STATIC_BLOCAGE_SUR_PAS_TROP_GRAND > 5) {
3096 STATIC_BLOCAGE_SUR_PAS_TROP_GRAND = 0;
3097 return IntWalk_PasTropGrand;
3100 if(Status == IntWalk_OK) {
3101 STATIC_BLOCAGE_SUR_PAS_TROP_GRAND=0;
3102 //-- Try to increase the step
3106 else { //-- CurrentVector > vector*0.5
3107 if (FlecheCourante > fleche) { //-- Current step too Great
3108 Ratio = fleche/FlecheCourante;
3109 pasuv[0] = Ratio*pasuv[0];
3110 pasuv[1] = Ratio*pasuv[1];
3111 pasuv[2] = Ratio*pasuv[2];
3112 pasuv[3] = Ratio*pasuv[3];
3113 //if(++STATIC_BLOCAGE_SUR_PAS_TROP_GRAND > 5) {
3114 // STATIC_BLOCAGE_SUR_PAS_TROP_GRAND = 0;
3115 return IntWalk_PasTropGrand;
3118 else { //-- vector/2 < CurrentVector <= vector
3119 Ratio = 0.75 * (fleche / FlecheCourante);
3123 if(Status != IntWalk_PointConfondu)
3125 //Here, aCosBetweenTangent >= 0.0 definitely.
3128 Brief algorithm description.
3129 We have two (not-coincindent) intersection point (P1 and P2). In every point,
3130 vector of tangent (to the intersection curve) is known (vectors T1 and T2).
3131 Basing on these data, we create osculating circle.
3133 * - arc of osculating circle
3140 Let me draw your attention to the following facts:
3141 1. Vectors T1 and T2 direct FROM (not TO) points P1 and P2. Therefore,
3142 one of previously computed vector should be reversed.
3144 In this case, the absolute (!) value of the deflection between the arc of
3145 the osculating circle and the P1P2 segment can be computed as follows:
3146 e = d*(1-sin(B/2))/(2*cos(B/2)), (1)
3147 where d is the length of P1P2 segment, B is the angle between vectors T1 and T2.
3155 Later, the normal state of algorithm work is (as we apply)
3156 tolconf/2 <= e <= tolconf.
3157 In this case, we keep previous step.
3159 If e < tolconf/2 then the local curvature of the intersection curve is small.
3160 As result, the step should be increased.
3162 If e > tolconf then the step is too big. Therefore, we should decrease one.
3164 Condition (1) is equivalent to
3165 sin(B/2) = 1 - 2/(1+(d/(2*e))^2) = Fs(e),
3166 cos(B) = 1 - 2*Fs(e)^2 = Fd(e),
3167 where Fs(e)and Fd(e) are some function with parameter "deflection".
3169 Let mean that Fs(e) is decreasing function. Fd(e) is increasing function,
3170 in the range, where Fs(e) > 0.0 (i.e. when e < d/2).
3172 Now, let substitute required deflection (tolconf or tolconf/2) to z. Then
3173 it is necessary to check if e < z or if e > z.
3175 In this case, it is enough to comapare Fs(e) and Fs(z).
3176 At that Fs(e) > 0 because sin(B/2) > 0 always.
3178 Therefore, if Fs(z) < 0.0 then Fs(e) > Fs(z) ==> e < z definitely.
3179 If Fs(z) > 0.0 then we can compare Fs(z)^2 and Fs(e)^2 or, in substance,
3180 values Fd(e) and Fd(z). If Fd(e) > Fd(z) then e > z and vice versa.
3183 //Fd(e) is already known (Fd(e) == -aCosBetweenTangent)
3185 const Standard_Real anInvSqAbsArcDeflMax = 0.25*aSqDist/(tolconf*tolconf);
3186 const Standard_Real aSinB2Max = 1.0 - 2.0/(1.0 + anInvSqAbsArcDeflMax);
3188 if(aSinB2Max >= 0.0 && (aCosBetweenTangent <= 2.0 * aSinB2Max * aSinB2Max - 1.0))
3189 {//Real deflection is greater or equal than tolconf
3190 Status = IntWalk_PasTropGrand;
3193 {//Real deflection is less than tolconf
3194 const Standard_Real anInvSqAbsArcDeflMin = 4.0*anInvSqAbsArcDeflMax;
3195 const Standard_Real aSinB2Min = 1.0 - 2.0/(1.0 + anInvSqAbsArcDeflMin);
3197 if((aSinB2Min < 0.0) || (aCosBetweenTangent >= 2.0 * aSinB2Min * aSinB2Min - 1.0))
3198 {//Real deflection is less than tolconf/2.0
3199 Status = IntWalk_StepTooSmall;
3203 if(Status == IntWalk_PasTropGrand)
3205 pasuv[0]*=0.5; pasuv[1]*=0.5; pasuv[2]*=0.5; pasuv[3]*=0.5;
3209 if(Status == IntWalk_StepTooSmall)
3211 pasuv[0] = Max(pasuv[0], AbsDu1);
3212 pasuv[1] = Max(pasuv[1], AbsDv1);
3213 pasuv[2] = Max(pasuv[2], AbsDu2);
3214 pasuv[3] = Max(pasuv[3], AbsDv2);
3216 pasInit[0] = Max(pasInit[0], AbsDu1);
3217 pasInit[1] = Max(pasInit[1], AbsDv1);
3218 pasInit[2] = Max(pasInit[2], AbsDu2);
3219 pasInit[3] = Max(pasInit[3], AbsDv2);
3225 pasuv[0] = Max(myStepMin[0],Min(Min(Ratio*AbsDu1,pasuv[0]),pasInit[0]));
3226 pasuv[1] = Max(myStepMin[1],Min(Min(Ratio*AbsDv1,pasuv[1]),pasInit[1]));
3227 pasuv[2] = Max(myStepMin[2],Min(Min(Ratio*AbsDu2,pasuv[2]),pasInit[2]));
3228 pasuv[3] = Max(myStepMin[3],Min(Min(Ratio*AbsDv2,pasuv[3]),pasInit[3]));
3230 if(Status == IntWalk_OK) STATIC_BLOCAGE_SUR_PAS_TROP_GRAND=0;
3234 Standard_Boolean IntWalk_PWalking::
3235 TestArret(const Standard_Boolean DejaReparti,
3236 TColStd_Array1OfReal& Param,
3237 IntImp_ConstIsoparametric& ChoixIso)
3240 // test if the point of intersection set by these parameters remains in the
3241 // natural domain of each square.
3242 // if the point outpasses reframe to find the best iso (border)
3243 // that intersects easiest the other square
3244 // otherwise test if closed line is present
3247 Standard_Real Uvd[4],Uvf[4],Epsuv[4],Duv[4],Uvp[4],dv,dv2,ParC[4];
3248 Standard_Real DPc,DPb;
3249 Standard_Integer i = 0, k = 0;
3254 previousPoint.Parameters(Uvp[0],Uvp[1],Uvp[2],Uvp[3]);
3256 Standard_Real SolParam[4];
3257 myIntersectionOn2S.Point().Parameters(SolParam[0],SolParam[1],SolParam[2],SolParam[3]);
3259 Standard_Boolean Trouve = Standard_False;
3261 Uvd[0]=Um1; Uvf[0]=UM1; Uvd[1]=Vm1; Uvf[1]=VM1;
3262 Uvd[2]=Um2; Uvf[2]=UM2; Uvd[3]=Vm2; Uvf[3]=VM2;
3264 Standard_Integer im1;
3265 for ( i = 1,im1 = 0;i<=4;i++,im1++) {
3272 if (Param(i) < (Uvd[im1]-Epsuv[im1]) ||
3273 SolParam[im1] < (Uvd[im1]-Epsuv[im1])) //-- Current ----- Bound Inf ----- Previous
3275 Trouve = Standard_True; //--
3276 DPc = Uvp[im1]-Param(i); //-- Previous - Current
3277 DPb = Uvp[im1]-Uvd[im1]; //-- Previous - Bound Inf
3278 ParC[im1] = Uvd[im1]; //-- ParamCorrige
3279 dv = Param(k)-Uvp[k-1]; //-- Current - Previous (other Direction)
3281 if(dv2>RealEpsilon()) { //-- Progress at the other Direction ?
3282 Duv[im1] = DPc*DPb + dv2;
3283 Duv[im1] = Duv[im1]*Duv[im1]/(DPc*DPc+dv2)/(DPb*DPb+dv2);
3286 Duv[im1]=-1.0; //-- If no progress, do not change
3287 } //-- the choice of iso
3289 else if (Param(i) > (Uvf[im1] + Epsuv[im1]) ||
3290 SolParam[im1] > (Uvf[im1] + Epsuv[im1]))//-- Previous ----- Bound Sup ----- Current
3292 Trouve = Standard_True; //--
3293 DPc = Param(i)-Uvp[im1]; //-- Current - Previous
3294 DPb = Uvf[im1]-Uvp[im1]; //-- Bound Sup - Previous
3295 ParC[im1] = Uvf[im1]; //-- Param Corrige
3296 dv = Param(k)-Uvp[k-1]; //-- Current - Previous (other Direction)
3298 if(dv2>RealEpsilon()) { //-- Progress in other Direction ?
3299 Duv[im1] = DPc*DPb + dv2;
3300 Duv[im1] = Duv[im1]*Duv[im1]/(DPc*DPc+dv2)/(DPb*DPb+dv2);
3303 Duv[im1]=-1.0; //-- If no progress, do not change
3304 } //-- the choice of iso
3313 //--------------------------------------------------
3314 //-- One of Parameters u1,v1,u2,v2 is outside of --
3315 //-- the natural limits. --
3316 //-- Find the best direction of --
3317 //-- progress and reframe the parameters. --
3318 //--------------------------------------------------
3319 Standard_Real ddv = -1.0;
3321 for (i=0;i<=3;i++) {
3322 Param(i+1) = ParC[i];
3329 ChoixIso = ChoixRef[k];
3332 if((ParC[0]<=Uvd[0]+Epsuv[0]) || (ParC[0]>=Uvf[0]-Epsuv[0])) {
3333 ChoixIso = IntImp_UIsoparametricOnCaro1;
3335 else if((ParC[1]<=Uvd[1]+Epsuv[1]) || (ParC[1]>=Uvf[1]-Epsuv[1])) {
3336 ChoixIso = IntImp_VIsoparametricOnCaro1;
3338 else if((ParC[2]<=Uvd[2]+Epsuv[2]) || (ParC[2]>=Uvf[2]-Epsuv[2])) {
3339 ChoixIso = IntImp_UIsoparametricOnCaro2;
3341 else if((ParC[3]<=Uvd[3]+Epsuv[3]) || (ParC[3]>=Uvf[3]-Epsuv[3])) {
3342 ChoixIso = IntImp_VIsoparametricOnCaro2;
3345 close = Standard_False;
3346 return Standard_True;
3350 if (!DejaReparti) { // find if line closed
3353 const IntSurf_PntOn2S& POn2S1=line->Value(1);
3355 POn2S1.ParametersOnS1(u,v);
3356 gp_Pnt2d P1uvS1(u,v);
3357 previousPoint.ParametersOnS1(u,v);
3358 gp_Pnt2d PrevuvS1(u,v);
3359 myIntersectionOn2S.Point().ParametersOnS1(u,v);
3360 gp_Pnt2d myIntersuvS1(u,v);
3361 Standard_Boolean close2dS1 = (P1uvS1.XY()-PrevuvS1.XY())*
3362 (P1uvS1.XY()-myIntersuvS1.XY()) < 0.0;
3364 POn2S1.ParametersOnS2(u,v);
3365 gp_Pnt2d P1uvS2(u,v);
3366 previousPoint.ParametersOnS2(u,v);
3367 gp_Pnt2d PrevuvS2(u,v);
3368 myIntersectionOn2S.Point().ParametersOnS2(u,v);
3369 gp_Pnt2d myIntersuvS2(u,v);
3370 Standard_Boolean close2dS2 = (P1uvS2.XY()-PrevuvS2.XY())*
3371 (P1uvS2.XY()-myIntersuvS2.XY()) < 0.0;
3373 close = close2dS1 && close2dS2;
3376 else return Standard_False;