0031668: Visualization - WebGL sample doesn't work on Emscripten 1.39
[occt.git] / src / ApproxInt / ApproxInt_ImpPrmSvSurfaces.gxx
CommitLineData
b311480e 1// Created on: 1993-03-17
2// Created by: Laurent BUCHARD
3// Copyright (c) 1993-1999 Matra Datavision
973c2be1 4// Copyright (c) 1999-2014 OPEN CASCADE SAS
b311480e 5//
973c2be1 6// This file is part of Open CASCADE Technology software library.
b311480e 7//
d5f74e42 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
973c2be1 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.
b311480e 13//
973c2be1 14// Alternatively, this file may be used under the terms of Open CASCADE
15// commercial license or contractual agreement.
7fd59977 16
2c26a53d 17#include <IntSurf_PntOn2S.hxx>
7fd59977 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
4e14c88f 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);
7fd59977 37
4e14c88f 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}
7fd59977 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;
7fd59977 360 this->Compute(tu1,tv1,tu2,tv2,aP,aT,aTS1,aTS2);
7fd59977 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}
4e14c88f 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
7fd59977 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
2c26a53d 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;
7fd59977 497 Standard_Real binfu,bsupu,binfv,bsupv;
4e14c88f 498 binfu = ThePSurfaceTool::FirstUParameter(aPSurf);
499 binfv = ThePSurfaceTool::FirstVParameter(aPSurf);
500 bsupu = ThePSurfaceTool::LastUParameter(aPSurf);
501 bsupv = ThePSurfaceTool::LastVParameter(aPSurf);
7fd59977 502 BornInf(1) = binfu; BornSup(1) = bsupu;
503 BornInf(2) = binfv; BornSup(2) = bsupv;
2c26a53d 504 Standard_Real TranslationU = 0., TranslationV = 0.;
7fd59977 505
2c26a53d 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 }
7fd59977 515
516 Standard_Real PourTesterU = X(1);
517 Standard_Real PourTesterV = X(2);
2c26a53d 518
519 math_FunctionSetRoot Rsnld(MyZerImpFunc);
520 Rsnld.SetTolerance(Tolerance);
7fd59977 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));
7fd59977 528
4e14c88f 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);
7fd59977 558 }
4e14c88f 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);
7fd59977 580 }
4e14c88f 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;
7fd59977 606 }
4e14c88f 607 else
608 {
609 aNormalPrm.Divide(sqrt(aSQMagnPrm));
7fd59977 610 }
4e14c88f 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));
7fd59977 630 }
631
4e14c88f 632 if(isImpSingular && isPrmSingular)
633 {
634 //All is OK. All abnormal cases were processed above.
7fd59977 635
4e14c88f 636 MyTguv1 = Tguv1;
637 MyTguv2 = Tguv2;
638
639 MyIsTangent=Standard_True;
640 return MyIsTangent;
7fd59977 641 }
4e14c88f 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 }
7fd59977 689 }
4e14c88f 690
691 MyTguv1 = Tguv1;
692 MyTguv2 = Tguv2;
693
7fd59977 694 MyIsTangent=Standard_True;
4e14c88f 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;
7fd59977 711 }
712 else {
7fd59977 713 //-- cout<<" ApproxInt_ImpImpSvSurfaces.gxx : Distance apres recadrage Trop Grande "<<endl;
714
715 MyIsTangent=Standard_False;
716 MyHasBeenComputed = MyHasBeenComputedbis = Standard_False;
4e14c88f 717 return Standard_False;
7fd59977 718 }
719 }
720 else {
721 MyIsTangent = Standard_False;
722 MyHasBeenComputed = MyHasBeenComputedbis = Standard_False;
4e14c88f 723 return Standard_False;
7fd59977 724 }
725}
726
2c26a53d 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();
7fd59977 742
2c26a53d 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//--------------------------------------------------------------------------------
7fd59977 809
2c26a53d 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();
7fd59977 824
2c26a53d 825 math_Vector F(1,1);
7fd59977 826
2c26a53d 827 TranslationU = 0.0;
828 TranslationV = 0.0;
7fd59977 829
2c26a53d 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}