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