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