0025019: Command "bsection" in Test Harness with flag build pcurve on second shape...
[occt.git] / src / IntTools / IntTools_FaceFace.cxx
1 // Created on: 2000-11-23
2 // Created by: Michael KLOKOV
3 // Copyright (c) 2000-2014 OPEN CASCADE SAS
4 //
5 // This file is part of Open CASCADE Technology software library.
6 //
7 // This library is free software; you can redistribute it and/or modify it under
8 // the terms of the GNU Lesser General Public License version 2.1 as published
9 // by the Free Software Foundation, with special exception defined in the file
10 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
11 // distribution for complete text of the license and disclaimer of any warranty.
12 //
13 // Alternatively, this file may be used under the terms of Open CASCADE
14 // commercial license or contractual agreement.
15
16 #include <IntTools_FaceFace.ixx>
17
18 #include <Precision.hxx>
19
20 #include <TColStd_HArray1OfReal.hxx>
21 #include <TColStd_Array1OfReal.hxx>
22 #include <TColStd_Array1OfInteger.hxx>
23 #include <TColStd_SequenceOfReal.hxx>
24 #include <TColStd_ListOfInteger.hxx>
25 #include <TColStd_ListIteratorOfListOfInteger.hxx>
26 #include <TColStd_Array1OfListOfInteger.hxx>
27
28 #include <gp_Lin2d.hxx>
29 #include <gp_Ax22d.hxx>
30 #include <gp_Circ2d.hxx>
31 #include <gp_Torus.hxx>
32 #include <gp_Cylinder.hxx>
33
34 #include <Bnd_Box.hxx>
35
36 #include <TColgp_HArray1OfPnt2d.hxx>
37 #include <TColgp_SequenceOfPnt2d.hxx>
38 #include <TColgp_Array1OfPnt.hxx>
39 #include <TColgp_Array1OfPnt2d.hxx>
40
41 #include <IntAna_QuadQuadGeo.hxx>
42
43 #include <IntSurf_PntOn2S.hxx>
44 #include <IntSurf_LineOn2S.hxx>
45 #include <IntSurf_PntOn2S.hxx>
46 #include <IntSurf_ListOfPntOn2S.hxx>
47 #include <IntRes2d_Domain.hxx>
48 #include <ProjLib_Plane.hxx>
49
50 #include <IntPatch_GLine.hxx>
51 #include <IntPatch_RLine.hxx>
52 #include <IntPatch_WLine.hxx>
53 #include <IntPatch_ALine.hxx>
54 #include <IntPatch_ALineToWLine.hxx>
55
56 #include <ElSLib.hxx>
57 #include <ElCLib.hxx>
58
59 #include <Extrema_ExtCC.hxx>
60 #include <Extrema_POnCurv.hxx>
61 #include <BndLib_AddSurface.hxx>
62
63 #include <Adaptor3d_SurfacePtr.hxx>
64 #include <Adaptor2d_HLine2d.hxx>
65
66 #include <GeomAbs_SurfaceType.hxx>
67 #include <GeomAbs_CurveType.hxx>
68
69 #include <Geom_Surface.hxx>
70 #include <Geom_Line.hxx>
71 #include <Geom_Circle.hxx>
72 #include <Geom_Ellipse.hxx>
73 #include <Geom_Parabola.hxx>
74 #include <Geom_Hyperbola.hxx>
75 #include <Geom_TrimmedCurve.hxx>
76 #include <Geom_BSplineCurve.hxx>
77 #include <Geom_RectangularTrimmedSurface.hxx>
78 #include <Geom_OffsetSurface.hxx>
79 #include <Geom_Curve.hxx>
80 #include <Geom_Conic.hxx>
81
82 #include <Geom2d_TrimmedCurve.hxx>
83 #include <Geom2d_BSplineCurve.hxx>
84 #include <Geom2d_Line.hxx>
85 #include <Geom2d_Curve.hxx>
86 #include <Geom2d_Circle.hxx>
87
88 #include <Geom2dAPI_InterCurveCurve.hxx>
89 #include <Geom2dInt_GInter.hxx>
90 #include <GeomAdaptor_Curve.hxx>
91 #include <GeomAdaptor_HSurface.hxx>
92 #include <GeomAdaptor_Surface.hxx>
93 #include <GeomLib_CheckBSplineCurve.hxx>
94 #include <GeomLib_Check2dBSplineCurve.hxx>
95
96 #include <GeomInt_WLApprox.hxx>
97 #include <GeomProjLib.hxx>
98 #include <GeomAPI_ProjectPointOnSurf.hxx>
99 #include <Geom2dAdaptor_Curve.hxx>
100 #include <TopoDS.hxx>
101 #include <TopoDS_Edge.hxx>
102 #include <TopExp_Explorer.hxx>
103
104 #include <BRep_Tool.hxx>
105 #include <BRepTools.hxx>
106 #include <BRepAdaptor_Surface.hxx>
107
108 #include <IntTools_Curve.hxx>
109 #include <IntTools_Tools.hxx>
110 #include <IntTools_Tools.hxx>
111 #include <IntTools_TopolTool.hxx>
112 #include <IntTools_PntOnFace.hxx>
113 #include <IntTools_PntOn2Faces.hxx>
114 #include <IntTools_Context.hxx>
115 #include <IntSurf_ListIteratorOfListOfPntOn2S.hxx>
116
117 static
118   void RefineVector(gp_Vec2d& aV2D);
119 #ifdef DEB_DUMPWLINE
120 static
121   void DumpWLine(const Handle(IntPatch_WLine)& aWLine);
122 #endif
123 //
124 static
125   void TolR3d(const TopoDS_Face& ,
126               const TopoDS_Face& ,
127               Standard_Real& );
128 static 
129   Handle(Geom_Curve) MakeBSpline (const Handle(IntPatch_WLine)&,
130                                   const Standard_Integer,
131                                   const Standard_Integer);
132
133 static 
134   void Parameters(const Handle(GeomAdaptor_HSurface)&,
135                   const Handle(GeomAdaptor_HSurface)&,
136                   const gp_Pnt&,
137                   Standard_Real&,
138                   Standard_Real&,
139                   Standard_Real&,
140                   Standard_Real&);
141
142 static 
143   void BuildPCurves (Standard_Real f,Standard_Real l,Standard_Real& Tol,
144                      const Handle (Geom_Surface)& S,
145                      const Handle (Geom_Curve)&   C,
146                      Handle (Geom2d_Curve)& C2d);
147
148 static 
149   void CorrectSurfaceBoundaries(const TopoDS_Face&  theFace,
150                                 const Standard_Real theTolerance,
151                                 Standard_Real&      theumin,
152                                 Standard_Real&      theumax, 
153                                 Standard_Real&      thevmin, 
154                                 Standard_Real&      thevmax);
155 static
156   Standard_Boolean NotUseSurfacesForApprox
157           (const TopoDS_Face& aF1,
158            const TopoDS_Face& aF2,
159            const Handle(IntPatch_WLine)& WL,
160            const Standard_Integer ifprm,
161            const Standard_Integer ilprm);
162
163 static 
164   Handle(IntPatch_WLine) ComputePurgedWLine(const Handle(IntPatch_WLine)& theWLine);
165
166 static 
167   Handle(Geom2d_BSplineCurve) MakeBSpline2d(const Handle(IntPatch_WLine)& theWLine,
168                                             const Standard_Integer                         ideb,
169                                             const Standard_Integer                         ifin,
170                                             const Standard_Boolean                         onFirst);
171
172 static 
173   Standard_Boolean DecompositionOfWLine(const Handle(IntPatch_WLine)& theWLine,
174                                         const Handle(GeomAdaptor_HSurface)&            theSurface1, 
175                                         const Handle(GeomAdaptor_HSurface)&            theSurface2,
176                                         const TopoDS_Face&                             theFace1,
177                                         const TopoDS_Face&                             theFace2,
178                                         const IntTools_LineConstructor&                theLConstructor,
179                                         const Standard_Boolean                         theAvoidLConstructor,
180                                         IntPatch_SequenceOfLine&                       theNewLines,
181                                         Standard_Real&                                 theReachedTol3d,
182                                         const Handle(IntTools_Context)& );
183
184 static 
185   Standard_Boolean ParameterOutOfBoundary(const Standard_Real       theParameter, 
186                                           const Handle(Geom_Curve)& theCurve, 
187                                           const TopoDS_Face&        theFace1, 
188                                           const TopoDS_Face&        theFace2,
189                                           const Standard_Real       theOtherParameter,
190                                           const Standard_Boolean    bIncreasePar,
191                                           Standard_Real&            theNewParameter,
192                                           const Handle(IntTools_Context)& );
193
194 static 
195   Standard_Boolean IsCurveValid(Handle(Geom2d_Curve)& thePCurve);
196
197 static 
198   Standard_Boolean IsPointOnBoundary(const Standard_Real theParameter,
199                                      const Standard_Real theFirstBoundary,
200                                      const Standard_Real theSecondBoundary,
201                                      const Standard_Real theResolution,
202                                      Standard_Boolean&   IsOnFirstBoundary);
203 static
204   Standard_Boolean FindPoint(const gp_Pnt2d&     theFirstPoint,
205                              const gp_Pnt2d&     theLastPoint,
206                              const Standard_Real theUmin, 
207                              const Standard_Real theUmax,
208                              const Standard_Real theVmin,
209                              const Standard_Real theVmax,
210                              gp_Pnt2d&           theNewPoint);
211
212
213 static 
214   Standard_Integer ComputeTangentZones( const Handle(GeomAdaptor_HSurface)& theSurface1,
215                                        const Handle(GeomAdaptor_HSurface)&  theSurface2,
216                                        const TopoDS_Face&                   theFace1,
217                                        const TopoDS_Face&                   theFace2,
218                                        Handle(TColgp_HArray1OfPnt2d)&       theResultOnS1,
219                                        Handle(TColgp_HArray1OfPnt2d)&       theResultOnS2,
220                                        Handle(TColStd_HArray1OfReal)&       theResultRadius,
221                                        const Handle(IntTools_Context)& );
222
223 static
224   Standard_Boolean FindPoint(const gp_Pnt2d&     theFirstPoint,
225                              const gp_Pnt2d&     theLastPoint,
226                              const Standard_Real theUmin, 
227                              const Standard_Real theUmax,
228                              const Standard_Real theVmin,
229                              const Standard_Real theVmax,
230                              const gp_Pnt2d&     theTanZoneCenter,
231                              const Standard_Real theZoneRadius,
232                              Handle(GeomAdaptor_HSurface) theGASurface,
233                              gp_Pnt2d&           theNewPoint);
234
235 static
236   Standard_Boolean IsInsideTanZone(const gp_Pnt2d&     thePoint,
237                                    const gp_Pnt2d&     theTanZoneCenter,
238                                    const Standard_Real theZoneRadius,
239                                    Handle(GeomAdaptor_HSurface) theGASurface);
240
241 static
242   gp_Pnt2d AdjustByNeighbour(const gp_Pnt2d&     theaNeighbourPoint,
243                              const gp_Pnt2d&     theOriginalPoint,
244                              Handle(GeomAdaptor_HSurface) theGASurface);
245 static
246   Standard_Boolean  ApproxWithPCurves(const gp_Cylinder& theCyl, 
247                                       const gp_Sphere& theSph);
248
249 static void  PerformPlanes(const Handle(GeomAdaptor_HSurface)& theS1, 
250                            const Handle(GeomAdaptor_HSurface)& theS2, 
251                            const Standard_Real TolAng, 
252                            const Standard_Real TolTang, 
253                            const Standard_Boolean theApprox1,
254                            const Standard_Boolean theApprox2,
255                            IntTools_SequenceOfCurves& theSeqOfCurve, 
256                            Standard_Boolean& theTangentFaces);
257
258 static Standard_Boolean ClassifyLin2d(const Handle(GeomAdaptor_HSurface)& theS, 
259                                       const gp_Lin2d& theLin2d, 
260                                       const Standard_Real theTol,
261                                       Standard_Real& theP1, 
262                                       Standard_Real& theP2);
263 //
264 static
265   void ApproxParameters(const Handle(GeomAdaptor_HSurface)& aHS1,
266                         const Handle(GeomAdaptor_HSurface)& aHS2,
267                         Standard_Integer& iDegMin,
268                         Standard_Integer& iNbIter,
269                         Standard_Integer& iDegMax);
270
271 static
272   void Tolerances(const Handle(GeomAdaptor_HSurface)& aHS1,
273                   const Handle(GeomAdaptor_HSurface)& aHS2,
274                   Standard_Real& aTolTang);
275
276 static
277   Standard_Boolean SortTypes(const GeomAbs_SurfaceType aType1,
278                              const GeomAbs_SurfaceType aType2);
279 static
280   Standard_Integer IndexType(const GeomAbs_SurfaceType aType);
281
282 //
283 static
284   Standard_Real MaxSquareDistance (const Standard_Real aT,
285                                    const Handle(Geom_Curve)& aC3D,
286                                    const Handle(Geom2d_Curve)& aC2D1,
287                                    const Handle(Geom2d_Curve)& aC2D2,
288                                    const Handle(GeomAdaptor_HSurface) myHS1,
289                                    const Handle(GeomAdaptor_HSurface) myHS2,
290                                    const TopoDS_Face& aF1,
291                                    const TopoDS_Face& aF2,
292                                    const Handle(IntTools_Context)& aCtx);
293
294 static
295   Standard_Boolean CheckPCurve(const Handle(Geom2d_Curve)& aPC, 
296                                const TopoDS_Face& aFace);
297
298 //
299 static
300   Standard_Real FindMaxSquareDistance (const Standard_Real aA,
301                                        const Standard_Real aB,
302                                        const Standard_Real aEps,
303                                        const Handle(Geom_Curve)& aC3D,
304                                        const Handle(Geom2d_Curve)& aC2D1,
305                                        const Handle(Geom2d_Curve)& aC2D2,
306                                        const Handle(GeomAdaptor_HSurface)& myHS1,
307                                        const Handle(GeomAdaptor_HSurface)& myHS2,
308                                        const TopoDS_Face& aF1,
309                                        const TopoDS_Face& aF2,
310                                        const Handle(IntTools_Context)& aCtx);
311
312 //=======================================================================
313 //function : 
314 //purpose  : 
315 //=======================================================================
316 IntTools_FaceFace::IntTools_FaceFace()
317 {
318   myIsDone=Standard_False;
319   myTangentFaces=Standard_False;
320   //
321   myHS1 = new GeomAdaptor_HSurface ();
322   myHS2 = new GeomAdaptor_HSurface ();
323   myTolReached2d=0.; 
324   myTolReached3d=0.;
325   SetParameters(Standard_True, Standard_True, Standard_True, 1.e-07);
326   
327 }
328 //=======================================================================
329 //function : SetContext
330 //purpose  : 
331 //======================================================================= 
332 void IntTools_FaceFace::SetContext(const Handle(IntTools_Context)& aContext)
333 {
334   myContext=aContext;
335 }
336 //=======================================================================
337 //function : Context
338 //purpose  : 
339 //======================================================================= 
340 const Handle(IntTools_Context)& IntTools_FaceFace::Context()const
341 {
342   return myContext;
343 }
344 //=======================================================================
345 //function : Face1
346 //purpose  : 
347 //======================================================================= 
348 const TopoDS_Face& IntTools_FaceFace::Face1() const
349 {
350   return myFace1;
351 }
352 //=======================================================================
353 //function : Face2
354 //purpose  : 
355 //======================================================================= 
356 const TopoDS_Face& IntTools_FaceFace::Face2() const
357 {
358   return myFace2;
359 }
360 //=======================================================================
361 //function : TangentFaces
362 //purpose  : 
363 //======================================================================= 
364 Standard_Boolean IntTools_FaceFace::TangentFaces() const
365 {
366   return myTangentFaces;
367 }
368 //=======================================================================
369 //function : Points
370 //purpose  : 
371 //======================================================================= 
372 const IntTools_SequenceOfPntOn2Faces& IntTools_FaceFace::Points() const
373 {
374   return myPnts;
375 }
376 //=======================================================================
377 //function : IsDone
378 //purpose  : 
379 //======================================================================= 
380 Standard_Boolean IntTools_FaceFace::IsDone() const
381 {
382   return myIsDone;
383 }
384 //=======================================================================
385 //function : TolReached3d
386 //purpose  : 
387 //=======================================================================
388 Standard_Real IntTools_FaceFace::TolReached3d() const
389 {
390   return myTolReached3d;
391 }
392 //=======================================================================
393 //function : Lines
394 //purpose  : return lines of intersection
395 //=======================================================================
396 const IntTools_SequenceOfCurves& IntTools_FaceFace::Lines() const
397 {
398   StdFail_NotDone_Raise_if
399     (!myIsDone,
400      "IntTools_FaceFace::Lines() => myIntersector NOT DONE");
401   return mySeqOfCurve;
402 }
403 //=======================================================================
404 //function : TolReached2d
405 //purpose  : 
406 //=======================================================================
407 Standard_Real IntTools_FaceFace::TolReached2d() const
408 {
409   return myTolReached2d;
410 }
411 // =======================================================================
412 // function: SetParameters
413 //
414 // =======================================================================
415 void IntTools_FaceFace::SetParameters(const Standard_Boolean ToApproxC3d,
416                                       const Standard_Boolean ToApproxC2dOnS1,
417                                       const Standard_Boolean ToApproxC2dOnS2,
418                                       const Standard_Real ApproximationTolerance) 
419 {
420   myApprox = ToApproxC3d;
421   myApprox1 = ToApproxC2dOnS1;
422   myApprox2 = ToApproxC2dOnS2;
423   myTolApprox = ApproximationTolerance;
424 }
425 //=======================================================================
426 //function : SetList
427 //purpose  : 
428 //=======================================================================
429 void IntTools_FaceFace::SetList(IntSurf_ListOfPntOn2S& aListOfPnts)
430 {
431   myListOfPnts = aListOfPnts;  
432 }
433
434
435 static Standard_Boolean isTreatAnalityc(const TopoDS_Face& theF1,
436                                         const TopoDS_Face& theF2)
437 {
438   const Standard_Real Tolang = 1.e-8;
439   const Standard_Real aTolF1=BRep_Tool::Tolerance(theF1);
440   const Standard_Real aTolF2=BRep_Tool::Tolerance(theF2);
441   const Standard_Real aTolSum = aTolF1 + aTolF2;
442   Standard_Real aHigh = 0.0;
443
444   const BRepAdaptor_Surface aBAS1(theF1), aBAS2(theF2);
445   const GeomAbs_SurfaceType aType1=aBAS1.GetType();
446   const GeomAbs_SurfaceType aType2=aBAS2.GetType();
447   
448   gp_Pln aS1;
449   gp_Cylinder aS2;
450   if(aType1 == GeomAbs_Plane)
451   {
452     aS1=aBAS1.Plane();
453   }
454   else if(aType2 == GeomAbs_Plane)
455   {
456     aS1=aBAS2.Plane();
457   }
458   else
459   {
460     return Standard_True;
461   }
462
463   if(aType1 == GeomAbs_Cylinder)
464   {
465     aS2=aBAS1.Cylinder();
466     const Standard_Real VMin = aBAS1.FirstVParameter();
467     const Standard_Real VMax = aBAS1.LastVParameter();
468
469     if( Precision::IsNegativeInfinite(VMin) ||
470         Precision::IsPositiveInfinite(VMax))
471           return Standard_True;
472     else
473       aHigh = VMax - VMin;
474   }
475   else if(aType2 == GeomAbs_Cylinder)
476   {
477     aS2=aBAS2.Cylinder();
478
479     const Standard_Real VMin = aBAS2.FirstVParameter();
480     const Standard_Real VMax = aBAS2.LastVParameter();
481
482     if( Precision::IsNegativeInfinite(VMin) ||
483         Precision::IsPositiveInfinite(VMax))
484           return Standard_True;
485     else
486       aHigh = VMax - VMin;
487   }
488   else
489   {
490     return Standard_True;
491   }
492
493   IntAna_QuadQuadGeo inter;
494   inter.Perform(aS1,aS2,Tolang,aTolSum, aHigh);
495   if(inter.TypeInter() == IntAna_Ellipse)
496   {
497     const gp_Elips anEl = inter.Ellipse(1);
498     const Standard_Real aMajorR = anEl.MajorRadius();
499     const Standard_Real aMinorR = anEl.MinorRadius();
500     
501     return (aMajorR < 100000.0 * aMinorR);
502   }
503   else
504   {
505     return inter.IsDone();
506   }
507 }
508 //=======================================================================
509 //function : Perform
510 //purpose  : intersect surfaces of the faces
511 //=======================================================================
512   void IntTools_FaceFace::Perform(const TopoDS_Face& aF1,
513                                 const TopoDS_Face& aF2)
514 {
515   Standard_Boolean RestrictLine = Standard_False, hasCone = Standard_False;
516   
517   if (myContext.IsNull()) {
518     myContext=new IntTools_Context;
519   }
520
521   mySeqOfCurve.Clear();
522   myTolReached2d=0.;
523   myTolReached3d=0.;
524   myIsDone = Standard_False;
525   myNbrestr=0;//?
526
527   myFace1=aF1;
528   myFace2=aF2;
529
530   const BRepAdaptor_Surface aBAS1(myFace1, Standard_False);
531   const BRepAdaptor_Surface aBAS2(myFace2, Standard_False);
532   GeomAbs_SurfaceType aType1=aBAS1.GetType();
533   GeomAbs_SurfaceType aType2=aBAS2.GetType();
534
535   const Standard_Boolean bReverse=SortTypes(aType1, aType2);
536   if (bReverse)
537   {
538     myFace1=aF2;
539     myFace2=aF1;
540     aType1=aBAS2.GetType();
541     aType2=aBAS1.GetType();
542
543     if (myListOfPnts.Extent())
544     {
545       Standard_Real aU1,aV1,aU2,aV2;
546       IntSurf_ListIteratorOfListOfPntOn2S aItP2S;
547       //
548       aItP2S.Initialize(myListOfPnts);
549       for (; aItP2S.More(); aItP2S.Next())
550       {
551         IntSurf_PntOn2S& aP2S=aItP2S.Value();
552         aP2S.Parameters(aU1,aV1,aU2,aV2);
553         aP2S.SetValue(aU2,aV2,aU1,aV1);
554       }
555     }
556     //
557     Standard_Boolean anAproxTmp = myApprox1;
558     myApprox1 = myApprox2;
559     myApprox2 = anAproxTmp;
560   }
561
562
563   const Handle(Geom_Surface) S1=BRep_Tool::Surface(myFace1);
564   const Handle(Geom_Surface) S2=BRep_Tool::Surface(myFace2);
565
566   const Standard_Real aTolF1=BRep_Tool::Tolerance(myFace1);
567   const Standard_Real aTolF2=BRep_Tool::Tolerance(myFace2);
568
569   Standard_Real TolArc = aTolF1 + aTolF2;
570   Standard_Real TolTang = TolArc;
571
572   const Standard_Boolean isFace1Quad = (aType1 == GeomAbs_Cylinder ||
573                                         aType1 == GeomAbs_Cone ||
574                                         aType1 == GeomAbs_Torus);
575
576   const Standard_Boolean isFace2Quad = (aType2 == GeomAbs_Cylinder ||
577                                         aType2 == GeomAbs_Cone ||
578                                         aType2 == GeomAbs_Torus);
579
580   if(aType1==GeomAbs_Plane && aType2==GeomAbs_Plane)  {
581     Standard_Real umin, umax, vmin, vmax;
582     Standard_Real dU, dV;
583     //
584     BRepTools::UVBounds(myFace1, umin, umax, vmin, vmax);
585     dU=0.1*(umax-umin);
586     dV=0.1*(vmax-vmin);
587     umin=umin-dU;
588     umax=umax+dU;
589     vmin=vmin-dV;
590     vmax=vmax+dV;
591     myHS1->ChangeSurface().Load(S1, umin, umax, vmin, vmax);
592     //
593     BRepTools::UVBounds(myFace2, umin, umax, vmin, vmax);
594     dU=0.1*(umax-umin);
595     dV=0.1*(vmax-vmin);
596     umin=umin-dU;
597     umax=umax+dU;
598     vmin=vmin-dV;
599     vmax=vmax+dV;
600     myHS2->ChangeSurface().Load(S2, umin, umax, vmin, vmax);
601     //
602     Standard_Real TolAng = 1.e-8;
603     //
604     PerformPlanes(myHS1, myHS2, TolAng, TolTang, myApprox1, myApprox2, 
605                   mySeqOfCurve, myTangentFaces);
606     //
607     myIsDone = Standard_True;
608     
609     if(!myTangentFaces) {
610       const Standard_Integer NbLinPP = mySeqOfCurve.Length();
611       if(NbLinPP) {
612         Standard_Real aTolFMax;
613         myTolReached3d = 1.e-7;
614         aTolFMax=Max(aTolF1, aTolF2);
615         if (aTolFMax>myTolReached3d) {
616           myTolReached3d=aTolFMax;
617         }
618         //
619         myTolReached2d = myTolReached3d;
620
621         if (bReverse) {
622           Handle(Geom2d_Curve) aC2D1, aC2D2;
623           const Standard_Integer aNbLin = mySeqOfCurve.Length();
624           for (Standard_Integer i = 1; i <= aNbLin; ++i) {
625             IntTools_Curve& aIC=mySeqOfCurve(i);
626             aC2D1=aIC.FirstCurve2d();
627             aC2D2=aIC.SecondCurve2d();
628             aIC.SetFirstCurve2d(aC2D2);
629             aIC.SetSecondCurve2d(aC2D1);
630           }
631         }
632       }
633     }
634     return;
635   }//if(aType1==GeomAbs_Plane && aType2==GeomAbs_Plane){
636
637   if ((aType1==GeomAbs_Plane) && isFace2Quad)
638   {
639     Standard_Real dU, dV;
640
641     // F1
642     Standard_Real umin, umax, vmin, vmax;
643     BRepTools::UVBounds(myFace1, umin, umax, vmin, vmax);
644
645     dU=0.1*(umax-umin);
646     dV=0.1*(vmax-vmin);
647     umin=umin-dU;
648     umax=umax+dU;
649     vmin=vmin-dV;
650     vmax=vmax+dV;
651     myHS1->ChangeSurface().Load(S1, umin, umax, vmin, vmax);
652     // F2
653     BRepTools::UVBounds(myFace2, umin, umax, vmin, vmax);
654     CorrectSurfaceBoundaries(myFace2, (aTolF1 + aTolF2) * 2., umin, umax, vmin, vmax);
655     myHS2->ChangeSurface().Load(S2, umin, umax, vmin, vmax);
656     //
657     if( aType2==GeomAbs_Cone ) { 
658       TolArc = 0.0001; 
659       hasCone = Standard_True; 
660     }
661   }
662   else if ((aType2==GeomAbs_Plane) && isFace1Quad)
663   {
664     Standard_Real dU, dV;
665
666     //F1
667     Standard_Real umin, umax, vmin, vmax;
668     BRepTools::UVBounds(myFace1, umin, umax, vmin, vmax);
669     CorrectSurfaceBoundaries(myFace1, (aTolF1 + aTolF2) * 2., umin, umax, vmin, vmax);
670     myHS1->ChangeSurface().Load(S1, umin, umax, vmin, vmax);
671     // F2
672     BRepTools::UVBounds(myFace2, umin, umax, vmin, vmax);
673     dU=0.1*(umax-umin);
674     dV=0.1*(vmax-vmin);
675     umin=umin-dU;
676     umax=umax+dU;
677     vmin=vmin-dV;
678     vmax=vmax+dV;
679     myHS2->ChangeSurface().Load(S2, umin, umax, vmin, vmax);
680     //
681     if( aType1==GeomAbs_Cone ) {
682       TolArc = 0.0001; 
683       hasCone = Standard_True; 
684     }
685   }
686   else
687   {
688     Standard_Real umin, umax, vmin, vmax;
689     BRepTools::UVBounds(myFace1, umin, umax, vmin, vmax);
690     CorrectSurfaceBoundaries(myFace1, (aTolF1 + aTolF2) * 2., umin, umax, vmin, vmax);
691     myHS1->ChangeSurface().Load(S1, umin, umax, vmin, vmax);
692     BRepTools::UVBounds(myFace2, umin, umax, vmin, vmax);
693     CorrectSurfaceBoundaries(myFace2, (aTolF1 + aTolF2) * 2., umin, umax, vmin, vmax);
694     myHS2->ChangeSurface().Load(S2, umin, umax, vmin, vmax);
695   }
696
697   const Handle(IntTools_TopolTool) dom1 = new IntTools_TopolTool(myHS1);
698   const Handle(IntTools_TopolTool) dom2 = new IntTools_TopolTool(myHS2);
699
700   myLConstruct.Load(dom1, dom2, myHS1, myHS2);
701   
702
703   Tolerances(myHS1, myHS2, TolTang);
704
705   {
706     const Standard_Real UVMaxStep = 0.001;
707     const Standard_Real Deflection = (hasCone) ? 0.085 : 0.1;
708   myIntersector.SetTolerances(TolArc, TolTang, UVMaxStep, Deflection); 
709   }
710   
711   if((myHS1->IsUClosed() && !myHS1->IsUPeriodic()) || 
712      (myHS1->IsVClosed() && !myHS1->IsVPeriodic()) ||
713      (myHS2->IsUClosed() && !myHS2->IsUPeriodic()) || 
714      (myHS2->IsVClosed() && !myHS2->IsVPeriodic()))
715   {
716     RestrictLine = Standard_True;
717   }
718   //
719   if((aType1 != GeomAbs_BSplineSurface) &&
720       (aType1 != GeomAbs_BezierSurface)  &&
721      (aType1 != GeomAbs_OtherSurface)  &&
722      (aType2 != GeomAbs_BSplineSurface) &&
723       (aType2 != GeomAbs_BezierSurface)  &&
724      (aType2 != GeomAbs_OtherSurface))
725   {
726     RestrictLine = Standard_True;
727
728     if ((aType1 == GeomAbs_Torus) ||
729         (aType2 == GeomAbs_Torus))
730     {
731       myListOfPnts.Clear();
732     }
733   }
734
735   //
736   if(!RestrictLine)
737   {
738     TopExp_Explorer aExp;
739     for(Standard_Integer i = 0; (!RestrictLine) && (i < 2); i++)
740     {
741       const TopoDS_Face& aF=(!i) ? myFace1 : myFace2;
742       aExp.Init(aF, TopAbs_EDGE);
743       for(; aExp.More(); aExp.Next())
744       {
745         const TopoDS_Edge& aE=TopoDS::Edge(aExp.Current());
746
747         if(BRep_Tool::Degenerated(aE))
748         {
749           RestrictLine = Standard_True;
750           break;
751         }
752       }
753     }
754   }
755
756   const Standard_Boolean isGeomInt = isTreatAnalityc(aF1, aF2);
757   myIntersector.Perform(myHS1, dom1, myHS2, dom2, TolArc, TolTang, 
758                                   myListOfPnts, RestrictLine, isGeomInt);
759
760   myIsDone = myIntersector.IsDone();
761
762   if (myIsDone)
763   {
764     myTangentFaces=myIntersector.TangentFaces();
765     if (myTangentFaces) {
766       return;
767     }
768     //
769     if(RestrictLine) {
770       myListOfPnts.Clear(); // to use LineConstructor
771     }
772     //
773     const Standard_Integer aNbLin = myIntersector.NbLines();
774     for (Standard_Integer i=1; i <= aNbLin; ++i) {
775       MakeCurve(i, dom1, dom2);
776     }
777     //
778     ComputeTolReached3d();
779     //
780     if (bReverse) {
781       Handle(Geom2d_Curve) aC2D1, aC2D2;
782       //
783       const Standard_Integer aNbLin=mySeqOfCurve.Length();
784       for (Standard_Integer i=1; i<=aNbLin; ++i)
785       {
786         IntTools_Curve& aIC=mySeqOfCurve(i);
787         aC2D1=aIC.FirstCurve2d();
788         aC2D2=aIC.SecondCurve2d();
789         aIC.SetFirstCurve2d(aC2D2);
790         aIC.SetSecondCurve2d(aC2D1);
791       }
792     }
793
794     // Points
795     Standard_Real U1,V1,U2,V2;
796     IntTools_PntOnFace aPntOnF1, aPntOnF2;
797     IntTools_PntOn2Faces aPntOn2Faces;
798     //
799     const Standard_Integer aNbPnts = myIntersector.NbPnts();
800     for (Standard_Integer i=1; i <= aNbPnts; ++i)
801     {
802       const IntSurf_PntOn2S& aISPnt=myIntersector.Point(i).PntOn2S();
803       const gp_Pnt& aPnt=aISPnt.Value();
804       aISPnt.Parameters(U1,V1,U2,V2);
805       aPntOnF1.Init(myFace1, aPnt, U1, V1);
806       aPntOnF2.Init(myFace2, aPnt, U2, V2);
807       //
808       if (!bReverse)
809       {
810         aPntOn2Faces.SetP1(aPntOnF1);
811         aPntOn2Faces.SetP2(aPntOnF2);
812       }
813       else
814       {
815         aPntOn2Faces.SetP2(aPntOnF1);
816         aPntOn2Faces.SetP1(aPntOnF2);
817       }
818
819       myPnts.Append(aPntOn2Faces);
820     }
821   }
822 }
823
824 //=======================================================================
825 //function :ComputeTolReached3d 
826 //purpose  : 
827 //=======================================================================
828   void IntTools_FaceFace::ComputeTolReached3d()
829 {
830   Standard_Boolean bCase1;
831   Standard_Integer aNbLin, i;
832   GeomAbs_SurfaceType aType1, aType2;
833   //
834   aNbLin=myIntersector.NbLines();
835   if (!aNbLin) {
836     return;
837   }
838   //
839   aType1=myHS1->Surface().GetType();
840   aType2=myHS2->Surface().GetType();
841   //
842   bCase1=((aType1==GeomAbs_Plane && aType2==GeomAbs_SurfaceOfExtrusion) ||
843           (aType2==GeomAbs_Plane && aType1==GeomAbs_SurfaceOfExtrusion));
844   //
845   if (aType1==GeomAbs_Cylinder && aType2==GeomAbs_Cylinder) {
846     if (aNbLin==2){ 
847       Handle(IntPatch_Line) aIL1, aIL2;
848       IntPatch_IType aTL1, aTL2;
849       //
850       aIL1=myIntersector.Line(1);
851       aIL2=myIntersector.Line(2);
852       aTL1=aIL1->ArcType();
853       aTL2=aIL2->ArcType();
854       if (aTL1==IntPatch_Lin && aTL2==IntPatch_Lin) {
855         Standard_Real aD, aDTresh, dTol;
856         gp_Lin aL1, aL2;
857         //
858         dTol=1.e-8;
859         aDTresh=1.5e-6;
860         //
861         aL1=Handle(IntPatch_GLine)::DownCast(aIL1)->Line();
862         aL2=Handle(IntPatch_GLine)::DownCast(aIL2)->Line();
863         aD=aL1.Distance(aL2);
864         aD=0.5*aD;
865         if (aD<aDTresh) {
866           myTolReached3d=aD+dTol;
867         }
868         return;
869       }
870     }
871     //ZZ
872     if (aNbLin) {// Check the distances
873       Standard_Integer  aNbP, j ;
874       Standard_Real aT1, aT2, dT, aD2, aD2Max, aEps, aT11, aT12;
875       //
876       aD2Max=0.;
877       aNbP=10;
878       aNbLin=mySeqOfCurve.Length();
879       //
880       for (i=1; i<=aNbLin; ++i) {
881         const IntTools_Curve& aIC=mySeqOfCurve(i);
882         const Handle(Geom_Curve)& aC3D=aIC.Curve();
883         const Handle(Geom2d_Curve)& aC2D1=aIC.FirstCurve2d();
884         const Handle(Geom2d_Curve)& aC2D2=aIC.SecondCurve2d();
885         //
886         if (aC3D.IsNull()) {
887           continue;
888         }
889         const Handle(Geom_BSplineCurve)& aBC=
890           Handle(Geom_BSplineCurve)::DownCast(aC3D);
891         if (aBC.IsNull()) {
892           continue;
893         }
894         //
895         aT1=aBC->FirstParameter();
896         aT2=aBC->LastParameter();
897         //
898         aEps=0.01*(aT2-aT1);
899         dT=(aT2-aT1)/aNbP;
900         for (j=1; j<aNbP; ++j) {
901           aT11=aT1+j*dT;
902           aT12=aT11+dT;
903           aD2=FindMaxSquareDistance(aT11, aT12, aEps, aC3D, aC2D1, aC2D2,
904                                     myHS1, myHS2, myFace1, myFace2, myContext);
905           if (aD2>aD2Max) {
906             aD2Max=aD2;
907           }
908         }
909       }//for (i=1; i<=aNbLin; ++i) {
910       //
911       myTolReached3d=sqrt(aD2Max);
912     }// if (aNbLin) 
913   }// if (aType1==GeomAbs_Cylinder && aType2==GeomAbs_Cylinder) {
914   //
915   //904/G3 f
916   else if (aType1==GeomAbs_Plane && aType2==GeomAbs_Plane) {
917     Standard_Real aTolF1, aTolF2, aTolFMax, aTolTresh;
918     //
919     aTolTresh=1.e-7;
920     //
921     aTolF1 = BRep_Tool::Tolerance(myFace1);
922     aTolF2 = BRep_Tool::Tolerance(myFace2);
923     aTolFMax=Max(aTolF1, aTolF2);
924     //
925     if (aTolFMax>aTolTresh) {
926       myTolReached3d=aTolFMax;
927     }
928   }//if (aType1==GeomAbs_Plane && aType2==GeomAbs_Plane) {
929   //t
930   //IFV Bug OCC20297 
931   else if((aType1 == GeomAbs_Cylinder && aType2 == GeomAbs_Plane) ||
932           (aType2 == GeomAbs_Cylinder && aType1 == GeomAbs_Plane)) {
933     if(aNbLin == 1) {
934       const Handle(IntPatch_Line)& aIL1 = myIntersector.Line(1);
935       if(aIL1->ArcType() == IntPatch_Circle) {
936         gp_Circ aCir = Handle(IntPatch_GLine)::DownCast(aIL1)->Circle();
937         gp_XYZ aCirDir = aCir.Axis().Direction().XYZ();
938         gp_XYZ aPlDir;
939         gp_Pln aPln;
940         if(aType1 == GeomAbs_Plane) {
941           aPln = myHS1->Surface().Plane();
942         }
943         else {
944           aPln = myHS2->Surface().Plane();
945         }
946         aPlDir = aPln.Axis().Direction().XYZ();
947         Standard_Real cs = aCirDir*aPlDir;
948         if(cs < 0.) aPlDir.Reverse();
949         Standard_Real eps = 1.e-14;
950         if(!aPlDir.IsEqual(aCirDir, eps)) {
951           Standard_Integer aNbP = 11;
952           Standard_Real dt = 2.*M_PI / (aNbP - 1), t;
953           for(t = 0.; t < 2.*M_PI; t += dt) {
954             Standard_Real d = aPln.Distance(ElCLib::Value(t, aCir)); 
955             if(myTolReached3d < d) myTolReached3d = d;
956           }
957           myTolReached3d *= 1.1;
958         }
959       } //aIL1->ArcType() == IntPatch_Circle
960     } //aNbLin == 1
961   } // aType1 == GeomAbs_Cylinder && aType2 == GeomAbs_Plane) 
962   //End IFV Bug OCC20297
963   //
964   else if ((aType1==GeomAbs_Plane && aType2==GeomAbs_Torus) ||
965            (aType2==GeomAbs_Plane && aType1==GeomAbs_Torus)) {
966     aNbLin=mySeqOfCurve.Length();
967     if (aNbLin!=1) {
968       return;
969     }
970     //
971     Standard_Integer aNbP;
972     Standard_Real aT, aT1, aT2, dT, aUT, aVT, aUP, aVP;
973     Standard_Real aDP, aDT, aDmax;
974     gp_Pln aPln;
975     gp_Torus aTorus;
976     gp_Pnt aP, aPP, aPT;
977     //
978     const IntTools_Curve& aIC=mySeqOfCurve(1);
979     const Handle(Geom_Curve)& aC3D=aIC.Curve();
980     const Handle(Geom_BSplineCurve)& aBS=
981       Handle(Geom_BSplineCurve)::DownCast(aC3D);
982     if (aBS.IsNull()) {
983       return;
984     }
985     //
986     aT1=aBS->FirstParameter();
987     aT2=aBS->LastParameter();
988     //
989     aPln  =(aType1==GeomAbs_Plane) ? myHS1->Plane() : myHS2->Plane();
990     aTorus=(aType1==GeomAbs_Plane) ? myHS2->Torus() : myHS1->Torus();
991     //
992     aDmax=-1.;
993     aNbP=11;
994     dT=(aT2-aT1)/(aNbP-1);
995     for (i=0; i<aNbP; ++i) {
996       aT=aT1+i*dT;
997       if (i==aNbP-1) {
998         aT=aT2;
999       }
1000       //
1001       aC3D->D0(aT, aP);
1002       //
1003       ElSLib::Parameters(aPln, aP, aUP, aVP);
1004       aPP=ElSLib::Value(aUP, aVP, aPln);
1005       aDP=aP.SquareDistance(aPP);
1006       if (aDP>aDmax) {
1007         aDmax=aDP;
1008       }
1009       //
1010       ElSLib::Parameters(aTorus, aP, aUT, aVT);
1011       aPT=ElSLib::Value(aUT, aVT, aTorus);
1012       aDT=aP.SquareDistance(aPT);
1013       if (aDT>aDmax) {
1014         aDmax=aDT;
1015       }
1016     }
1017     //
1018     if (aDmax > myTolReached3d*myTolReached3d) {
1019       myTolReached3d=sqrt(aDmax);
1020       myTolReached3d=1.1*myTolReached3d;
1021     }
1022   }// if ((aType1==GeomAbs_Plane && aType2==GeomAbs_Torus) ||
1023   //
1024   else if ((aType1==GeomAbs_SurfaceOfRevolution && aType2==GeomAbs_Cylinder) ||
1025            (aType2==GeomAbs_SurfaceOfRevolution && aType1==GeomAbs_Cylinder)) {
1026     Standard_Integer j, aNbP;
1027     Standard_Real aT, aT1, aT2, dT, aD2max, aD2;
1028     //
1029     aNbLin=mySeqOfCurve.Length();
1030     aD2max=0.;
1031     aNbP=11;
1032     //
1033     for (i=1; i<=aNbLin; ++i) {
1034       const IntTools_Curve& aIC=mySeqOfCurve(i);
1035       const Handle(Geom_Curve)& aC3D=aIC.Curve();
1036       const Handle(Geom2d_Curve)& aC2D1=aIC.FirstCurve2d();
1037       const Handle(Geom2d_Curve)& aC2D2=aIC.SecondCurve2d();
1038       //
1039       if (aC3D.IsNull()) {
1040         continue;
1041       }
1042       const Handle(Geom_BSplineCurve)& aBC=
1043         Handle(Geom_BSplineCurve)::DownCast(aC3D);
1044       if (aBC.IsNull()) {
1045         return;
1046       }
1047       //
1048       aT1=aBC->FirstParameter();
1049       aT2=aBC->LastParameter();
1050       //
1051       dT=(aT2-aT1)/(aNbP-1);
1052       for (j=0; j<aNbP; ++j) {
1053         aT=aT1+j*dT;
1054         if (j==aNbP-1) {
1055           aT=aT2;
1056         }
1057         //
1058         aD2=MaxSquareDistance(aT, aC3D, aC2D1, aC2D2,
1059                               myHS1, myHS2, myFace1, myFace2, myContext);
1060         if (aD2>aD2max) {
1061           aD2max=aD2;
1062         }
1063       }//for (j=0; j<aNbP; ++j) {
1064      
1065     }//for (i=1; i<=aNbLin; ++i) {
1066     //
1067     aD2=myTolReached3d*myTolReached3d;
1068     if (aD2max > aD2) {
1069       myTolReached3d=sqrt(aD2max);
1070     }
1071   }//if((aType1==GeomAbs_SurfaceOfRevolution ...
1072   else if ((aType1==GeomAbs_Plane && aType2==GeomAbs_Sphere) ||
1073            (aType2==GeomAbs_Plane && aType1==GeomAbs_Sphere)) {
1074     Standard_Integer  j, aNbP;
1075     Standard_Real aT1, aT2, dT, aD2max, aD2, aEps, aT11, aT12;
1076     //
1077     aNbLin=mySeqOfCurve.Length();
1078     aD2max=0.;
1079     aNbP=10;
1080     //
1081     for (i=1; i<=aNbLin; ++i) {
1082       const IntTools_Curve& aIC=mySeqOfCurve(i);
1083       const Handle(Geom_Curve)& aC3D=aIC.Curve();
1084       const Handle(Geom2d_Curve)& aC2D1=aIC.FirstCurve2d();
1085       const Handle(Geom2d_Curve)& aC2D2=aIC.SecondCurve2d();
1086       //
1087       const Handle(Geom2d_BSplineCurve)& aBC2D1=
1088         Handle(Geom2d_BSplineCurve)::DownCast(aC2D1);
1089       const Handle(Geom2d_BSplineCurve)& aBC2D2=
1090         Handle(Geom2d_BSplineCurve)::DownCast(aC2D2);
1091       //
1092       if (aBC2D1.IsNull() && aBC2D2.IsNull()) {
1093         return;
1094       }
1095       //
1096       if (!aBC2D1.IsNull()) {
1097         aT1=aBC2D1->FirstParameter();
1098         aT2=aBC2D1->LastParameter();
1099       }
1100       else {
1101         aT1=aBC2D2->FirstParameter();
1102         aT2=aBC2D2->LastParameter();
1103       }
1104       //
1105       aEps=0.01*(aT2-aT1);
1106       dT=(aT2-aT1)/aNbP;
1107       for (j=0; j<aNbP; ++j) {
1108         aT11=aT1+j*dT;
1109         aT12=aT11+dT;
1110         if (j==aNbP-1) {
1111           aT12=aT2;
1112         }
1113         //
1114         aD2=FindMaxSquareDistance(aT11, aT12, aEps, aC3D, aC2D1, aC2D2,
1115                                   myHS1, myHS2, myFace1, myFace2, myContext);
1116         if (aD2>aD2max) {
1117           aD2max=aD2;
1118         }
1119       }//for (j=0; j<aNbP; ++j) {
1120      
1121     }//for (i=1; i<=aNbLin; ++i) {
1122     //
1123     aD2=myTolReached3d*myTolReached3d;
1124     if (aD2max > aD2) {
1125       myTolReached3d=sqrt(aD2max);
1126     }
1127   }//else if ((aType1==GeomAbs_Plane && aType2==GeomAbs_Sphere) ...
1128   else if (!myApprox || bCase1) {
1129   //else if (!myApprox) {
1130     Standard_Integer aNbP, j;
1131     Standard_Real aT1, aT2, dT, aD2, aD2Max, aEps, aT11, aT12;
1132     //
1133     aD2Max=0.;
1134     aNbLin=mySeqOfCurve.Length();
1135     //
1136     for (i=1; i<=aNbLin; ++i) {
1137       const IntTools_Curve& aIC=mySeqOfCurve(i);
1138       const Handle(Geom_Curve)& aC3D=aIC.Curve();
1139       const Handle(Geom2d_Curve)& aC2D1=aIC.FirstCurve2d();
1140       const Handle(Geom2d_Curve)& aC2D2=aIC.SecondCurve2d();
1141       //
1142       if (aC3D.IsNull()) {
1143         continue;
1144 }
1145       const Handle(Geom_BSplineCurve)& aBC=
1146         Handle(Geom_BSplineCurve)::DownCast(aC3D);
1147       if (aBC.IsNull()) {
1148         continue;
1149       }
1150       //
1151       aT1=aBC->FirstParameter();
1152       aT2=aBC->LastParameter();
1153       //
1154       aEps=0.0001*(aT2-aT1);
1155       aNbP=11;
1156       dT=(aT2-aT1)/aNbP;
1157       for (j=1; j<aNbP-1; ++j) {
1158         aT11=aT1+j*dT;
1159         aT12=aT11+dT;
1160         aD2=FindMaxSquareDistance(aT11, aT12, aEps, aC3D, aC2D1, aC2D2,
1161                                   myHS1, myHS2, myFace1, myFace2, myContext);
1162         if (aD2>aD2Max) {
1163           aD2Max=aD2;
1164         }
1165       }
1166     }//for (i=1; i<=aNbLin; ++i) {
1167     myTolReached3d=sqrt(aD2Max);
1168   }
1169 }
1170 //=======================================================================
1171 //function : MakeCurve
1172 //purpose  : 
1173 //=======================================================================
1174   void IntTools_FaceFace::MakeCurve(const Standard_Integer Index,
1175                                     const Handle(Adaptor3d_TopolTool)& dom1,
1176                                     const Handle(Adaptor3d_TopolTool)& dom2) 
1177 {
1178   Standard_Boolean bDone, rejectSurface, reApprox, bAvoidLineConstructor;
1179   Standard_Boolean ok, bPCurvesOk;
1180   Standard_Integer i, j, aNbParts;
1181   Standard_Real fprm, lprm;
1182   Standard_Real Tolpc;
1183   Handle(IntPatch_Line) L;
1184   IntPatch_IType typl;
1185   Handle(Geom_Curve) newc;
1186   //
1187   const Standard_Real TOLCHECK   =0.0000001;
1188   const Standard_Real TOLANGCHECK=0.1;
1189   //
1190   rejectSurface = Standard_False;
1191   reApprox = Standard_False;
1192   //
1193   bPCurvesOk = Standard_True;
1194  
1195  reapprox:;
1196   
1197   Tolpc = myTolApprox;
1198   bAvoidLineConstructor = Standard_False;
1199   L = myIntersector.Line(Index);
1200   typl = L->ArcType();
1201   //
1202   if(typl==IntPatch_Walking) {
1203     Handle(IntPatch_Line) anewL;
1204     //
1205     const Handle(IntPatch_WLine)& aWLine=
1206       Handle(IntPatch_WLine)::DownCast(L);
1207     //DumpWLine(aWLine);
1208
1209     anewL = ComputePurgedWLine(aWLine);
1210     if(anewL.IsNull()) {
1211       return;
1212     }
1213     L = anewL;
1214     
1215     //const Handle(IntPatch_WLine)& aWLineX = Handle(IntPatch_WLine)::DownCast(L);
1216     //DumpWLine(aWLineX);
1217
1218     //
1219     if(!myListOfPnts.IsEmpty()) {
1220       bAvoidLineConstructor = Standard_True;
1221     }
1222
1223     Standard_Integer nbp = aWLine->NbPnts();
1224     const IntSurf_PntOn2S& p1 = aWLine->Point(1);
1225     const IntSurf_PntOn2S& p2 = aWLine->Point(nbp);
1226
1227     const gp_Pnt& P1 = p1.Value();
1228     const gp_Pnt& P2 = p2.Value();
1229
1230     if(P1.SquareDistance(P2) < 1.e-14) {
1231       bAvoidLineConstructor = Standard_False;
1232     }
1233   }
1234   //
1235   // Line Constructor
1236   if(!bAvoidLineConstructor) {
1237     myLConstruct.Perform(L);
1238     //
1239     bDone=myLConstruct.IsDone();
1240     aNbParts=myLConstruct.NbParts();
1241     if (!bDone|| !aNbParts) {
1242       return;
1243     }
1244   }
1245   // Do the Curve
1246   
1247   
1248   typl=L->ArcType();
1249   switch (typl) {
1250   //########################################  
1251   // Line, Parabola, Hyperbola
1252   //########################################  
1253   case IntPatch_Lin:
1254   case IntPatch_Parabola: 
1255   case IntPatch_Hyperbola: {
1256     if (typl == IntPatch_Lin) {
1257       newc = 
1258         new Geom_Line (Handle(IntPatch_GLine)::DownCast(L)->Line());
1259     }
1260
1261     else if (typl == IntPatch_Parabola) {
1262       newc = 
1263         new Geom_Parabola(Handle(IntPatch_GLine)::DownCast(L)->Parabola());
1264     }
1265     
1266     else if (typl == IntPatch_Hyperbola) {
1267       newc = 
1268         new Geom_Hyperbola (Handle(IntPatch_GLine)::DownCast(L)->Hyperbola());
1269     }
1270     //
1271     // myTolReached3d
1272     if (typl == IntPatch_Lin) {
1273       TolR3d (myFace1, myFace2, myTolReached3d);
1274     }
1275     //
1276     aNbParts=myLConstruct.NbParts();
1277     for (i=1; i<=aNbParts; i++) {
1278       myLConstruct.Part(i, fprm, lprm);
1279       
1280       if (!Precision::IsNegativeInfinite(fprm) && 
1281           !Precision::IsPositiveInfinite(lprm)) {
1282         //
1283         IntTools_Curve aCurve;
1284         //
1285         Handle(Geom_TrimmedCurve) aCT3D=new Geom_TrimmedCurve(newc, fprm, lprm);
1286         aCurve.SetCurve(aCT3D);
1287         if (typl == IntPatch_Parabola) {
1288           Standard_Real aTolF1, aTolF2, aTolBase;
1289           
1290           aTolF1 = BRep_Tool::Tolerance(myFace1);
1291           aTolF2 = BRep_Tool::Tolerance(myFace2);
1292           aTolBase=aTolF1+aTolF2;
1293           myTolReached3d=IntTools_Tools::CurveTolerance(aCT3D, aTolBase);
1294         }
1295         //
1296         aCurve.SetCurve(new Geom_TrimmedCurve(newc, fprm, lprm));
1297         if(myApprox1) { 
1298           Handle (Geom2d_Curve) C2d;
1299           BuildPCurves(fprm, lprm, Tolpc, myHS1->ChangeSurface().Surface(), newc, C2d);
1300           if(Tolpc>myTolReached2d || myTolReached2d==0.) { 
1301             myTolReached2d=Tolpc;
1302           }
1303             //            
1304             aCurve.SetFirstCurve2d(new Geom2d_TrimmedCurve(C2d,fprm,lprm));
1305           }
1306           else { 
1307             Handle(Geom2d_BSplineCurve) H1;
1308             //
1309             aCurve.SetFirstCurve2d(H1);
1310           }
1311         
1312         if(myApprox2) { 
1313           Handle (Geom2d_Curve) C2d;
1314           BuildPCurves(fprm,lprm,Tolpc,myHS2->ChangeSurface().Surface(),newc,C2d);
1315           if(Tolpc>myTolReached2d || myTolReached2d==0.) { 
1316             myTolReached2d=Tolpc;
1317           }
1318           //
1319           aCurve.SetSecondCurve2d(new Geom2d_TrimmedCurve(C2d,fprm,lprm));
1320           }
1321         else { 
1322           Handle(Geom2d_BSplineCurve) H1;
1323           //
1324           aCurve.SetSecondCurve2d(H1);
1325         }
1326         mySeqOfCurve.Append(aCurve);
1327       } // end of if (!Precision::IsNegativeInfinite(fprm) &&  !Precision::IsPositiveInfinite(lprm))
1328       else {
1329         //  on regarde si on garde
1330         //
1331         Standard_Boolean bFNIt, bLPIt;
1332         Standard_Real aTestPrm, dT=100.;
1333
1334         bFNIt=Precision::IsNegativeInfinite(fprm);
1335         bLPIt=Precision::IsPositiveInfinite(lprm);
1336         
1337         aTestPrm=0.;
1338         
1339         if (bFNIt && !bLPIt) {
1340           aTestPrm=lprm-dT;
1341         }
1342         else if (!bFNIt && bLPIt) {
1343           aTestPrm=fprm+dT;
1344         }
1345         
1346         gp_Pnt ptref(newc->Value(aTestPrm));
1347         //
1348         GeomAbs_SurfaceType typS1 = myHS1->GetType();
1349         GeomAbs_SurfaceType typS2 = myHS2->GetType();
1350         if( typS1 == GeomAbs_SurfaceOfExtrusion ||
1351             typS1 == GeomAbs_OffsetSurface ||
1352             typS1 == GeomAbs_SurfaceOfRevolution ||
1353             typS2 == GeomAbs_SurfaceOfExtrusion ||
1354             typS2 == GeomAbs_OffsetSurface ||
1355             typS2 == GeomAbs_SurfaceOfRevolution) 
1356         {
1357           Handle(Geom2d_BSplineCurve) H1;
1358           mySeqOfCurve.Append(IntTools_Curve(newc, H1, H1));
1359           continue;
1360         }
1361
1362         Standard_Real u1, v1, u2, v2, Tol;
1363         
1364         Tol = Precision::Confusion();
1365         Parameters(myHS1, myHS2, ptref,  u1, v1, u2, v2);
1366         ok = (dom1->Classify(gp_Pnt2d(u1, v1), Tol) != TopAbs_OUT);
1367         if(ok) { 
1368           ok = (dom2->Classify(gp_Pnt2d(u2,v2),Tol) != TopAbs_OUT); 
1369         }
1370         if (ok) {
1371           Handle(Geom2d_BSplineCurve) H1;
1372           mySeqOfCurve.Append(IntTools_Curve(newc, H1, H1));
1373         }
1374       }
1375     }// end of for (i=1; i<=myLConstruct.NbParts(); i++)
1376   }// case IntPatch_Lin:  case IntPatch_Parabola:  case IntPatch_Hyperbola:
1377     break;
1378
1379   //########################################  
1380   // Circle and Ellipse
1381   //########################################  
1382   case IntPatch_Circle: 
1383   case IntPatch_Ellipse: {
1384
1385     if (typl == IntPatch_Circle) {
1386       newc = new Geom_Circle
1387         (Handle(IntPatch_GLine)::DownCast(L)->Circle());
1388     }
1389     else { //IntPatch_Ellipse
1390       newc = new Geom_Ellipse
1391         (Handle(IntPatch_GLine)::DownCast(L)->Ellipse());
1392     }
1393     //
1394     // myTolReached3d
1395     TolR3d (myFace1, myFace2, myTolReached3d);
1396     //
1397     aNbParts=myLConstruct.NbParts();
1398     //
1399     Standard_Real aPeriod, aNul;
1400     TColStd_SequenceOfReal aSeqFprm,  aSeqLprm;
1401     
1402     aNul=0.;
1403     aPeriod=M_PI+M_PI;
1404
1405     for (i=1; i<=aNbParts; i++) {
1406       myLConstruct.Part(i, fprm, lprm);
1407
1408       if (fprm < aNul && lprm > aNul) {
1409         // interval that goes through 0. is divided on two intervals;
1410         while (fprm<aNul || fprm>aPeriod) fprm=fprm+aPeriod;
1411         while (lprm<aNul || lprm>aPeriod) lprm=lprm+aPeriod;
1412         //
1413         if((aPeriod - fprm) > Tolpc) {
1414           aSeqFprm.Append(fprm);
1415           aSeqLprm.Append(aPeriod);
1416         }
1417         else {
1418           gp_Pnt P1 = newc->Value(fprm);
1419           gp_Pnt P2 = newc->Value(aPeriod);
1420           Standard_Real aTolDist = BRep_Tool::Tolerance(myFace1) + BRep_Tool::Tolerance(myFace2);
1421           aTolDist = (myTolReached3d > aTolDist) ? myTolReached3d : aTolDist;
1422
1423           if(P1.Distance(P2) > aTolDist) {
1424             Standard_Real anewpar = fprm;
1425
1426             if(ParameterOutOfBoundary(fprm, newc, myFace1, myFace2, 
1427                                       lprm, Standard_False, anewpar, myContext)) {
1428               fprm = anewpar;
1429             }
1430             aSeqFprm.Append(fprm);
1431             aSeqLprm.Append(aPeriod);
1432           }
1433         }
1434
1435         //
1436         if((lprm - aNul) > Tolpc) {
1437           aSeqFprm.Append(aNul);
1438           aSeqLprm.Append(lprm);
1439         }
1440         else {
1441           gp_Pnt P1 = newc->Value(aNul);
1442           gp_Pnt P2 = newc->Value(lprm);
1443           Standard_Real aTolDist = BRep_Tool::Tolerance(myFace1) + BRep_Tool::Tolerance(myFace2);
1444           aTolDist = (myTolReached3d > aTolDist) ? myTolReached3d : aTolDist;
1445
1446           if(P1.Distance(P2) > aTolDist) {
1447             Standard_Real anewpar = lprm;
1448
1449             if(ParameterOutOfBoundary(lprm, newc, myFace1, myFace2, 
1450                                       fprm, Standard_True, anewpar, myContext)) {
1451               lprm = anewpar;
1452             }
1453             aSeqFprm.Append(aNul);
1454             aSeqLprm.Append(lprm);
1455           }
1456         }
1457       }
1458       else {
1459         // usual interval 
1460         aSeqFprm.Append(fprm);
1461         aSeqLprm.Append(lprm);
1462       }
1463     }
1464     
1465     //
1466     aNbParts=aSeqFprm.Length();
1467     for (i=1; i<=aNbParts; i++) {
1468       fprm=aSeqFprm(i);
1469       lprm=aSeqLprm(i);
1470       //
1471       Standard_Real aRealEpsilon=RealEpsilon();
1472       if (Abs(fprm) > aRealEpsilon || Abs(lprm-2.*M_PI) > aRealEpsilon) {
1473         //==============================================
1474         ////
1475         IntTools_Curve aCurve;
1476         Handle(Geom_TrimmedCurve) aTC3D=new Geom_TrimmedCurve(newc,fprm,lprm);
1477         aCurve.SetCurve(aTC3D);
1478         fprm=aTC3D->FirstParameter();
1479         lprm=aTC3D->LastParameter ();
1480         ////         
1481         if (typl == IntPatch_Circle || typl == IntPatch_Ellipse) {//// 
1482           if(myApprox1) { 
1483             Handle (Geom2d_Curve) C2d;
1484             BuildPCurves(fprm,lprm,Tolpc,myHS1->ChangeSurface().Surface(),newc,C2d);
1485             if(Tolpc>myTolReached2d || myTolReached2d==0) { 
1486               myTolReached2d=Tolpc;
1487             }
1488             //
1489             aCurve.SetFirstCurve2d(C2d);
1490           }
1491           else { //// 
1492             Handle(Geom2d_BSplineCurve) H1;
1493             aCurve.SetFirstCurve2d(H1);
1494           }
1495
1496
1497           if(myApprox2) { 
1498             Handle (Geom2d_Curve) C2d;
1499             BuildPCurves(fprm,lprm,Tolpc,myHS2->ChangeSurface().Surface(),newc,C2d);
1500             if(Tolpc>myTolReached2d || myTolReached2d==0) { 
1501               myTolReached2d=Tolpc;
1502             }
1503             //
1504             aCurve.SetSecondCurve2d(C2d);
1505           }
1506           else { 
1507             Handle(Geom2d_BSplineCurve) H1;
1508             aCurve.SetSecondCurve2d(H1);
1509           }
1510         }
1511         
1512         else { 
1513           Handle(Geom2d_BSplineCurve) H1;
1514           aCurve.SetFirstCurve2d(H1);
1515           aCurve.SetSecondCurve2d(H1);
1516         }
1517         mySeqOfCurve.Append(aCurve);
1518           //==============================================        
1519       } //if (Abs(fprm) > RealEpsilon() || Abs(lprm-2.*M_PI) > RealEpsilon())
1520
1521       else {
1522         //  on regarde si on garde
1523         //
1524         if (aNbParts==1) {
1525 //           if (Abs(fprm) < RealEpsilon() &&  Abs(lprm-2.*M_PI) < RealEpsilon()) {
1526           if (Abs(fprm) <= aRealEpsilon && Abs(lprm-2.*M_PI) <= aRealEpsilon) {
1527             IntTools_Curve aCurve;
1528             Handle(Geom_TrimmedCurve) aTC3D=new Geom_TrimmedCurve(newc,fprm,lprm);
1529             aCurve.SetCurve(aTC3D);
1530             fprm=aTC3D->FirstParameter();
1531             lprm=aTC3D->LastParameter ();
1532             
1533             if(myApprox1) { 
1534               Handle (Geom2d_Curve) C2d;
1535               BuildPCurves(fprm,lprm,Tolpc,myHS1->ChangeSurface().Surface(),newc,C2d);
1536               if(Tolpc>myTolReached2d || myTolReached2d==0) { 
1537                 myTolReached2d=Tolpc;
1538               }
1539               //
1540               aCurve.SetFirstCurve2d(C2d);
1541             }
1542             else { //// 
1543               Handle(Geom2d_BSplineCurve) H1;
1544               aCurve.SetFirstCurve2d(H1);
1545             }
1546
1547             if(myApprox2) { 
1548               Handle (Geom2d_Curve) C2d;
1549               BuildPCurves(fprm,lprm,Tolpc,myHS2->ChangeSurface().Surface(),newc,C2d);
1550               if(Tolpc>myTolReached2d || myTolReached2d==0) { 
1551                 myTolReached2d=Tolpc;
1552               }
1553               //
1554               aCurve.SetSecondCurve2d(C2d);
1555             }
1556             else { 
1557               Handle(Geom2d_BSplineCurve) H1;
1558               aCurve.SetSecondCurve2d(H1);
1559             }
1560             mySeqOfCurve.Append(aCurve);
1561             break;
1562           }
1563         }
1564         //
1565         Standard_Real aTwoPIdiv17, u1, v1, u2, v2, Tol;
1566
1567         aTwoPIdiv17=2.*M_PI/17.;
1568
1569         for (j=0; j<=17; j++) {
1570           gp_Pnt ptref (newc->Value (j*aTwoPIdiv17));
1571           Tol = Precision::Confusion();
1572
1573           Parameters(myHS1, myHS2, ptref, u1, v1, u2, v2);
1574           ok = (dom1->Classify(gp_Pnt2d(u1,v1),Tol) != TopAbs_OUT);
1575           if(ok) { 
1576             ok = (dom2->Classify(gp_Pnt2d(u2,v2),Tol) != TopAbs_OUT);
1577           }
1578           if (ok) {
1579             IntTools_Curve aCurve;
1580             aCurve.SetCurve(newc);
1581             //==============================================
1582             if (typl == IntPatch_Circle || typl == IntPatch_Ellipse) {
1583               
1584               if(myApprox1) { 
1585                 Handle (Geom2d_Curve) C2d;
1586                 BuildPCurves(fprm, lprm, Tolpc, myHS1->ChangeSurface().Surface(), newc, C2d);
1587                 if(Tolpc>myTolReached2d || myTolReached2d==0) { 
1588                   myTolReached2d=Tolpc;
1589                 }
1590                 // 
1591                 aCurve.SetFirstCurve2d(C2d);
1592               }
1593               else { 
1594                 Handle(Geom2d_BSplineCurve) H1;
1595                 aCurve.SetFirstCurve2d(H1);
1596               }
1597                 
1598               if(myApprox2) { 
1599                 Handle (Geom2d_Curve) C2d;
1600                 BuildPCurves(fprm, lprm, Tolpc,myHS2->ChangeSurface().Surface(), newc, C2d);
1601                 if(Tolpc>myTolReached2d || myTolReached2d==0) { 
1602                   myTolReached2d=Tolpc;
1603                 }
1604                 //                 
1605                 aCurve.SetSecondCurve2d(C2d);
1606               }
1607                 
1608               else { 
1609                 Handle(Geom2d_BSplineCurve) H1;
1610                 aCurve.SetSecondCurve2d(H1);
1611               }
1612             }//  end of if (typl == IntPatch_Circle || typl == IntPatch_Ellipse)
1613              
1614             else { 
1615               Handle(Geom2d_BSplineCurve) H1;
1616               //         
1617               aCurve.SetFirstCurve2d(H1);
1618               aCurve.SetSecondCurve2d(H1);
1619             }
1620             //==============================================        
1621             //
1622             mySeqOfCurve.Append(aCurve);
1623             break;
1624
1625             }//  end of if (ok) {
1626           }//  end of for (Standard_Integer j=0; j<=17; j++)
1627         }//  end of else { on regarde si on garde
1628       }// for (i=1; i<=myLConstruct.NbParts(); i++)
1629     }// IntPatch_Circle: IntPatch_Ellipse:
1630     break;
1631     
1632   case IntPatch_Analytic: {
1633     IntSurf_Quadric quad1,quad2;
1634     GeomAbs_SurfaceType typs = myHS1->Surface().GetType();
1635     
1636     switch (typs) {
1637       case GeomAbs_Plane:
1638         quad1.SetValue(myHS1->Surface().Plane());
1639         break;
1640       case GeomAbs_Cylinder:
1641         quad1.SetValue(myHS1->Surface().Cylinder());
1642         break;
1643       case GeomAbs_Cone:
1644         quad1.SetValue(myHS1->Surface().Cone());
1645         break;
1646       case GeomAbs_Sphere:
1647         quad1.SetValue(myHS1->Surface().Sphere());
1648         break;
1649     case GeomAbs_Torus:
1650       quad1.SetValue(myHS1->Surface().Torus());
1651       break;
1652       default:
1653         Standard_ConstructionError::Raise("GeomInt_IntSS::MakeCurve 1");
1654       }
1655       
1656     typs = myHS2->Surface().GetType();
1657     
1658     switch (typs) {
1659       case GeomAbs_Plane:
1660         quad2.SetValue(myHS2->Surface().Plane());
1661         break;
1662       case GeomAbs_Cylinder:
1663         quad2.SetValue(myHS2->Surface().Cylinder());
1664         break;
1665       case GeomAbs_Cone:
1666         quad2.SetValue(myHS2->Surface().Cone());
1667         break;
1668       case GeomAbs_Sphere:
1669         quad2.SetValue(myHS2->Surface().Sphere());
1670         break;
1671     case GeomAbs_Torus:
1672       quad2.SetValue(myHS2->Surface().Torus());
1673       break;
1674       default:
1675         Standard_ConstructionError::Raise("GeomInt_IntSS::MakeCurve 2");
1676       }
1677     //
1678     //=========
1679     IntPatch_ALineToWLine convert (quad1, quad2);
1680       
1681     if (!myApprox) {
1682       aNbParts=myLConstruct.NbParts();
1683       for (i=1; i<=aNbParts; i++) {
1684         myLConstruct.Part(i, fprm, lprm);
1685         Handle(IntPatch_WLine) WL = 
1686           convert.MakeWLine(Handle(IntPatch_ALine)::DownCast(L), fprm, lprm);
1687         //
1688         Handle(Geom2d_BSplineCurve) H1;
1689         Handle(Geom2d_BSplineCurve) H2;
1690
1691         if(myApprox1) {
1692           H1 = MakeBSpline2d(WL, 1, WL->NbPnts(), Standard_True);
1693         }
1694         
1695         if(myApprox2) {
1696           H2 = MakeBSpline2d(WL, 1, WL->NbPnts(), Standard_False);
1697         }
1698         //          
1699         mySeqOfCurve.Append(IntTools_Curve(MakeBSpline(WL,1,WL->NbPnts()), H1, H2));
1700       }
1701     } // if (!myApprox)
1702
1703     else { // myApprox=TRUE
1704       GeomInt_WLApprox theapp3d;
1705       // 
1706       Standard_Real tol2d = myTolApprox;
1707       //         
1708       theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, Standard_True);
1709       
1710       aNbParts=myLConstruct.NbParts();
1711       for (i=1; i<=aNbParts; i++) {
1712         myLConstruct.Part(i, fprm, lprm);
1713         Handle(IntPatch_WLine) WL = 
1714           convert.MakeWLine(Handle(IntPatch_ALine):: DownCast(L),fprm,lprm);
1715
1716         theapp3d.Perform(myHS1,myHS2,WL,Standard_True,myApprox1,myApprox2, 1, WL->NbPnts());
1717
1718         if (!theapp3d.IsDone()) {
1719           //
1720           Handle(Geom2d_BSplineCurve) H1;
1721           Handle(Geom2d_BSplineCurve) H2;
1722
1723           if(myApprox1) {
1724             H1 = MakeBSpline2d(WL, 1, WL->NbPnts(), Standard_True);
1725           }
1726           
1727           if(myApprox2) {
1728             H2 = MakeBSpline2d(WL, 1, WL->NbPnts(), Standard_False);
1729           }
1730           //          
1731           mySeqOfCurve.Append(IntTools_Curve(MakeBSpline(WL,1,WL->NbPnts()), H1, H2));
1732         }
1733
1734         else {
1735           if(myApprox1 || myApprox2) { 
1736             if( theapp3d.TolReached2d()>myTolReached2d || myTolReached2d==0) { 
1737               myTolReached2d = theapp3d.TolReached2d();
1738             }
1739           }
1740           
1741           if( theapp3d.TolReached3d()>myTolReached3d || myTolReached3d==0) { 
1742             myTolReached3d = theapp3d.TolReached3d();
1743           }
1744
1745           Standard_Integer aNbMultiCurves, nbpoles;
1746           aNbMultiCurves=theapp3d.NbMultiCurves();
1747           for (j=1; j<=aNbMultiCurves; j++) {
1748             const AppParCurves_MultiBSpCurve& mbspc = theapp3d.Value(j);
1749             nbpoles = mbspc.NbPoles();
1750             
1751             TColgp_Array1OfPnt tpoles(1, nbpoles);
1752             mbspc.Curve(1, tpoles);
1753             Handle(Geom_BSplineCurve) BS=new Geom_BSplineCurve(tpoles,
1754                                                                mbspc.Knots(),
1755                                                                mbspc.Multiplicities(),
1756                                                                mbspc.Degree());
1757             
1758             GeomLib_CheckBSplineCurve Check(BS,TOLCHECK,TOLANGCHECK);
1759             Check.FixTangent(Standard_True,Standard_True);
1760             // 
1761             IntTools_Curve aCurve;
1762             aCurve.SetCurve(BS);
1763             
1764             if(myApprox1) { 
1765               TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
1766               mbspc.Curve(2,tpoles2d);
1767               Handle(Geom2d_BSplineCurve) BS2=new Geom2d_BSplineCurve(tpoles2d,
1768                                                                       mbspc.Knots(),
1769                                                                       mbspc.Multiplicities(),
1770                                                                       mbspc.Degree());
1771
1772               GeomLib_Check2dBSplineCurve newCheck(BS2,TOLCHECK,TOLANGCHECK);
1773               newCheck.FixTangent(Standard_True,Standard_True);
1774               //                 
1775               aCurve.SetFirstCurve2d(BS2);
1776             }
1777             else {
1778               Handle(Geom2d_BSplineCurve) H1;
1779               aCurve.SetFirstCurve2d(H1);
1780             }
1781             
1782             if(myApprox2) { 
1783               TColgp_Array1OfPnt2d tpoles2d(1, nbpoles);
1784               Standard_Integer TwoOrThree;
1785               TwoOrThree=myApprox1 ? 3 : 2;
1786               mbspc.Curve(TwoOrThree, tpoles2d);
1787               Handle(Geom2d_BSplineCurve) BS2 =new Geom2d_BSplineCurve(tpoles2d,
1788                                                                        mbspc.Knots(),
1789                                                                        mbspc.Multiplicities(),
1790                                                                        mbspc.Degree());
1791                 
1792               GeomLib_Check2dBSplineCurve newCheck(BS2,TOLCHECK,TOLANGCHECK);
1793               newCheck.FixTangent(Standard_True,Standard_True);
1794               //         
1795               aCurve.SetSecondCurve2d(BS2);
1796             }
1797             else { 
1798               Handle(Geom2d_BSplineCurve) H2;
1799               aCurve.SetSecondCurve2d(H2);
1800             }
1801             // 
1802             mySeqOfCurve.Append(aCurve);
1803
1804           }// for (j=1; j<=aNbMultiCurves; j++) {
1805         }// else from if (!theapp3d.IsDone())
1806       }// for (i=1; i<=aNbParts; i++) {
1807     }// else { // myApprox=TRUE
1808   }// case IntPatch_Analytic:
1809     break;
1810
1811   case IntPatch_Walking:{
1812     Handle(IntPatch_WLine) WL = 
1813       Handle(IntPatch_WLine)::DownCast(L);
1814     //
1815     Standard_Integer ifprm, ilprm;
1816     //
1817     if (!myApprox) {
1818       aNbParts = 1;
1819       if(!bAvoidLineConstructor){
1820         aNbParts=myLConstruct.NbParts();
1821       }
1822       for (i=1; i<=aNbParts; ++i) {
1823         Handle(Geom2d_BSplineCurve) H1, H2;
1824         Handle(Geom_Curve) aBSp;
1825         //
1826         if(bAvoidLineConstructor) {
1827           ifprm = 1;
1828           ilprm = WL->NbPnts();
1829         }
1830         else {
1831           myLConstruct.Part(i, fprm, lprm);
1832           ifprm=(Standard_Integer)fprm;
1833           ilprm=(Standard_Integer)lprm;
1834         }
1835         //
1836         if(myApprox1) {
1837           H1 = MakeBSpline2d(WL, ifprm, ilprm, Standard_True);
1838         }
1839         //
1840         if(myApprox2) {
1841           H2 = MakeBSpline2d(WL, ifprm, ilprm, Standard_False);
1842         }
1843         //           
1844         aBSp=MakeBSpline(WL, ifprm, ilprm);
1845         IntTools_Curve aIC(aBSp, H1, H2);
1846         mySeqOfCurve.Append(aIC);
1847       }// for (i=1; i<=aNbParts; ++i) {
1848     }// if (!myApprox) {
1849     //
1850     else { // X
1851       Standard_Boolean bIsDecomposited;
1852       Standard_Integer nbiter, aNbSeqOfL;
1853       Standard_Real tol2d, aTolApproxImp;
1854       IntPatch_SequenceOfLine aSeqOfL;
1855       GeomInt_WLApprox theapp3d;
1856       Approx_ParametrizationType aParType = Approx_ChordLength;
1857       //
1858       Standard_Boolean anApprox1 = myApprox1;
1859       Standard_Boolean anApprox2 = myApprox2;
1860       //
1861       aTolApproxImp=1.e-5;
1862       tol2d = myTolApprox;
1863
1864       GeomAbs_SurfaceType typs1, typs2;
1865       typs1 = myHS1->Surface().GetType();
1866       typs2 = myHS2->Surface().GetType();
1867       Standard_Boolean anWithPC = Standard_True;
1868
1869       if(typs1 == GeomAbs_Cylinder && typs2 == GeomAbs_Sphere) {
1870         anWithPC = 
1871           ApproxWithPCurves(myHS1->Surface().Cylinder(), myHS2->Surface().Sphere());
1872       }
1873       else if (typs1 == GeomAbs_Sphere && typs2 == GeomAbs_Cylinder) {
1874         anWithPC = 
1875           ApproxWithPCurves(myHS2->Surface().Cylinder(), myHS1->Surface().Sphere());
1876       }
1877       //
1878       if(!anWithPC) {
1879         myTolApprox = aTolApproxImp;//1.e-5; 
1880         anApprox1 = Standard_False;
1881         anApprox2 = Standard_False;
1882         //         
1883         tol2d = myTolApprox;
1884       }
1885         
1886       if(myHS1 == myHS2) { 
1887         theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, Standard_False, aParType);
1888         rejectSurface = Standard_True;
1889       }
1890       else { 
1891         if(reApprox && !rejectSurface)
1892           theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, Standard_False, aParType);
1893         else {
1894           Standard_Integer iDegMax, iDegMin, iNbIter;
1895           //
1896           ApproxParameters(myHS1, myHS2, iDegMin, iDegMax, iNbIter);
1897           theapp3d.SetParameters(myTolApprox, tol2d, iDegMin, iDegMax, iNbIter, Standard_True, aParType);
1898         }
1899       }
1900       //
1901       Standard_Real aReachedTol = Precision::Confusion();
1902       bIsDecomposited=DecompositionOfWLine(WL,
1903                                            myHS1, 
1904                                            myHS2, 
1905                                            myFace1, 
1906                                            myFace2, 
1907                                            myLConstruct, 
1908                                            bAvoidLineConstructor, 
1909                                            aSeqOfL, 
1910                                            aReachedTol,
1911                                            myContext);
1912       if ( bIsDecomposited && ( myTolReached3d < aReachedTol ) ) {
1913         myTolReached3d = aReachedTol;
1914       }
1915       //
1916       aNbSeqOfL=aSeqOfL.Length();
1917       //
1918       if (bIsDecomposited) {
1919         nbiter=aNbSeqOfL;
1920       }
1921       else {
1922         nbiter=1;
1923         aNbParts=1;
1924         if (!bAvoidLineConstructor) {
1925           aNbParts=myLConstruct.NbParts();
1926           nbiter=aNbParts;
1927         }
1928       }
1929       //
1930       for(i = 1; i <= nbiter; ++i) {
1931         if(bIsDecomposited) {
1932           WL = Handle(IntPatch_WLine)::DownCast(aSeqOfL.Value(i));
1933           ifprm = 1;
1934           ilprm = WL->NbPnts();
1935         }
1936         else {
1937           if(bAvoidLineConstructor) {
1938             ifprm = 1;
1939             ilprm = WL->NbPnts();
1940           }
1941           else {
1942             myLConstruct.Part(i, fprm, lprm);
1943             ifprm = (Standard_Integer)fprm;
1944             ilprm = (Standard_Integer)lprm;
1945           }
1946         }
1947         //-- lbr : 
1948         //-- Si une des surfaces est un plan , on approxime en 2d
1949         //-- sur cette surface et on remonte les points 2d en 3d.
1950         if(typs1 == GeomAbs_Plane) { 
1951           theapp3d.Perform(myHS1, myHS2, WL, Standard_False,Standard_True, myApprox2,ifprm,ilprm);
1952         }          
1953         else if(typs2 == GeomAbs_Plane) { 
1954           theapp3d.Perform(myHS1,myHS2,WL,Standard_False,myApprox1,Standard_True,ifprm,ilprm);
1955         }
1956         else { 
1957           //
1958           if (myHS1 != myHS2){
1959             if ((typs1==GeomAbs_BezierSurface || typs1==GeomAbs_BSplineSurface) &&
1960                 (typs2==GeomAbs_BezierSurface || typs2==GeomAbs_BSplineSurface)) {
1961              
1962               theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, Standard_True, aParType);
1963               
1964               Standard_Boolean bUseSurfaces;
1965               bUseSurfaces=NotUseSurfacesForApprox(myFace1, myFace2, WL, ifprm,  ilprm);
1966               if (bUseSurfaces) {
1967                 // ######
1968                 rejectSurface = Standard_True;
1969                 // ######
1970                 theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, Standard_False, aParType);
1971               }
1972             }
1973           }
1974           //
1975           theapp3d.Perform(myHS1,myHS2,WL,Standard_True,anApprox1,anApprox2,ifprm,ilprm);
1976         }
1977           //           
1978         if (!theapp3d.IsDone()) {
1979           Handle(Geom2d_BSplineCurve) H1;
1980           Handle(Geom2d_BSplineCurve) H2;
1981           //           
1982           Handle(Geom_Curve) aBSp=MakeBSpline(WL,ifprm, ilprm);
1983           //
1984           if(myApprox1) {
1985             H1 = MakeBSpline2d(WL, ifprm, ilprm, Standard_True);
1986           }
1987           //
1988           if(myApprox2) {
1989             H2 = MakeBSpline2d(WL, ifprm, ilprm, Standard_False);
1990           }
1991           //           
1992           IntTools_Curve aIC(aBSp, H1, H2);
1993           mySeqOfCurve.Append(aIC);
1994         }
1995         
1996         else {
1997           if(myApprox1 || myApprox2 || (typs1==GeomAbs_Plane || typs2==GeomAbs_Plane)) { 
1998             if( theapp3d.TolReached2d()>myTolReached2d || myTolReached2d==0.) { 
1999               myTolReached2d = theapp3d.TolReached2d();
2000             }
2001           }
2002           if(typs1==GeomAbs_Plane || typs2==GeomAbs_Plane) { 
2003             myTolReached3d = myTolReached2d;
2004             //
2005             if (typs1==GeomAbs_Torus || typs2==GeomAbs_Torus) {
2006               if (myTolReached3d<1.e-6) {
2007                 myTolReached3d = theapp3d.TolReached3d();
2008                 myTolReached3d=1.e-6;
2009               }
2010             }
2011           }
2012           else  if( theapp3d.TolReached3d()>myTolReached3d || myTolReached3d==0.) { 
2013             myTolReached3d = theapp3d.TolReached3d();
2014           }
2015           
2016           Standard_Integer aNbMultiCurves, nbpoles;
2017           aNbMultiCurves=theapp3d.NbMultiCurves(); 
2018           for (j=1; j<=aNbMultiCurves; j++) {
2019             if(typs1 == GeomAbs_Plane) {
2020               const AppParCurves_MultiBSpCurve& mbspc = theapp3d.Value(j);
2021               nbpoles = mbspc.NbPoles();
2022               
2023               TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
2024               TColgp_Array1OfPnt   tpoles(1,nbpoles);
2025               
2026               mbspc.Curve(1,tpoles2d);
2027               const gp_Pln&  Pln = myHS1->Surface().Plane();
2028               //
2029               Standard_Integer ik; 
2030               for(ik = 1; ik<= nbpoles; ik++) { 
2031                 tpoles.SetValue(ik,
2032                                 ElSLib::Value(tpoles2d.Value(ik).X(),
2033                                               tpoles2d.Value(ik).Y(),
2034                                               Pln));
2035               }
2036               //
2037               Handle(Geom_BSplineCurve) BS = 
2038                 new Geom_BSplineCurve(tpoles,
2039                                       mbspc.Knots(),
2040                                       mbspc.Multiplicities(),
2041                                       mbspc.Degree());
2042               GeomLib_CheckBSplineCurve Check(BS,TOLCHECK,TOLANGCHECK);
2043               Check.FixTangent(Standard_True, Standard_True);
2044               //         
2045               IntTools_Curve aCurve;
2046               aCurve.SetCurve(BS);
2047
2048               if(myApprox1) { 
2049                 Handle(Geom2d_BSplineCurve) BS1 = 
2050                   new Geom2d_BSplineCurve(tpoles2d,
2051                                           mbspc.Knots(),
2052                                           mbspc.Multiplicities(),
2053                                           mbspc.Degree());
2054                 GeomLib_Check2dBSplineCurve Check1(BS1,TOLCHECK,TOLANGCHECK);
2055                 Check1.FixTangent(Standard_True,Standard_True);
2056                 //
2057                 // ############################################
2058                 if(!rejectSurface && !reApprox) {
2059                   Standard_Boolean isValid = IsCurveValid(BS1);
2060                   if(!isValid) {
2061                     reApprox = Standard_True;
2062                     goto reapprox;
2063                   }
2064                 }
2065                 // ############################################
2066                 aCurve.SetFirstCurve2d(BS1);
2067               }
2068               else {
2069                 Handle(Geom2d_BSplineCurve) H1;
2070                 aCurve.SetFirstCurve2d(H1);
2071               }
2072
2073               if(myApprox2) { 
2074                 mbspc.Curve(2, tpoles2d);
2075                 
2076                 Handle(Geom2d_BSplineCurve) BS2 = new Geom2d_BSplineCurve(tpoles2d,
2077                                                                           mbspc.Knots(),
2078                                                                           mbspc.Multiplicities(),
2079                                                                           mbspc.Degree());
2080                 GeomLib_Check2dBSplineCurve newCheck(BS2,TOLCHECK,TOLANGCHECK);
2081                 newCheck.FixTangent(Standard_True,Standard_True);
2082                 
2083                 // ###########################################
2084                 if(!rejectSurface && !reApprox) {
2085                   Standard_Boolean isValid = IsCurveValid(BS2);
2086                   if(!isValid) {
2087                     reApprox = Standard_True;
2088                     goto reapprox;
2089                   }
2090                 }
2091                 // ###########################################
2092                 // 
2093                 aCurve.SetSecondCurve2d(BS2);
2094               }
2095               else { 
2096                 Handle(Geom2d_BSplineCurve) H2;
2097                 //                 
2098                   aCurve.SetSecondCurve2d(H2);
2099               }
2100               //
2101               mySeqOfCurve.Append(aCurve);
2102             }//if(typs1 == GeomAbs_Plane) {
2103             
2104             else if(typs2 == GeomAbs_Plane) { 
2105               const AppParCurves_MultiBSpCurve& mbspc = theapp3d.Value(j);
2106               nbpoles = mbspc.NbPoles();
2107               
2108               TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
2109               TColgp_Array1OfPnt   tpoles(1,nbpoles);
2110               mbspc.Curve((myApprox1==Standard_True)? 2 : 1,tpoles2d);
2111               const gp_Pln&  Pln = myHS2->Surface().Plane();
2112               //
2113               Standard_Integer ik; 
2114               for(ik = 1; ik<= nbpoles; ik++) { 
2115                 tpoles.SetValue(ik,
2116                                 ElSLib::Value(tpoles2d.Value(ik).X(),
2117                                               tpoles2d.Value(ik).Y(),
2118                                               Pln));
2119                 
2120               }
2121               //
2122               Handle(Geom_BSplineCurve) BS=new Geom_BSplineCurve(tpoles,
2123                                                                  mbspc.Knots(),
2124                                                                  mbspc.Multiplicities(),
2125                                                                  mbspc.Degree());
2126               GeomLib_CheckBSplineCurve Check(BS,TOLCHECK,TOLANGCHECK);
2127               Check.FixTangent(Standard_True,Standard_True);
2128               //         
2129               IntTools_Curve aCurve;
2130               aCurve.SetCurve(BS);
2131
2132               if(myApprox2) {
2133                 Handle(Geom2d_BSplineCurve) BS1=new Geom2d_BSplineCurve(tpoles2d,
2134                                                                         mbspc.Knots(),
2135                                                                         mbspc.Multiplicities(),
2136                                                                         mbspc.Degree());
2137                 GeomLib_Check2dBSplineCurve Check1(BS1,TOLCHECK,TOLANGCHECK);
2138                 Check1.FixTangent(Standard_True,Standard_True);
2139                 //         
2140                 // ###########################################
2141                 if(!rejectSurface && !reApprox) {
2142                   Standard_Boolean isValid = IsCurveValid(BS1);
2143                   if(!isValid) {
2144                     reApprox = Standard_True;
2145                     goto reapprox;
2146                   }
2147                 }
2148                 // ###########################################
2149                 bPCurvesOk = CheckPCurve(BS1, myFace2);
2150                 aCurve.SetSecondCurve2d(BS1);
2151               }
2152               else {
2153                 Handle(Geom2d_BSplineCurve) H2;
2154                 aCurve.SetSecondCurve2d(H2);
2155               }
2156               
2157               if(myApprox1) { 
2158                 mbspc.Curve(1,tpoles2d);
2159                 Handle(Geom2d_BSplineCurve) BS2=new Geom2d_BSplineCurve(tpoles2d,
2160                                                                         mbspc.Knots(),
2161                                                                         mbspc.Multiplicities(),
2162                                                                         mbspc.Degree());
2163                 GeomLib_Check2dBSplineCurve Check2(BS2,TOLCHECK,TOLANGCHECK);
2164                 Check2.FixTangent(Standard_True,Standard_True);
2165                 //
2166                 // ###########################################
2167                 if(!rejectSurface && !reApprox) {
2168                   Standard_Boolean isValid = IsCurveValid(BS2);
2169                   if(!isValid) {
2170                     reApprox = Standard_True;
2171                     goto reapprox;
2172                   }
2173                 }
2174                 // ###########################################
2175                 bPCurvesOk = bPCurvesOk && CheckPCurve(BS2, myFace1);
2176                 aCurve.SetFirstCurve2d(BS2);
2177               }
2178               else { 
2179                 Handle(Geom2d_BSplineCurve) H1;
2180                 //                 
2181                 aCurve.SetFirstCurve2d(H1);
2182               }
2183               //
2184               //if points of the pcurves are out of the faces bounds
2185               //create 3d and 2d curves without approximation
2186               if (!bPCurvesOk) {
2187                 Handle(Geom2d_BSplineCurve) H1, H2;
2188                 bPCurvesOk = Standard_True;
2189                 //           
2190                 Handle(Geom_Curve) aBSp=MakeBSpline(WL,ifprm, ilprm);
2191                 
2192                 if(myApprox1) {
2193                   H1 = MakeBSpline2d(WL, ifprm, ilprm, Standard_True);
2194                   bPCurvesOk = CheckPCurve(H1, myFace1);
2195                 }
2196                 
2197                 if(myApprox2) {
2198                   H2 = MakeBSpline2d(WL, ifprm, ilprm, Standard_False);
2199                   bPCurvesOk = bPCurvesOk && CheckPCurve(H2, myFace2);
2200                 }
2201                 //
2202                 //if pcurves created without approximation are out of the 
2203                 //faces bounds, use approximated 3d and 2d curves
2204                 if (bPCurvesOk) {
2205                   IntTools_Curve aIC(aBSp, H1, H2);
2206                   mySeqOfCurve.Append(aIC);
2207                 } else {
2208                   mySeqOfCurve.Append(aCurve);
2209                 }
2210               } else {
2211                 mySeqOfCurve.Append(aCurve);
2212               }
2213             }// else if(typs2 == GeomAbs_Plane)
2214             //
2215             else { //typs2 != GeomAbs_Plane && typs1 != GeomAbs_Plane
2216               Standard_Boolean bIsValid1, bIsValid2;
2217               Handle(Geom_BSplineCurve) BS;
2218               Handle(Geom2d_BSplineCurve) aH2D;        
2219               IntTools_Curve aCurve; 
2220               //
2221               bIsValid1=Standard_True;
2222               bIsValid2=Standard_True;
2223               //
2224               const AppParCurves_MultiBSpCurve& mbspc = theapp3d.Value(j);
2225               nbpoles = mbspc.NbPoles();
2226               TColgp_Array1OfPnt tpoles(1,nbpoles);
2227               mbspc.Curve(1,tpoles);
2228               BS=new Geom_BSplineCurve(tpoles,
2229                                                                  mbspc.Knots(),
2230                                                                  mbspc.Multiplicities(),
2231                                                                  mbspc.Degree());
2232               GeomLib_CheckBSplineCurve Check(BS,TOLCHECK,TOLANGCHECK);
2233               Check.FixTangent(Standard_True,Standard_True);
2234               //                 
2235               aCurve.SetCurve(BS);
2236               aCurve.SetFirstCurve2d(aH2D);
2237               aCurve.SetSecondCurve2d(aH2D);
2238               //
2239               if(myApprox1) { 
2240                 if(anApprox1) {
2241                   Handle(Geom2d_BSplineCurve) BS1;
2242                   TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
2243                   mbspc.Curve(2,tpoles2d);
2244                   //
2245                   BS1=new Geom2d_BSplineCurve(tpoles2d,
2246                                                                         mbspc.Knots(),
2247                                                                         mbspc.Multiplicities(),
2248                                                                         mbspc.Degree());
2249                   GeomLib_Check2dBSplineCurve newCheck(BS1,TOLCHECK,TOLANGCHECK);
2250                   newCheck.FixTangent(Standard_True,Standard_True);
2251                   //         
2252                   if (!reApprox) {
2253                     bIsValid1=CheckPCurve(BS1, myFace1);
2254                   }
2255                   //
2256                   aCurve.SetFirstCurve2d(BS1);
2257                 }
2258                 else {
2259                   Handle(Geom2d_BSplineCurve) BS1;
2260                   fprm = BS->FirstParameter();
2261                   lprm = BS->LastParameter();
2262
2263                   Handle(Geom2d_Curve) C2d;
2264                   Standard_Real aTol = myTolApprox;
2265                   BuildPCurves(fprm, lprm, aTol, myHS1->ChangeSurface().Surface(), BS, C2d);
2266                   BS1 = Handle(Geom2d_BSplineCurve)::DownCast(C2d);
2267                   aCurve.SetFirstCurve2d(BS1);
2268                 }
2269               } // if(myApprox1) { 
2270                 //                 
2271               if(myApprox2) { 
2272                 if(anApprox2) {
2273                   Handle(Geom2d_BSplineCurve) BS2;
2274                   TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
2275                   mbspc.Curve((myApprox1==Standard_True)? 3 : 2,tpoles2d);
2276                   BS2=new Geom2d_BSplineCurve(tpoles2d,
2277                                                                         mbspc.Knots(),
2278                                                                         mbspc.Multiplicities(),
2279                                                                         mbspc.Degree());
2280                   GeomLib_Check2dBSplineCurve newCheck(BS2,TOLCHECK,TOLANGCHECK);
2281                   newCheck.FixTangent(Standard_True,Standard_True);
2282                 //                 
2283                   if (!reApprox) {
2284                     bIsValid2=CheckPCurve(BS2, myFace2);        
2285                   }
2286                   aCurve.SetSecondCurve2d(BS2);
2287                 }
2288                 else {
2289                   Handle(Geom2d_BSplineCurve) BS2;
2290                   fprm = BS->FirstParameter();
2291                   lprm = BS->LastParameter();
2292
2293                   Handle(Geom2d_Curve) C2d;
2294                   Standard_Real aTol = myTolApprox;
2295                   BuildPCurves(fprm, lprm, aTol, myHS2->ChangeSurface().Surface(), BS, C2d);
2296                   BS2 = Handle(Geom2d_BSplineCurve)::DownCast(C2d);
2297                   aCurve.SetSecondCurve2d(BS2);
2298                 }
2299               } //if(myApprox2) { 
2300               if (!bIsValid1 || !bIsValid2) {
2301                 myTolApprox=aTolApproxImp;//1.e-5;
2302                 tol2d = myTolApprox;
2303                 reApprox = Standard_True;
2304                 goto reapprox;
2305               }
2306                 //                 
2307               mySeqOfCurve.Append(aCurve);
2308             }
2309           }
2310         }
2311       }
2312     }// else { // X
2313   }// case IntPatch_Walking:{
2314     break;
2315     
2316   case IntPatch_Restriction: 
2317     break;
2318   default:
2319     break;
2320
2321   }
2322 }
2323
2324 //=======================================================================
2325 //function : BuildPCurves
2326 //purpose  : 
2327 //=======================================================================
2328 void BuildPCurves (Standard_Real f,
2329                    Standard_Real l,
2330                    Standard_Real& Tol,
2331                    const Handle (Geom_Surface)& S,
2332                    const Handle (Geom_Curve)&   C,
2333                    Handle (Geom2d_Curve)& C2d)
2334 {
2335   if (!C2d.IsNull()) {
2336     return;
2337   }
2338   //
2339   Standard_Real umin,umax,vmin,vmax;
2340   // 
2341   S->Bounds(umin, umax, vmin, vmax);
2342   // in class ProjLib_Function the range of parameters is shrank by 1.e-09
2343   if((l - f) > 2.e-09) {
2344     C2d = GeomProjLib::Curve2d(C,f,l,S,umin,umax,vmin,vmax,Tol);
2345     //
2346     if (C2d.IsNull()) {
2347       // proj. a circle that goes through the pole on a sphere to the sphere     
2348       Tol += Precision::Confusion();
2349       C2d = GeomProjLib::Curve2d(C,f,l,S,Tol);
2350     }
2351   }
2352   else {
2353     if((l - f) > Epsilon(Abs(f))) {
2354       GeomAPI_ProjectPointOnSurf aProjector1, aProjector2;
2355       gp_Pnt P1 = C->Value(f);
2356       gp_Pnt P2 = C->Value(l);
2357       aProjector1.Init(P1, S);
2358       aProjector2.Init(P2, S);
2359       
2360       if(aProjector1.IsDone() && aProjector2.IsDone()) {
2361         Standard_Real U=0., V=0.;
2362         aProjector1.LowerDistanceParameters(U, V);
2363         gp_Pnt2d p1(U, V);
2364         
2365         aProjector2.LowerDistanceParameters(U, V);
2366         gp_Pnt2d p2(U, V);
2367         
2368         if(p1.Distance(p2) > gp::Resolution()) {
2369           TColgp_Array1OfPnt2d poles(1,2);
2370           TColStd_Array1OfReal knots(1,2);
2371           TColStd_Array1OfInteger mults(1,2);
2372           poles(1) = p1;
2373           poles(2) = p2;
2374           knots(1) = f;
2375           knots(2) = l;
2376           mults(1) = mults(2) = 2;
2377           
2378           C2d = new Geom2d_BSplineCurve(poles,knots,mults,1);
2379           
2380           // compute reached tolerance.begin
2381           gp_Pnt PMid = C->Value((f + l) * 0.5);
2382           aProjector1.Perform(PMid);
2383           
2384           if(aProjector1.IsDone()) {
2385             aProjector1.LowerDistanceParameters(U, V);
2386             gp_Pnt2d pmidproj(U, V);
2387             gp_Pnt2d pmidcurve2d = C2d->Value((f + l) * 0.5);
2388             Standard_Real adist = pmidcurve2d.Distance(pmidproj);
2389             Tol = (adist > Tol) ? adist : Tol;
2390           }
2391           // compute reached tolerance.end
2392         }
2393       }
2394     }
2395   }
2396   //
2397   if (S->IsUPeriodic() && !C2d.IsNull()) {
2398     // Recadre dans le domaine UV de la face
2399     Standard_Real aTm, U0, aEps, period, du, U0x;
2400     Standard_Boolean bAdjust;
2401     //
2402     aEps = Precision::PConfusion();
2403     period = S->UPeriod();
2404     //
2405     aTm = .5*(f + l);
2406     gp_Pnt2d pm = C2d->Value(aTm);
2407     U0 = pm.X();
2408     //
2409     bAdjust = 
2410       IntTools_Tools::AdjustPeriodic(U0, umin, umax, period, U0x, du, aEps);
2411     if (bAdjust) {
2412       gp_Vec2d T1(du, 0.);
2413       C2d->Translate(T1);
2414     }
2415   }
2416 }
2417
2418 //=======================================================================
2419 //function : Parameters
2420 //purpose  : 
2421 //=======================================================================
2422  void Parameters(const Handle(GeomAdaptor_HSurface)& HS1,
2423                  const Handle(GeomAdaptor_HSurface)& HS2,
2424                  const gp_Pnt& Ptref,
2425                  Standard_Real& U1,
2426                  Standard_Real& V1,
2427                  Standard_Real& U2,
2428                  Standard_Real& V2)
2429 {
2430
2431   IntSurf_Quadric quad1,quad2;
2432   GeomAbs_SurfaceType typs = HS1->Surface().GetType();
2433
2434   switch (typs) {
2435   case GeomAbs_Plane:
2436     quad1.SetValue(HS1->Surface().Plane());
2437     break;
2438   case GeomAbs_Cylinder:
2439     quad1.SetValue(HS1->Surface().Cylinder());
2440     break;
2441   case GeomAbs_Cone:
2442     quad1.SetValue(HS1->Surface().Cone());
2443     break;
2444   case GeomAbs_Sphere:
2445     quad1.SetValue(HS1->Surface().Sphere());
2446     break;
2447   case GeomAbs_Torus:
2448     quad1.SetValue(HS1->Surface().Torus());
2449     break;
2450   default:
2451     Standard_ConstructionError::Raise("GeomInt_IntSS::MakeCurve");
2452   }
2453   
2454   typs = HS2->Surface().GetType();
2455   switch (typs) {
2456   case GeomAbs_Plane:
2457     quad2.SetValue(HS2->Surface().Plane());
2458     break;
2459   case GeomAbs_Cylinder:
2460     quad2.SetValue(HS2->Surface().Cylinder());
2461     break;
2462   case GeomAbs_Cone:
2463     quad2.SetValue(HS2->Surface().Cone());
2464     break;
2465   case GeomAbs_Sphere:
2466     quad2.SetValue(HS2->Surface().Sphere());
2467     break;
2468   case GeomAbs_Torus:
2469     quad2.SetValue(HS2->Surface().Torus());
2470     break;
2471   default:
2472     Standard_ConstructionError::Raise("GeomInt_IntSS::MakeCurve");
2473   }
2474
2475   quad1.Parameters(Ptref,U1,V1);
2476   quad2.Parameters(Ptref,U2,V2);
2477 }
2478
2479 //=======================================================================
2480 //function : MakeBSpline
2481 //purpose  : 
2482 //=======================================================================
2483 Handle(Geom_Curve) MakeBSpline  (const Handle(IntPatch_WLine)& WL,
2484                                  const Standard_Integer ideb,
2485                                  const Standard_Integer ifin)
2486 {
2487   Standard_Integer i,nbpnt = ifin-ideb+1;
2488   TColgp_Array1OfPnt poles(1,nbpnt);
2489   TColStd_Array1OfReal knots(1,nbpnt);
2490   TColStd_Array1OfInteger mults(1,nbpnt);
2491   Standard_Integer ipidebm1;
2492   for(i=1,ipidebm1=i+ideb-1; i<=nbpnt;ipidebm1++, i++) {
2493     poles(i) = WL->Point(ipidebm1).Value();
2494     mults(i) = 1;
2495     knots(i) = i-1;
2496   }
2497   mults(1) = mults(nbpnt) = 2;
2498   return
2499     new Geom_BSplineCurve(poles,knots,mults,1);
2500 }
2501 //
2502
2503 //=======================================================================
2504 //function : MakeBSpline2d
2505 //purpose  : 
2506 //=======================================================================
2507 Handle(Geom2d_BSplineCurve) MakeBSpline2d(const Handle(IntPatch_WLine)& theWLine,
2508                                           const Standard_Integer ideb,
2509                                           const Standard_Integer ifin,
2510                                           const Standard_Boolean onFirst)
2511 {
2512   Standard_Integer i, nbpnt = ifin-ideb+1;
2513   TColgp_Array1OfPnt2d poles(1,nbpnt);
2514   TColStd_Array1OfReal knots(1,nbpnt);
2515   TColStd_Array1OfInteger mults(1,nbpnt);
2516   Standard_Integer ipidebm1;
2517
2518   for(i = 1, ipidebm1 = i+ideb-1; i <= nbpnt; ipidebm1++, i++) {
2519       Standard_Real U, V;
2520       if(onFirst)
2521         theWLine->Point(ipidebm1).ParametersOnS1(U, V);
2522       else
2523         theWLine->Point(ipidebm1).ParametersOnS2(U, V);
2524       poles(i).SetCoord(U, V);
2525       mults(i) = 1;
2526       knots(i) = i-1;
2527     }
2528     mults(1) = mults(nbpnt) = 2;
2529
2530   return new Geom2d_BSplineCurve(poles,knots,mults,1);
2531 }
2532 //=======================================================================
2533 //function : PrepareLines3D
2534 //purpose  : 
2535 //=======================================================================
2536   void IntTools_FaceFace::PrepareLines3D(const Standard_Boolean bToSplit)
2537 {
2538   Standard_Integer i, aNbCurves;
2539   GeomAbs_SurfaceType aType1, aType2;
2540   IntTools_SequenceOfCurves aNewCvs;
2541   //
2542   // 1. Treatment closed  curves
2543   aNbCurves=mySeqOfCurve.Length();
2544   for (i=1; i<=aNbCurves; ++i) {
2545     const IntTools_Curve& aIC=mySeqOfCurve(i);
2546     //
2547     if (bToSplit) {
2548       Standard_Integer j, aNbC;
2549       IntTools_SequenceOfCurves aSeqCvs;
2550       //
2551       aNbC=IntTools_Tools::SplitCurve(aIC, aSeqCvs);
2552       if (aNbC) {
2553         for (j=1; j<=aNbC; ++j) {
2554           const IntTools_Curve& aICNew=aSeqCvs(j);
2555           aNewCvs.Append(aICNew);
2556         }
2557       }
2558       else {
2559         aNewCvs.Append(aIC);
2560       }
2561     }
2562     else {
2563       aNewCvs.Append(aIC);
2564     }
2565   }
2566   //
2567   // 2. Plane\Cone intersection when we had 4 curves
2568   aType1=myHS1->GetType();
2569   aType2=myHS2->GetType();
2570   aNbCurves=aNewCvs.Length();
2571   //
2572   if ((aType1==GeomAbs_Plane && aType2==GeomAbs_Cone) ||
2573       (aType2==GeomAbs_Plane && aType1==GeomAbs_Cone)) {
2574     if (aNbCurves==4) {
2575       GeomAbs_CurveType aCType1;
2576       //
2577       aCType1=aNewCvs(1).Type();
2578       if (aCType1==GeomAbs_Line) {
2579         IntTools_SequenceOfCurves aSeqIn, aSeqOut;
2580         //
2581         for (i=1; i<=aNbCurves; ++i) {
2582           const IntTools_Curve& aIC=aNewCvs(i);
2583           aSeqIn.Append(aIC);
2584         }
2585         //
2586         IntTools_Tools::RejectLines(aSeqIn, aSeqOut);
2587         //
2588         aNewCvs.Clear();
2589         aNbCurves=aSeqOut.Length(); 
2590         for (i=1; i<=aNbCurves; ++i) {
2591           const IntTools_Curve& aIC=aSeqOut(i);
2592           aNewCvs.Append(aIC);
2593         }
2594       }
2595     }
2596   }// if ((aType1==GeomAbs_Plane && aType2==GeomAbs_Cone)...
2597   //
2598   // 3. Fill  mySeqOfCurve
2599   mySeqOfCurve.Clear();
2600   aNbCurves=aNewCvs.Length();
2601   for (i=1; i<=aNbCurves; ++i) {
2602     const IntTools_Curve& aIC=aNewCvs(i);
2603     mySeqOfCurve.Append(aIC);
2604   }
2605 }
2606 //=======================================================================
2607 //function : CorrectSurfaceBoundaries
2608 //purpose  : 
2609 //=======================================================================
2610  void CorrectSurfaceBoundaries(const TopoDS_Face&  theFace,
2611                               const Standard_Real theTolerance,
2612                               Standard_Real&      theumin,
2613                               Standard_Real&      theumax, 
2614                               Standard_Real&      thevmin, 
2615                               Standard_Real&      thevmax) 
2616 {
2617   Standard_Boolean enlarge, isuperiodic, isvperiodic;
2618   Standard_Real uinf, usup, vinf, vsup, delta;
2619   GeomAbs_SurfaceType aType;
2620   Handle(Geom_Surface) aSurface;
2621   //
2622   aSurface = BRep_Tool::Surface(theFace);
2623   aSurface->Bounds(uinf, usup, vinf, vsup);
2624   delta = theTolerance;
2625   enlarge = Standard_False;
2626   //
2627   GeomAdaptor_Surface anAdaptorSurface(aSurface);
2628   //
2629   if(aSurface->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface))) {
2630     Handle(Geom_Surface) aBasisSurface = 
2631       (Handle(Geom_RectangularTrimmedSurface)::DownCast(aSurface))->BasisSurface();
2632     
2633     if(aBasisSurface->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface)) ||
2634        aBasisSurface->IsKind(STANDARD_TYPE(Geom_OffsetSurface))) {
2635       return;
2636     }
2637   }
2638   //
2639   if(aSurface->IsKind(STANDARD_TYPE(Geom_OffsetSurface))) {
2640     Handle(Geom_Surface) aBasisSurface = 
2641       (Handle(Geom_OffsetSurface)::DownCast(aSurface))->BasisSurface();
2642     
2643     if(aBasisSurface->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface)) ||
2644        aBasisSurface->IsKind(STANDARD_TYPE(Geom_OffsetSurface))) {
2645       return;
2646     }
2647   }
2648   //
2649   isuperiodic = anAdaptorSurface.IsUPeriodic();
2650   isvperiodic = anAdaptorSurface.IsVPeriodic();
2651   //
2652   aType=anAdaptorSurface.GetType();
2653   if((aType==GeomAbs_BezierSurface) ||
2654      (aType==GeomAbs_BSplineSurface) ||
2655      (aType==GeomAbs_SurfaceOfExtrusion) ||
2656      (aType==GeomAbs_SurfaceOfRevolution) ||
2657      (aType==GeomAbs_Cylinder)) {
2658     enlarge=Standard_True;
2659   }
2660   //
2661   if(!isuperiodic && enlarge) {
2662
2663     if((theumin - uinf) > delta )
2664       theumin -= delta;
2665     else {
2666       theumin = uinf;
2667     }
2668
2669     if((usup - theumax) > delta )
2670       theumax += delta;
2671     else
2672       theumax = usup;
2673   }
2674   //
2675   if(!isvperiodic && enlarge) {
2676     if((thevmin - vinf) > delta ) {
2677       thevmin -= delta;
2678     }
2679     else { 
2680       thevmin = vinf;
2681     }
2682     if((vsup - thevmax) > delta ) {
2683       thevmax += delta;
2684     }
2685     else {
2686       thevmax = vsup;
2687     }
2688   }
2689   //
2690   {
2691     Standard_Integer aNbP;
2692     Standard_Real aXP, dXfact, aXmid, aX1, aX2, aTolPA;
2693     //
2694     aTolPA=Precision::Angular();
2695     // U
2696     if (isuperiodic) {
2697       aXP=anAdaptorSurface.UPeriod();
2698       dXfact=theumax-theumin;
2699       if (dXfact-aTolPA>aXP) {
2700         aXmid=0.5*(theumax+theumin);
2701         aNbP=RealToInt(aXmid/aXP);
2702         if (aXmid<0.) {
2703           aNbP=aNbP-1;
2704         }
2705         aX1=aNbP*aXP;
2706         if (theumin>aTolPA) {
2707           aX1=theumin+aNbP*aXP;
2708         }
2709         aX2=aX1+aXP;
2710         if (theumin<aX1) {
2711           theumin=aX1;
2712         }
2713         if (theumax>aX2) {
2714           theumax=aX2;
2715         }
2716       }
2717     }
2718     // V
2719     if (isvperiodic) {
2720       aXP=anAdaptorSurface.VPeriod();
2721       dXfact=thevmax-thevmin;
2722       if (dXfact-aTolPA>aXP) {
2723         aXmid=0.5*(thevmax+thevmin);
2724         aNbP=RealToInt(aXmid/aXP);
2725         if (aXmid<0.) {
2726           aNbP=aNbP-1;
2727         }
2728         aX1=aNbP*aXP;
2729         if (thevmin>aTolPA) {
2730           aX1=thevmin+aNbP*aXP;
2731         }
2732         aX2=aX1+aXP;
2733         if (thevmin<aX1) {
2734           thevmin=aX1;
2735         }
2736         if (thevmax>aX2) {
2737           thevmax=aX2;
2738         }
2739       }
2740     }
2741   }
2742   //
2743   if(isuperiodic || isvperiodic) {
2744     Standard_Boolean correct = Standard_False;
2745     Standard_Boolean correctU = Standard_False;
2746     Standard_Boolean correctV = Standard_False;
2747     Bnd_Box2d aBox;
2748     TopExp_Explorer anExp;
2749
2750     for(anExp.Init(theFace, TopAbs_EDGE); anExp.More(); anExp.Next()) {
2751       if(BRep_Tool::IsClosed(TopoDS::Edge(anExp.Current()), theFace)) {
2752         correct = Standard_True;
2753         Standard_Real f, l;
2754         TopoDS_Edge anEdge = TopoDS::Edge(anExp.Current());
2755         
2756         for(Standard_Integer i = 0; i < 2; i++) {
2757           if(i==0) {
2758             anEdge.Orientation(TopAbs_FORWARD);
2759           }
2760           else {
2761             anEdge.Orientation(TopAbs_REVERSED);
2762           }
2763           Handle(Geom2d_Curve) aCurve = BRep_Tool::CurveOnSurface(anEdge, theFace, f, l);
2764           
2765           if(aCurve.IsNull()) {
2766             correct = Standard_False;
2767             break;
2768           }
2769           Handle(Geom2d_Line) aLine = Handle(Geom2d_Line)::DownCast(aCurve);
2770
2771           if(aLine.IsNull()) {
2772             correct = Standard_False;
2773             break;
2774           }
2775           gp_Dir2d anUDir(1., 0.);
2776           gp_Dir2d aVDir(0., 1.);
2777           Standard_Real anAngularTolerance = Precision::Angular();
2778
2779           correctU = correctU || aLine->Position().Direction().IsParallel(aVDir, anAngularTolerance);
2780           correctV = correctV || aLine->Position().Direction().IsParallel(anUDir, anAngularTolerance);
2781           
2782           gp_Pnt2d pp1 = aCurve->Value(f);
2783           aBox.Add(pp1);
2784           gp_Pnt2d pp2 = aCurve->Value(l);
2785           aBox.Add(pp2);
2786         }
2787         if(!correct)
2788           break;
2789       }
2790     }
2791
2792     if(correct) {
2793       Standard_Real umin, vmin, umax, vmax;
2794       aBox.Get(umin, vmin, umax, vmax);
2795
2796       if(isuperiodic && correctU) {
2797         
2798         if(theumin < umin)
2799           theumin = umin;
2800         
2801         if(theumax > umax) {
2802           theumax = umax;
2803         }
2804       }
2805       if(isvperiodic && correctV) {
2806         
2807         if(thevmin < vmin)
2808           thevmin = vmin;
2809         if(thevmax > vmax)
2810           thevmax = vmax;
2811       }
2812     }
2813   }
2814 }
2815 //
2816 //
2817 // The block is dedicated to determine whether WLine [ifprm, ilprm]
2818 // crosses the degenerated zone on each given surface or not.
2819 // If Yes -> We will not use info about surfaces during approximation
2820 // because inside degenerated zone of the surface the approx. algo.
2821 // uses wrong values of normal, etc., and resulting curve will have
2822 // oscillations that we would not like to have. 
2823  
2824
2825
2826 static
2827   Standard_Boolean IsDegeneratedZone(const gp_Pnt2d& aP2d,
2828                                      const Handle(Geom_Surface)& aS,
2829                                      const Standard_Integer iDir);
2830 static
2831   Standard_Boolean IsPointInDegeneratedZone(const IntSurf_PntOn2S& aP2S,
2832                                             const TopoDS_Face& aF1,
2833                                             const TopoDS_Face& aF2);
2834 //=======================================================================
2835 //function :  NotUseSurfacesForApprox
2836 //purpose  : 
2837 //=======================================================================
2838 Standard_Boolean NotUseSurfacesForApprox(const TopoDS_Face& aF1,
2839                                          const TopoDS_Face& aF2,
2840                                          const Handle(IntPatch_WLine)& WL,
2841                                          const Standard_Integer ifprm,
2842                                          const Standard_Integer ilprm)
2843 {
2844   Standard_Boolean bPInDZ;
2845
2846   Handle(IntSurf_LineOn2S) aLineOn2S=WL->Curve();
2847   
2848   const IntSurf_PntOn2S& aP2Sfprm=aLineOn2S->Value(ifprm);
2849   bPInDZ=IsPointInDegeneratedZone(aP2Sfprm, aF1, aF2);
2850   if (bPInDZ) {
2851     return bPInDZ;
2852   }
2853
2854   const IntSurf_PntOn2S& aP2Slprm=aLineOn2S->Value(ilprm);
2855   bPInDZ=IsPointInDegeneratedZone(aP2Slprm, aF1, aF2);
2856   
2857   return bPInDZ;
2858 }
2859 //=======================================================================
2860 //function : IsPointInDegeneratedZone
2861 //purpose  : 
2862 //=======================================================================
2863 Standard_Boolean IsPointInDegeneratedZone(const IntSurf_PntOn2S& aP2S,
2864                                           const TopoDS_Face& aF1,
2865                                           const TopoDS_Face& aF2)
2866                                           
2867 {
2868   Standard_Boolean bFlag=Standard_True;
2869   Standard_Real US11, US12, VS11, VS12, US21, US22, VS21, VS22;
2870   Standard_Real U1, V1, U2, V2, aDelta, aD;
2871   gp_Pnt2d aP2d;
2872
2873   Handle(Geom_Surface)aS1 = BRep_Tool::Surface(aF1);
2874   aS1->Bounds(US11, US12, VS11, VS12);
2875   GeomAdaptor_Surface aGAS1(aS1);
2876
2877   Handle(Geom_Surface)aS2 = BRep_Tool::Surface(aF2);
2878   aS1->Bounds(US21, US22, VS21, VS22);
2879   GeomAdaptor_Surface aGAS2(aS2);
2880   //
2881   //const gp_Pnt& aP=aP2S.Value();
2882   aP2S.Parameters(U1, V1, U2, V2);
2883   //
2884   aDelta=1.e-7;
2885   // Check on Surf 1
2886   aD=aGAS1.UResolution(aDelta);
2887   aP2d.SetCoord(U1, V1);
2888   if (fabs(U1-US11) < aD) {
2889     bFlag=IsDegeneratedZone(aP2d, aS1, 1);
2890     if (bFlag) {
2891       return bFlag;
2892     }
2893   }
2894   if (fabs(U1-US12) < aD) {
2895     bFlag=IsDegeneratedZone(aP2d, aS1, 1);
2896     if (bFlag) {
2897       return bFlag;
2898     }
2899   }
2900   aD=aGAS1.VResolution(aDelta);
2901   if (fabs(V1-VS11) < aDelta) {
2902     bFlag=IsDegeneratedZone(aP2d, aS1, 2);
2903     if (bFlag) {
2904       return bFlag;
2905     }
2906   }
2907   if (fabs(V1-VS12) < aDelta) {
2908     bFlag=IsDegeneratedZone(aP2d, aS1, 2);
2909     if (bFlag) {
2910       return bFlag;
2911     }
2912   }
2913   // Check on Surf 2
2914   aD=aGAS2.UResolution(aDelta);
2915   aP2d.SetCoord(U2, V2);
2916   if (fabs(U2-US21) < aDelta) {
2917     bFlag=IsDegeneratedZone(aP2d, aS2, 1);
2918     if (bFlag) {
2919       return bFlag;
2920     }
2921   }
2922   if (fabs(U2-US22) < aDelta) {
2923     bFlag=IsDegeneratedZone(aP2d, aS2, 1);
2924     if (bFlag) {
2925       return bFlag;
2926     }
2927   }
2928   aD=aGAS2.VResolution(aDelta);
2929   if (fabs(V2-VS21) < aDelta) {
2930     bFlag=IsDegeneratedZone(aP2d, aS2, 2);
2931     if (bFlag) {  
2932       return bFlag;
2933     }
2934   }
2935   if (fabs(V2-VS22) < aDelta) {
2936     bFlag=IsDegeneratedZone(aP2d, aS2, 2);
2937     if (bFlag) {
2938       return bFlag;
2939     }
2940   }
2941   return !bFlag;
2942 }
2943
2944 //=======================================================================
2945 //function : IsDegeneratedZone
2946 //purpose  : 
2947 //=======================================================================
2948 Standard_Boolean IsDegeneratedZone(const gp_Pnt2d& aP2d,
2949                                    const Handle(Geom_Surface)& aS,
2950                                    const Standard_Integer iDir)
2951 {
2952   Standard_Boolean bFlag=Standard_True;
2953   Standard_Real US1, US2, VS1, VS2, dY, dX, d1, d2, dD;
2954   Standard_Real aXm, aYm, aXb, aYb, aXe, aYe;
2955   aS->Bounds(US1, US2, VS1, VS2); 
2956
2957   gp_Pnt aPm, aPb, aPe;
2958   
2959   aXm=aP2d.X();
2960   aYm=aP2d.Y();
2961   
2962   aS->D0(aXm, aYm, aPm); 
2963   
2964   dX=1.e-5;
2965   dY=1.e-5;
2966   dD=1.e-12;
2967
2968   if (iDir==1) {
2969     aXb=aXm;
2970     aXe=aXm;
2971     aYb=aYm-dY;
2972     if (aYb < VS1) {
2973       aYb=VS1;
2974     }
2975     aYe=aYm+dY;
2976     if (aYe > VS2) {
2977       aYe=VS2;
2978     }
2979     aS->D0(aXb, aYb, aPb);
2980     aS->D0(aXe, aYe, aPe);
2981     
2982     d1=aPm.Distance(aPb);
2983     d2=aPm.Distance(aPe);
2984     if (d1 < dD && d2 < dD) {
2985       return bFlag;
2986     }
2987     return !bFlag;
2988   }
2989   //
2990   else if (iDir==2) {
2991     aYb=aYm;
2992     aYe=aYm;
2993     aXb=aXm-dX;
2994     if (aXb < US1) {
2995       aXb=US1;
2996     }
2997     aXe=aXm+dX;
2998     if (aXe > US2) {
2999       aXe=US2;
3000     }
3001     aS->D0(aXb, aYb, aPb);
3002     aS->D0(aXe, aYe, aPe);
3003     
3004     d1=aPm.Distance(aPb);
3005     d2=aPm.Distance(aPe);
3006     if (d1 < dD && d2 < dD) {
3007       return bFlag;
3008     }
3009     return !bFlag;
3010   }
3011   return !bFlag;
3012 }
3013
3014 //=========================================================================
3015 // static function : ComputePurgedWLine
3016 // purpose : Removes equal points (leave one of equal points) from theWLine
3017 //           and recompute vertex parameters.
3018 //           Returns new WLine or null WLine if the number
3019 //           of the points is less than 2.
3020 //=========================================================================
3021 Handle(IntPatch_WLine) ComputePurgedWLine(const Handle(IntPatch_WLine)& theWLine) {
3022  
3023   Standard_Integer i, k, v, nb, nbvtx;
3024   Handle(IntPatch_WLine) aResult;
3025   nbvtx = theWLine->NbVertex();
3026   nb = theWLine->NbPnts();
3027   if (nb==2) {
3028     const IntSurf_PntOn2S& p1 = theWLine->Point(1);
3029     const IntSurf_PntOn2S& p2 = theWLine->Point(2);
3030     if(p1.Value().IsEqual(p2.Value(), gp::Resolution())) {
3031       return aResult;
3032     }
3033   }
3034   //
3035   Handle(IntPatch_WLine) aLocalWLine;
3036   Handle(IntPatch_WLine) aTmpWLine = theWLine;
3037   Handle(IntSurf_LineOn2S) aLineOn2S = new IntSurf_LineOn2S();
3038   aLocalWLine = new IntPatch_WLine(aLineOn2S, Standard_False);
3039   for(i = 1; i <= nb; i++) {
3040     aLineOn2S->Add(theWLine->Point(i));
3041   }
3042
3043   for(v = 1; v <= nbvtx; v++) {
3044     aLocalWLine->AddVertex(theWLine->Vertex(v));
3045   }
3046   
3047   for(i = 1; i <= aLineOn2S->NbPoints(); i++) {
3048     Standard_Integer aStartIndex = i + 1;
3049     Standard_Integer anEndIndex = i + 5;
3050     nb = aLineOn2S->NbPoints();
3051     anEndIndex = (anEndIndex > nb) ? nb : anEndIndex;
3052
3053     if((aStartIndex > nb) || (anEndIndex <= 1)) {
3054       continue;
3055     }
3056     k = aStartIndex;
3057
3058     while(k <= anEndIndex) {
3059       
3060       if(i != k) {
3061         IntSurf_PntOn2S p1 = aLineOn2S->Value(i);
3062         IntSurf_PntOn2S p2 = aLineOn2S->Value(k);
3063         
3064         if(p1.Value().IsEqual(p2.Value(), gp::Resolution())) {
3065           aTmpWLine = aLocalWLine;
3066           aLocalWLine = new IntPatch_WLine(aLineOn2S, Standard_False);
3067
3068           for(v = 1; v <= aTmpWLine->NbVertex(); v++) {
3069             IntPatch_Point aVertex = aTmpWLine->Vertex(v);
3070             Standard_Integer avertexindex = (Standard_Integer)aVertex.ParameterOnLine();
3071
3072             if(avertexindex >= k) {
3073               aVertex.SetParameter(aVertex.ParameterOnLine() - 1.);
3074             }
3075             aLocalWLine->AddVertex(aVertex);
3076           }
3077           aLineOn2S->RemovePoint(k);
3078           anEndIndex--;
3079           continue;
3080         }
3081       }
3082       k++;
3083     }
3084   }
3085
3086   if(aLineOn2S->NbPoints() > 1) {
3087     aResult = aLocalWLine;
3088   }
3089   return aResult;
3090 }
3091
3092 //=======================================================================
3093 //function : TolR3d
3094 //purpose  : 
3095 //=======================================================================
3096 void TolR3d(const TopoDS_Face& aF1,
3097             const TopoDS_Face& aF2,
3098             Standard_Real& myTolReached3d)
3099 {
3100   Standard_Real aTolF1, aTolF2, aTolFMax, aTolTresh;
3101       
3102   aTolTresh=2.999999e-3;
3103   aTolF1 = BRep_Tool::Tolerance(aF1);
3104   aTolF2 = BRep_Tool::Tolerance(aF2);
3105   aTolFMax=Max(aTolF1, aTolF2);
3106   
3107   if (aTolFMax>aTolTresh) {
3108     myTolReached3d=aTolFMax;
3109   }
3110 }
3111 //=======================================================================
3112 //function : IsPointOnBoundary
3113 //purpose  : 
3114 //=======================================================================
3115 Standard_Boolean IsPointOnBoundary(const Standard_Real theParameter,
3116                                    const Standard_Real theFirstBoundary,
3117                                    const Standard_Real theSecondBoundary,
3118                                    const Standard_Real theResolution,
3119                                    Standard_Boolean&   IsOnFirstBoundary) 
3120 {
3121   Standard_Boolean bRet;
3122   Standard_Integer i;
3123   Standard_Real adist;
3124   //
3125   bRet=Standard_False;
3126   for(i = 0; i < 2; ++i) {
3127     IsOnFirstBoundary = (i == 0);
3128     if (IsOnFirstBoundary) {
3129       adist = fabs(theParameter - theFirstBoundary);
3130     }
3131     else {
3132       adist = fabs(theParameter - theSecondBoundary);
3133     }
3134     if(adist < theResolution) {
3135       return !bRet;
3136     }
3137   }
3138   return bRet;
3139 }
3140 // ------------------------------------------------------------------------------------------------
3141 // static function: FindPoint
3142 // purpose:
3143 // ------------------------------------------------------------------------------------------------
3144 Standard_Boolean FindPoint(const gp_Pnt2d&     theFirstPoint,
3145                            const gp_Pnt2d&     theLastPoint,
3146                            const Standard_Real theUmin, 
3147                            const Standard_Real theUmax,
3148                            const Standard_Real theVmin,
3149                            const Standard_Real theVmax,
3150                            gp_Pnt2d&           theNewPoint) {
3151   
3152   gp_Vec2d aVec(theFirstPoint, theLastPoint);
3153   Standard_Integer i = 0, j = 0;
3154
3155   for(i = 0; i < 4; i++) {
3156     gp_Vec2d anOtherVec;
3157     gp_Vec2d anOtherVecNormal;
3158     gp_Pnt2d aprojpoint = theLastPoint;    
3159
3160     if((i % 2) == 0) {
3161       anOtherVec.SetX(0.);
3162       anOtherVec.SetY(1.);
3163       anOtherVecNormal.SetX(1.);
3164       anOtherVecNormal.SetY(0.);
3165
3166       if(i < 2)
3167         aprojpoint.SetX(theUmin);
3168       else
3169         aprojpoint.SetX(theUmax);
3170     }
3171     else {
3172       anOtherVec.SetX(1.);
3173       anOtherVec.SetY(0.);
3174       anOtherVecNormal.SetX(0.);
3175       anOtherVecNormal.SetY(1.);
3176
3177       if(i < 2)
3178         aprojpoint.SetY(theVmin);
3179       else
3180         aprojpoint.SetY(theVmax);
3181     }
3182     gp_Vec2d anormvec = aVec;
3183     anormvec.Normalize();
3184     RefineVector(anormvec);
3185     Standard_Real adot1 = anormvec.Dot(anOtherVecNormal);
3186
3187     if(fabs(adot1) < Precision::Angular())
3188       continue;
3189     Standard_Real adist = 0.;
3190     Standard_Boolean bIsOut = Standard_False;
3191
3192     if((i % 2) == 0) {
3193       adist = (i < 2) ? fabs(theLastPoint.X() - theUmin) : fabs(theLastPoint.X() - theUmax);
3194       bIsOut = (i < 2) ? (theLastPoint.X() < theUmin) : (theLastPoint.X() > theUmax);
3195     }
3196     else {
3197       adist = (i < 2) ? fabs(theLastPoint.Y() - theVmin) : fabs(theLastPoint.Y() - theVmax);
3198       bIsOut = (i < 2) ? (theLastPoint.Y() < theVmin) : (theLastPoint.Y() > theVmax);
3199     }
3200     Standard_Real anoffset = adist * anOtherVec.Dot(anormvec) / adot1;
3201
3202     for(j = 0; j < 2; j++) {
3203       anoffset = (j == 0) ? anoffset : -anoffset;
3204       gp_Pnt2d acurpoint(aprojpoint.XY() + (anOtherVec.XY()*anoffset));
3205       gp_Vec2d acurvec(theLastPoint, acurpoint);
3206       if ( bIsOut )
3207         acurvec.Reverse();
3208
3209       Standard_Real aDotX, anAngleX;
3210       //
3211       aDotX = aVec.Dot(acurvec);
3212       anAngleX = aVec.Angle(acurvec);
3213       //
3214       if(aDotX > 0. && fabs(anAngleX) < Precision::PConfusion()) {
3215         if((i % 2) == 0) {
3216           if((acurpoint.Y() >= theVmin) &&
3217              (acurpoint.Y() <= theVmax)) {
3218             theNewPoint = acurpoint;
3219             return Standard_True;
3220           }
3221         }
3222         else {
3223           if((acurpoint.X() >= theUmin) &&
3224              (acurpoint.X() <= theUmax)) {
3225             theNewPoint = acurpoint;
3226             return Standard_True;
3227           }
3228         }
3229       }
3230     }
3231   }
3232   return Standard_False;
3233 }
3234
3235
3236 // ------------------------------------------------------------------------------------------------
3237 // static function: FindPoint
3238 // purpose: Find point on the boundary of radial tangent zone
3239 // ------------------------------------------------------------------------------------------------
3240 Standard_Boolean FindPoint(const gp_Pnt2d&     theFirstPoint,
3241                            const gp_Pnt2d&     theLastPoint,
3242                            const Standard_Real theUmin, 
3243                            const Standard_Real theUmax,
3244                            const Standard_Real theVmin,
3245                            const Standard_Real theVmax,
3246                            const gp_Pnt2d&     theTanZoneCenter,
3247                            const Standard_Real theZoneRadius,
3248                            Handle(GeomAdaptor_HSurface) theGASurface,
3249                            gp_Pnt2d&           theNewPoint) {
3250   theNewPoint = theLastPoint;
3251
3252   if ( !IsInsideTanZone( theLastPoint, theTanZoneCenter, theZoneRadius, theGASurface) )
3253     return Standard_False;
3254
3255   Standard_Real aUResolution = theGASurface->UResolution( theZoneRadius );
3256   Standard_Real aVResolution = theGASurface->VResolution( theZoneRadius );
3257
3258   Standard_Real aRadius = ( aUResolution < aVResolution ) ? aUResolution : aVResolution;
3259   gp_Ax22d anAxis( theTanZoneCenter, gp_Dir2d(1, 0), gp_Dir2d(0, 1) );
3260   gp_Circ2d aCircle( anAxis, aRadius );
3261   
3262   //
3263   gp_Vec2d aDir( theLastPoint.XY() - theFirstPoint.XY() );
3264   Standard_Real aLength = aDir.Magnitude();
3265   if ( aLength <= gp::Resolution() )
3266     return Standard_False;
3267   gp_Lin2d aLine( theFirstPoint, aDir );
3268
3269   //
3270   Handle(Geom2d_Line) aCLine = new Geom2d_Line( aLine );
3271   Handle(Geom2d_TrimmedCurve) aC1 = new Geom2d_TrimmedCurve( aCLine, 0, aLength );
3272   Handle(Geom2d_Circle) aC2 = new Geom2d_Circle( aCircle );
3273
3274   Standard_Real aTol = aRadius * 0.001;
3275   aTol = ( aTol < Precision::PConfusion() ) ? Precision::PConfusion() : aTol;
3276
3277   Geom2dAPI_InterCurveCurve anIntersector;
3278   anIntersector.Init( aC1, aC2, aTol );
3279
3280   if ( anIntersector.NbPoints() == 0 )
3281     return Standard_False;
3282
3283   Standard_Boolean aFound = Standard_False;
3284   Standard_Real aMinDist = aLength * aLength;
3285   Standard_Integer i = 0;
3286   for ( i = 1; i <= anIntersector.NbPoints(); i++ ) {
3287     gp_Pnt2d aPInt = anIntersector.Point( i );
3288     if ( aPInt.SquareDistance( theFirstPoint ) < aMinDist ) {
3289       if ( ( aPInt.X() >= theUmin ) && ( aPInt.X() <= theUmax ) &&
3290            ( aPInt.Y() >= theVmin ) && ( aPInt.Y() <= theVmax ) ) {
3291         theNewPoint = aPInt;
3292         aFound = Standard_True;
3293       }
3294     }
3295   }
3296
3297   return aFound;
3298 }
3299
3300 // ------------------------------------------------------------------------------------------------
3301 // static function: IsInsideTanZone
3302 // purpose: Check if point is inside a radial tangent zone
3303 // ------------------------------------------------------------------------------------------------
3304 Standard_Boolean IsInsideTanZone(const gp_Pnt2d&     thePoint,
3305                                  const gp_Pnt2d&     theTanZoneCenter,
3306                                  const Standard_Real theZoneRadius,
3307                                  Handle(GeomAdaptor_HSurface) theGASurface) {
3308
3309   Standard_Real aUResolution = theGASurface->UResolution( theZoneRadius );
3310   Standard_Real aVResolution = theGASurface->VResolution( theZoneRadius );
3311   Standard_Real aRadiusSQR = ( aUResolution < aVResolution ) ? aUResolution : aVResolution;
3312   aRadiusSQR *= aRadiusSQR;
3313   if ( thePoint.SquareDistance( theTanZoneCenter ) <= aRadiusSQR )
3314     return Standard_True;
3315   return Standard_False;
3316 }
3317
3318 // ------------------------------------------------------------------------------------------------
3319 // static function: CheckTangentZonesExist
3320 // purpose: Check if tangent zone exists
3321 // ------------------------------------------------------------------------------------------------
3322 Standard_Boolean CheckTangentZonesExist( const Handle(GeomAdaptor_HSurface)& theSurface1,
3323                                         const Handle(GeomAdaptor_HSurface)&  theSurface2 ) 
3324 {
3325   if ( ( theSurface1->GetType() != GeomAbs_Torus ) ||
3326       ( theSurface2->GetType() != GeomAbs_Torus ) )
3327     return Standard_False;
3328
3329   gp_Torus aTor1 = theSurface1->Torus();
3330   gp_Torus aTor2 = theSurface2->Torus();
3331
3332   if ( aTor1.Location().Distance( aTor2.Location() ) > Precision::Confusion() )
3333     return Standard_False;
3334
3335   if ( ( fabs( aTor1.MajorRadius() - aTor2.MajorRadius() ) > Precision::Confusion() ) ||
3336        ( fabs( aTor1.MinorRadius() - aTor2.MinorRadius() ) > Precision::Confusion() ) )
3337     return Standard_False;
3338
3339   if ( ( aTor1.MajorRadius() < aTor1.MinorRadius() ) ||
3340        ( aTor2.MajorRadius() < aTor2.MinorRadius() ) )
3341     return Standard_False;
3342   return Standard_True;
3343 }
3344
3345 // ------------------------------------------------------------------------------------------------
3346 // static function: ComputeTangentZones
3347 // purpose: 
3348 // ------------------------------------------------------------------------------------------------
3349 Standard_Integer ComputeTangentZones( const Handle(GeomAdaptor_HSurface)& theSurface1,
3350                                      const Handle(GeomAdaptor_HSurface)&  theSurface2,
3351                                      const TopoDS_Face&                   theFace1,
3352                                      const TopoDS_Face&                   theFace2,
3353                                      Handle(TColgp_HArray1OfPnt2d)&       theResultOnS1,
3354                                      Handle(TColgp_HArray1OfPnt2d)&       theResultOnS2,
3355                                      Handle(TColStd_HArray1OfReal)&       theResultRadius,
3356                                      const Handle(IntTools_Context)& aContext)
3357 {
3358   Standard_Integer aResult = 0;
3359   if ( !CheckTangentZonesExist( theSurface1, theSurface2 ) )
3360     return aResult;
3361
3362
3363   TColgp_SequenceOfPnt2d aSeqResultS1, aSeqResultS2;
3364   TColStd_SequenceOfReal aSeqResultRad;
3365
3366   gp_Torus aTor1 = theSurface1->Torus();
3367   gp_Torus aTor2 = theSurface2->Torus();
3368
3369   gp_Ax2 anax1( aTor1.Location(), aTor1.Axis().Direction() );
3370   gp_Ax2 anax2( aTor2.Location(), aTor2.Axis().Direction() );
3371   Standard_Integer j = 0;
3372
3373   for ( j = 0; j < 2; j++ ) {
3374     Standard_Real aCoef = ( j == 0 ) ? -1 : 1;
3375     Standard_Real aRadius1 = fabs(aTor1.MajorRadius() + aCoef * aTor1.MinorRadius());
3376     Standard_Real aRadius2 = fabs(aTor2.MajorRadius() + aCoef * aTor2.MinorRadius());
3377
3378     gp_Circ aCircle1( anax1, aRadius1 );
3379     gp_Circ aCircle2( anax2, aRadius2 );
3380
3381     // roughly compute radius of tangent zone for perpendicular case
3382     Standard_Real aCriteria = Precision::Confusion() * 0.5;
3383
3384     Standard_Real aT1 = aCriteria;
3385     Standard_Real