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