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