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