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