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