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