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