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