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