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