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         //  on regarde si on garde
 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 { on regarde si on garde
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         //-- lbr :
1466         //-- Si une des surfaces est un plan , on approxime en 2d
1467         //-- sur cette surface et on remonte les points 2d en 3d.
1468         if(typs1 == GeomAbs_Plane) {
1469           theapp3d.Perform(myHS1, myHS2, WL, Standard_False,Standard_True, myApprox2,ifprm,ilprm);
1470         }
1471         else if(typs2 == GeomAbs_Plane) {
1472           theapp3d.Perform(myHS1,myHS2,WL,Standard_False,myApprox1,Standard_True,ifprm,ilprm);
1473         }
1474         else {
1475           //
1476           if (myHS1 != myHS2){
1477             if ((typs1==GeomAbs_BezierSurface || typs1==GeomAbs_BSplineSurface) &&
1478                 (typs2==GeomAbs_BezierSurface || typs2==GeomAbs_BSplineSurface)) {
1479
1480               theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, Standard_True, aParType);
1481
1482               Standard_Boolean bUseSurfaces;
1483               bUseSurfaces=NotUseSurfacesForApprox(myFace1, myFace2, WL, ifprm,  ilprm);
1484               if (bUseSurfaces) {
1485                 // ######
1486                 rejectSurface = Standard_True;
1487                 // ######
1488                 theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, Standard_False, aParType);
1489               }
1490             }
1491           }
1492           //
1493           theapp3d.Perform(myHS1,myHS2,WL,Standard_True,anApprox1,anApprox2,ifprm,ilprm);
1494         }
1495
1496         if (!theapp3d.IsDone()) {
1497           //
1498           Handle(Geom2d_BSplineCurve) H1;
1499           //
1500           Handle(Geom_Curve) aBSp=MakeBSpline(WL,ifprm, ilprm);
1501           Handle(Geom2d_BSplineCurve) H2;
1502
1503           if(myApprox1) {
1504             H1 = MakeBSpline2d(WL, ifprm, ilprm, Standard_True);
1505           }
1506
1507           if(myApprox2) {
1508             H2 = MakeBSpline2d(WL, ifprm, ilprm, Standard_False);
1509           }
1510           //
1511           IntTools_Curve aIC(aBSp, H1, H2);
1512           mySeqOfCurve.Append(aIC);
1513         }
1514
1515         else {
1516           if(myApprox1 || myApprox2 || (typs1==GeomAbs_Plane || typs2==GeomAbs_Plane)) {
1517             if( theapp3d.TolReached2d()>myTolReached2d || myTolReached2d==0.) {
1518               myTolReached2d = theapp3d.TolReached2d();
1519             }
1520           }
1521           if(typs1==GeomAbs_Plane || typs2==GeomAbs_Plane) {
1522             myTolReached3d = myTolReached2d;
1523             //
1524             if (typs1==GeomAbs_Torus || typs2==GeomAbs_Torus) {
1525               if (myTolReached3d<1.e-6) {
1526                 myTolReached3d = theapp3d.TolReached3d();
1527                 myTolReached3d=1.e-6;
1528               }
1529             }
1530             //
1531           }
1532           else  if( theapp3d.TolReached3d()>myTolReached3d || myTolReached3d==0.) {
1533             myTolReached3d = theapp3d.TolReached3d();
1534           }
1535
1536           Standard_Integer aNbMultiCurves, nbpoles;
1537           aNbMultiCurves=theapp3d.NbMultiCurves();
1538           for (j=1; j<=aNbMultiCurves; j++) {
1539             if(typs1 == GeomAbs_Plane) {
1540               const AppParCurves_MultiBSpCurve& mbspc = theapp3d.Value(j);
1541               nbpoles = mbspc.NbPoles();
1542
1543               TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
1544               TColgp_Array1OfPnt   tpoles(1,nbpoles);
1545
1546               mbspc.Curve(1,tpoles2d);
1547               const gp_Pln&  Pln = myHS1->Surface().Plane();
1548               //
1549               Standard_Integer ik;
1550               for(ik = 1; ik<= nbpoles; ik++) {
1551                 tpoles.SetValue(ik,
1552                                 ElSLib::Value(tpoles2d.Value(ik).X(),
1553                                               tpoles2d.Value(ik).Y(),
1554                                               Pln));
1555               }
1556               //
1557               Handle(Geom_BSplineCurve) BS =
1558                 new Geom_BSplineCurve(tpoles,
1559                                       mbspc.Knots(),
1560                                       mbspc.Multiplicities(),
1561                                       mbspc.Degree());
1562               GeomLib_CheckBSplineCurve Check(BS,TOLCHECK,TOLANGCHECK);
1563               Check.FixTangent(Standard_True, Standard_True);
1564               //
1565               IntTools_Curve aCurve;
1566               aCurve.SetCurve(BS);
1567
1568               if(myApprox1) {
1569                 Handle(Geom2d_BSplineCurve) BS1 =
1570                   new Geom2d_BSplineCurve(tpoles2d,
1571                                           mbspc.Knots(),
1572                                           mbspc.Multiplicities(),
1573                                           mbspc.Degree());
1574                 GeomLib_Check2dBSplineCurve Check1(BS1,TOLCHECK,TOLANGCHECK);
1575                 Check1.FixTangent(Standard_True,Standard_True);
1576                 //
1577                 // ############################################
1578                 if(!rejectSurface && !reApprox) {
1579                   Standard_Boolean isValid = IsCurveValid(BS1);
1580                   if(!isValid) {
1581                     reApprox = Standard_True;
1582                     goto reapprox;
1583                   }
1584                 }
1585                 // ############################################
1586                 aCurve.SetFirstCurve2d(BS1);
1587               }
1588               else {
1589                 Handle(Geom2d_BSplineCurve) H1;
1590                 aCurve.SetFirstCurve2d(H1);
1591               }
1592
1593               if(myApprox2) {
1594                 mbspc.Curve(2, tpoles2d);
1595
1596                 Handle(Geom2d_BSplineCurve) BS2 = new Geom2d_BSplineCurve(tpoles2d,
1597                                                                           mbspc.Knots(),
1598                                                                           mbspc.Multiplicities(),
1599                                                                           mbspc.Degree());
1600                 GeomLib_Check2dBSplineCurve newCheck(BS2,TOLCHECK,TOLANGCHECK);
1601                 newCheck.FixTangent(Standard_True,Standard_True);
1602
1603                 // ###########################################
1604                 if(!rejectSurface && !reApprox) {
1605                   Standard_Boolean isValid = IsCurveValid(BS2);
1606                   if(!isValid) {
1607                     reApprox = Standard_True;
1608                     goto reapprox;
1609                   }
1610                 }
1611                 // ###########################################
1612                 //
1613                 aCurve.SetSecondCurve2d(BS2);
1614               }
1615               else {
1616                 Handle(Geom2d_BSplineCurve) H2;
1617                 //
1618                   aCurve.SetSecondCurve2d(H2);
1619               }
1620               //
1621               mySeqOfCurve.Append(aCurve);
1622             }
1623
1624             else if(typs2 == GeomAbs_Plane) {
1625               const AppParCurves_MultiBSpCurve& mbspc = theapp3d.Value(j);
1626               nbpoles = mbspc.NbPoles();
1627
1628               TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
1629               TColgp_Array1OfPnt   tpoles(1,nbpoles);
1630               mbspc.Curve((myApprox1==Standard_True)? 2 : 1,tpoles2d);
1631               const gp_Pln&  Pln = myHS2->Surface().Plane();
1632               //
1633               Standard_Integer ik;
1634               for(ik = 1; ik<= nbpoles; ik++) {
1635                 tpoles.SetValue(ik,
1636                                 ElSLib::Value(tpoles2d.Value(ik).X(),
1637                                               tpoles2d.Value(ik).Y(),
1638                                               Pln));
1639
1640               }
1641               //
1642               Handle(Geom_BSplineCurve) BS=new Geom_BSplineCurve(tpoles,
1643                                                                  mbspc.Knots(),
1644                                                                  mbspc.Multiplicities(),
1645                                                                  mbspc.Degree());
1646               GeomLib_CheckBSplineCurve Check(BS,TOLCHECK,TOLANGCHECK);
1647               Check.FixTangent(Standard_True,Standard_True);
1648               //
1649               IntTools_Curve aCurve;
1650               aCurve.SetCurve(BS);
1651
1652               if(myApprox2) {
1653                 Handle(Geom2d_BSplineCurve) BS1=new Geom2d_BSplineCurve(tpoles2d,
1654                                                                         mbspc.Knots(),
1655                                                                         mbspc.Multiplicities(),
1656                                                                         mbspc.Degree());
1657                 GeomLib_Check2dBSplineCurve Check1(BS1,TOLCHECK,TOLANGCHECK);
1658                 Check1.FixTangent(Standard_True,Standard_True);
1659                 //
1660                 // ###########################################
1661                 if(!rejectSurface && !reApprox) {
1662                   Standard_Boolean isValid = IsCurveValid(BS1);
1663                   if(!isValid) {
1664                     reApprox = Standard_True;
1665                     goto reapprox;
1666                   }
1667                 }
1668                 // ###########################################
1669                 aCurve.SetSecondCurve2d(BS1);
1670               }
1671               else {
1672                 Handle(Geom2d_BSplineCurve) H2;
1673                 aCurve.SetSecondCurve2d(H2);
1674               }
1675
1676               if(myApprox1) {
1677                 mbspc.Curve(1,tpoles2d);
1678                 Handle(Geom2d_BSplineCurve) BS2=new Geom2d_BSplineCurve(tpoles2d,
1679                                                                         mbspc.Knots(),
1680                                                                         mbspc.Multiplicities(),
1681                                                                         mbspc.Degree());
1682                 GeomLib_Check2dBSplineCurve Check2(BS2,TOLCHECK,TOLANGCHECK);
1683                 Check2.FixTangent(Standard_True,Standard_True);
1684                 //
1685                 // ###########################################
1686                 if(!rejectSurface && !reApprox) {
1687                   Standard_Boolean isValid = IsCurveValid(BS2);
1688                   if(!isValid) {
1689                     reApprox = Standard_True;
1690                     goto reapprox;
1691                   }
1692                 }
1693                 // ###########################################
1694                 aCurve.SetFirstCurve2d(BS2);
1695               }
1696               else {
1697                 Handle(Geom2d_BSplineCurve) H1;
1698                 //
1699                 aCurve.SetFirstCurve2d(H1);
1700               }
1701               //
1702               mySeqOfCurve.Append(aCurve);
1703             }
1704             else {
1705               const AppParCurves_MultiBSpCurve& mbspc = theapp3d.Value(j);
1706               nbpoles = mbspc.NbPoles();
1707               TColgp_Array1OfPnt tpoles(1,nbpoles);
1708               mbspc.Curve(1,tpoles);
1709               Handle(Geom_BSplineCurve) BS=new Geom_BSplineCurve(tpoles,
1710                                                                  mbspc.Knots(),
1711                                                                  mbspc.Multiplicities(),
1712                                                                  mbspc.Degree());
1713               GeomLib_CheckBSplineCurve Check(BS,TOLCHECK,TOLANGCHECK);
1714               Check.FixTangent(Standard_True,Standard_True);
1715               //
1716               IntTools_Curve aCurve;
1717               aCurve.SetCurve(BS);
1718
1719               if(myApprox1) {
1720                 if(anApprox1) {
1721                   TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
1722                   mbspc.Curve(2,tpoles2d);
1723                   Handle(Geom2d_BSplineCurve) BS1=new Geom2d_BSplineCurve(tpoles2d,
1724                                                                         mbspc.Knots(),
1725                                                                         mbspc.Multiplicities(),
1726                                                                         mbspc.Degree());
1727                   GeomLib_Check2dBSplineCurve newCheck(BS1,TOLCHECK,TOLANGCHECK);
1728                   newCheck.FixTangent(Standard_True,Standard_True);
1729                   //
1730                   aCurve.SetFirstCurve2d(BS1);
1731                 }
1732                 else {
1733                   Handle(Geom2d_BSplineCurve) BS1;
1734                   fprm = BS->FirstParameter();
1735                   lprm = BS->LastParameter();
1736
1737                   Handle(Geom2d_Curve) C2d;
1738                   Standard_Real aTol = myTolApprox;
1739                   BuildPCurves(fprm, lprm, aTol, myHS1->ChangeSurface().Surface(), BS, C2d);
1740                   BS1 = Handle(Geom2d_BSplineCurve)::DownCast(C2d);
1741                   aCurve.SetFirstCurve2d(BS1);
1742                 }
1743
1744               }
1745               else {
1746                 Handle(Geom2d_BSplineCurve) H1;
1747                 //
1748                 aCurve.SetFirstCurve2d(H1);
1749               }
1750               if(myApprox2) {
1751                 if(anApprox2) {
1752                   TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
1753                   mbspc.Curve((myApprox1==Standard_True)? 3 : 2,tpoles2d);
1754                   Handle(Geom2d_BSplineCurve) BS2=new Geom2d_BSplineCurve(tpoles2d,
1755                                                                         mbspc.Knots(),
1756                                                                         mbspc.Multiplicities(),
1757                                                                         mbspc.Degree());
1758                   GeomLib_Check2dBSplineCurve newCheck(BS2,TOLCHECK,TOLANGCHECK);
1759                   newCheck.FixTangent(Standard_True,Standard_True);
1760                 //
1761                   aCurve.SetSecondCurve2d(BS2);
1762                 }
1763                 else {
1764                   Handle(Geom2d_BSplineCurve) BS2;
1765                   fprm = BS->FirstParameter();
1766                   lprm = BS->LastParameter();
1767
1768                   Handle(Geom2d_Curve) C2d;
1769                   Standard_Real aTol = myTolApprox;
1770                   BuildPCurves(fprm, lprm, aTol, myHS2->ChangeSurface().Surface(), BS, C2d);
1771                   BS2 = Handle(Geom2d_BSplineCurve)::DownCast(C2d);
1772                   aCurve.SetSecondCurve2d(BS2);
1773                 }
1774
1775               }
1776               else {
1777                 Handle(Geom2d_BSplineCurve) H2;
1778                 //
1779                 aCurve.SetSecondCurve2d(H2);
1780               }
1781               //
1782               mySeqOfCurve.Append(aCurve);
1783             }
1784           }
1785         }
1786       }
1787     }// else { // X
1788   }// case IntPatch_Walking:{
1789     break;
1790
1791   case IntPatch_Restriction:
1792     break;
1793
1794   }
1795 }
1796
1797 //=======================================================================
1798 //function : BuildPCurves
1799 //purpose  :
1800 //=======================================================================
1801  void BuildPCurves (Standard_Real f,
1802                     Standard_Real l,
1803                     Standard_Real& Tol,
1804                     const Handle (Geom_Surface)& S,
1805                     const Handle (Geom_Curve)&   C,
1806                     Handle (Geom2d_Curve)& C2d)
1807 {
1808
1809   Standard_Real umin,umax,vmin,vmax;
1810   //
1811
1812   if (C2d.IsNull()) {
1813
1814     // in class ProjLib_Function the range of parameters is shrank by 1.e-09
1815     if((l - f) > 2.e-09) {
1816       C2d = GeomProjLib::Curve2d(C,f,l,S,Tol);
1817       //
1818       if (C2d.IsNull()) {
1819         // proj. a circle that goes through the pole on a sphere to the sphere
1820         Tol=Tol+1.e-7;
1821         C2d = GeomProjLib::Curve2d(C,f,l,S,Tol);
1822       }
1823     }
1824     else {
1825       if((l - f) > Epsilon(Abs(f))) {
1826         GeomAPI_ProjectPointOnSurf aProjector1, aProjector2;
1827         gp_Pnt P1 = C->Value(f);
1828         gp_Pnt P2 = C->Value(l);
1829         aProjector1.Init(P1, S);
1830         aProjector2.Init(P2, S);
1831
1832         if(aProjector1.IsDone() && aProjector2.IsDone()) {
1833           Standard_Real U=0., V=0.;
1834           aProjector1.LowerDistanceParameters(U, V);
1835           gp_Pnt2d p1(U, V);
1836
1837           aProjector2.LowerDistanceParameters(U, V);
1838           gp_Pnt2d p2(U, V);
1839
1840           if(p1.Distance(p2) > gp::Resolution()) {
1841             TColgp_Array1OfPnt2d poles(1,2);
1842             TColStd_Array1OfReal knots(1,2);
1843             TColStd_Array1OfInteger mults(1,2);
1844             poles(1) = p1;
1845             poles(2) = p2;
1846             knots(1) = f;
1847             knots(2) = l;
1848             mults(1) = mults(2) = 2;
1849
1850             C2d = new Geom2d_BSplineCurve(poles,knots,mults,1);
1851
1852             // compute reached tolerance.begin
1853             gp_Pnt PMid = C->Value((f + l) * 0.5);
1854             aProjector1.Perform(PMid);
1855
1856             if(aProjector1.IsDone()) {
1857               aProjector1.LowerDistanceParameters(U, V);
1858               gp_Pnt2d pmidproj(U, V);
1859               gp_Pnt2d pmidcurve2d = C2d->Value((f + l) * 0.5);
1860               Standard_Real adist = pmidcurve2d.Distance(pmidproj);
1861               Tol = (adist > Tol) ? adist : Tol;
1862             }
1863             // compute reached tolerance.end
1864           }
1865         }
1866       }
1867     }
1868     //
1869     S->Bounds(umin, umax, vmin, vmax);
1870
1871     if (S->IsUPeriodic() && !C2d.IsNull()) {
1872       // Recadre dans le domaine UV de la face
1873       Standard_Real period, U0, du, aEps;
1874
1875       du =0.0;
1876       aEps=Precision::PConfusion();
1877       period = S->UPeriod();
1878       gp_Pnt2d Pf = C2d->Value(f);
1879       U0=Pf.X();
1880       //
1881       gp_Pnt2d Pl = C2d->Value(l);
1882
1883       U0 = Min(Pl.X(), U0);
1884 //       while(U0-umin<aEps) {
1885       while(U0-umin<-aEps) {
1886         U0+=period;
1887         du+=period;
1888       }
1889       //
1890       while(U0-umax>aEps) {
1891         U0-=period;
1892         du-=period;
1893       }
1894       if (du != 0) {
1895         gp_Vec2d T1(du,0.);
1896         C2d->Translate(T1);
1897       }
1898     }
1899   }
1900   if (C2d.IsNull()) {
1901     BOPTColStd_Dump::PrintMessage("BuildPCurves()=> Echec ProjLib\n");
1902   }
1903 }
1904
1905 //=======================================================================
1906 //function : Parameters
1907 //purpose  :
1908 //=======================================================================
1909  void Parameters(const Handle(GeomAdaptor_HSurface)& HS1,
1910                  const Handle(GeomAdaptor_HSurface)& HS2,
1911                  const gp_Pnt& Ptref,
1912                  Standard_Real& U1,
1913                  Standard_Real& V1,
1914                  Standard_Real& U2,
1915                  Standard_Real& V2)
1916 {
1917
1918   IntSurf_Quadric quad1,quad2;
1919   GeomAbs_SurfaceType typs = HS1->Surface().GetType();
1920
1921   switch (typs) {
1922   case GeomAbs_Plane:
1923     quad1.SetValue(HS1->Surface().Plane());
1924     break;
1925   case GeomAbs_Cylinder:
1926     quad1.SetValue(HS1->Surface().Cylinder());
1927     break;
1928   case GeomAbs_Cone:
1929     quad1.SetValue(HS1->Surface().Cone());
1930     break;
1931   case GeomAbs_Sphere:
1932     quad1.SetValue(HS1->Surface().Sphere());
1933     break;
1934   default:
1935     Standard_ConstructionError::Raise("GeomInt_IntSS::MakeCurve");
1936   }
1937
1938   typs = HS2->Surface().GetType();
1939   switch (typs) {
1940   case GeomAbs_Plane:
1941     quad2.SetValue(HS2->Surface().Plane());
1942     break;
1943   case GeomAbs_Cylinder:
1944     quad2.SetValue(HS2->Surface().Cylinder());
1945     break;
1946   case GeomAbs_Cone:
1947     quad2.SetValue(HS2->Surface().Cone());
1948     break;
1949   case GeomAbs_Sphere:
1950     quad2.SetValue(HS2->Surface().Sphere());
1951     break;
1952   default:
1953     Standard_ConstructionError::Raise("GeomInt_IntSS::MakeCurve");
1954   }
1955
1956   quad1.Parameters(Ptref,U1,V1);
1957   quad2.Parameters(Ptref,U2,V2);
1958 }
1959
1960 //=======================================================================
1961 //function : MakeBSpline
1962 //purpose  :
1963 //=======================================================================
1964 Handle(Geom_Curve) MakeBSpline  (const Handle(IntPatch_WLine)& WL,
1965                                  const Standard_Integer ideb,
1966                                  const Standard_Integer ifin)
1967 {
1968   Standard_Integer i,nbpnt = ifin-ideb+1;
1969   TColgp_Array1OfPnt poles(1,nbpnt);
1970   TColStd_Array1OfReal knots(1,nbpnt);
1971   TColStd_Array1OfInteger mults(1,nbpnt);
1972   Standard_Integer ipidebm1;
1973   for(i=1,ipidebm1=i+ideb-1; i<=nbpnt;ipidebm1++, i++) {
1974     poles(i) = WL->Point(ipidebm1).Value();
1975     mults(i) = 1;
1976     knots(i) = i-1;
1977   }
1978   mults(1) = mults(nbpnt) = 2;
1979   return
1980     new Geom_BSplineCurve(poles,knots,mults,1);
1981 }
1982 //
1983
1984 //=======================================================================
1985 //function : MakeBSpline2d
1986 //purpose  :
1987 //=======================================================================
1988 Handle(Geom2d_BSplineCurve) MakeBSpline2d(const Handle(IntPatch_WLine)& theWLine,
1989                                           const Standard_Integer ideb,
1990                                           const Standard_Integer ifin,
1991                                           const Standard_Boolean onFirst)
1992 {
1993   Standard_Integer i, nbpnt = ifin-ideb+1;
1994   TColgp_Array1OfPnt2d poles(1,nbpnt);
1995   TColStd_Array1OfReal knots(1,nbpnt);
1996   TColStd_Array1OfInteger mults(1,nbpnt);
1997   Standard_Integer ipidebm1;
1998
1999   for(i = 1, ipidebm1 = i+ideb-1; i <= nbpnt; ipidebm1++, i++) {
2000       Standard_Real U, V;
2001       if(onFirst)
2002         theWLine->Point(ipidebm1).ParametersOnS1(U, V);
2003       else
2004         theWLine->Point(ipidebm1).ParametersOnS2(U, V);
2005       poles(i).SetCoord(U, V);
2006       mults(i) = 1;
2007       knots(i) = i-1;
2008     }
2009     mults(1) = mults(nbpnt) = 2;
2010
2011   return new Geom2d_BSplineCurve(poles,knots,mults,1);
2012 }
2013 //=======================================================================
2014 //function : PrepareLines3D
2015 //purpose  :
2016 //=======================================================================
2017   void IntTools_FaceFace::PrepareLines3D()
2018 {
2019   Standard_Integer i, aNbCurves, j, aNbNewCurves;
2020   IntTools_SequenceOfCurves aNewCvs;
2021
2022   //
2023   // 1. Treatment of periodic and closed  curves
2024   aNbCurves=mySeqOfCurve.Length();
2025   for (i=1; i<=aNbCurves; i++) {
2026     const IntTools_Curve& aIC=mySeqOfCurve(i);
2027     // DEBUG
2028     // const Handle(Geom_Curve)&   aC3D =aIC.Curve();
2029     // const Handle(Geom2d_Curve)& aC2D1=aIC.FirstCurve2d();
2030     // const Handle(Geom2d_Curve)& aC2D2=aIC.SecondCurve2d();
2031     //
2032     IntTools_SequenceOfCurves aSeqCvs;
2033     aNbNewCurves=IntTools_Tools::SplitCurve(aIC, aSeqCvs);
2034
2035     if (aNbNewCurves) {
2036       for (j=1; j<=aNbNewCurves; j++) {
2037         const IntTools_Curve& aICNew=aSeqCvs(j);
2038         aNewCvs.Append(aICNew);
2039       }
2040     }
2041     //
2042     else {
2043       aNewCvs.Append(aIC);
2044     }
2045   }
2046   //
2047   // 2. Plane\Cone intersection when we had 4 curves
2048   GeomAbs_SurfaceType aType1, aType2;
2049   BRepAdaptor_Surface aBS1, aBS2;
2050
2051   aBS1.Initialize(myFace1);
2052   aType1=aBS1.GetType();
2053
2054   aBS2.Initialize(myFace2);
2055   aType2=aBS2.GetType();
2056
2057   if ((aType1==GeomAbs_Plane && aType2==GeomAbs_Cone) ||
2058       (aType2==GeomAbs_Plane && aType1==GeomAbs_Cone)) {
2059     aNbCurves=aNewCvs.Length();
2060     if (aNbCurves==4) {
2061       GeomAbs_CurveType aCType1=aNewCvs(1).Type();
2062       if (aCType1==GeomAbs_Line) {
2063         IntTools_SequenceOfCurves aSeqIn, aSeqOut;
2064         //
2065         for (i=1; i<=aNbCurves; i++) {
2066           const IntTools_Curve& aIC=aNewCvs(i);
2067           aSeqIn.Append(aIC);
2068         }
2069         //
2070         IntTools_Tools::RejectLines(aSeqIn, aSeqOut);
2071         //
2072         aNewCvs.Clear();
2073         aNbCurves=aSeqOut.Length();
2074         for (i=1; i<=aNbCurves; i++) {
2075           const IntTools_Curve& aIC=aSeqOut(i);
2076           aNewCvs.Append(aIC);
2077         }
2078         //
2079       }
2080     }
2081   }// end of if ((aType1==GeomAbs_Plane && ...
2082   //
2083   // 3. Fill  mySeqOfCurve
2084   mySeqOfCurve.Clear();
2085   aNbCurves=aNewCvs.Length();
2086   for (i=1; i<=aNbCurves; i++) {
2087     const IntTools_Curve& aIC=aNewCvs(i);
2088     mySeqOfCurve.Append(aIC);
2089   }
2090
2091 }
2092
2093
2094 //=======================================================================
2095 //function : CorrectSurfaceBoundaries
2096 //purpose  :
2097 //=======================================================================
2098  void CorrectSurfaceBoundaries(const TopoDS_Face&  theFace,
2099                                const Standard_Real theTolerance,
2100                                Standard_Real&      theumin,
2101                                Standard_Real&      theumax,
2102                                Standard_Real&      thevmin,
2103                                Standard_Real&      thevmax)
2104 {
2105   Standard_Boolean enlarge, isuperiodic, isvperiodic;
2106   Standard_Real uinf, usup, vinf, vsup, delta;
2107   GeomAbs_SurfaceType aType;
2108   Handle(Geom_Surface) aSurface;
2109   //
2110   aSurface = BRep_Tool::Surface(theFace);
2111   aSurface->Bounds(uinf, usup, vinf, vsup);
2112   delta = theTolerance;
2113   enlarge = Standard_False;
2114   //
2115   GeomAdaptor_Surface anAdaptorSurface(aSurface);
2116   //
2117   if(aSurface->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface))) {
2118     Handle(Geom_Surface) aBasisSurface =
2119       (Handle(Geom_RectangularTrimmedSurface)::DownCast(aSurface))->BasisSurface();
2120
2121     if(aBasisSurface->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface)) ||
2122        aBasisSurface->IsKind(STANDARD_TYPE(Geom_OffsetSurface))) {
2123       return;
2124     }
2125   }
2126   //
2127   if(aSurface->IsKind(STANDARD_TYPE(Geom_OffsetSurface))) {
2128     Handle(Geom_Surface) aBasisSurface =
2129       (Handle(Geom_OffsetSurface)::DownCast(aSurface))->BasisSurface();
2130
2131     if(aBasisSurface->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface)) ||
2132        aBasisSurface->IsKind(STANDARD_TYPE(Geom_OffsetSurface))) {
2133       return;
2134     }
2135   }
2136   //
2137   isuperiodic = anAdaptorSurface.IsUPeriodic();
2138   isvperiodic = anAdaptorSurface.IsVPeriodic();
2139   //
2140   aType=anAdaptorSurface.GetType();
2141   if((aType==GeomAbs_BezierSurface) ||
2142      (aType==GeomAbs_BSplineSurface) ||
2143      (aType==GeomAbs_SurfaceOfExtrusion) ||
2144      (aType==GeomAbs_SurfaceOfRevolution)) {
2145     enlarge=Standard_True;
2146   }
2147   //
2148   if(!isuperiodic && enlarge) {
2149
2150     if((theumin - uinf) > delta )
2151       theumin -= delta;
2152     else {
2153       theumin = uinf;
2154     }
2155
2156     if((usup - theumax) > delta )
2157       theumax += delta;
2158     else
2159       theumax = usup;
2160   }
2161   //
2162   if(!isvperiodic && enlarge) {
2163     if((thevmin - vinf) > delta ) {
2164       thevmin -= delta;
2165     }
2166     else {
2167       thevmin = vinf;
2168     }
2169     if((vsup - thevmax) > delta ) {
2170       thevmax += delta;
2171     }
2172     else {
2173       thevmax = vsup;
2174     }
2175   }
2176   //
2177   {
2178     Standard_Integer aNbP;
2179     Standard_Real aXP, dXfact, aXmid, aX1, aX2, aTolPA;
2180     //
2181     aTolPA=Precision::Angular();
2182     // U
2183     if (isuperiodic) {
2184       aXP=anAdaptorSurface.UPeriod();
2185       dXfact=theumax-theumin;
2186       if (dXfact-aTolPA>aXP) {
2187         aXmid=0.5*(theumax+theumin);
2188         aNbP=RealToInt(aXmid/aXP);
2189         if (aXmid<0.) {
2190           aNbP=aNbP-1;
2191         }
2192         aX1=aNbP*aXP;
2193         aX2=aX1+aXP;
2194         if (theumin<aX1) {
2195           theumin=aX1;
2196         }
2197         if (theumax>aX2) {
2198           theumax=aX2;
2199         }
2200       }
2201     }
2202     // V
2203     if (isvperiodic) {
2204       aXP=anAdaptorSurface.VPeriod();
2205       dXfact=thevmax-thevmin;
2206       if (dXfact-aTolPA>aXP) {
2207         aXmid=0.5*(thevmax+thevmin);
2208         aNbP=RealToInt(aXmid/aXP);
2209         if (aXmid<0.) {
2210           aNbP=aNbP-1;
2211         }
2212         aX1=aNbP*aXP;
2213         aX2=aX1+aXP;
2214         if (thevmin<aX1) {
2215           thevmin=aX1;
2216         }
2217         if (thevmax>aX2) {
2218           thevmax=aX2;
2219         }
2220       }
2221     }
2222   }
2223   //
2224   if(isuperiodic || isvperiodic) {
2225     Standard_Boolean correct = Standard_False;
2226     Standard_Boolean correctU = Standard_False;
2227     Standard_Boolean correctV = Standard_False;
2228     Bnd_Box2d aBox;
2229     TopExp_Explorer anExp;
2230
2231     for(anExp.Init(theFace, TopAbs_EDGE); anExp.More(); anExp.Next()) {
2232       if(BRep_Tool::IsClosed(TopoDS::Edge(anExp.Current()), theFace)) {
2233         correct = Standard_True;
2234         Standard_Real f, l;
2235         TopoDS_Edge anEdge = TopoDS::Edge(anExp.Current());
2236
2237         for(Standard_Integer i = 0; i < 2; i++) {
2238           if(i==0) {
2239             anEdge.Orientation(TopAbs_FORWARD);
2240           }
2241           else {
2242             anEdge.Orientation(TopAbs_REVERSED);
2243           }
2244           Handle(Geom2d_Curve) aCurve = BRep_Tool::CurveOnSurface(anEdge, theFace, f, l);
2245
2246           if(aCurve.IsNull()) {
2247             correct = Standard_False;
2248             break;
2249           }
2250           Handle(Geom2d_Line) aLine = Handle(Geom2d_Line)::DownCast(aCurve);
2251
2252           if(aLine.IsNull()) {
2253             correct = Standard_False;
2254             break;
2255           }
2256           gp_Dir2d anUDir(1., 0.);
2257           gp_Dir2d aVDir(0., 1.);
2258           Standard_Real anAngularTolerance = Precision::Angular();
2259
2260           correctU = correctU || aLine->Position().Direction().IsParallel(aVDir, anAngularTolerance);
2261           correctV = correctV || aLine->Position().Direction().IsParallel(anUDir, anAngularTolerance);
2262
2263           gp_Pnt2d pp1 = aCurve->Value(f);
2264           aBox.Add(pp1);
2265           gp_Pnt2d pp2 = aCurve->Value(l);
2266           aBox.Add(pp2);
2267         }
2268         if(!correct)
2269           break;
2270       }
2271     }
2272
2273     if(correct) {
2274       Standard_Real umin, vmin, umax, vmax;
2275       aBox.Get(umin, vmin, umax, vmax);
2276
2277       if(isuperiodic && correctU) {
2278
2279         if(theumin < umin)
2280           theumin = umin;
2281
2282         if(theumax > umax) {
2283           theumax = umax;
2284         }
2285       }
2286       if(isvperiodic && correctV) {
2287
2288         if(thevmin < vmin)
2289           thevmin = vmin;
2290         if(thevmax > vmax)
2291           thevmax = vmax;
2292       }
2293     }
2294   }
2295 }
2296 //
2297 //
2298 // The block is dedicated to determine whether WLine [ifprm, ilprm]
2299 // crosses the degenerated zone on each given surface or not.
2300 // If Yes -> We will not use info about surfaces during approximation
2301 // because inside degenerated zone of the surface the approx. alogo.
2302 // uses wrong values of normal, etc., and resulting curve will have
2303 // oscillations that we would not like to have.
2304 //                                               PKV Tue Feb 12 2002
2305
2306
2307 static
2308   Standard_Boolean IsDegeneratedZone(const gp_Pnt2d& aP2d,
2309                                      const Handle(Geom_Surface)& aS,
2310                                      const Standard_Integer iDir);
2311 static
2312   Standard_Boolean IsPointInDegeneratedZone(const IntSurf_PntOn2S& aP2S,
2313                                             const TopoDS_Face& aF1,
2314                                             const TopoDS_Face& aF2);
2315 //=======================================================================
2316 //function :  NotUseSurfacesForApprox
2317 //purpose  :
2318 //=======================================================================
2319 Standard_Boolean NotUseSurfacesForApprox(const TopoDS_Face& aF1,
2320                                          const TopoDS_Face& aF2,
2321                                          const Handle(IntPatch_WLine)& WL,
2322                                          const Standard_Integer ifprm,
2323                                          const Standard_Integer ilprm)
2324 {
2325   Standard_Boolean bPInDZ;
2326
2327   Handle(IntSurf_LineOn2S) aLineOn2S=WL->Curve();
2328
2329   const IntSurf_PntOn2S& aP2Sfprm=aLineOn2S->Value(ifprm);
2330   bPInDZ=IsPointInDegeneratedZone(aP2Sfprm, aF1, aF2);
2331   if (bPInDZ) {
2332     return bPInDZ;
2333   }
2334
2335   const IntSurf_PntOn2S& aP2Slprm=aLineOn2S->Value(ilprm);
2336   bPInDZ=IsPointInDegeneratedZone(aP2Slprm, aF1, aF2);
2337
2338   return bPInDZ;
2339 }
2340 //=======================================================================
2341 //function : IsPointInDegeneratedZone
2342 //purpose  :
2343 //=======================================================================
2344 Standard_Boolean IsPointInDegeneratedZone(const IntSurf_PntOn2S& aP2S,
2345                                           const TopoDS_Face& aF1,
2346                                           const TopoDS_Face& aF2)
2347
2348 {
2349   Standard_Boolean bFlag=Standard_True;
2350   Standard_Real US11, US12, VS11, VS12, US21, US22, VS21, VS22;
2351   Standard_Real U1, V1, U2, V2, aDelta, aD;
2352   gp_Pnt2d aP2d;
2353
2354   Handle(Geom_Surface)aS1 = BRep_Tool::Surface(aF1);
2355   aS1->Bounds(US11, US12, VS11, VS12);
2356   GeomAdaptor_Surface aGAS1(aS1);
2357
2358   Handle(Geom_Surface)aS2 = BRep_Tool::Surface(aF2);
2359   aS1->Bounds(US21, US22, VS21, VS22);
2360   GeomAdaptor_Surface aGAS2(aS2);
2361   //
2362   //const gp_Pnt& aP=aP2S.Value();
2363   aP2S.Parameters(U1, V1, U2, V2);
2364   //
2365   aDelta=1.e-7;
2366   // Check on Surf 1
2367   aD=aGAS1.UResolution(aDelta);
2368   aP2d.SetCoord(U1, V1);
2369   if (fabs(U1-US11) < aD) {
2370     bFlag=IsDegeneratedZone(aP2d, aS1, 1);
2371     if (bFlag) {
2372       return bFlag;
2373     }
2374   }
2375   if (fabs(U1-US12) < aD) {
2376     bFlag=IsDegeneratedZone(aP2d, aS1, 1);
2377     if (bFlag) {
2378       return bFlag;
2379     }
2380   }
2381   aD=aGAS1.VResolution(aDelta);
2382   if (fabs(V1-VS11) < aDelta) {
2383     bFlag=IsDegeneratedZone(aP2d, aS1, 2);
2384     if (bFlag) {
2385       return bFlag;
2386     }
2387   }
2388   if (fabs(V1-VS12) < aDelta) {
2389     bFlag=IsDegeneratedZone(aP2d, aS1, 2);
2390     if (bFlag) {
2391       return bFlag;
2392     }
2393   }
2394   // Check on Surf 2
2395   aD=aGAS2.UResolution(aDelta);
2396   aP2d.SetCoord(U2, V2);
2397   if (fabs(U2-US21) < aDelta) {
2398     bFlag=IsDegeneratedZone(aP2d, aS2, 1);
2399     if (bFlag) {
2400       return bFlag;
2401     }
2402   }
2403   if (fabs(U2-US22) < aDelta) {
2404     bFlag=IsDegeneratedZone(aP2d, aS2, 1);
2405     if (bFlag) {
2406       return bFlag;
2407     }
2408   }
2409   aD=aGAS2.VResolution(aDelta);
2410   if (fabs(V2-VS21) < aDelta) {
2411     bFlag=IsDegeneratedZone(aP2d, aS2, 2);
2412     if (bFlag) {
2413       return bFlag;
2414     }
2415   }
2416   if (fabs(V2-VS22) < aDelta) {
2417     bFlag=IsDegeneratedZone(aP2d, aS2, 2);
2418     if (bFlag) {
2419       return bFlag;
2420     }
2421   }
2422   return !bFlag;
2423 }
2424
2425 //=======================================================================
2426 //function : IsDegeneratedZone
2427 //purpose  :
2428 //=======================================================================
2429 Standard_Boolean IsDegeneratedZone(const gp_Pnt2d& aP2d,
2430                                    const Handle(Geom_Surface)& aS,
2431                                    const Standard_Integer iDir)
2432 {
2433   Standard_Boolean bFlag=Standard_True;
2434   Standard_Real US1, US2, VS1, VS2, dY, dX, d1, d2, dD;
2435   Standard_Real aXm, aYm, aXb, aYb, aXe, aYe;
2436   aS->Bounds(US1, US2, VS1, VS2);
2437
2438   gp_Pnt aPm, aPb, aPe;
2439
2440   aXm=aP2d.X();
2441   aYm=aP2d.Y();
2442
2443   aS->D0(aXm, aYm, aPm);
2444
2445   dX=1.e-5;
2446   dY=1.e-5;
2447   dD=1.e-12;
2448
2449   if (iDir==1) {
2450     aXb=aXm;
2451     aXe=aXm;
2452     aYb=aYm-dY;
2453     if (aYb < VS1) {
2454       aYb=VS1;
2455     }
2456     aYe=aYm+dY;
2457     if (aYe > VS2) {
2458       aYe=VS2;
2459     }
2460     aS->D0(aXb, aYb, aPb);
2461     aS->D0(aXe, aYe, aPe);
2462
2463     d1=aPm.Distance(aPb);
2464     d2=aPm.Distance(aPe);
2465     if (d1 < dD && d2 < dD) {
2466       return bFlag;
2467     }
2468     return !bFlag;
2469   }
2470   //
2471   else if (iDir==2) {
2472     aYb=aYm;
2473     aYe=aYm;
2474     aXb=aXm-dX;
2475     if (aXb < US1) {
2476       aXb=US1;
2477     }
2478     aXe=aXm+dX;
2479     if (aXe > US2) {
2480       aXe=US2;
2481     }
2482     aS->D0(aXb, aYb, aPb);
2483     aS->D0(aXe, aYe, aPe);
2484
2485     d1=aPm.Distance(aPb);
2486     d2=aPm.Distance(aPe);
2487     if (d1 < dD && d2 < dD) {
2488       return bFlag;
2489     }
2490     return !bFlag;
2491   }
2492   return !bFlag;
2493 }
2494
2495 //=========================================================================
2496 // static function : ComputePurgedWLine
2497 // purpose : Removes equal points (leave one of equal points) from theWLine
2498 //           and recompute vertex parameters.
2499 //           Returns new WLine or null WLine if the number
2500 //           of the points is less than 2.
2501 //=========================================================================
2502 Handle(IntPatch_WLine) ComputePurgedWLine(const Handle(IntPatch_WLine)& theWLine) {
2503   Handle(IntPatch_WLine) aResult;
2504   Handle(IntPatch_WLine) aLocalWLine;
2505   Handle(IntPatch_WLine) aTmpWLine = theWLine;
2506
2507   Handle(IntSurf_LineOn2S) aLineOn2S = new IntSurf_LineOn2S();
2508   aLocalWLine = new IntPatch_WLine(aLineOn2S, Standard_False);
2509   Standard_Integer i, k, v, nb, nbvtx;
2510   nbvtx = theWLine->NbVertex();
2511   nb = theWLine->NbPnts();
2512
2513   for(i = 1; i <= nb; i++) {
2514     aLineOn2S->Add(theWLine->Point(i));
2515   }
2516
2517   for(v = 1; v <= nbvtx; v++) {
2518     aLocalWLine->AddVertex(theWLine->Vertex(v));
2519   }
2520
2521   for(i = 1; i <= aLineOn2S->NbPoints(); i++) {
2522     Standard_Integer aStartIndex = i + 1;
2523     Standard_Integer anEndIndex = i + 5;
2524     nb = aLineOn2S->NbPoints();
2525     anEndIndex = (anEndIndex > nb) ? nb : anEndIndex;
2526
2527     if((aStartIndex >= nb) || (anEndIndex <= 1)) {
2528       continue;
2529     }
2530     k = aStartIndex;
2531
2532     while(k <= anEndIndex) {
2533
2534       if(i != k) {
2535         IntSurf_PntOn2S p1 = aLineOn2S->Value(i);
2536         IntSurf_PntOn2S p2 = aLineOn2S->Value(k);
2537
2538         if(p1.Value().IsEqual(p2.Value(), gp::Resolution())) {
2539           aTmpWLine = aLocalWLine;
2540           aLocalWLine = new IntPatch_WLine(aLineOn2S, Standard_False);
2541
2542           for(v = 1; v <= aTmpWLine->NbVertex(); v++) {
2543             IntPatch_Point aVertex = aTmpWLine->Vertex(v);
2544             Standard_Integer avertexindex = (Standard_Integer)aVertex.ParameterOnLine();
2545
2546             if(avertexindex >= k) {
2547               aVertex.SetParameter(aVertex.ParameterOnLine() - 1.);
2548             }
2549             aLocalWLine->AddVertex(aVertex);
2550           }
2551           aLineOn2S->RemovePoint(k);
2552           anEndIndex--;
2553           continue;
2554         }
2555       }
2556       k++;
2557     }
2558   }
2559
2560   if(aLineOn2S->NbPoints() > 1) {
2561     aResult = aLocalWLine;
2562   }
2563   return aResult;
2564 }
2565
2566 //=======================================================================
2567 //function : TolR3d
2568 //purpose  :
2569 //=======================================================================
2570 void TolR3d(const TopoDS_Face& aF1,
2571             const TopoDS_Face& aF2,
2572             Standard_Real& myTolReached3d)
2573 {
2574   Standard_Real aTolF1, aTolF2, aTolFMax, aTolTresh;
2575
2576   aTolTresh=2.999999e-3;
2577   aTolF1 = BRep_Tool::Tolerance(aF1);
2578   aTolF2 = BRep_Tool::Tolerance(aF2);
2579   aTolFMax=Max(aTolF1, aTolF2);
2580
2581   if (aTolFMax>aTolTresh) {
2582     myTolReached3d=aTolFMax;
2583   }
2584 }
2585 //=======================================================================
2586 //function : AdjustPeriodic
2587 //purpose  :
2588 //=======================================================================
2589 Standard_Real AdjustPeriodic(const Standard_Real theParameter,
2590                              const Standard_Real parmin,
2591                              const Standard_Real parmax,
2592                              const Standard_Real thePeriod,
2593                              Standard_Real&      theOffset)
2594 {
2595   Standard_Real aresult;
2596   //
2597   theOffset = 0.;
2598   aresult = theParameter;
2599   while(aresult < parmin) {
2600     aresult += thePeriod;
2601     theOffset += thePeriod;
2602   }
2603
2604   while(aresult > parmax) {
2605     aresult -= thePeriod;
2606     theOffset -= thePeriod;
2607   }
2608   return aresult;
2609 }
2610 //=======================================================================
2611 //function : IsPointOnBoundary
2612 //purpose  :
2613 //=======================================================================
2614 Standard_Boolean IsPointOnBoundary(const Standard_Real theParameter,
2615                                    const Standard_Real theFirstBoundary,
2616                                    const Standard_Real theSecondBoundary,
2617                                    const Standard_Real theResolution,
2618                                    Standard_Boolean&   IsOnFirstBoundary)
2619 {
2620   Standard_Boolean bRet;
2621   Standard_Integer i;
2622   Standard_Real adist;
2623   //
2624   bRet=Standard_False;
2625   for(i = 0; i < 2; ++i) {
2626     IsOnFirstBoundary = (i == 0);
2627     if (IsOnFirstBoundary) {
2628       adist = fabs(theParameter - theFirstBoundary);
2629     }
2630     else {
2631       adist = fabs(theParameter - theSecondBoundary);
2632     }
2633     if(adist < theResolution) {
2634       return !bRet;
2635     }
2636   }
2637   return bRet;
2638 }
2639 // ------------------------------------------------------------------------------------------------
2640 // static function: FindPoint
2641 // purpose:
2642 // ------------------------------------------------------------------------------------------------
2643 Standard_Boolean FindPoint(const gp_Pnt2d&     theFirstPoint,
2644                            const gp_Pnt2d&     theLastPoint,
2645                            const Standard_Real theUmin,
2646                            const Standard_Real theUmax,
2647                            const Standard_Real theVmin,
2648                            const Standard_Real theVmax,
2649                            gp_Pnt2d&           theNewPoint) {
2650
2651   gp_Vec2d aVec(theFirstPoint, theLastPoint);
2652   Standard_Integer i = 0, j = 0;
2653
2654   for(i = 0; i < 4; i++) {
2655     gp_Vec2d anOtherVec;
2656     gp_Vec2d anOtherVecNormal;
2657     gp_Pnt2d aprojpoint = theLastPoint;
2658
2659     if((i % 2) == 0) {
2660       anOtherVec.SetX(0.);
2661       anOtherVec.SetY(1.);
2662       anOtherVecNormal.SetX(1.);
2663       anOtherVecNormal.SetY(0.);
2664
2665       if(i < 2)
2666         aprojpoint.SetX(theUmin);
2667       else
2668         aprojpoint.SetX(theUmax);
2669     }
2670     else {
2671       anOtherVec.SetX(1.);
2672       anOtherVec.SetY(0.);
2673       anOtherVecNormal.SetX(0.);
2674       anOtherVecNormal.SetY(1.);
2675
2676       if(i < 2)
2677         aprojpoint.SetY(theVmin);
2678       else
2679         aprojpoint.SetY(theVmax);
2680     }
2681     gp_Vec2d anormvec = aVec;
2682     anormvec.Normalize();
2683     Standard_Real adot1 = anormvec.Dot(anOtherVecNormal);
2684
2685     if(fabs(adot1) < Precision::Angular())
2686       continue;
2687     Standard_Real adist = 0.;
2688     Standard_Boolean bIsOut = Standard_False;
2689
2690     if((i % 2) == 0) {
2691       adist = (i < 2) ? fabs(theLastPoint.X() - theUmin) : fabs(theLastPoint.X() - theUmax);
2692       bIsOut = (i < 2) ? (theLastPoint.X() < theUmin) : (theLastPoint.X() > theUmax);
2693     }
2694     else {
2695       adist = (i < 2) ? fabs(theLastPoint.Y() - theVmin) : fabs(theLastPoint.Y() - theVmax);
2696       bIsOut = (i < 2) ? (theLastPoint.Y() < theVmin) : (theLastPoint.Y() > theVmax);
2697     }
2698     Standard_Real anoffset = adist * anOtherVec.Dot(anormvec) / adot1;
2699
2700     for(j = 0; j < 2; j++) {
2701       anoffset = (j == 0) ? anoffset : -anoffset;
2702       gp_Pnt2d acurpoint(aprojpoint.XY() + (anOtherVec.XY()*anoffset));
2703       gp_Vec2d acurvec(theLastPoint, acurpoint);
2704       if ( bIsOut )
2705         acurvec.Reverse();
2706
2707       if((aVec.Dot(acurvec) > 0.) &&
2708          (aVec.Angle(acurvec) < Precision::PConfusion())) {
2709         if((i % 2) == 0) {
2710           if((acurpoint.Y() >= theVmin) &&
2711              (acurpoint.Y() <= theVmax)) {
2712             theNewPoint = acurpoint;
2713             return Standard_True;
2714           }
2715         }
2716         else {
2717           if((acurpoint.X() >= theUmin) &&
2718              (acurpoint.X() <= theUmax)) {
2719             theNewPoint = acurpoint;
2720             return Standard_True;
2721           }
2722         }
2723       }
2724     }
2725   }
2726   return Standard_False;
2727 }
2728
2729
2730 // ------------------------------------------------------------------------------------------------
2731 // static function: FindPoint
2732 // purpose: Find point on the boundary of radial tangent zone
2733 // ------------------------------------------------------------------------------------------------
2734 Standard_Boolean FindPoint(const gp_Pnt2d&     theFirstPoint,
2735                            const gp_Pnt2d&     theLastPoint,
2736                            const Standard_Real theUmin,
2737                            const Standard_Real theUmax,
2738                            const Standard_Real theVmin,
2739                            const Standard_Real theVmax,
2740                            const gp_Pnt2d&     theTanZoneCenter,
2741                            const Standard_Real theZoneRadius,
2742                            Handle(GeomAdaptor_HSurface) theGASurface,
2743                            gp_Pnt2d&           theNewPoint) {
2744   theNewPoint = theLastPoint;
2745
2746   if ( !IsInsideTanZone( theLastPoint, theTanZoneCenter, theZoneRadius, theGASurface) )
2747     return Standard_False;
2748
2749   Standard_Real aUResolution = theGASurface->UResolution( theZoneRadius );
2750   Standard_Real aVResolution = theGASurface->VResolution( theZoneRadius );
2751
2752   Standard_Real aRadius = ( aUResolution < aVResolution ) ? aUResolution : aVResolution;
2753   gp_Ax22d anAxis( theTanZoneCenter, gp_Dir2d(1, 0), gp_Dir2d(0, 1) );
2754   gp_Circ2d aCircle( anAxis, aRadius );
2755
2756   //
2757   gp_Vec2d aDir( theLastPoint.XY() - theFirstPoint.XY() );
2758   Standard_Real aLength = aDir.Magnitude();
2759   if ( aLength <= gp::Resolution() )
2760     return Standard_False;
2761   gp_Lin2d aLine( theFirstPoint, aDir );
2762
2763   //
2764   Handle(Geom2d_Line) aCLine = new Geom2d_Line( aLine );
2765   Handle(Geom2d_TrimmedCurve) aC1 = new Geom2d_TrimmedCurve( aCLine, 0, aLength );
2766   Handle(Geom2d_Circle) aC2 = new Geom2d_Circle( aCircle );
2767
2768   Standard_Real aTol = aRadius * 0.001;
2769   aTol = ( aTol < Precision::PConfusion() ) ? Precision::PConfusion() : aTol;
2770
2771   Geom2dAPI_InterCurveCurve anIntersector;
2772   anIntersector.Init( aC1, aC2, aTol );
2773
2774   if ( anIntersector.NbPoints() == 0 )
2775     return Standard_False;
2776
2777   Standard_Boolean aFound = Standard_False;
2778   Standard_Real aMinDist = aLength * aLength;
2779   Standard_Integer i = 0;
2780   for ( i = 1; i <= anIntersector.NbPoints(); i++ ) {
2781     gp_Pnt2d aPInt = anIntersector.Point( i );
2782     if ( aPInt.SquareDistance( theFirstPoint ) < aMinDist ) {
2783       if ( ( aPInt.X() >= theUmin ) && ( aPInt.X() <= theUmax ) &&
2784            ( aPInt.Y() >= theVmin ) && ( aPInt.Y() <= theVmax ) ) {
2785         theNewPoint = aPInt;
2786         aFound = Standard_True;
2787       }
2788     }
2789   }
2790
2791   return aFound;
2792 }
2793
2794 // ------------------------------------------------------------------------------------------------
2795 // static function: IsInsideTanZone
2796 // purpose: Check if point is inside a radial tangent zone
2797 // ------------------------------------------------------------------------------------------------
2798 Standard_Boolean IsInsideTanZone(const gp_Pnt2d&     thePoint,
2799                                  const gp_Pnt2d&     theTanZoneCenter,
2800                                  const Standard_Real theZoneRadius,
2801                                  Handle(GeomAdaptor_HSurface) theGASurface) {
2802
2803   Standard_Real aUResolution = theGASurface->UResolution( theZoneRadius );
2804   Standard_Real aVResolution = theGASurface->VResolution( theZoneRadius );
2805   Standard_Real aRadiusSQR = ( aUResolution < aVResolution ) ? aUResolution : aVResolution;
2806   aRadiusSQR *= aRadiusSQR;
2807   if ( thePoint.SquareDistance( theTanZoneCenter ) <= aRadiusSQR )
2808     return Standard_True;
2809   return Standard_False;
2810 }
2811
2812 // ------------------------------------------------------------------------------------------------
2813 // static function: CheckTangentZonesExist
2814 // purpose: Check if tangent zone exists
2815 // ------------------------------------------------------------------------------------------------
2816 Standard_Boolean CheckTangentZonesExist( const Handle(GeomAdaptor_HSurface)& theSurface1,
2817                                         const Handle(GeomAdaptor_HSurface)&  theSurface2 )
2818 {
2819   if ( ( theSurface1->GetType() != GeomAbs_Torus ) ||
2820       ( theSurface2->GetType() != GeomAbs_Torus ) )
2821     return Standard_False;
2822
2823   IntTools_Context aContext;
2824
2825   gp_Torus aTor1 = theSurface1->Torus();
2826   gp_Torus aTor2 = theSurface2->Torus();
2827
2828   if ( aTor1.Location().Distance( aTor2.Location() ) > Precision::Confusion() )
2829     return Standard_False;
2830
2831   if ( ( fabs( aTor1.MajorRadius() - aTor2.MajorRadius() ) > Precision::Confusion() ) ||
2832        ( fabs( aTor1.MinorRadius() - aTor2.MinorRadius() ) > Precision::Confusion() ) )
2833     return Standard_False;
2834
2835   if ( ( aTor1.MajorRadius() < aTor1.MinorRadius() ) ||
2836        ( aTor2.MajorRadius() < aTor2.MinorRadius() ) )
2837     return Standard_False;
2838   return Standard_True;
2839 }
2840
2841 // ------------------------------------------------------------------------------------------------
2842 // static function: ComputeTangentZones
2843 // purpose:
2844 // ------------------------------------------------------------------------------------------------
2845 Standard_Integer ComputeTangentZones( const Handle(GeomAdaptor_HSurface)& theSurface1,
2846                                      const Handle(GeomAdaptor_HSurface)&  theSurface2,
2847                                      const TopoDS_Face&                   theFace1,
2848                                      const TopoDS_Face&                   theFace2,
2849                                      Handle(TColgp_HArray1OfPnt2d)&       theResultOnS1,
2850                                      Handle(TColgp_HArray1OfPnt2d)&       theResultOnS2,
2851                                      Handle(TColStd_HArray1OfReal)&       theResultRadius) {
2852   Standard_Integer aResult = 0;
2853   if ( !CheckTangentZonesExist( theSurface1, theSurface2 ) )
2854     return aResult;
2855
2856   IntTools_Context aContext;
2857
2858   TColgp_SequenceOfPnt2d aSeqResultS1, aSeqResultS2;
2859   TColStd_SequenceOfReal aSeqResultRad;
2860
2861   gp_Torus aTor1 = theSurface1->Torus();
2862   gp_Torus aTor2 = theSurface2->Torus();
2863
2864   gp_Ax2 anax1( aTor1.Location(), aTor1.Axis().Direction() );
2865   gp_Ax2 anax2( aTor2.Location(), aTor2.Axis().Direction() );
2866   Standard_Integer j = 0;
2867
2868   for ( j = 0; j < 2; j++ ) {
2869     Standard_Real aCoef = ( j == 0 ) ? -1 : 1;
2870     Standard_Real aRadius1 = fabs(aTor1.MajorRadius() + aCoef * aTor1.MinorRadius());
2871     Standard_Real aRadius2 = fabs(aTor2.MajorRadius() + aCoef * aTor2.MinorRadius());
2872
2873     gp_Circ aCircle1( anax1, aRadius1 );
2874     gp_Circ aCircle2( anax2, aRadius2 );
2875
2876     // roughly compute radius of tangent zone for perpendicular case
2877     Standard_Real aCriteria = Precision::Confusion() * 0.5;
2878
2879     Standard_Real aT1 = aCriteria;
2880     Standard_Real aT2 = aCriteria;
2881     if ( j == 0 ) {
2882       // internal tangency
2883       Standard_Real aR = ( aRadius1 > aTor2.MinorRadius() ) ? aRadius1 : aTor2.MinorRadius();
2884       //aT1 = aCriteria * aCriteria + aR * aR - ( aR - aCriteria ) * ( aR - aCriteria );
2885       aT1 = 2. * aR * aCriteria;
2886       aT2 = aT1;
2887     }
2888     else {
2889       // external tangency
2890       Standard_Real aRb = ( aRadius1 > aTor2.MinorRadius() ) ? aRadius1 : aTor2.MinorRadius();
2891       Standard_Real aRm = ( aRadius1 < aTor2.MinorRadius() ) ? aRadius1 : aTor2.MinorRadius();
2892       Standard_Real aDelta = aRb - aCriteria;
2893       aDelta *= aDelta;
2894       aDelta -= aRm * aRm;
2895       aDelta /= 2. * (aRb - aRm);
2896       aDelta -= 0.5 * (aRb - aRm);
2897
2898       aT1 = 2. * aRm * (aRm - aDelta);
2899       aT2 = aT1;
2900     }
2901     aCriteria = ( aT1 > aT2) ? aT1 : aT2;
2902     if ( aCriteria > 0 )
2903       aCriteria = sqrt( aCriteria );
2904
2905     if ( aCriteria > 0.5 * aTor1.MinorRadius() ) {
2906       // too big zone -> drop to minimum
2907       aCriteria = Precision::Confusion();
2908     }
2909
2910     GeomAdaptor_Curve aC1( new Geom_Circle(aCircle1) );
2911     GeomAdaptor_Curve aC2( new Geom_Circle(aCircle2) );
2912     Extrema_ExtCC anExtrema(aC1, aC2, 0, 2.*Standard_PI, 0, 2.*Standard_PI,
2913                             Precision::PConfusion(), Precision::PConfusion());
2914
2915     if ( anExtrema.IsDone() ) {
2916
2917       Standard_Integer i = 0;
2918       for ( i = 1; i <= anExtrema.NbExt(); i++ ) {
2919         if ( anExtrema.SquareDistance(i) > aCriteria * aCriteria )
2920           continue;
2921
2922         Extrema_POnCurv P1, P2;
2923         anExtrema.Points( i, P1, P2 );
2924
2925         Standard_Boolean bFoundResult = Standard_True;
2926         gp_Pnt2d pr1, pr2;
2927
2928         Standard_Integer surfit = 0;
2929         for ( surfit = 0; surfit < 2; surfit++ ) {
2930           GeomAPI_ProjectPointOnSurf& aProjector = (surfit == 0) ? aContext.ProjPS(theFace1) : aContext.ProjPS(theFace2);
2931
2932           gp_Pnt aP3d = (surfit == 0) ? P1.Value() : P2.Value();
2933           aProjector.Perform(aP3d);
2934
2935           if(!aProjector.IsDone())
2936             bFoundResult = Standard_False;
2937           else {
2938             if(aProjector.LowerDistance() > aCriteria) {
2939               bFoundResult = Standard_False;
2940             }
2941             else {
2942               Standard_Real foundU = 0, foundV = 0;
2943               aProjector.LowerDistanceParameters(foundU, foundV);
2944               if ( surfit == 0 )
2945                 pr1 = gp_Pnt2d( foundU, foundV );
2946               else
2947                 pr2 = gp_Pnt2d( foundU, foundV );
2948             }
2949           }
2950         }
2951         if ( bFoundResult ) {
2952           aSeqResultS1.Append( pr1 );
2953           aSeqResultS2.Append( pr2 );
2954           aSeqResultRad.Append( aCriteria );
2955
2956           // torus is u and v periodic
2957           const Standard_Real twoPI = Standard_PI + Standard_PI;
2958           Standard_Real arr1tmp[2] = {pr1.X(), pr1.Y()};
2959           Standard_Real arr2tmp[2] = {pr2.X(), pr2.Y()};
2960
2961           // iteration on period bounds
2962           for ( Standard_Integer k1 = 0; k1 < 2; k1++ ) {
2963             Standard_Real aBound = ( k1 == 0 ) ? 0 : twoPI;
2964             Standard_Real aShift = ( k1 == 0 ) ? twoPI : -twoPI;
2965
2966             // iteration on surfaces
2967             for ( Standard_Integer k2 = 0; k2 < 2; k2++ ) {
2968               Standard_Real* arr1 = ( k2 == 0 ) ? arr1tmp : arr2tmp;
2969               Standard_Real* arr2 = ( k2 != 0 ) ? arr1tmp : arr2tmp;
2970               TColgp_SequenceOfPnt2d& aSeqS1 = ( k2 == 0 ) ? aSeqResultS1 : aSeqResultS2;
2971               TColgp_SequenceOfPnt2d& aSeqS2 = ( k2 != 0 ) ? aSeqResultS1 : aSeqResultS2;
2972
2973               if (fabs(arr1[0] - aBound) < Precision::PConfusion()) {
2974                 aSeqS1.Append( gp_Pnt2d( arr1[0] + aShift, arr1[1] ) );
2975                 aSeqS2.Append( gp_Pnt2d( arr2[0], arr2[1] ) );
2976                 aSeqResultRad.Append( aCriteria );
2977               }
2978               if (fabs(arr1[1] - aBound) < Precision::PConfusion()) {
2979                 aSeqS1.Append( gp_Pnt2d( arr1[0], arr1[1] + aShift) );
2980                 aSeqS2.Append( gp_Pnt2d( arr2[0], arr2[1] ) );
2981                 aSeqResultRad.Append( aCriteria );
2982               }
2983             }
2984           } //
2985         }
2986       }
2987     }
2988   }
2989   aResult = aSeqResultRad.Length();
2990
2991   if ( aResult > 0 ) {
2992     theResultOnS1 = new TColgp_HArray1OfPnt2d( 1, aResult );
2993     theResultOnS2 = new TColgp_HArray1OfPnt2d( 1, aResult );
2994     theResultRadius = new TColStd_HArray1OfReal( 1, aResult );
2995
2996     for ( Standard_Integer i = 1 ; i <= aResult; i++ ) {
2997       theResultOnS1->SetValue( i, aSeqResultS1.Value(i) );
2998       theResultOnS2->SetValue( i, aSeqResultS2.Value(i) );
2999       theResultRadius->SetValue( i, aSeqResultRad.Value(i) );
3000     }
3001   }
3002   return aResult;
3003 }
3004
3005 // ------------------------------------------------------------------------------------------------
3006 // static function: AdjustByNeighbour
3007 // purpose:
3008 // ------------------------------------------------------------------------------------------------
3009 gp_Pnt2d AdjustByNeighbour(const gp_Pnt2d&     theaNeighbourPoint,
3010                            const gp_Pnt2d&     theOriginalPoint,
3011                            Handle(GeomAdaptor_HSurface) theGASurface) {
3012
3013   gp_Pnt2d ap1 = theaNeighbourPoint;
3014   gp_Pnt2d ap2 = theOriginalPoint;
3015
3016   if ( theGASurface->IsUPeriodic() ) {
3017     Standard_Real aPeriod     = theGASurface->UPeriod();
3018     gp_Pnt2d aPTest = ap2;
3019     Standard_Real aSqDistMin = 1.e+100;
3020
3021     for ( Standard_Integer pIt = -1; pIt <= 1; pIt++) {
3022       aPTest.SetX( theOriginalPoint.X() + aPeriod * pIt );
3023       Standard_Real dd = ap1.SquareDistance( aPTest );
3024
3025       if ( dd < aSqDistMin ) {
3026         ap2 = aPTest;
3027         aSqDistMin = dd;
3028       }
3029     }
3030   }
3031   if ( theGASurface->IsVPeriodic() ) {
3032     Standard_Real aPeriod     = theGASurface->VPeriod();
3033     gp_Pnt2d aPTest = ap2;
3034     Standard_Real aSqDistMin = 1.e+100;
3035
3036     for ( Standard_Integer pIt = -1; pIt <= 1; pIt++) {
3037       aPTest.SetY( theOriginalPoint.Y() + aPeriod * pIt );
3038       Standard_Real dd = ap1.SquareDistance( aPTest );
3039
3040       if ( dd < aSqDistMin ) {
3041         ap2 = aPTest;
3042         aSqDistMin = dd;
3043       }
3044     }
3045   }
3046   return ap2;
3047 }
3048
3049 // ------------------------------------------------------------------------------------------------
3050 //function: DecompositionOfWLine
3051 // purpose:
3052 // ------------------------------------------------------------------------------------------------
3053 Standard_Boolean DecompositionOfWLine(const Handle(IntPatch_WLine)& theWLine,
3054                                       const Handle(GeomAdaptor_HSurface)&            theSurface1,
3055                                       const Handle(GeomAdaptor_HSurface)&            theSurface2,
3056                                       const TopoDS_Face&                             theFace1,
3057                                       const TopoDS_Face&                             theFace2,
3058                                       const IntTools_LineConstructor&                theLConstructor,
3059                                       const Standard_Boolean                         theAvoidLConstructor,
3060                                       IntPatch_SequenceOfLine&                       theNewLines,
3061                                       Standard_Real&                                 theReachedTol3d) {
3062
3063   Standard_Boolean bRet, bAvoidLineConstructor;
3064   Standard_Integer aNbPnts, aNbParts;
3065   //
3066   bRet=Standard_False;
3067   aNbPnts=theWLine->NbPnts();
3068   bAvoidLineConstructor=theAvoidLConstructor;
3069   //
3070   if(!aNbPnts){
3071     return bRet;
3072   }
3073   if (!bAvoidLineConstructor) {
3074     aNbParts=theLConstructor.NbParts();
3075     if (!aNbParts) {
3076       return bRet;
3077     }
3078   }
3079   //
3080   Standard_Boolean bIsPrevPointOnBoundary, bIsPointOnBoundary, bIsCurrentPointOnBoundary;
3081   Standard_Integer nblines, pit, i, j;
3082   Standard_Real aTol;
3083   TColStd_Array1OfListOfInteger anArrayOfLines(1, aNbPnts);
3084   TColStd_Array1OfInteger       anArrayOfLineType(1, aNbPnts);
3085   TColStd_ListOfInteger aListOfPointIndex;
3086   IntTools_Context aContext;
3087
3088   Handle(TColgp_HArray1OfPnt2d) aTanZoneS1;
3089   Handle(TColgp_HArray1OfPnt2d) aTanZoneS2;
3090   Handle(TColStd_HArray1OfReal) aTanZoneRadius;
3091   Standard_Integer aNbZone = ComputeTangentZones( theSurface1, theSurface2, theFace1, theFace2,
3092                                                  aTanZoneS1, aTanZoneS2, aTanZoneRadius );
3093
3094   //
3095   nblines=0;
3096   aTol=Precision::Confusion();
3097   aTol=0.5*aTol;
3098   bIsPrevPointOnBoundary=Standard_False;
3099   bIsPointOnBoundary=Standard_False;
3100   //
3101   // 1. ...
3102   //
3103   // Points
3104   for(pit = 1; pit <= aNbPnts; ++pit) {
3105     Standard_Boolean bIsOnFirstBoundary, isperiodic;
3106     Standard_Real aResolution, aPeriod, alowerboundary, aupperboundary, U, V;
3107     Standard_Real aParameter, anoffset, anAdjustPar;
3108     Standard_Real umin, umax, vmin, vmax;
3109     //
3110     bIsCurrentPointOnBoundary = Standard_False;
3111     const IntSurf_PntOn2S& aPoint = theWLine->Point(pit);
3112     //
3113     // Surface
3114     for(i = 0; i < 2; ++i) {
3115       Handle(GeomAdaptor_HSurface) aGASurface = (!i) ? theSurface1 : theSurface2;
3116       aGASurface->ChangeSurface().Surface()->Bounds(umin, umax, vmin, vmax);
3117       if(!i) {
3118         aPoint.ParametersOnS1(U, V);
3119       }
3120       else {
3121         aPoint.ParametersOnS2(U, V);
3122       }
3123       // U, V
3124       for(j = 0; j < 2; j++) {
3125         isperiodic = (!j) ? aGASurface->IsUPeriodic() : aGASurface->IsVPeriodic();
3126         if(!isperiodic){
3127           continue;
3128         }
3129         //
3130         if (!j) {
3131           aResolution=aGASurface->UResolution(aTol);
3132           aPeriod=aGASurface->UPeriod();
3133           alowerboundary=umin;
3134           aupperboundary=umax;
3135           aParameter=U;
3136         }
3137         else {
3138           aResolution=aGASurface->VResolution(aTol);
3139           aPeriod=aGASurface->VPeriod();
3140           alowerboundary=vmin;
3141           aupperboundary=vmax;
3142           aParameter=V;
3143         }
3144
3145         anoffset = 0.;
3146         anAdjustPar = AdjustPeriodic(aParameter,
3147                                      alowerboundary,
3148                                      aupperboundary,
3149                                      aPeriod,
3150                                      anoffset);
3151         //
3152         bIsOnFirstBoundary = Standard_True;// ?
3153         bIsPointOnBoundary=
3154           IsPointOnBoundary(anAdjustPar,
3155                             alowerboundary,
3156                             aupperboundary,
3157                             aResolution,
3158                             bIsOnFirstBoundary);
3159         //
3160         if(bIsPointOnBoundary) {
3161           bIsCurrentPointOnBoundary = Standard_True;
3162           break;
3163         }
3164         else {
3165           // check if a point belong to a tangent zone. Begin
3166           Standard_Integer zIt = 0;
3167           for ( zIt = 1; zIt <= aNbZone; zIt++ ) {
3168             gp_Pnt2d aPZone = (i == 0) ? aTanZoneS1->Value(zIt) : aTanZoneS2->Value(zIt);
3169             Standard_Real aZoneRadius = aTanZoneRadius->Value(zIt);
3170
3171             if ( IsInsideTanZone(gp_Pnt2d( U, V ), aPZone, aZoneRadius, aGASurface ) ) {
3172               // set boundary flag to split the curve by a tangent zone
3173               bIsPointOnBoundary = Standard_True;
3174               bIsCurrentPointOnBoundary = Standard_True;
3175               if ( theReachedTol3d < aZoneRadius ) {
3176                 theReachedTol3d = aZoneRadius;
3177               }
3178               break;
3179             }
3180           }
3181         }
3182       }//for(j = 0; j < 2; j++) {
3183
3184       if(bIsCurrentPointOnBoundary){
3185         break;
3186       }
3187     }//for(i = 0; i < 2; ++i) {
3188     //
3189     if((bIsCurrentPointOnBoundary != bIsPrevPointOnBoundary)) {
3190       if(!aListOfPointIndex.IsEmpty()) {
3191         nblines++;
3192         anArrayOfLines.SetValue(nblines, aListOfPointIndex);
3193         anArrayOfLineType.SetValue(nblines, bIsPrevPointOnBoundary);
3194         aListOfPointIndex.Clear();
3195       }
3196       bIsPrevPointOnBoundary = bIsCurrentPointOnBoundary;
3197     }
3198     aListOfPointIndex.Append(pit);
3199   } //for(pit = 1; pit <= aNbPnts; ++pit) {
3200   //
3201   if(!aListOfPointIndex.IsEmpty()) {
3202     nblines++;
3203     anArrayOfLines.SetValue(nblines, aListOfPointIndex);
3204     anArrayOfLineType.SetValue(nblines, bIsPrevPointOnBoundary);
3205     aListOfPointIndex.Clear();
3206   }
3207   //
3208   if(nblines<=1) {
3209     return bRet; //Standard_False;
3210   }
3211   //
3212   //
3213   // 2. Correct wlines.begin
3214   TColStd_Array1OfListOfInteger anArrayOfLineEnds(1, nblines);
3215   Handle(IntSurf_LineOn2S) aSeqOfPntOn2S = new IntSurf_LineOn2S();
3216   //
3217   for(i = 1; i <= nblines; i++) {
3218     if(anArrayOfLineType.Value(i) != 0) {
3219       continue;
3220     }
3221     const TColStd_ListOfInteger& aListOfIndex = anArrayOfLines.Value(i);
3222     if(aListOfIndex.Extent() < 2) {
3223       continue;
3224     }
3225     TColStd_ListOfInteger aListOfFLIndex;
3226
3227     for(j = 0; j < 2; j++) {
3228       Standard_Integer aneighbourindex = (j == 0) ? (i - 1) : (i + 1);
3229
3230       if((aneighbourindex < 1) || (aneighbourindex > nblines))
3231         continue;
3232
3233       if(anArrayOfLineType.Value(aneighbourindex) == 0)
3234         continue;
3235       const TColStd_ListOfInteger& aNeighbour = anArrayOfLines.Value(aneighbourindex);
3236       Standard_Integer anIndex = (j == 0) ? aNeighbour.Last() : aNeighbour.First();
3237       const IntSurf_PntOn2S& aPoint = theWLine->Point(anIndex);
3238
3239       IntSurf_PntOn2S aNewP = aPoint;
3240
3241       for(Standard_Integer surfit = 0; surfit < 2; surfit++) {
3242
3243         Handle(GeomAdaptor_HSurface) aGASurface = (surfit == 0) ? theSurface1 : theSurface2;
3244         Standard_Real umin=0., umax=0., vmin=0., vmax=0.;
3245         aGASurface->ChangeSurface().Surface()->Bounds(umin, umax, vmin, vmax);
3246         Standard_Real U=0., V=0.;
3247
3248         if(surfit == 0)
3249           aNewP.ParametersOnS1(U, V);
3250         else
3251           aNewP.ParametersOnS2(U, V);
3252         Standard_Integer nbboundaries = 0;
3253
3254         Standard_Boolean bIsNearBoundary = Standard_False;
3255         Standard_Integer aZoneIndex = 0;
3256         Standard_Integer bIsUBoundary = Standard_False; // use if nbboundaries == 1
3257         Standard_Integer bIsFirstBoundary = Standard_False; // use if nbboundaries == 1
3258
3259
3260         for(Standard_Integer parit = 0; parit < 2; parit++) {
3261           Standard_Boolean isperiodic = (parit == 0) ? aGASurface->IsUPeriodic() : aGASurface->IsVPeriodic();
3262
3263           Standard_Real aResolution = (parit == 0) ? aGASurface->UResolution(aTol) : aGASurface->VResolution(aTol);
3264           Standard_Real alowerboundary = (parit == 0) ? umin : vmin;
3265           Standard_Real aupperboundary = (parit == 0) ? umax : vmax;
3266
3267           Standard_Real aParameter = (parit == 0) ? U : V;
3268           Standard_Boolean bIsOnFirstBoundary = Standard_True;
3269
3270           if(!isperiodic) {
3271             bIsPointOnBoundary=
3272               IsPointOnBoundary(aParameter, alowerboundary, aupperboundary, aResolution, bIsOnFirstBoundary);
3273             if(bIsPointOnBoundary) {
3274               bIsUBoundary = (parit == 0);
3275               bIsFirstBoundary = bIsOnFirstBoundary;
3276               nbboundaries++;
3277             }
3278           }
3279           else {
3280             Standard_Real aPeriod     = (parit == 0) ? aGASurface->UPeriod() : aGASurface->VPeriod();
3281             Standard_Real anoffset = 0.;
3282             Standard_Real anAdjustPar = AdjustPeriodic(aParameter, alowerboundary, aupperboundary, aPeriod, anoffset);
3283
3284             bIsPointOnBoundary=
3285               IsPointOnBoundary(anAdjustPar, alowerboundary, aupperboundary, aResolution, bIsOnFirstBoundary);
3286             if(bIsPointOnBoundary) {
3287               bIsUBoundary = (parit == 0);
3288               bIsFirstBoundary = bIsOnFirstBoundary;
3289               nbboundaries++;
3290             }
3291             else {
3292               //check neighbourhood of boundary
3293               Standard_Real anEpsilon = aResolution * 100.;
3294               Standard_Real aPart = ( aupperboundary - alowerboundary ) * 0.1;
3295               anEpsilon = ( anEpsilon > aPart ) ? aPart : anEpsilon;
3296
3297               bIsNearBoundary = IsPointOnBoundary(anAdjustPar, alowerboundary, aupperboundary,
3298                                                   anEpsilon, bIsOnFirstBoundary);
3299
3300             }
3301           }
3302         }
3303
3304         // check if a point belong to a tangent zone. Begin
3305         for ( Standard_Integer zIt = 1; zIt <= aNbZone; zIt++ ) {
3306           gp_Pnt2d aPZone = (surfit == 0) ? aTanZoneS1->Value(zIt) : aTanZoneS2->Value(zIt);
3307           Standard_Real aZoneRadius = aTanZoneRadius->Value(zIt);
3308
3309           Standard_Integer aneighbourpointindex1 = (j == 0) ? aListOfIndex.First() : aListOfIndex.Last();
3310           const IntSurf_PntOn2S& aNeighbourPoint = theWLine->Point(aneighbourpointindex1);
3311           Standard_Real nU1, nV1;
3312
3313           if(surfit == 0)
3314             aNeighbourPoint.ParametersOnS1(nU1, nV1);
3315           else
3316             aNeighbourPoint.ParametersOnS2(nU1, nV1);
3317           gp_Pnt2d ap1(nU1, nV1);
3318           gp_Pnt2d ap2 = AdjustByNeighbour( ap1, gp_Pnt2d( U, V ), aGASurface );
3319
3320
3321           if ( IsInsideTanZone( ap2, aPZone, aZoneRadius, aGASurface ) ) {
3322             aZoneIndex = zIt;
3323             bIsNearBoundary = Standard_True;
3324             if ( theReachedTol3d < aZoneRadius ) {
3325               theReachedTol3d = aZoneRadius;
3326             }
3327           }
3328         }
3329         // check if a point belong to a tangent zone. End
3330         Standard_Boolean bComputeLineEnd = Standard_False;
3331
3332         if(nbboundaries == 2) {
3333           //xf
3334           //bComputeLineEnd = Standard_True;
3335           //xt
3336         }
3337         else if(nbboundaries == 1) {
3338           Standard_Boolean isperiodic = (bIsUBoundary) ? aGASurface->IsUPeriodic() : aGASurface->IsVPeriodic();
3339
3340           if(isperiodic) {
3341             Standard_Real alowerboundary = (bIsUBoundary) ? umin : vmin;
3342             Standard_Real aupperboundary = (bIsUBoundary) ? umax : vmax;
3343             Standard_Real aPeriod     = (bIsUBoundary) ? aGASurface->UPeriod() : aGASurface->VPeriod();
3344             Standard_Real aParameter = (bIsUBoundary) ? U : V;
3345             Standard_Real anoffset = 0.;
3346             Standard_Real anAdjustPar = AdjustPeriodic(aParameter, alowerboundary, aupperboundary, aPeriod, anoffset);
3347
3348             Standard_Real adist = (bIsFirstBoundary) ? fabs(anAdjustPar - alowerboundary) : fabs(anAdjustPar - aupperboundary);
3349             Standard_Real anotherPar = (bIsFirstBoundary) ? (aupperboundary - adist) : (alowerboundary + adist);
3350             anotherPar += anoffset;
3351             Standard_Integer aneighbourpointindex = (j == 0) ? aListOfIndex.First() : aListOfIndex.Last();
3352             const IntSurf_PntOn2S& aNeighbourPoint = theWLine->Point(aneighbourpointindex);
3353             Standard_Real nU1, nV1;
3354
3355             if(surfit == 0)
3356               aNeighbourPoint.ParametersOnS1(nU1, nV1);
3357             else
3358               aNeighbourPoint.ParametersOnS2(nU1, nV1);
3359
3360             Standard_Real adist1 = (bIsUBoundary) ? fabs(nU1 - U) : fabs(nV1 - V);
3361             Standard_Real adist2 = (bIsUBoundary) ? fabs(nU1 - anotherPar) : fabs(nV1 - anotherPar);
3362             bComputeLineEnd = Standard_True;
3363             Standard_Boolean bCheckAngle1 = Standard_False;
3364             Standard_Boolean bCheckAngle2 = Standard_False;
3365             gp_Vec2d aNewVec;
3366             Standard_Real anewU = (bIsUBoundary) ? anotherPar : U;
3367             Standard_Real anewV = (bIsUBoundary) ? V : anotherPar;
3368
3369             if(((adist1 - adist2) > Precision::PConfusion()) &&
3370                (adist2 < (aPeriod / 4.))) {
3371               bCheckAngle1 = Standard_True;
3372               aNewVec = gp_Vec2d(gp_Pnt2d(nU1, nV1), gp_Pnt2d(anewU, anewV));
3373
3374               if(aNewVec.SquareMagnitude() < (gp::Resolution() * gp::Resolution())) {
3375                 aNewP.SetValue((surfit == 0), anewU, anewV);
3376                 bCheckAngle1 = Standard_False;
3377               }
3378             }
3379             else if(adist1 < (aPeriod / 4.)) {
3380               bCheckAngle2 = Standard_True;
3381               aNewVec = gp_Vec2d(gp_Pnt2d(nU1, nV1), gp_Pnt2d(U, V));
3382
3383               if(aNewVec.SquareMagnitude() < (gp::Resolution() * gp::Resolution())) {
3384                 bCheckAngle2 = Standard_False;
3385               }
3386             }
3387
3388             if(bCheckAngle1 || bCheckAngle2) {
3389               // assume there are at least two points in line (see "if" above)
3390               Standard_Integer anindexother = aneighbourpointindex;
3391
3392               while((anindexother <= aListOfIndex.Last()) && (anindexother >= aListOfIndex.First())) {
3393                 anindexother = (j == 0) ? (anindexother + 1) : (anindexother - 1);
3394                 const IntSurf_PntOn2S& aPrevNeighbourPoint = theWLine->Point(anindexother);
3395                 Standard_Real nU2, nV2;
3396
3397                 if(surfit == 0)
3398                   aPrevNeighbourPoint.ParametersOnS1(nU2, nV2);
3399                 else
3400                   aPrevNeighbourPoint.ParametersOnS2(nU2, nV2);
3401                 gp_Vec2d aVecOld(gp_Pnt2d(nU2, nV2), gp_Pnt2d(nU1, nV1));
3402
3403                 if(aVecOld.SquareMagnitude() <= (gp::Resolution() * gp::Resolution())) {
3404                   continue;
3405                 }
3406                 else {
3407                   Standard_Real anAngle = aNewVec.Angle(aVecOld);
3408
3409                   if((fabs(anAngle) < (Standard_PI * 0.25)) && (aNewVec.Dot(aVecOld) > 0.)) {
3410
3411                     if(bCheckAngle1) {
3412                       Standard_Real U1, U2, V1, V2;
3413                       IntSurf_PntOn2S atmppoint = aNewP;
3414                       atmppoint.SetValue((surfit == 0), anewU, anewV);
3415                       atmppoint.Parameters(U1, V1, U2, V2);
3416                       gp_Pnt P1 = theSurface1->Value(U1, V1);
3417                       gp_Pnt P2 = theSurface2->Value(U2, V2);
3418                       gp_Pnt P0 = aPoint.Value();
3419
3420                       if(P0.IsEqual(P1, aTol) &&
3421                          P0.IsEqual(P2, aTol) &&
3422                          P1.IsEqual(P2, aTol)) {
3423                         bComputeLineEnd = Standard_False;
3424                         aNewP.SetValue((surfit == 0), anewU, anewV);
3425                       }
3426                     }
3427
3428                     if(bCheckAngle2) {
3429                       bComputeLineEnd = Standard_False;
3430                     }
3431                   }
3432                   break;
3433                 }
3434               } // end while(anindexother...)
3435             }
3436           }
3437         }
3438         else if ( bIsNearBoundary ) {
3439           bComputeLineEnd = Standard_True;
3440         }
3441
3442         if(bComputeLineEnd) {
3443
3444           gp_Pnt2d anewpoint;
3445           Standard_Boolean found = Standard_False;
3446
3447           if ( bIsNearBoundary ) {
3448             // re-compute point near natural boundary or near tangent zone
3449             Standard_Real u1, v1, u2, v2;
3450             aNewP.Parameters( u1, v1, u2, v2 );
3451             if(surfit == 0)
3452               anewpoint = gp_Pnt2d( u1, v1 );
3453             else
3454               anewpoint = gp_Pnt2d( u2, v2 );
3455
3456             Standard_Integer aneighbourpointindex1 = (j == 0) ? aListOfIndex.First() : aListOfIndex.Last();
3457             const IntSurf_PntOn2S& aNeighbourPoint = theWLine->Point(aneighbourpointindex1);
3458             Standard_Real nU1, nV1;
3459
3460             if(surfit == 0)
3461               aNeighbourPoint.ParametersOnS1(nU1, nV1);
3462             else
3463               aNeighbourPoint.ParametersOnS2(nU1, nV1);
3464             gp_Pnt2d ap1(nU1, nV1);
3465             gp_Pnt2d ap2;
3466
3467
3468             if ( aZoneIndex ) {
3469               // exclude point from a tangent zone
3470               anewpoint = AdjustByNeighbour( ap1, anewpoint, aGASurface );
3471               gp_Pnt2d aPZone = (surfit == 0) ? aTanZoneS1->Value(aZoneIndex) : aTanZoneS2->Value(aZoneIndex);
3472               Standard_Real aZoneRadius = aTanZoneRadius->Value(aZoneIndex);
3473
3474               if ( FindPoint(ap1, anewpoint, umin, umax, vmin, vmax,
3475                              aPZone, aZoneRadius, aGASurface, ap2) ) {
3476                 anewpoint = ap2;
3477                 found = Standard_True;
3478               }
3479             }
3480             else if ( aGASurface->IsUPeriodic() || aGASurface->IsVPeriodic() ) {
3481               // re-compute point near boundary if shifted on a period
3482               ap2 = AdjustByNeighbour( ap1, anewpoint, aGASurface );
3483
3484               if ( ( ap2.X() < umin ) || ( ap2.X() > umax ) ||
3485                   ( ap2.Y() < vmin ) || ( ap2.Y() > vmax ) ) {
3486                 found = FindPoint(ap1, ap2, umin, umax, vmin, vmax, anewpoint);
3487               }
3488               else {
3489                 anewpoint = ap2;
3490                 aNewP.SetValue( (surfit == 0), anewpoint.X(), anewpoint.Y() );
3491               }
3492             }
3493           }
3494           else {
3495
3496             Standard_Integer aneighbourpointindex1 = (j == 0) ? aListOfIndex.First() : aListOfIndex.Last();
3497             const IntSurf_PntOn2S& aNeighbourPoint = theWLine->Point(aneighbourpointindex1);
3498             Standard_Real nU1, nV1;
3499
3500             if(surfit == 0)
3501               aNeighbourPoint.ParametersOnS1(nU1, nV1);
3502             else
3503               aNeighbourPoint.ParametersOnS2(nU1, nV1);
3504             gp_Pnt2d ap1(nU1, nV1);
3505             gp_Pnt2d ap2(nU1, nV1);
3506             Standard_Integer aneighbourpointindex2 = aneighbourpointindex1;
3507
3508             while((aneighbourpointindex2 <= aListOfIndex.Last()) && (aneighbourpointindex2 >= aListOfIndex.First())) {
3509               aneighbourpointindex2 = (j == 0) ? (aneighbourpointindex2 + 1) : (aneighbourpointindex2 - 1);
3510               const IntSurf_PntOn2S& aPrevNeighbourPoint = theWLine->Point(aneighbourpointindex2);
3511               Standard_Real nU2, nV2;
3512
3513               if(surfit == 0)
3514                 aPrevNeighbourPoint.ParametersOnS1(nU2, nV2);
3515               else
3516                 aPrevNeighbourPoint.ParametersOnS2(nU2, nV2);
3517               ap2.SetX(nU2);
3518               ap2.SetY(nV2);
3519
3520               if(ap1.SquareDistance(ap2) > (gp::Resolution() * gp::Resolution())) {
3521                 break;
3522               }
3523             }
3524             found = FindPoint(ap2, ap1, umin, umax, vmin, vmax, anewpoint);
3525           }
3526
3527           if(found) {
3528             // check point
3529             Standard_Real aCriteria = BRep_Tool::Tolerance(theFace1) + BRep_Tool::Tolerance(theFace2);
3530             GeomAPI_ProjectPointOnSurf& aProjector = (surfit == 0) ? aContext.ProjPS(theFace2) : aContext.ProjPS(theFace1);
3531             Handle(GeomAdaptor_HSurface) aSurface = (surfit == 0) ? theSurface1 : theSurface2;
3532
3533             Handle(GeomAdaptor_HSurface) aSurfaceOther = (surfit == 0) ? theSurface2 : theSurface1;
3534
3535             gp_Pnt aP3d = aSurface->Value(anewpoint.X(), anewpoint.Y());
3536             aProjector.Perform(aP3d);
3537
3538             if(aProjector.IsDone()) {
3539               if(aProjector.LowerDistance() < aCriteria) {
3540                 Standard_Real foundU = U, foundV = V;
3541                 aProjector.LowerDistanceParameters(foundU, foundV);
3542
3543                 //Correction of projected coordinates. Begin
3544                 //Note, it may be shifted on a period
3545                 Standard_Integer aneindex1 = (j == 0) ? aListOfIndex.First() : aListOfIndex.Last();
3546                 const IntSurf_PntOn2S& aNeighbourPoint = theWLine->Point(aneindex1);
3547                 Standard_Real nUn, nVn;
3548
3549                 if(surfit == 0)
3550                   aNeighbourPoint.ParametersOnS2(nUn, nVn);
3551                 else
3552                   aNeighbourPoint.ParametersOnS1(nUn, nVn);
3553                 gp_Pnt2d aNeighbour2d(nUn, nVn);
3554                 gp_Pnt2d anAdjustedPoint = AdjustByNeighbour( aNeighbour2d, gp_Pnt2d(foundU, foundV), aSurfaceOther );
3555                 foundU = anAdjustedPoint.X();
3556                 foundV = anAdjustedPoint.Y();
3557
3558                 if ( ( anAdjustedPoint.X() < umin ) && ( anAdjustedPoint.X() > umax ) &&
3559                     ( anAdjustedPoint.Y() < vmin ) && ( anAdjustedPoint.Y() > vmax ) ) {
3560                   // attempt to roughly re-compute point
3561                   foundU = ( foundU < umin ) ? umin : foundU;
3562                   foundU = ( foundU > umax ) ? umax : foundU;
3563                   foundV = ( foundV < vmin ) ? vmin : foundV;
3564                   foundV = ( foundV > vmax ) ? vmax : foundV;
3565
3566                   GeomAPI_ProjectPointOnSurf& aProjector2 = (surfit == 0) ? aContext.ProjPS(theFace1) : aContext.ProjPS(theFace2);
3567
3568                   aP3d = aSurfaceOther->Value(foundU, foundV);
3569                   aProjector2.Perform(aP3d);
3570
3571                   if(aProjector2.IsDone()) {
3572                     if(aProjector2.LowerDistance() < aCriteria) {
3573                       Standard_Real foundU2 = anewpoint.X(), foundV2 = anewpoint.Y();
3574                       aProjector2.LowerDistanceParameters(foundU2, foundV2);
3575                       anewpoint.SetX(foundU2);
3576                       anewpoint.SetY(foundV2);
3577                     }
3578                   }
3579                 }
3580                 //Correction of projected coordinates. End
3581
3582                 if(surfit == 0)
3583                   aNewP.SetValue(aP3d, anewpoint.X(), anewpoint.Y(), foundU, foundV);
3584                 else
3585                   aNewP.SetValue(aP3d, foundU, foundV, anewpoint.X(), anewpoint.Y());
3586               }
3587             }
3588           }
3589         }
3590       }
3591       aSeqOfPntOn2S->Add(aNewP);
3592       aListOfFLIndex.Append(aSeqOfPntOn2S->NbPoints());
3593     }
3594     anArrayOfLineEnds.SetValue(i, aListOfFLIndex);
3595   }
3596   // Correct wlines.end
3597
3598   // Split wlines.begin
3599   Standard_Integer nbiter;
3600   //
3601   nbiter=1;
3602   if (!bAvoidLineConstructor) {
3603     nbiter=theLConstructor.NbParts();
3604   }
3605   //
3606   for(j = 1; j <= nbiter; ++j) {
3607     Standard_Real fprm, lprm;
3608     Standard_Integer ifprm, ilprm;
3609     //
3610     if(bAvoidLineConstructor) {
3611       ifprm = 1;
3612       ilprm = theWLine->NbPnts();
3613     }
3614     else {
3615       theLConstructor.Part(j, fprm, lprm);
3616       ifprm = (Standard_Integer)fprm;
3617       ilprm = (Standard_Integer)lprm;
3618     }
3619
3620     Handle(IntSurf_LineOn2S) aLineOn2S = new IntSurf_LineOn2S();
3621     //
3622     for(i = 1; i <= nblines; i++) {
3623       if(anArrayOfLineType.Value(i) != 0) {
3624         continue;
3625       }
3626       const TColStd_ListOfInteger& aListOfIndex = anArrayOfLines.Value(i);
3627
3628       if(aListOfIndex.Extent() < 2) {
3629         continue;
3630       }
3631       const TColStd_ListOfInteger& aListOfFLIndex = anArrayOfLineEnds.Value(i);
3632       Standard_Boolean bhasfirstpoint = (aListOfFLIndex.Extent() == 2);
3633       Standard_Boolean bhaslastpoint = (aListOfFLIndex.Extent() == 2);
3634
3635       if(!bhasfirstpoint && !aListOfFLIndex.IsEmpty()) {
3636         bhasfirstpoint = (i != 1);
3637       }
3638
3639       if(!bhaslastpoint && !aListOfFLIndex.IsEmpty()) {
3640         bhaslastpoint = (i != nblines);
3641       }
3642       Standard_Boolean bIsFirstInside = ((ifprm >= aListOfIndex.First()) && (ifprm <= aListOfIndex.Last()));
3643       Standard_Boolean bIsLastInside =  ((ilprm >= aListOfIndex.First()) && (ilprm <= aListOfIndex.Last()));
3644
3645       if(!bIsFirstInside && !bIsLastInside) {
3646         if((ifprm < aListOfIndex.First()) && (ilprm > aListOfIndex.Last())) {
3647           // append whole line, and boundaries if neccesary
3648           if(bhasfirstpoint) {
3649             const IntSurf_PntOn2S& aP = aSeqOfPntOn2S->Value(aListOfFLIndex.First());
3650             aLineOn2S->Add(aP);
3651           }
3652           TColStd_ListIteratorOfListOfInteger anIt(aListOfIndex);
3653
3654           for(; anIt.More(); anIt.Next()) {
3655             const IntSurf_PntOn2S& aP = theWLine->Point(anIt.Value());
3656             aLineOn2S->Add(aP);
3657           }
3658
3659           if(bhaslastpoint) {
3660             const IntSurf_PntOn2S& aP = aSeqOfPntOn2S->Value(aListOfFLIndex.Last());
3661             aLineOn2S->Add(aP);
3662           }
3663
3664           // check end of split line (end is almost always)
3665           Standard_Integer aneighbour = i + 1;
3666           Standard_Boolean bIsEndOfLine = Standard_True;
3667
3668           if(aneighbour <= nblines) {
3669             const TColStd_ListOfInteger& aListOfNeighbourIndex = anArrayOfLines.Value(aneighbour);
3670
3671             if((anArrayOfLineType.Value(aneighbour) != 0) &&
3672                (aListOfNeighbourIndex.IsEmpty())) {
3673               bIsEndOfLine = Standard_False;
3674             }
3675           }
3676
3677           if(bIsEndOfLine) {
3678             if(aLineOn2S->NbPoints() > 1) {
3679               Handle(IntPatch_WLine) aNewWLine =
3680                 new IntPatch_WLine(aLineOn2S, Standard_False);
3681               theNewLines.Append(aNewWLine);
3682             }
3683             aLineOn2S = new IntSurf_LineOn2S();
3684           }
3685         }
3686         continue;
3687       }
3688       // end if(!bIsFirstInside && !bIsLastInside)
3689
3690       if(bIsFirstInside && bIsLastInside) {
3691         // append inside points between ifprm and ilprm
3692         TColStd_ListIteratorOfListOfInteger anIt(aListOfIndex);
3693
3694         for(; anIt.More(); anIt.Next()) {
3695           if((anIt.Value() < ifprm) || (anIt.Value() > ilprm))
3696             continue;
3697           const IntSurf_PntOn2S& aP = theWLine->Point(anIt.Value());
3698           aLineOn2S->Add(aP);
3699         }
3700       }
3701       else {
3702
3703         if(bIsFirstInside) {
3704           // append points from ifprm to last point + boundary point
3705           TColStd_ListIteratorOfListOfInteger anIt(aListOfIndex);
3706
3707           for(; anIt.More(); anIt.Next()) {
3708             if(anIt.Value() < ifprm)
3709               continue;
3710             const IntSurf_PntOn2S& aP = theWLine->Point(anIt.Value());
3711             aLineOn2S->Add(aP);
3712           }
3713
3714           if(bhaslastpoint) {
3715             const IntSurf_PntOn2S& aP = aSeqOfPntOn2S->Value(aListOfFLIndex.Last());
3716             aLineOn2S->Add(aP);
3717           }
3718           // check end of split line (end is almost always)
3719           Standard_Integer aneighbour = i + 1;
3720           Standard_Boolean bIsEndOfLine = Standard_True;
3721
3722           if(aneighbour <= nblines) {
3723             const TColStd_ListOfInteger& aListOfNeighbourIndex = anArrayOfLines.Value(aneighbour);
3724
3725             if((anArrayOfLineType.Value(aneighbour) != 0) &&
3726                (aListOfNeighbourIndex.IsEmpty())) {
3727               bIsEndOfLine = Standard_False;
3728             }
3729           }
3730
3731           if(bIsEndOfLine) {
3732             if(aLineOn2S->NbPoints() > 1) {
3733               Handle(IntPatch_WLine) aNewWLine =
3734                 new IntPatch_WLine(aLineOn2S, Standard_False);
3735               theNewLines.Append(aNewWLine);
3736             }
3737             aLineOn2S = new IntSurf_LineOn2S();
3738           }
3739         }
3740         // end if(bIsFirstInside)
3741
3742         if(bIsLastInside) {
3743           // append points from first boundary point to ilprm
3744           if(bhasfirstpoint) {
3745             const IntSurf_PntOn2S& aP = aSeqOfPntOn2S->Value(aListOfFLIndex.First());
3746             aLineOn2S->Add(aP);
3747           }
3748           TColStd_ListIteratorOfListOfInteger anIt(aListOfIndex);
3749
3750           for(; anIt.More(); anIt.Next()) {
3751             if(anIt.Value() > ilprm)
3752               continue;
3753             const IntSurf_PntOn2S& aP = theWLine->Point(anIt.Value());
3754             aLineOn2S->Add(aP);
3755           }
3756         }
3757         //end if(bIsLastInside)
3758       }
3759     }
3760
3761     if(aLineOn2S->NbPoints() > 1) {
3762       Handle(IntPatch_WLine) aNewWLine =
3763         new IntPatch_WLine(aLineOn2S, Standard_False);
3764       theNewLines.Append(aNewWLine);
3765     }
3766   }
3767   // Split wlines.end
3768
3769   return Standard_True;
3770 }
3771
3772 // ------------------------------------------------------------------------------------------------
3773 // static function: ParameterOutOfBoundary
3774 // purpose:         Computes a new parameter for given curve. The corresponding 2d points
3775 //                  does not lay on any boundary of given faces
3776 // ------------------------------------------------------------------------------------------------
3777 Standard_Boolean ParameterOutOfBoundary(const Standard_Real       theParameter,
3778                                         const Handle(Geom_Curve)& theCurve,
3779                                         const TopoDS_Face&        theFace1,
3780                                         const TopoDS_Face&        theFace2,
3781                                         const Standard_Real       theOtherParameter,
3782                                         const Standard_Boolean    bIncreasePar,
3783                                         Standard_Real&            theNewParameter) {
3784   Standard_Boolean bIsComputed = Standard_False;
3785   theNewParameter = theParameter;
3786
3787   IntTools_Context aContext;
3788   Standard_Real acurpar = theParameter;
3789   TopAbs_State aState = TopAbs_ON;
3790   Standard_Integer iter = 0;
3791   Standard_Real asumtol = BRep_Tool::Tolerance(theFace1) + BRep_Tool::Tolerance(theFace2);
3792   Standard_Real adelta = asumtol * 0.1;
3793   adelta = (adelta < Precision::Confusion()) ? Precision::Confusion() : adelta;
3794   Handle(Geom_Surface) aSurf1 = BRep_Tool::Surface(theFace1);
3795   Handle(Geom_Surface) aSurf2 = BRep_Tool::Surface(theFace2);
3796
3797   Standard_Real u1, u2, v1, v2;
3798
3799   GeomAPI_ProjectPointOnSurf aPrj1;
3800   aSurf1->Bounds(u1, u2, v1, v2);
3801   aPrj1.Init(aSurf1, u1, u2, v1, v2);
3802
3803   GeomAPI_ProjectPointOnSurf aPrj2;
3804   aSurf2->Bounds(u1, u2, v1, v2);
3805   aPrj2.Init(aSurf2, u1, u2, v1, v2);
3806
3807   while(aState == TopAbs_ON) {
3808     if(bIncreasePar)
3809       acurpar += adelta;
3810     else
3811       acurpar -= adelta;
3812     gp_Pnt aPCurrent = theCurve->Value(acurpar);
3813     aPrj1.Perform(aPCurrent);
3814     Standard_Real U=0., V=0.;
3815
3816     if(aPrj1.IsDone()) {
3817       aPrj1.LowerDistanceParameters(U, V);
3818       aState = aContext.StatePointFace(theFace1, gp_Pnt2d(U, V));
3819     }
3820
3821     if(aState != TopAbs_ON) {
3822       aPrj2.Perform(aPCurrent);
3823
3824       if(aPrj2.IsDone()) {
3825         aPrj2.LowerDistanceParameters(U, V);
3826         aState = aContext.StatePointFace(theFace2, gp_Pnt2d(U, V));
3827       }
3828     }
3829
3830     if(iter > 11) {
3831       break;
3832     }
3833     iter++;
3834   }
3835
3836   if(iter <= 11) {
3837     theNewParameter = acurpar;
3838     bIsComputed = Standard_True;
3839
3840     if(bIncreasePar) {
3841       if(acurpar >= theOtherParameter)
3842         theNewParameter = theOtherParameter;
3843     }
3844     else {
3845       if(acurpar <= theOtherParameter)
3846         theNewParameter = theOtherParameter;
3847     }
3848   }
3849   return bIsComputed;
3850 }
3851
3852 Standard_Boolean IsCurveValid(Handle(Geom2d_Curve)& thePCurve)
3853 {
3854   if(thePCurve.IsNull())
3855     return Standard_False;
3856
3857   Standard_Real tolint = 1.e-10;
3858   Geom2dAdaptor_Curve PCA;
3859   IntRes2d_Domain PCD;
3860   Geom2dInt_GInter PCI;
3861
3862   Standard_Real pf = 0., pl = 0.;
3863   gp_Pnt2d pntf, pntl;
3864
3865   if(!thePCurve->IsClosed() && !thePCurve->IsPeriodic()) {
3866     pf = thePCurve->FirstParameter();
3867     pl = thePCurve->LastParameter();
3868     pntf = thePCurve->Value(pf);
3869     pntl = thePCurve->Value(pl);
3870     PCA.Load(thePCurve);
3871     if(!PCA.IsPeriodic()) {
3872       if(PCA.FirstParameter() > pf) pf = PCA.FirstParameter();
3873       if(PCA.LastParameter()  < pl) pl = PCA.LastParameter();
3874     }
3875     PCD.SetValues(pntf,pf,tolint,pntl,pl,tolint);
3876     PCI.Perform(PCA,PCD,tolint,tolint);
3877     if(PCI.IsDone())
3878       if(PCI.NbPoints() > 0) {
3879         return Standard_False;
3880       }
3881   }
3882
3883   return Standard_True;
3884 }
3885
3886 //=======================================================================
3887 //static function : ApproxWithPCurves
3888 //purpose  : for bug 20964 only
3889 //=======================================================================
3890
3891 Standard_Boolean ApproxWithPCurves(const gp_Cylinder& theCyl,
3892                                    const gp_Sphere& theSph)
3893 {
3894   Standard_Boolean bRes = Standard_True;
3895   Standard_Real R1 = theCyl.Radius(), R2 = theSph.Radius();
3896
3897   if(R1 < 2.*R2) return bRes;
3898
3899   gp_Lin anCylAx(theCyl.Axis());
3900
3901   Standard_Real aDist = anCylAx.Distance(theSph.Location());
3902   Standard_Real aDRel = Abs(aDist - R1)/R2;
3903
3904   if(aDRel > .2) return bRes;
3905
3906   Standard_Real par = ElCLib::Parameter(anCylAx, theSph.Location());
3907   gp_Pnt aP = ElCLib::Value(par, anCylAx);
3908   gp_Vec aV(aP, theSph.Location());
3909
3910   Standard_Real dd = aV.Dot(theSph.Position().XDirection());
3911
3912   if(aDist < R1 && dd > 0.) return Standard_False;
3913   if(aDist > R1 && dd < 0.) return Standard_False;
3914
3915
3916   return bRes;
3917 }
3918
3919 //=======================================================================
3920 //function : PerformPlanes
3921 //purpose  :
3922 //=======================================================================
3923
3924 void  PerformPlanes(const Handle(GeomAdaptor_HSurface)& theS1,
3925                     const Handle(GeomAdaptor_HSurface)& theS2,
3926                     const Standard_Real TolAng,
3927                     const Standard_Real TolTang,
3928                     const Standard_Boolean theApprox1,
3929                     const Standard_Boolean theApprox2,
3930                     IntTools_SequenceOfCurves& theSeqOfCurve,
3931                     Standard_Boolean& theTangentFaces)
3932 {
3933
3934   gp_Pln aPln1 = theS1->Surface().Plane();
3935   gp_Pln aPln2 = theS2->Surface().Plane();
3936
3937   IntAna_QuadQuadGeo aPlnInter(aPln1, aPln2, TolAng, TolTang);
3938
3939   if(!aPlnInter.IsDone()) {
3940     theTangentFaces = Standard_False;
3941     return;
3942   }
3943
3944   IntAna_ResultType aResType = aPlnInter.TypeInter();
3945
3946   if(aResType == IntAna_Same) {
3947     theTangentFaces = Standard_True;
3948     return;
3949   }
3950
3951   theTangentFaces = Standard_False;
3952
3953   if(aResType == IntAna_Empty) {
3954     return;
3955   }
3956
3957   gp_Lin aLin = aPlnInter.Line(1);
3958
3959   ProjLib_Plane aProj;
3960
3961   aProj.Init(aPln1);
3962   aProj.Project(aLin);
3963   gp_Lin2d aLin2d1 = aProj.Line();
3964   //
3965   aProj.Init(aPln2);
3966   aProj.Project(aLin);
3967   gp_Lin2d aLin2d2 = aProj.Line();
3968   //
3969   //classify line2d1 relatively first plane
3970   Standard_Real P11, P12;
3971   Standard_Boolean IsCrossed = ClassifyLin2d(theS1, aLin2d1, TolTang, P11, P12);
3972   if(!IsCrossed) return;
3973   //classify line2d2 relatively second plane
3974   Standard_Real P21, P22;
3975   IsCrossed = ClassifyLin2d(theS2, aLin2d2, TolTang, P21, P22);
3976   if(!IsCrossed) return;
3977
3978   //Analysis of parametric intervals: must have common part
3979
3980   if(P21 >= P12) return;
3981   if(P22 <= P11) return;
3982
3983   Standard_Real pmin, pmax;
3984   pmin = Max(P11, P21);
3985   pmax = Min(P12, P22);
3986
3987   if(pmax - pmin <= TolTang) return;
3988
3989   Handle(Geom_Line) aGLin = new Geom_Line(aLin);
3990
3991   IntTools_Curve aCurve;
3992   Handle(Geom_TrimmedCurve) aGTLin = new Geom_TrimmedCurve(aGLin, pmin, pmax);
3993
3994   aCurve.SetCurve(aGTLin);
3995
3996   if(theApprox1) {
3997     Handle(Geom2d_Line) C2d = new Geom2d_Line(aLin2d1);
3998     aCurve.SetFirstCurve2d(new Geom2d_TrimmedCurve(C2d, pmin, pmax));
3999   }
4000   else {
4001     Handle(Geom2d_Curve) H1;
4002     aCurve.SetFirstCurve2d(H1);
4003   }
4004   if(theApprox2) {
4005     Handle(Geom2d_Line) C2d = new Geom2d_Line(aLin2d2);
4006     aCurve.SetSecondCurve2d(new Geom2d_TrimmedCurve(C2d, pmin, pmax));
4007   }
4008   else {
4009     Handle(Geom2d_Curve) H1;
4010     aCurve.SetFirstCurve2d(H1);
4011   }
4012
4013   theSeqOfCurve.Append(aCurve);
4014
4015 }
4016
4017 //=======================================================================
4018 //function : ClassifyLin2d
4019 //purpose  :
4020 //=======================================================================
4021 static inline Standard_Boolean INTER(const Standard_Real d1,
4022                                      const Standard_Real d2,
4023                                      const Standard_Real tol)
4024 {
4025   return (d1 > tol && d2 < -tol) ||
4026          (d1 < -tol && d2 > tol) ||
4027          ((d1 <= tol && d1 >= -tol) && (d2 > tol || d2 < -tol)) ||
4028          ((d2 <= tol && d2 >= -tol) && (d1 > tol || d1 < -tol));
4029 }
4030 static inline Standard_Boolean COINC(const Standard_Real d1,
4031                                      const Standard_Real d2,
4032                                      const Standard_Real tol)
4033 {
4034   return (d1 <= tol && d1 >= -tol) && (d2 <= tol && d2 >= -tol);
4035 }
4036 Standard_Boolean ClassifyLin2d(const Handle(GeomAdaptor_HSurface)& theS,
4037                                const gp_Lin2d& theLin2d,
4038                                const Standard_Real theTol,
4039                                Standard_Real& theP1,
4040                                Standard_Real& theP2)
4041
4042 {
4043   Standard_Real xmin, xmax, ymin, ymax, d1, d2, A, B, C;
4044   Standard_Real par[2];
4045   Standard_Integer nbi = 0;
4046
4047   xmin = theS->Surface().FirstUParameter();
4048   xmax = theS->Surface().LastUParameter();
4049   ymin = theS->Surface().FirstVParameter();
4050   ymax = theS->Surface().LastVParameter();
4051
4052   theLin2d.Coefficients(A, B, C);
4053
4054   //xmin, ymin <-> xmin, ymax
4055   d1 = A*xmin + B*ymin + C;
4056   d2 = A*xmin + B*ymax + C;
4057
4058   if(INTER(d1, d2, theTol)) {
4059     //Intersection with boundary
4060     Standard_Real y = -(C + A*xmin)/B;
4061     par[nbi] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmin, y));
4062     nbi++;
4063   }
4064   else if (COINC(d1, d2, theTol)) {
4065     //Coincidence with boundary
4066     par[0] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmin, ymin));
4067     par[1] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmin, ymax));
4068     nbi = 2;
4069   }
4070
4071   if(nbi == 2) {
4072
4073     if(fabs(par[0]-par[1]) > theTol) {
4074       theP1 = Min(par[0], par[1]);
4075       theP2 = Max(par[0], par[1]);
4076       return Standard_True;
4077     }
4078     else return Standard_False;
4079
4080   }
4081
4082   //xmin, ymax <-> xmax, ymax
4083   d1 = d2;
4084   d2 = A*xmax + B*ymax + C;
4085
4086   if(d1 > theTol || d1 < -theTol) {//to avoid checking of
4087                                    //coincidence with the same point
4088     if(INTER(d1, d2, theTol)) {
4089       Standard_Real x = -(C + B*ymax)/A;
4090       par[nbi] = ElCLib::Parameter(theLin2d, gp_Pnt2d(x, ymax));
4091       nbi++;
4092     }
4093     else if (COINC(d1, d2, theTol)) {
4094       par[0] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmin, ymax));
4095       par[1] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmax, ymax));
4096       nbi = 2;
4097     }
4098   }
4099
4100   if(nbi == 2) {
4101
4102     if(fabs(par[0]-par[1]) > theTol) {
4103       theP1 = Min(par[0], par[1]);
4104       theP2 = Max(par[0], par[1]);
4105       return Standard_True;
4106     }
4107     else return Standard_False;
4108
4109   }
4110
4111   //xmax, ymax <-> xmax, ymin
4112   d1 = d2;
4113   d2 = A*xmax + B*ymin + C;
4114
4115   if(d1 > theTol || d1 < -theTol) {
4116     if(INTER(d1, d2, theTol)) {
4117       Standard_Real y = -(C + A*xmax)/B;
4118       par[nbi] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmax, y));
4119       nbi++;
4120     }
4121     else if (COINC(d1, d2, theTol)) {
4122       par[0] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmax, ymax));
4123       par[1] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmax, ymin));
4124       nbi = 2;
4125     }
4126   }
4127
4128   if(nbi == 2) {
4129     if(fabs(par[0]-par[1]) > theTol) {
4130       theP1 = Min(par[0], par[1]);
4131       theP2 = Max(par[0], par[1]);
4132       return Standard_True;
4133     }
4134     else return Standard_False;
4135   }
4136
4137   //xmax, ymin <-> xmin, ymin
4138   d1 = d2;
4139   d2 = A*xmin + B*ymin + C;
4140
4141   if(d1 > theTol || d1 < -theTol) {
4142     if(INTER(d1, d2, theTol)) {
4143       Standard_Real x = -(C + B*ymin)/A;
4144       par[nbi] = ElCLib::Parameter(theLin2d, gp_Pnt2d(x, ymin));
4145       nbi++;
4146     }
4147     else if (COINC(d1, d2, theTol)) {
4148       par[0] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmax, ymin));
4149       par[1] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmin, ymin));
4150       nbi = 2;
4151     }
4152   }
4153
4154   if(nbi == 2) {
4155     if(fabs(par[0]-par[1]) > theTol) {
4156       theP1 = Min(par[0], par[1]);
4157       theP2 = Max(par[0], par[1]);
4158       return Standard_True;
4159     }
4160     else return Standard_False;
4161   }
4162
4163   return Standard_False;
4164
4165 }
4166 //
4167 //modified by NIZNHY-PKV Thu Sep 23 08:37:32 2010f
4168 //=======================================================================
4169 //function : ApproxParameters
4170 //purpose  :
4171 //=======================================================================
4172 void ApproxParameters(const Handle(GeomAdaptor_HSurface)& aHS1,
4173                       const Handle(GeomAdaptor_HSurface)& aHS2,
4174                       Standard_Integer& iDegMin,
4175                       Standard_Integer& iDegMax,
4176                       Standard_Integer& iNbIterMax,
4177                       Standard_Boolean& anWithPC
4178                       )
4179 {
4180   GeomAbs_SurfaceType aTS1, aTS2;
4181   //
4182   iDegMin=4;
4183   iDegMax=8;
4184   iNbIterMax=0;
4185   //
4186   aTS1=aHS1->Surface().GetType();
4187   aTS2=aHS2->Surface().GetType();
4188   //
4189   // Cylinder/Torus
4190   if ((aTS1==GeomAbs_Cylinder && aTS2==GeomAbs_Torus) ||
4191       (aTS2==GeomAbs_Cylinder && aTS1==GeomAbs_Torus)) {
4192     Standard_Real aRC, aRT, dR, aPC, aDPC, aScPr;
4193     gp_Cylinder aCylinder;
4194     gp_Torus aTorus;
4195     //
4196     aPC=Precision::Confusion();
4197     //
4198     aCylinder=(aTS1==GeomAbs_Cylinder)?
4199       aHS1->Surface().Cylinder() :
4200       aHS2->Surface().Cylinder();
4201     aTorus=(aTS1==GeomAbs_Torus)?
4202       aHS1->Surface().Torus() :
4203       aHS2->Surface().Torus();
4204     //
4205     aRC=aCylinder.Radius();
4206     aRT=aTorus.MinorRadius();
4207     dR=aRC-aRT;
4208     if (dR<0.) {
4209       dR=-dR;
4210     }
4211     //
4212     if (dR<aPC) {
4213       const gp_Ax3& aPosT=aTorus.Position();
4214       const gp_Dir& aDirT=aPosT.Direction();
4215       //
4216       const gp_Ax1& aAxisC=aCylinder.Axis();
4217       const gp_Dir& aDirC=aAxisC.Direction();
4218       //
4219       aScPr=aDirT*aDirC;
4220       if (aScPr<0.){
4221         aScPr=-aScPr;
4222       }
4223       //
4224       if (aScPr<aPC) {
4225         gp_Pln aPlnT(aPosT);
4226         const gp_Pnt& aLocC=aCylinder.Location();
4227
4228         aDPC=aPlnT.SquareDistance(aLocC);
4229         if (aDPC<2.5e-9) {
4230           iDegMax=6;
4231           anWithPC=Standard_False;
4232           iNbIterMax=10;
4233         }
4234       }
4235     }
4236   }
4237 }
4238 //modified by NIZNHY-PKV Thu Sep 23 08:37:39 2010t