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