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