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