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