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