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