src/IntTools/IntTools_FaceFace.cxx

   1 // Created on: 2000-11-23
   2 // Created by: Michael KLOKOV
   3 // Copyright (c) 2000-2014 OPEN CASCADE SAS
   4 //
   5 // This file is part of Open CASCADE Technology software library.
   6 //
   7 // This library is free software; you can redistribute it and/or modify it under
   8 // the terms of the GNU Lesser General Public License version 2.1 as published
   9 // by the Free Software Foundation, with special exception defined in the file
  10 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
  11 // distribution for complete text of the license and disclaimer of any warranty.
  12 //
  13 // Alternatively, this file may be used under the terms of Open CASCADE
  14 // commercial license or contractual agreement.
  15
  16 #include <IntTools_FaceFace.ixx>
  17
  18 #include <Precision.hxx>
  19
  20 #include <TColStd_HArray1OfReal.hxx>
  21 #include <TColStd_Array1OfReal.hxx>
  22 #include <TColStd_Array1OfInteger.hxx>
  23 #include <TColStd_SequenceOfReal.hxx>
  24 #include <TColStd_ListOfInteger.hxx>
  25 #include <TColStd_ListIteratorOfListOfInteger.hxx>
  26 #include <TColStd_Array1OfListOfInteger.hxx>
  27
  28 #include <gp_Lin2d.hxx>
  29 #include <gp_Ax22d.hxx>
  30 #include <gp_Circ2d.hxx>
  31 #include <gp_Torus.hxx>
  32 #include <gp_Cylinder.hxx>
  33
  34 #include <Bnd_Box.hxx>
  35
  36 #include <TColgp_HArray1OfPnt2d.hxx>
  37 #include <TColgp_SequenceOfPnt2d.hxx>
  38 #include <TColgp_Array1OfPnt.hxx>
  39 #include <TColgp_Array1OfPnt2d.hxx>
  40
  41 #include <IntAna_QuadQuadGeo.hxx>
  42
  43 #include <IntSurf_PntOn2S.hxx>
  44 #include <IntSurf_LineOn2S.hxx>
  45 #include <IntSurf_PntOn2S.hxx>
  46 #include <IntSurf_ListOfPntOn2S.hxx>
  47 #include <IntRes2d_Domain.hxx>
  48 #include <ProjLib_Plane.hxx>
  49
  50 #include <IntPatch_GLine.hxx>
  51 #include <IntPatch_RLine.hxx>
  52 #include <IntPatch_WLine.hxx>
  53 #include <IntPatch_ALine.hxx>
  54 #include <IntPatch_ALineToWLine.hxx>
  55
  56 #include <ElSLib.hxx>
  57 #include <ElCLib.hxx>
  58
  59 #include <Extrema_ExtCC.hxx>
  60 #include <Extrema_POnCurv.hxx>
  61 #include <BndLib_AddSurface.hxx>
  62
  63 #include <Adaptor3d_SurfacePtr.hxx>
  64 #include <Adaptor2d_HLine2d.hxx>
  65
  66 #include <GeomAbs_SurfaceType.hxx>
  67 #include <GeomAbs_CurveType.hxx>
  68
  69 #include <Geom_Surface.hxx>
  70 #include <Geom_Line.hxx>
  71 #include <Geom_Circle.hxx>
  72 #include <Geom_Ellipse.hxx>
  73 #include <Geom_Parabola.hxx>
  74 #include <Geom_Hyperbola.hxx>
  75 #include <Geom_TrimmedCurve.hxx>
  76 #include <Geom_BSplineCurve.hxx>
  77 #include <Geom_RectangularTrimmedSurface.hxx>
  78 #include <Geom_OffsetSurface.hxx>
  79 #include <Geom_Curve.hxx>
  80 #include <Geom_Conic.hxx>
  81
  82 #include <Geom2d_TrimmedCurve.hxx>
  83 #include <Geom2d_BSplineCurve.hxx>
  84 #include <Geom2d_Line.hxx>
  85 #include <Geom2d_Curve.hxx>
  86 #include <Geom2d_Circle.hxx>
  87
  88 #include <Geom2dAPI_InterCurveCurve.hxx>
  89 #include <Geom2dInt_GInter.hxx>
  90 #include <Geom2dAdaptor.hxx>
  91 #include <GeomAdaptor_Curve.hxx>
  92 #include <GeomAdaptor_HSurface.hxx>
  93 #include <GeomAdaptor_Surface.hxx>
  94 #include <GeomLib_CheckBSplineCurve.hxx>
  95 #include <GeomLib_Check2dBSplineCurve.hxx>
  96
  97 #include <GeomInt_WLApprox.hxx>
  98 #include <GeomProjLib.hxx>
  99 #include <GeomAPI_ProjectPointOnSurf.hxx>
 100 #include <Geom2dAdaptor_Curve.hxx>
 101 #include <TopoDS.hxx>
 102 #include <TopoDS_Edge.hxx>
 103 #include <TopExp_Explorer.hxx>
 104
 105 #include <BRep_Tool.hxx>
 106 #include <BRepTools.hxx>
 107 #include <BRepAdaptor_Surface.hxx>
 108
 109 #include <IntTools_Curve.hxx>
 110 #include <IntTools_Tools.hxx>
 111 #include <IntTools_Tools.hxx>
 112 #include <IntTools_TopolTool.hxx>
 113 #include <IntTools_PntOnFace.hxx>
 114 #include <IntTools_PntOn2Faces.hxx>
 115 #include <IntTools_Context.hxx>
 116 #include <IntSurf_ListIteratorOfListOfPntOn2S.hxx>
 117 #include <GeomInt.hxx>
 118
 119 #include <Approx_CurveOnSurface.hxx>
 120 #include <GeomAdaptor.hxx>
 121 #include <GeomInt_IntSS.hxx>
 122
 123 static
 124   void RefineVector(gp_Vec2d& aV2D);
 125 #ifdef OCCT_DEBUG_DUMPWLINE
 126 static
 127   void DumpWLine(const Handle(IntPatch_WLine)& aWLine);
 128 #endif
 129 //
 130 static
 131   void TolR3d(const TopoDS_Face& ,
 132               const TopoDS_Face& ,
 133               Standard_Real& );
 134 static
 135   Handle(Geom_Curve) MakeBSpline (const Handle(IntPatch_WLine)&,
 136                                   const Standard_Integer,
 137                                   const Standard_Integer);
 138
 139 static
 140   void Parameters(const Handle(GeomAdaptor_HSurface)&,
 141                   const Handle(GeomAdaptor_HSurface)&,
 142                   const gp_Pnt&,
 143                   Standard_Real&,
 144                   Standard_Real&,
 145                   Standard_Real&,
 146                   Standard_Real&);
 147
 148 static
 149   void CorrectSurfaceBoundaries(const TopoDS_Face&  theFace,
 150                                 const Standard_Real theTolerance,
 151                                 Standard_Real&      theumin,
 152                                 Standard_Real&      theumax,
 153                                 Standard_Real&      thevmin,
 154                                 Standard_Real&      thevmax);
 155
 156 static
 157   Standard_Boolean NotUseSurfacesForApprox
 158           (const TopoDS_Face& aF1,
 159            const TopoDS_Face& aF2,
 160            const Handle(IntPatch_WLine)& WL,
 161            const Standard_Integer ifprm,
 162            const Standard_Integer ilprm);
 163
 164 static
 165   Handle(IntPatch_WLine) ComputePurgedWLine(const Handle(IntPatch_WLine)& theWLine);
 166
 167 static
 168   Handle(Geom2d_BSplineCurve) MakeBSpline2d(const Handle(IntPatch_WLine)& theWLine,
 169                                             const Standard_Integer                         ideb,
 170                                             const Standard_Integer                         ifin,
 171                                             const Standard_Boolean                         onFirst);
 172
 173 static
 174   Standard_Boolean DecompositionOfWLine(const Handle(IntPatch_WLine)& theWLine,
 175                                         const Handle(GeomAdaptor_HSurface)&            theSurface1,
 176                                         const Handle(GeomAdaptor_HSurface)&            theSurface2,
 177                                         const TopoDS_Face&                             theFace1,
 178                                         const TopoDS_Face&                             theFace2,
 179                                         const GeomInt_LineConstructor&                 theLConstructor,
 180                                         const Standard_Boolean                         theAvoidLConstructor,
 181                                         IntPatch_SequenceOfLine&                       theNewLines,
 182                                         Standard_Real&                                 theReachedTol3d,
 183                                         const Handle(IntTools_Context)& );
 184
 185 static
 186   Standard_Boolean ParameterOutOfBoundary(const Standard_Real       theParameter,
 187                                           const Handle(Geom_Curve)& theCurve,
 188                                           const TopoDS_Face&        theFace1,
 189                                           const TopoDS_Face&        theFace2,
 190                                           const Standard_Real       theOtherParameter,
 191                                           const Standard_Boolean    bIncreasePar,
 192                                           Standard_Real&            theNewParameter,
 193                                           const Handle(IntTools_Context)& );
 194
 195 static
 196   Standard_Boolean IsCurveValid(Handle(Geom2d_Curve)& thePCurve);
 197
 198 static
 199   Standard_Boolean IsPointOnBoundary(const Standard_Real theParameter,
 200                                      const Standard_Real theFirstBoundary,
 201                                      const Standard_Real theSecondBoundary,
 202                                      const Standard_Real theResolution,
 203                                      Standard_Boolean&   IsOnFirstBoundary);
 204 static
 205   Standard_Boolean FindPoint(const gp_Pnt2d&     theFirstPoint,
 206                              const gp_Pnt2d&     theLastPoint,
 207                              const Standard_Real theUmin,
 208                              const Standard_Real theUmax,
 209                              const Standard_Real theVmin,
 210                              const Standard_Real theVmax,
 211                              gp_Pnt2d&           theNewPoint);
 212
 213
 214 static
 215   Standard_Integer ComputeTangentZones( const Handle(GeomAdaptor_HSurface)& theSurface1,
 216                                        const Handle(GeomAdaptor_HSurface)&  theSurface2,
 217                                        const TopoDS_Face&                   theFace1,
 218                                        const TopoDS_Face&                   theFace2,
 219                                        Handle(TColgp_HArray1OfPnt2d)&       theResultOnS1,
 220                                        Handle(TColgp_HArray1OfPnt2d)&       theResultOnS2,
 221                                        Handle(TColStd_HArray1OfReal)&       theResultRadius,
 222                                        const Handle(IntTools_Context)& );
 223
 224 static
 225   Standard_Boolean FindPoint(const gp_Pnt2d&     theFirstPoint,
 226                              const gp_Pnt2d&     theLastPoint,
 227                              const Standard_Real theUmin,
 228                              const Standard_Real theUmax,
 229                              const Standard_Real theVmin,
 230                              const Standard_Real theVmax,
 231                              const gp_Pnt2d&     theTanZoneCenter,
 232                              const Standard_Real theZoneRadius,
 233                              Handle(GeomAdaptor_HSurface) theGASurface,
 234                              gp_Pnt2d&           theNewPoint);
 235
 236 static
 237   Standard_Boolean IsInsideTanZone(const gp_Pnt2d&     thePoint,
 238                                    const gp_Pnt2d&     theTanZoneCenter,
 239                                    const Standard_Real theZoneRadius,
 240                                    Handle(GeomAdaptor_HSurface) theGASurface);
 241
 242 static
 243   gp_Pnt2d AdjustByNeighbour(const gp_Pnt2d&     theaNeighbourPoint,
 244                              const gp_Pnt2d&     theOriginalPoint,
 245                              Handle(GeomAdaptor_HSurface) theGASurface);
 246 static
 247   Standard_Boolean  ApproxWithPCurves(const gp_Cylinder& theCyl,
 248                                       const gp_Sphere& theSph);
 249
 250 static void  PerformPlanes(const Handle(GeomAdaptor_HSurface)& theS1,
 251                            const Handle(GeomAdaptor_HSurface)& theS2,
 252                            const Standard_Real TolAng,
 253                            const Standard_Real TolTang,
 254                            const Standard_Boolean theApprox1,
 255                            const Standard_Boolean theApprox2,
 256                            IntTools_SequenceOfCurves& theSeqOfCurve,
 257                            Standard_Boolean& theTangentFaces);
 258
 259 static Standard_Boolean ClassifyLin2d(const Handle(GeomAdaptor_HSurface)& theS,
 260                                       const gp_Lin2d& theLin2d,
 261                                       const Standard_Real theTol,
 262                                       Standard_Real& theP1,
 263                                       Standard_Real& theP2);
 264 //
 265 static
 266   void ApproxParameters(const Handle(GeomAdaptor_HSurface)& aHS1,
 267                         const Handle(GeomAdaptor_HSurface)& aHS2,
 268                         Standard_Integer& iDegMin,
 269                         Standard_Integer& iNbIter,
 270                         Standard_Integer& iDegMax);
 271
 272 static
 273   void Tolerances(const Handle(GeomAdaptor_HSurface)& aHS1,
 274                   const Handle(GeomAdaptor_HSurface)& aHS2,
 275                   Standard_Real& aTolTang);
 276
 277 static
 278   Standard_Boolean SortTypes(const GeomAbs_SurfaceType aType1,
 279                              const GeomAbs_SurfaceType aType2);
 280 static
 281   Standard_Integer IndexType(const GeomAbs_SurfaceType aType);
 282
 283 //
 284 static
 285   Standard_Boolean CheckPCurve(const Handle(Geom2d_Curve)& aPC,
 286                                const TopoDS_Face& aFace);
 287
 288 //
 289 static
 290   Standard_Real MaxDistance(const Handle(Geom_Curve)& theC,
 291                             const Standard_Real aT,
 292                             GeomAPI_ProjectPointOnSurf& theProjPS);
 293 static
 294   Standard_Real FindMaxDistance(const Handle(Geom_Curve)& theC,
 295                                 const Standard_Real theFirst,
 296                                 const Standard_Real theLast,
 297                                 GeomAPI_ProjectPointOnSurf& theProjPS,
 298                                 const Standard_Real theEps);
 299 static
 300   Standard_Real FindMaxDistance(const Handle(Geom_Curve)& theCurve,
 301                                 const Standard_Real theFirst,
 302                                 const Standard_Real theLast,
 303                                 const TopoDS_Face& theFace,
 304                                 const Handle(IntTools_Context)& theContext);
 305
 306 static
 307   void CorrectPlaneBoundaries(Standard_Real& aUmin,
 308                               Standard_Real& aUmax,
 309                               Standard_Real& aVmin,
 310                               Standard_Real& aVmax);
 311
 312 //=======================================================================
 313 //function :
 314 //purpose  :
 315 //=======================================================================
 316 IntTools_FaceFace::IntTools_FaceFace()
 317 {
 318   myIsDone=Standard_False;
 319   myTangentFaces=Standard_False;
 320   //
 321   myHS1 = new GeomAdaptor_HSurface ();
 322   myHS2 = new GeomAdaptor_HSurface ();
 323   myTolReached2d=0.;
 324   myTolReached3d=0.;
 325   SetParameters(Standard_True, Standard_True, Standard_True, 1.e-07);
 326
 327 }
 328 //=======================================================================
 329 //function : SetContext
 330 //purpose  :
 331 //=======================================================================
 332 void IntTools_FaceFace::SetContext(const Handle(IntTools_Context)& aContext)
 333 {
 334   myContext=aContext;
 335 }
 336 //=======================================================================
 337 //function : Context
 338 //purpose  :
 339 //=======================================================================
 340 const Handle(IntTools_Context)& IntTools_FaceFace::Context()const
 341 {
 342   return myContext;
 343 }
 344 //=======================================================================
 345 //function : Face1
 346 //purpose  :
 347 //=======================================================================
 348 const TopoDS_Face& IntTools_FaceFace::Face1() const
 349 {
 350   return myFace1;
 351 }
 352 //=======================================================================
 353 //function : Face2
 354 //purpose  :
 355 //=======================================================================
 356 const TopoDS_Face& IntTools_FaceFace::Face2() const
 357 {
 358   return myFace2;
 359 }
 360 //=======================================================================
 361 //function : TangentFaces
 362 //purpose  :
 363 //=======================================================================
 364 Standard_Boolean IntTools_FaceFace::TangentFaces() const
 365 {
 366   return myTangentFaces;
 367 }
 368 //=======================================================================
 369 //function : Points
 370 //purpose  :
 371 //=======================================================================
 372 const IntTools_SequenceOfPntOn2Faces& IntTools_FaceFace::Points() const
 373 {
 374   return myPnts;
 375 }
 376 //=======================================================================
 377 //function : IsDone
 378 //purpose  :
 379 //=======================================================================
 380 Standard_Boolean IntTools_FaceFace::IsDone() const
 381 {
 382   return myIsDone;
 383 }
 384 //=======================================================================
 385 //function : TolReached3d
 386 //purpose  :
 387 //=======================================================================
 388 Standard_Real IntTools_FaceFace::TolReached3d() const
 389 {
 390   return myTolReached3d;
 391 }
 392 //=======================================================================
 393 //function : Lines
 394 //purpose  : return lines of intersection
 395 //=======================================================================
 396 const IntTools_SequenceOfCurves& IntTools_FaceFace::Lines() const
 397 {
 398   StdFail_NotDone_Raise_if
 399     (!myIsDone,
 400      "IntTools_FaceFace::Lines() => myIntersector NOT DONE");
 401   return mySeqOfCurve;
 402 }
 403 //=======================================================================
 404 //function : TolReached2d
 405 //purpose  :
 406 //=======================================================================
 407 Standard_Real IntTools_FaceFace::TolReached2d() const
 408 {
 409   return myTolReached2d;
 410 }
 411 // =======================================================================
 412 // function: SetParameters
 413 //
 414 // =======================================================================
 415 void IntTools_FaceFace::SetParameters(const Standard_Boolean ToApproxC3d,
 416                                       const Standard_Boolean ToApproxC2dOnS1,
 417                                       const Standard_Boolean ToApproxC2dOnS2,
 418                                       const Standard_Real ApproximationTolerance)
 419 {
 420   myApprox = ToApproxC3d;
 421   myApprox1 = ToApproxC2dOnS1;
 422   myApprox2 = ToApproxC2dOnS2;
 423   myTolApprox = ApproximationTolerance;
 424 }
 425 //=======================================================================
 426 //function : SetList
 427 //purpose  :
 428 //=======================================================================
 429 void IntTools_FaceFace::SetList(IntSurf_ListOfPntOn2S& aListOfPnts)
 430 {
 431   myListOfPnts = aListOfPnts;
 432 }
 433
 434
 435 static Standard_Boolean isTreatAnalityc(const TopoDS_Face& theF1,
 436                                         const TopoDS_Face& theF2)
 437 {
 438   const Standard_Real Tolang = 1.e-8;
 439   const Standard_Real aTolF1=BRep_Tool::Tolerance(theF1);
 440   const Standard_Real aTolF2=BRep_Tool::Tolerance(theF2);
 441   const Standard_Real aTolSum = aTolF1 + aTolF2;
 442   Standard_Real aHigh = 0.0;
 443
 444   const BRepAdaptor_Surface aBAS1(theF1), aBAS2(theF2);
 445   const GeomAbs_SurfaceType aType1=aBAS1.GetType();
 446   const GeomAbs_SurfaceType aType2=aBAS2.GetType();
 447
 448   gp_Pln aS1;
 449   gp_Cylinder aS2;
 450   if(aType1 == GeomAbs_Plane)
 451   {
 452     aS1=aBAS1.Plane();
 453   }
 454   else if(aType2 == GeomAbs_Plane)
 455   {
 456     aS1=aBAS2.Plane();
 457   }
 458   else
 459   {
 460     return Standard_True;
 461   }
 462
 463   if(aType1 == GeomAbs_Cylinder)
 464   {
 465     aS2=aBAS1.Cylinder();
 466     const Standard_Real VMin = aBAS1.FirstVParameter();
 467     const Standard_Real VMax = aBAS1.LastVParameter();
 468
 469     if( Precision::IsNegativeInfinite(VMin) ||
 470         Precision::IsPositiveInfinite(VMax))
 471           return Standard_True;
 472     else
 473       aHigh = VMax - VMin;
 474   }
 475   else if(aType2 == GeomAbs_Cylinder)
 476   {
 477     aS2=aBAS2.Cylinder();
 478
 479     const Standard_Real VMin = aBAS2.FirstVParameter();
 480     const Standard_Real VMax = aBAS2.LastVParameter();
 481
 482     if( Precision::IsNegativeInfinite(VMin) ||
 483         Precision::IsPositiveInfinite(VMax))
 484           return Standard_True;
 485     else
 486       aHigh = VMax - VMin;
 487   }
 488   else
 489   {
 490     return Standard_True;
 491   }
 492
 493   IntAna_QuadQuadGeo inter;
 494   inter.Perform(aS1,aS2,Tolang,aTolSum, aHigh);
 495   if(inter.TypeInter() == IntAna_Ellipse)
 496   {
 497     const gp_Elips anEl = inter.Ellipse(1);
 498     const Standard_Real aMajorR = anEl.MajorRadius();
 499     const Standard_Real aMinorR = anEl.MinorRadius();
 500
 501     return (aMajorR < 100000.0 * aMinorR);
 502   }
 503   else
 504   {
 505     return inter.IsDone();
 506   }
 507 }
 508 //=======================================================================
 509 //function : Perform
 510 //purpose  : intersect surfaces of the faces
 511 //=======================================================================
 512 void IntTools_FaceFace::Perform(const TopoDS_Face& aF1,
 513                                 const TopoDS_Face& aF2)
 514 {
 515   Standard_Boolean RestrictLine = Standard_False, hasCone = Standard_False;
 516
 517   if (myContext.IsNull()) {
 518     myContext=new IntTools_Context;
 519   }
 520
 521   mySeqOfCurve.Clear();
 522   myTolReached2d=0.;
 523   myTolReached3d=0.;
 524   myIsDone = Standard_False;
 525   myNbrestr=0;//?
 526
 527   myFace1=aF1;
 528   myFace2=aF2;
 529
 530   const BRepAdaptor_Surface aBAS1(myFace1, Standard_False);
 531   const BRepAdaptor_Surface aBAS2(myFace2, Standard_False);
 532   GeomAbs_SurfaceType aType1=aBAS1.GetType();
 533   GeomAbs_SurfaceType aType2=aBAS2.GetType();
 534
 535   const Standard_Boolean bReverse=SortTypes(aType1, aType2);
 536   if (bReverse)
 537   {
 538     myFace1=aF2;
 539     myFace2=aF1;
 540     aType1=aBAS2.GetType();
 541     aType2=aBAS1.GetType();
 542
 543     if (myListOfPnts.Extent())
 544     {
 545       Standard_Real aU1,aV1,aU2,aV2;
 546       IntSurf_ListIteratorOfListOfPntOn2S aItP2S;
 547       //
 548       aItP2S.Initialize(myListOfPnts);
 549       for (; aItP2S.More(); aItP2S.Next())
 550       {
 551         IntSurf_PntOn2S& aP2S=aItP2S.Value();
 552         aP2S.Parameters(aU1,aV1,aU2,aV2);
 553         aP2S.SetValue(aU2,aV2,aU1,aV1);
 554       }
 555     }
 556     //
 557     Standard_Boolean anAproxTmp = myApprox1;
 558     myApprox1 = myApprox2;
 559     myApprox2 = anAproxTmp;
 560   }
 561
 562
 563   const Handle(Geom_Surface) S1=BRep_Tool::Surface(myFace1);
 564   const Handle(Geom_Surface) S2=BRep_Tool::Surface(myFace2);
 565
 566   const Standard_Real aTolF1=BRep_Tool::Tolerance(myFace1);
 567   const Standard_Real aTolF2=BRep_Tool::Tolerance(myFace2);
 568
 569   Standard_Real TolArc = aTolF1 + aTolF2;
 570   Standard_Real TolTang = TolArc;
 571
 572   const Standard_Boolean isFace1Quad = (aType1 == GeomAbs_Cylinder ||
 573                                         aType1 == GeomAbs_Cone ||
 574                                         aType1 == GeomAbs_Torus);
 575
 576   const Standard_Boolean isFace2Quad = (aType2 == GeomAbs_Cylinder ||
 577                                         aType2 == GeomAbs_Cone ||
 578                                         aType2 == GeomAbs_Torus);
 579
 580   if(aType1==GeomAbs_Plane && aType2==GeomAbs_Plane)  {
 581     Standard_Real umin, umax, vmin, vmax;
 582     //
 583     BRepTools::UVBounds(myFace1, umin, umax, vmin, vmax);
 584     CorrectPlaneBoundaries(umin, umax, vmin, vmax);
 585     myHS1->ChangeSurface().Load(S1, umin, umax, vmin, vmax);
 586     //
 587     BRepTools::UVBounds(myFace2, umin, umax, vmin, vmax);
 588     CorrectPlaneBoundaries(umin, umax, vmin, vmax);
 589     myHS2->ChangeSurface().Load(S2, umin, umax, vmin, vmax);
 590     //
 591     Standard_Real TolAng = 1.e-8;
 592     //
 593     PerformPlanes(myHS1, myHS2, TolAng, TolTang, myApprox1, myApprox2,
 594                   mySeqOfCurve, myTangentFaces);
 595     //
 596     myIsDone = Standard_True;
 597
 598     if(!myTangentFaces) {
 599       const Standard_Integer NbLinPP = mySeqOfCurve.Length();
 600       if(NbLinPP) {
 601         Standard_Real aTolFMax;
 602         myTolReached3d = 1.e-7;
 603         aTolFMax=Max(aTolF1, aTolF2);
 604         if (aTolFMax>myTolReached3d) {
 605           myTolReached3d=aTolFMax;
 606         }
 607         //
 608         myTolReached2d = myTolReached3d;
 609
 610         if (bReverse) {
 611           Handle(Geom2d_Curve) aC2D1, aC2D2;
 612           const Standard_Integer aNbLin = mySeqOfCurve.Length();
 613           for (Standard_Integer i = 1; i <= aNbLin; ++i) {
 614             IntTools_Curve& aIC=mySeqOfCurve(i);
 615             aC2D1=aIC.FirstCurve2d();
 616             aC2D2=aIC.SecondCurve2d();
 617             aIC.SetFirstCurve2d(aC2D2);
 618             aIC.SetSecondCurve2d(aC2D1);
 619           }
 620         }
 621       }
 622     }
 623     return;
 624   }//if(aType1==GeomAbs_Plane && aType2==GeomAbs_Plane){
 625
 626   if ((aType1==GeomAbs_Plane) && isFace2Quad)
 627   {
 628     Standard_Real umin, umax, vmin, vmax;
 629     // F1
 630     BRepTools::UVBounds(myFace1, umin, umax, vmin, vmax);
 631     CorrectPlaneBoundaries(umin, umax, vmin, vmax);
 632     myHS1->ChangeSurface().Load(S1, umin, umax, vmin, vmax);
 633     // F2
 634     BRepTools::UVBounds(myFace2, umin, umax, vmin, vmax);
 635     CorrectSurfaceBoundaries(myFace2, (aTolF1 + aTolF2) * 2., umin, umax, vmin, vmax);
 636     myHS2->ChangeSurface().Load(S2, umin, umax, vmin, vmax);
 637     //
 638     if( aType2==GeomAbs_Cone ) {
 639       TolArc = 0.0001;
 640       hasCone = Standard_True;
 641     }
 642   }
 643   else if ((aType2==GeomAbs_Plane) && isFace1Quad)
 644   {
 645     Standard_Real umin, umax, vmin, vmax;
 646     //F1
 647     BRepTools::UVBounds(myFace1, umin, umax, vmin, vmax);
 648     CorrectSurfaceBoundaries(myFace1, (aTolF1 + aTolF2) * 2., umin, umax, vmin, vmax);
 649     myHS1->ChangeSurface().Load(S1, umin, umax, vmin, vmax);
 650     // F2
 651     BRepTools::UVBounds(myFace2, umin, umax, vmin, vmax);
 652     CorrectPlaneBoundaries(umin, umax, vmin, vmax);
 653     myHS2->ChangeSurface().Load(S2, umin, umax, vmin, vmax);
 654     //
 655     if( aType1==GeomAbs_Cone ) {
 656       TolArc = 0.0001;
 657       hasCone = Standard_True;
 658     }
 659   }
 660   else
 661   {
 662     Standard_Real umin, umax, vmin, vmax;
 663     BRepTools::UVBounds(myFace1, umin, umax, vmin, vmax);
 664     CorrectSurfaceBoundaries(myFace1, (aTolF1 + aTolF2) * 2., umin, umax, vmin, vmax);
 665     myHS1->ChangeSurface().Load(S1, umin, umax, vmin, vmax);
 666     BRepTools::UVBounds(myFace2, umin, umax, vmin, vmax);
 667     CorrectSurfaceBoundaries(myFace2, (aTolF1 + aTolF2) * 2., umin, umax, vmin, vmax);
 668     myHS2->ChangeSurface().Load(S2, umin, umax, vmin, vmax);
 669   }
 670
 671   const Handle(IntTools_TopolTool) dom1 = new IntTools_TopolTool(myHS1);
 672   const Handle(IntTools_TopolTool) dom2 = new IntTools_TopolTool(myHS2);
 673
 674   myLConstruct.Load(dom1, dom2, myHS1, myHS2);
 675
 676
 677   Tolerances(myHS1, myHS2, TolTang);
 678
 679   {
 680     const Standard_Real UVMaxStep = 0.001;
 681     const Standard_Real Deflection = (hasCone) ? 0.085 : 0.1;
 682   myIntersector.SetTolerances(TolArc, TolTang, UVMaxStep, Deflection);
 683   }
 684
 685   if((myHS1->IsUClosed() && !myHS1->IsUPeriodic()) ||
 686      (myHS1->IsVClosed() && !myHS1->IsVPeriodic()) ||
 687      (myHS2->IsUClosed() && !myHS2->IsUPeriodic()) ||
 688      (myHS2->IsVClosed() && !myHS2->IsVPeriodic()))
 689   {
 690     RestrictLine = Standard_True;
 691   }
 692   //
 693   if((aType1 != GeomAbs_BSplineSurface) &&
 694       (aType1 != GeomAbs_BezierSurface)  &&
 695      (aType1 != GeomAbs_OtherSurface)  &&
 696      (aType2 != GeomAbs_BSplineSurface) &&
 697       (aType2 != GeomAbs_BezierSurface)  &&
 698      (aType2 != GeomAbs_OtherSurface))
 699   {
 700     RestrictLine = Standard_True;
 701
 702     if ((aType1 == GeomAbs_Torus) ||
 703         (aType2 == GeomAbs_Torus))
 704     {
 705       myListOfPnts.Clear();
 706     }
 707   }
 708
 709   //
 710   if(!RestrictLine)
 711   {
 712     TopExp_Explorer aExp;
 713     for(Standard_Integer i = 0; (!RestrictLine) && (i < 2); i++)
 714     {
 715       const TopoDS_Face& aF=(!i) ? myFace1 : myFace2;
 716       aExp.Init(aF, TopAbs_EDGE);
 717       for(; aExp.More(); aExp.Next())
 718       {
 719         const TopoDS_Edge& aE=TopoDS::Edge(aExp.Current());
 720
 721         if(BRep_Tool::Degenerated(aE))
 722         {
 723           RestrictLine = Standard_True;
 724           break;
 725         }
 726       }
 727     }
 728   }
 729
 730   const Standard_Boolean isGeomInt = isTreatAnalityc(aF1, aF2);
 731   myIntersector.Perform(myHS1, dom1, myHS2, dom2, TolArc, TolTang,
 732                                   myListOfPnts, RestrictLine, isGeomInt);
 733
 734   myIsDone = myIntersector.IsDone();
 735
 736   if (myIsDone)
 737   {
 738     myTangentFaces=myIntersector.TangentFaces();
 739     if (myTangentFaces) {
 740       return;
 741     }
 742     //
 743     if(RestrictLine) {
 744       myListOfPnts.Clear(); // to use LineConstructor
 745     }
 746     //
 747     const Standard_Integer aNbLin = myIntersector.NbLines();
 748     for (Standard_Integer i=1; i <= aNbLin; ++i) {
 749       MakeCurve(i, dom1, dom2);
 750     }
 751     //
 752     ComputeTolReached3d();
 753     //
 754     if (bReverse) {
 755       Handle(Geom2d_Curve) aC2D1, aC2D2;
 756       //
 757       const Standard_Integer aNbLin=mySeqOfCurve.Length();
 758       for (Standard_Integer i=1; i<=aNbLin; ++i)
 759       {
 760         IntTools_Curve& aIC=mySeqOfCurve(i);
 761         aC2D1=aIC.FirstCurve2d();
 762         aC2D2=aIC.SecondCurve2d();
 763         aIC.SetFirstCurve2d(aC2D2);
 764         aIC.SetSecondCurve2d(aC2D1);
 765       }
 766     }
 767
 768     // Points
 769     Standard_Real U1,V1,U2,V2;
 770     IntTools_PntOnFace aPntOnF1, aPntOnF2;
 771     IntTools_PntOn2Faces aPntOn2Faces;
 772     //
 773     const Standard_Integer aNbPnts = myIntersector.NbPnts();
 774     for (Standard_Integer i=1; i <= aNbPnts; ++i)
 775     {
 776       const IntSurf_PntOn2S& aISPnt=myIntersector.Point(i).PntOn2S();
 777       const gp_Pnt& aPnt=aISPnt.Value();
 778       aISPnt.Parameters(U1,V1,U2,V2);
 779       aPntOnF1.Init(myFace1, aPnt, U1, V1);
 780       aPntOnF2.Init(myFace2, aPnt, U2, V2);
 781       //
 782       if (!bReverse)
 783       {
 784         aPntOn2Faces.SetP1(aPntOnF1);
 785         aPntOn2Faces.SetP2(aPntOnF2);
 786       }
 787       else
 788       {
 789         aPntOn2Faces.SetP2(aPntOnF1);
 790         aPntOn2Faces.SetP1(aPntOnF2);
 791       }
 792
 793       myPnts.Append(aPntOn2Faces);
 794     }
 795   }
 796 }
 797
 798 //=======================================================================
 799 //function : ComputeTolerance
 800 //purpose  :
 801 //=======================================================================
 802 Standard_Real IntTools_FaceFace::ComputeTolerance()
 803 {
 804   Standard_Integer i, j, aNbLin;
 805   Standard_Real aFirst, aLast, aD, aDMax, aT;
 806   Handle(Geom_Surface) aS1, aS2;
 807   //
 808   aDMax = 0;
 809   aNbLin = mySeqOfCurve.Length();
 810   //
 811   aS1 = myHS1->ChangeSurface().Surface();
 812   aS2 = myHS2->ChangeSurface().Surface();
 813   //
 814   for (i = 1; i <= aNbLin; ++i) {
 815     const IntTools_Curve& aIC = mySeqOfCurve(i);
 816     const Handle(Geom_Curve)& aC3D = aIC.Curve();
 817     if (aC3D.IsNull()) {
 818       continue;
 819     }
 820     //
 821     aFirst = aC3D->FirstParameter();
 822     aLast  = aC3D->LastParameter();
 823     //
 824     const Handle(Geom2d_Curve)& aC2D1 = aIC.FirstCurve2d();
 825     const Handle(Geom2d_Curve)& aC2D2 = aIC.SecondCurve2d();
 826     //
 827     for (j = 0; j < 2; ++j) {
 828       const Handle(Geom2d_Curve)& aC2D = !j ? aC2D1 : aC2D2;
 829       const Handle(Geom_Surface)& aS = !j ? aS1 : aS2;
 830       //
 831       if (!aC2D.IsNull()) {
 832         if (IntTools_Tools::ComputeTolerance
 833             (aC3D, aC2D, aS, aFirst, aLast, aD, aT)) {
 834           if (aD > aDMax) {
 835             aDMax = aD;
 836           }
 837         }
 838       }
 839       //
 840       const TopoDS_Face& aF = !i ? myFace1 : myFace2;
 841       aD = FindMaxDistance(aC3D, aFirst, aLast, aF, myContext);
 842       if (aD > aDMax) {
 843         aDMax = aD;
 844       }
 845     }
 846   }
 847   //
 848   return aDMax;
 849 }
 850
 851 //=======================================================================
 852 //function :ComputeTolReached3d
 853 //purpose  :
 854 //=======================================================================
 855   void IntTools_FaceFace::ComputeTolReached3d()
 856 {
 857   Standard_Integer aNbLin, i;
 858   GeomAbs_SurfaceType aType1, aType2;
 859   //
 860   aNbLin=myIntersector.NbLines();
 861   if (!aNbLin) {
 862     return;
 863   }
 864   //
 865   aType1=myHS1->Surface().GetType();
 866   aType2=myHS2->Surface().GetType();
 867   //
 868   if (aType1==GeomAbs_Cylinder && aType2==GeomAbs_Cylinder) {
 869     if (aNbLin==2){
 870       Handle(IntPatch_Line) aIL1, aIL2;
 871       IntPatch_IType aTL1, aTL2;
 872       //
 873       aIL1=myIntersector.Line(1);
 874       aIL2=myIntersector.Line(2);
 875       aTL1=aIL1->ArcType();
 876       aTL2=aIL2->ArcType();
 877       if (aTL1==IntPatch_Lin && aTL2==IntPatch_Lin) {
 878         Standard_Real aD, aDTresh, dTol;
 879         gp_Lin aL1, aL2;
 880         //
 881         dTol=1.e-8;
 882         aDTresh=1.5e-6;
 883         //
 884         aL1=Handle(IntPatch_GLine)::DownCast(aIL1)->Line();
 885         aL2=Handle(IntPatch_GLine)::DownCast(aIL2)->Line();
 886         aD=aL1.Distance(aL2);
 887         aD=0.5*aD;
 888         if (aD<aDTresh) {
 889           myTolReached3d=aD+dTol;
 890         }
 891         return;
 892       }
 893     }
 894     //ZZ
 895     if (aNbLin) {// Check the distances
 896       Standard_Real aDMax;
 897       //
 898       aDMax = ComputeTolerance();
 899       if (aDMax > 0.) {
 900         myTolReached3d = aDMax;
 901       }
 902     }// if (aNbLin)
 903   }// if (aType1==GeomAbs_Cylinder && aType2==GeomAbs_Cylinder) {
 904   //
 905   //904/G3 f
 906   else if (aType1==GeomAbs_Plane && aType2==GeomAbs_Plane) {
 907     Standard_Real aTolF1, aTolF2, aTolFMax, aTolTresh;
 908     //
 909     aTolTresh=1.e-7;
 910     //
 911     aTolF1 = BRep_Tool::Tolerance(myFace1);
 912     aTolF2 = BRep_Tool::Tolerance(myFace2);
 913     aTolFMax=Max(aTolF1, aTolF2);
 914     //
 915     if (aTolFMax>aTolTresh) {
 916       myTolReached3d=aTolFMax;
 917     }
 918   }//if (aType1==GeomAbs_Plane && aType2==GeomAbs_Plane) {
 919   //t
 920   //IFV Bug OCC20297
 921   else if((aType1 == GeomAbs_Cylinder && aType2 == GeomAbs_Plane) ||
 922           (aType2 == GeomAbs_Cylinder && aType1 == GeomAbs_Plane)) {
 923     if(aNbLin == 1) {
 924       const Handle(IntPatch_Line)& aIL1 = myIntersector.Line(1);
 925       if(aIL1->ArcType() == IntPatch_Circle) {
 926         gp_Circ aCir = Handle(IntPatch_GLine)::DownCast(aIL1)->Circle();
 927         gp_XYZ aCirDir = aCir.Axis().Direction().XYZ();
 928         gp_XYZ aPlDir;
 929         gp_Pln aPln;
 930         if(aType1 == GeomAbs_Plane) {
 931           aPln = myHS1->Surface().Plane();
 932         }
 933         else {
 934           aPln = myHS2->Surface().Plane();
 935         }
 936         aPlDir = aPln.Axis().Direction().XYZ();
 937         Standard_Real cs = aCirDir*aPlDir;
 938         if(cs < 0.) aPlDir.Reverse();
 939         Standard_Real eps = 1.e-14;
 940         if(!aPlDir.IsEqual(aCirDir, eps)) {
 941           Standard_Integer aNbP = 11;
 942           Standard_Real dt = 2.*M_PI / (aNbP - 1), t;
 943           for(t = 0.; t < 2.*M_PI; t += dt) {
 944             Standard_Real d = aPln.Distance(ElCLib::Value(t, aCir));
 945             if(myTolReached3d < d) myTolReached3d = d;
 946           }
 947           myTolReached3d *= 1.1;
 948         }
 949       } //aIL1->ArcType() == IntPatch_Circle
 950     } //aNbLin == 1
 951   } // aType1 == GeomAbs_Cylinder && aType2 == GeomAbs_Plane)
 952   //End IFV Bug OCC20297
 953   //
 954   else if ((aType1==GeomAbs_Plane && aType2==GeomAbs_Torus) ||
 955            (aType2==GeomAbs_Plane && aType1==GeomAbs_Torus)) {
 956     aNbLin=mySeqOfCurve.Length();
 957     if (aNbLin!=1) {
 958       return;
 959     }
 960     //
 961     Standard_Integer aNbP;
 962     Standard_Real aT, aT1, aT2, dT, aUT, aVT, aUP, aVP;
 963     Standard_Real aDP, aDT, aDmax;
 964     gp_Pln aPln;
 965     gp_Torus aTorus;
 966     gp_Pnt aP, aPP, aPT;
 967     //
 968     const IntTools_Curve& aIC=mySeqOfCurve(1);
 969     const Handle(Geom_Curve)& aC3D=aIC.Curve();
 970     const Handle(Geom_BSplineCurve)& aBS=
 971       Handle(Geom_BSplineCurve)::DownCast(aC3D);
 972     if (aBS.IsNull()) {
 973       return;
 974     }
 975     //
 976     aT1=aBS->FirstParameter();
 977     aT2=aBS->LastParameter();
 978     //
 979     aPln  =(aType1==GeomAbs_Plane) ? myHS1->Plane() : myHS2->Plane();
 980     aTorus=(aType1==GeomAbs_Plane) ? myHS2->Torus() : myHS1->Torus();
 981     //
 982     aDmax=-1.;
 983     aNbP=11;
 984     dT=(aT2-aT1)/(aNbP-1);
 985     for (i=0; i<aNbP; ++i) {
 986       aT=aT1+i*dT;
 987       if (i==aNbP-1) {
 988         aT=aT2;
 989       }
 990       //
 991       aC3D->D0(aT, aP);
 992       //
 993       ElSLib::Parameters(aPln, aP, aUP, aVP);
 994       aPP=ElSLib::Value(aUP, aVP, aPln);
 995       aDP=aP.SquareDistance(aPP);
 996       if (aDP>aDmax) {
 997         aDmax=aDP;
 998       }
 999       //
1000       ElSLib::Parameters(aTorus, aP, aUT, aVT);
1001       aPT=ElSLib::Value(aUT, aVT, aTorus);
1002       aDT=aP.SquareDistance(aPT);
1003       if (aDT>aDmax) {
1004         aDmax=aDT;
1005       }
1006     }
1007     //
1008     if (aDmax > myTolReached3d*myTolReached3d) {
1009       myTolReached3d=sqrt(aDmax);
1010       myTolReached3d=1.1*myTolReached3d;
1011     }
1012   }// if ((aType1==GeomAbs_Plane && aType2==GeomAbs_Torus) ||
1013   //
1014   else if ((aType1==GeomAbs_SurfaceOfRevolution && aType2==GeomAbs_Cylinder) ||
1015            (aType2==GeomAbs_SurfaceOfRevolution && aType1==GeomAbs_Cylinder) ||
1016            (aType1==GeomAbs_Plane && aType2==GeomAbs_Sphere) ||
1017            (aType2==GeomAbs_Plane && aType1==GeomAbs_Sphere) ||
1018            (aType1==GeomAbs_Plane && aType2==GeomAbs_SurfaceOfExtrusion) ||
1019            (aType2==GeomAbs_Plane && aType1==GeomAbs_SurfaceOfExtrusion) ||
1020            (aType1==GeomAbs_Plane && aType2==GeomAbs_BSplineSurface) ||
1021            (aType2==GeomAbs_Plane && aType1==GeomAbs_BSplineSurface) ||
1022            !myApprox) {
1023     //
1024     Standard_Real aDMax;
1025     //
1026     aDMax = ComputeTolerance();
1027     if (aDMax > myTolReached3d) {
1028       myTolReached3d = aDMax;
1029     }
1030   }
1031 }
1032
1033 //=======================================================================
1034 //function : MakeCurve
1035 //purpose  :
1036 //=======================================================================
1037   void IntTools_FaceFace::MakeCurve(const Standard_Integer Index,
1038                                     const Handle(Adaptor3d_TopolTool)& dom1,
1039                                     const Handle(Adaptor3d_TopolTool)& dom2)
1040 {
1041   Standard_Boolean bDone, rejectSurface, reApprox, bAvoidLineConstructor;
1042   Standard_Boolean ok, bPCurvesOk;
1043   Standard_Integer i, j, aNbParts;
1044   Standard_Real fprm, lprm;
1045   Standard_Real Tolpc;
1046   Handle(IntPatch_Line) L;
1047   IntPatch_IType typl;
1048   Handle(Geom_Curve) newc;
1049   //
1050   const Standard_Real TOLCHECK   =0.0000001;
1051   const Standard_Real TOLANGCHECK=0.1;
1052   //
1053   rejectSurface = Standard_False;
1054   reApprox = Standard_False;
1055   //
1056   bPCurvesOk = Standard_True;
1057
1058  reapprox:;
1059
1060   Tolpc = myTolApprox;
1061   bAvoidLineConstructor = Standard_False;
1062   L = myIntersector.Line(Index);
1063   typl = L->ArcType();
1064   //
1065   if(typl==IntPatch_Walking) {
1066     Handle(IntPatch_Line) anewL;
1067     //
1068     const Handle(IntPatch_WLine)& aWLine=
1069       Handle(IntPatch_WLine)::DownCast(L);
1070     //DumpWLine(aWLine);
1071
1072     anewL = ComputePurgedWLine(aWLine);
1073     if(anewL.IsNull()) {
1074       return;
1075     }
1076     L = anewL;
1077
1078     //const Handle(IntPatch_WLine)& aWLineX = Handle(IntPatch_WLine)::DownCast(L);
1079     //DumpWLine(aWLineX);
1080
1081     //
1082     if(!myListOfPnts.IsEmpty()) {
1083       bAvoidLineConstructor = Standard_True;
1084     }
1085
1086     Standard_Integer nbp = aWLine->NbPnts();
1087     const IntSurf_PntOn2S& p1 = aWLine->Point(1);
1088     const IntSurf_PntOn2S& p2 = aWLine->Point(nbp);
1089
1090     const gp_Pnt& P1 = p1.Value();
1091     const gp_Pnt& P2 = p2.Value();
1092
1093     if(P1.SquareDistance(P2) < 1.e-14) {
1094       bAvoidLineConstructor = Standard_False;
1095     }
1096   }
1097
1098   typl=L->ArcType();
1099
1100   //
1101   // Line Constructor
1102   if(!bAvoidLineConstructor) {
1103     myLConstruct.Perform(L);
1104     //
1105     bDone=myLConstruct.IsDone();
1106     if(!bDone)
1107     {
1108       return;
1109     }
1110
1111     if(typl != IntPatch_Restriction)
1112     {
1113       aNbParts=myLConstruct.NbParts();
1114       if (aNbParts <= 0)
1115       {
1116         return;
1117       }
1118     }
1119   }
1120   // Do the Curve
1121
1122
1123   switch (typl) {
1124   //########################################
1125   // Line, Parabola, Hyperbola
1126   //########################################
1127   case IntPatch_Lin:
1128   case IntPatch_Parabola:
1129   case IntPatch_Hyperbola: {
1130     if (typl == IntPatch_Lin) {
1131       newc =
1132         new Geom_Line (Handle(IntPatch_GLine)::DownCast(L)->Line());
1133     }
1134
1135     else if (typl == IntPatch_Parabola) {
1136       newc =
1137         new Geom_Parabola(Handle(IntPatch_GLine)::DownCast(L)->Parabola());
1138     }
1139
1140     else if (typl == IntPatch_Hyperbola) {
1141       newc =
1142         new Geom_Hyperbola (Handle(IntPatch_GLine)::DownCast(L)->Hyperbola());
1143     }
1144     //
1145     // myTolReached3d
1146     if (typl == IntPatch_Lin) {
1147       TolR3d (myFace1, myFace2, myTolReached3d);
1148     }
1149     //
1150     aNbParts=myLConstruct.NbParts();
1151     for (i=1; i<=aNbParts; i++) {
1152       Standard_Boolean bFNIt, bLPIt;
1153       //
1154       myLConstruct.Part(i, fprm, lprm);
1155         //
1156       bFNIt=Precision::IsNegativeInfinite(fprm);
1157       bLPIt=Precision::IsPositiveInfinite(lprm);
1158       //
1159       if (!bFNIt && !bLPIt) {
1160         //
1161         IntTools_Curve aCurve;
1162         //
1163         Handle(Geom_TrimmedCurve) aCT3D=new Geom_TrimmedCurve(newc, fprm, lprm);
1164         aCurve.SetCurve(aCT3D);
1165         if (typl == IntPatch_Parabola) {
1166           Standard_Real aTolF1, aTolF2, aTolBase;
1167
1168           aTolF1 = BRep_Tool::Tolerance(myFace1);
1169           aTolF2 = BRep_Tool::Tolerance(myFace2);
1170           aTolBase=aTolF1+aTolF2;
1171           myTolReached3d=IntTools_Tools::CurveTolerance(aCT3D, aTolBase);
1172         }
1173         //
1174         aCurve.SetCurve(new Geom_TrimmedCurve(newc, fprm, lprm));
1175         if(myApprox1) {
1176           Handle (Geom2d_Curve) C2d;
1177           GeomInt_IntSS::BuildPCurves(fprm, lprm, Tolpc,
1178                 myHS1->ChangeSurface().Surface(), newc, C2d);
1179           if(Tolpc>myTolReached2d || myTolReached2d==0.) {
1180             myTolReached2d=Tolpc;
1181           }
1182             //
1183             aCurve.SetFirstCurve2d(new Geom2d_TrimmedCurve(C2d,fprm,lprm));
1184           }
1185           else {
1186             Handle(Geom2d_BSplineCurve) H1;
1187             //
1188             aCurve.SetFirstCurve2d(H1);
1189           }
1190         //
1191         if(myApprox2) {
1192           Handle (Geom2d_Curve) C2d;
1193           GeomInt_IntSS::BuildPCurves(fprm, lprm, Tolpc,
1194                     myHS2->ChangeSurface().Surface(), newc, C2d);
1195           if(Tolpc>myTolReached2d || myTolReached2d==0.) {
1196             myTolReached2d=Tolpc;
1197           }
1198           //
1199           aCurve.SetSecondCurve2d(new Geom2d_TrimmedCurve(C2d,fprm,lprm));
1200           }
1201         else {
1202           Handle(Geom2d_BSplineCurve) H1;
1203           //
1204           aCurve.SetSecondCurve2d(H1);
1205         }
1206         mySeqOfCurve.Append(aCurve);
1207       } //if (!bFNIt && !bLPIt) {
1208       else {
1209         //  on regarde si on garde
1210         //
1211         Standard_Real aTestPrm, dT=100.;
1212         //
1213         aTestPrm=0.;
1214         if (bFNIt && !bLPIt) {
1215           aTestPrm=lprm-dT;
1216         }
1217         else if (!bFNIt && bLPIt) {
1218           aTestPrm=fprm+dT;
1219         }
1220         else {
1221           // i.e, if (bFNIt && bLPIt)
1222           aTestPrm=IntTools_Tools::IntermediatePoint(-dT, dT);
1223         }
1224         //
1225         gp_Pnt ptref(newc->Value(aTestPrm));
1226         //
1227         GeomAbs_SurfaceType typS1 = myHS1->GetType();
1228         GeomAbs_SurfaceType typS2 = myHS2->GetType();
1229         if( typS1 == GeomAbs_SurfaceOfExtrusion ||
1230             typS1 == GeomAbs_OffsetSurface ||
1231             typS1 == GeomAbs_SurfaceOfRevolution ||
1232             typS2 == GeomAbs_SurfaceOfExtrusion ||
1233             typS2 == GeomAbs_OffsetSurface ||
1234             typS2 == GeomAbs_SurfaceOfRevolution) {
1235           Handle(Geom2d_BSplineCurve) H1;
1236           mySeqOfCurve.Append(IntTools_Curve(newc, H1, H1));
1237           continue;
1238         }
1239
1240         Standard_Real u1, v1, u2, v2, Tol;
1241
1242         Tol = Precision::Confusion();
1243         Parameters(myHS1, myHS2, ptref,  u1, v1, u2, v2);
1244         ok = (dom1->Classify(gp_Pnt2d(u1, v1), Tol) != TopAbs_OUT);
1245         if(ok) {
1246           ok = (dom2->Classify(gp_Pnt2d(u2,v2),Tol) != TopAbs_OUT);
1247         }
1248         if (ok) {
1249           Handle(Geom2d_BSplineCurve) H1;
1250           mySeqOfCurve.Append(IntTools_Curve(newc, H1, H1));
1251         }
1252       }
1253     }// for (i=1; i<=aNbParts; i++) {
1254   }// case IntPatch_Lin:  case IntPatch_Parabola:  case IntPatch_Hyperbola:
1255     break;
1256
1257   //########################################
1258   // Circle and Ellipse
1259   //########################################
1260   case IntPatch_Circle:
1261   case IntPatch_Ellipse: {
1262
1263     if (typl == IntPatch_Circle) {
1264       newc = new Geom_Circle
1265         (Handle(IntPatch_GLine)::DownCast(L)->Circle());
1266     }
1267     else { //IntPatch_Ellipse
1268       newc = new Geom_Ellipse
1269         (Handle(IntPatch_GLine)::DownCast(L)->Ellipse());
1270     }
1271     //
1272     // myTolReached3d
1273     TolR3d (myFace1, myFace2, myTolReached3d);
1274     //
1275     aNbParts=myLConstruct.NbParts();
1276     //
1277     Standard_Real aPeriod, aNul;
1278     TColStd_SequenceOfReal aSeqFprm,  aSeqLprm;
1279
1280     aNul=0.;
1281     aPeriod=M_PI+M_PI;
1282
1283     for (i=1; i<=aNbParts; i++) {
1284       myLConstruct.Part(i, fprm, lprm);
1285
1286       if (fprm < aNul && lprm > aNul) {
1287         // interval that goes through 0. is divided on two intervals;
1288         while (fprm<aNul || fprm>aPeriod) fprm=fprm+aPeriod;
1289         while (lprm<aNul || lprm>aPeriod) lprm=lprm+aPeriod;
1290         //
1291         if((aPeriod - fprm) > Tolpc) {
1292           aSeqFprm.Append(fprm);
1293           aSeqLprm.Append(aPeriod);
1294         }
1295         else {
1296           gp_Pnt P1 = newc->Value(fprm);
1297           gp_Pnt P2 = newc->Value(aPeriod);
1298           Standard_Real aTolDist = BRep_Tool::Tolerance(myFace1) + BRep_Tool::Tolerance(myFace2);
1299           aTolDist = (myTolReached3d > aTolDist) ? myTolReached3d : aTolDist;
1300
1301           if(P1.Distance(P2) > aTolDist) {
1302             Standard_Real anewpar = fprm;
1303
1304             if(ParameterOutOfBoundary(fprm, newc, myFace1, myFace2,
1305                                       lprm, Standard_False, anewpar, myContext)) {
1306               fprm = anewpar;
1307             }
1308             aSeqFprm.Append(fprm);
1309             aSeqLprm.Append(aPeriod);
1310           }
1311         }
1312
1313         //
1314         if((lprm - aNul) > Tolpc) {
1315           aSeqFprm.Append(aNul);
1316           aSeqLprm.Append(lprm);
1317         }
1318         else {
1319           gp_Pnt P1 = newc->Value(aNul);
1320           gp_Pnt P2 = newc->Value(lprm);
1321           Standard_Real aTolDist = BRep_Tool::Tolerance(myFace1) + BRep_Tool::Tolerance(myFace2);
1322           aTolDist = (myTolReached3d > aTolDist) ? myTolReached3d : aTolDist;
1323
1324           if(P1.Distance(P2) > aTolDist) {
1325             Standard_Real anewpar = lprm;
1326
1327             if(ParameterOutOfBoundary(lprm, newc, myFace1, myFace2,
1328                                       fprm, Standard_True, anewpar, myContext)) {
1329               lprm = anewpar;
1330             }
1331             aSeqFprm.Append(aNul);
1332             aSeqLprm.Append(lprm);
1333           }
1334         }
1335       }
1336       else {
1337         // usual interval
1338         aSeqFprm.Append(fprm);
1339         aSeqLprm.Append(lprm);
1340       }
1341     }
1342
1343     //
1344     aNbParts=aSeqFprm.Length();
1345     for (i=1; i<=aNbParts; i++) {
1346       fprm=aSeqFprm(i);
1347       lprm=aSeqLprm(i);
1348       //
1349       Standard_Real aRealEpsilon=RealEpsilon();
1350       if (Abs(fprm) > aRealEpsilon || Abs(lprm-2.*M_PI) > aRealEpsilon) {
1351         //==============================================
1352         ////
1353         IntTools_Curve aCurve;
1354         Handle(Geom_TrimmedCurve) aTC3D=new Geom_TrimmedCurve(newc,fprm,lprm);
1355         aCurve.SetCurve(aTC3D);
1356         fprm=aTC3D->FirstParameter();
1357         lprm=aTC3D->LastParameter ();
1358         ////
1359         if (typl == IntPatch_Circle || typl == IntPatch_Ellipse) {////
1360           if(myApprox1) {
1361             Handle (Geom2d_Curve) C2d;
1362             GeomInt_IntSS::BuildPCurves(fprm, lprm, Tolpc,
1363                         myHS1->ChangeSurface().Surface(), newc, C2d);
1364             if(Tolpc>myTolReached2d || myTolReached2d==0) {
1365               myTolReached2d=Tolpc;
1366             }
1367             //
1368             aCurve.SetFirstCurve2d(C2d);
1369           }
1370           else { ////
1371             Handle(Geom2d_BSplineCurve) H1;
1372             aCurve.SetFirstCurve2d(H1);
1373           }
1374
1375
1376           if(myApprox2) {
1377             Handle (Geom2d_Curve) C2d;
1378             GeomInt_IntSS::BuildPCurves(fprm,lprm,Tolpc,
1379                     myHS2->ChangeSurface().Surface(),newc,C2d);
1380             if(Tolpc>myTolReached2d || myTolReached2d==0) {
1381               myTolReached2d=Tolpc;
1382             }
1383             //
1384             aCurve.SetSecondCurve2d(C2d);
1385           }
1386           else {
1387             Handle(Geom2d_BSplineCurve) H1;
1388             aCurve.SetSecondCurve2d(H1);
1389           }
1390         }
1391
1392         else {
1393           Handle(Geom2d_BSplineCurve) H1;
1394           aCurve.SetFirstCurve2d(H1);
1395           aCurve.SetSecondCurve2d(H1);
1396         }
1397         mySeqOfCurve.Append(aCurve);
1398           //==============================================
1399       } //if (Abs(fprm) > RealEpsilon() || Abs(lprm-2.*M_PI) > RealEpsilon())
1400
1401       else {
1402         //  on regarde si on garde
1403         //
1404         if (aNbParts==1) {
1405 //           if (Abs(fprm) < RealEpsilon() &&  Abs(lprm-2.*M_PI) < RealEpsilon()) {
1406           if (Abs(fprm) <= aRealEpsilon && Abs(lprm-2.*M_PI) <= aRealEpsilon) {
1407             IntTools_Curve aCurve;
1408             Handle(Geom_TrimmedCurve) aTC3D=new Geom_TrimmedCurve(newc,fprm,lprm);
1409             aCurve.SetCurve(aTC3D);
1410             fprm=aTC3D->FirstParameter();
1411             lprm=aTC3D->LastParameter ();
1412
1413             if(myApprox1) {
1414               Handle (Geom2d_Curve) C2d;
1415               GeomInt_IntSS::BuildPCurves(fprm,lprm,Tolpc,
1416                     myHS1->ChangeSurface().Surface(),newc,C2d);
1417               if(Tolpc>myTolReached2d || myTolReached2d==0) {
1418                 myTolReached2d=Tolpc;
1419               }
1420               //
1421               aCurve.SetFirstCurve2d(C2d);
1422             }
1423             else { ////
1424               Handle(Geom2d_BSplineCurve) H1;
1425               aCurve.SetFirstCurve2d(H1);
1426             }
1427
1428             if(myApprox2) {
1429               Handle (Geom2d_Curve) C2d;
1430               GeomInt_IntSS::BuildPCurves(fprm,lprm,Tolpc,
1431                     myHS2->ChangeSurface().Surface(),newc,C2d);
1432               if(Tolpc>myTolReached2d || myTolReached2d==0) {
1433                 myTolReached2d=Tolpc;
1434               }
1435               //
1436               aCurve.SetSecondCurve2d(C2d);
1437             }
1438             else {
1439               Handle(Geom2d_BSplineCurve) H1;
1440               aCurve.SetSecondCurve2d(H1);
1441             }
1442             mySeqOfCurve.Append(aCurve);
1443             break;
1444           }
1445         }
1446         //
1447         Standard_Real aTwoPIdiv17, u1, v1, u2, v2, Tol;
1448
1449         aTwoPIdiv17=2.*M_PI/17.;
1450
1451         for (j=0; j<=17; j++) {
1452           gp_Pnt ptref (newc->Value (j*aTwoPIdiv17));
1453           Tol = Precision::Confusion();
1454
1455           Parameters(myHS1, myHS2, ptref, u1, v1, u2, v2);
1456           ok = (dom1->Classify(gp_Pnt2d(u1,v1),Tol) != TopAbs_OUT);
1457           if(ok) {
1458             ok = (dom2->Classify(gp_Pnt2d(u2,v2),Tol) != TopAbs_OUT);
1459           }
1460           if (ok) {
1461             IntTools_Curve aCurve;
1462             aCurve.SetCurve(newc);
1463             //==============================================
1464             if (typl == IntPatch_Circle || typl == IntPatch_Ellipse) {
1465
1466               if(myApprox1) {
1467                 Handle (Geom2d_Curve) C2d;
1468                 GeomInt_IntSS::BuildPCurves(fprm, lprm, Tolpc,
1469                         myHS1->ChangeSurface().Surface(), newc, C2d);
1470                 if(Tolpc>myTolReached2d || myTolReached2d==0) {
1471                   myTolReached2d=Tolpc;
1472                 }
1473                 //
1474                 aCurve.SetFirstCurve2d(C2d);
1475               }
1476               else {
1477                 Handle(Geom2d_BSplineCurve) H1;
1478                 aCurve.SetFirstCurve2d(H1);
1479               }
1480
1481               if(myApprox2) {
1482                 Handle (Geom2d_Curve) C2d;
1483                 GeomInt_IntSS::BuildPCurves(fprm, lprm, Tolpc,
1484                         myHS2->ChangeSurface().Surface(), newc, C2d);
1485                 if(Tolpc>myTolReached2d || myTolReached2d==0) {
1486                   myTolReached2d=Tolpc;
1487                 }
1488                 //
1489                 aCurve.SetSecondCurve2d(C2d);
1490               }
1491
1492               else {
1493                 Handle(Geom2d_BSplineCurve) H1;
1494                 aCurve.SetSecondCurve2d(H1);
1495               }
1496             }//  end of if (typl == IntPatch_Circle || typl == IntPatch_Ellipse)
1497
1498             else {
1499               Handle(Geom2d_BSplineCurve) H1;
1500               //
1501               aCurve.SetFirstCurve2d(H1);
1502               aCurve.SetSecondCurve2d(H1);
1503             }
1504             //==============================================
1505             //
1506             mySeqOfCurve.Append(aCurve);
1507             break;
1508
1509             }//  end of if (ok) {
1510           }//  end of for (Standard_Integer j=0; j<=17; j++)
1511         }//  end of else { on regarde si on garde
1512       }// for (i=1; i<=myLConstruct.NbParts(); i++)
1513     }// IntPatch_Circle: IntPatch_Ellipse:
1514     break;
1515
1516   case IntPatch_Analytic: {
1517     IntSurf_Quadric quad1,quad2;
1518     GeomAbs_SurfaceType typs = myHS1->Surface().GetType();
1519
1520     switch (typs) {
1521       case GeomAbs_Plane:
1522         quad1.SetValue(myHS1->Surface().Plane());
1523         break;
1524       case GeomAbs_Cylinder:
1525         quad1.SetValue(myHS1->Surface().Cylinder());
1526         break;
1527       case GeomAbs_Cone:
1528         quad1.SetValue(myHS1->Surface().Cone());
1529         break;
1530       case GeomAbs_Sphere:
1531         quad1.SetValue(myHS1->Surface().Sphere());
1532         break;
1533     case GeomAbs_Torus:
1534       quad1.SetValue(myHS1->Surface().Torus());
1535       break;
1536       default:
1537         Standard_ConstructionError::Raise("GeomInt_IntSS::MakeCurve 1");
1538       }
1539
1540     typs = myHS2->Surface().GetType();
1541
1542     switch (typs) {
1543       case GeomAbs_Plane:
1544         quad2.SetValue(myHS2->Surface().Plane());
1545         break;
1546       case GeomAbs_Cylinder:
1547         quad2.SetValue(myHS2->Surface().Cylinder());
1548         break;
1549       case GeomAbs_Cone:
1550         quad2.SetValue(myHS2->Surface().Cone());
1551         break;
1552       case GeomAbs_Sphere:
1553         quad2.SetValue(myHS2->Surface().Sphere());
1554         break;
1555     case GeomAbs_Torus:
1556       quad2.SetValue(myHS2->Surface().Torus());
1557       break;
1558       default:
1559         Standard_ConstructionError::Raise("GeomInt_IntSS::MakeCurve 2");
1560       }
1561     //
1562     //=========
1563     IntPatch_ALineToWLine convert (quad1, quad2);
1564
1565     if (!myApprox) {
1566       aNbParts=myLConstruct.NbParts();
1567       for (i=1; i<=aNbParts; i++) {
1568         myLConstruct.Part(i, fprm, lprm);
1569         Handle(IntPatch_WLine) WL =
1570           convert.MakeWLine(Handle(IntPatch_ALine)::DownCast(L), fprm, lprm);
1571         //
1572         Handle(Geom2d_BSplineCurve) H1;
1573         Handle(Geom2d_BSplineCurve) H2;
1574
1575         if(myApprox1) {
1576           H1 = MakeBSpline2d(WL, 1, WL->NbPnts(), Standard_True);
1577         }
1578
1579         if(myApprox2) {
1580           H2 = MakeBSpline2d(WL, 1, WL->NbPnts(), Standard_False);
1581         }
1582         //
1583         mySeqOfCurve.Append(IntTools_Curve(MakeBSpline(WL,1,WL->NbPnts()), H1, H2));
1584       }
1585     } // if (!myApprox)
1586
1587     else { // myApprox=TRUE
1588       GeomInt_WLApprox theapp3d;
1589       //
1590       Standard_Real tol2d = myTolApprox;
1591       //
1592       theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, Standard_True);
1593
1594       aNbParts=myLConstruct.NbParts();
1595       for (i=1; i<=aNbParts; i++) {
1596         myLConstruct.Part(i, fprm, lprm);
1597         Handle(IntPatch_WLine) WL =
1598           convert.MakeWLine(Handle(IntPatch_ALine):: DownCast(L),fprm,lprm);
1599
1600         theapp3d.Perform(myHS1,myHS2,WL,Standard_True,myApprox1,myApprox2, 1, WL->NbPnts());
1601
1602         if (!theapp3d.IsDone()) {
1603           //
1604           Handle(Geom2d_BSplineCurve) H1;
1605           Handle(Geom2d_BSplineCurve) H2;
1606
1607           if(myApprox1) {
1608             H1 = MakeBSpline2d(WL, 1, WL->NbPnts(), Standard_True);
1609           }
1610
1611           if(myApprox2) {
1612             H2 = MakeBSpline2d(WL, 1, WL->NbPnts(), Standard_False);
1613           }
1614           //
1615           mySeqOfCurve.Append(IntTools_Curve(MakeBSpline(WL,1,WL->NbPnts()), H1, H2));
1616         }
1617
1618         else {
1619           if(myApprox1 || myApprox2) {
1620             if( theapp3d.TolReached2d()>myTolReached2d || myTolReached2d==0) {
1621               myTolReached2d = theapp3d.TolReached2d();
1622             }
1623           }
1624
1625           if( theapp3d.TolReached3d()>myTolReached3d || myTolReached3d==0) {
1626             myTolReached3d = theapp3d.TolReached3d();
1627           }
1628
1629           Standard_Integer aNbMultiCurves, nbpoles;
1630           aNbMultiCurves=theapp3d.NbMultiCurves();
1631           for (j=1; j<=aNbMultiCurves; j++) {
1632             const AppParCurves_MultiBSpCurve& mbspc = theapp3d.Value(j);
1633             nbpoles = mbspc.NbPoles();
1634
1635             TColgp_Array1OfPnt tpoles(1, nbpoles);
1636             mbspc.Curve(1, tpoles);
1637             Handle(Geom_BSplineCurve) BS=new Geom_BSplineCurve(tpoles,
1638                                                                mbspc.Knots(),
1639                                                                mbspc.Multiplicities(),
1640                                                                mbspc.Degree());
1641
1642             GeomLib_CheckBSplineCurve Check(BS,TOLCHECK,TOLANGCHECK);
1643             Check.FixTangent(Standard_True,Standard_True);
1644             //
1645             IntTools_Curve aCurve;
1646             aCurve.SetCurve(BS);
1647
1648             if(myApprox1) {
1649               TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
1650               mbspc.Curve(2,tpoles2d);
1651               Handle(Geom2d_BSplineCurve) BS2=new Geom2d_BSplineCurve(tpoles2d,
1652                                                                       mbspc.Knots(),
1653                                                                       mbspc.Multiplicities(),
1654                                                                       mbspc.Degree());
1655
1656               GeomLib_Check2dBSplineCurve newCheck(BS2,TOLCHECK,TOLANGCHECK);
1657               newCheck.FixTangent(Standard_True,Standard_True);
1658               //
1659               aCurve.SetFirstCurve2d(BS2);
1660             }
1661             else {
1662               Handle(Geom2d_BSplineCurve) H1;
1663               aCurve.SetFirstCurve2d(H1);
1664             }
1665
1666             if(myApprox2) {
1667               TColgp_Array1OfPnt2d tpoles2d(1, nbpoles);
1668               Standard_Integer TwoOrThree;
1669               TwoOrThree=myApprox1 ? 3 : 2;
1670               mbspc.Curve(TwoOrThree, tpoles2d);
1671               Handle(Geom2d_BSplineCurve) BS2 =new Geom2d_BSplineCurve(tpoles2d,
1672                                                                        mbspc.Knots(),
1673                                                                        mbspc.Multiplicities(),
1674                                                                        mbspc.Degree());
1675
1676               GeomLib_Check2dBSplineCurve newCheck(BS2,TOLCHECK,TOLANGCHECK);
1677               newCheck.FixTangent(Standard_True,Standard_True);
1678               //
1679               aCurve.SetSecondCurve2d(BS2);
1680             }
1681             else {
1682               Handle(Geom2d_BSplineCurve) H2;
1683               aCurve.SetSecondCurve2d(H2);
1684             }
1685             //
1686             mySeqOfCurve.Append(aCurve);
1687
1688           }// for (j=1; j<=aNbMultiCurves; j++) {
1689         }// else from if (!theapp3d.IsDone())
1690       }// for (i=1; i<=aNbParts; i++) {
1691     }// else { // myApprox=TRUE
1692   }// case IntPatch_Analytic:
1693     break;
1694
1695   case IntPatch_Walking:{
1696     Handle(IntPatch_WLine) WL =
1697       Handle(IntPatch_WLine)::DownCast(L);
1698     //
1699     Standard_Integer ifprm, ilprm;
1700     //
1701     if (!myApprox) {
1702       aNbParts = 1;
1703       if(!bAvoidLineConstructor){
1704         aNbParts=myLConstruct.NbParts();
1705       }
1706       for (i=1; i<=aNbParts; ++i) {
1707         Handle(Geom2d_BSplineCurve) H1, H2;
1708         Handle(Geom_Curve) aBSp;
1709         //
1710         if(bAvoidLineConstructor) {
1711           ifprm = 1;
1712           ilprm = WL->NbPnts();
1713         }
1714         else {
1715           myLConstruct.Part(i, fprm, lprm);
1716           ifprm=(Standard_Integer)fprm;
1717           ilprm=(Standard_Integer)lprm;
1718         }
1719         //
1720         if(myApprox1) {
1721           H1 = MakeBSpline2d(WL, ifprm, ilprm, Standard_True);
1722         }
1723         //
1724         if(myApprox2) {
1725           H2 = MakeBSpline2d(WL, ifprm, ilprm, Standard_False);
1726         }
1727         //
1728         aBSp=MakeBSpline(WL, ifprm, ilprm);
1729         IntTools_Curve aIC(aBSp, H1, H2);
1730         mySeqOfCurve.Append(aIC);
1731       }// for (i=1; i<=aNbParts; ++i) {
1732     }// if (!myApprox) {
1733     //
1734     else { // X
1735       Standard_Boolean bIsDecomposited;
1736       Standard_Integer nbiter, aNbSeqOfL;
1737       Standard_Real tol2d, aTolApproxImp;
1738       IntPatch_SequenceOfLine aSeqOfL;
1739       GeomInt_WLApprox theapp3d;
1740       Approx_ParametrizationType aParType = Approx_ChordLength;
1741       //
1742       Standard_Boolean anApprox1 = myApprox1;
1743       Standard_Boolean anApprox2 = myApprox2;
1744       //
1745       aTolApproxImp=1.e-5;
1746       tol2d = myTolApprox;
1747
1748       GeomAbs_SurfaceType typs1, typs2;
1749       typs1 = myHS1->Surface().GetType();
1750       typs2 = myHS2->Surface().GetType();
1751       Standard_Boolean anWithPC = Standard_True;
1752
1753       if(typs1 == GeomAbs_Cylinder && typs2 == GeomAbs_Sphere) {
1754         anWithPC =
1755           ApproxWithPCurves(myHS1->Surface().Cylinder(), myHS2->Surface().Sphere());
1756       }
1757       else if (typs1 == GeomAbs_Sphere && typs2 == GeomAbs_Cylinder) {
1758         anWithPC =
1759           ApproxWithPCurves(myHS2->Surface().Cylinder(), myHS1->Surface().Sphere());
1760       }
1761       //
1762       if(!anWithPC) {
1763         myTolApprox = aTolApproxImp;//1.e-5;
1764         anApprox1 = Standard_False;
1765         anApprox2 = Standard_False;
1766         //
1767         tol2d = myTolApprox;
1768       }
1769
1770       if(myHS1 == myHS2) {
1771         theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, Standard_False, aParType);
1772         rejectSurface = Standard_True;
1773       }
1774       else {
1775         if(reApprox && !rejectSurface)
1776           theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, Standard_False, aParType);
1777         else {
1778           Standard_Integer iDegMax, iDegMin, iNbIter;
1779           //
1780           ApproxParameters(myHS1, myHS2, iDegMin, iDegMax, iNbIter);
1781           theapp3d.SetParameters(myTolApprox, tol2d, iDegMin, iDegMax, iNbIter, Standard_True, aParType);
1782         }
1783       }
1784       //
1785       Standard_Real aReachedTol = Precision::Confusion();
1786       bIsDecomposited=DecompositionOfWLine(WL,
1787                                            myHS1,
1788                                            myHS2,
1789                                            myFace1,
1790                                            myFace2,
1791                                            myLConstruct,
1792                                            bAvoidLineConstructor,
1793                                            aSeqOfL,
1794                                            aReachedTol,
1795                                            myContext);
1796       if ( bIsDecomposited && ( myTolReached3d < aReachedTol ) ) {
1797         myTolReached3d = aReachedTol;
1798       }
1799       //
1800       aNbSeqOfL=aSeqOfL.Length();
1801       //
1802       if (bIsDecomposited) {
1803         nbiter=aNbSeqOfL;
1804       }
1805       else {
1806         nbiter=1;
1807         aNbParts=1;
1808         if (!bAvoidLineConstructor) {
1809           aNbParts=myLConstruct.NbParts();
1810           nbiter=aNbParts;
1811         }
1812       }
1813       //
1814       for(i = 1; i <= nbiter; ++i) {
1815         if(bIsDecomposited) {
1816           WL = Handle(IntPatch_WLine)::DownCast(aSeqOfL.Value(i));
1817           ifprm = 1;
1818           ilprm = WL->NbPnts();
1819         }
1820         else {
1821           if(bAvoidLineConstructor) {
1822             ifprm = 1;
1823             ilprm = WL->NbPnts();
1824           }
1825           else {
1826             myLConstruct.Part(i, fprm, lprm);
1827             ifprm = (Standard_Integer)fprm;
1828             ilprm = (Standard_Integer)lprm;
1829           }
1830         }
1831         //-- lbr :
1832         //-- Si une des surfaces est un plan , on approxime en 2d
1833         //-- sur cette surface et on remonte les points 2d en 3d.
1834         if(typs1 == GeomAbs_Plane) {
1835           theapp3d.Perform(myHS1, myHS2, WL, Standard_False,Standard_True, myApprox2,ifprm,ilprm);
1836         }
1837         else if(typs2 == GeomAbs_Plane) {
1838           theapp3d.Perform(myHS1,myHS2,WL,Standard_False,myApprox1,Standard_True,ifprm,ilprm);
1839         }
1840         else {
1841           //
1842           if (myHS1 != myHS2){
1843             if ((typs1==GeomAbs_BezierSurface || typs1==GeomAbs_BSplineSurface) &&
1844                 (typs2==GeomAbs_BezierSurface || typs2==GeomAbs_BSplineSurface)) {
1845
1846               theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, Standard_True, aParType);
1847
1848               Standard_Boolean bUseSurfaces;
1849               bUseSurfaces=NotUseSurfacesForApprox(myFace1, myFace2, WL, ifprm,  ilprm);
1850               if (bUseSurfaces) {
1851                 // ######
1852                 rejectSurface = Standard_True;
1853                 // ######
1854                 theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, Standard_False, aParType);
1855               }
1856             }
1857           }
1858           //
1859           theapp3d.Perform(myHS1,myHS2,WL,Standard_True,anApprox1,anApprox2,ifprm,ilprm);
1860         }
1861           //
1862         if (!theapp3d.IsDone()) {
1863           Handle(Geom2d_BSplineCurve) H1;
1864           Handle(Geom2d_BSplineCurve) H2;
1865           //
1866           Handle(Geom_Curve) aBSp=MakeBSpline(WL,ifprm, ilprm);
1867           //
1868           if(myApprox1) {
1869             H1 = MakeBSpline2d(WL, ifprm, ilprm, Standard_True);
1870           }
1871           //
1872           if(myApprox2) {
1873             H2 = MakeBSpline2d(WL, ifprm, ilprm, Standard_False);
1874           }
1875           //
1876           IntTools_Curve aIC(aBSp, H1, H2);
1877           mySeqOfCurve.Append(aIC);
1878         }
1879
1880         else {
1881           if(myApprox1 || myApprox2 || (typs1==GeomAbs_Plane || typs2==GeomAbs_Plane)) {
1882             if( theapp3d.TolReached2d()>myTolReached2d || myTolReached2d==0.) {
1883               myTolReached2d = theapp3d.TolReached2d();
1884             }
1885           }
1886           if(typs1==GeomAbs_Plane || typs2==GeomAbs_Plane) {
1887             myTolReached3d = myTolReached2d;
1888             //
1889             if (typs1==GeomAbs_Torus || typs2==GeomAbs_Torus) {
1890               if (myTolReached3d<1.e-6) {
1891                 myTolReached3d = theapp3d.TolReached3d();
1892                 myTolReached3d=1.e-6;
1893               }
1894             }
1895           }
1896           else  if( theapp3d.TolReached3d()>myTolReached3d || myTolReached3d==0.) {
1897             myTolReached3d = theapp3d.TolReached3d();
1898           }
1899
1900           Standard_Integer aNbMultiCurves, nbpoles;
1901           aNbMultiCurves=theapp3d.NbMultiCurves();
1902           for (j=1; j<=aNbMultiCurves; j++) {
1903             if(typs1 == GeomAbs_Plane) {
1904               const AppParCurves_MultiBSpCurve& mbspc = theapp3d.Value(j);
1905               nbpoles = mbspc.NbPoles();
1906
1907               TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
1908               TColgp_Array1OfPnt   tpoles(1,nbpoles);
1909
1910               mbspc.Curve(1,tpoles2d);
1911               const gp_Pln&  Pln = myHS1->Surface().Plane();
1912               //
1913               Standard_Integer ik;
1914               for(ik = 1; ik<= nbpoles; ik++) {
1915                 tpoles.SetValue(ik,
1916                                 ElSLib::Value(tpoles2d.Value(ik).X(),
1917                                               tpoles2d.Value(ik).Y(),
1918                                               Pln));
1919               }
1920               //
1921               Handle(Geom_BSplineCurve) BS =
1922                 new Geom_BSplineCurve(tpoles,
1923                                       mbspc.Knots(),
1924                                       mbspc.Multiplicities(),
1925                                       mbspc.Degree());
1926               GeomLib_CheckBSplineCurve Check(BS,TOLCHECK,TOLANGCHECK);
1927               Check.FixTangent(Standard_True, Standard_True);
1928               //
1929               IntTools_Curve aCurve;
1930               aCurve.SetCurve(BS);
1931
1932               if(myApprox1) {
1933                 Handle(Geom2d_BSplineCurve) BS1 =
1934                   new Geom2d_BSplineCurve(tpoles2d,
1935                                           mbspc.Knots(),
1936                                           mbspc.Multiplicities(),
1937                                           mbspc.Degree());
1938                 GeomLib_Check2dBSplineCurve Check1(BS1,TOLCHECK,TOLANGCHECK);
1939                 Check1.FixTangent(Standard_True,Standard_True);
1940                 //
1941                 // ############################################
1942                 if(!rejectSurface && !reApprox) {
1943                   Standard_Boolean isValid = IsCurveValid(BS1);
1944                   if(!isValid) {
1945                     reApprox = Standard_True;
1946                     goto reapprox;
1947                   }
1948                 }
1949                 // ############################################
1950                 aCurve.SetFirstCurve2d(BS1);
1951               }
1952               else {
1953                 Handle(Geom2d_BSplineCurve) H1;
1954                 aCurve.SetFirstCurve2d(H1);
1955               }
1956
1957               if(myApprox2) {
1958                 mbspc.Curve(2, tpoles2d);
1959
1960                 Handle(Geom2d_BSplineCurve) BS2 = new Geom2d_BSplineCurve(tpoles2d,
1961                                                                           mbspc.Knots(),
1962                                                                           mbspc.Multiplicities(),
1963                                                                           mbspc.Degree());
1964                 GeomLib_Check2dBSplineCurve newCheck(BS2,TOLCHECK,TOLANGCHECK);
1965                 newCheck.FixTangent(Standard_True,Standard_True);
1966
1967                 // ###########################################
1968                 if(!rejectSurface && !reApprox) {
1969                   Standard_Boolean isValid = IsCurveValid(BS2);
1970                   if(!isValid) {
1971                     reApprox = Standard_True;
1972                     goto reapprox;
1973                   }
1974                 }
1975                 // ###########################################
1976                 //
1977                 aCurve.SetSecondCurve2d(BS2);
1978               }
1979               else {
1980                 Handle(Geom2d_BSplineCurve) H2;
1981                 //
1982                   aCurve.SetSecondCurve2d(H2);
1983               }
1984               //
1985               mySeqOfCurve.Append(aCurve);
1986             }//if(typs1 == GeomAbs_Plane) {
1987
1988             else if(typs2 == GeomAbs_Plane) {
1989               const AppParCurves_MultiBSpCurve& mbspc = theapp3d.Value(j);
1990               nbpoles = mbspc.NbPoles();
1991
1992               TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
1993               TColgp_Array1OfPnt   tpoles(1,nbpoles);
1994               mbspc.Curve((myApprox1==Standard_True)? 2 : 1,tpoles2d);
1995               const gp_Pln&  Pln = myHS2->Surface().Plane();
1996               //
1997               Standard_Integer ik;
1998               for(ik = 1; ik<= nbpoles; ik++) {
1999                 tpoles.SetValue(ik,
2000                                 ElSLib::Value(tpoles2d.Value(ik).X(),
2001                                               tpoles2d.Value(ik).Y(),
2002                                               Pln));
2003
2004               }
2005               //
2006               Handle(Geom_BSplineCurve) BS=new Geom_BSplineCurve(tpoles,
2007                                                                  mbspc.Knots(),
2008                                                                  mbspc.Multiplicities(),
2009                                                                  mbspc.Degree());
2010               GeomLib_CheckBSplineCurve Check(BS,TOLCHECK,TOLANGCHECK);
2011               Check.FixTangent(Standard_True,Standard_True);
2012               //
2013               IntTools_Curve aCurve;
2014               aCurve.SetCurve(BS);
2015
2016               if(myApprox2) {
2017                 Handle(Geom2d_BSplineCurve) BS1=new Geom2d_BSplineCurve(tpoles2d,
2018                                                                         mbspc.Knots(),
2019                                                                         mbspc.Multiplicities(),
2020                                                                         mbspc.Degree());
2021                 GeomLib_Check2dBSplineCurve Check1(BS1,TOLCHECK,TOLANGCHECK);
2022                 Check1.FixTangent(Standard_True,Standard_True);
2023                 //
2024                 // ###########################################
2025                 if(!rejectSurface && !reApprox) {
2026                   Standard_Boolean isValid = IsCurveValid(BS1);
2027                   if(!isValid) {
2028                     reApprox = Standard_True;
2029                     goto reapprox;
2030                   }
2031                 }
2032                 // ###########################################
2033                 bPCurvesOk = CheckPCurve(BS1, myFace2);
2034                 aCurve.SetSecondCurve2d(BS1);
2035               }
2036               else {
2037                 Handle(Geom2d_BSplineCurve) H2;
2038                 aCurve.SetSecondCurve2d(H2);
2039               }
2040
2041               if(myApprox1) {
2042                 mbspc.Curve(1,tpoles2d);
2043                 Handle(Geom2d_BSplineCurve) BS2=new Geom2d_BSplineCurve(tpoles2d,
2044                                                                         mbspc.Knots(),
2045                                                                         mbspc.Multiplicities(),
2046                                                                         mbspc.Degree());
2047                 GeomLib_Check2dBSplineCurve Check2(BS2,TOLCHECK,TOLANGCHECK);
2048                 Check2.FixTangent(Standard_True,Standard_True);
2049                 //
2050                 // ###########################################
2051                 if(!rejectSurface && !reApprox) {
2052                   Standard_Boolean isValid = IsCurveValid(BS2);
2053                   if(!isValid) {
2054                     reApprox = Standard_True;
2055                     goto reapprox;
2056                   }
2057                 }
2058                 // ###########################################
2059                 bPCurvesOk = bPCurvesOk && CheckPCurve(BS2, myFace1);
2060                 aCurve.SetFirstCurve2d(BS2);
2061               }
2062               else {
2063                 Handle(Geom2d_BSplineCurve) H1;
2064                 //
2065                 aCurve.SetFirstCurve2d(H1);
2066               }
2067               //
2068               //if points of the pcurves are out of the faces bounds
2069               //create 3d and 2d curves without approximation
2070               if (!bPCurvesOk) {
2071                 Handle(Geom2d_BSplineCurve) H1, H2;
2072                 bPCurvesOk = Standard_True;
2073                 //
2074                 Handle(Geom_Curve) aBSp=MakeBSpline(WL,ifprm, ilprm);
2075
2076                 if(myApprox1) {
2077                   H1 = MakeBSpline2d(WL, ifprm, ilprm, Standard_True);
2078                   bPCurvesOk = CheckPCurve(H1, myFace1);
2079                 }
2080
2081                 if(myApprox2) {
2082                   H2 = MakeBSpline2d(WL, ifprm, ilprm, Standard_False);
2083                   bPCurvesOk = bPCurvesOk && CheckPCurve(H2, myFace2);
2084                 }
2085                 //
2086                 //if pcurves created without approximation are out of the
2087                 //faces bounds, use approximated 3d and 2d curves
2088                 if (bPCurvesOk) {
2089                   IntTools_Curve aIC(aBSp, H1, H2);
2090                   mySeqOfCurve.Append(aIC);
2091                 } else {
2092                   mySeqOfCurve.Append(aCurve);
2093                 }
2094               } else {
2095                 mySeqOfCurve.Append(aCurve);
2096               }
2097             }// else if(typs2 == GeomAbs_Plane)
2098             //
2099             else { //typs2 != GeomAbs_Plane && typs1 != GeomAbs_Plane
2100               Standard_Boolean bIsValid1, bIsValid2;
2101               Handle(Geom_BSplineCurve) BS;
2102               Handle(Geom2d_BSplineCurve) aH2D;
2103               IntTools_Curve aCurve;
2104               //
2105               bIsValid1=Standard_True;
2106               bIsValid2=Standard_True;
2107               //
2108               const AppParCurves_MultiBSpCurve& mbspc = theapp3d.Value(j);
2109               nbpoles = mbspc.NbPoles();
2110               TColgp_Array1OfPnt tpoles(1,nbpoles);
2111               mbspc.Curve(1,tpoles);
2112               BS=new Geom_BSplineCurve(tpoles,
2113                                                                  mbspc.Knots(),
2114                                                                  mbspc.Multiplicities(),
2115                                                                  mbspc.Degree());
2116               GeomLib_CheckBSplineCurve Check(BS,TOLCHECK,TOLANGCHECK);
2117               Check.FixTangent(Standard_True,Standard_True);
2118               //
2119               aCurve.SetCurve(BS);
2120               aCurve.SetFirstCurve2d(aH2D);
2121               aCurve.SetSecondCurve2d(aH2D);
2122               //
2123               if(myApprox1) {
2124                 if(anApprox1) {
2125                   Handle(Geom2d_BSplineCurve) BS1;
2126                   TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
2127                   mbspc.Curve(2,tpoles2d);
2128                   //
2129                   BS1=new Geom2d_BSplineCurve(tpoles2d,
2130                                                                         mbspc.Knots(),
2131                                                                         mbspc.Multiplicities(),
2132                                                                         mbspc.Degree());
2133                   GeomLib_Check2dBSplineCurve newCheck(BS1,TOLCHECK,TOLANGCHECK);
2134                   newCheck.FixTangent(Standard_True,Standard_True);
2135                   //
2136                   if (!reApprox) {
2137                     bIsValid1=CheckPCurve(BS1, myFace1);
2138                   }
2139                   //
2140                   aCurve.SetFirstCurve2d(BS1);
2141                 }
2142                 else {
2143                   Handle(Geom2d_BSplineCurve) BS1;
2144                   fprm = BS->FirstParameter();
2145                   lprm = BS->LastParameter();
2146
2147                   Handle(Geom2d_Curve) C2d;
2148                   Standard_Real aTol = myTolApprox;
2149                   GeomInt_IntSS::BuildPCurves(fprm, lprm, aTol,
2150                             myHS1->ChangeSurface().Surface(), BS, C2d);
2151                   BS1 = Handle(Geom2d_BSplineCurve)::DownCast(C2d);
2152                   aCurve.SetFirstCurve2d(BS1);
2153                 }
2154               } // if(myApprox1) {
2155                 //
2156               if(myApprox2) {
2157                 if(anApprox2) {
2158                   Handle(Geom2d_BSplineCurve) BS2;
2159                   TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
2160                   mbspc.Curve((myApprox1==Standard_True)? 3 : 2,tpoles2d);
2161                   BS2=new Geom2d_BSplineCurve(tpoles2d,
2162                                                                         mbspc.Knots(),
2163                                                                         mbspc.Multiplicities(),
2164                                                                         mbspc.Degree());
2165                   GeomLib_Check2dBSplineCurve newCheck(BS2,TOLCHECK,TOLANGCHECK);
2166                   newCheck.FixTangent(Standard_True,Standard_True);
2167                 //
2168                   if (!reApprox) {
2169                     bIsValid2=CheckPCurve(BS2, myFace2);
2170                   }
2171                   aCurve.SetSecondCurve2d(BS2);
2172                 }
2173                 else {
2174                   Handle(Geom2d_BSplineCurve) BS2;
2175                   fprm = BS->FirstParameter();
2176                   lprm = BS->LastParameter();
2177
2178                   Handle(Geom2d_Curve) C2d;
2179                   Standard_Real aTol = myTolApprox;
2180                   GeomInt_IntSS::BuildPCurves(fprm, lprm, aTol,
2181                             myHS2->ChangeSurface().Surface(), BS, C2d);
2182                   BS2 = Handle(Geom2d_BSplineCurve)::DownCast(C2d);
2183                   aCurve.SetSecondCurve2d(BS2);
2184                 }
2185               } //if(myApprox2) {
2186               if (!bIsValid1 || !bIsValid2) {
2187                 myTolApprox=aTolApproxImp;//1.e-5;
2188                 tol2d = myTolApprox;
2189                 reApprox = Standard_True;
2190                 goto reapprox;
2191               }
2192                 //
2193               mySeqOfCurve.Append(aCurve);
2194             }
2195           }
2196         }
2197       }
2198     }// else { // X
2199   }// case IntPatch_Walking:{
2200     break;
2201
2202   case IntPatch_Restriction:
2203     {
2204       GeomAbs_SurfaceType typS1 = myHS1->Surface().GetType();
2205       GeomAbs_SurfaceType typS2 = myHS2->Surface().GetType();
2206       Standard_Boolean isAnalS1 = Standard_False;
2207       switch (typS1)
2208       {
2209       case GeomAbs_Plane:
2210       case GeomAbs_Cylinder:
2211       case GeomAbs_Sphere:
2212       case GeomAbs_Cone:
2213       case GeomAbs_Torus: isAnalS1 = Standard_True; break;
2214       default: break;
2215       }
2216
2217       Standard_Integer isAnalS2 = Standard_False;
2218       switch (typS2)
2219       {
2220       case GeomAbs_Plane:
2221       case GeomAbs_Cylinder:
2222       case GeomAbs_Sphere:
2223       case GeomAbs_Cone:
2224       case GeomAbs_Torus: isAnalS2 = Standard_True; break;
2225       default: break;
2226       }
2227
2228       Handle(IntPatch_RLine) RL =
2229         Handle(IntPatch_RLine)::DownCast(L);
2230       Handle(Geom_Curve) aC3d;
2231       Handle(Geom2d_Curve) aC2d1, aC2d2;
2232       Standard_Real aTolReached;
2233       GeomInt_IntSS::TreatRLine(RL, myHS1, myHS2, aC3d,
2234                                   aC2d1, aC2d2, aTolReached);
2235
2236       if(aC3d.IsNull())
2237         break;
2238
2239       Bnd_Box2d aBox1, aBox2;
2240
2241       const Standard_Real aU1f = myHS1->FirstUParameter(),
2242                           aV1f = myHS1->FirstVParameter(),
2243                           aU1l = myHS1->LastUParameter(),
2244                           aV1l = myHS1->LastVParameter();
2245       const Standard_Real aU2f = myHS2->FirstUParameter(),
2246                           aV2f = myHS2->FirstVParameter(),
2247                           aU2l = myHS2->LastUParameter(),
2248                           aV2l = myHS2->LastVParameter();
2249
2250       aBox1.Add(gp_Pnt2d(aU1f, aV1f));
2251       aBox1.Add(gp_Pnt2d(aU1l, aV1l));
2252       aBox2.Add(gp_Pnt2d(aU2f, aV2f));
2253       aBox2.Add(gp_Pnt2d(aU2l, aV2l));
2254
2255       GeomInt_VectorOfReal anArrayOfParameters;
2256
2257       //We consider here that the intersection line is same-parameter-line
2258       anArrayOfParameters.Append(aC3d->FirstParameter());
2259       anArrayOfParameters.Append(aC3d->LastParameter());
2260
2261       GeomInt_IntSS::
2262         TrimILineOnSurfBoundaries(aC2d1, aC2d2, aBox1, aBox2, anArrayOfParameters);
2263
2264       const Standard_Integer aNbIntersSolutionsm1 = anArrayOfParameters.Length() - 1;
2265
2266       //Trim RLine found.
2267       for(Standard_Integer anInd = 0; anInd < aNbIntersSolutionsm1; anInd++)
2268       {
2269         const Standard_Real aParF = anArrayOfParameters(anInd),
2270                             aParL = anArrayOfParameters(anInd+1);
2271
2272         if((aParL - aParF) <= Precision::PConfusion())
2273           continue;
2274
2275         const Standard_Real aPar = 0.5*(aParF + aParL);
2276         gp_Pnt2d aPt;
2277
2278         Handle(Geom2d_Curve) aCurv2d1, aCurv2d2;
2279         if(!aC2d1.IsNull())
2280         {
2281           aC2d1->D0(aPar, aPt);
2282
2283           if(aBox1.IsOut(aPt))
2284             continue;
2285
2286           if(myApprox1)
2287             aCurv2d1 = new Geom2d_TrimmedCurve(aC2d1, aParF, aParL);
2288         }
2289
2290         if(!aC2d2.IsNull())
2291         {
2292           aC2d2->D0(aPar, aPt);
2293
2294           if(aBox2.IsOut(aPt))
2295             continue;
2296
2297           if(myApprox2)
2298             aCurv2d2 = new Geom2d_TrimmedCurve(aC2d2, aParF, aParL);
2299         }
2300
2301         Handle(Geom_Curve) aCurv3d = new Geom_TrimmedCurve(aC3d, aParF, aParL);
2302
2303         IntTools_Curve aIC(aCurv3d, aCurv2d1, aCurv2d2);
2304         mySeqOfCurve.Append(aIC);
2305       }
2306     }
2307     break;
2308   default:
2309     break;
2310
2311   }
2312 }
2313
2314 //=======================================================================
2315 //function : Parameters
2316 //purpose  :
2317 //=======================================================================
2318  void Parameters(const Handle(GeomAdaptor_HSurface)& HS1,
2319                  const Handle(GeomAdaptor_HSurface)& HS2,
2320                  const gp_Pnt& Ptref,
2321                  Standard_Real& U1,
2322                  Standard_Real& V1,
2323                  Standard_Real& U2,
2324                  Standard_Real& V2)
2325 {
2326
2327   IntSurf_Quadric quad1,quad2;
2328   GeomAbs_SurfaceType typs = HS1->Surface().GetType();
2329
2330   switch (typs) {
2331   case GeomAbs_Plane:
2332     quad1.SetValue(HS1->Surface().Plane());
2333     break;
2334   case GeomAbs_Cylinder:
2335     quad1.SetValue(HS1->Surface().Cylinder());
2336     break;
2337   case GeomAbs_Cone:
2338     quad1.SetValue(HS1->Surface().Cone());
2339     break;
2340   case GeomAbs_Sphere:
2341     quad1.SetValue(HS1->Surface().Sphere());
2342     break;
2343   case GeomAbs_Torus:
2344     quad1.SetValue(HS1->Surface().Torus());
2345     break;
2346   default:
2347     Standard_ConstructionError::Raise("GeomInt_IntSS::MakeCurve");
2348   }
2349
2350   typs = HS2->Surface().GetType();
2351   switch (typs) {
2352   case GeomAbs_Plane:
2353     quad2.SetValue(HS2->Surface().Plane());
2354     break;
2355   case GeomAbs_Cylinder:
2356     quad2.SetValue(HS2->Surface().Cylinder());
2357     break;
2358   case GeomAbs_Cone:
2359     quad2.SetValue(HS2->Surface().Cone());
2360     break;
2361   case GeomAbs_Sphere:
2362     quad2.SetValue(HS2->Surface().Sphere());
2363     break;
2364   case GeomAbs_Torus:
2365     quad2.SetValue(HS2->Surface().Torus());
2366     break;
2367   default:
2368     Standard_ConstructionError::Raise("GeomInt_IntSS::MakeCurve");
2369   }
2370
2371   quad1.Parameters(Ptref,U1,V1);
2372   quad2.Parameters(Ptref,U2,V2);
2373 }
2374
2375 //=======================================================================
2376 //function : MakeBSpline
2377 //purpose  :
2378 //=======================================================================
2379 Handle(Geom_Curve) MakeBSpline  (const Handle(IntPatch_WLine)& WL,
2380                                  const Standard_Integer ideb,
2381                                  const Standard_Integer ifin)
2382 {
2383   Standard_Integer i,nbpnt = ifin-ideb+1;
2384   TColgp_Array1OfPnt poles(1,nbpnt);
2385   TColStd_Array1OfReal knots(1,nbpnt);
2386   TColStd_Array1OfInteger mults(1,nbpnt);
2387   Standard_Integer ipidebm1;
2388   for(i=1,ipidebm1=i+ideb-1; i<=nbpnt;ipidebm1++, i++) {
2389     poles(i) = WL->Point(ipidebm1).Value();
2390     mults(i) = 1;
2391     knots(i) = i-1;
2392   }
2393   mults(1) = mults(nbpnt) = 2;
2394   return
2395     new Geom_BSplineCurve(poles,knots,mults,1);
2396 }
2397 //
2398
2399 //=======================================================================
2400 //function : MakeBSpline2d
2401 //purpose  :
2402 //=======================================================================
2403 Handle(Geom2d_BSplineCurve) MakeBSpline2d(const Handle(IntPatch_WLine)& theWLine,
2404                                           const Standard_Integer ideb,
2405                                           const Standard_Integer ifin,
2406                                           const Standard_Boolean onFirst)
2407 {
2408   Standard_Integer i, nbpnt = ifin-ideb+1;
2409   TColgp_Array1OfPnt2d poles(1,nbpnt);
2410   TColStd_Array1OfReal knots(1,nbpnt);
2411   TColStd_Array1OfInteger mults(1,nbpnt);
2412   Standard_Integer ipidebm1;
2413
2414   for(i = 1, ipidebm1 = i+ideb-1; i <= nbpnt; ipidebm1++, i++) {
2415       Standard_Real U, V;
2416       if(onFirst)
2417         theWLine->Point(ipidebm1).ParametersOnS1(U, V);
2418       else
2419         theWLine->Point(ipidebm1).ParametersOnS2(U, V);
2420       poles(i).SetCoord(U, V);
2421       mults(i) = 1;
2422       knots(i) = i-1;
2423     }
2424     mults(1) = mults(nbpnt) = 2;
2425
2426   return new Geom2d_BSplineCurve(poles,knots,mults,1);
2427 }
2428 //=======================================================================
2429 //function : PrepareLines3D
2430 //purpose  :
2431 //=======================================================================
2432   void IntTools_FaceFace::PrepareLines3D(const Standard_Boolean bToSplit)
2433 {
2434   Standard_Integer i, aNbCurves;
2435   GeomAbs_SurfaceType aType1, aType2;
2436   IntTools_SequenceOfCurves aNewCvs;
2437   //
2438   // 1. Treatment closed  curves
2439   aNbCurves=mySeqOfCurve.Length();
2440   for (i=1; i<=aNbCurves; ++i) {
2441     const IntTools_Curve& aIC=mySeqOfCurve(i);
2442     //
2443     if (bToSplit) {
2444       Standard_Integer j, aNbC;
2445       IntTools_SequenceOfCurves aSeqCvs;
2446       //
2447       aNbC=IntTools_Tools::SplitCurve(aIC, aSeqCvs);
2448       if (aNbC) {
2449         for (j=1; j<=aNbC; ++j) {
2450           const IntTools_Curve& aICNew=aSeqCvs(j);
2451           aNewCvs.Append(aICNew);
2452         }
2453       }
2454       else {
2455         aNewCvs.Append(aIC);
2456       }
2457     }
2458     else {
2459       aNewCvs.Append(aIC);
2460     }
2461   }
2462   //
2463   // 2. Plane\Cone intersection when we had 4 curves
2464   aType1=myHS1->GetType();
2465   aType2=myHS2->GetType();
2466   aNbCurves=aNewCvs.Length();
2467   //
2468   if ((aType1==GeomAbs_Plane && aType2==GeomAbs_Cone) ||
2469       (aType2==GeomAbs_Plane && aType1==GeomAbs_Cone)) {
2470     if (aNbCurves==4) {
2471       GeomAbs_CurveType aCType1;
2472       //
2473       aCType1=aNewCvs(1).Type();
2474       if (aCType1==GeomAbs_Line) {
2475         IntTools_SequenceOfCurves aSeqIn, aSeqOut;
2476         //
2477         for (i=1; i<=aNbCurves; ++i) {
2478           const IntTools_Curve& aIC=aNewCvs(i);
2479           aSeqIn.Append(aIC);
2480         }
2481         //
2482         IntTools_Tools::RejectLines(aSeqIn, aSeqOut);
2483         //
2484         aNewCvs.Clear();
2485         aNbCurves=aSeqOut.Length();
2486         for (i=1; i<=aNbCurves; ++i) {
2487           const IntTools_Curve& aIC=aSeqOut(i);
2488           aNewCvs.Append(aIC);
2489         }
2490       }
2491     }
2492   }// if ((aType1==GeomAbs_Plane && aType2==GeomAbs_Cone)...
2493   //
2494   // 3. Fill  mySeqOfCurve
2495   mySeqOfCurve.Clear();
2496   aNbCurves=aNewCvs.Length();
2497   for (i=1; i<=aNbCurves; ++i) {
2498     const IntTools_Curve& aIC=aNewCvs(i);
2499     mySeqOfCurve.Append(aIC);
2500   }
2501 }
2502 //=======================================================================
2503 //function : CorrectSurfaceBoundaries
2504 //purpose  :
2505 //=======================================================================
2506  void CorrectSurfaceBoundaries(const TopoDS_Face&  theFace,
2507                               const Standard_Real theTolerance,
2508                               Standard_Real&      theumin,
2509                               Standard_Real&      theumax,
2510                               Standard_Real&      thevmin,
2511                               Standard_Real&      thevmax)
2512 {
2513   Standard_Boolean enlarge, isuperiodic, isvperiodic;
2514   Standard_Real uinf, usup, vinf, vsup, delta;
2515   GeomAbs_SurfaceType aType;
2516   Handle(Geom_Surface) aSurface;
2517   //
2518   aSurface = BRep_Tool::Surface(theFace);
2519   aSurface->Bounds(uinf, usup, vinf, vsup);
2520   delta = theTolerance;
2521   enlarge = Standard_False;
2522   //
2523   GeomAdaptor_Surface anAdaptorSurface(aSurface);
2524   //
2525   if(aSurface->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface))) {
2526     Handle(Geom_Surface) aBasisSurface =
2527       (Handle(Geom_RectangularTrimmedSurface)::DownCast(aSurface))->BasisSurface();
2528
2529     if(aBasisSurface->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface)) ||
2530        aBasisSurface->IsKind(STANDARD_TYPE(Geom_OffsetSurface))) {
2531       return;
2532     }
2533   }
2534   //
2535   if(aSurface->IsKind(STANDARD_TYPE(Geom_OffsetSurface))) {
2536     Handle(Geom_Surface) aBasisSurface =
2537       (Handle(Geom_OffsetSurface)::DownCast(aSurface))->BasisSurface();
2538
2539     if(aBasisSurface->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface)) ||
2540        aBasisSurface->IsKind(STANDARD_TYPE(Geom_OffsetSurface))) {
2541       return;
2542     }
2543   }
2544   //
2545   isuperiodic = anAdaptorSurface.IsUPeriodic();
2546   isvperiodic = anAdaptorSurface.IsVPeriodic();
2547   //
2548   aType=anAdaptorSurface.GetType();
2549   if((aType==GeomAbs_BezierSurface) ||
2550      (aType==GeomAbs_BSplineSurface) ||
2551      (aType==GeomAbs_SurfaceOfExtrusion) ||
2552      (aType==GeomAbs_SurfaceOfRevolution) ||
2553      (aType==GeomAbs_Cylinder)) {
2554     enlarge=Standard_True;
2555   }
2556   //
2557   if(!isuperiodic && enlarge) {
2558
2559     if(!Precision::IsInfinite(theumin) &&
2560         ((theumin - uinf) > delta))
2561       theumin -= delta;
2562     else {
2563       theumin = uinf;
2564     }
2565
2566     if(!Precision::IsInfinite(theumax) &&
2567         ((usup - theumax) > delta))
2568       theumax += delta;
2569     else
2570       theumax = usup;
2571   }
2572   //
2573   if(!isvperiodic && enlarge) {
2574     if(!Precision::IsInfinite(thevmin) &&
2575         ((thevmin - vinf) > delta)) {
2576       thevmin -= delta;
2577     }
2578     else {
2579       thevmin = vinf;
2580     }
2581     if(!Precision::IsInfinite(thevmax) &&
2582         ((vsup - thevmax) > delta)) {
2583       thevmax += delta;
2584     }
2585     else {
2586       thevmax = vsup;
2587     }
2588   }
2589   //
2590   {
2591     Standard_Integer aNbP;
2592     Standard_Real aXP, dXfact, aXmid, aX1, aX2, aTolPA;
2593     //
2594     aTolPA=Precision::Angular();
2595     // U
2596     if (isuperiodic) {
2597       aXP=anAdaptorSurface.UPeriod();
2598       dXfact=theumax-theumin;
2599       if (dXfact-aTolPA>aXP) {
2600         aXmid=0.5*(theumax+theumin);
2601         aNbP=RealToInt(aXmid/aXP);
2602         if (aXmid<0.) {
2603           aNbP=aNbP-1;
2604         }
2605         aX1=aNbP*aXP;
2606         if (theumin>aTolPA) {
2607           aX1=theumin+aNbP*aXP;
2608         }
2609         aX2=aX1+aXP;
2610         if (theumin<aX1) {
2611           theumin=aX1;
2612         }
2613         if (theumax>aX2) {
2614           theumax=aX2;
2615         }
2616       }
2617     }
2618     // V
2619     if (isvperiodic) {
2620       aXP=anAdaptorSurface.VPeriod();
2621       dXfact=thevmax-thevmin;
2622       if (dXfact-aTolPA>aXP) {
2623         aXmid=0.5*(thevmax+thevmin);
2624         aNbP=RealToInt(aXmid/aXP);
2625         if (aXmid<0.) {
2626           aNbP=aNbP-1;
2627         }
2628         aX1=aNbP*aXP;
2629         if (thevmin>aTolPA) {
2630           aX1=thevmin+aNbP*aXP;
2631         }
2632         aX2=aX1+aXP;
2633         if (thevmin<aX1) {
2634           thevmin=aX1;
2635         }
2636         if (thevmax>aX2) {
2637           thevmax=aX2;
2638         }
2639       }
2640     }
2641   }
2642   //
2643   if(isuperiodic || isvperiodic) {
2644     Standard_Boolean correct = Standard_False;
2645     Standard_Boolean correctU = Standard_False;
2646     Standard_Boolean correctV = Standard_False;
2647     Bnd_Box2d aBox;
2648     TopExp_Explorer anExp;
2649
2650     for(anExp.Init(theFace, TopAbs_EDGE); anExp.More(); anExp.Next()) {
2651       if(BRep_Tool::IsClosed(TopoDS::Edge(anExp.Current()), theFace)) {
2652         correct = Standard_True;
2653         Standard_Real f, l;
2654         TopoDS_Edge anEdge = TopoDS::Edge(anExp.Current());
2655
2656         for(Standard_Integer i = 0; i < 2; i++) {
2657           if(i==0) {
2658             anEdge.Orientation(TopAbs_FORWARD);
2659           }
2660           else {
2661             anEdge.Orientation(TopAbs_REVERSED);
2662           }
2663           Handle(Geom2d_Curve) aCurve = BRep_Tool::CurveOnSurface(anEdge, theFace, f, l);
2664
2665           if(aCurve.IsNull()) {
2666             correct = Standard_False;
2667             break;
2668           }
2669           Handle(Geom2d_Line) aLine = Handle(Geom2d_Line)::DownCast(aCurve);
2670
2671           if(aLine.IsNull()) {
2672             correct = Standard_False;
2673             break;
2674           }
2675           gp_Dir2d anUDir(1., 0.);
2676           gp_Dir2d aVDir(0., 1.);
2677           Standard_Real anAngularTolerance = Precision::Angular();
2678
2679           correctU = correctU || aLine->Position().Direction().IsParallel(aVDir, anAngularTolerance);
2680           correctV = correctV || aLine->Position().Direction().IsParallel(anUDir, anAngularTolerance);
2681
2682           gp_Pnt2d pp1 = aCurve->Value(f);
2683           aBox.Add(pp1);
2684           gp_Pnt2d pp2 = aCurve->Value(l);
2685           aBox.Add(pp2);
2686         }
2687         if(!correct)
2688           break;
2689       }
2690     }
2691
2692     if(correct) {
2693       Standard_Real umin, vmin, umax, vmax;
2694       aBox.Get(umin, vmin, umax, vmax);
2695
2696       if(isuperiodic && correctU) {
2697
2698         if(theumin < umin)
2699           theumin = umin;
2700
2701         if(theumax > umax) {
2702           theumax = umax;
2703         }
2704       }
2705       if(isvperiodic && correctV) {
2706
2707         if(thevmin < vmin)
2708           thevmin = vmin;
2709         if(thevmax > vmax)
2710           thevmax = vmax;
2711       }
2712     }
2713   }
2714 }
2715 //
2716 //
2717 // The block is dedicated to determine whether WLine [ifprm, ilprm]
2718 // crosses the degenerated zone on each given surface or not.
2719 // If Yes -> We will not use info about surfaces during approximation
2720 // because inside degenerated zone of the surface the approx. algo.
2721 // uses wrong values of normal, etc., and resulting curve will have
2722 // oscillations that we would not like to have.
2723
2724
2725
2726 static
2727   Standard_Boolean IsDegeneratedZone(const gp_Pnt2d& aP2d,
2728                                      const Handle(Geom_Surface)& aS,
2729                                      const Standard_Integer iDir);
2730 static
2731   Standard_Boolean IsPointInDegeneratedZone(const IntSurf_PntOn2S& aP2S,
2732                                             const TopoDS_Face& aF1,
2733                                             const TopoDS_Face& aF2);
2734 //=======================================================================
2735 //function :  NotUseSurfacesForApprox
2736 //purpose  :
2737 //=======================================================================
2738 Standard_Boolean NotUseSurfacesForApprox(const TopoDS_Face& aF1,
2739                                          const TopoDS_Face& aF2,
2740                                          const Handle(IntPatch_WLine)& WL,
2741                                          const Standard_Integer ifprm,
2742                                          const Standard_Integer ilprm)
2743 {
2744   Standard_Boolean bPInDZ;
2745
2746   Handle(IntSurf_LineOn2S) aLineOn2S=WL->Curve();
2747
2748   const IntSurf_PntOn2S& aP2Sfprm=aLineOn2S->Value(ifprm);
2749   bPInDZ=IsPointInDegeneratedZone(aP2Sfprm, aF1, aF2);
2750   if (bPInDZ) {
2751     return bPInDZ;
2752   }
2753
2754   const IntSurf_PntOn2S& aP2Slprm=aLineOn2S->Value(ilprm);
2755   bPInDZ=IsPointInDegeneratedZone(aP2Slprm, aF1, aF2);
2756
2757   return bPInDZ;
2758 }
2759 //=======================================================================
2760 //function : IsPointInDegeneratedZone
2761 //purpose  :
2762 //=======================================================================
2763 Standard_Boolean IsPointInDegeneratedZone(const IntSurf_PntOn2S& aP2S,
2764                                           const TopoDS_Face& aF1,
2765                                           const TopoDS_Face& aF2)
2766
2767 {
2768   Standard_Boolean bFlag=Standard_True;
2769   Standard_Real US11, US12, VS11, VS12, US21, US22, VS21, VS22;
2770   Standard_Real U1, V1, U2, V2, aDelta, aD;
2771   gp_Pnt2d aP2d;
2772
2773   Handle(Geom_Surface)aS1 = BRep_Tool::Surface(aF1);
2774   aS1->Bounds(US11, US12, VS11, VS12);
2775   GeomAdaptor_Surface aGAS1(aS1);
2776
2777   Handle(Geom_Surface)aS2 = BRep_Tool::Surface(aF2);
2778   aS1->Bounds(US21, US22, VS21, VS22);
2779   GeomAdaptor_Surface aGAS2(aS2);
2780   //
2781   //const gp_Pnt& aP=aP2S.Value();
2782   aP2S.Parameters(U1, V1, U2, V2);
2783   //
2784   aDelta=1.e-7;
2785   // Check on Surf 1
2786   aD=aGAS1.UResolution(aDelta);
2787   aP2d.SetCoord(U1, V1);
2788   if (fabs(U1-US11) < aD) {
2789     bFlag=IsDegeneratedZone(aP2d, aS1, 1);
2790     if (bFlag) {
2791       return bFlag;
2792     }
2793   }
2794   if (fabs(U1-US12) < aD) {
2795     bFlag=IsDegeneratedZone(aP2d, aS1, 1);
2796     if (bFlag) {
2797       return bFlag;
2798     }
2799   }
2800   aD=aGAS1.VResolution(aDelta);
2801   if (fabs(V1-VS11) < aDelta) {
2802     bFlag=IsDegeneratedZone(aP2d, aS1, 2);
2803     if (bFlag) {
2804       return bFlag;
2805     }
2806   }
2807   if (fabs(V1-VS12) < aDelta) {
2808     bFlag=IsDegeneratedZone(aP2d, aS1, 2);
2809     if (bFlag) {
2810       return bFlag;
2811     }
2812   }
2813   // Check on Surf 2
2814   aD=aGAS2.UResolution(aDelta);
2815   aP2d.SetCoord(U2, V2);
2816   if (fabs(U2-US21) < aDelta) {
2817     bFlag=IsDegeneratedZone(aP2d, aS2, 1);
2818     if (bFlag) {
2819       return bFlag;
2820     }
2821   }
2822   if (fabs(U2-US22) < aDelta) {
2823     bFlag=IsDegeneratedZone(aP2d, aS2, 1);
2824     if (bFlag) {
2825       return bFlag;
2826     }
2827   }
2828   aD=aGAS2.VResolution(aDelta);
2829   if (fabs(V2-VS21) < aDelta) {
2830     bFlag=IsDegeneratedZone(aP2d, aS2, 2);
2831     if (bFlag) {
2832       return bFlag;
2833     }
2834   }
2835   if (fabs(V2-VS22) < aDelta) {
2836     bFlag=IsDegeneratedZone(aP2d, aS2, 2);
2837     if (bFlag) {
2838       return bFlag;
2839     }
2840   }
2841   return !bFlag;
2842 }
2843
2844 //=======================================================================
2845 //function : IsDegeneratedZone
2846 //purpose  :
2847 //=======================================================================
2848 Standard_Boolean IsDegeneratedZone(const gp_Pnt2d& aP2d,
2849                                    const Handle(Geom_Surface)& aS,
2850                                    const Standard_Integer iDir)
2851 {
2852   Standard_Boolean bFlag=Standard_True;
2853   Standard_Real US1, US2, VS1, VS2, dY, dX, d1, d2, dD;
2854   Standard_Real aXm, aYm, aXb, aYb, aXe, aYe;
2855   aS->Bounds(US1, US2, VS1, VS2);
2856
2857   gp_Pnt aPm, aPb, aPe;
2858
2859   aXm=aP2d.X();
2860   aYm=aP2d.Y();
2861
2862   aS->D0(aXm, aYm, aPm);
2863
2864   dX=1.e-5;
2865   dY=1.e-5;
2866   dD=1.e-12;
2867
2868   if (iDir==1) {
2869     aXb=aXm;
2870     aXe=aXm;
2871     aYb=aYm-dY;
2872     if (aYb < VS1) {
2873       aYb=VS1;
2874     }
2875     aYe=aYm+dY;
2876     if (aYe > VS2) {
2877       aYe=VS2;
2878     }
2879     aS->D0(aXb, aYb, aPb);
2880     aS->D0(aXe, aYe, aPe);
2881
2882     d1=aPm.Distance(aPb);
2883     d2=aPm.Distance(aPe);
2884     if (d1 < dD && d2 < dD) {
2885       return bFlag;
2886     }
2887     return !bFlag;
2888   }
2889   //
2890   else if (iDir==2) {
2891     aYb=aYm;
2892     aYe=aYm;
2893     aXb=aXm-dX;
2894     if (aXb < US1) {
2895       aXb=US1;
2896     }
2897     aXe=aXm+dX;
2898     if (aXe > US2) {
2899       aXe=US2;
2900     }
2901     aS->D0(aXb, aYb, aPb);
2902     aS->D0(aXe, aYe, aPe);
2903
2904     d1=aPm.Distance(aPb);
2905     d2=aPm.Distance(aPe);
2906     if (d1 < dD && d2 < dD) {
2907       return bFlag;
2908     }
2909     return !bFlag;
2910   }
2911   return !bFlag;
2912 }
2913
2914 //=========================================================================
2915 // static function : ComputePurgedWLine
2916 // purpose : Removes equal points (leave one of equal points) from theWLine
2917 //           and recompute vertex parameters.
2918 //           Returns new WLine or null WLine if the number
2919 //           of the points is less than 2.
2920 //=========================================================================
2921 Handle(IntPatch_WLine) ComputePurgedWLine(const Handle(IntPatch_WLine)& theWLine) {
2922
2923   Standard_Integer i, k, v, nb, nbvtx;
2924   Handle(IntPatch_WLine) aResult;
2925   nbvtx = theWLine->NbVertex();
2926   nb = theWLine->NbPnts();
2927   if (nb==2) {
2928     const IntSurf_PntOn2S& p1 = theWLine->Point(1);
2929     const IntSurf_PntOn2S& p2 = theWLine->Point(2);
2930     if(p1.Value().IsEqual(p2.Value(), gp::Resolution())) {
2931       return aResult;
2932     }
2933   }
2934   //
2935   Handle(IntPatch_WLine) aLocalWLine;
2936   Handle(IntPatch_WLine) aTmpWLine = theWLine;
2937   Handle(IntSurf_LineOn2S) aLineOn2S = new IntSurf_LineOn2S();
2938   aLocalWLine = new IntPatch_WLine(aLineOn2S, Standard_False);
2939   for(i = 1; i <= nb; i++) {
2940     aLineOn2S->Add(theWLine->Point(i));
2941   }
2942
2943   for(v = 1; v <= nbvtx; v++) {
2944     aLocalWLine->AddVertex(theWLine->Vertex(v));
2945   }
2946
2947   for(i = 1; i <= aLineOn2S->NbPoints(); i++) {
2948     Standard_Integer aStartIndex = i + 1;
2949     Standard_Integer anEndIndex = i + 5;
2950     nb = aLineOn2S->NbPoints();
2951     anEndIndex = (anEndIndex > nb) ? nb : anEndIndex;
2952
2953     if((aStartIndex > nb) || (anEndIndex <= 1)) {
2954       continue;
2955     }
2956     k = aStartIndex;
2957
2958     while(k <= anEndIndex) {
2959
2960       if(i != k) {
2961         IntSurf_PntOn2S p1 = aLineOn2S->Value(i);
2962         IntSurf_PntOn2S p2 = aLineOn2S->Value(k);
2963
2964         if(p1.Value().IsEqual(p2.Value(), gp::Resolution())) {
2965           aTmpWLine = aLocalWLine;
2966           aLocalWLine = new IntPatch_WLine(aLineOn2S, Standard_False);
2967
2968           for(v = 1; v <= aTmpWLine->NbVertex(); v++) {
2969             IntPatch_Point aVertex = aTmpWLine->Vertex(v);
2970             Standard_Integer avertexindex = (Standard_Integer)aVertex.ParameterOnLine();
2971
2972             if(avertexindex >= k) {
2973               aVertex.SetParameter(aVertex.ParameterOnLine() - 1.);
2974             }
2975             aLocalWLine->AddVertex(aVertex);
2976           }
2977           aLineOn2S->RemovePoint(k);
2978           anEndIndex--;
2979           continue;
2980         }
2981       }
2982       k++;
2983     }
2984   }
2985
2986   if(aLineOn2S->NbPoints() > 1) {
2987     aResult = aLocalWLine;
2988   }
2989   return aResult;
2990 }
2991
2992 //=======================================================================
2993 //function : TolR3d
2994 //purpose  :
2995 //=======================================================================
2996 void TolR3d(const TopoDS_Face& aF1,
2997             const TopoDS_Face& aF2,
2998             Standard_Real& myTolReached3d)
2999 {
3000   Standard_Real aTolF1, aTolF2, aTolFMax, aTolTresh;
3001
3002   aTolTresh=2.999999e-3;
3003   aTolF1 = BRep_Tool::Tolerance(aF1);
3004   aTolF2 = BRep_Tool::Tolerance(aF2);
3005   aTolFMax=Max(aTolF1, aTolF2);
3006
3007   if (aTolFMax>aTolTresh) {
3008     myTolReached3d=aTolFMax;
3009   }
3010 }
3011 //=======================================================================
3012 //function : IsPointOnBoundary
3013 //purpose  :
3014 //=======================================================================
3015 Standard_Boolean IsPointOnBoundary(const Standard_Real theParameter,
3016                                    const Standard_Real theFirstBoundary,
3017                                    const Standard_Real theSecondBoundary,
3018                                    const Standard_Real theResolution,
3019                                    Standard_Boolean&   IsOnFirstBoundary)
3020 {
3021   Standard_Boolean bRet;
3022   Standard_Integer i;
3023   Standard_Real adist;
3024   //
3025   bRet=Standard_False;
3026   for(i = 0; i < 2; ++i) {
3027     IsOnFirstBoundary = (i == 0);
3028     if (IsOnFirstBoundary) {
3029       adist = fabs(theParameter - theFirstBoundary);
3030     }
3031     else {
3032       adist = fabs(theParameter - theSecondBoundary);
3033     }
3034     if(adist < theResolution) {
3035       return !bRet;
3036     }
3037   }
3038   return bRet;
3039 }
3040 // ------------------------------------------------------------------------------------------------
3041 // static function: FindPoint
3042 // purpose:
3043 // ------------------------------------------------------------------------------------------------
3044 Standard_Boolean FindPoint(const gp_Pnt2d&     theFirstPoint,
3045                            const gp_Pnt2d&     theLastPoint,
3046                            const Standard_Real theUmin,
3047                            const Standard_Real theUmax,
3048                            const Standard_Real theVmin,
3049                            const Standard_Real theVmax,
3050                            gp_Pnt2d&           theNewPoint) {
3051
3052   gp_Vec2d aVec(theFirstPoint, theLastPoint);
3053   Standard_Integer i = 0, j = 0;
3054
3055   for(i = 0; i < 4; i++) {
3056     gp_Vec2d anOtherVec;
3057     gp_Vec2d anOtherVecNormal;
3058     gp_Pnt2d aprojpoint = theLastPoint;
3059
3060     if((i % 2) == 0) {
3061       anOtherVec.SetX(0.);
3062       anOtherVec.SetY(1.);
3063       anOtherVecNormal.SetX(1.);
3064       anOtherVecNormal.SetY(0.);
3065
3066       if(i < 2)
3067         aprojpoint.SetX(theUmin);
3068       else
3069         aprojpoint.SetX(theUmax);
3070     }
3071     else {
3072       anOtherVec.SetX(1.);
3073       anOtherVec.SetY(0.);
3074       anOtherVecNormal.SetX(0.);
3075       anOtherVecNormal.SetY(1.);
3076
3077       if(i < 2)
3078         aprojpoint.SetY(theVmin);
3079       else
3080         aprojpoint.SetY(theVmax);
3081     }
3082     gp_Vec2d anormvec = aVec;
3083     anormvec.Normalize();
3084     RefineVector(anormvec);
3085     Standard_Real adot1 = anormvec.Dot(anOtherVecNormal);
3086
3087     if(fabs(adot1) < Precision::Angular())
3088       continue;
3089     Standard_Real adist = 0.;
3090     Standard_Boolean bIsOut = Standard_False;
3091
3092     if((i % 2) == 0) {
3093       adist = (i < 2) ? fabs(theLastPoint.X() - theUmin) : fabs(theLastPoint.X() - theUmax);
3094       bIsOut = (i < 2) ? (theLastPoint.X() < theUmin) : (theLastPoint.X() > theUmax);
3095     }
3096     else {
3097       adist = (i < 2) ? fabs(theLastPoint.Y() - theVmin) : fabs(theLastPoint.Y() - theVmax);
3098       bIsOut = (i < 2) ? (theLastPoint.Y() < theVmin) : (theLastPoint.Y() > theVmax);
3099     }
3100     Standard_Real anoffset = adist * anOtherVec.Dot(anormvec) / adot1;
3101
3102     for(j = 0; j < 2; j++) {
3103       anoffset = (j == 0) ? anoffset : -anoffset;
3104       gp_Pnt2d acurpoint(aprojpoint.XY() + (anOtherVec.XY()*anoffset));
3105       gp_Vec2d acurvec(theLastPoint, acurpoint);
3106       if ( bIsOut )
3107         acurvec.Reverse();
3108
3109       Standard_Real aDotX, anAngleX;
3110       //
3111       aDotX = aVec.Dot(acurvec);
3112       anAngleX = aVec.Angle(acurvec);
3113       //
3114       if(aDotX > 0. && fabs(anAngleX) < Precision::PConfusion()) {
3115         if((i % 2) == 0) {
3116           if((acurpoint.Y() >= theVmin) &&
3117              (acurpoint.Y() <= theVmax)) {
3118             theNewPoint = acurpoint;
3119             return Standard_True;
3120           }
3121         }
3122         else {
3123           if((acurpoint.X() >= theUmin) &&
3124              (acurpoint.X() <= theUmax)) {
3125             theNewPoint = acurpoint;
3126             return Standard_True;
3127           }
3128         }
3129       }
3130     }
3131   }
3132   return Standard_False;
3133 }
3134
3135
3136 // ------------------------------------------------------------------------------------------------
3137 // static function: FindPoint
3138 // purpose: Find point on the boundary of radial tangent zone
3139 // ------------------------------------------------------------------------------------------------
3140 Standard_Boolean FindPoint(const gp_Pnt2d&     theFirstPoint,
3141                            const gp_Pnt2d&     theLastPoint,
3142                            const Standard_Real theUmin,
3143                            const Standard_Real theUmax,
3144                            const Standard_Real theVmin,
3145                            const Standard_Real theVmax,
3146                            const gp_Pnt2d&     theTanZoneCenter,
3147                            const Standard_Real theZoneRadius,
3148                            Handle(GeomAdaptor_HSurface) theGASurface,
3149                            gp_Pnt2d&           theNewPoint) {
3150   theNewPoint = theLastPoint;
3151
3152   if ( !IsInsideTanZone( theLastPoint, theTanZoneCenter, theZoneRadius, theGASurface) )
3153     return Standard_False;
3154
3155   Standard_Real aUResolution = theGASurface->UResolution( theZoneRadius );
3156   Standard_Real aVResolution = theGASurface->VResolution( theZoneRadius );
3157
3158   Standard_Real aRadius = ( aUResolution < aVResolution ) ? aUResolution : aVResolution;
3159   gp_Ax22d anAxis( theTanZoneCenter, gp_Dir2d(1, 0), gp_Dir2d(0, 1) );
3160   gp_Circ2d aCircle( anAxis, aRadius );
3161
3162   //
3163   gp_Vec2d aDir( theLastPoint.XY() - theFirstPoint.XY() );
3164   Standard_Real aLength = aDir.Magnitude();
3165   if ( aLength <= gp::Resolution() )
3166     return Standard_False;
3167   gp_Lin2d aLine( theFirstPoint, aDir );
3168
3169   //
3170   Handle(Geom2d_Line) aCLine = new Geom2d_Line( aLine );
3171   Handle(Geom2d_TrimmedCurve) aC1 = new Geom2d_TrimmedCurve( aCLine, 0, aLength );
3172   Handle(Geom2d_Circle) aC2 = new Geom2d_Circle( aCircle );
3173
3174   Standard_Real aTol = aRadius * 0.001;
3175   aTol = ( aTol < Precision::PConfusion() ) ? Precision::PConfusion() : aTol;
3176
3177   Geom2dAPI_InterCurveCurve anIntersector;
3178   anIntersector.Init( aC1, aC2, aTol );
3179
3180   if ( anIntersector.NbPoints() == 0 )
3181     return Standard_False;
3182
3183   Standard_Boolean aFound = Standard_False;
3184   Standard_Real aMinDist = aLength * aLength;
3185   Standard_Integer i = 0;
3186   for ( i = 1; i <= anIntersector.NbPoints(); i++ ) {
3187     gp_Pnt2d aPInt = anIntersector.Point( i );
3188     if ( aPInt.SquareDistance( theFirstPoint ) < aMinDist ) {
3189       if ( ( aPInt.X() >= theUmin ) && ( aPInt.X() <= theUmax ) &&
3190            ( aPInt.Y() >= theVmin ) && ( aPInt.Y() <= theVmax ) ) {
3191         theNewPoint = aPInt;
3192         aFound = Standard_True;
3193       }
3194     }
3195   }
3196
3197   return aFound;
3198 }
3199
3200 // ------------------------------------------------------------------------------------------------
3201 // static function: IsInsideTanZone
3202 // purpose: Check if point is inside a radial tangent zone
3203 // ------------------------------------------------------------------------------------------------
3204 Standard_Boolean IsInsideTanZone(const gp_Pnt2d&     thePoint,
3205                                  const gp_Pnt2d&     theTanZoneCenter,
3206                                  const Standard_Real theZoneRadius,
3207                                  Handle(GeomAdaptor_HSurface) theGASurface) {
3208
3209   Standard_Real aUResolution = theGASurface->UResolution( theZoneRadius );
3210   Standard_Real aVResolution = theGASurface->VResolution( theZoneRadius );
3211   Standard_Real aRadiusSQR = ( aUResolution < aVResolution ) ? aUResolution : aVResolution;
3212   aRadiusSQR *= aRadiusSQR;
3213   if ( thePoint.SquareDistance( theTanZoneCenter ) <= aRadiusSQR )
3214     return Standard_True;
3215   return Standard_False;
3216 }
3217
3218 // ------------------------------------------------------------------------------------------------
3219 // static function: CheckTangentZonesExist
3220 // purpose: Check if tangent zone exists
3221 // ------------------------------------------------------------------------------------------------
3222 Standard_Boolean CheckTangentZonesExist( const Handle(GeomAdaptor_HSurface)& theSurface1,
3223                                         const Handle(GeomAdaptor_HSurface)&  theSurface2 )
3224 {
3225   if ( ( theSurface1->GetType() != GeomAbs_Torus ) ||
3226       ( theSurface2->GetType() != GeomAbs_Torus ) )
3227     return Standard_False;
3228
3229   gp_Torus aTor1 = theSurface1->Torus();
3230   gp_Torus aTor2 = theSurface2->Torus();
3231
3232   if ( aTor1.Location().Distance( aTor2.Location() ) > Precision::Confusion() )
3233     return Standard_False;
3234
3235   if ( ( fabs( aTor1.MajorRadius() - aTor2.MajorRadius() ) > Precision::Confusion() ) ||
3236        ( fabs( aTor1.MinorRadius() - aTor2.MinorRadius() ) > Precision::Confusion() ) )
3237     return Standard_False;
3238
3239   if ( ( aTor1.MajorRadius() < aTor1.MinorRadius() ) ||
3240        ( aTor2.MajorRadius() < aTor2.MinorRadius() ) )
3241     return Standard_False;
3242   return Standard_True;
3243 }
3244
3245 // ------------------------------------------------------------------------------------------------
3246 // static function: ComputeTangentZones
3247 // purpose:
3248 // ------------------------------------------------------------------------------------------------
3249 Standard_Integer ComputeTangentZones( const Handle(GeomAdaptor_HSurface)& theSurface1,
3250                                      const Handle(GeomAdaptor_HSurface)&  theSurface2,
3251                                      const TopoDS_Face&                   theFace1,
3252                                      const TopoDS_Face&                   theFace2,
3253                                      Handle(TColgp_HArray1OfPnt2d)&       theResultOnS1,
3254                                      Handle(TColgp_HArray1OfPnt2d)&       theResultOnS2,
3255                                      Handle(TColStd_HArray1OfReal)&       theResultRadius,
3256                                      const Handle(IntTools_Context)& aContext)
3257 {
3258   Standard_Integer aResult = 0;
3259   if ( !CheckTangentZonesExist( theSurface1, theSurface2 ) )
3260     return aResult;
3261
3262
3263   TColgp_SequenceOfPnt2d aSeqResultS1, aSeqResultS2;
3264   TColStd_SequenceOfReal aSeqResultRad;
3265
3266   gp_Torus aTor1 = theSurface1->Torus();
3267   gp_Torus aTor2 = theSurface2->Torus();
3268
3269   gp_Ax2 anax1( aTor1.Location(), aTor1.Axis().Direction() );
3270   gp_Ax2 anax2( aTor2.Location(), aTor2.Axis().Direction() );
3271   Standard_Integer j = 0;
3272
3273   for ( j = 0; j < 2; j++ ) {
3274     Standard_Real aCoef = ( j == 0 ) ? -1 : 1;
3275     Standard_Real aRadius1 = fabs(aTor1.MajorRadius() + aCoef * aTor1.MinorRadius());
3276     Standard_Real aRadius2 = fabs(aTor2.MajorRadius() + aCoef * aTor2.MinorRadius());
3277
3278     gp_Circ aCircle1( anax1, aRadius1 );
3279     gp_Circ aCircle2( anax2, aRadius2 );
3280
3281     // roughly compute radius of tangent zone for perpendicular case
3282     Standard_Real aCriteria = Precision::Confusion() * 0.5;
3283
3284     Standard_Real aT1 = aCriteria;
3285     Standard_Real aT2 = aCriteria;
3286     if ( j == 0 ) {
3287       // internal tangency
3288       Standard_Real aR = ( aRadius1 > aTor2.MinorRadius() ) ? aRadius1 : aTor2.MinorRadius();
3289       //aT1 = aCriteria * aCriteria + aR * aR - ( aR - aCriteria ) * ( aR - aCriteria );
3290       aT1 = 2. * aR * aCriteria;
3291       aT2 = aT1;
3292     }
3293     else {
3294       // external tangency
3295       Standard_Real aRb = ( aRadius1 > aTor2.MinorRadius() ) ? aRadius1 : aTor2.MinorRadius();
3296       Standard_Real aRm = ( aRadius1 < aTor2.MinorRadius() ) ? aRadius1 : aTor2.MinorRadius();
3297       Standard_Real aDelta = aRb - aCriteria;
3298       aDelta *= aDelta;
3299       aDelta -= aRm * aRm;
3300       aDelta /= 2. * (aRb - aRm);
3301       aDelta -= 0.5 * (aRb - aRm);
3302
3303       aT1 = 2. * aRm * (aRm - aDelta);
3304       aT2 = aT1;
3305     }
3306     aCriteria = ( aT1 > aT2) ? aT1 : aT2;
3307     if ( aCriteria > 0 )
3308       aCriteria = sqrt( aCriteria );
3309
3310     if ( aCriteria > 0.5 * aTor1.MinorRadius() ) {
3311       // too big zone -> drop to minimum
3312       aCriteria = Precision::Confusion();
3313     }
3314
3315     GeomAdaptor_Curve aC1( new Geom_Circle(aCircle1) );
3316     GeomAdaptor_Curve aC2( new Geom_Circle(aCircle2) );
3317     Extrema_ExtCC anExtrema(aC1, aC2, 0, 2. * M_PI, 0, 2. * M_PI,
3318                             Precision::PConfusion(), Precision::PConfusion());
3319
3320     if ( anExtrema.IsDone() ) {
3321
3322       Standard_Integer i = 0;
3323       for ( i = 1; i <= anExtrema.NbExt(); i++ ) {
3324         if ( anExtrema.SquareDistance(i) > aCriteria * aCriteria )
3325           continue;
3326
3327         Extrema_POnCurv P1, P2;
3328         anExtrema.Points( i, P1, P2 );
3329
3330         Standard_Boolean bFoundResult = Standard_True;
3331         gp_Pnt2d pr1, pr2;
3332
3333         Standard_Integer surfit = 0;
3334         for ( surfit = 0; surfit < 2; surfit++ ) {
3335           GeomAPI_ProjectPointOnSurf& aProjector =
3336             (surfit == 0) ? aContext->ProjPS(theFace1) : aContext->ProjPS(theFace2);
3337
3338           gp_Pnt aP3d = (surfit == 0) ? P1.Value() : P2.Value();
3339           aProjector.Perform(aP3d);
3340
3341           if(!aProjector.IsDone())
3342             bFoundResult = Standard_False;
3343           else {
3344             if(aProjector.LowerDistance() > aCriteria) {
3345               bFoundResult = Standard_False;
3346             }
3347             else {
3348               Standard_Real foundU = 0, foundV = 0;
3349               aProjector.LowerDistanceParameters(foundU, foundV);
3350               if ( surfit == 0 )
3351                 pr1 = gp_Pnt2d( foundU, foundV );
3352               else
3353                 pr2 = gp_Pnt2d( foundU, foundV );
3354             }
3355           }
3356         }
3357         if ( bFoundResult ) {
3358           aSeqResultS1.Append( pr1 );
3359           aSeqResultS2.Append( pr2 );
3360           aSeqResultRad.Append( aCriteria );
3361
3362           // torus is u and v periodic
3363           const Standard_Real twoPI = M_PI + M_PI;
3364           Standard_Real arr1tmp[2] = {pr1.X(), pr1.Y()};
3365           Standard_Real arr2tmp[2] = {pr2.X(), pr2.Y()};
3366
3367           // iteration on period bounds
3368           for ( Standard_Integer k1 = 0; k1 < 2; k1++ ) {
3369             Standard_Real aBound = ( k1 == 0 ) ? 0 : twoPI;
3370             Standard_Real aShift = ( k1 == 0 ) ? twoPI : -twoPI;
3371
3372             // iteration on surfaces
3373             for ( Standard_Integer k2 = 0; k2 < 2; k2++ ) {
3374               Standard_Real* arr1 = ( k2 == 0 ) ? arr1tmp : arr2tmp;
3375               Standard_Real* arr2 = ( k2 != 0 ) ? arr1tmp : arr2tmp;
3376               TColgp_SequenceOfPnt2d& aSeqS1 = ( k2 == 0 ) ? aSeqResultS1 : aSeqResultS2;
3377               TColgp_SequenceOfPnt2d& aSeqS2 = ( k2 != 0 ) ? aSeqResultS1 : aSeqResultS2;
3378
3379               if (fabs(arr1[0] - aBound) < Precision::PConfusion()) {
3380                 aSeqS1.Append( gp_Pnt2d( arr1[0] + aShift, arr1[1] ) );
3381                 aSeqS2.Append( gp_Pnt2d( arr2[0], arr2[1] ) );
3382                 aSeqResultRad.Append( aCriteria );
3383               }
3384               if (fabs(arr1[1] - aBound) < Precision::PConfusion()) {
3385                 aSeqS1.Append( gp_Pnt2d( arr1[0], arr1[1] + aShift) );
3386                 aSeqS2.Append( gp_Pnt2d( arr2[0], arr2[1] ) );
3387                 aSeqResultRad.Append( aCriteria );
3388               }
3389             }
3390           } //
3391         }
3392       }
3393     }
3394   }
3395   aResult = aSeqResultRad.Length();
3396
3397   if ( aResult > 0 ) {
3398     theResultOnS1 = new TColgp_HArray1OfPnt2d( 1, aResult );
3399     theResultOnS2 = new TColgp_HArray1OfPnt2d( 1, aResult );
3400     theResultRadius = new TColStd_HArray1OfReal( 1, aResult );
3401
3402     for ( Standard_Integer i = 1 ; i <= aResult; i++ ) {
3403       theResultOnS1->SetValue( i, aSeqResultS1.Value(i) );
3404       theResultOnS2->SetValue( i, aSeqResultS2.Value(i) );
3405       theResultRadius->SetValue( i, aSeqResultRad.Value(i) );
3406     }
3407   }
3408   return aResult;
3409 }
3410
3411 // ------------------------------------------------------------------------------------------------
3412 // static function: AdjustByNeighbour
3413 // purpose:
3414 // ------------------------------------------------------------------------------------------------
3415 gp_Pnt2d AdjustByNeighbour(const gp_Pnt2d&     theaNeighbourPoint,
3416                            const gp_Pnt2d&     theOriginalPoint,
3417                            Handle(GeomAdaptor_HSurface) theGASurface) {
3418
3419   gp_Pnt2d ap1 = theaNeighbourPoint;
3420   gp_Pnt2d ap2 = theOriginalPoint;
3421
3422   if ( theGASurface->IsUPeriodic() ) {
3423     Standard_Real aPeriod     = theGASurface->UPeriod();
3424     gp_Pnt2d aPTest = ap2;
3425     Standard_Real aSqDistMin = 1.e+100;
3426
3427     for ( Standard_Integer pIt = -1; pIt <= 1; pIt++) {
3428       aPTest.SetX( theOriginalPoint.X() + aPeriod * pIt );
3429       Standard_Real dd = ap1.SquareDistance( aPTest );
3430
3431       if ( dd < aSqDistMin ) {
3432         ap2 = aPTest;
3433         aSqDistMin = dd;
3434       }
3435     }
3436   }
3437   if ( theGASurface->IsVPeriodic() ) {
3438     Standard_Real aPeriod     = theGASurface->VPeriod();
3439     gp_Pnt2d aPTest = ap2;
3440     Standard_Real aSqDistMin = 1.e+100;
3441
3442     for ( Standard_Integer pIt = -1; pIt <= 1; pIt++) {
3443       aPTest.SetY( theOriginalPoint.Y() + aPeriod * pIt );
3444       Standard_Real dd = ap1.SquareDistance( aPTest );
3445
3446       if ( dd < aSqDistMin ) {
3447         ap2 = aPTest;
3448         aSqDistMin = dd;
3449       }
3450     }
3451   }
3452   return ap2;
3453 }
3454
3455 // ------------------------------------------------------------------------------------------------
3456 //function: DecompositionOfWLine
3457 // purpose:
3458 // ------------------------------------------------------------------------------------------------
3459 Standard_Boolean DecompositionOfWLine(const Handle(IntPatch_WLine)& theWLine,
3460                                       const Handle(GeomAdaptor_HSurface)&            theSurface1,
3461                                       const Handle(GeomAdaptor_HSurface)&            theSurface2,
3462                                       const TopoDS_Face&                             theFace1,
3463                                       const TopoDS_Face&                             theFace2,
3464                                       const GeomInt_LineConstructor&                 theLConstructor,
3465                                       const Standard_Boolean                         theAvoidLConstructor,
3466                                       IntPatch_SequenceOfLine&                       theNewLines,
3467                                       Standard_Real&                                 theReachedTol3d,
3468                                       const Handle(IntTools_Context)& aContext)
3469 {
3470
3471   Standard_Boolean bRet, bAvoidLineConstructor;
3472   Standard_Integer aNbPnts, aNbParts;
3473   //
3474   bRet=Standard_False;
3475   aNbPnts=theWLine->NbPnts();
3476   bAvoidLineConstructor=theAvoidLConstructor;
3477   //
3478   if(!aNbPnts){
3479     return bRet;
3480   }
3481   if (!bAvoidLineConstructor) {
3482     aNbParts=theLConstructor.NbParts();
3483     if (!aNbParts) {
3484       return bRet;
3485     }
3486   }
3487   //
3488   Standard_Boolean bIsPrevPointOnBoundary, bIsPointOnBoundary, bIsCurrentPointOnBoundary;
3489   Standard_Integer nblines, pit, i, j;
3490   Standard_Real aTol;
3491   TColStd_Array1OfListOfInteger anArrayOfLines(1, aNbPnts);
3492   TColStd_Array1OfInteger       anArrayOfLineType(1, aNbPnts);
3493   TColStd_ListOfInteger aListOfPointIndex;
3494
3495   Handle(TColgp_HArray1OfPnt2d) aTanZoneS1;
3496   Handle(TColgp_HArray1OfPnt2d) aTanZoneS2;
3497   Handle(TColStd_HArray1OfReal) aTanZoneRadius;
3498   Standard_Integer aNbZone = ComputeTangentZones( theSurface1, theSurface2, theFace1, theFace2,
3499                                                  aTanZoneS1, aTanZoneS2, aTanZoneRadius, aContext);
3500
3501   //
3502   nblines=0;
3503   aTol=Precision::Confusion();
3504   aTol=0.5*aTol;
3505   bIsPrevPointOnBoundary=Standard_False;
3506   bIsPointOnBoundary=Standard_False;
3507   //
3508   // 1. ...
3509   //
3510   // Points
3511   for(pit = 1; pit <= aNbPnts; ++pit) {
3512     Standard_Boolean bIsOnFirstBoundary, isperiodic;
3513     Standard_Real aResolution, aPeriod, alowerboundary, aupperboundary, U, V;
3514     Standard_Real aParameter, anoffset, anAdjustPar;
3515     Standard_Real umin, umax, vmin, vmax;
3516     //
3517     bIsCurrentPointOnBoundary = Standard_False;
3518     const IntSurf_PntOn2S& aPoint = theWLine->Point(pit);
3519     //
3520     // Surface
3521     for(i = 0; i < 2; ++i) {
3522       Handle(GeomAdaptor_HSurface) aGASurface = (!i) ? theSurface1 : theSurface2;
3523       aGASurface->ChangeSurface().Surface()->Bounds(umin, umax, vmin, vmax);
3524       if(!i) {
3525         aPoint.ParametersOnS1(U, V);
3526       }
3527       else {
3528         aPoint.ParametersOnS2(U, V);
3529       }
3530       // U, V
3531       for(j = 0; j < 2; j++) {
3532         isperiodic = (!j) ? aGASurface->IsUPeriodic() : aGASurface->IsVPeriodic();
3533         if(!isperiodic){
3534           continue;
3535         }
3536         //
3537         if (!j) {
3538           aResolution=aGASurface->UResolution(aTol);
3539           aPeriod=aGASurface->UPeriod();
3540           alowerboundary=umin;
3541           aupperboundary=umax;
3542           aParameter=U;
3543         }
3544         else {
3545           aResolution=aGASurface->VResolution(aTol);
3546           aPeriod=aGASurface->VPeriod();
3547           alowerboundary=vmin;
3548           aupperboundary=vmax;
3549           aParameter=V;
3550         }
3551
3552         GeomInt::AdjustPeriodic(aParameter,
3553                                        alowerboundary,
3554                                        aupperboundary,
3555                                        aPeriod,
3556                                        anAdjustPar,
3557                                        anoffset);
3558         //
3559         bIsOnFirstBoundary = Standard_True;// ?
3560         bIsPointOnBoundary=
3561           IsPointOnBoundary(anAdjustPar,
3562                             alowerboundary,
3563                             aupperboundary,
3564                             aResolution,
3565                             bIsOnFirstBoundary);
3566         //
3567         if(bIsPointOnBoundary) {
3568           bIsCurrentPointOnBoundary = Standard_True;
3569           break;
3570         }
3571         else {
3572           // check if a point belong to a tangent zone. Begin
3573           Standard_Integer zIt = 0;
3574           for ( zIt = 1; zIt <= aNbZone; zIt++ ) {
3575             gp_Pnt2d aPZone = (i == 0) ? aTanZoneS1->Value(zIt) : aTanZoneS2->Value(zIt);
3576             Standard_Real aZoneRadius = aTanZoneRadius->Value(zIt);
3577
3578             if ( IsInsideTanZone(gp_Pnt2d( U, V ), aPZone, aZoneRadius, aGASurface ) ) {
3579               // set boundary flag to split the curve by a tangent zone
3580               bIsPointOnBoundary = Standard_True;
3581               bIsCurrentPointOnBoundary = Standard_True;
3582               if ( theReachedTol3d < aZoneRadius ) {
3583                 theReachedTol3d = aZoneRadius;
3584               }
3585               break;
3586             }
3587           }
3588         }
3589       }//for(j = 0; j < 2; j++) {
3590
3591       if(bIsCurrentPointOnBoundary){
3592         break;
3593       }
3594     }//for(i = 0; i < 2; ++i) {
3595     //
3596     if((bIsCurrentPointOnBoundary != bIsPrevPointOnBoundary)) {
3597       if(!aListOfPointIndex.IsEmpty()) {
3598         nblines++;
3599         anArrayOfLines.SetValue(nblines, aListOfPointIndex);
3600         anArrayOfLineType.SetValue(nblines, bIsPrevPointOnBoundary);
3601         aListOfPointIndex.Clear();
3602       }
3603       bIsPrevPointOnBoundary = bIsCurrentPointOnBoundary;
3604     }
3605     aListOfPointIndex.Append(pit);
3606   } //for(pit = 1; pit <= aNbPnts; ++pit) {
3607   //
3608   if(!aListOfPointIndex.IsEmpty()) {
3609     nblines++;
3610     anArrayOfLines.SetValue(nblines, aListOfPointIndex);
3611     anArrayOfLineType.SetValue(nblines, bIsPrevPointOnBoundary);
3612     aListOfPointIndex.Clear();
3613   }
3614   //
3615   if(nblines<=1) {
3616     return bRet; //Standard_False;
3617   }
3618   //
3619   //
3620   // 2. Correct wlines.begin
3621   TColStd_Array1OfListOfInteger anArrayOfLineEnds(1, nblines);
3622   Handle(IntSurf_LineOn2S) aSeqOfPntOn2S = new IntSurf_LineOn2S();
3623   //
3624   for(i = 1; i <= nblines; i++) {
3625     if(anArrayOfLineType.Value(i) != 0) {
3626       continue;
3627     }
3628     const TColStd_ListOfInteger& aListOfIndex = anArrayOfLines.Value(i);
3629     TColStd_ListOfInteger aListOfFLIndex;
3630
3631     for(j = 0; j < 2; j++) {
3632       Standard_Integer aneighbourindex = (j == 0) ? (i - 1) : (i + 1);
3633
3634       if((aneighbourindex < 1) || (aneighbourindex > nblines))
3635         continue;
3636
3637       if(anArrayOfLineType.Value(aneighbourindex) == 0)
3638         continue;
3639       const TColStd_ListOfInteger& aNeighbour = anArrayOfLines.Value(aneighbourindex);
3640       Standard_Integer anIndex = (j == 0) ? aNeighbour.Last() : aNeighbour.First();
3641       const IntSurf_PntOn2S& aPoint = theWLine->Point(anIndex);
3642
3643       IntSurf_PntOn2S aNewP = aPoint;
3644       if(aListOfIndex.Extent() < 2) {
3645         aSeqOfPntOn2S->Add(aNewP);
3646         aListOfFLIndex.Append(aSeqOfPntOn2S->NbPoints());
3647         continue;
3648       }
3649       //
3650       Standard_Integer iFirst = aListOfIndex.First();
3651       Standard_Integer iLast  = aListOfIndex.Last();
3652       //
3653       for(Standard_Integer surfit = 0; surfit < 2; surfit++) {
3654
3655         Handle(GeomAdaptor_HSurface) aGASurface = (surfit == 0) ? theSurface1 : theSurface2;
3656         Standard_Real umin=0., umax=0., vmin=0., vmax=0.;
3657         aGASurface->ChangeSurface().Surface()->Bounds(umin, umax, vmin, vmax);
3658         Standard_Real U=0., V=0.;
3659
3660         if(surfit == 0)
3661           aNewP.ParametersOnS1(U, V);
3662         else
3663           aNewP.ParametersOnS2(U, V);
3664         Standard_Integer nbboundaries = 0;
3665
3666         Standard_Boolean bIsNearBoundary = Standard_False;
3667         Standard_Integer aZoneIndex = 0;
3668         Standard_Integer bIsUBoundary = Standard_False; // use if nbboundaries == 1
3669         Standard_Integer bIsFirstBoundary = Standard_False; // use if nbboundaries == 1
3670
3671
3672         for(Standard_Integer parit = 0; parit < 2; parit++) {
3673           Standard_Boolean isperiodic = (parit == 0) ? aGASurface->IsUPeriodic() : aGASurface->IsVPeriodic();
3674
3675           Standard_Real aResolution = (parit == 0) ? aGASurface->UResolution(aTol) : aGASurface->VResolution(aTol);
3676           Standard_Real alowerboundary = (parit == 0) ? umin : vmin;
3677           Standard_Real aupperboundary = (parit == 0) ? umax : vmax;
3678
3679           Standard_Real aParameter = (parit == 0) ? U : V;
3680           Standard_Boolean bIsOnFirstBoundary = Standard_True;
3681
3682           if(!isperiodic) {
3683             bIsPointOnBoundary=
3684               IsPointOnBoundary(aParameter, alowerboundary, aupperboundary, aResolution, bIsOnFirstBoundary);
3685             if(bIsPointOnBoundary) {
3686               bIsUBoundary = (parit == 0);
3687               bIsFirstBoundary = bIsOnFirstBoundary;
3688               nbboundaries++;
3689             }
3690           }
3691           else {
3692             Standard_Real aPeriod     = (parit == 0) ? aGASurface->UPeriod() : aGASurface->VPeriod();
3693             Standard_Real anoffset, anAdjustPar;
3694             GeomInt::AdjustPeriodic(aParameter, alowerboundary, aupperboundary,
3695                                            aPeriod, anAdjustPar, anoffset);
3696
3697             bIsPointOnBoundary=
3698               IsPointOnBoundary(anAdjustPar, alowerboundary, aupperboundary, aResolution, bIsOnFirstBoundary);
3699             if(bIsPointOnBoundary) {
3700               bIsUBoundary = (parit == 0);
3701               bIsFirstBoundary = bIsOnFirstBoundary;
3702               nbboundaries++;
3703             }
3704             else {
3705             //check neighbourhood of boundary
3706             Standard_Real anEpsilon = aResolution * 100.;
3707             Standard_Real aPart = ( aupperboundary - alowerboundary ) * 0.1;
3708             anEpsilon = ( anEpsilon > aPart ) ? aPart : anEpsilon;
3709
3710               bIsNearBoundary = IsPointOnBoundary(anAdjustPar, alowerboundary, aupperboundary,
3711                                                   anEpsilon, bIsOnFirstBoundary);
3712
3713             }
3714           }
3715         }
3716
3717         // check if a point belong to a tangent zone. Begin
3718         for ( Standard_Integer zIt = 1; zIt <= aNbZone; zIt++ ) {
3719           gp_Pnt2d aPZone = (surfit == 0) ? aTanZoneS1->Value(zIt) : aTanZoneS2->Value(zIt);
3720           Standard_Real aZoneRadius = aTanZoneRadius->Value(zIt);
3721
3722           Standard_Integer aneighbourpointindex1 = (j == 0) ? iFirst : iLast;
3723           const IntSurf_PntOn2S& aNeighbourPoint = theWLine->Point(aneighbourpointindex1);
3724           Standard_Real nU1, nV1;
3725
3726           if(surfit == 0)
3727             aNeighbourPoint.ParametersOnS1(nU1, nV1);
3728           else
3729             aNeighbourPoint.ParametersOnS2(nU1, nV1);
3730           gp_Pnt2d ap1(nU1, nV1);
3731           gp_Pnt2d ap2 = AdjustByNeighbour( ap1, gp_Pnt2d( U, V ), aGASurface );
3732
3733
3734           if ( IsInsideTanZone( ap2, aPZone, aZoneRadius, aGASurface ) ) {
3735             aZoneIndex = zIt;
3736             bIsNearBoundary = Standard_True;
3737             if ( theReachedTol3d < aZoneRadius ) {
3738               theReachedTol3d = aZoneRadius;
3739             }
3740           }
3741         }
3742         // check if a point belong to a tangent zone. End
3743         Standard_Boolean bComputeLineEnd = Standard_False;
3744
3745         if(nbboundaries == 2) {
3746           //xf
3747           bComputeLineEnd = Standard_True;
3748           //xt
3749         }
3750         else if(nbboundaries == 1) {
3751           Standard_Boolean isperiodic = (bIsUBoundary) ? aGASurface->IsUPeriodic() : aGASurface->IsVPeriodic();
3752
3753           if(isperiodic) {
3754             Standard_Real alowerboundary = (bIsUBoundary) ? umin : vmin;
3755             Standard_Real aupperboundary = (bIsUBoundary) ? umax : vmax;
3756             Standard_Real aPeriod     = (bIsUBoundary) ? aGASurface->UPeriod() : aGASurface->VPeriod();
3757             Standard_Real aParameter = (bIsUBoundary) ? U : V;
3758             Standard_Real anoffset, anAdjustPar;
3759             GeomInt::AdjustPeriodic(aParameter, alowerboundary, aupperboundary,
3760                                            aPeriod, anAdjustPar, anoffset);
3761
3762             Standard_Real adist = (bIsFirstBoundary) ? fabs(anAdjustPar - alowerboundary) : fabs(anAdjustPar - aupperboundary);
3763             Standard_Real anotherPar = (bIsFirstBoundary) ? (aupperboundary - adist) : (alowerboundary + adist);
3764             anotherPar += anoffset;
3765             Standard_Integer aneighbourpointindex = (j == 0) ? aListOfIndex.First() : aListOfIndex.Last();
3766             const IntSurf_PntOn2S& aNeighbourPoint = theWLine->Point(aneighbourpointindex);
3767             Standard_Real nU1, nV1;
3768
3769             if(surfit == 0)
3770               aNeighbourPoint.ParametersOnS1(nU1, nV1);
3771             else
3772               aNeighbourPoint.ParametersOnS2(nU1, nV1);
3773
3774             Standard_Real adist1 = (bIsUBoundary) ? fabs(nU1 - U) : fabs(nV1 - V);
3775             Standard_Real adist2 = (bIsUBoundary) ? fabs(nU1 - anotherPar) : fabs(nV1 - anotherPar);
3776             bComputeLineEnd = Standard_True;
3777             Standard_Boolean bCheckAngle1 = Standard_False;
3778             Standard_Boolean bCheckAngle2 = Standard_False;
3779             gp_Vec2d aNewVec;
3780             Standard_Real anewU = (bIsUBoundary) ? anotherPar : U;
3781             Standard_Real anewV = (bIsUBoundary) ? V : anotherPar;
3782
3783             if(((adist1 - adist2) > Precision::PConfusion()) &&
3784                (adist2 < (aPeriod / 4.))) {
3785               bCheckAngle1 = Standard_True;
3786               aNewVec = gp_Vec2d(gp_Pnt2d(nU1, nV1), gp_Pnt2d(anewU, anewV));
3787
3788               if(aNewVec.SquareMagnitude() < gp::Resolution()) {
3789                 aNewP.SetValue((surfit == 0), anewU, anewV);
3790                 bCheckAngle1 = Standard_False;
3791               }
3792             }
3793             else if(adist1 < (aPeriod / 4.)) {
3794               bCheckAngle2 = Standard_True;
3795               aNewVec = gp_Vec2d(gp_Pnt2d(nU1, nV1), gp_Pnt2d(U, V));
3796
3797               if(aNewVec.SquareMagnitude() < gp::Resolution()) {
3798                 bCheckAngle2 = Standard_False;
3799               }
3800             }
3801
3802             if(bCheckAngle1 || bCheckAngle2) {
3803               // assume there are at least two points in line (see "if" above)
3804               Standard_Integer anindexother = aneighbourpointindex;
3805
3806               while((anindexother <= iLast) && (anindexother >= iFirst)) {
3807                 anindexother = (j == 0) ? (anindexother + 1) : (anindexother - 1);
3808                 const IntSurf_PntOn2S& aPrevNeighbourPoint = theWLine->Point(anindexother);
3809                 Standard_Real nU2, nV2;
3810
3811                 if(surfit == 0)
3812                   aPrevNeighbourPoint.ParametersOnS1(nU2, nV2);
3813                 else
3814                   aPrevNeighbourPoint.ParametersOnS2(nU2, nV2);
3815                 gp_Vec2d aVecOld(gp_Pnt2d(nU2, nV2), gp_Pnt2d(nU1, nV1));
3816
3817                 if(aVecOld.SquareMagnitude() <= gp::Resolution()) {
3818                   continue;
3819                 }
3820                 else {
3821                   Standard_Real anAngle = aNewVec.Angle(aVecOld);
3822
3823                   if((fabs(anAngle) < (M_PI * 0.25)) && (aNewVec.Dot(aVecOld) > 0.)) {
3824
3825                     if(bCheckAngle1) {
3826                       Standard_Real U1, U2, V1, V2;
3827                       IntSurf_PntOn2S atmppoint = aNewP;
3828                       atmppoint.SetValue((surfit == 0), anewU, anewV);
3829                       atmppoint.Parameters(U1, V1, U2, V2);
3830                       gp_Pnt P1 = theSurface1->Value(U1, V1);
3831                       gp_Pnt P2 = theSurface2->Value(U2, V2);
3832                       gp_Pnt P0 = aPoint.Value();
3833
3834                       if(P0.IsEqual(P1, aTol) &&
3835                          P0.IsEqual(P2, aTol) &&
3836                          P1.IsEqual(P2, aTol)) {
3837                         bComputeLineEnd = Standard_False;
3838                         aNewP.SetValue((surfit == 0), anewU, anewV);
3839                       }
3840                     }
3841
3842                     if(bCheckAngle2) {
3843                       bComputeLineEnd = Standard_False;
3844                     }
3845                   }
3846                   break;
3847                 }
3848               } // end while(anindexother...)
3849             }
3850           }
3851         }
3852         else if ( bIsNearBoundary ) {
3853           bComputeLineEnd = Standard_True;
3854         }
3855
3856         if(bComputeLineEnd) {
3857
3858           gp_Pnt2d anewpoint;
3859           Standard_Boolean found = Standard_False;
3860
3861           if ( bIsNearBoundary ) {
3862             // re-compute point near natural boundary or near tangent zone
3863             Standard_Real u1, v1, u2, v2;
3864             aNewP.Parameters( u1, v1, u2, v2 );
3865             if(surfit == 0)
3866               anewpoint = gp_Pnt2d( u1, v1 );
3867             else
3868               anewpoint = gp_Pnt2d( u2, v2 );
3869
3870             Standard_Integer aneighbourpointindex1 = (j == 0) ? iFirst : iLast;
3871             const IntSurf_PntOn2S& aNeighbourPoint = theWLine->Point(aneighbourpointindex1);
3872             Standard_Real nU1, nV1;
3873
3874             if(surfit == 0)
3875               aNeighbourPoint.ParametersOnS1(nU1, nV1);
3876             else
3877               aNeighbourPoint.ParametersOnS2(nU1, nV1);
3878             gp_Pnt2d ap1(nU1, nV1);
3879             gp_Pnt2d ap2;
3880
3881
3882             if ( aZoneIndex ) {
3883               // exclude point from a tangent zone
3884               anewpoint = AdjustByNeighbour( ap1, anewpoint, aGASurface );
3885               gp_Pnt2d aPZone = (surfit == 0) ? aTanZoneS1->Value(aZoneIndex) : aTanZoneS2->Value(aZoneIndex);
3886               Standard_Real aZoneRadius = aTanZoneRadius->Value(aZoneIndex);
3887
3888               if ( FindPoint(ap1, anewpoint, umin, umax, vmin, vmax,
3889                              aPZone, aZoneRadius, aGASurface, ap2) ) {
3890                 anewpoint = ap2;
3891                 found = Standard_True;
3892               }
3893             }
3894             else if ( aGASurface->IsUPeriodic() || aGASurface->IsVPeriodic() ) {
3895               // re-compute point near boundary if shifted on a period
3896               ap2 = AdjustByNeighbour( ap1, anewpoint, aGASurface );
3897
3898               if ( ( ap2.X() < umin ) || ( ap2.X() > umax ) ||
3899                   ( ap2.Y() < vmin ) || ( ap2.Y() > vmax ) ) {
3900                 found = FindPoint(ap1, ap2, umin, umax, vmin, vmax, anewpoint);
3901               }
3902               else {
3903                 anewpoint = ap2;
3904                 aNewP.SetValue( (surfit == 0), anewpoint.X(), anewpoint.Y() );
3905               }
3906             }
3907           }
3908           else {
3909
3910             Standard_Integer aneighbourpointindex1 = (j == 0) ? iFirst : iLast;
3911             const IntSurf_PntOn2S& aNeighbourPoint = theWLine->Point(aneighbourpointindex1);
3912             Standard_Real nU1, nV1;
3913
3914             if(surfit == 0)
3915               aNeighbourPoint.ParametersOnS1(nU1, nV1);
3916             else
3917               aNeighbourPoint.ParametersOnS2(nU1, nV1);
3918             gp_Pnt2d ap1(nU1, nV1);
3919             gp_Pnt2d ap2(nU1, nV1);
3920             Standard_Integer aneighbourpointindex2 = aneighbourpointindex1;
3921
3922             while((aneighbourpointindex2 <= iLast) && (aneighbourpointindex2 >= iFirst)) {
3923               aneighbourpointindex2 = (j == 0) ? (aneighbourpointindex2 + 1) : (aneighbourpointindex2 - 1);
3924               const IntSurf_PntOn2S& aPrevNeighbourPoint = theWLine->Point(aneighbourpointindex2);
3925               Standard_Real nU2, nV2;
3926
3927               if(surfit == 0)
3928                 aPrevNeighbourPoint.ParametersOnS1(nU2, nV2);
3929               else
3930                 aPrevNeighbourPoint.ParametersOnS2(nU2, nV2);
3931               ap2.SetX(nU2);
3932               ap2.SetY(nV2);
3933
3934               if(ap1.SquareDistance(ap2) > gp::Resolution()) {
3935                 break;
3936               }
3937             }
3938             found = FindPoint(ap2, ap1, umin, umax, vmin, vmax, anewpoint);
3939           }
3940
3941           if(found) {
3942             // check point
3943             Standard_Real aCriteria = BRep_Tool::Tolerance(theFace1) + BRep_Tool::Tolerance(theFace2);
3944             GeomAPI_ProjectPointOnSurf& aProjector =
3945               (surfit == 0) ? aContext->ProjPS(theFace2) : aContext->ProjPS(theFace1);
3946             Handle(GeomAdaptor_HSurface) aSurface = (surfit == 0) ? theSurface1 : theSurface2;
3947
3948             Handle(GeomAdaptor_HSurface) aSurfaceOther = (surfit == 0) ? theSurface2 : theSurface1;
3949
3950             gp_Pnt aP3d = aSurface->Value(anewpoint.X(), anewpoint.Y());
3951             aProjector.Perform(aP3d);
3952
3953             if(aProjector.IsDone()) {
3954               if(aProjector.LowerDistance() < aCriteria) {
3955                 Standard_Real foundU = U, foundV = V;
3956                 aProjector.LowerDistanceParameters(foundU, foundV);
3957
3958                 //Correction of projected coordinates. Begin
3959                 //Note, it may be shifted on a period
3960                 Standard_Integer aneindex1 = (j == 0) ? iFirst : iLast;
3961                 const IntSurf_PntOn2S& aNeighbourPoint = theWLine->Point(aneindex1);
3962                 Standard_Real nUn, nVn;
3963
3964                 if(surfit == 0)
3965                   aNeighbourPoint.ParametersOnS2(nUn, nVn);
3966                 else
3967                   aNeighbourPoint.ParametersOnS1(nUn, nVn);
3968                 gp_Pnt2d aNeighbour2d(nUn, nVn);
3969                 gp_Pnt2d anAdjustedPoint = AdjustByNeighbour( aNeighbour2d, gp_Pnt2d(foundU, foundV), aSurfaceOther );
3970                 foundU = anAdjustedPoint.X();
3971                 foundV = anAdjustedPoint.Y();
3972
3973                 if ( ( anAdjustedPoint.X() < umin ) && ( anAdjustedPoint.X() > umax ) &&
3974                     ( anAdjustedPoint.Y() < vmin ) && ( anAdjustedPoint.Y() > vmax ) ) {
3975                   // attempt to roughly re-compute point
3976                   foundU = ( foundU < umin ) ? umin : foundU;
3977                   foundU = ( foundU > umax ) ? umax : foundU;
3978                   foundV = ( foundV < vmin ) ? vmin : foundV;
3979                   foundV = ( foundV > vmax ) ? vmax : foundV;
3980
3981                   GeomAPI_ProjectPointOnSurf& aProjector2 =
3982                     (surfit == 0) ? aContext->ProjPS(theFace1) : aContext->ProjPS(theFace2);
3983
3984                   aP3d = aSurfaceOther->Value(foundU, foundV);
3985                   aProjector2.Perform(aP3d);
3986
3987                   if(aProjector2.IsDone()) {
3988                     if(aProjector2.LowerDistance() < aCriteria) {
3989                       Standard_Real foundU2 = anewpoint.X(), foundV2 = anewpoint.Y();
3990                       aProjector2.LowerDistanceParameters(foundU2, foundV2);
3991                       anewpoint.SetX(foundU2);
3992                       anewpoint.SetY(foundV2);
3993                     }
3994                   }
3995                 }
3996                 //Correction of projected coordinates. End
3997
3998                 if(surfit == 0)
3999                   aNewP.SetValue(aP3d, anewpoint.X(), anewpoint.Y(), foundU, foundV);
4000                 else
4001                   aNewP.SetValue(aP3d, foundU, foundV, anewpoint.X(), anewpoint.Y());
4002               }
4003             }
4004           }
4005         }
4006       }
4007       aSeqOfPntOn2S->Add(aNewP);
4008       aListOfFLIndex.Append(aSeqOfPntOn2S->NbPoints());
4009     }
4010     anArrayOfLineEnds.SetValue(i, aListOfFLIndex);
4011   }
4012   // Correct wlines.end
4013
4014   // Split wlines.begin
4015   Standard_Integer nbiter;
4016   //
4017   nbiter=1;
4018   if (!bAvoidLineConstructor) {
4019     nbiter=theLConstructor.NbParts();
4020   }
4021   //
4022   for(j = 1; j <= nbiter; ++j) {
4023     Standard_Real fprm, lprm;
4024     Standard_Integer ifprm, ilprm;
4025     //
4026     if(bAvoidLineConstructor) {
4027       ifprm = 1;
4028       ilprm = theWLine->NbPnts();
4029     }
4030     else {
4031       theLConstructor.Part(j, fprm, lprm);
4032       ifprm = (Standard_Integer)fprm;
4033       ilprm = (Standard_Integer)lprm;
4034     }
4035
4036     Handle(IntSurf_LineOn2S) aLineOn2S = new IntSurf_LineOn2S();
4037     //
4038     for(i = 1; i <= nblines; i++) {
4039       if(anArrayOfLineType.Value(i) != 0) {
4040         continue;
4041       }
4042       const TColStd_ListOfInteger& aListOfIndex = anArrayOfLines.Value(i);
4043       const TColStd_ListOfInteger& aListOfFLIndex = anArrayOfLineEnds.Value(i);
4044       Standard_Boolean bhasfirstpoint = (aListOfFLIndex.Extent() == 2);
4045       Standard_Boolean bhaslastpoint = (aListOfFLIndex.Extent() == 2);
4046
4047       if(!bhasfirstpoint && !aListOfFLIndex.IsEmpty()) {
4048         bhasfirstpoint = (i != 1);
4049       }
4050
4051       if(!bhaslastpoint && !aListOfFLIndex.IsEmpty()) {
4052         bhaslastpoint = (i != nblines);
4053       }
4054
4055       Standard_Integer iFirst = aListOfIndex.First();
4056       Standard_Integer iLast  = aListOfIndex.Last();
4057       Standard_Boolean bIsFirstInside = ((ifprm >= iFirst) && (ifprm <= iLast));
4058       Standard_Boolean bIsLastInside =  ((ilprm >= iFirst) && (ilprm <= iLast));
4059
4060       if(!bIsFirstInside && !bIsLastInside) {
4061         if((ifprm < iFirst) && (ilprm > iLast)) {
4062           // append whole line, and boundaries if neccesary
4063           if(bhasfirstpoint) {
4064             pit = aListOfFLIndex.First();
4065             const IntSurf_PntOn2S& aP = aSeqOfPntOn2S->Value(pit);
4066             aLineOn2S->Add(aP);
4067           }
4068           TColStd_ListIteratorOfListOfInteger anIt(aListOfIndex);
4069
4070           for(; anIt.More(); anIt.Next()) {
4071             pit = anIt.Value();
4072             const IntSurf_PntOn2S& aP = theWLine->Point(pit);
4073             aLineOn2S->Add(aP);
4074           }
4075
4076           if(bhaslastpoint) {
4077             pit = aListOfFLIndex.Last();
4078             const IntSurf_PntOn2S& aP = aSeqOfPntOn2S->Value(pit);
4079             aLineOn2S->Add(aP);
4080           }
4081
4082           // check end of split line (end is almost always)
4083           Standard_Integer aneighbour = i + 1;
4084           Standard_Boolean bIsEndOfLine = Standard_True;
4085
4086           if(aneighbour <= nblines) {
4087             const TColStd_ListOfInteger& aListOfNeighbourIndex = anArrayOfLines.Value(aneighbour);
4088
4089             if((anArrayOfLineType.Value(aneighbour) != 0) &&
4090                (aListOfNeighbourIndex.IsEmpty())) {
4091               bIsEndOfLine = Standard_False;
4092             }
4093           }
4094
4095           if(bIsEndOfLine) {
4096             if(aLineOn2S->NbPoints() > 1) {
4097               Handle(IntPatch_WLine) aNewWLine =
4098                 new IntPatch_WLine(aLineOn2S, Standard_False);
4099               theNewLines.Append(aNewWLine);
4100             }
4101             aLineOn2S = new IntSurf_LineOn2S();
4102           }
4103         }
4104         continue;
4105       }
4106       // end if(!bIsFirstInside && !bIsLastInside)
4107
4108       if(bIsFirstInside && bIsLastInside) {
4109         // append inside points between ifprm and ilprm
4110         TColStd_ListIteratorOfListOfInteger anIt(aListOfIndex);
4111
4112         for(; anIt.More(); anIt.Next()) {
4113           pit = anIt.Value();
4114           if((pit < ifprm) || (pit > ilprm))
4115             continue;
4116           const IntSurf_PntOn2S& aP = theWLine->Point(pit);
4117           aLineOn2S->Add(aP);
4118         }
4119       }
4120       else {
4121
4122         if(bIsFirstInside) {
4123           // append points from ifprm to last point + boundary point
4124           TColStd_ListIteratorOfListOfInteger anIt(aListOfIndex);
4125
4126           for(; anIt.More(); anIt.Next()) {
4127             pit = anIt.Value();
4128             if(pit < ifprm)
4129               continue;
4130             const IntSurf_PntOn2S& aP = theWLine->Point(pit);
4131             aLineOn2S->Add(aP);
4132           }
4133
4134           if(bhaslastpoint) {
4135             pit = aListOfFLIndex.Last();
4136             const IntSurf_PntOn2S& aP = aSeqOfPntOn2S->Value(pit);
4137             aLineOn2S->Add(aP);
4138           }
4139           // check end of split line (end is almost always)
4140           Standard_Integer aneighbour = i + 1;
4141           Standard_Boolean bIsEndOfLine = Standard_True;
4142
4143           if(aneighbour <= nblines) {
4144             const TColStd_ListOfInteger& aListOfNeighbourIndex = anArrayOfLines.Value(aneighbour);
4145
4146             if((anArrayOfLineType.Value(aneighbour) != 0) &&
4147                (aListOfNeighbourIndex.IsEmpty())) {
4148               bIsEndOfLine = Standard_False;
4149             }
4150           }
4151
4152           if(bIsEndOfLine) {
4153             if(aLineOn2S->NbPoints() > 1) {
4154               Handle(IntPatch_WLine) aNewWLine =
4155                 new IntPatch_WLine(aLineOn2S, Standard_False);
4156               theNewLines.Append(aNewWLine);
4157             }
4158             aLineOn2S = new IntSurf_LineOn2S();
4159           }
4160         }
4161         // end if(bIsFirstInside)
4162
4163         if(bIsLastInside) {
4164           // append points from first boundary point to ilprm
4165           if(bhasfirstpoint) {
4166             pit = aListOfFLIndex.First();
4167             const IntSurf_PntOn2S& aP = aSeqOfPntOn2S->Value(pit);
4168             aLineOn2S->Add(aP);
4169           }
4170           TColStd_ListIteratorOfListOfInteger anIt(aListOfIndex);
4171
4172           for(; anIt.More(); anIt.Next()) {
4173             pit = anIt.Value();
4174             if(pit > ilprm)
4175               continue;
4176             const IntSurf_PntOn2S& aP = theWLine->Point(pit);
4177             aLineOn2S->Add(aP);
4178           }
4179         }
4180         //end if(bIsLastInside)
4181       }
4182     }
4183
4184     if(aLineOn2S->NbPoints() > 1) {
4185       Handle(IntPatch_WLine) aNewWLine =
4186         new IntPatch_WLine(aLineOn2S, Standard_False);
4187       theNewLines.Append(aNewWLine);
4188     }
4189   }
4190   // Split wlines.end
4191
4192   return Standard_True;
4193 }
4194
4195 // ------------------------------------------------------------------------------------------------
4196 // static function: ParameterOutOfBoundary
4197 // purpose:         Computes a new parameter for given curve. The corresponding 2d points
4198 //                  does not lay on any boundary of given faces
4199 // ------------------------------------------------------------------------------------------------
4200 Standard_Boolean ParameterOutOfBoundary(const Standard_Real       theParameter,
4201                                         const Handle(Geom_Curve)& theCurve,
4202                                         const TopoDS_Face&        theFace1,
4203                                         const TopoDS_Face&        theFace2,
4204                                         const Standard_Real       theOtherParameter,
4205                                         const Standard_Boolean    bIncreasePar,
4206                                         Standard_Real&            theNewParameter,
4207                                         const Handle(IntTools_Context)& aContext)
4208 {
4209   Standard_Boolean bIsComputed = Standard_False;
4210   theNewParameter = theParameter;
4211
4212   Standard_Real acurpar = theParameter;
4213   TopAbs_State aState = TopAbs_ON;
4214   Standard_Integer iter = 0;
4215   Standard_Real asumtol = BRep_Tool::Tolerance(theFace1) + BRep_Tool::Tolerance(theFace2);
4216   Standard_Real adelta = asumtol * 0.1;
4217   adelta = (adelta < Precision::Confusion()) ? Precision::Confusion() : adelta;
4218   Handle(Geom_Surface) aSurf1 = BRep_Tool::Surface(theFace1);
4219   Handle(Geom_Surface) aSurf2 = BRep_Tool::Surface(theFace2);
4220
4221   Standard_Real u1, u2, v1, v2;
4222
4223   GeomAPI_ProjectPointOnSurf aPrj1;
4224   aSurf1->Bounds(u1, u2, v1, v2);
4225   aPrj1.Init(aSurf1, u1, u2, v1, v2);
4226
4227   GeomAPI_ProjectPointOnSurf aPrj2;
4228   aSurf2->Bounds(u1, u2, v1, v2);
4229   aPrj2.Init(aSurf2, u1, u2, v1, v2);
4230
4231   while(aState == TopAbs_ON) {
4232     if(bIncreasePar)
4233       acurpar += adelta;
4234     else
4235       acurpar -= adelta;
4236     gp_Pnt aPCurrent = theCurve->Value(acurpar);
4237     aPrj1.Perform(aPCurrent);
4238     Standard_Real U=0., V=0.;
4239
4240     if(aPrj1.IsDone()) {
4241       aPrj1.LowerDistanceParameters(U, V);
4242       aState = aContext->StatePointFace(theFace1, gp_Pnt2d(U, V));
4243     }
4244
4245     if(aState != TopAbs_ON) {
4246       aPrj2.Perform(aPCurrent);
4247
4248       if(aPrj2.IsDone()) {
4249         aPrj2.LowerDistanceParameters(U, V);
4250         aState = aContext->StatePointFace(theFace2, gp_Pnt2d(U, V));
4251       }
4252     }
4253
4254     if(iter > 11) {
4255       break;
4256     }
4257     iter++;
4258   }
4259
4260   if(iter <= 11) {
4261     theNewParameter = acurpar;
4262     bIsComputed = Standard_True;
4263
4264     if(bIncreasePar) {
4265       if(acurpar >= theOtherParameter)
4266         theNewParameter = theOtherParameter;
4267     }
4268     else {
4269       if(acurpar <= theOtherParameter)
4270         theNewParameter = theOtherParameter;
4271     }
4272   }
4273   return bIsComputed;
4274 }
4275
4276 //=======================================================================
4277 //function : IsCurveValid
4278 //purpose  :
4279 //=======================================================================
4280 Standard_Boolean IsCurveValid(Handle(Geom2d_Curve)& thePCurve)
4281 {
4282   if(thePCurve.IsNull())
4283     return Standard_False;
4284
4285   Standard_Real tolint = 1.e-10;
4286   Geom2dAdaptor_Curve PCA;
4287   IntRes2d_Domain PCD;
4288   Geom2dInt_GInter PCI;
4289
4290   Standard_Real pf = 0., pl = 0.;
4291   gp_Pnt2d pntf, pntl;
4292
4293   if(!thePCurve->IsClosed() && !thePCurve->IsPeriodic()) {
4294     pf = thePCurve->FirstParameter();
4295     pl = thePCurve->LastParameter();
4296     pntf = thePCurve->Value(pf);
4297     pntl = thePCurve->Value(pl);
4298     PCA.Load(thePCurve);
4299     if(!PCA.IsPeriodic()) {
4300       if(PCA.FirstParameter() > pf) pf = PCA.FirstParameter();
4301       if(PCA.LastParameter()  < pl) pl = PCA.LastParameter();
4302     }
4303     PCD.SetValues(pntf,pf,tolint,pntl,pl,tolint);
4304     PCI.Perform(PCA,PCD,tolint,tolint);
4305     if(PCI.IsDone())
4306       if(PCI.NbPoints() > 0) {
4307         return Standard_False;
4308       }
4309   }
4310
4311   return Standard_True;
4312 }
4313
4314 //=======================================================================
4315 //static function : ApproxWithPCurves
4316 //purpose  : for bug 20964 only
4317 //=======================================================================
4318 Standard_Boolean ApproxWithPCurves(const gp_Cylinder& theCyl,
4319                                    const gp_Sphere& theSph)
4320 {
4321   Standard_Boolean bRes = Standard_True;
4322   Standard_Real R1 = theCyl.Radius(), R2 = theSph.Radius();
4323   //
4324   {
4325     Standard_Real aD2, aRc2, aEps;
4326     gp_Pnt aApexSph;
4327     //
4328     aEps=1.E-7;
4329     aRc2=R1*R1;
4330     //
4331     const gp_Ax3& aAx3Sph=theSph.Position();
4332     const gp_Pnt& aLocSph=aAx3Sph.Location();
4333     const gp_Dir& aDirSph=aAx3Sph.Direction();
4334     //
4335     const gp_Ax1& aAx1Cyl=theCyl.Axis();
4336     gp_Lin aLinCyl(aAx1Cyl);
4337     //
4338     aApexSph.SetXYZ(aLocSph.XYZ()+R2*aDirSph.XYZ());
4339     aD2=aLinCyl.SquareDistance(aApexSph);
4340     if (fabs(aD2-aRc2)<aEps) {
4341       return !bRes;
4342     }
4343     //
4344     aApexSph.SetXYZ(aLocSph.XYZ()-R2*aDirSph.XYZ());
4345     aD2=aLinCyl.SquareDistance(aApexSph);
4346     if (fabs(aD2-aRc2)<aEps) {
4347       return !bRes;
4348     }
4349   }
4350   //
4351
4352   if(R1 < 2.*R2) {
4353     return bRes;
4354   }
4355   gp_Lin anCylAx(theCyl.Axis());
4356
4357   Standard_Real aDist = anCylAx.Distance(theSph.Location());
4358   Standard_Real aDRel = Abs(aDist - R1)/R2;
4359
4360   if(aDRel > .2) return bRes;
4361
4362   Standard_Real par = ElCLib::Parameter(anCylAx, theSph.Location());
4363   gp_Pnt aP = ElCLib::Value(par, anCylAx);
4364   gp_Vec aV(aP, theSph.Location());
4365
4366   Standard_Real dd = aV.Dot(theSph.Position().XDirection());
4367
4368   if(aDist < R1 && dd > 0.) return Standard_False;
4369   if(aDist > R1 && dd < 0.) return Standard_False;
4370
4371
4372   return bRes;
4373 }
4374 //=======================================================================
4375 //function : PerformPlanes
4376 //purpose  :
4377 //=======================================================================
4378 void  PerformPlanes(const Handle(GeomAdaptor_HSurface)& theS1,
4379                     const Handle(GeomAdaptor_HSurface)& theS2,
4380                     const Standard_Real TolAng,
4381                     const Standard_Real TolTang,
4382                     const Standard_Boolean theApprox1,
4383                     const Standard_Boolean theApprox2,
4384                     IntTools_SequenceOfCurves& theSeqOfCurve,
4385                     Standard_Boolean& theTangentFaces)
4386 {
4387
4388   gp_Pln aPln1 = theS1->Surface().Plane();
4389   gp_Pln aPln2 = theS2->Surface().Plane();
4390
4391   IntAna_QuadQuadGeo aPlnInter(aPln1, aPln2, TolAng, TolTang);
4392
4393   if(!aPlnInter.IsDone()) {
4394     theTangentFaces = Standard_False;
4395     return;
4396   }
4397
4398   IntAna_ResultType aResType = aPlnInter.TypeInter();
4399
4400   if(aResType == IntAna_Same) {
4401     theTangentFaces = Standard_True;
4402     return;
4403   }
4404
4405   theTangentFaces = Standard_False;
4406
4407   if(aResType == IntAna_Empty) {
4408     return;
4409   }
4410
4411   gp_Lin aLin = aPlnInter.Line(1);
4412
4413   ProjLib_Plane aProj;
4414
4415   aProj.Init(aPln1);
4416   aProj.Project(aLin);
4417   gp_Lin2d aLin2d1 = aProj.Line();
4418   //
4419   aProj.Init(aPln2);
4420   aProj.Project(aLin);
4421   gp_Lin2d aLin2d2 = aProj.Line();
4422   //
4423   //classify line2d1 relatively first plane
4424   Standard_Real P11, P12;
4425   Standard_Boolean IsCrossed = ClassifyLin2d(theS1, aLin2d1, TolTang, P11, P12);
4426   if(!IsCrossed) return;
4427   //classify line2d2 relatively second plane
4428   Standard_Real P21, P22;
4429   IsCrossed = ClassifyLin2d(theS2, aLin2d2, TolTang, P21, P22);
4430   if(!IsCrossed) return;
4431
4432   //Analysis of parametric intervals: must have common part
4433
4434   if(P21 >= P12) return;
4435   if(P22 <= P11) return;
4436
4437   Standard_Real pmin, pmax;
4438   pmin = Max(P11, P21);
4439   pmax = Min(P12, P22);
4440
4441   if(pmax - pmin <= TolTang) return;
4442
4443   Handle(Geom_Line) aGLin = new Geom_Line(aLin);
4444
4445   IntTools_Curve aCurve;
4446   Handle(Geom_TrimmedCurve) aGTLin = new Geom_TrimmedCurve(aGLin, pmin, pmax);
4447
4448   aCurve.SetCurve(aGTLin);
4449
4450   if(theApprox1) {
4451     Handle(Geom2d_Line) C2d = new Geom2d_Line(aLin2d1);
4452     aCurve.SetFirstCurve2d(new Geom2d_TrimmedCurve(C2d, pmin, pmax));
4453   }
4454   else {
4455     Handle(Geom2d_Curve) H1;
4456     aCurve.SetFirstCurve2d(H1);
4457   }
4458   if(theApprox2) {
4459     Handle(Geom2d_Line) C2d = new Geom2d_Line(aLin2d2);
4460     aCurve.SetSecondCurve2d(new Geom2d_TrimmedCurve(C2d, pmin, pmax));
4461   }
4462   else {
4463     Handle(Geom2d_Curve) H1;
4464     aCurve.SetFirstCurve2d(H1);
4465   }
4466
4467   theSeqOfCurve.Append(aCurve);
4468
4469 }
4470
4471 //=======================================================================
4472 //function : ClassifyLin2d
4473 //purpose  :
4474 //=======================================================================
4475 static inline Standard_Boolean INTER(const Standard_Real d1,
4476                                      const Standard_Real d2,
4477                                      const Standard_Real tol)
4478 {
4479   return (d1 > tol && d2 < -tol) ||
4480          (d1 < -tol && d2 > tol) ||
4481          ((d1 <= tol && d1 >= -tol) && (d2 > tol || d2 < -tol)) ||
4482          ((d2 <= tol && d2 >= -tol) && (d1 > tol || d1 < -tol));
4483 }
4484 static inline Standard_Boolean COINC(const Standard_Real d1,
4485                                      const Standard_Real d2,
4486                                      const Standard_Real tol)
4487 {
4488   return (d1 <= tol && d1 >= -tol) && (d2 <= tol && d2 >= -tol);
4489 }
4490 Standard_Boolean ClassifyLin2d(const Handle(GeomAdaptor_HSurface)& theS,
4491                                const gp_Lin2d& theLin2d,
4492                                const Standard_Real theTol,
4493                                Standard_Real& theP1,
4494                                Standard_Real& theP2)
4495
4496 {
4497   Standard_Real xmin, xmax, ymin, ymax, d1, d2, A, B, C;
4498   Standard_Real par[2];
4499   Standard_Integer nbi = 0;
4500
4501   xmin = theS->Surface().FirstUParameter();
4502   xmax = theS->Surface().LastUParameter();
4503   ymin = theS->Surface().FirstVParameter();
4504   ymax = theS->Surface().LastVParameter();
4505
4506   theLin2d.Coefficients(A, B, C);
4507
4508   //xmin, ymin <-> xmin, ymax
4509   d1 = A*xmin + B*ymin + C;
4510   d2 = A*xmin + B*ymax + C;
4511
4512   if(INTER(d1, d2, theTol)) {
4513     //Intersection with boundary
4514     Standard_Real y = -(C + A*xmin)/B;
4515     par[nbi] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmin, y));
4516     nbi++;
4517   }
4518   else if (COINC(d1, d2, theTol)) {
4519     //Coincidence with boundary
4520     par[0] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmin, ymin));
4521     par[1] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmin, ymax));
4522     nbi = 2;
4523   }
4524
4525   if(nbi == 2) {
4526
4527     if(fabs(par[0]-par[1]) > theTol) {
4528       theP1 = Min(par[0], par[1]);
4529       theP2 = Max(par[0], par[1]);
4530       return Standard_True;
4531     }
4532     else return Standard_False;
4533
4534   }
4535
4536   //xmin, ymax <-> xmax, ymax
4537   d1 = d2;
4538   d2 = A*xmax + B*ymax + C;
4539
4540   if(d1 > theTol || d1 < -theTol) {//to avoid checking of
4541                                    //coincidence with the same point
4542     if(INTER(d1, d2, theTol)) {
4543       Standard_Real x = -(C + B*ymax)/A;
4544       par[nbi] = ElCLib::Parameter(theLin2d, gp_Pnt2d(x, ymax));
4545       nbi++;
4546     }
4547     else if (COINC(d1, d2, theTol)) {
4548       par[0] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmin, ymax));
4549       par[1] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmax, ymax));
4550       nbi = 2;
4551     }
4552   }
4553
4554   if(nbi == 2) {
4555
4556     if(fabs(par[0]-par[1]) > theTol) {
4557       theP1 = Min(par[0], par[1]);
4558       theP2 = Max(par[0], par[1]);
4559       return Standard_True;
4560     }
4561     else return Standard_False;
4562
4563   }
4564
4565   //xmax, ymax <-> xmax, ymin
4566   d1 = d2;
4567   d2 = A*xmax + B*ymin + C;
4568
4569   if(d1 > theTol || d1 < -theTol) {
4570     if(INTER(d1, d2, theTol)) {
4571       Standard_Real y = -(C + A*xmax)/B;
4572       par[nbi] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmax, y));
4573       nbi++;
4574     }
4575     else if (COINC(d1, d2, theTol)) {
4576       par[0] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmax, ymax));
4577       par[1] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmax, ymin));
4578       nbi = 2;
4579     }
4580   }
4581
4582   if(nbi == 2) {
4583     if(fabs(par[0]-par[1]) > theTol) {
4584       theP1 = Min(par[0], par[1]);
4585       theP2 = Max(par[0], par[1]);
4586       return Standard_True;
4587     }
4588     else return Standard_False;
4589   }
4590
4591   //xmax, ymin <-> xmin, ymin
4592   d1 = d2;
4593   d2 = A*xmin + B*ymin + C;
4594
4595   if(d1 > theTol || d1 < -theTol) {
4596     if(INTER(d1, d2, theTol)) {
4597       Standard_Real x = -(C + B*ymin)/A;
4598       par[nbi] = ElCLib::Parameter(theLin2d, gp_Pnt2d(x, ymin));
4599       nbi++;
4600     }
4601     else if (COINC(d1, d2, theTol)) {
4602       par[0] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmax, ymin));
4603       par[1] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmin, ymin));
4604       nbi = 2;
4605     }
4606   }
4607
4608   if(nbi == 2) {
4609     if(fabs(par[0]-par[1]) > theTol) {
4610       theP1 = Min(par[0], par[1]);
4611       theP2 = Max(par[0], par[1]);
4612       return Standard_True;
4613     }
4614     else return Standard_False;
4615   }
4616
4617   return Standard_False;
4618
4619 }
4620 //
4621 //=======================================================================
4622 //function : ApproxParameters
4623 //purpose  :
4624 //=======================================================================
4625 void ApproxParameters(const Handle(GeomAdaptor_HSurface)& aHS1,
4626                       const Handle(GeomAdaptor_HSurface)& aHS2,
4627                       Standard_Integer& iDegMin,
4628                       Standard_Integer& iDegMax,
4629                       Standard_Integer& iNbIter)
4630
4631 {
4632   GeomAbs_SurfaceType aTS1, aTS2;
4633
4634   //
4635   iNbIter=0;
4636   iDegMin=4;
4637   iDegMax=8;
4638   //
4639   aTS1=aHS1->Surface().GetType();
4640   aTS2=aHS2->Surface().GetType();
4641   //
4642   // Cylinder/Torus
4643   if ((aTS1==GeomAbs_Cylinder && aTS2==GeomAbs_Torus) ||
4644       (aTS2==GeomAbs_Cylinder && aTS1==GeomAbs_Torus)) {
4645     Standard_Real aRC, aRT, dR, aPC;
4646     gp_Cylinder aCylinder;
4647     gp_Torus aTorus;
4648     //
4649     aPC=Precision::Confusion();
4650     //
4651     aCylinder=(aTS1==GeomAbs_Cylinder)? aHS1->Surface().Cylinder() : aHS2->Surface().Cylinder();
4652     aTorus=(aTS1==GeomAbs_Torus)? aHS1->Surface().Torus() : aHS2->Surface().Torus();
4653     //
4654     aRC=aCylinder.Radius();
4655     aRT=aTorus.MinorRadius();
4656     dR=aRC-aRT;
4657     if (dR<0.) {
4658       dR=-dR;
4659     }
4660     //
4661     if (dR<aPC) {
4662       iDegMax=6;
4663     }
4664   }
4665   if (aTS1==GeomAbs_Cylinder && aTS2==GeomAbs_Cylinder) {
4666     iNbIter=1;
4667   }
4668 }
4669 //=======================================================================
4670 //function : Tolerances
4671 //purpose  :
4672 //=======================================================================
4673 void Tolerances(const Handle(GeomAdaptor_HSurface)& aHS1,
4674                 const Handle(GeomAdaptor_HSurface)& aHS2,
4675                 Standard_Real& aTolTang)
4676 {
4677   GeomAbs_SurfaceType aTS1, aTS2;
4678   //
4679   aTS1=aHS1->Surface().GetType();
4680   aTS2=aHS2->Surface().GetType();
4681   //
4682   // Cylinder/Torus
4683   if ((aTS1==GeomAbs_Cylinder && aTS2==GeomAbs_Torus) ||
4684       (aTS2==GeomAbs_Cylinder && aTS1==GeomAbs_Torus)) {
4685     Standard_Real aRC, aRT, dR, aPC;
4686     gp_Cylinder aCylinder;
4687     gp_Torus aTorus;
4688     //
4689     aPC=Precision::Confusion();
4690     //
4691     aCylinder=(aTS1==GeomAbs_Cylinder)? aHS1->Surface().Cylinder() : aHS2->Surface().Cylinder();
4692     aTorus=(aTS1==GeomAbs_Torus)? aHS1->Surface().Torus() : aHS2->Surface().Torus();
4693     //
4694     aRC=aCylinder.Radius();
4695     aRT=aTorus.MinorRadius();
4696     dR=aRC-aRT;
4697     if (dR<0.) {
4698       dR=-dR;
4699     }
4700     //
4701     if (dR<aPC) {
4702       aTolTang=0.1*aTolTang;
4703     }
4704   }
4705 }
4706 //=======================================================================
4707 //function : SortTypes
4708 //purpose  :
4709 //=======================================================================
4710 Standard_Boolean SortTypes(const GeomAbs_SurfaceType aType1,
4711                            const GeomAbs_SurfaceType aType2)
4712 {
4713   Standard_Boolean bRet;
4714   Standard_Integer aI1, aI2;
4715   //
4716   bRet=Standard_False;
4717   //
4718   aI1=IndexType(aType1);
4719   aI2=IndexType(aType2);
4720   if (aI1<aI2){
4721     bRet=!bRet;
4722   }
4723   return bRet;
4724 }
4725 //=======================================================================
4726 //function : IndexType
4727 //purpose  :
4728 //=======================================================================
4729 Standard_Integer IndexType(const GeomAbs_SurfaceType aType)
4730 {
4731   Standard_Integer aIndex;
4732   //
4733   aIndex=11;
4734   //
4735   if (aType==GeomAbs_Plane) {
4736     aIndex=0;
4737   }
4738   else if (aType==GeomAbs_Cylinder) {
4739     aIndex=1;
4740   }
4741   else if (aType==GeomAbs_Cone) {
4742     aIndex=2;
4743   }
4744   else if (aType==GeomAbs_Sphere) {
4745     aIndex=3;
4746   }
4747   else if (aType==GeomAbs_Torus) {
4748     aIndex=4;
4749   }
4750   else if (aType==GeomAbs_BezierSurface) {
4751     aIndex=5;
4752   }
4753   else if (aType==GeomAbs_BSplineSurface) {
4754     aIndex=6;
4755   }
4756   else if (aType==GeomAbs_SurfaceOfRevolution) {
4757     aIndex=7;
4758   }
4759   else if (aType==GeomAbs_SurfaceOfExtrusion) {
4760     aIndex=8;
4761   }
4762   else if (aType==GeomAbs_OffsetSurface) {
4763     aIndex=9;
4764   }
4765   else if (aType==GeomAbs_OtherSurface) {
4766     aIndex=10;
4767   }
4768   return aIndex;
4769 }
4770 #ifdef OCCT_DEBUG_DUMPWLINE
4771 //=======================================================================
4772 //function : DumpWLine
4773 //purpose  :
4774 //=======================================================================
4775 void DumpWLine(const Handle(IntPatch_WLine)& aWLine)
4776 {
4777   Standard_Integer i, aNbPnts;
4778   Standard_Real aX, aY, aZ, aU1, aV1, aU2, aV2;
4779   //
4780   printf(" *WLine\n");
4781   aNbPnts=aWLine->NbPnts();
4782   for (i=1; i<=aNbPnts; ++i) {
4783     const IntSurf_PntOn2S aPntOn2S=aWLine->Point(i);
4784     const gp_Pnt& aP3D=aPntOn2S.Value();
4785     aP3D.Coord(aX, aY, aZ);
4786     aPntOn2S.Parameters(aU1, aV1, aU2, aV2);
4787     //
4788     printf("point p_%d %lf %lf %lf\n", i, aX, aY, aZ);
4789     //printf("point p_%d %20.15lf %20.15lf %20.15lf %20.15lf %20.15lf %20.15lf %20.15lf\n",
4790         //   i, aX, aY, aZ, aU1, aV1, aU2, aV2);
4791   }
4792 }
4793 #endif
4794 //=======================================================================
4795 //function : RefineVector
4796 //purpose  :
4797 //=======================================================================
4798 void RefineVector(gp_Vec2d& aV2D)
4799 {
4800   Standard_Integer k,m;
4801   Standard_Real aC[2], aEps, aR1, aR2, aNum;
4802   //
4803   aEps=RealEpsilon();
4804   aR1=1.-aEps;
4805   aR2=1.+aEps;
4806   //
4807   aV2D.Coord(aC[0], aC[1]);
4808   //
4809   for (k=0; k<2; ++k) {
4810     m=(k+1)%2;
4811     aNum=fabs(aC[k]);
4812     if (aNum>aR1 && aNum<aR2) {
4813       if (aC[k]<0.) {
4814         aC[k]=-1.;
4815       }
4816       else {
4817         aC[k]=1.;
4818       }
4819       aC[m]=0.;
4820       break;
4821     }
4822   }
4823   aV2D.SetCoord(aC[0], aC[1]);
4824 }
4825 //=======================================================================
4826 // Function : FindMaxDistance
4827 // purpose :
4828 //=======================================================================
4829 Standard_Real FindMaxDistance(const Handle(Geom_Curve)& theCurve,
4830                               const Standard_Real theFirst,
4831                               const Standard_Real theLast,
4832                               const TopoDS_Face& theFace,
4833                               const Handle(IntTools_Context)& theContext)
4834 {
4835   Standard_Integer aNbS;
4836   Standard_Real aT1, aT2, aDt, aD, aDMax, anEps;
4837   //
4838   aNbS = 11;
4839   aDt = (theLast - theFirst) / aNbS;
4840   aDMax = 0.;
4841   anEps = 1.e-4 * aDt;
4842   //
4843   GeomAPI_ProjectPointOnSurf& aProjPS = theContext->ProjPS(theFace);
4844   aT2 = theFirst;
4845   for (;;) {
4846     aT1 = aT2;
4847     aT2 += aDt;
4848     //
4849     if (aT2 > theLast) {
4850       break;
4851     }
4852     //
4853     aD = FindMaxDistance(theCurve, aT1, aT2, aProjPS, anEps);
4854     if (aD > aDMax) {
4855       aDMax = aD;
4856     }
4857   }
4858   //
4859   return aDMax;
4860 }
4861 //=======================================================================
4862 // Function : FindMaxDistance
4863 // purpose :
4864 //=======================================================================
4865 Standard_Real FindMaxDistance(const Handle(Geom_Curve)& theC,
4866                               const Standard_Real theFirst,
4867                               const Standard_Real theLast,
4868                               GeomAPI_ProjectPointOnSurf& theProjPS,
4869                               const Standard_Real theEps)
4870 {
4871   Standard_Real aA, aB, aCf, aX, aX1, aX2, aF1, aF2, aF;
4872   //
4873   aCf = 0.61803398874989484820458683436564;//(sqrt(5.)-1)/2.;
4874   aA = theFirst;
4875   aB = theLast;
4876   //
4877   aX1 = aB - aCf * (aB - aA);
4878   aF1 = MaxDistance(theC, aX1, theProjPS);
4879   aX2 = aA + aCf * (aB - aA);
4880   aF2 = MaxDistance(theC, aX2, theProjPS);
4881   //
4882   for (;;) {
4883     if ((aB - aA) < theEps) {
4884       break;
4885     }
4886     //
4887     if (aF1 > aF2) {
4888       aB = aX2;
4889       aX2 = aX1;
4890       aF2 = aF1;
4891       aX1 = aB - aCf * (aB - aA);
4892       aF1 = MaxDistance(theC, aX1, theProjPS);
4893     }
4894     else {
4895       aA = aX1;
4896       aX1 = aX2;
4897       aF1 = aF2;
4898       aX2 = aA + aCf * (aB - aA);
4899       aF2 = MaxDistance(theC, aX2, theProjPS);
4900     }
4901   }
4902   //
4903   aX = 0.5 * (aA + aB);
4904   aF = MaxDistance(theC, aX, theProjPS);
4905   //
4906   if (aF1 > aF) {
4907     aF = aF1;
4908   }
4909   //
4910   if (aF2 > aF) {
4911     aF = aF2;
4912   }
4913   //
4914   return aF;
4915 }
4916 //=======================================================================
4917 // Function : MaxDistance
4918 // purpose :
4919 //=======================================================================
4920 Standard_Real MaxDistance(const Handle(Geom_Curve)& theC,
4921                           const Standard_Real aT,
4922                           GeomAPI_ProjectPointOnSurf& theProjPS)
4923 {
4924   Standard_Real aD;
4925   gp_Pnt aP;
4926   //
4927   theC->D0(aT, aP);
4928   theProjPS.Perform(aP);
4929   aD = theProjPS.NbPoints() ? theProjPS.LowerDistance() : 0.;
4930   //
4931   return aD;
4932 }
4933
4934 //=======================================================================
4935 //function : CheckPCurve
4936 //purpose  : Checks if points of the pcurve are out of the face bounds.
4937 //=======================================================================
4938   Standard_Boolean CheckPCurve(const Handle(Geom2d_Curve)& aPC,
4939                                const TopoDS_Face& aFace)
4940 {
4941   const Standard_Integer NPoints = 23;
4942   Standard_Integer i;
4943   Standard_Real umin,umax,vmin,vmax;
4944
4945   BRepTools::UVBounds(aFace, umin, umax, vmin, vmax);
4946   Standard_Real tolU = Max ((umax-umin)*0.01, Precision::Confusion());
4947   Standard_Real tolV = Max ((vmax-vmin)*0.01, Precision::Confusion());
4948   Standard_Real fp = aPC->FirstParameter();
4949   Standard_Real lp = aPC->LastParameter();
4950
4951
4952   // adjust domain for periodic surfaces
4953   TopLoc_Location aLoc;
4954   Handle(Geom_Surface) aSurf = BRep_Tool::Surface(aFace, aLoc);
4955   if (aSurf->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface))) {
4956     aSurf = (Handle(Geom_RectangularTrimmedSurface)::DownCast(aSurf))->BasisSurface();
4957   }
4958   gp_Pnt2d pnt = aPC->Value((fp+lp)/2);
4959   Standard_Real u,v;
4960   pnt.Coord(u,v);
4961   //
4962   if (aSurf->IsUPeriodic()) {
4963     Standard_Real aPer = aSurf->UPeriod();
4964     Standard_Integer nshift = (Standard_Integer) ((u-umin)/aPer);
4965     if (u < umin+aPer*nshift) nshift--;
4966     umin += aPer*nshift;
4967     umax += aPer*nshift;
4968   }
4969   if (aSurf->IsVPeriodic()) {
4970     Standard_Real aPer = aSurf->VPeriod();
4971     Standard_Integer nshift = (Standard_Integer) ((v-vmin)/aPer);
4972     if (v < vmin+aPer*nshift) nshift--;
4973     vmin += aPer*nshift;
4974     vmax += aPer*nshift;
4975   }
4976   //
4977   //--------------------------------------------------------
4978   Standard_Boolean bRet;
4979   Standard_Integer j, aNbIntervals;
4980   Standard_Real aT, dT;
4981   gp_Pnt2d aP2D;
4982   //
4983   Geom2dAdaptor_Curve aGAC(aPC);
4984   aNbIntervals=aGAC.NbIntervals(GeomAbs_CN);
4985   //
4986   TColStd_Array1OfReal aTI(1, aNbIntervals+1);
4987   aGAC.Intervals(aTI,GeomAbs_CN);
4988   //
4989   bRet=Standard_False;
4990   //
4991   aT=aGAC.FirstParameter();
4992   for (j=1; j<=aNbIntervals; ++j) {
4993     dT=(aTI(j+1)-aTI(j))/NPoints;
4994     //
4995     for (i=1; i<NPoints; i++) {
4996       aT=aT+dT;
4997       aGAC.D0(aT, aP2D);
4998       aP2D.Coord(u,v);
4999     if (umin-u > tolU || u-umax > tolU ||
5000           vmin-v > tolV || v-vmax > tolV) {
5001         return bRet;
5002   }
5003 }
5004   }
5005   return !bRet;
5006 }
5007 //=======================================================================
5008 //function : CorrectPlaneBoundaries
5009 //purpose  :
5010 //=======================================================================
5011  void CorrectPlaneBoundaries(Standard_Real& aUmin,
5012                              Standard_Real& aUmax,
5013                              Standard_Real& aVmin,
5014                              Standard_Real& aVmax)
5015 {
5016   if (!(Precision::IsInfinite(aUmin) ||
5017         Precision::IsInfinite(aUmax))) {
5018     Standard_Real dU;
5019     //
5020     dU=0.1*(aUmax-aUmin);
5021     aUmin=aUmin-dU;
5022     aUmax=aUmax+dU;
5023   }
5024   if (!(Precision::IsInfinite(aVmin) ||
5025         Precision::IsInfinite(aVmax))) {
5026     Standard_Real dV;
5027     //
5028     dV=0.1*(aVmax-aVmin);
5029     aVmin=aVmin-dV;
5030     aVmax=aVmax+dV;
5031   }
5032 }