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