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