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