bd1e4f1dab0941a3725f0e0cd34484e414e892d9
[occt.git] / src / IntTools / IntTools_FaceFace.cxx
1 // Created on: 2000-11-23
2 // Created by: Michael KLOKOV
3 // Copyright (c) 2000-2012 OPEN CASCADE SAS
4 //
5 // The content of this file is subject to the Open CASCADE Technology Public
6 // License Version 6.5 (the "License"). You may not use the content of this file
7 // except in compliance with the License. Please obtain a copy of the License
8 // at http://www.opencascade.org and read it completely before using this file.
9 //
10 // The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
11 // main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
12 //
13 // The Original Code and all software distributed under the License is
14 // distributed on an "AS IS" basis, without warranty of any kind, and the
15 // Initial Developer hereby disclaims all such warranties, including without
16 // limitation, any warranties of merchantability, fitness for a particular
17 // purpose or non-infringement. Please see the License for the specific terms
18 // and conditions governing the rights and limitations under the License.
19
20
21
22 #include <IntTools_FaceFace.ixx>
23
24 #include <Precision.hxx>
25
26 #include <TColStd_HArray1OfReal.hxx>
27 #include <TColStd_Array1OfReal.hxx>
28 #include <TColStd_Array1OfInteger.hxx>
29 #include <TColStd_SequenceOfReal.hxx>
30 #include <TColStd_ListOfInteger.hxx>
31 #include <TColStd_ListIteratorOfListOfInteger.hxx>
32 #include <TColStd_Array1OfListOfInteger.hxx>
33
34 #include <gp_Lin2d.hxx>
35 #include <gp_Ax22d.hxx>
36 #include <gp_Circ2d.hxx>
37 #include <gp_Torus.hxx>
38 #include <gp_Cylinder.hxx>
39
40 #include <Bnd_Box.hxx>
41
42 #include <TColgp_HArray1OfPnt2d.hxx>
43 #include <TColgp_SequenceOfPnt2d.hxx>
44 #include <TColgp_Array1OfPnt.hxx>
45 #include <TColgp_Array1OfPnt2d.hxx>
46
47 #include <IntAna_QuadQuadGeo.hxx>
48
49 #include <IntSurf_PntOn2S.hxx>
50 #include <IntSurf_LineOn2S.hxx>
51 #include <IntSurf_PntOn2S.hxx>
52 #include <IntSurf_ListOfPntOn2S.hxx>
53 #include <IntRes2d_Domain.hxx>
54 #include <ProjLib_Plane.hxx>
55
56 #include <IntPatch_GLine.hxx>
57 #include <IntPatch_RLine.hxx>
58 #include <IntPatch_WLine.hxx>
59 #include <IntPatch_ALine.hxx>
60 #include <IntPatch_ALineToWLine.hxx>
61
62 #include <ElSLib.hxx>
63 #include <ElCLib.hxx>
64
65 #include <Extrema_ExtCC.hxx>
66 #include <Extrema_POnCurv.hxx>
67 #include <BndLib_AddSurface.hxx>
68
69 #include <Adaptor3d_SurfacePtr.hxx>
70 #include <Adaptor2d_HLine2d.hxx>
71
72 #include <GeomAbs_SurfaceType.hxx>
73 #include <GeomAbs_CurveType.hxx>
74
75 #include <Geom_Surface.hxx>
76 #include <Geom_Line.hxx>
77 #include <Geom_Circle.hxx>
78 #include <Geom_Ellipse.hxx>
79 #include <Geom_Parabola.hxx>
80 #include <Geom_Hyperbola.hxx>
81 #include <Geom_TrimmedCurve.hxx>
82 #include <Geom_BSplineCurve.hxx>
83 #include <Geom_RectangularTrimmedSurface.hxx>
84 #include <Geom_OffsetSurface.hxx>
85 #include <Geom_Curve.hxx>
86 #include <Geom_Conic.hxx>
87
88 #include <Geom2d_TrimmedCurve.hxx>
89 #include <Geom2d_BSplineCurve.hxx>
90 #include <Geom2d_Line.hxx>
91 #include <Geom2d_Curve.hxx>
92 #include <Geom2d_Circle.hxx>
93
94 #include <Geom2dAPI_InterCurveCurve.hxx>
95 #include <Geom2dInt_GInter.hxx>
96 #include <GeomAdaptor_Curve.hxx>
97 #include <GeomAdaptor_HSurface.hxx>
98 #include <GeomAdaptor_Surface.hxx>
99 #include <GeomLib_CheckBSplineCurve.hxx>
100 #include <GeomLib_Check2dBSplineCurve.hxx>
101
102 #include <GeomInt_WLApprox.hxx>
103 #include <GeomProjLib.hxx>
104 #include <GeomAPI_ProjectPointOnSurf.hxx>
105 #include <Geom2dAdaptor_Curve.hxx>
106 #include <TopoDS.hxx>
107 #include <TopoDS_Edge.hxx>
108 #include <TopExp_Explorer.hxx>
109
110 #include <BRep_Tool.hxx>
111 #include <BRepTools.hxx>
112 #include <BRepAdaptor_Surface.hxx>
113
114 #include <IntTools_Curve.hxx>
115 #include <IntTools_Tools.hxx>
116 #include <IntTools_Tools.hxx>
117 #include <IntTools_TopolTool.hxx>
118 #include <IntTools_PntOnFace.hxx>
119 #include <IntTools_PntOn2Faces.hxx>
120 #include <BOPInt_Context.hxx>
121 #include <IntSurf_ListIteratorOfListOfPntOn2S.hxx>
122
123 static
124   void RefineVector(gp_Vec2d& aV2D);
125 #ifdef DEB_DUMPWLINE
126 static
127   void DumpWLine(const Handle(IntPatch_WLine)& aWLine);
128 #endif
129 //
130 static
131   void TolR3d(const TopoDS_Face& ,
132               const TopoDS_Face& ,
133               Standard_Real& );
134 static 
135   Handle(Geom_Curve) MakeBSpline (const Handle(IntPatch_WLine)&,
136                                   const Standard_Integer,
137                                   const Standard_Integer);
138
139 static 
140   void Parameters(const Handle(GeomAdaptor_HSurface)&,
141                   const Handle(GeomAdaptor_HSurface)&,
142                   const gp_Pnt&,
143                   Standard_Real&,
144                   Standard_Real&,
145                   Standard_Real&,
146                   Standard_Real&);
147
148 static 
149   void BuildPCurves (Standard_Real f,Standard_Real l,Standard_Real& Tol,
150                      const Handle (Geom_Surface)& S,
151                      const Handle (Geom_Curve)&   C,
152                      Handle (Geom2d_Curve)& C2d);
153
154 static 
155   void CorrectSurfaceBoundaries(const TopoDS_Face&  theFace,
156                                 const Standard_Real theTolerance,
157                                 Standard_Real&      theumin,
158                                 Standard_Real&      theumax, 
159                                 Standard_Real&      thevmin, 
160                                 Standard_Real&      thevmax);
161 static
162   Standard_Boolean NotUseSurfacesForApprox
163           (const TopoDS_Face& aF1,
164            const TopoDS_Face& aF2,
165            const Handle(IntPatch_WLine)& WL,
166            const Standard_Integer ifprm,
167            const Standard_Integer ilprm);
168
169 static 
170   Handle(IntPatch_WLine) ComputePurgedWLine(const Handle(IntPatch_WLine)& theWLine);
171
172 static 
173   Standard_Real AdjustPeriodic(const Standard_Real theParameter,
174                                const Standard_Real parmin,
175                                const Standard_Real parmax,
176                                const Standard_Real thePeriod,
177                                Standard_Real&      theOffset);
178
179 static 
180   Handle(Geom2d_BSplineCurve) MakeBSpline2d(const Handle(IntPatch_WLine)& theWLine,
181                                             const Standard_Integer                         ideb,
182                                             const Standard_Integer                         ifin,
183                                             const Standard_Boolean                         onFirst);
184
185 static 
186   Standard_Boolean DecompositionOfWLine(const Handle(IntPatch_WLine)& theWLine,
187                                         const Handle(GeomAdaptor_HSurface)&            theSurface1, 
188                                         const Handle(GeomAdaptor_HSurface)&            theSurface2,
189                                         const TopoDS_Face&                             theFace1,
190                                         const TopoDS_Face&                             theFace2,
191                                         const IntTools_LineConstructor&                theLConstructor,
192                                         const Standard_Boolean                         theAvoidLConstructor,
193                                         IntPatch_SequenceOfLine&                       theNewLines,
194                                         Standard_Real&                                 theReachedTol3d,
195                                         const Handle(BOPInt_Context)& );
196
197 static 
198   Standard_Boolean ParameterOutOfBoundary(const Standard_Real       theParameter, 
199                                           const Handle(Geom_Curve)& theCurve, 
200                                           const TopoDS_Face&        theFace1, 
201                                           const TopoDS_Face&        theFace2,
202                                           const Standard_Real       theOtherParameter,
203                                           const Standard_Boolean    bIncreasePar,
204                                           Standard_Real&            theNewParameter,
205                                           const Handle(BOPInt_Context)& );
206
207 static 
208   Standard_Boolean IsCurveValid(Handle(Geom2d_Curve)& thePCurve);
209
210 static 
211   Standard_Boolean IsPointOnBoundary(const Standard_Real theParameter,
212                                      const Standard_Real theFirstBoundary,
213                                      const Standard_Real theSecondBoundary,
214                                      const Standard_Real theResolution,
215                                      Standard_Boolean&   IsOnFirstBoundary);
216 static
217   Standard_Boolean FindPoint(const gp_Pnt2d&     theFirstPoint,
218                              const gp_Pnt2d&     theLastPoint,
219                              const Standard_Real theUmin, 
220                              const Standard_Real theUmax,
221                              const Standard_Real theVmin,
222                              const Standard_Real theVmax,
223                              gp_Pnt2d&           theNewPoint);
224
225
226 static 
227   Standard_Integer ComputeTangentZones( const Handle(GeomAdaptor_HSurface)& theSurface1,
228                                        const Handle(GeomAdaptor_HSurface)&  theSurface2,
229                                        const TopoDS_Face&                   theFace1,
230                                        const TopoDS_Face&                   theFace2,
231                                        Handle(TColgp_HArray1OfPnt2d)&       theResultOnS1,
232                                        Handle(TColgp_HArray1OfPnt2d)&       theResultOnS2,
233                                        Handle(TColStd_HArray1OfReal)&       theResultRadius,
234                                        const Handle(BOPInt_Context)& );
235
236 static
237   Standard_Boolean FindPoint(const gp_Pnt2d&     theFirstPoint,
238                              const gp_Pnt2d&     theLastPoint,
239                              const Standard_Real theUmin, 
240                              const Standard_Real theUmax,
241                              const Standard_Real theVmin,
242                              const Standard_Real theVmax,
243                              const gp_Pnt2d&     theTanZoneCenter,
244                              const Standard_Real theZoneRadius,
245                              Handle(GeomAdaptor_HSurface) theGASurface,
246                              gp_Pnt2d&           theNewPoint);
247
248 static
249   Standard_Boolean IsInsideTanZone(const gp_Pnt2d&     thePoint,
250                                    const gp_Pnt2d&     theTanZoneCenter,
251                                    const Standard_Real theZoneRadius,
252                                    Handle(GeomAdaptor_HSurface) theGASurface);
253
254 static
255   gp_Pnt2d AdjustByNeighbour(const gp_Pnt2d&     theaNeighbourPoint,
256                              const gp_Pnt2d&     theOriginalPoint,
257                              Handle(GeomAdaptor_HSurface) theGASurface);
258 static
259   Standard_Boolean  ApproxWithPCurves(const gp_Cylinder& theCyl, 
260                                       const gp_Sphere& theSph);
261
262 static void  PerformPlanes(const Handle(GeomAdaptor_HSurface)& theS1, 
263                            const Handle(GeomAdaptor_HSurface)& theS2, 
264                            const Standard_Real TolAng, 
265                            const Standard_Real TolTang, 
266                            const Standard_Boolean theApprox1,
267                            const Standard_Boolean theApprox2,
268                            IntTools_SequenceOfCurves& theSeqOfCurve, 
269                            Standard_Boolean& theTangentFaces);
270
271 static Standard_Boolean ClassifyLin2d(const Handle(GeomAdaptor_HSurface)& theS, 
272                                       const gp_Lin2d& theLin2d, 
273                                       const Standard_Real theTol,
274                                       Standard_Real& theP1, 
275                                       Standard_Real& theP2);
276 //
277 static
278   void ApproxParameters(const Handle(GeomAdaptor_HSurface)& aHS1,
279                         const Handle(GeomAdaptor_HSurface)& aHS2,
280                         Standard_Integer& iDegMin,
281                         Standard_Integer& iNbIter,
282                         Standard_Integer& iDegMax);
283
284 static
285   void Tolerances(const Handle(GeomAdaptor_HSurface)& aHS1,
286                   const Handle(GeomAdaptor_HSurface)& aHS2,
287                   Standard_Real& aTolTang);
288
289 static
290   Standard_Boolean SortTypes(const GeomAbs_SurfaceType aType1,
291                              const GeomAbs_SurfaceType aType2);
292 static
293   Standard_Integer IndexType(const GeomAbs_SurfaceType aType);
294
295 //
296 static
297   Standard_Real MaxSquareDistance (const Standard_Real aT,
298                                    const Handle(Geom_Curve)& aC3D,
299                                    const Handle(Geom2d_Curve)& aC2D1,
300                                    const Handle(Geom2d_Curve)& aC2D2,
301                                    const Handle(GeomAdaptor_HSurface) myHS1,
302                                    const Handle(GeomAdaptor_HSurface) myHS2,
303                                    const TopoDS_Face& aF1,
304                                    const TopoDS_Face& aF2,
305                                    const Handle(BOPInt_Context)& aCtx);
306
307 static
308   Standard_Boolean CheckPCurve(const Handle(Geom2d_Curve)& aPC, 
309                                const TopoDS_Face& aFace);
310
311 //
312 static
313   Standard_Real FindMaxSquareDistance (const Standard_Real aA,
314                                        const Standard_Real aB,
315                                        const Standard_Real aEps,
316                                        const Handle(Geom_Curve)& aC3D,
317                                        const Handle(Geom2d_Curve)& aC2D1,
318                                        const Handle(Geom2d_Curve)& aC2D2,
319                                        const Handle(GeomAdaptor_HSurface)& myHS1,
320                                        const Handle(GeomAdaptor_HSurface)& myHS2,
321                                        const TopoDS_Face& aF1,
322                                        const TopoDS_Face& aF2,
323                                        const Handle(BOPInt_Context)& aCtx);
324
325 //=======================================================================
326 //function : 
327 //purpose  : 
328 //=======================================================================
329 IntTools_FaceFace::IntTools_FaceFace()
330 {
331   myIsDone=Standard_False;
332   myTangentFaces=Standard_False;
333   //
334   myHS1 = new GeomAdaptor_HSurface ();
335   myHS2 = new GeomAdaptor_HSurface ();
336   myTolReached2d=0.; 
337   myTolReached3d=0.;
338   SetParameters(Standard_True, Standard_True, Standard_True, 1.e-07);
339   
340 }
341 //=======================================================================
342 //function : SetContext
343 //purpose  : 
344 //======================================================================= 
345 void IntTools_FaceFace::SetContext(const Handle(BOPInt_Context)& aContext)
346 {
347   myContext=aContext;
348 }
349 //=======================================================================
350 //function : Context
351 //purpose  : 
352 //======================================================================= 
353 const Handle(BOPInt_Context)& IntTools_FaceFace::Context()const
354 {
355   return myContext;
356 }
357 //=======================================================================
358 //function : Face1
359 //purpose  : 
360 //======================================================================= 
361 const TopoDS_Face& IntTools_FaceFace::Face1() const
362 {
363   return myFace1;
364 }
365 //=======================================================================
366 //function : Face2
367 //purpose  : 
368 //======================================================================= 
369 const TopoDS_Face& IntTools_FaceFace::Face2() const
370 {
371   return myFace2;
372 }
373 //=======================================================================
374 //function : TangentFaces
375 //purpose  : 
376 //======================================================================= 
377 Standard_Boolean IntTools_FaceFace::TangentFaces() const
378 {
379   return myTangentFaces;
380 }
381 //=======================================================================
382 //function : Points
383 //purpose  : 
384 //======================================================================= 
385 const IntTools_SequenceOfPntOn2Faces& IntTools_FaceFace::Points() const
386 {
387   return myPnts;
388 }
389 //=======================================================================
390 //function : IsDone
391 //purpose  : 
392 //======================================================================= 
393 Standard_Boolean IntTools_FaceFace::IsDone() const
394 {
395   return myIsDone;
396 }
397 //=======================================================================
398 //function : TolReached3d
399 //purpose  : 
400 //=======================================================================
401 Standard_Real IntTools_FaceFace::TolReached3d() const
402 {
403   return myTolReached3d;
404 }
405 //=======================================================================
406 //function : Lines
407 //purpose  : return lines of intersection
408 //=======================================================================
409 const IntTools_SequenceOfCurves& IntTools_FaceFace::Lines() const
410 {
411   StdFail_NotDone_Raise_if
412     (!myIsDone,
413      "IntTools_FaceFace::Lines() => myIntersector NOT DONE");
414   return mySeqOfCurve;
415 }
416 //=======================================================================
417 //function : TolReached2d
418 //purpose  : 
419 //=======================================================================
420 Standard_Real IntTools_FaceFace::TolReached2d() const
421 {
422   return myTolReached2d;
423 }
424 // =======================================================================
425 // function: SetParameters
426 //
427 // =======================================================================
428 void IntTools_FaceFace::SetParameters(const Standard_Boolean ToApproxC3d,
429                                       const Standard_Boolean ToApproxC2dOnS1,
430                                       const Standard_Boolean ToApproxC2dOnS2,
431                                       const Standard_Real ApproximationTolerance) 
432 {
433   myApprox = ToApproxC3d;
434   myApprox1 = ToApproxC2dOnS1;
435   myApprox2 = ToApproxC2dOnS2;
436   myTolApprox = ApproximationTolerance;
437 }
438 //=======================================================================
439 //function : SetList
440 //purpose  : 
441 //=======================================================================
442 void IntTools_FaceFace::SetList(IntSurf_ListOfPntOn2S& aListOfPnts)
443 {
444   myListOfPnts = aListOfPnts;  
445 }
446
447
448 static Standard_Boolean isTreatAnalityc(const TopoDS_Face& theF1,
449                                         const TopoDS_Face& theF2)
450 {
451   const Standard_Real Tolang = 1.e-8;
452   const Standard_Real aTolF1=BRep_Tool::Tolerance(theF1);
453   const Standard_Real aTolF2=BRep_Tool::Tolerance(theF2);
454   const Standard_Real aTolSum = aTolF1 + aTolF2;
455   Standard_Real aHigh = 0.0;
456
457   const BRepAdaptor_Surface aBAS1(theF1), aBAS2(theF2);
458   const GeomAbs_SurfaceType aType1=aBAS1.GetType();
459   const GeomAbs_SurfaceType aType2=aBAS2.GetType();
460   
461   gp_Pln aS1;
462   gp_Cylinder aS2;
463   if(aType1 == GeomAbs_Plane)
464   {
465     aS1=aBAS1.Plane();
466   }
467   else if(aType2 == GeomAbs_Plane)
468   {
469     aS1=aBAS2.Plane();
470   }
471   else
472   {
473     return Standard_True;
474   }
475
476   if(aType1 == GeomAbs_Cylinder)
477   {
478     aS2=aBAS1.Cylinder();
479     const Standard_Real VMin = aBAS1.FirstVParameter();
480     const Standard_Real VMax = aBAS1.LastVParameter();
481
482     if( Precision::IsNegativeInfinite(VMin) ||
483         Precision::IsPositiveInfinite(VMax))
484           return Standard_True;
485     else
486       aHigh = VMax - VMin;
487   }
488   else if(aType2 == GeomAbs_Cylinder)
489   {
490     aS2=aBAS2.Cylinder();
491
492     const Standard_Real VMin = aBAS2.FirstVParameter();
493     const Standard_Real VMax = aBAS2.LastVParameter();
494
495     if( Precision::IsNegativeInfinite(VMin) ||
496         Precision::IsPositiveInfinite(VMax))
497           return Standard_True;
498     else
499       aHigh = VMax - VMin;
500   }
501   else
502   {
503     return Standard_True;
504   }
505
506   IntAna_QuadQuadGeo inter;
507   inter.Perform(aS1,aS2,Tolang,aTolSum, aHigh);
508   if(inter.TypeInter() == IntAna_Ellipse)
509   {
510     const gp_Elips anEl = inter.Ellipse(1);
511     const Standard_Real aMajorR = anEl.MajorRadius();
512     const Standard_Real aMinorR = anEl.MinorRadius();
513     
514     return (aMajorR < 100000.0 * aMinorR);
515   }
516   else
517   {
518     return inter.IsDone();
519   }
520 }
521
522
523
524 //=======================================================================
525 //function : Perform
526 //purpose  : intersect surfaces of the faces
527 //=======================================================================
528   void IntTools_FaceFace::Perform(const TopoDS_Face& aF1,
529                                   const TopoDS_Face& aF2)
530 {
531   Standard_Boolean RestrictLine = Standard_False, hasCone = Standard_False;
532   
533   if (myContext.IsNull()) {
534     myContext=new BOPInt_Context;
535   }
536
537   mySeqOfCurve.Clear();
538   myTolReached2d=0.;
539   myTolReached3d=0.;
540   myIsDone = Standard_False;
541   myNbrestr=0;//?
542
543   myFace1=aF1;
544   myFace2=aF2;
545
546   const BRepAdaptor_Surface aBAS1(myFace1, Standard_False);
547   const BRepAdaptor_Surface aBAS2(myFace2, Standard_False);
548   GeomAbs_SurfaceType aType1=aBAS1.GetType();
549   GeomAbs_SurfaceType aType2=aBAS2.GetType();
550
551   const Standard_Boolean bReverse=SortTypes(aType1, aType2);
552   if (bReverse)
553   {
554     myFace1=aF2;
555     myFace2=aF1;
556     aType1=aBAS2.GetType();
557     aType2=aBAS1.GetType();
558
559     if (myListOfPnts.Extent())
560     {
561       Standard_Real aU1,aV1,aU2,aV2;
562       IntSurf_ListIteratorOfListOfPntOn2S aItP2S;
563       //
564       aItP2S.Initialize(myListOfPnts);
565       for (; aItP2S.More(); aItP2S.Next())
566       {
567         IntSurf_PntOn2S& aP2S=aItP2S.Value();
568         aP2S.Parameters(aU1,aV1,aU2,aV2);
569         aP2S.SetValue(aU2,aV2,aU1,aV1);
570       }
571     }
572   }
573
574
575   const Handle(Geom_Surface) S1=BRep_Tool::Surface(myFace1);
576   const Handle(Geom_Surface) S2=BRep_Tool::Surface(myFace2);
577
578   const Standard_Real aTolF1=BRep_Tool::Tolerance(myFace1);
579   const Standard_Real aTolF2=BRep_Tool::Tolerance(myFace2);
580
581   Standard_Real TolArc = aTolF1 + aTolF2;
582   Standard_Real TolTang = TolArc;
583
584   const Standard_Boolean isFace1Quad = (aType1 == GeomAbs_Cylinder ||
585                                         aType1 == GeomAbs_Cone ||
586                                         aType1 == GeomAbs_Torus);
587
588   const Standard_Boolean isFace2Quad = (aType2 == GeomAbs_Cylinder ||
589                                         aType2 == GeomAbs_Cone ||
590                                         aType2 == GeomAbs_Torus);
591
592   if(aType1==GeomAbs_Plane && aType2==GeomAbs_Plane)
593   {
594     Standard_Real umin, umax, vmin, vmax;
595     BRepTools::UVBounds(myFace1, umin, umax, vmin, vmax);
596     myHS1->ChangeSurface().Load(S1, umin, umax, vmin, vmax);
597     //
598     BRepTools::UVBounds(myFace2, umin, umax, vmin, vmax);
599     myHS2->ChangeSurface().Load(S2, umin, umax, vmin, vmax);
600     Standard_Real TolAng = 1.e-8;
601
602     PerformPlanes(myHS1, myHS2, TolAng, TolTang, myApprox1, myApprox2,
603                                           mySeqOfCurve, myTangentFaces);
604
605     myIsDone = Standard_True;
606     
607     if(!myTangentFaces)
608     {
609       const Standard_Integer NbLinPP = mySeqOfCurve.Length();
610       if(NbLinPP)
611       {
612         Standard_Real aTolFMax;
613         myTolReached3d = 1.e-7;
614         aTolFMax=Max(aTolF1, aTolF2);
615         if (aTolFMax>myTolReached3d)
616         {
617           myTolReached3d=aTolFMax;
618         }
619
620         myTolReached2d = myTolReached3d;
621
622         if (bReverse)
623         {
624           Handle(Geom2d_Curve) aC2D1, aC2D2;
625           const Standard_Integer aNbLin = mySeqOfCurve.Length();
626           for (Standard_Integer i = 1; i <= aNbLin; ++i)
627           {
628             IntTools_Curve& aIC=mySeqOfCurve(i);
629             aC2D1=aIC.FirstCurve2d();
630             aC2D2=aIC.SecondCurve2d();
631             aIC.SetFirstCurve2d(aC2D2);
632             aIC.SetSecondCurve2d(aC2D1);
633           }
634         }
635       }
636     }
637
638     return;
639   }//if(aType1==GeomAbs_Plane && aType2==GeomAbs_Plane){
640
641   if ((aType1==GeomAbs_Plane) && isFace2Quad)
642   {
643     Standard_Real dU, dV;
644
645     // F1
646     Standard_Real umin, umax, vmin, vmax;
647     BRepTools::UVBounds(myFace1, umin, umax, vmin, vmax);
648
649     dU=0.1*(umax-umin);
650     dV=0.1*(vmax-vmin);
651     umin=umin-dU;
652     umax=umax+dU;
653     vmin=vmin-dV;
654     vmax=vmax+dV;
655     myHS1->ChangeSurface().Load(S1, umin, umax, vmin, vmax);
656     // F2
657     BRepTools::UVBounds(myFace2, umin, umax, vmin, vmax);
658     CorrectSurfaceBoundaries(myFace2, (aTolF1 + aTolF2) * 2., umin, umax, vmin, vmax);
659     myHS2->ChangeSurface().Load(S2, umin, umax, vmin, vmax);
660     //
661     if( aType2==GeomAbs_Cone ) { 
662       TolArc = 0.0001; 
663       hasCone = Standard_True; 
664     }
665   }
666   else if ((aType2==GeomAbs_Plane) && isFace1Quad)
667   {
668     Standard_Real dU, dV;
669
670     //F1
671     Standard_Real umin, umax, vmin, vmax;
672     BRepTools::UVBounds(myFace1, umin, umax, vmin, vmax);
673     CorrectSurfaceBoundaries(myFace1, (aTolF1 + aTolF2) * 2., umin, umax, vmin, vmax);
674     myHS1->ChangeSurface().Load(S1, umin, umax, vmin, vmax);
675     // F2
676     BRepTools::UVBounds(myFace2, umin, umax, vmin, vmax);
677     dU=0.1*(umax-umin);
678     dV=0.1*(vmax-vmin);
679     umin=umin-dU;
680     umax=umax+dU;
681     vmin=vmin-dV;
682     vmax=vmax+dV;
683     myHS2->ChangeSurface().Load(S2, umin, umax, vmin, vmax);
684     //
685     if( aType1==GeomAbs_Cone ) {
686       TolArc = 0.0001; 
687       hasCone = Standard_True; 
688     }
689   }
690   else
691   {
692     Standard_Real umin, umax, vmin, vmax;
693     BRepTools::UVBounds(myFace1, umin, umax, vmin, vmax);
694     CorrectSurfaceBoundaries(myFace1, (aTolF1 + aTolF2) * 2., umin, umax, vmin, vmax);
695     myHS1->ChangeSurface().Load(S1, umin, umax, vmin, vmax);
696     BRepTools::UVBounds(myFace2, umin, umax, vmin, vmax);
697     CorrectSurfaceBoundaries(myFace2, (aTolF1 + aTolF2) * 2., umin, umax, vmin, vmax);
698     myHS2->ChangeSurface().Load(S2, umin, umax, vmin, vmax);
699   }
700
701   const Handle(IntTools_TopolTool) dom1 = new IntTools_TopolTool(myHS1);
702   const Handle(IntTools_TopolTool) dom2 = new IntTools_TopolTool(myHS2);
703
704   myLConstruct.Load(dom1, dom2, myHS1, myHS2);
705   
706
707   Tolerances(myHS1, myHS2, TolTang);
708
709   {
710     const Standard_Real UVMaxStep = 0.001;
711     const Standard_Real Deflection = (hasCone) ? 0.085 : 0.1;
712     myIntersector.SetTolerances(TolArc, TolTang, UVMaxStep, Deflection);
713   }
714   
715   if((myHS1->IsUClosed() && !myHS1->IsUPeriodic()) || 
716      (myHS1->IsVClosed() && !myHS1->IsVPeriodic()) ||
717      (myHS2->IsUClosed() && !myHS2->IsUPeriodic()) || 
718      (myHS2->IsVClosed() && !myHS2->IsVPeriodic()))
719   {
720     RestrictLine = Standard_True;
721   }
722   //
723   if((aType1 != GeomAbs_BSplineSurface) &&
724      (aType1 != GeomAbs_BezierSurface) &&
725      (aType1 != GeomAbs_OtherSurface)  &&
726      (aType2 != GeomAbs_BSplineSurface) &&
727      (aType2 != GeomAbs_BezierSurface) &&
728      (aType2 != GeomAbs_OtherSurface))
729   {
730     RestrictLine = Standard_True;
731
732     if ((aType1 == GeomAbs_Torus) ||
733         (aType2 == GeomAbs_Torus))
734     {
735       myListOfPnts.Clear();
736     }
737   }
738
739   //
740   if(!RestrictLine)
741   {
742     TopExp_Explorer aExp;
743     for(Standard_Integer i = 0; (!RestrictLine) && (i < 2); i++)
744     {
745       const TopoDS_Face& aF=(!i) ? myFace1 : myFace2;
746       aExp.Init(aF, TopAbs_EDGE);
747       for(; aExp.More(); aExp.Next())
748       {
749         const TopoDS_Edge& aE=TopoDS::Edge(aExp.Current());
750
751         if(BRep_Tool::Degenerated(aE))
752         {
753           RestrictLine = Standard_True;
754           break;
755         }
756       }
757     }
758   }
759
760   const Standard_Boolean isGeomInt = isTreatAnalityc(aF1, aF2);
761   myIntersector.Perform(myHS1, dom1, myHS2, dom2, TolArc, TolTang, 
762                                   myListOfPnts, RestrictLine, isGeomInt);
763
764   myIsDone = myIntersector.IsDone();
765
766   if (myIsDone)
767   {
768     myTangentFaces=myIntersector.TangentFaces();
769     if (myTangentFaces) {
770       return;
771     }
772     //
773     if(RestrictLine) {
774       myListOfPnts.Clear(); // to use LineConstructor
775     }
776     //
777     const Standard_Integer aNbLin = myIntersector.NbLines();
778     for (Standard_Integer i=1; i <= aNbLin; ++i) {
779       MakeCurve(i, dom1, dom2);
780     }
781     //
782     ComputeTolReached3d();
783     //
784     if (bReverse) {
785       Handle(Geom2d_Curve) aC2D1, aC2D2;
786       //
787       const Standard_Integer aNbLin=mySeqOfCurve.Length();
788       for (Standard_Integer i=1; i<=aNbLin; ++i)
789       {
790         IntTools_Curve& aIC=mySeqOfCurve(i);
791         aC2D1=aIC.FirstCurve2d();
792         aC2D2=aIC.SecondCurve2d();
793         aIC.SetFirstCurve2d(aC2D2);
794         aIC.SetSecondCurve2d(aC2D1);
795       }
796     }
797
798     // Points
799     Standard_Real U1,V1,U2,V2;
800     IntTools_PntOnFace aPntOnF1, aPntOnF2;
801     IntTools_PntOn2Faces aPntOn2Faces;
802     //
803     const Standard_Integer aNbPnts = myIntersector.NbPnts();
804     for (Standard_Integer i=1; i <= aNbPnts; ++i)
805     {
806       const IntSurf_PntOn2S& aISPnt=myIntersector.Point(i).PntOn2S();
807       const gp_Pnt& aPnt=aISPnt.Value();
808       aISPnt.Parameters(U1,V1,U2,V2);
809       aPntOnF1.Init(myFace1, aPnt, U1, V1);
810       aPntOnF2.Init(myFace2, aPnt, U2, V2);
811       //
812       if (!bReverse)
813       {
814         aPntOn2Faces.SetP1(aPntOnF1);
815         aPntOn2Faces.SetP2(aPntOnF2);
816       }
817       else
818       {
819         aPntOn2Faces.SetP2(aPntOnF1);
820         aPntOn2Faces.SetP1(aPntOnF2);
821       }
822
823       myPnts.Append(aPntOn2Faces);
824     }
825   }
826 }
827
828 //=======================================================================
829 //function :ComputeTolReached3d 
830 //purpose  : 
831 //=======================================================================
832   void IntTools_FaceFace::ComputeTolReached3d()
833 {
834   Standard_Boolean bCase1;
835   Standard_Integer aNbLin, i;
836   GeomAbs_SurfaceType aType1, aType2;
837   //
838   aNbLin=myIntersector.NbLines();
839   if (!aNbLin) {
840     return;
841   }
842   //
843   aType1=myHS1->Surface().GetType();
844   aType2=myHS2->Surface().GetType();
845   //
846   bCase1=((aType1==GeomAbs_Plane && aType2==GeomAbs_SurfaceOfExtrusion) ||
847           (aType2==GeomAbs_Plane && aType1==GeomAbs_SurfaceOfExtrusion));
848   //
849   if (aType1==GeomAbs_Cylinder && aType2==GeomAbs_Cylinder) {
850     if (aNbLin==2){ 
851       Handle(IntPatch_Line) aIL1, aIL2;
852       IntPatch_IType aTL1, aTL2;
853       //
854       aIL1=myIntersector.Line(1);
855       aIL2=myIntersector.Line(2);
856       aTL1=aIL1->ArcType();
857       aTL2=aIL2->ArcType();
858       if (aTL1==IntPatch_Lin && aTL2==IntPatch_Lin) {
859         Standard_Real aD, aDTresh, dTol;
860         gp_Lin aL1, aL2;
861         //
862         dTol=1.e-8;
863         aDTresh=1.5e-6;
864         //
865         aL1=Handle(IntPatch_GLine)::DownCast(aIL1)->Line();
866         aL2=Handle(IntPatch_GLine)::DownCast(aIL2)->Line();
867         aD=aL1.Distance(aL2);
868         aD=0.5*aD;
869         if (aD<aDTresh) {
870           myTolReached3d=aD+dTol;
871         }
872         return;
873       }
874     }
875     //ZZ
876     if (aNbLin) {// Check the distances
877       Standard_Integer  aNbP, j ;
878       Standard_Real aT1, aT2, dT, aD2, aD2Max, aEps, aT11, aT12;
879       //
880       aD2Max=0.;
881       aNbP=10;
882       aNbLin=mySeqOfCurve.Length();
883       //
884       for (i=1; i<=aNbLin; ++i) {
885         const IntTools_Curve& aIC=mySeqOfCurve(i);
886         const Handle(Geom_Curve)& aC3D=aIC.Curve();
887         const Handle(Geom2d_Curve)& aC2D1=aIC.FirstCurve2d();
888         const Handle(Geom2d_Curve)& aC2D2=aIC.SecondCurve2d();
889         //
890         if (aC3D.IsNull()) {
891           continue;
892         }
893         const Handle(Geom_BSplineCurve)& aBC=
894           Handle(Geom_BSplineCurve)::DownCast(aC3D);
895         if (aBC.IsNull()) {
896           continue;
897         }
898         //
899         aT1=aBC->FirstParameter();
900         aT2=aBC->LastParameter();
901         //
902         aEps=0.01*(aT2-aT1);
903         dT=(aT2-aT1)/aNbP;
904         for (j=1; j<aNbP; ++j) {
905           aT11=aT1+j*dT;
906           aT12=aT11+dT;
907           aD2=FindMaxSquareDistance(aT11, aT12, aEps, aC3D, aC2D1, aC2D2,
908                                     myHS1, myHS2, myFace1, myFace2, myContext);
909           if (aD2>aD2Max) {
910             aD2Max=aD2;
911           }
912         }
913       }//for (i=1; i<=aNbLin; ++i) {
914       //
915       myTolReached3d=sqrt(aD2Max);
916     }// if (aNbLin) 
917   }// if (aType1==GeomAbs_Cylinder && aType2==GeomAbs_Cylinder) {
918   //
919   //904/G3 f
920   else if (aType1==GeomAbs_Plane && aType2==GeomAbs_Plane) {
921     Standard_Real aTolF1, aTolF2, aTolFMax, aTolTresh;
922     //
923     aTolTresh=1.e-7;
924     //
925     aTolF1 = BRep_Tool::Tolerance(myFace1);
926     aTolF2 = BRep_Tool::Tolerance(myFace2);
927     aTolFMax=Max(aTolF1, aTolF2);
928     //
929     if (aTolFMax>aTolTresh) {
930       myTolReached3d=aTolFMax;
931     }
932   }//if (aType1==GeomAbs_Plane && aType2==GeomAbs_Plane) {
933   //t
934   //IFV Bug OCC20297 
935   else if((aType1 == GeomAbs_Cylinder && aType2 == GeomAbs_Plane) ||
936           (aType2 == GeomAbs_Cylinder && aType1 == GeomAbs_Plane)) {
937     if(aNbLin == 1) {
938       const Handle(IntPatch_Line)& aIL1 = myIntersector.Line(1);
939       if(aIL1->ArcType() == IntPatch_Circle) {
940         gp_Circ aCir = Handle(IntPatch_GLine)::DownCast(aIL1)->Circle();
941         gp_XYZ aCirDir = aCir.Axis().Direction().XYZ();
942         gp_XYZ aPlDir;
943         gp_Pln aPln;
944         if(aType1 == GeomAbs_Plane) {
945           aPln = myHS1->Surface().Plane();
946         }
947         else {
948           aPln = myHS2->Surface().Plane();
949         }
950         aPlDir = aPln.Axis().Direction().XYZ();
951         Standard_Real cs = aCirDir*aPlDir;
952         if(cs < 0.) aPlDir.Reverse();
953         Standard_Real eps = 1.e-14;
954         if(!aPlDir.IsEqual(aCirDir, eps)) {
955           Standard_Integer aNbP = 11;
956           Standard_Real dt = 2.*M_PI / (aNbP - 1), t;
957           for(t = 0.; t < 2.*M_PI; t += dt) {
958             Standard_Real d = aPln.Distance(ElCLib::Value(t, aCir)); 
959             if(myTolReached3d < d) myTolReached3d = d;
960           }
961           myTolReached3d *= 1.1;
962         }
963       } //aIL1->ArcType() == IntPatch_Circle
964     } //aNbLin == 1
965   } // aType1 == GeomAbs_Cylinder && aType2 == GeomAbs_Plane) 
966   //End IFV Bug OCC20297
967   //
968   else if ((aType1==GeomAbs_Plane && aType2==GeomAbs_Torus) ||
969            (aType2==GeomAbs_Plane && aType1==GeomAbs_Torus)) {
970     aNbLin=mySeqOfCurve.Length();
971     if (aNbLin!=1) {
972       return;
973     }
974     //
975     Standard_Integer aNbP;
976     Standard_Real aT, aT1, aT2, dT, aUT, aVT, aUP, aVP;
977     Standard_Real aDP, aDT, aDmax;
978     gp_Pln aPln;
979     gp_Torus aTorus;
980     gp_Pnt aP, aPP, aPT;
981     //
982     const IntTools_Curve& aIC=mySeqOfCurve(1);
983     const Handle(Geom_Curve)& aC3D=aIC.Curve();
984     const Handle(Geom_BSplineCurve)& aBS=
985       Handle(Geom_BSplineCurve)::DownCast(aC3D);
986     if (aBS.IsNull()) {
987       return;
988     }
989     //
990     aT1=aBS->FirstParameter();
991     aT2=aBS->LastParameter();
992     //
993     aPln  =(aType1==GeomAbs_Plane) ? myHS1->Plane() : myHS2->Plane();
994     aTorus=(aType1==GeomAbs_Plane) ? myHS2->Torus() : myHS1->Torus();
995     //
996     aDmax=-1.;
997     aNbP=11;
998     dT=(aT2-aT1)/(aNbP-1);
999     for (i=0; i<aNbP; ++i) {
1000       aT=aT1+i*dT;
1001       if (i==aNbP-1) {
1002         aT=aT2;
1003       }
1004       //
1005       aC3D->D0(aT, aP);
1006       //
1007       ElSLib::Parameters(aPln, aP, aUP, aVP);
1008       aPP=ElSLib::Value(aUP, aVP, aPln);
1009       aDP=aP.SquareDistance(aPP);
1010       if (aDP>aDmax) {
1011         aDmax=aDP;
1012       }
1013       //
1014       ElSLib::Parameters(aTorus, aP, aUT, aVT);
1015       aPT=ElSLib::Value(aUT, aVT, aTorus);
1016       aDT=aP.SquareDistance(aPT);
1017       if (aDT>aDmax) {
1018         aDmax=aDT;
1019       }
1020     }
1021     //
1022     if (aDmax > myTolReached3d*myTolReached3d) {
1023       myTolReached3d=sqrt(aDmax);
1024       myTolReached3d=1.1*myTolReached3d;
1025     }
1026   }// if ((aType1==GeomAbs_Plane && aType2==GeomAbs_Torus) ||
1027   //
1028   else if ((aType1==GeomAbs_SurfaceOfRevolution && aType2==GeomAbs_Cylinder) ||
1029            (aType2==GeomAbs_SurfaceOfRevolution && aType1==GeomAbs_Cylinder)) {
1030     Standard_Integer j, aNbP;
1031     Standard_Real aT, aT1, aT2, dT, aD2max, aD2;
1032     //
1033     aNbLin=mySeqOfCurve.Length();
1034     aD2max=0.;
1035     aNbP=11;
1036     //
1037     for (i=1; i<=aNbLin; ++i) {
1038       const IntTools_Curve& aIC=mySeqOfCurve(i);
1039       const Handle(Geom_Curve)& aC3D=aIC.Curve();
1040       const Handle(Geom2d_Curve)& aC2D1=aIC.FirstCurve2d();
1041       const Handle(Geom2d_Curve)& aC2D2=aIC.SecondCurve2d();
1042       //
1043       if (aC3D.IsNull()) {
1044         continue;
1045       }
1046       const Handle(Geom_BSplineCurve)& aBC=
1047         Handle(Geom_BSplineCurve)::DownCast(aC3D);
1048       if (aBC.IsNull()) {
1049         return;
1050       }
1051       //
1052       aT1=aBC->FirstParameter();
1053       aT2=aBC->LastParameter();
1054       //
1055       dT=(aT2-aT1)/(aNbP-1);
1056       for (j=0; j<aNbP; ++j) {
1057         aT=aT1+j*dT;
1058         if (j==aNbP-1) {
1059           aT=aT2;
1060         }
1061         //
1062         aD2=MaxSquareDistance(aT, aC3D, aC2D1, aC2D2,
1063                               myHS1, myHS2, myFace1, myFace2, myContext);
1064         if (aD2>aD2max) {
1065           aD2max=aD2;
1066         }
1067       }//for (j=0; j<aNbP; ++j) {
1068      
1069     }//for (i=1; i<=aNbLin; ++i) {
1070     //
1071     aD2=myTolReached3d*myTolReached3d;
1072     if (aD2max > aD2) {
1073       myTolReached3d=sqrt(aD2max);
1074     }
1075   }//if((aType1==GeomAbs_SurfaceOfRevolution ...
1076   else if ((aType1==GeomAbs_Plane && aType2==GeomAbs_Sphere) ||
1077            (aType2==GeomAbs_Plane && aType1==GeomAbs_Sphere)) {
1078     Standard_Integer  j, aNbP;
1079     Standard_Real aT1, aT2, dT, aD2max, aD2, aEps, aT11, aT12;
1080     //
1081     aNbLin=mySeqOfCurve.Length();
1082     aD2max=0.;
1083     aNbP=10;
1084     //
1085     for (i=1; i<=aNbLin; ++i) {
1086       const IntTools_Curve& aIC=mySeqOfCurve(i);
1087       const Handle(Geom_Curve)& aC3D=aIC.Curve();
1088       const Handle(Geom2d_Curve)& aC2D1=aIC.FirstCurve2d();
1089       const Handle(Geom2d_Curve)& aC2D2=aIC.SecondCurve2d();
1090       //
1091       const Handle(Geom2d_BSplineCurve)& aBC2D1=
1092         Handle(Geom2d_BSplineCurve)::DownCast(aC2D1);
1093       const Handle(Geom2d_BSplineCurve)& aBC2D2=
1094         Handle(Geom2d_BSplineCurve)::DownCast(aC2D2);
1095       //
1096       if (aBC2D1.IsNull() && aBC2D2.IsNull()) {
1097         return;
1098       }
1099       //
1100       if (!aBC2D1.IsNull()) {
1101         aT1=aBC2D1->FirstParameter();
1102         aT2=aBC2D1->LastParameter();
1103       }
1104       else {
1105         aT1=aBC2D2->FirstParameter();
1106         aT2=aBC2D2->LastParameter();
1107       }
1108       //
1109       aEps=0.01*(aT2-aT1);
1110       dT=(aT2-aT1)/aNbP;
1111       for (j=0; j<aNbP; ++j) {
1112         aT11=aT1+j*dT;
1113         aT12=aT11+dT;
1114         if (j==aNbP-1) {
1115           aT12=aT2;
1116         }
1117         //
1118         aD2=FindMaxSquareDistance(aT11, aT12, aEps, aC3D, aC2D1, aC2D2,
1119                                   myHS1, myHS2, myFace1, myFace2, myContext);
1120         if (aD2>aD2max) {
1121           aD2max=aD2;
1122         }
1123       }//for (j=0; j<aNbP; ++j) {
1124      
1125     }//for (i=1; i<=aNbLin; ++i) {
1126     //
1127     aD2=myTolReached3d*myTolReached3d;
1128     if (aD2max > aD2) {
1129       myTolReached3d=sqrt(aD2max);
1130     }
1131   }//else if ((aType1==GeomAbs_Plane && aType2==GeomAbs_Sphere) ...
1132   else if (!myApprox || bCase1) {
1133   //else if (!myApprox) {
1134     Standard_Integer aNbP, j;
1135     Standard_Real aT1, aT2, dT, aD2, aD2Max, aEps, aT11, aT12;
1136     //
1137     aD2Max=0.;
1138     aNbLin=mySeqOfCurve.Length();
1139     //
1140     for (i=1; i<=aNbLin; ++i) {
1141       const IntTools_Curve& aIC=mySeqOfCurve(i);
1142       const Handle(Geom_Curve)& aC3D=aIC.Curve();
1143       const Handle(Geom2d_Curve)& aC2D1=aIC.FirstCurve2d();
1144       const Handle(Geom2d_Curve)& aC2D2=aIC.SecondCurve2d();
1145       //
1146       if (aC3D.IsNull()) {
1147         continue;
1148       }
1149       const Handle(Geom_BSplineCurve)& aBC=
1150         Handle(Geom_BSplineCurve)::DownCast(aC3D);
1151       if (aBC.IsNull()) {
1152         continue;
1153       }
1154       //
1155       aT1=aBC->FirstParameter();
1156       aT2=aBC->LastParameter();
1157       //
1158       aEps=0.0001*(aT2-aT1);
1159       aNbP=11;
1160       dT=(aT2-aT1)/aNbP;
1161       for (j=1; j<aNbP-1; ++j) {
1162         aT11=aT1+j*dT;
1163         aT12=aT11+dT;
1164         aD2=FindMaxSquareDistance(aT11, aT12, aEps, aC3D, aC2D1, aC2D2,
1165                                   myHS1, myHS2, myFace1, myFace2, myContext);
1166         if (aD2>aD2Max) {
1167           aD2Max=aD2;
1168         }
1169       }
1170     }//for (i=1; i<=aNbLin; ++i) {
1171     myTolReached3d=sqrt(aD2Max);
1172   }
1173 }
1174 //=======================================================================
1175 //function : MakeCurve
1176 //purpose  : 
1177 //=======================================================================
1178   void IntTools_FaceFace::MakeCurve(const Standard_Integer Index,
1179                                     const Handle(Adaptor3d_TopolTool)& dom1,
1180                                     const Handle(Adaptor3d_TopolTool)& dom2) 
1181 {
1182   Standard_Boolean bDone, rejectSurface, reApprox, bAvoidLineConstructor,
1183                    bPCurvesOk;
1184   Standard_Boolean ok;
1185   Standard_Integer i, j, aNbParts;
1186   Standard_Real fprm, lprm;
1187   Standard_Real Tolpc;
1188   Handle(IntPatch_Line) L;
1189   IntPatch_IType typl;
1190   Handle(Geom_Curve) newc;
1191   //
1192   const Standard_Real TOLCHECK   =0.0000001;
1193   const Standard_Real TOLANGCHECK=0.1;
1194   //
1195   rejectSurface = Standard_False;
1196   reApprox = Standard_False;
1197   //
1198   bPCurvesOk = Standard_True;
1199  
1200  reapprox:;
1201   
1202   Tolpc = myTolApprox;
1203   bAvoidLineConstructor = Standard_False;
1204   L = myIntersector.Line(Index);
1205   typl = L->ArcType();
1206   //
1207   if(typl==IntPatch_Walking) {
1208     Handle(IntPatch_Line) anewL;
1209     //
1210     const Handle(IntPatch_WLine)& aWLine=
1211       Handle(IntPatch_WLine)::DownCast(L);
1212 #ifdef DEB_DUMPWLINE
1213     DumpWLine(aWLine);
1214 #endif
1215     anewL = ComputePurgedWLine(aWLine);
1216     if(anewL.IsNull()) {
1217       return;
1218     }
1219     L = anewL;
1220 #ifdef DEB_DUMPWLINE
1221       const Handle(IntPatch_WLine)& aWLineX = Handle(IntPatch_WLine)::DownCast(L);
1222       DumpWLine(aWLineX);
1223 #endif
1224     //
1225     if(!myListOfPnts.IsEmpty()) {
1226       bAvoidLineConstructor = Standard_True;
1227     }
1228
1229     Standard_Integer nbp = aWLine->NbPnts();
1230     const IntSurf_PntOn2S& p1 = aWLine->Point(1);
1231     const IntSurf_PntOn2S& p2 = aWLine->Point(nbp);
1232
1233     const gp_Pnt& P1 = p1.Value();
1234     const gp_Pnt& P2 = p2.Value();
1235
1236     if(P1.SquareDistance(P2) < 1.e-14) {
1237       bAvoidLineConstructor = Standard_False;
1238     }
1239
1240   }
1241   //
1242   // Line Constructor
1243   if(!bAvoidLineConstructor) {
1244     myLConstruct.Perform(L);
1245     //
1246     bDone=myLConstruct.IsDone();
1247     aNbParts=myLConstruct.NbParts();
1248     if (!bDone|| !aNbParts) {
1249       return;
1250     }
1251   }
1252   // Do the Curve
1253   
1254   
1255   typl=L->ArcType();
1256   switch (typl) {
1257   //########################################  
1258   // Line, Parabola, Hyperbola
1259   //########################################  
1260   case IntPatch_Lin:
1261   case IntPatch_Parabola: 
1262   case IntPatch_Hyperbola: {
1263     if (typl == IntPatch_Lin) {
1264       newc = 
1265         new Geom_Line (Handle(IntPatch_GLine)::DownCast(L)->Line());
1266     }
1267
1268     else if (typl == IntPatch_Parabola) {
1269       newc = 
1270         new Geom_Parabola(Handle(IntPatch_GLine)::DownCast(L)->Parabola());
1271     }
1272     
1273     else if (typl == IntPatch_Hyperbola) {
1274       newc = 
1275         new Geom_Hyperbola (Handle(IntPatch_GLine)::DownCast(L)->Hyperbola());
1276     }
1277     //
1278     // myTolReached3d
1279     if (typl == IntPatch_Lin) {
1280       TolR3d (myFace1, myFace2, myTolReached3d);
1281     }
1282     //
1283     aNbParts=myLConstruct.NbParts();
1284     for (i=1; i<=aNbParts; i++) {
1285       myLConstruct.Part(i, fprm, lprm);
1286       
1287       if (!Precision::IsNegativeInfinite(fprm) && 
1288           !Precision::IsPositiveInfinite(lprm)) {
1289         //
1290         IntTools_Curve aCurve;
1291         //
1292         Handle(Geom_TrimmedCurve) aCT3D=new Geom_TrimmedCurve(newc, fprm, lprm);
1293         aCurve.SetCurve(aCT3D);
1294         if (typl == IntPatch_Parabola) {
1295           Standard_Real aTolF1, aTolF2, aTolBase;
1296           
1297           aTolF1 = BRep_Tool::Tolerance(myFace1);
1298           aTolF2 = BRep_Tool::Tolerance(myFace2);
1299           aTolBase=aTolF1+aTolF2;
1300           myTolReached3d=IntTools_Tools::CurveTolerance(aCT3D, aTolBase);
1301         }
1302         //
1303         aCurve.SetCurve(new Geom_TrimmedCurve(newc, fprm, lprm));
1304         if(myApprox1) { 
1305           Handle (Geom2d_Curve) C2d;
1306           BuildPCurves(fprm, lprm, Tolpc, myHS1->ChangeSurface().Surface(), newc, C2d);
1307           if(Tolpc>myTolReached2d || myTolReached2d==0.) { 
1308             myTolReached2d=Tolpc;
1309           }
1310             //     
1311             aCurve.SetFirstCurve2d(new Geom2d_TrimmedCurve(C2d,fprm,lprm));
1312           }
1313           else { 
1314             Handle(Geom2d_BSplineCurve) H1;
1315             //
1316             aCurve.SetFirstCurve2d(H1);
1317           }
1318         
1319         if(myApprox2) { 
1320           Handle (Geom2d_Curve) C2d;
1321           BuildPCurves(fprm,lprm,Tolpc,myHS2->ChangeSurface().Surface(),newc,C2d);
1322           if(Tolpc>myTolReached2d || myTolReached2d==0.) { 
1323             myTolReached2d=Tolpc;
1324           }
1325           //
1326           aCurve.SetSecondCurve2d(new Geom2d_TrimmedCurve(C2d,fprm,lprm));
1327           }
1328         else { 
1329           Handle(Geom2d_BSplineCurve) H1;
1330           //
1331           aCurve.SetSecondCurve2d(H1);
1332         }
1333         mySeqOfCurve.Append(aCurve);
1334       } // end of if (!Precision::IsNegativeInfinite(fprm) &&  !Precision::IsPositiveInfinite(lprm))
1335
1336       else {
1337         //  on regarde si on garde
1338         //
1339         Standard_Boolean bFNIt, bLPIt;
1340         Standard_Real aTestPrm, dT=100.;
1341
1342         bFNIt=Precision::IsNegativeInfinite(fprm);
1343         bLPIt=Precision::IsPositiveInfinite(lprm);
1344         
1345         aTestPrm=0.;
1346         
1347         if (bFNIt && !bLPIt) {
1348           aTestPrm=lprm-dT;
1349         }
1350         else if (!bFNIt && bLPIt) {
1351           aTestPrm=fprm+dT;
1352         }
1353         
1354         gp_Pnt ptref(newc->Value(aTestPrm));
1355         //
1356
1357         Standard_Real u1, v1, u2, v2, Tol;
1358         
1359         Tol = Precision::Confusion();
1360         Parameters(myHS1, myHS2, ptref,  u1, v1, u2, v2);
1361         ok = (dom1->Classify(gp_Pnt2d(u1, v1), Tol) != TopAbs_OUT);
1362         if(ok) { 
1363           ok = (dom2->Classify(gp_Pnt2d(u2,v2),Tol) != TopAbs_OUT); 
1364         }
1365         if (ok) {
1366           Handle(Geom2d_BSplineCurve) H1;
1367           mySeqOfCurve.Append(IntTools_Curve(newc, H1, H1));
1368         }
1369       }
1370     }// end of for (i=1; i<=myLConstruct.NbParts(); i++)
1371   }// case IntPatch_Lin:  case IntPatch_Parabola:  case IntPatch_Hyperbola:
1372     break;
1373
1374   //########################################  
1375   // Circle and Ellipse
1376   //########################################  
1377   case IntPatch_Circle: 
1378   case IntPatch_Ellipse: {
1379
1380     if (typl == IntPatch_Circle) {
1381       newc = new Geom_Circle
1382         (Handle(IntPatch_GLine)::DownCast(L)->Circle());
1383     }
1384     else { //IntPatch_Ellipse
1385       newc = new Geom_Ellipse
1386         (Handle(IntPatch_GLine)::DownCast(L)->Ellipse());
1387     }
1388     //
1389     // myTolReached3d
1390     TolR3d (myFace1, myFace2, myTolReached3d);
1391     //
1392     aNbParts=myLConstruct.NbParts();
1393     //
1394     Standard_Real aPeriod, aNul;
1395     TColStd_SequenceOfReal aSeqFprm,  aSeqLprm;
1396     
1397     aNul=0.;
1398     aPeriod=M_PI+M_PI;
1399
1400     for (i=1; i<=aNbParts; i++) {
1401       myLConstruct.Part(i, fprm, lprm);
1402
1403       if (fprm < aNul && lprm > aNul) {
1404         // interval that goes through 0. is divided on two intervals;
1405         while (fprm<aNul || fprm>aPeriod) fprm=fprm+aPeriod;
1406         while (lprm<aNul || lprm>aPeriod) lprm=lprm+aPeriod;
1407         //
1408         if((aPeriod - fprm) > Tolpc) {
1409           aSeqFprm.Append(fprm);
1410           aSeqLprm.Append(aPeriod);
1411         }
1412         else {
1413           gp_Pnt P1 = newc->Value(fprm);
1414           gp_Pnt P2 = newc->Value(aPeriod);
1415           Standard_Real aTolDist = BRep_Tool::Tolerance(myFace1) + BRep_Tool::Tolerance(myFace2);
1416           aTolDist = (myTolReached3d > aTolDist) ? myTolReached3d : aTolDist;
1417
1418           if(P1.Distance(P2) > aTolDist) {
1419             Standard_Real anewpar = fprm;
1420
1421             if(ParameterOutOfBoundary(fprm, newc, myFace1, myFace2, lprm, Standard_False, anewpar, myContext)) {
1422               fprm = anewpar;
1423             }
1424             aSeqFprm.Append(fprm);
1425             aSeqLprm.Append(aPeriod);
1426           }
1427         }
1428
1429         //
1430         if((lprm - aNul) > Tolpc) {
1431           aSeqFprm.Append(aNul);
1432           aSeqLprm.Append(lprm);
1433         }
1434         else {
1435           gp_Pnt P1 = newc->Value(aNul);
1436           gp_Pnt P2 = newc->Value(lprm);
1437           Standard_Real aTolDist = BRep_Tool::Tolerance(myFace1) + BRep_Tool::Tolerance(myFace2);
1438           aTolDist = (myTolReached3d > aTolDist) ? myTolReached3d : aTolDist;
1439
1440           if(P1.Distance(P2) > aTolDist) {
1441             Standard_Real anewpar = lprm;
1442
1443             if(ParameterOutOfBoundary(lprm, newc, myFace1, myFace2, fprm, Standard_True, anewpar, myContext)) {
1444               lprm = anewpar;
1445             }
1446             aSeqFprm.Append(aNul);
1447             aSeqLprm.Append(lprm);
1448           }
1449         }
1450       }
1451       else {
1452         // usual interval 
1453         aSeqFprm.Append(fprm);
1454         aSeqLprm.Append(lprm);
1455       }
1456     }
1457     
1458     //
1459     aNbParts=aSeqFprm.Length();
1460     for (i=1; i<=aNbParts; i++) {
1461       fprm=aSeqFprm(i);
1462       lprm=aSeqLprm(i);
1463       //
1464       Standard_Real aRealEpsilon=RealEpsilon();
1465       if (Abs(fprm) > aRealEpsilon || Abs(lprm-2.*M_PI) > aRealEpsilon) {
1466         //==============================================
1467         ////
1468         IntTools_Curve aCurve;
1469         Handle(Geom_TrimmedCurve) aTC3D=new Geom_TrimmedCurve(newc,fprm,lprm);
1470         aCurve.SetCurve(aTC3D);
1471         fprm=aTC3D->FirstParameter();
1472         lprm=aTC3D->LastParameter ();
1473         ////    
1474         if (typl == IntPatch_Circle || typl == IntPatch_Ellipse) {//// 
1475           if(myApprox1) { 
1476             Handle (Geom2d_Curve) C2d;
1477             BuildPCurves(fprm,lprm,Tolpc,myHS1->ChangeSurface().Surface(),newc,C2d);
1478             if(Tolpc>myTolReached2d || myTolReached2d==0) { 
1479               myTolReached2d=Tolpc;
1480             }
1481             //
1482             aCurve.SetFirstCurve2d(C2d);
1483           }
1484           else { //// 
1485             Handle(Geom2d_BSplineCurve) H1;
1486             aCurve.SetFirstCurve2d(H1);
1487           }
1488
1489
1490           if(myApprox2) { 
1491             Handle (Geom2d_Curve) C2d;
1492             BuildPCurves(fprm,lprm,Tolpc,myHS2->ChangeSurface().Surface(),newc,C2d);
1493             if(Tolpc>myTolReached2d || myTolReached2d==0) { 
1494               myTolReached2d=Tolpc;
1495             }
1496             //
1497             aCurve.SetSecondCurve2d(C2d);
1498           }
1499           else { 
1500             Handle(Geom2d_BSplineCurve) H1;
1501             aCurve.SetSecondCurve2d(H1);
1502           }
1503         }
1504         
1505         else { 
1506           Handle(Geom2d_BSplineCurve) H1;
1507           aCurve.SetFirstCurve2d(H1);
1508           aCurve.SetSecondCurve2d(H1);
1509         }
1510         mySeqOfCurve.Append(aCurve);
1511           //==============================================      
1512       } //if (Abs(fprm) > RealEpsilon() || Abs(lprm-2.*M_PI) > RealEpsilon())
1513
1514       else {
1515         //  on regarde si on garde
1516         //
1517         if (aNbParts==1) {
1518 //        if (Abs(fprm) < RealEpsilon() &&  Abs(lprm-2.*M_PI) < RealEpsilon()) {
1519           if (Abs(fprm) <= aRealEpsilon && Abs(lprm-2.*M_PI) <= aRealEpsilon) {
1520             IntTools_Curve aCurve;
1521             Handle(Geom_TrimmedCurve) aTC3D=new Geom_TrimmedCurve(newc,fprm,lprm);
1522             aCurve.SetCurve(aTC3D);
1523             fprm=aTC3D->FirstParameter();
1524             lprm=aTC3D->LastParameter ();
1525             
1526             if(myApprox1) { 
1527               Handle (Geom2d_Curve) C2d;
1528               BuildPCurves(fprm,lprm,Tolpc,myHS1->ChangeSurface().Surface(),newc,C2d);
1529               if(Tolpc>myTolReached2d || myTolReached2d==0) { 
1530                 myTolReached2d=Tolpc;
1531               }
1532               //
1533               aCurve.SetFirstCurve2d(C2d);
1534             }
1535             else { //// 
1536               Handle(Geom2d_BSplineCurve) H1;
1537               aCurve.SetFirstCurve2d(H1);
1538             }
1539
1540             if(myApprox2) { 
1541               Handle (Geom2d_Curve) C2d;
1542               BuildPCurves(fprm,lprm,Tolpc,myHS2->ChangeSurface().Surface(),newc,C2d);
1543               if(Tolpc>myTolReached2d || myTolReached2d==0) { 
1544                 myTolReached2d=Tolpc;
1545               }
1546               //
1547               aCurve.SetSecondCurve2d(C2d);
1548             }
1549             else { 
1550               Handle(Geom2d_BSplineCurve) H1;
1551               aCurve.SetSecondCurve2d(H1);
1552             }
1553             mySeqOfCurve.Append(aCurve);
1554             break;
1555           }
1556         }
1557         //
1558         Standard_Real aTwoPIdiv17, u1, v1, u2, v2, Tol;
1559
1560         aTwoPIdiv17=2.*M_PI/17.;
1561
1562         for (j=0; j<=17; j++) {
1563           gp_Pnt ptref (newc->Value (j*aTwoPIdiv17));
1564           Tol = Precision::Confusion();
1565
1566           Parameters(myHS1, myHS2, ptref, u1, v1, u2, v2);
1567           ok = (dom1->Classify(gp_Pnt2d(u1,v1),Tol) != TopAbs_OUT);
1568           if(ok) { 
1569             ok = (dom2->Classify(gp_Pnt2d(u2,v2),Tol) != TopAbs_OUT);
1570           }
1571           if (ok) {
1572             IntTools_Curve aCurve;
1573             aCurve.SetCurve(newc);
1574             //==============================================
1575             if (typl == IntPatch_Circle || typl == IntPatch_Ellipse) {
1576               
1577               if(myApprox1) { 
1578                 Handle (Geom2d_Curve) C2d;
1579                 BuildPCurves(fprm, lprm, Tolpc, myHS1->ChangeSurface().Surface(), newc, C2d);
1580                 if(Tolpc>myTolReached2d || myTolReached2d==0) { 
1581                   myTolReached2d=Tolpc;
1582                 }
1583                 // 
1584                 aCurve.SetFirstCurve2d(C2d);
1585               }
1586               else { 
1587                 Handle(Geom2d_BSplineCurve) H1;
1588                 aCurve.SetFirstCurve2d(H1);
1589               }
1590                 
1591               if(myApprox2) { 
1592                 Handle (Geom2d_Curve) C2d;
1593                 BuildPCurves(fprm, lprm, Tolpc,myHS2->ChangeSurface().Surface(), newc, C2d);
1594                 if(Tolpc>myTolReached2d || myTolReached2d==0) { 
1595                   myTolReached2d=Tolpc;
1596                 }
1597                 //              
1598                 aCurve.SetSecondCurve2d(C2d);
1599               }
1600                 
1601               else { 
1602                 Handle(Geom2d_BSplineCurve) H1;
1603                 aCurve.SetSecondCurve2d(H1);
1604               }
1605             }//  end of if (typl == IntPatch_Circle || typl == IntPatch_Ellipse)
1606              
1607             else { 
1608               Handle(Geom2d_BSplineCurve) H1;
1609               //        
1610               aCurve.SetFirstCurve2d(H1);
1611               aCurve.SetSecondCurve2d(H1);
1612             }
1613             //==============================================    
1614             //
1615             mySeqOfCurve.Append(aCurve);
1616             break;
1617
1618             }//  end of if (ok) {
1619           }//  end of for (Standard_Integer j=0; j<=17; j++)
1620         }//  end of else { on regarde si on garde
1621       }// for (i=1; i<=myLConstruct.NbParts(); i++)
1622     }// IntPatch_Circle: IntPatch_Ellipse:
1623     break;
1624     
1625   case IntPatch_Analytic: {
1626     IntSurf_Quadric quad1,quad2;
1627     GeomAbs_SurfaceType typs = myHS1->Surface().GetType();
1628     
1629     switch (typs) {
1630       case GeomAbs_Plane:
1631         quad1.SetValue(myHS1->Surface().Plane());
1632         break;
1633       case GeomAbs_Cylinder:
1634         quad1.SetValue(myHS1->Surface().Cylinder());
1635         break;
1636       case GeomAbs_Cone:
1637         quad1.SetValue(myHS1->Surface().Cone());
1638         break;
1639       case GeomAbs_Sphere:
1640         quad1.SetValue(myHS1->Surface().Sphere());
1641         break;
1642       default:
1643         Standard_ConstructionError::Raise("GeomInt_IntSS::MakeCurve 1");
1644       }
1645       
1646     typs = myHS2->Surface().GetType();
1647     
1648     switch (typs) {
1649       case GeomAbs_Plane:
1650         quad2.SetValue(myHS2->Surface().Plane());
1651         break;
1652       case GeomAbs_Cylinder:
1653         quad2.SetValue(myHS2->Surface().Cylinder());
1654         break;
1655       case GeomAbs_Cone:
1656         quad2.SetValue(myHS2->Surface().Cone());
1657         break;
1658       case GeomAbs_Sphere:
1659         quad2.SetValue(myHS2->Surface().Sphere());
1660         break;
1661       default:
1662         Standard_ConstructionError::Raise("GeomInt_IntSS::MakeCurve 2");
1663       }
1664     //
1665     //=========
1666     IntPatch_ALineToWLine convert (quad1, quad2);
1667       
1668     if (!myApprox) {
1669       aNbParts=myLConstruct.NbParts();
1670       for (i=1; i<=aNbParts; i++) {
1671         myLConstruct.Part(i, fprm, lprm);
1672         Handle(IntPatch_WLine) WL = 
1673           convert.MakeWLine(Handle(IntPatch_ALine)::DownCast(L), fprm, lprm);
1674         //
1675         Handle(Geom2d_BSplineCurve) H1;
1676         Handle(Geom2d_BSplineCurve) H2;
1677
1678         if(myApprox1) {
1679           H1 = MakeBSpline2d(WL, 1, WL->NbPnts(), Standard_True);
1680         }
1681         
1682         if(myApprox2) {
1683           H2 = MakeBSpline2d(WL, 1, WL->NbPnts(), Standard_False);
1684         }
1685         //       
1686         mySeqOfCurve.Append(IntTools_Curve(MakeBSpline(WL,1,WL->NbPnts()), H1, H2));
1687       }
1688     } // if (!myApprox)
1689
1690     else { // myApprox=TRUE
1691       GeomInt_WLApprox theapp3d;
1692       // 
1693       Standard_Real tol2d = myTolApprox;
1694       //        
1695       theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, Standard_True);
1696       
1697       aNbParts=myLConstruct.NbParts();
1698       for (i=1; i<=aNbParts; i++) {
1699         myLConstruct.Part(i, fprm, lprm);
1700         Handle(IntPatch_WLine) WL = 
1701           convert.MakeWLine(Handle(IntPatch_ALine):: DownCast(L),fprm,lprm);
1702
1703         theapp3d.Perform(myHS1,myHS2,WL,Standard_True,myApprox1,myApprox2, 1, WL->NbPnts());
1704
1705         if (!theapp3d.IsDone()) {
1706           //
1707           Handle(Geom2d_BSplineCurve) H1;
1708           Handle(Geom2d_BSplineCurve) H2;
1709
1710           if(myApprox1) {
1711             H1 = MakeBSpline2d(WL, 1, WL->NbPnts(), Standard_True);
1712           }
1713           
1714           if(myApprox2) {
1715             H2 = MakeBSpline2d(WL, 1, WL->NbPnts(), Standard_False);
1716           }
1717           //     
1718           mySeqOfCurve.Append(IntTools_Curve(MakeBSpline(WL,1,WL->NbPnts()), H1, H2));
1719         }
1720
1721         else {
1722           if(myApprox1 || myApprox2) { 
1723             if( theapp3d.TolReached2d()>myTolReached2d || myTolReached2d==0) { 
1724               myTolReached2d = theapp3d.TolReached2d();
1725             }
1726           }
1727           
1728           if( theapp3d.TolReached3d()>myTolReached3d || myTolReached3d==0) { 
1729             myTolReached3d = theapp3d.TolReached3d();
1730           }
1731
1732           Standard_Integer aNbMultiCurves, nbpoles;
1733           aNbMultiCurves=theapp3d.NbMultiCurves();
1734           for (j=1; j<=aNbMultiCurves; j++) {
1735             const AppParCurves_MultiBSpCurve& mbspc = theapp3d.Value(j);
1736             nbpoles = mbspc.NbPoles();
1737             
1738             TColgp_Array1OfPnt tpoles(1, nbpoles);
1739             mbspc.Curve(1, tpoles);
1740             Handle(Geom_BSplineCurve) BS=new Geom_BSplineCurve(tpoles,
1741                                                                mbspc.Knots(),
1742                                                                mbspc.Multiplicities(),
1743                                                                mbspc.Degree());
1744             
1745             GeomLib_CheckBSplineCurve Check(BS,TOLCHECK,TOLANGCHECK);
1746             Check.FixTangent(Standard_True,Standard_True);
1747             // 
1748             IntTools_Curve aCurve;
1749             aCurve.SetCurve(BS);
1750             
1751             if(myApprox1) { 
1752               TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
1753               mbspc.Curve(2,tpoles2d);
1754               Handle(Geom2d_BSplineCurve) BS2=new Geom2d_BSplineCurve(tpoles2d,
1755                                                                       mbspc.Knots(),
1756                                                                       mbspc.Multiplicities(),
1757                                                                       mbspc.Degree());
1758
1759               GeomLib_Check2dBSplineCurve newCheck(BS2,TOLCHECK,TOLANGCHECK);
1760               newCheck.FixTangent(Standard_True,Standard_True);
1761               //                
1762               aCurve.SetFirstCurve2d(BS2);
1763             }
1764             else {
1765               Handle(Geom2d_BSplineCurve) H1;
1766               aCurve.SetFirstCurve2d(H1);
1767             }
1768             
1769             if(myApprox2) { 
1770               TColgp_Array1OfPnt2d tpoles2d(1, nbpoles);
1771               Standard_Integer TwoOrThree;
1772               TwoOrThree=myApprox1 ? 3 : 2;
1773               mbspc.Curve(TwoOrThree, tpoles2d);
1774               Handle(Geom2d_BSplineCurve) BS2 =new Geom2d_BSplineCurve(tpoles2d,
1775                                                                        mbspc.Knots(),
1776                                                                        mbspc.Multiplicities(),
1777                                                                        mbspc.Degree());
1778                 
1779               GeomLib_Check2dBSplineCurve newCheck(BS2,TOLCHECK,TOLANGCHECK);
1780               newCheck.FixTangent(Standard_True,Standard_True);
1781               //        
1782               aCurve.SetSecondCurve2d(BS2);
1783             }
1784             else { 
1785               Handle(Geom2d_BSplineCurve) H2;
1786               aCurve.SetSecondCurve2d(H2);
1787             }
1788             // 
1789             mySeqOfCurve.Append(aCurve);
1790
1791           }// for (j=1; j<=aNbMultiCurves; j++) {
1792         }// else from if (!theapp3d.IsDone())
1793       }// for (i=1; i<=aNbParts; i++) {
1794     }// else { // myApprox=TRUE
1795   }// case IntPatch_Analytic:
1796     break;
1797
1798   case IntPatch_Walking:{
1799     Handle(IntPatch_WLine) WL = 
1800       Handle(IntPatch_WLine)::DownCast(L);
1801     //
1802     Standard_Integer ifprm, ilprm;
1803     //
1804     if (!myApprox) {
1805       aNbParts = 1;
1806       if(!bAvoidLineConstructor){
1807         aNbParts=myLConstruct.NbParts();
1808       }
1809       for (i=1; i<=aNbParts; ++i) {
1810         Handle(Geom2d_BSplineCurve) H1, H2;
1811         Handle(Geom_Curve) aBSp;
1812         //
1813         if(bAvoidLineConstructor) {
1814           ifprm = 1;
1815           ilprm = WL->NbPnts();
1816         }
1817         else {
1818           myLConstruct.Part(i, fprm, lprm);
1819           ifprm=(Standard_Integer)fprm;
1820           ilprm=(Standard_Integer)lprm;
1821         }
1822         //
1823         if(myApprox1) {
1824           H1 = MakeBSpline2d(WL, ifprm, ilprm, Standard_True);
1825         }
1826         //
1827         if(myApprox2) {
1828           H2 = MakeBSpline2d(WL, ifprm, ilprm, Standard_False);
1829         }
1830         //        
1831         aBSp=MakeBSpline(WL, ifprm, ilprm);
1832         IntTools_Curve aIC(aBSp, H1, H2);
1833         mySeqOfCurve.Append(aIC);
1834       }// for (i=1; i<=aNbParts; ++i) {
1835     }// if (!myApprox) {
1836     //
1837     else { // X
1838       Standard_Boolean bIsDecomposited;
1839       Standard_Integer nbiter, aNbSeqOfL;
1840       Standard_Real tol2d;
1841       IntPatch_SequenceOfLine aSeqOfL;
1842       GeomInt_WLApprox theapp3d;
1843       Approx_ParametrizationType aParType = Approx_ChordLength;
1844       //
1845       Standard_Boolean anApprox1 = myApprox1;
1846       Standard_Boolean anApprox2 = myApprox2;
1847
1848       tol2d = myTolApprox;
1849
1850       GeomAbs_SurfaceType typs1, typs2;
1851       typs1 = myHS1->Surface().GetType();
1852       typs2 = myHS2->Surface().GetType();
1853       Standard_Boolean anWithPC = Standard_True;
1854
1855       if(typs1 == GeomAbs_Cylinder && typs2 == GeomAbs_Sphere) {
1856         anWithPC = 
1857           ApproxWithPCurves(myHS1->Surface().Cylinder(), myHS2->Surface().Sphere());
1858       }
1859       else if (typs1 == GeomAbs_Sphere && typs2 == GeomAbs_Cylinder) {
1860         anWithPC = 
1861           ApproxWithPCurves(myHS2->Surface().Cylinder(), myHS1->Surface().Sphere());
1862       }
1863       if(!anWithPC) {
1864         //aParType = Approx_Centripetal; 
1865         myTolApprox = 1.e-5; 
1866         anApprox1 = Standard_False;
1867         anApprox2 = Standard_False;
1868         //      
1869         tol2d = myTolApprox;
1870       }
1871         
1872       if(myHS1 == myHS2) { 
1873         //
1874         theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, Standard_False, aParType);
1875         rejectSurface = Standard_True;
1876       }
1877       else { 
1878         if(reApprox && !rejectSurface)
1879           theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, Standard_False, aParType);
1880         else {
1881           Standard_Integer iDegMax, iDegMin, iNbIter;
1882           //
1883           ApproxParameters(myHS1, myHS2, iDegMin, iDegMax, iNbIter);
1884           theapp3d.SetParameters(myTolApprox, tol2d, iDegMin, iDegMax, iNbIter, Standard_True, aParType);
1885           //
1886         }
1887       }
1888       //
1889       Standard_Real aReachedTol = Precision::Confusion();
1890       bIsDecomposited=DecompositionOfWLine(WL,
1891                                            myHS1, 
1892                                            myHS2, 
1893                                            myFace1, 
1894                                            myFace2, 
1895                                            myLConstruct, 
1896                                            bAvoidLineConstructor, 
1897                                            aSeqOfL, 
1898                                            aReachedTol,
1899                                            myContext);
1900       if ( bIsDecomposited && ( myTolReached3d < aReachedTol ) )
1901         myTolReached3d = aReachedTol;
1902
1903       //
1904       aNbSeqOfL=aSeqOfL.Length();
1905       //
1906       if (bIsDecomposited) {
1907         nbiter=aNbSeqOfL;
1908       }
1909       else {
1910         nbiter=1;
1911         aNbParts=1;
1912         if (!bAvoidLineConstructor) {
1913           aNbParts=myLConstruct.NbParts();
1914           nbiter=aNbParts;
1915         }
1916       }
1917       //
1918       // nbiter=(bIsDecomposited) ? aSeqOfL.Length() :
1919       //   ((bAvoidLineConstructor) ? 1 :aNbParts);
1920       //
1921       for(i = 1; i <= nbiter; ++i) {
1922         if(bIsDecomposited) {
1923           WL = Handle(IntPatch_WLine)::DownCast(aSeqOfL.Value(i));
1924           ifprm = 1;
1925           ilprm = WL->NbPnts();
1926         }
1927         else {
1928           if(bAvoidLineConstructor) {
1929             ifprm = 1;
1930             ilprm = WL->NbPnts();
1931           }
1932           else {
1933             myLConstruct.Part(i, fprm, lprm);
1934             ifprm = (Standard_Integer)fprm;
1935             ilprm = (Standard_Integer)lprm;
1936           }
1937         }
1938         //-- lbr : 
1939         //-- Si une des surfaces est un plan , on approxime en 2d
1940         //-- sur cette surface et on remonte les points 2d en 3d.
1941         if(typs1 == GeomAbs_Plane) { 
1942           theapp3d.Perform(myHS1, myHS2, WL, Standard_False,Standard_True, myApprox2,ifprm,ilprm);
1943         }         
1944         else if(typs2 == GeomAbs_Plane) { 
1945           theapp3d.Perform(myHS1,myHS2,WL,Standard_False,myApprox1,Standard_True,ifprm,ilprm);
1946         }
1947         else { 
1948           //
1949           if (myHS1 != myHS2){
1950             if ((typs1==GeomAbs_BezierSurface || typs1==GeomAbs_BSplineSurface) &&
1951                 (typs2==GeomAbs_BezierSurface || typs2==GeomAbs_BSplineSurface)) {
1952              
1953               theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, Standard_True, aParType);
1954               
1955               Standard_Boolean bUseSurfaces;
1956               bUseSurfaces=NotUseSurfacesForApprox(myFace1, myFace2, WL, ifprm,  ilprm);
1957               if (bUseSurfaces) {
1958                 // ######
1959                 rejectSurface = Standard_True;
1960                 // ######
1961                 theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, Standard_False, aParType);
1962               }
1963             }
1964           }
1965           //
1966           theapp3d.Perform(myHS1,myHS2,WL,Standard_True,anApprox1,anApprox2,ifprm,ilprm);
1967         }
1968          
1969         if (!theapp3d.IsDone()) {
1970           //      
1971           Handle(Geom2d_BSplineCurve) H1;
1972           //      
1973           Handle(Geom_Curve) aBSp=MakeBSpline(WL,ifprm, ilprm);
1974           Handle(Geom2d_BSplineCurve) H2;
1975
1976           if(myApprox1) {
1977             H1 = MakeBSpline2d(WL, ifprm, ilprm, Standard_True);
1978           }
1979           
1980           if(myApprox2) {
1981             H2 = MakeBSpline2d(WL, ifprm, ilprm, Standard_False);
1982           }
1983           //      
1984           IntTools_Curve aIC(aBSp, H1, H2);
1985           mySeqOfCurve.Append(aIC);
1986         }
1987         
1988         else {
1989           if(myApprox1 || myApprox2 || (typs1==GeomAbs_Plane || typs2==GeomAbs_Plane)) { 
1990             if( theapp3d.TolReached2d()>myTolReached2d || myTolReached2d==0.) { 
1991               myTolReached2d = theapp3d.TolReached2d();
1992             }
1993           }
1994           if(typs1==GeomAbs_Plane || typs2==GeomAbs_Plane) { 
1995             myTolReached3d = myTolReached2d;
1996             //
1997             if (typs1==GeomAbs_Torus || typs2==GeomAbs_Torus) {
1998               if (myTolReached3d<1.e-6) {
1999                 myTolReached3d = theapp3d.TolReached3d();
2000                 myTolReached3d=1.e-6;
2001               }
2002             }
2003             //
2004           }
2005           else  if( theapp3d.TolReached3d()>myTolReached3d || myTolReached3d==0.) { 
2006             myTolReached3d = theapp3d.TolReached3d();
2007           }
2008           
2009           Standard_Integer aNbMultiCurves, nbpoles;
2010           aNbMultiCurves=theapp3d.NbMultiCurves(); 
2011           for (j=1; j<=aNbMultiCurves; j++) {
2012             if(typs1 == GeomAbs_Plane) {
2013               const AppParCurves_MultiBSpCurve& mbspc = theapp3d.Value(j);
2014               nbpoles = mbspc.NbPoles();
2015               
2016               TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
2017               TColgp_Array1OfPnt   tpoles(1,nbpoles);
2018               
2019               mbspc.Curve(1,tpoles2d);
2020               const gp_Pln&  Pln = myHS1->Surface().Plane();
2021               //
2022               Standard_Integer ik; 
2023               for(ik = 1; ik<= nbpoles; ik++) { 
2024                 tpoles.SetValue(ik,
2025                                 ElSLib::Value(tpoles2d.Value(ik).X(),
2026                                               tpoles2d.Value(ik).Y(),
2027                                               Pln));
2028               }
2029               //
2030               Handle(Geom_BSplineCurve) BS = 
2031                 new Geom_BSplineCurve(tpoles,
2032                                       mbspc.Knots(),
2033                                       mbspc.Multiplicities(),
2034                                       mbspc.Degree());
2035               GeomLib_CheckBSplineCurve Check(BS,TOLCHECK,TOLANGCHECK);
2036               Check.FixTangent(Standard_True, Standard_True);
2037               //        
2038               IntTools_Curve aCurve;
2039               aCurve.SetCurve(BS);
2040
2041               if(myApprox1) { 
2042                 Handle(Geom2d_BSplineCurve) BS1 = 
2043                   new Geom2d_BSplineCurve(tpoles2d,
2044                                           mbspc.Knots(),
2045                                           mbspc.Multiplicities(),
2046                                           mbspc.Degree());
2047                 GeomLib_Check2dBSplineCurve Check1(BS1,TOLCHECK,TOLANGCHECK);
2048                 Check1.FixTangent(Standard_True,Standard_True);
2049                 //
2050                 // ############################################
2051                 if(!rejectSurface && !reApprox) {
2052                   Standard_Boolean isValid = IsCurveValid(BS1);
2053                   if(!isValid) {
2054                     reApprox = Standard_True;
2055                     goto reapprox;
2056                   }
2057                 }
2058                 // ############################################
2059                 aCurve.SetFirstCurve2d(BS1);
2060               }
2061               else {
2062                 Handle(Geom2d_BSplineCurve) H1;
2063                 aCurve.SetFirstCurve2d(H1);
2064               }
2065
2066               if(myApprox2) { 
2067                 mbspc.Curve(2, tpoles2d);
2068                 
2069                 Handle(Geom2d_BSplineCurve) BS2 = new Geom2d_BSplineCurve(tpoles2d,
2070                                                                           mbspc.Knots(),
2071                                                                           mbspc.Multiplicities(),
2072                                                                           mbspc.Degree());
2073                 GeomLib_Check2dBSplineCurve newCheck(BS2,TOLCHECK,TOLANGCHECK);
2074                 newCheck.FixTangent(Standard_True,Standard_True);
2075                 
2076                 // ###########################################
2077                 if(!rejectSurface && !reApprox) {
2078                   Standard_Boolean isValid = IsCurveValid(BS2);
2079                   if(!isValid) {
2080                     reApprox = Standard_True;
2081                     goto reapprox;
2082                   }
2083                 }
2084                 // ###########################################
2085                 // 
2086                 aCurve.SetSecondCurve2d(BS2);
2087               }
2088               else { 
2089                 Handle(Geom2d_BSplineCurve) H2;
2090                 //              
2091                   aCurve.SetSecondCurve2d(H2);
2092               }
2093               //
2094               mySeqOfCurve.Append(aCurve);
2095             }
2096             
2097             else if(typs2 == GeomAbs_Plane) { 
2098               const AppParCurves_MultiBSpCurve& mbspc = theapp3d.Value(j);
2099               nbpoles = mbspc.NbPoles();
2100               
2101               TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
2102               TColgp_Array1OfPnt   tpoles(1,nbpoles);
2103               mbspc.Curve((myApprox1==Standard_True)? 2 : 1,tpoles2d);
2104               const gp_Pln&  Pln = myHS2->Surface().Plane();
2105               //
2106               Standard_Integer ik; 
2107               for(ik = 1; ik<= nbpoles; ik++) { 
2108                 tpoles.SetValue(ik,
2109                                 ElSLib::Value(tpoles2d.Value(ik).X(),
2110                                               tpoles2d.Value(ik).Y(),
2111                                               Pln));
2112                 
2113               }
2114               //
2115               Handle(Geom_BSplineCurve) BS=new Geom_BSplineCurve(tpoles,
2116                                                                  mbspc.Knots(),
2117                                                                  mbspc.Multiplicities(),
2118                                                                  mbspc.Degree());
2119               GeomLib_CheckBSplineCurve Check(BS,TOLCHECK,TOLANGCHECK);
2120               Check.FixTangent(Standard_True,Standard_True);
2121               //        
2122               IntTools_Curve aCurve;
2123               aCurve.SetCurve(BS);
2124
2125               if(myApprox2) {
2126                 Handle(Geom2d_BSplineCurve) BS1=new Geom2d_BSplineCurve(tpoles2d,
2127                                                                         mbspc.Knots(),
2128                                                                         mbspc.Multiplicities(),
2129                                                                         mbspc.Degree());
2130                 GeomLib_Check2dBSplineCurve Check1(BS1,TOLCHECK,TOLANGCHECK);
2131                 Check1.FixTangent(Standard_True,Standard_True);
2132                 //      
2133                 // ###########################################
2134                 if(!rejectSurface && !reApprox) {
2135                   Standard_Boolean isValid = IsCurveValid(BS1);
2136                   if(!isValid) {
2137                     reApprox = Standard_True;
2138                     goto reapprox;
2139                   }
2140                 }
2141                 // ###########################################
2142                 bPCurvesOk = CheckPCurve(BS1, myFace2);
2143                 aCurve.SetSecondCurve2d(BS1);
2144               }
2145               else {
2146                 Handle(Geom2d_BSplineCurve) H2;
2147                 aCurve.SetSecondCurve2d(H2);
2148               }
2149               
2150               if(myApprox1) { 
2151                 mbspc.Curve(1,tpoles2d);
2152                 Handle(Geom2d_BSplineCurve) BS2=new Geom2d_BSplineCurve(tpoles2d,
2153                                                                         mbspc.Knots(),
2154                                                                         mbspc.Multiplicities(),
2155                                                                         mbspc.Degree());
2156                 GeomLib_Check2dBSplineCurve Check2(BS2,TOLCHECK,TOLANGCHECK);
2157                 Check2.FixTangent(Standard_True,Standard_True);
2158                 //
2159                 // ###########################################
2160                 if(!rejectSurface && !reApprox) {
2161                   Standard_Boolean isValid = IsCurveValid(BS2);
2162                   if(!isValid) {
2163                     reApprox = Standard_True;
2164                     goto reapprox;
2165                   }
2166                 }
2167                 // ###########################################
2168                 bPCurvesOk = bPCurvesOk && CheckPCurve(BS2, myFace1);
2169                 aCurve.SetFirstCurve2d(BS2);
2170               }
2171               else { 
2172                 Handle(Geom2d_BSplineCurve) H1;
2173                 //              
2174                 aCurve.SetFirstCurve2d(H1);
2175               }
2176               //
2177               //if points of the pcurves are out of the faces bounds
2178               //create 3d and 2d curves without approximation
2179               if (!bPCurvesOk) {
2180                 Handle(Geom2d_BSplineCurve) H1, H2;
2181                 bPCurvesOk = Standard_True;
2182                 //        
2183                 Handle(Geom_Curve) aBSp=MakeBSpline(WL,ifprm, ilprm);
2184                 
2185                 if(myApprox1) {
2186                   H1 = MakeBSpline2d(WL, ifprm, ilprm, Standard_True);
2187                   bPCurvesOk = CheckPCurve(H1, myFace1);
2188                 }
2189                 
2190                 if(myApprox2) {
2191                   H2 = MakeBSpline2d(WL, ifprm, ilprm, Standard_False);
2192                   bPCurvesOk = bPCurvesOk && CheckPCurve(H2, myFace2);
2193                 }
2194                 //
2195                 //if pcurves created without approximation are out of the 
2196                 //faces bounds, use approximated 3d and 2d curves
2197                 if (bPCurvesOk) {
2198                   IntTools_Curve aIC(aBSp, H1, H2);
2199                   mySeqOfCurve.Append(aIC);
2200                 } else {
2201                   mySeqOfCurve.Append(aCurve);
2202                 }
2203               } else {
2204                 mySeqOfCurve.Append(aCurve);
2205               }
2206             }
2207             else { 
2208               const AppParCurves_MultiBSpCurve& mbspc = theapp3d.Value(j);
2209               nbpoles = mbspc.NbPoles();
2210               TColgp_Array1OfPnt tpoles(1,nbpoles);
2211               mbspc.Curve(1,tpoles);
2212               Handle(Geom_BSplineCurve) BS=new Geom_BSplineCurve(tpoles,
2213                                                                  mbspc.Knots(),
2214                                                                  mbspc.Multiplicities(),
2215                                                                  mbspc.Degree());
2216               GeomLib_CheckBSplineCurve Check(BS,TOLCHECK,TOLANGCHECK);
2217               Check.FixTangent(Standard_True,Standard_True);
2218               //                
2219               IntTools_Curve aCurve;
2220               aCurve.SetCurve(BS);
2221               
2222               if(myApprox1) { 
2223                 if(anApprox1) {
2224                   TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
2225                   mbspc.Curve(2,tpoles2d);
2226                   Handle(Geom2d_BSplineCurve) BS1=new Geom2d_BSplineCurve(tpoles2d,
2227                                                                         mbspc.Knots(),
2228                                                                         mbspc.Multiplicities(),
2229                                                                         mbspc.Degree());
2230                   GeomLib_Check2dBSplineCurve newCheck(BS1,TOLCHECK,TOLANGCHECK);
2231                   newCheck.FixTangent(Standard_True,Standard_True);
2232                   //    
2233                   aCurve.SetFirstCurve2d(BS1);
2234                 }
2235                 else {
2236                   Handle(Geom2d_BSplineCurve) BS1;
2237                   fprm = BS->FirstParameter();
2238                   lprm = BS->LastParameter();
2239
2240                   Handle(Geom2d_Curve) C2d;
2241                   Standard_Real aTol = myTolApprox;
2242                   BuildPCurves(fprm, lprm, aTol, myHS1->ChangeSurface().Surface(), BS, C2d);
2243                   BS1 = Handle(Geom2d_BSplineCurve)::DownCast(C2d);
2244                   aCurve.SetFirstCurve2d(BS1);
2245                 }
2246                 
2247               }
2248               else {
2249                 Handle(Geom2d_BSplineCurve) H1;
2250                 //              
2251                 aCurve.SetFirstCurve2d(H1);
2252               }
2253               if(myApprox2) { 
2254                 if(anApprox2) {
2255                   TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
2256                   mbspc.Curve((myApprox1==Standard_True)? 3 : 2,tpoles2d);
2257                   Handle(Geom2d_BSplineCurve) BS2=new Geom2d_BSplineCurve(tpoles2d,
2258                                                                         mbspc.Knots(),
2259                                                                         mbspc.Multiplicities(),
2260                                                                         mbspc.Degree());
2261                   GeomLib_Check2dBSplineCurve newCheck(BS2,TOLCHECK,TOLANGCHECK);
2262                   newCheck.FixTangent(Standard_True,Standard_True);
2263                 //              
2264                   aCurve.SetSecondCurve2d(BS2);
2265                 }
2266                 else {
2267                   Handle(Geom2d_BSplineCurve) BS2;
2268                   fprm = BS->FirstParameter();
2269                   lprm = BS->LastParameter();
2270
2271                   Handle(Geom2d_Curve) C2d;
2272                   Standard_Real aTol = myTolApprox;
2273                   BuildPCurves(fprm, lprm, aTol, myHS2->ChangeSurface().Surface(), BS, C2d);
2274                   BS2 = Handle(Geom2d_BSplineCurve)::DownCast(C2d);
2275                   aCurve.SetSecondCurve2d(BS2);
2276                 }
2277                 
2278               }
2279               else { 
2280                 Handle(Geom2d_BSplineCurve) H2;
2281                 //              
2282                 aCurve.SetSecondCurve2d(H2);
2283               }
2284               //
2285               mySeqOfCurve.Append(aCurve);
2286             }
2287           }
2288         }
2289       }
2290     }// else { // X
2291   }// case IntPatch_Walking:{
2292     break;
2293     
2294   case IntPatch_Restriction: 
2295     break;
2296
2297   }
2298 }
2299
2300 //=======================================================================
2301 //function : BuildPCurves
2302 //purpose  : 
2303 //=======================================================================
2304  void BuildPCurves (Standard_Real f,
2305                     Standard_Real l,
2306                     Standard_Real& Tol,
2307                     const Handle (Geom_Surface)& S,
2308                     const Handle (Geom_Curve)&   C,
2309                     Handle (Geom2d_Curve)& C2d)
2310 {
2311
2312   Standard_Real umin,umax,vmin,vmax;
2313   // 
2314
2315   if (C2d.IsNull()) {
2316
2317     // in class ProjLib_Function the range of parameters is shrank by 1.e-09
2318     if((l - f) > 2.e-09) {
2319       C2d = GeomProjLib::Curve2d(C,f,l,S,Tol);
2320       //
2321       if (C2d.IsNull()) {
2322         // proj. a circle that goes through the pole on a sphere to the sphere     
2323         Tol=Tol+1.e-7;
2324         C2d = GeomProjLib::Curve2d(C,f,l,S,Tol);
2325       }
2326     }
2327     else {
2328       if((l - f) > Epsilon(Abs(f))) {
2329         GeomAPI_ProjectPointOnSurf aProjector1, aProjector2;
2330         gp_Pnt P1 = C->Value(f);
2331         gp_Pnt P2 = C->Value(l);
2332         aProjector1.Init(P1, S);
2333         aProjector2.Init(P2, S);
2334
2335         if(aProjector1.IsDone() && aProjector2.IsDone()) {
2336           Standard_Real U=0., V=0.;
2337           aProjector1.LowerDistanceParameters(U, V);
2338           gp_Pnt2d p1(U, V);
2339
2340           aProjector2.LowerDistanceParameters(U, V);
2341           gp_Pnt2d p2(U, V);
2342
2343           if(p1.Distance(p2) > gp::Resolution()) {
2344             TColgp_Array1OfPnt2d poles(1,2);
2345             TColStd_Array1OfReal knots(1,2);
2346             TColStd_Array1OfInteger mults(1,2);
2347             poles(1) = p1;
2348             poles(2) = p2;
2349             knots(1) = f;
2350             knots(2) = l;
2351             mults(1) = mults(2) = 2;
2352
2353             C2d = new Geom2d_BSplineCurve(poles,knots,mults,1);
2354
2355             // compute reached tolerance.begin
2356             gp_Pnt PMid = C->Value((f + l) * 0.5);
2357             aProjector1.Perform(PMid);
2358
2359             if(aProjector1.IsDone()) {
2360               aProjector1.LowerDistanceParameters(U, V);
2361               gp_Pnt2d pmidproj(U, V);
2362               gp_Pnt2d pmidcurve2d = C2d->Value((f + l) * 0.5);
2363               Standard_Real adist = pmidcurve2d.Distance(pmidproj);
2364               Tol = (adist > Tol) ? adist : Tol;
2365             }
2366             // compute reached tolerance.end
2367           }
2368         }
2369       }
2370     }
2371     //
2372     S->Bounds(umin, umax, vmin, vmax);
2373
2374     if (S->IsUPeriodic() && !C2d.IsNull()) {
2375       // Recadre dans le domaine UV de la face
2376       Standard_Real period, U0, du, aEps; 
2377       
2378       du =0.0;
2379       aEps=Precision::PConfusion();
2380       period = S->UPeriod();
2381       gp_Pnt2d Pf = C2d->Value(f);
2382       U0=Pf.X();
2383       //
2384       gp_Pnt2d Pl = C2d->Value(l);
2385       
2386       U0 = Min(Pl.X(), U0);
2387 //       while(U0-umin<aEps) { 
2388       while(U0-umin<-aEps) { 
2389         U0+=period;
2390         du+=period;
2391       }
2392       //
2393       while(U0-umax>aEps) { 
2394         U0-=period;
2395         du-=period;
2396       }
2397       if (du != 0) {
2398         gp_Vec2d T1(du,0.);
2399         C2d->Translate(T1);
2400       }
2401     }
2402   }
2403 }
2404
2405 //=======================================================================
2406 //function : Parameters
2407 //purpose  : 
2408 //=======================================================================
2409  void Parameters(const Handle(GeomAdaptor_HSurface)& HS1,
2410                  const Handle(GeomAdaptor_HSurface)& HS2,
2411                  const gp_Pnt& Ptref,
2412                  Standard_Real& U1,
2413                  Standard_Real& V1,
2414                  Standard_Real& U2,
2415                  Standard_Real& V2)
2416 {
2417
2418   IntSurf_Quadric quad1,quad2;
2419   GeomAbs_SurfaceType typs = HS1->Surface().GetType();
2420
2421   switch (typs) {
2422   case GeomAbs_Plane:
2423     quad1.SetValue(HS1->Surface().Plane());
2424     break;
2425   case GeomAbs_Cylinder:
2426     quad1.SetValue(HS1->Surface().Cylinder());
2427     break;
2428   case GeomAbs_Cone:
2429     quad1.SetValue(HS1->Surface().Cone());
2430     break;
2431   case GeomAbs_Sphere:
2432     quad1.SetValue(HS1->Surface().Sphere());
2433     break;
2434   default:
2435     Standard_ConstructionError::Raise("GeomInt_IntSS::MakeCurve");
2436   }
2437   
2438   typs = HS2->Surface().GetType();
2439   switch (typs) {
2440   case GeomAbs_Plane:
2441     quad2.SetValue(HS2->Surface().Plane());
2442     break;
2443   case GeomAbs_Cylinder:
2444     quad2.SetValue(HS2->Surface().Cylinder());
2445     break;
2446   case GeomAbs_Cone:
2447     quad2.SetValue(HS2->Surface().Cone());
2448     break;
2449   case GeomAbs_Sphere:
2450     quad2.SetValue(HS2->Surface().Sphere());
2451     break;
2452   default:
2453     Standard_ConstructionError::Raise("GeomInt_IntSS::MakeCurve");
2454   }
2455
2456   quad1.Parameters(Ptref,U1,V1);
2457   quad2.Parameters(Ptref,U2,V2);
2458 }
2459
2460 //=======================================================================
2461 //function : MakeBSpline
2462 //purpose  : 
2463 //=======================================================================
2464 Handle(Geom_Curve) MakeBSpline  (const Handle(IntPatch_WLine)& WL,
2465                                  const Standard_Integer ideb,
2466                                  const Standard_Integer ifin)
2467 {
2468   Standard_Integer i,nbpnt = ifin-ideb+1;
2469   TColgp_Array1OfPnt poles(1,nbpnt);
2470   TColStd_Array1OfReal knots(1,nbpnt);
2471   TColStd_Array1OfInteger mults(1,nbpnt);
2472   Standard_Integer ipidebm1;
2473   for(i=1,ipidebm1=i+ideb-1; i<=nbpnt;ipidebm1++, i++) {
2474     poles(i) = WL->Point(ipidebm1).Value();
2475     mults(i) = 1;
2476     knots(i) = i-1;
2477   }
2478   mults(1) = mults(nbpnt) = 2;
2479   return
2480     new Geom_BSplineCurve(poles,knots,mults,1);
2481 }
2482 //
2483
2484 //=======================================================================
2485 //function : MakeBSpline2d
2486 //purpose  : 
2487 //=======================================================================
2488 Handle(Geom2d_BSplineCurve) MakeBSpline2d(const Handle(IntPatch_WLine)& theWLine,
2489                                           const Standard_Integer ideb,
2490                                           const Standard_Integer ifin,
2491                                           const Standard_Boolean onFirst)
2492 {
2493   Standard_Integer i, nbpnt = ifin-ideb+1;
2494   TColgp_Array1OfPnt2d poles(1,nbpnt);
2495   TColStd_Array1OfReal knots(1,nbpnt);
2496   TColStd_Array1OfInteger mults(1,nbpnt);
2497   Standard_Integer ipidebm1;
2498
2499   for(i = 1, ipidebm1 = i+ideb-1; i <= nbpnt; ipidebm1++, i++) {
2500       Standard_Real U, V;
2501       if(onFirst)
2502         theWLine->Point(ipidebm1).ParametersOnS1(U, V);
2503       else
2504         theWLine->Point(ipidebm1).ParametersOnS2(U, V);
2505       poles(i).SetCoord(U, V);
2506       mults(i) = 1;
2507       knots(i) = i-1;
2508     }
2509     mults(1) = mults(nbpnt) = 2;
2510
2511   return new Geom2d_BSplineCurve(poles,knots,mults,1);
2512 }
2513 //=======================================================================
2514 //function : PrepareLines3D
2515 //purpose  : 
2516 //=======================================================================
2517   void IntTools_FaceFace::PrepareLines3D(const Standard_Boolean bToSplit)
2518 {
2519   Standard_Integer i, aNbCurves;
2520   GeomAbs_SurfaceType aType1, aType2;
2521   IntTools_SequenceOfCurves aNewCvs;
2522   //
2523   // 1. Treatment closed  curves
2524   aNbCurves=mySeqOfCurve.Length();
2525   for (i=1; i<=aNbCurves; ++i) {
2526     const IntTools_Curve& aIC=mySeqOfCurve(i);
2527     //
2528     if (bToSplit) {
2529       Standard_Integer j, aNbC;
2530       IntTools_SequenceOfCurves aSeqCvs;
2531       //
2532       aNbC=IntTools_Tools::SplitCurve(aIC, aSeqCvs);
2533       if (aNbC) {
2534         for (j=1; j<=aNbC; ++j) {
2535           const IntTools_Curve& aICNew=aSeqCvs(j);
2536           aNewCvs.Append(aICNew);
2537         }
2538       }
2539       else {
2540         aNewCvs.Append(aIC);
2541       }
2542     }
2543     else {
2544       aNewCvs.Append(aIC);
2545     }
2546   }
2547   //
2548   // 2. Plane\Cone intersection when we had 4 curves
2549   aType1=myHS1->GetType();
2550   aType2=myHS2->GetType();
2551   aNbCurves=aNewCvs.Length();
2552   //
2553   if ((aType1==GeomAbs_Plane && aType2==GeomAbs_Cone) ||
2554       (aType2==GeomAbs_Plane && aType1==GeomAbs_Cone)) {
2555     if (aNbCurves==4) {
2556       GeomAbs_CurveType aCType1;
2557       //
2558       aCType1=aNewCvs(1).Type();
2559       if (aCType1==GeomAbs_Line) {
2560         IntTools_SequenceOfCurves aSeqIn, aSeqOut;
2561         //
2562         for (i=1; i<=aNbCurves; ++i) {
2563           const IntTools_Curve& aIC=aNewCvs(i);
2564           aSeqIn.Append(aIC);
2565         }
2566         //
2567         IntTools_Tools::RejectLines(aSeqIn, aSeqOut);
2568         //
2569         aNewCvs.Clear();
2570         aNbCurves=aSeqOut.Length(); 
2571         for (i=1; i<=aNbCurves; ++i) {
2572           const IntTools_Curve& aIC=aSeqOut(i);
2573           aNewCvs.Append(aIC);
2574         }
2575       }
2576     }
2577   }// if ((aType1==GeomAbs_Plane && aType2==GeomAbs_Cone)...
2578   //
2579   // 3. Fill  mySeqOfCurve
2580   mySeqOfCurve.Clear();
2581   aNbCurves=aNewCvs.Length();
2582   for (i=1; i<=aNbCurves; ++i) {
2583     const IntTools_Curve& aIC=aNewCvs(i);
2584     mySeqOfCurve.Append(aIC);
2585   }
2586 }
2587 //=======================================================================
2588 //function : CorrectSurfaceBoundaries
2589 //purpose  : 
2590 //=======================================================================
2591  void CorrectSurfaceBoundaries(const TopoDS_Face&  theFace,
2592                                const Standard_Real theTolerance,
2593                                Standard_Real&      theumin,
2594                                Standard_Real&      theumax, 
2595                                Standard_Real&      thevmin, 
2596                                Standard_Real&      thevmax) 
2597 {
2598   Standard_Boolean enlarge, isuperiodic, isvperiodic;
2599   Standard_Real uinf, usup, vinf, vsup, delta;
2600   GeomAbs_SurfaceType aType;
2601   Handle(Geom_Surface) aSurface;
2602   //
2603   aSurface = BRep_Tool::Surface(theFace);
2604   aSurface->Bounds(uinf, usup, vinf, vsup);
2605   delta = theTolerance;
2606   enlarge = Standard_False;
2607   //
2608   GeomAdaptor_Surface anAdaptorSurface(aSurface);
2609   //
2610   if(aSurface->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface))) {
2611     Handle(Geom_Surface) aBasisSurface = 
2612       (Handle(Geom_RectangularTrimmedSurface)::DownCast(aSurface))->BasisSurface();
2613     
2614     if(aBasisSurface->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface)) ||
2615        aBasisSurface->IsKind(STANDARD_TYPE(Geom_OffsetSurface))) {
2616       return;
2617     }
2618   }
2619   //
2620   if(aSurface->IsKind(STANDARD_TYPE(Geom_OffsetSurface))) {
2621     Handle(Geom_Surface) aBasisSurface = 
2622       (Handle(Geom_OffsetSurface)::DownCast(aSurface))->BasisSurface();
2623     
2624     if(aBasisSurface->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface)) ||
2625        aBasisSurface->IsKind(STANDARD_TYPE(Geom_OffsetSurface))) {
2626       return;
2627     }
2628   }
2629   //
2630   isuperiodic = anAdaptorSurface.IsUPeriodic();
2631   isvperiodic = anAdaptorSurface.IsVPeriodic();
2632   //
2633   aType=anAdaptorSurface.GetType();
2634   if((aType==GeomAbs_BezierSurface) ||
2635      (aType==GeomAbs_BSplineSurface) ||
2636      (aType==GeomAbs_SurfaceOfExtrusion) ||
2637      (aType==GeomAbs_SurfaceOfRevolution)) {
2638     enlarge=Standard_True;
2639   }
2640   //
2641   if(!isuperiodic && enlarge) {
2642
2643     if((theumin - uinf) > delta )
2644       theumin -= delta;
2645     else {
2646       theumin = uinf;
2647     }
2648
2649     if((usup - theumax) > delta )
2650       theumax += delta;
2651     else
2652       theumax = usup;
2653   }
2654   //
2655   if(!isvperiodic && enlarge) {
2656     if((thevmin - vinf) > delta ) {
2657       thevmin -= delta;
2658     }
2659     else { 
2660       thevmin = vinf;
2661     }
2662     if((vsup - thevmax) > delta ) {
2663       thevmax += delta;
2664     }
2665     else {
2666       thevmax = vsup;
2667     }
2668   }
2669   //
2670   {
2671     Standard_Integer aNbP;
2672     Standard_Real aXP, dXfact, aXmid, aX1, aX2, aTolPA;
2673     //
2674     aTolPA=Precision::Angular();
2675     // U
2676     if (isuperiodic) {
2677       aXP=anAdaptorSurface.UPeriod();
2678       dXfact=theumax-theumin;
2679       if (dXfact-aTolPA>aXP) {
2680         aXmid=0.5*(theumax+theumin);
2681         aNbP=RealToInt(aXmid/aXP);
2682         if (aXmid<0.) {
2683           aNbP=aNbP-1;
2684         }
2685         aX1=aNbP*aXP;
2686         if (theumin>aTolPA) {
2687           aX1=theumin+aNbP*aXP;
2688         }
2689         aX2=aX1+aXP;
2690         if (theumin<aX1) {
2691           theumin=aX1;
2692         }
2693         if (theumax>aX2) {
2694           theumax=aX2;
2695         }
2696       }
2697     }
2698     // V
2699     if (isvperiodic) {
2700       aXP=anAdaptorSurface.VPeriod();
2701       dXfact=thevmax-thevmin;
2702       if (dXfact-aTolPA>aXP) {
2703         aXmid=0.5*(thevmax+thevmin);
2704         aNbP=RealToInt(aXmid/aXP);
2705         if (aXmid<0.) {
2706           aNbP=aNbP-1;
2707         }
2708         aX1=aNbP*aXP;
2709         if (thevmin>aTolPA) {
2710           aX1=thevmin+aNbP*aXP;
2711         }
2712         aX2=aX1+aXP;
2713         if (thevmin<aX1) {
2714           thevmin=aX1;
2715         }
2716         if (thevmax>aX2) {
2717           thevmax=aX2;
2718         }
2719       }
2720     }
2721   }
2722   //
2723   if(isuperiodic || isvperiodic) {
2724     Standard_Boolean correct = Standard_False;
2725     Standard_Boolean correctU = Standard_False;
2726     Standard_Boolean correctV = Standard_False;
2727     Bnd_Box2d aBox;
2728     TopExp_Explorer anExp;
2729
2730     for(anExp.Init(theFace, TopAbs_EDGE); anExp.More(); anExp.Next()) {
2731       if(BRep_Tool::IsClosed(TopoDS::Edge(anExp.Current()), theFace)) {
2732         correct = Standard_True;
2733         Standard_Real f, l;
2734         TopoDS_Edge anEdge = TopoDS::Edge(anExp.Current());
2735         
2736         for(Standard_Integer i = 0; i < 2; i++) {
2737           if(i==0) {
2738             anEdge.Orientation(TopAbs_FORWARD);
2739           }
2740           else {
2741             anEdge.Orientation(TopAbs_REVERSED);
2742           }
2743           Handle(Geom2d_Curve) aCurve = BRep_Tool::CurveOnSurface(anEdge, theFace, f, l);
2744           
2745           if(aCurve.IsNull()) {
2746             correct = Standard_False;
2747             break;
2748           }
2749           Handle(Geom2d_Line) aLine = Handle(Geom2d_Line)::DownCast(aCurve);
2750
2751           if(aLine.IsNull()) {
2752             correct = Standard_False;
2753             break;
2754           }
2755           gp_Dir2d anUDir(1., 0.);
2756           gp_Dir2d aVDir(0., 1.);
2757           Standard_Real anAngularTolerance = Precision::Angular();
2758
2759           correctU = correctU || aLine->Position().Direction().IsParallel(aVDir, anAngularTolerance);
2760           correctV = correctV || aLine->Position().Direction().IsParallel(anUDir, anAngularTolerance);
2761           
2762           gp_Pnt2d pp1 = aCurve->Value(f);
2763           aBox.Add(pp1);
2764           gp_Pnt2d pp2 = aCurve->Value(l);
2765           aBox.Add(pp2);
2766         }
2767         if(!correct)
2768           break;
2769       }
2770     }
2771
2772     if(correct) {
2773       Standard_Real umin, vmin, umax, vmax;
2774       aBox.Get(umin, vmin, umax, vmax);
2775
2776       if(isuperiodic && correctU) {
2777         
2778         if(theumin < umin)
2779           theumin = umin;
2780         
2781         if(theumax > umax) {
2782           theumax = umax;
2783         }
2784       }
2785       if(isvperiodic && correctV) {
2786         
2787         if(thevmin < vmin)
2788           thevmin = vmin;
2789         if(thevmax > vmax)
2790           thevmax = vmax;
2791       }
2792     }
2793   }
2794 }
2795 //
2796 //
2797 // The block is dedicated to determine whether WLine [ifprm, ilprm]
2798 // crosses the degenerated zone on each given surface or not.
2799 // If Yes -> We will not use info about surfaces during approximation
2800 // because inside degenerated zone of the surface the approx. algo.
2801 // uses wrong values of normal, etc., and resulting curve will have
2802 // oscillations that we would not like to have. 
2803  
2804
2805
2806 static
2807   Standard_Boolean IsDegeneratedZone(const gp_Pnt2d& aP2d,
2808                                      const Handle(Geom_Surface)& aS,
2809                                      const Standard_Integer iDir);
2810 static
2811   Standard_Boolean IsPointInDegeneratedZone(const IntSurf_PntOn2S& aP2S,
2812                                             const TopoDS_Face& aF1,
2813                                             const TopoDS_Face& aF2);
2814 //=======================================================================
2815 //function :  NotUseSurfacesForApprox
2816 //purpose  : 
2817 //=======================================================================
2818 Standard_Boolean NotUseSurfacesForApprox(const TopoDS_Face& aF1,
2819                                          const TopoDS_Face& aF2,
2820                                          const Handle(IntPatch_WLine)& WL,
2821                                          const Standard_Integer ifprm,
2822                                          const Standard_Integer ilprm)
2823 {
2824   Standard_Boolean bPInDZ;
2825
2826   Handle(IntSurf_LineOn2S) aLineOn2S=WL->Curve();
2827   
2828   const IntSurf_PntOn2S& aP2Sfprm=aLineOn2S->Value(ifprm);
2829   bPInDZ=IsPointInDegeneratedZone(aP2Sfprm, aF1, aF2);
2830   if (bPInDZ) {
2831     return bPInDZ;
2832   }
2833
2834   const IntSurf_PntOn2S& aP2Slprm=aLineOn2S->Value(ilprm);
2835   bPInDZ=IsPointInDegeneratedZone(aP2Slprm, aF1, aF2);
2836   
2837   return bPInDZ;
2838 }
2839 //=======================================================================
2840 //function : IsPointInDegeneratedZone
2841 //purpose  : 
2842 //=======================================================================
2843 Standard_Boolean IsPointInDegeneratedZone(const IntSurf_PntOn2S& aP2S,
2844                                           const TopoDS_Face& aF1,
2845                                           const TopoDS_Face& aF2)
2846                                           
2847 {
2848   Standard_Boolean bFlag=Standard_True;
2849   Standard_Real US11, US12, VS11, VS12, US21, US22, VS21, VS22;
2850   Standard_Real U1, V1, U2, V2, aDelta, aD;
2851   gp_Pnt2d aP2d;
2852
2853   Handle(Geom_Surface)aS1 = BRep_Tool::Surface(aF1);
2854   aS1->Bounds(US11, US12, VS11, VS12);
2855   GeomAdaptor_Surface aGAS1(aS1);
2856
2857   Handle(Geom_Surface)aS2 = BRep_Tool::Surface(aF2);
2858   aS1->Bounds(US21, US22, VS21, VS22);
2859   GeomAdaptor_Surface aGAS2(aS2);
2860   //
2861   //const gp_Pnt& aP=aP2S.Value();
2862   aP2S.Parameters(U1, V1, U2, V2);
2863   //
2864   aDelta=1.e-7;
2865   // Check on Surf 1
2866   aD=aGAS1.UResolution(aDelta);
2867   aP2d.SetCoord(U1, V1);
2868   if (fabs(U1-US11) < aD) {
2869     bFlag=IsDegeneratedZone(aP2d, aS1, 1);
2870     if (bFlag) {
2871       return bFlag;
2872     }
2873   }
2874   if (fabs(U1-US12) < aD) {
2875     bFlag=IsDegeneratedZone(aP2d, aS1, 1);
2876     if (bFlag) {
2877       return bFlag;
2878     }
2879   }
2880   aD=aGAS1.VResolution(aDelta);
2881   if (fabs(V1-VS11) < aDelta) {
2882     bFlag=IsDegeneratedZone(aP2d, aS1, 2);
2883     if (bFlag) {
2884       return bFlag;
2885     }
2886   }
2887   if (fabs(V1-VS12) < aDelta) {
2888     bFlag=IsDegeneratedZone(aP2d, aS1, 2);
2889     if (bFlag) {
2890       return bFlag;
2891     }
2892   }
2893   // Check on Surf 2
2894   aD=aGAS2.UResolution(aDelta);
2895   aP2d.SetCoord(U2, V2);
2896   if (fabs(U2-US21) < aDelta) {
2897     bFlag=IsDegeneratedZone(aP2d, aS2, 1);
2898     if (bFlag) {
2899       return bFlag;
2900     }
2901   }
2902   if (fabs(U2-US22) < aDelta) {
2903     bFlag=IsDegeneratedZone(aP2d, aS2, 1);
2904     if (bFlag) {
2905       return bFlag;
2906     }
2907   }
2908   aD=aGAS2.VResolution(aDelta);
2909   if (fabs(V2-VS21) < aDelta) {
2910     bFlag=IsDegeneratedZone(aP2d, aS2, 2);
2911     if (bFlag) {  
2912       return bFlag;
2913     }
2914   }
2915   if (fabs(V2-VS22) < aDelta) {
2916     bFlag=IsDegeneratedZone(aP2d, aS2, 2);
2917     if (bFlag) {
2918       return bFlag;
2919     }
2920   }
2921   return !bFlag;
2922 }
2923
2924 //=======================================================================
2925 //function : IsDegeneratedZone
2926 //purpose  : 
2927 //=======================================================================
2928 Standard_Boolean IsDegeneratedZone(const gp_Pnt2d& aP2d,
2929                                    const Handle(Geom_Surface)& aS,
2930                                    const Standard_Integer iDir)
2931 {
2932   Standard_Boolean bFlag=Standard_True;
2933   Standard_Real US1, US2, VS1, VS2, dY, dX, d1, d2, dD;
2934   Standard_Real aXm, aYm, aXb, aYb, aXe, aYe;
2935   aS->Bounds(US1, US2, VS1, VS2); 
2936
2937   gp_Pnt aPm, aPb, aPe;
2938   
2939   aXm=aP2d.X();
2940   aYm=aP2d.Y();
2941   
2942   aS->D0(aXm, aYm, aPm); 
2943   
2944   dX=1.e-5;
2945   dY=1.e-5;
2946   dD=1.e-12;
2947
2948   if (iDir==1) {
2949     aXb=aXm;
2950     aXe=aXm;
2951     aYb=aYm-dY;
2952     if (aYb < VS1) {
2953       aYb=VS1;
2954     }
2955     aYe=aYm+dY;
2956     if (aYe > VS2) {
2957       aYe=VS2;
2958     }
2959     aS->D0(aXb, aYb, aPb);
2960     aS->D0(aXe, aYe, aPe);
2961     
2962     d1=aPm.Distance(aPb);
2963     d2=aPm.Distance(aPe);
2964     if (d1 < dD && d2 < dD) {
2965       return bFlag;
2966     }
2967     return !bFlag;
2968   }
2969   //
2970   else if (iDir==2) {
2971     aYb=aYm;
2972     aYe=aYm;
2973     aXb=aXm-dX;
2974     if (aXb < US1) {
2975       aXb=US1;
2976     }
2977     aXe=aXm+dX;
2978     if (aXe > US2) {
2979       aXe=US2;
2980     }
2981     aS->D0(aXb, aYb, aPb);
2982     aS->D0(aXe, aYe, aPe);
2983     
2984     d1=aPm.Distance(aPb);
2985     d2=aPm.Distance(aPe);
2986     if (d1 < dD && d2 < dD) {
2987       return bFlag;
2988     }
2989     return !bFlag;
2990   }
2991   return !bFlag;
2992 }
2993
2994 //=========================================================================
2995 // static function : ComputePurgedWLine
2996 // purpose : Removes equal points (leave one of equal points) from theWLine
2997 //           and recompute vertex parameters.
2998 //           Returns new WLine or null WLine if the number
2999 //           of the points is less than 2.
3000 //=========================================================================
3001 Handle(IntPatch_WLine) ComputePurgedWLine(const Handle(IntPatch_WLine)& theWLine) {
3002  
3003   Standard_Integer i, k, v, nb, nbvtx;
3004   Handle(IntPatch_WLine) aResult;
3005   nbvtx = theWLine->NbVertex();
3006   nb = theWLine->NbPnts();
3007   if (nb==2) {
3008     const IntSurf_PntOn2S& p1 = theWLine->Point(1);
3009     const IntSurf_PntOn2S& p2 = theWLine->Point(2);
3010     if(p1.Value().IsEqual(p2.Value(), gp::Resolution())) {
3011       return aResult;
3012     }
3013   }
3014   //
3015   Handle(IntPatch_WLine) aLocalWLine;
3016   Handle(IntPatch_WLine) aTmpWLine = theWLine;
3017   Handle(IntSurf_LineOn2S) aLineOn2S = new IntSurf_LineOn2S();
3018   aLocalWLine = new IntPatch_WLine(aLineOn2S, Standard_False);
3019   for(i = 1; i <= nb; i++) {
3020     aLineOn2S->Add(theWLine->Point(i));
3021   }
3022
3023   for(v = 1; v <= nbvtx; v++) {
3024     aLocalWLine->AddVertex(theWLine->Vertex(v));
3025   }
3026   
3027   for(i = 1; i <= aLineOn2S->NbPoints(); i++) {
3028     Standard_Integer aStartIndex = i + 1;
3029     Standard_Integer anEndIndex = i + 5;
3030     nb = aLineOn2S->NbPoints();
3031     anEndIndex = (anEndIndex > nb) ? nb : anEndIndex;
3032
3033     if((aStartIndex > nb) || (anEndIndex <= 1)) {
3034       continue;
3035     }
3036     k = aStartIndex;
3037
3038     while(k <= anEndIndex) {
3039       
3040       if(i != k) {
3041         IntSurf_PntOn2S p1 = aLineOn2S->Value(i);
3042         IntSurf_PntOn2S p2 = aLineOn2S->Value(k);
3043         
3044         if(p1.Value().IsEqual(p2.Value(), gp::Resolution())) {
3045           aTmpWLine = aLocalWLine;
3046           aLocalWLine = new IntPatch_WLine(aLineOn2S, Standard_False);
3047
3048           for(v = 1; v <= aTmpWLine->NbVertex(); v++) {
3049             IntPatch_Point aVertex = aTmpWLine->Vertex(v);
3050             Standard_Integer avertexindex = (Standard_Integer)aVertex.ParameterOnLine();
3051
3052             if(avertexindex >= k) {
3053               aVertex.SetParameter(aVertex.ParameterOnLine() - 1.);
3054             }
3055             aLocalWLine->AddVertex(aVertex);
3056           }
3057           aLineOn2S->RemovePoint(k);
3058           anEndIndex--;
3059           continue;
3060         }
3061       }
3062       k++;
3063     }
3064   }
3065
3066   if(aLineOn2S->NbPoints() > 1) {
3067     aResult = aLocalWLine;
3068   }
3069   return aResult;
3070 }
3071
3072 //=======================================================================
3073 //function : TolR3d
3074 //purpose  : 
3075 //=======================================================================
3076 void TolR3d(const TopoDS_Face& aF1,
3077             const TopoDS_Face& aF2,
3078             Standard_Real& myTolReached3d)
3079 {
3080   Standard_Real aTolF1, aTolF2, aTolFMax, aTolTresh;
3081       
3082   aTolTresh=2.999999e-3;
3083   aTolF1 = BRep_Tool::Tolerance(aF1);
3084   aTolF2 = BRep_Tool::Tolerance(aF2);
3085   aTolFMax=Max(aTolF1, aTolF2);
3086   
3087   if (aTolFMax>aTolTresh) {
3088     myTolReached3d=aTolFMax;
3089   }
3090 }
3091 //=======================================================================
3092 //function : AdjustPeriodic
3093 //purpose  : 
3094 //=======================================================================
3095 Standard_Real AdjustPeriodic(const Standard_Real theParameter,
3096                              const Standard_Real parmin,
3097                              const Standard_Real parmax,
3098                              const Standard_Real thePeriod,
3099                              Standard_Real&      theOffset) 
3100 {
3101   Standard_Real aresult;
3102   //
3103   theOffset = 0.;
3104   aresult = theParameter;
3105   while(aresult < parmin) {
3106     aresult += thePeriod;
3107     theOffset += thePeriod;
3108   }
3109
3110   while(aresult > parmax) {
3111     aresult -= thePeriod;
3112     theOffset -= thePeriod;
3113   }
3114   return aresult;
3115 }
3116 //=======================================================================
3117 //function : IsPointOnBoundary
3118 //purpose  : 
3119 //=======================================================================
3120 Standard_Boolean IsPointOnBoundary(const Standard_Real theParameter,
3121                                    const Standard_Real theFirstBoundary,
3122                                    const Standard_Real theSecondBoundary,
3123                                    const Standard_Real theResolution,
3124                                    Standard_Boolean&   IsOnFirstBoundary) 
3125 {
3126   Standard_Boolean bRet;
3127   Standard_Integer i;
3128   Standard_Real adist;
3129   //
3130   bRet=Standard_False;
3131   for(i = 0; i < 2; ++i) {
3132     IsOnFirstBoundary = (i == 0);
3133     if (IsOnFirstBoundary) {
3134       adist = fabs(theParameter - theFirstBoundary);
3135     }
3136     else {
3137       adist = fabs(theParameter - theSecondBoundary);
3138     }
3139     if(adist < theResolution) {
3140       return !bRet;
3141     }
3142   }
3143   return bRet;
3144 }
3145 // ------------------------------------------------------------------------------------------------
3146 // static function: FindPoint
3147 // purpose:
3148 // ------------------------------------------------------------------------------------------------
3149 Standard_Boolean FindPoint(const gp_Pnt2d&     theFirstPoint,
3150                            const gp_Pnt2d&     theLastPoint,
3151                            const Standard_Real theUmin, 
3152                            const Standard_Real theUmax,
3153                            const Standard_Real theVmin,
3154                            const Standard_Real theVmax,
3155                            gp_Pnt2d&           theNewPoint) {
3156   
3157   gp_Vec2d aVec(theFirstPoint, theLastPoint);
3158   Standard_Integer i = 0, j = 0;
3159
3160   for(i = 0; i < 4; i++) {
3161     gp_Vec2d anOtherVec;
3162     gp_Vec2d anOtherVecNormal;
3163     gp_Pnt2d aprojpoint = theLastPoint;    
3164
3165     if((i % 2) == 0) {
3166       anOtherVec.SetX(0.);
3167       anOtherVec.SetY(1.);
3168       anOtherVecNormal.SetX(1.);
3169       anOtherVecNormal.SetY(0.);
3170
3171       if(i < 2)
3172         aprojpoint.SetX(theUmin);
3173       else
3174         aprojpoint.SetX(theUmax);
3175     }
3176     else {
3177       anOtherVec.SetX(1.);
3178       anOtherVec.SetY(0.);
3179       anOtherVecNormal.SetX(0.);
3180       anOtherVecNormal.SetY(1.);
3181
3182       if(i < 2)
3183         aprojpoint.SetY(theVmin);
3184       else
3185         aprojpoint.SetY(theVmax);
3186     }
3187     gp_Vec2d anormvec = aVec;
3188     anormvec.Normalize();
3189     RefineVector(anormvec);
3190     Standard_Real adot1 = anormvec.Dot(anOtherVecNormal);
3191
3192     if(fabs(adot1) < Precision::Angular())
3193       continue;
3194     Standard_Real adist = 0.;
3195     Standard_Boolean bIsOut = Standard_False;
3196
3197     if((i % 2) == 0) {
3198       adist = (i < 2) ? fabs(theLastPoint.X() - theUmin) : fabs(theLastPoint.X() - theUmax);
3199       bIsOut = (i < 2) ? (theLastPoint.X() < theUmin) : (theLastPoint.X() > theUmax);
3200     }
3201     else {
3202       adist = (i < 2) ? fabs(theLastPoint.Y() - theVmin) : fabs(theLastPoint.Y() - theVmax);
3203       bIsOut = (i < 2) ? (theLastPoint.Y() < theVmin) : (theLastPoint.Y() > theVmax);
3204     }
3205     Standard_Real anoffset = adist * anOtherVec.Dot(anormvec) / adot1;
3206
3207     for(j = 0; j < 2; j++) {
3208       anoffset = (j == 0) ? anoffset : -anoffset;
3209       gp_Pnt2d acurpoint(aprojpoint.XY() + (anOtherVec.XY()*anoffset));
3210       gp_Vec2d acurvec(theLastPoint, acurpoint);
3211       if ( bIsOut )
3212         acurvec.Reverse();
3213
3214       Standard_Real aDotX, anAngleX;
3215       //
3216       aDotX = aVec.Dot(acurvec);
3217       anAngleX = aVec.Angle(acurvec);
3218       //
3219       if(aDotX > 0. && fabs(anAngleX) < Precision::PConfusion()) {
3220         if((i % 2) == 0) {
3221           if((acurpoint.Y() >= theVmin) &&
3222              (acurpoint.Y() <= theVmax)) {
3223             theNewPoint = acurpoint;
3224             return Standard_True;
3225           }
3226         }
3227         else {
3228           if((acurpoint.X() >= theUmin) &&
3229              (acurpoint.X() <= theUmax)) {
3230             theNewPoint = acurpoint;
3231             return Standard_True;
3232           }
3233         }
3234       }
3235     }
3236   }
3237   return Standard_False;
3238 }
3239
3240
3241 // ------------------------------------------------------------------------------------------------
3242 // static function: FindPoint
3243 // purpose: Find point on the boundary of radial tangent zone
3244 // ------------------------------------------------------------------------------------------------
3245 Standard_Boolean FindPoint(const gp_Pnt2d&     theFirstPoint,
3246                            const gp_Pnt2d&     theLastPoint,
3247                            const Standard_Real theUmin, 
3248                            const Standard_Real theUmax,
3249                            const Standard_Real theVmin,
3250                            const Standard_Real theVmax,
3251                            const gp_Pnt2d&     theTanZoneCenter,
3252                            const Standard_Real theZoneRadius,
3253                            Handle(GeomAdaptor_HSurface) theGASurface,
3254                            gp_Pnt2d&           theNewPoint) {
3255   theNewPoint = theLastPoint;
3256
3257   if ( !IsInsideTanZone( theLastPoint, theTanZoneCenter, theZoneRadius, theGASurface) )
3258     return Standard_False;
3259
3260   Standard_Real aUResolution = theGASurface->UResolution( theZoneRadius );
3261   Standard_Real aVResolution = theGASurface->VResolution( theZoneRadius );
3262
3263   Standard_Real aRadius = ( aUResolution < aVResolution ) ? aUResolution : aVResolution;
3264   gp_Ax22d anAxis( theTanZoneCenter, gp_Dir2d(1, 0), gp_Dir2d(0, 1) );
3265   gp_Circ2d aCircle( anAxis, aRadius );
3266   
3267   //
3268   gp_Vec2d aDir( theLastPoint.XY() - theFirstPoint.XY() );
3269   Standard_Real aLength = aDir.Magnitude();
3270   if ( aLength <= gp::Resolution() )
3271     return Standard_False;
3272   gp_Lin2d aLine( theFirstPoint, aDir );
3273
3274   //
3275   Handle(Geom2d_Line) aCLine = new Geom2d_Line( aLine );
3276   Handle(Geom2d_TrimmedCurve) aC1 = new Geom2d_TrimmedCurve( aCLine, 0, aLength );
3277   Handle(Geom2d_Circle) aC2 = new Geom2d_Circle( aCircle );
3278
3279   Standard_Real aTol = aRadius * 0.001;
3280   aTol = ( aTol < Precision::PConfusion() ) ? Precision::PConfusion() : aTol;
3281
3282   Geom2dAPI_InterCurveCurve anIntersector;
3283   anIntersector.Init( aC1, aC2, aTol );
3284
3285   if ( anIntersector.NbPoints() == 0 )
3286     return Standard_False;
3287
3288   Standard_Boolean aFound = Standard_False;
3289   Standard_Real aMinDist = aLength * aLength;
3290   Standard_Integer i = 0;
3291   for ( i = 1; i <= anIntersector.NbPoints(); i++ ) {
3292     gp_Pnt2d aPInt = anIntersector.Point( i );
3293     if ( aPInt.SquareDistance( theFirstPoint ) < aMinDist ) {
3294       if ( ( aPInt.X() >= theUmin ) && ( aPInt.X() <= theUmax ) &&
3295            ( aPInt.Y() >= theVmin ) && ( aPInt.Y() <= theVmax ) ) {
3296         theNewPoint = aPInt;
3297         aFound = Standard_True;
3298       }
3299     }
3300   }
3301
3302   return aFound;
3303 }
3304
3305 // ------------------------------------------------------------------------------------------------
3306 // static function: IsInsideTanZone
3307 // purpose: Check if point is inside a radial tangent zone
3308 // ------------------------------------------------------------------------------------------------
3309 Standard_Boolean IsInsideTanZone(const gp_Pnt2d&     thePoint,
3310                                  const gp_Pnt2d&     theTanZoneCenter,
3311                                  const Standard_Real theZoneRadius,
3312                                  Handle(GeomAdaptor_HSurface) theGASurface) {
3313
3314   Standard_Real aUResolution = theGASurface->UResolution( theZoneRadius );
3315   Standard_Real aVResolution = theGASurface->VResolution( theZoneRadius );
3316   Standard_Real aRadiusSQR = ( aUResolution < aVResolution ) ? aUResolution : aVResolution;
3317   aRadiusSQR *= aRadiusSQR;
3318   if ( thePoint.SquareDistance( theTanZoneCenter ) <= aRadiusSQR )
3319     return Standard_True;
3320   return Standard_False;
3321 }
3322
3323 // ------------------------------------------------------------------------------------------------
3324 // static function: CheckTangentZonesExist
3325 // purpose: Check if tangent zone exists
3326 // ------------------------------------------------------------------------------------------------
3327 Standard_Boolean CheckTangentZonesExist( const Handle(GeomAdaptor_HSurface)& theSurface1,
3328                                         const Handle(GeomAdaptor_HSurface)&  theSurface2 ) 
3329 {
3330   if ( ( theSurface1->GetType() != GeomAbs_Torus ) ||
3331       ( theSurface2->GetType() != GeomAbs_Torus ) )
3332     return Standard_False;
3333
3334   gp_Torus aTor1 = theSurface1->Torus();
3335   gp_Torus aTor2 = theSurface2->Torus();
3336
3337   if ( aTor1.Location().Distance( aTor2.Location() ) > Precision::Confusion() )
3338     return Standard_False;
3339
3340   if ( ( fabs( aTor1.MajorRadius() - aTor2.MajorRadius() ) > Precision::Confusion() ) ||
3341        ( fabs( aTor1.MinorRadius() - aTor2.MinorRadius() ) > Precision::Confusion() ) )
3342     return Standard_False;
3343
3344   if ( ( aTor1.MajorRadius() < aTor1.MinorRadius() ) ||
3345        ( aTor2.MajorRadius() < aTor2.MinorRadius() ) )
3346     return Standard_False;
3347   return Standard_True;
3348 }
3349
3350 // ------------------------------------------------------------------------------------------------
3351 // static function: ComputeTangentZones
3352 // purpose: 
3353 // ------------------------------------------------------------------------------------------------
3354 Standard_Integer ComputeTangentZones( const Handle(GeomAdaptor_HSurface)& theSurface1,
3355                                      const Handle(GeomAdaptor_HSurface)&  theSurface2,
3356                                      const TopoDS_Face&                   theFace1,
3357                                      const TopoDS_Face&                   theFace2,
3358                                      Handle(TColgp_HArray1OfPnt2d)&       theResultOnS1,
3359                                      Handle(TColgp_HArray1OfPnt2d)&       theResultOnS2,
3360                                      Handle(TColStd_HArray1OfReal)&       theResultRadius,
3361                                      const Handle(BOPInt_Context)& aContext)
3362 {
3363   Standard_Integer aResult = 0;
3364   if ( !CheckTangentZonesExist( theSurface1, theSurface2 ) )
3365     return aResult;
3366
3367
3368   TColgp_SequenceOfPnt2d aSeqResultS1, aSeqResultS2;
3369   TColStd_SequenceOfReal aSeqResultRad;
3370
3371   gp_Torus aTor1 = theSurface1->Torus();
3372   gp_Torus aTor2 = theSurface2->Torus();
3373
3374   gp_Ax2 anax1( aTor1.Location(), aTor1.Axis().Direction() );
3375   gp_Ax2 anax2( aTor2.Location(), aTor2.Axis().Direction() );
3376   Standard_Integer j = 0;
3377
3378   for ( j = 0; j < 2; j++ ) {
3379     Standard_Real aCoef = ( j == 0 ) ? -1 : 1;
3380     Standard_Real aRadius1 = fabs(aTor1.MajorRadius() + aCoef * aTor1.MinorRadius());
3381     Standard_Real aRadius2 = fabs(aTor2.MajorRadius() + aCoef * aTor2.MinorRadius());
3382
3383     gp_Circ aCircle1( anax1, aRadius1 );
3384     gp_Circ aCircle2( anax2, aRadius2 );
3385
3386     // roughly compute radius of tangent zone for perpendicular case
3387     Standard_Real aCriteria = Precision::Confusion() * 0.5;
3388
3389     Standard_Real aT1 = aCriteria;
3390     Standard_Real aT2 = aCriteria;
3391     if ( j == 0 ) {
3392       // internal tangency
3393       Standard_Real aR = ( aRadius1 > aTor2.MinorRadius() ) ? aRadius1 : aTor2.MinorRadius();
3394       //aT1 = aCriteria * aCriteria + aR * aR - ( aR - aCriteria ) * ( aR - aCriteria );
3395       aT1 = 2. * aR * aCriteria;
3396       aT2 = aT1;
3397     }
3398     else {
3399       // external tangency
3400       Standard_Real aRb = ( aRadius1 > aTor2.MinorRadius() ) ? aRadius1 : aTor2.MinorRadius();
3401       Standard_Real aRm = ( aRadius1 < aTor2.MinorRadius() ) ? aRadius1 : aTor2.MinorRadius();
3402       Standard_Real aDelta = aRb - aCriteria;
3403       aDelta *= aDelta;
3404       aDelta -= aRm * aRm;
3405       aDelta /= 2. * (aRb - aRm);
3406       aDelta -= 0.5 * (aRb - aRm);
3407       
3408       aT1 = 2. * aRm * (aRm - aDelta);
3409       aT2 = aT1;
3410     }
3411     aCriteria = ( aT1 > aT2) ? aT1 : aT2;
3412     if ( aCriteria > 0 )
3413       aCriteria = sqrt( aCriteria );
3414
3415     if ( aCriteria > 0.5 * aTor1.MinorRadius() ) {
3416       // too big zone -> drop to minimum
3417       aCriteria = Precision::Confusion();
3418     }
3419
3420     GeomAdaptor_Curve aC1( new Geom_Circle(aCircle1) );
3421     GeomAdaptor_Curve aC2( new Geom_Circle(aCircle2) );
3422     Extrema_ExtCC anExtrema(aC1, aC2, 0, 2. * M_PI, 0, 2. * M_PI, 
3423                             Precision::PConfusion(), Precision::PConfusion());
3424         
3425     if ( anExtrema.IsDone() ) {
3426
3427       Standard_Integer i = 0;
3428       for ( i = 1; i <= anExtrema.NbExt(); i++ ) {
3429         if ( anExtrema.SquareDistance(i) > aCriteria * aCriteria )
3430           continue;
3431
3432         Extrema_POnCurv P1, P2;
3433         anExtrema.Points( i, P1, P2 );
3434
3435         Standard_Boolean bFoundResult = Standard_True;
3436         gp_Pnt2d pr1, pr2;
3437
3438         Standard_Integer surfit = 0;
3439         for ( surfit = 0; surfit < 2; surfit++ ) {
3440           GeomAPI_ProjectPointOnSurf& aProjector = 
3441             (surfit == 0) ? aContext->ProjPS(theFace1) : aContext->ProjPS(theFace2);
3442
3443           gp_Pnt aP3d = (surfit == 0) ? P1.Value() : P2.Value();
3444           aProjector.Perform(aP3d);
3445
3446           if(!aProjector.IsDone())
3447             bFoundResult = Standard_False;
3448           else {
3449             if(aProjector.LowerDistance() > aCriteria) {
3450               bFoundResult = Standard_False;
3451             }
3452             else {
3453               Standard_Real foundU = 0, foundV = 0;
3454               aProjector.LowerDistanceParameters(foundU, foundV);
3455               if ( surfit == 0 )
3456                 pr1 = gp_Pnt2d( foundU, foundV );
3457               else
3458                 pr2 = gp_Pnt2d( foundU, foundV );
3459             }
3460           }
3461         }
3462         if ( bFoundResult ) {
3463           aSeqResultS1.Append( pr1 );
3464           aSeqResultS2.Append( pr2 );
3465           aSeqResultRad.Append( aCriteria );
3466
3467           // torus is u and v periodic
3468           const Standard_Real twoPI = M_PI + M_PI;
3469           Standard_Real arr1tmp[2] = {pr1.X(), pr1.Y()};
3470           Standard_Real arr2tmp[2] = {pr2.X(), pr2.Y()};
3471
3472           // iteration on period bounds
3473           for ( Standard_Integer k1 = 0; k1 < 2; k1++ ) {
3474             Standard_Real aBound = ( k1 == 0 ) ? 0 : twoPI;
3475             Standard_Real aShift = ( k1 == 0 ) ? twoPI : -twoPI;
3476
3477             // iteration on surfaces
3478             for ( Standard_Integer k2 = 0; k2 < 2; k2++ ) {
3479               Standard_Real* arr1 = ( k2 == 0 ) ? arr1tmp : arr2tmp;
3480               Standard_Real* arr2 = ( k2 != 0 ) ? arr1tmp : arr2tmp;
3481               TColgp_SequenceOfPnt2d& aSeqS1 = ( k2 == 0 ) ? aSeqResultS1 : aSeqResultS2; 
3482               TColgp_SequenceOfPnt2d& aSeqS2 = ( k2 != 0 ) ? aSeqResultS1 : aSeqResultS2; 
3483
3484               if (fabs(arr1[0] - aBound) < Precision::PConfusion()) {
3485                 aSeqS1.Append( gp_Pnt2d( arr1[0] + aShift, arr1[1] ) );
3486                 aSeqS2.Append( gp_Pnt2d( arr2[0], arr2[1] ) );
3487                 aSeqResultRad.Append( aCriteria );
3488               }
3489               if (fabs(arr1[1] - aBound) < Precision::PConfusion()) {
3490                 aSeqS1.Append( gp_Pnt2d( arr1[0], arr1[1] + aShift) );
3491                 aSeqS2.Append( gp_Pnt2d( arr2[0], arr2[1] ) );
3492                 aSeqResultRad.Append( aCriteria );
3493               }
3494             }
3495           } //
3496         }
3497       }
3498     }
3499   }
3500   aResult = aSeqResultRad.Length();
3501
3502   if ( aResult > 0 ) {
3503     theResultOnS1 = new TColgp_HArray1OfPnt2d( 1, aResult );
3504     theResultOnS2 = new TColgp_HArray1OfPnt2d( 1, aResult );
3505     theResultRadius = new TColStd_HArray1OfReal( 1, aResult );
3506
3507     for ( Standard_Integer i = 1 ; i <= aResult; i++ ) {
3508       theResultOnS1->SetValue( i, aSeqResultS1.Value(i) );
3509       theResultOnS2->SetValue( i, aSeqResultS2.Value(i) );
3510       theResultRadius->SetValue( i, aSeqResultRad.Value(i) );
3511     }
3512   }
3513   return aResult;
3514 }
3515
3516 // ------------------------------------------------------------------------------------------------
3517 // static function: AdjustByNeighbour
3518 // purpose:
3519 // ------------------------------------------------------------------------------------------------
3520 gp_Pnt2d AdjustByNeighbour(const gp_Pnt2d&     theaNeighbourPoint,
3521                            const gp_Pnt2d&     theOriginalPoint,
3522                            Handle(GeomAdaptor_HSurface) theGASurface) {
3523   
3524   gp_Pnt2d ap1 = theaNeighbourPoint;
3525   gp_Pnt2d ap2 = theOriginalPoint;
3526
3527   if ( theGASurface->IsUPeriodic() ) {
3528     Standard_Real aPeriod     = theGASurface->UPeriod();
3529     gp_Pnt2d aPTest = ap2;
3530     Standard_Real aSqDistMin = 1.e+100;
3531
3532     for ( Standard_Integer pIt = -1; pIt <= 1; pIt++) {
3533       aPTest.SetX( theOriginalPoint.X() + aPeriod * pIt );
3534       Standard_Real dd = ap1.SquareDistance( aPTest );
3535
3536       if ( dd < aSqDistMin ) {
3537         ap2 = aPTest;
3538         aSqDistMin = dd;
3539       }
3540     }
3541   }
3542   if ( theGASurface->IsVPeriodic() ) {
3543     Standard_Real aPeriod     = theGASurface->VPeriod();
3544     gp_Pnt2d aPTest = ap2;
3545     Standard_Real aSqDistMin = 1.e+100;
3546
3547     for ( Standard_Integer pIt = -1; pIt <= 1; pIt++) {
3548       aPTest.SetY( theOriginalPoint.Y() + aPeriod * pIt );
3549       Standard_Real dd = ap1.SquareDistance( aPTest );
3550
3551       if ( dd < aSqDistMin ) {
3552         ap2 = aPTest;
3553         aSqDistMin = dd;
3554       }
3555     }
3556   }
3557   return ap2;
3558 }
3559
3560 // ------------------------------------------------------------------------------------------------
3561 //function: DecompositionOfWLine
3562 // purpose:
3563 // ------------------------------------------------------------------------------------------------
3564 Standard_Boolean DecompositionOfWLine(const Handle(IntPatch_WLine)& theWLine,
3565                                       const Handle(GeomAdaptor_HSurface)&            theSurface1, 
3566                                       const Handle(GeomAdaptor_HSurface)&            theSurface2,
3567                                       const TopoDS_Face&                             theFace1,
3568                                       const TopoDS_Face&                             theFace2,
3569                                       const IntTools_LineConstructor&                theLConstructor,
3570                                       const Standard_Boolean                         theAvoidLConstructor,
3571                                       IntPatch_SequenceOfLine&                       theNewLines,
3572                                       Standard_Real&                                 theReachedTol3d,
3573                                       const Handle(BOPInt_Context)& aContext) 
3574 {
3575
3576   Standard_Boolean bRet, bAvoidLineConstructor;
3577   Standard_Integer aNbPnts, aNbParts;
3578   //
3579   bRet=Standard_False;
3580   aNbPnts=theWLine->NbPnts();
3581   bAvoidLineConstructor=theAvoidLConstructor;
3582   //
3583   if(!aNbPnts){
3584     return bRet;
3585   }
3586   if (!bAvoidLineConstructor) {
3587     aNbParts=theLConstructor.NbParts();
3588     if (!aNbParts) {
3589       return bRet;
3590     }
3591   }
3592   //
3593   Standard_Boolean bIsPrevPointOnBoundary, bIsPointOnBoundary, bIsCurrentPointOnBoundary;
3594   Standard_Integer nblines, pit, i, j;
3595   Standard_Real aTol;
3596   TColStd_Array1OfListOfInteger anArrayOfLines(1, aNbPnts); 
3597   TColStd_Array1OfInteger       anArrayOfLineType(1, aNbPnts);
3598   TColStd_ListOfInteger aListOfPointIndex;
3599   
3600   Handle(TColgp_HArray1OfPnt2d) aTanZoneS1;
3601   Handle(TColgp_HArray1OfPnt2d) aTanZoneS2;
3602   Handle(TColStd_HArray1OfReal) aTanZoneRadius;
3603   Standard_Integer aNbZone = ComputeTangentZones( theSurface1, theSurface2, theFace1, theFace2,
3604                                                  aTanZoneS1, aTanZoneS2, aTanZoneRadius, aContext);
3605   
3606   //
3607   nblines=0;
3608   aTol=Precision::Confusion();
3609   aTol=0.5*aTol;
3610   bIsPrevPointOnBoundary=Standard_False;
3611   bIsPointOnBoundary=Standard_False;
3612   //
3613   // 1. ...
3614   //
3615   // Points
3616   for(pit = 1; pit <= aNbPnts; ++pit) {
3617     Standard_Boolean bIsOnFirstBoundary, isperiodic;
3618     Standard_Real aResolution, aPeriod, alowerboundary, aupperboundary, U, V;
3619     Standard_Real aParameter, anoffset, anAdjustPar;
3620     Standard_Real umin, umax, vmin, vmax;
3621     //
3622     bIsCurrentPointOnBoundary = Standard_False;
3623     const IntSurf_PntOn2S& aPoint = theWLine->Point(pit);
3624     //
3625     // Surface
3626     for(i = 0; i < 2; ++i) {
3627       Handle(GeomAdaptor_HSurface) aGASurface = (!i) ? theSurface1 : theSurface2;
3628       aGASurface->ChangeSurface().Surface()->Bounds(umin, umax, vmin, vmax);
3629       if(!i) {
3630         aPoint.ParametersOnS1(U, V);
3631       }
3632       else {
3633         aPoint.ParametersOnS2(U, V);
3634       }
3635       // U, V
3636       for(j = 0; j < 2; j++) {
3637         isperiodic = (!j) ? aGASurface->IsUPeriodic() : aGASurface->IsVPeriodic();
3638         if(!isperiodic){
3639           continue;
3640         }
3641         //
3642         if (!j) {
3643           aResolution=aGASurface->UResolution(aTol);
3644           aPeriod=aGASurface->UPeriod();
3645           alowerboundary=umin;
3646           aupperboundary=umax;
3647           aParameter=U;
3648         }
3649         else {
3650           aResolution=aGASurface->VResolution(aTol);
3651           aPeriod=aGASurface->VPeriod();
3652           alowerboundary=vmin;
3653           aupperboundary=vmax;
3654           aParameter=V;
3655         }
3656         
3657         anoffset = 0.;
3658         anAdjustPar = AdjustPeriodic(aParameter, 
3659                                      alowerboundary, 
3660                                      aupperboundary, 
3661                                      aPeriod, 
3662                                      anoffset);
3663         //
3664         bIsOnFirstBoundary = Standard_True;// ?
3665         bIsPointOnBoundary=
3666           IsPointOnBoundary(anAdjustPar, 
3667                             alowerboundary, 
3668                             aupperboundary,
3669                             aResolution, 
3670                             bIsOnFirstBoundary);
3671         //
3672         if(bIsPointOnBoundary) {
3673           bIsCurrentPointOnBoundary = Standard_True;
3674           break;
3675         }
3676         else {
3677           // check if a point belong to a tangent zone. Begin
3678           Standard_Integer zIt = 0;
3679           for ( zIt = 1; zIt <= aNbZone; zIt++ ) {
3680             gp_Pnt2d aPZone = (i == 0) ? aTanZoneS1->Value(zIt) : aTanZoneS2->Value(zIt);
3681             Standard_Real aZoneRadius = aTanZoneRadius->Value(zIt);
3682
3683             if ( IsInsideTanZone(gp_Pnt2d( U, V ), aPZone, aZoneRadius, aGASurface ) ) {
3684               // set boundary flag to split the curve by a tangent zone
3685               bIsPointOnBoundary = Standard_True;
3686               bIsCurrentPointOnBoundary = Standard_True;
3687               if ( theReachedTol3d < aZoneRadius ) {
3688                 theReachedTol3d = aZoneRadius;
3689               }
3690               break;
3691             }
3692           }
3693         }
3694       }//for(j = 0; j < 2; j++) {
3695
3696       if(bIsCurrentPointOnBoundary){
3697         break;
3698       }
3699     }//for(i = 0; i < 2; ++i) {
3700     //
3701     if((bIsCurrentPointOnBoundary != bIsPrevPointOnBoundary)) {
3702       if(!aListOfPointIndex.IsEmpty()) {
3703         nblines++;
3704         anArrayOfLines.SetValue(nblines, aListOfPointIndex);
3705         anArrayOfLineType.SetValue(nblines, bIsPrevPointOnBoundary);
3706         aListOfPointIndex.Clear();
3707       }
3708       bIsPrevPointOnBoundary = bIsCurrentPointOnBoundary;
3709     }
3710     aListOfPointIndex.Append(pit);
3711   } //for(pit = 1; pit <= aNbPnts; ++pit) {
3712   //
3713   if(!aListOfPointIndex.IsEmpty()) {
3714     nblines++;
3715     anArrayOfLines.SetValue(nblines, aListOfPointIndex);
3716     anArrayOfLineType.SetValue(nblines, bIsPrevPointOnBoundary);
3717     aListOfPointIndex.Clear();
3718   }
3719   //
3720   if(nblines<=1) {
3721     return bRet; //Standard_False;
3722   }
3723   //
3724   // 
3725   // 2. Correct wlines.begin
3726   TColStd_Array1OfListOfInteger anArrayOfLineEnds(1, nblines);
3727   Handle(IntSurf_LineOn2S) aSeqOfPntOn2S = new IntSurf_LineOn2S();
3728   //
3729   for(i = 1; i <= nblines; i++) {
3730     if(anArrayOfLineType.Value(i) != 0) {
3731       continue;
3732     }
3733     const TColStd_ListOfInteger& aListOfIndex = anArrayOfLines.Value(i);
3734     if(aListOfIndex.Extent() < 2) {
3735       continue;
3736     }
3737     TColStd_ListOfInteger aListOfFLIndex;
3738
3739     for(j = 0; j < 2; j++) {
3740       Standard_Integer aneighbourindex = (j == 0) ? (i - 1) : (i + 1);
3741
3742       if((aneighbourindex < 1) || (aneighbourindex > nblines))
3743         continue;
3744
3745       if(anArrayOfLineType.Value(aneighbourindex) == 0)
3746         continue;
3747       const TColStd_ListOfInteger& aNeighbour = anArrayOfLines.Value(aneighbourindex);
3748       Standard_Integer anIndex = (j == 0) ? aNeighbour.Last() : aNeighbour.First();
3749       const IntSurf_PntOn2S& aPoint = theWLine->Point(anIndex);
3750
3751       IntSurf_PntOn2S aNewP = aPoint;
3752       
3753       for(Standard_Integer surfit = 0; surfit < 2; surfit++) {
3754
3755         Handle(GeomAdaptor_HSurface) aGASurface = (surfit == 0) ? theSurface1 : theSurface2;
3756         Standard_Real umin=0., umax=0., vmin=0., vmax=0.;
3757         aGASurface->ChangeSurface().Surface()->Bounds(umin, umax, vmin, vmax);
3758         Standard_Real U=0., V=0.;
3759
3760         if(surfit == 0)
3761           aNewP.ParametersOnS1(U, V);
3762         else
3763           aNewP.ParametersOnS2(U, V);
3764         Standard_Integer nbboundaries = 0;
3765
3766         Standard_Boolean bIsNearBoundary = Standard_False;
3767         Standard_Integer aZoneIndex = 0;
3768         Standard_Integer bIsUBoundary = Standard_False; // use if nbboundaries == 1
3769         Standard_Integer bIsFirstBoundary = Standard_False; // use if nbboundaries == 1
3770         
3771
3772         for(Standard_Integer parit = 0; parit < 2; parit++) {
3773           Standard_Boolean isperiodic = (parit == 0) ? aGASurface->IsUPeriodic() : aGASurface->IsVPeriodic();
3774
3775           Standard_Real aResolution = (parit == 0) ? aGASurface->UResolution(aTol) : aGASurface->VResolution(aTol);
3776           Standard_Real alowerboundary = (parit == 0) ? umin : vmin;
3777           Standard_Real aupperboundary = (parit == 0) ? umax : vmax;
3778
3779           Standard_Real aParameter = (parit == 0) ? U : V;
3780           Standard_Boolean bIsOnFirstBoundary = Standard_True;
3781   
3782           if(!isperiodic) {
3783             bIsPointOnBoundary=
3784               IsPointOnBoundary(aParameter, alowerboundary, aupperboundary, aResolution, bIsOnFirstBoundary);
3785             if(bIsPointOnBoundary) {
3786               bIsUBoundary = (parit == 0);
3787               bIsFirstBoundary = bIsOnFirstBoundary;
3788               nbboundaries++;
3789             }
3790           }
3791           else {
3792             Standard_Real aPeriod     = (parit == 0) ? aGASurface->UPeriod() : aGASurface->VPeriod();
3793             Standard_Real anoffset = 0.;
3794             Standard_Real anAdjustPar = AdjustPeriodic(aParameter, alowerboundary, aupperboundary, aPeriod, anoffset);
3795
3796             bIsPointOnBoundary=
3797               IsPointOnBoundary(anAdjustPar, alowerboundary, aupperboundary, aResolution, bIsOnFirstBoundary);
3798             if(bIsPointOnBoundary) {
3799               bIsUBoundary = (parit == 0);
3800               bIsFirstBoundary = bIsOnFirstBoundary;
3801               nbboundaries++;
3802             }
3803             else {
3804               //check neighbourhood of boundary
3805               Standard_Real anEpsilon = aResolution * 100.;
3806               Standard_Real aPart = ( aupperboundary - alowerboundary ) * 0.1;
3807               anEpsilon = ( anEpsilon > aPart ) ? aPart : anEpsilon;
3808                 
3809               bIsNearBoundary = IsPointOnBoundary(anAdjustPar, alowerboundary, aupperboundary, 
3810                                                   anEpsilon, bIsOnFirstBoundary);
3811
3812             }
3813           }
3814         }
3815
3816         // check if a point belong to a tangent zone. Begin
3817         for ( Standard_Integer zIt = 1; zIt <= aNbZone; zIt++ ) {
3818           gp_Pnt2d aPZone = (surfit == 0) ? aTanZoneS1->Value(zIt) : aTanZoneS2->Value(zIt);
3819           Standard_Real aZoneRadius = aTanZoneRadius->Value(zIt);
3820
3821           Standard_Integer aneighbourpointindex1 = (j == 0) ? aListOfIndex.First() : aListOfIndex.Last();
3822           const IntSurf_PntOn2S& aNeighbourPoint = theWLine->Point(aneighbourpointindex1);
3823           Standard_Real nU1, nV1;
3824             
3825           if(surfit == 0)
3826             aNeighbourPoint.ParametersOnS1(nU1, nV1);
3827           else
3828             aNeighbourPoint.ParametersOnS2(nU1, nV1);
3829           gp_Pnt2d ap1(nU1, nV1);
3830           gp_Pnt2d ap2 = AdjustByNeighbour( ap1, gp_Pnt2d( U, V ), aGASurface );
3831
3832
3833           if ( IsInsideTanZone( ap2, aPZone, aZoneRadius, aGASurface ) ) {
3834             aZoneIndex = zIt;
3835             bIsNearBoundary = Standard_True;
3836             if ( theReachedTol3d < aZoneRadius ) {
3837               theReachedTol3d = aZoneRadius;
3838             }
3839           }
3840         }
3841         // check if a point belong to a tangent zone. End
3842         Standard_Boolean bComputeLineEnd = Standard_False;
3843
3844         if(nbboundaries == 2) {
3845           //xf
3846           bComputeLineEnd = Standard_True;
3847           //xt
3848         }
3849         else if(nbboundaries == 1) {
3850           Standard_Boolean isperiodic = (bIsUBoundary) ? aGASurface->IsUPeriodic() : aGASurface->IsVPeriodic();
3851
3852           if(isperiodic) {
3853             Standard_Real alowerboundary = (bIsUBoundary) ? umin : vmin;
3854             Standard_Real aupperboundary = (bIsUBoundary) ? umax : vmax;
3855             Standard_Real aPeriod     = (bIsUBoundary) ? aGASurface->UPeriod() : aGASurface->VPeriod();
3856             Standard_Real aParameter = (bIsUBoundary) ? U : V;
3857             Standard_Real anoffset = 0.;
3858             Standard_Real anAdjustPar = AdjustPeriodic(aParameter, alowerboundary, aupperboundary, aPeriod, anoffset);
3859
3860             Standard_Real adist = (bIsFirstBoundary) ? fabs(anAdjustPar - alowerboundary) : fabs(anAdjustPar - aupperboundary);
3861             Standard_Real anotherPar = (bIsFirstBoundary) ? (aupperboundary - adist) : (alowerboundary + adist);
3862             anotherPar += anoffset;
3863             Standard_Integer aneighbourpointindex = (j == 0) ? aListOfIndex.First() : aListOfIndex.Last();
3864             const IntSurf_PntOn2S& aNeighbourPoint = theWLine->Point(aneighbourpointindex);
3865             Standard_Real nU1, nV1;
3866
3867             if(surfit == 0)
3868               aNeighbourPoint.ParametersOnS1(nU1, nV1);
3869             else
3870               aNeighbourPoint.ParametersOnS2(nU1, nV1);
3871             
3872             Standard_Real adist1 = (bIsUBoundary) ? fabs(nU1 - U) : fabs(nV1 - V);
3873             Standard_Real adist2 = (bIsUBoundary) ? fabs(nU1 - anotherPar) : fabs(nV1 - anotherPar);
3874             bComputeLineEnd = Standard_True;
3875             Standard_Boolean bCheckAngle1 = Standard_False;
3876             Standard_Boolean bCheckAngle2 = Standard_False;
3877             gp_Vec2d aNewVec;
3878             Standard_Real anewU = (bIsUBoundary) ? anotherPar : U;
3879             Standard_Real anewV = (bIsUBoundary) ? V : anotherPar;
3880
3881             if(((adist1 - adist2) > Precision::PConfusion()) && 
3882                (adist2 < (aPeriod / 4.))) {
3883               bCheckAngle1 = Standard_True;
3884               aNewVec = gp_Vec2d(gp_Pnt2d(nU1, nV1), gp_Pnt2d(anewU, anewV));
3885
3886               if(aNewVec.SquareMagnitude() < (gp::Resolution() * gp::Resolution())) {
3887                 aNewP.SetValue((surfit == 0), anewU, anewV);
3888                 bCheckAngle1 = Standard_False;
3889               }
3890             }
3891             else if(adist1 < (aPeriod / 4.)) {
3892               bCheckAngle2 = Standard_True;
3893               aNewVec = gp_Vec2d(gp_Pnt2d(nU1, nV1), gp_Pnt2d(U, V));
3894
3895               if(aNewVec.SquareMagnitude() < (gp::Resolution() * gp::Resolution())) {
3896                 bCheckAngle2 = Standard_False;
3897               }
3898             }
3899
3900             if(bCheckAngle1 || bCheckAngle2) {
3901               // assume there are at least two points in line (see "if" above)
3902               Standard_Integer anindexother = aneighbourpointindex;
3903
3904               while((anindexother <= aListOfIndex.Last()) && (anindexother >= aListOfIndex.First())) {
3905                 anindexother = (j == 0) ? (anindexother + 1) : (anindexother - 1);
3906                 const IntSurf_PntOn2S& aPrevNeighbourPoint = theWLine->Point(anindexother);
3907                 Standard_Real nU2, nV2;
3908                 
3909                 if(surfit == 0)
3910                   aPrevNeighbourPoint.ParametersOnS1(nU2, nV2);
3911                 else
3912                   aPrevNeighbourPoint.ParametersOnS2(nU2, nV2);
3913                 gp_Vec2d aVecOld(gp_Pnt2d(nU2, nV2), gp_Pnt2d(nU1, nV1));
3914
3915                 if(aVecOld.SquareMagnitude() <= (gp::Resolution() * gp::Resolution())) {
3916                   continue;
3917                 }
3918                 else {
3919                   Standard_Real anAngle = aNewVec.Angle(aVecOld);
3920
3921                   if((fabs(anAngle) < (M_PI * 0.25)) && (aNewVec.Dot(aVecOld) > 0.)) {
3922
3923                     if(bCheckAngle1) {
3924                       Standard_Real U1, U2, V1, V2;
3925                       IntSurf_PntOn2S atmppoint = aNewP;
3926                       atmppoint.SetValue((surfit == 0), anewU, anewV);
3927                       atmppoint.Parameters(U1, V1, U2, V2);
3928                       gp_Pnt P1 = theSurface1->Value(U1, V1);