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