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 <GeomAdaptor_Curve.hxx>
  91 #include <GeomAdaptor_HSurface.hxx>
  92 #include <GeomAdaptor_Surface.hxx>
  93 #include <GeomLib_CheckBSplineCurve.hxx>
  94 #include <GeomLib_Check2dBSplineCurve.hxx>
  95
  96 #include <GeomInt_WLApprox.hxx>
  97 #include <GeomProjLib.hxx>
  98 #include <GeomAPI_ProjectPointOnSurf.hxx>
  99 #include <Geom2dAdaptor_Curve.hxx>
 100 #include <TopoDS.hxx>
 101 #include <TopoDS_Edge.hxx>
 102 #include <TopExp_Explorer.hxx>
 103
 104 #include <BRep_Tool.hxx>
 105 #include <BRepTools.hxx>
 106 #include <BRepAdaptor_Surface.hxx>
 107
 108 #include <IntTools_Curve.hxx>
 109 #include <IntTools_Tools.hxx>
 110 #include <IntTools_Tools.hxx>
 111 #include <IntTools_TopolTool.hxx>
 112 #include <IntTools_PntOnFace.hxx>
 113 #include <IntTools_PntOn2Faces.hxx>
 114 #include <IntTools_Context.hxx>
 115 #include <IntSurf_ListIteratorOfListOfPntOn2S.hxx>
 116 #include <GeomInt.hxx>
 117
 118 static
 119   void RefineVector(gp_Vec2d& aV2D);
 120 #ifdef OCCT_DEBUG_DUMPWLINE
 121 static
 122   void DumpWLine(const Handle(IntPatch_WLine)& aWLine);
 123 #endif
 124 //
 125 static
 126   void TolR3d(const TopoDS_Face& ,
 127               const TopoDS_Face& ,
 128               Standard_Real& );
 129 static
 130   Handle(Geom_Curve) MakeBSpline (const Handle(IntPatch_WLine)&,
 131                                   const Standard_Integer,
 132                                   const Standard_Integer);
 133
 134 static
 135   void Parameters(const Handle(GeomAdaptor_HSurface)&,
 136                   const Handle(GeomAdaptor_HSurface)&,
 137                   const gp_Pnt&,
 138                   Standard_Real&,
 139                   Standard_Real&,
 140                   Standard_Real&,
 141                   Standard_Real&);
 142
 143 static
 144   void BuildPCurves (Standard_Real f,Standard_Real l,Standard_Real& Tol,
 145                      const Handle (Geom_Surface)& S,
 146                      const Handle (Geom_Curve)&   C,
 147                      Handle (Geom2d_Curve)& C2d);
 148
 149 static
 150   void CorrectSurfaceBoundaries(const TopoDS_Face&  theFace,
 151                                 const Standard_Real theTolerance,
 152                                 Standard_Real&      theumin,
 153                                 Standard_Real&      theumax,
 154                                 Standard_Real&      thevmin,
 155                                 Standard_Real&      thevmax);
 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   // Line Constructor
1099   if(!bAvoidLineConstructor) {
1100     myLConstruct.Perform(L);
1101     //
1102     bDone=myLConstruct.IsDone();
1103     aNbParts=myLConstruct.NbParts();
1104     if (!bDone|| !aNbParts) {
1105       return;
1106     }
1107   }
1108   //
1109   // Do the Curve
1110   //
1111   switch (typl) {
1112   //########################################
1113   // Line, Parabola, Hyperbola
1114   //########################################
1115   case IntPatch_Lin:
1116   case IntPatch_Parabola:
1117   case IntPatch_Hyperbola: {
1118     if (typl == IntPatch_Lin) {
1119       newc =
1120         new Geom_Line (Handle(IntPatch_GLine)::DownCast(L)->Line());
1121     }
1122
1123     else if (typl == IntPatch_Parabola) {
1124       newc =
1125         new Geom_Parabola(Handle(IntPatch_GLine)::DownCast(L)->Parabola());
1126     }
1127
1128     else if (typl == IntPatch_Hyperbola) {
1129       newc =
1130         new Geom_Hyperbola (Handle(IntPatch_GLine)::DownCast(L)->Hyperbola());
1131     }
1132     //
1133     // myTolReached3d
1134     if (typl == IntPatch_Lin) {
1135       TolR3d (myFace1, myFace2, myTolReached3d);
1136     }
1137     //
1138     aNbParts=myLConstruct.NbParts();
1139     for (i=1; i<=aNbParts; i++) {
1140       Standard_Boolean bFNIt, bLPIt;
1141       //
1142       myLConstruct.Part(i, fprm, lprm);
1143       //
1144       bFNIt=Precision::IsNegativeInfinite(fprm);
1145       bLPIt=Precision::IsPositiveInfinite(lprm);
1146       //
1147       if (!bFNIt && !bLPIt) {
1148         //
1149         IntTools_Curve aCurve;
1150         //
1151         Handle(Geom_TrimmedCurve) aCT3D=new Geom_TrimmedCurve(newc, fprm, lprm);
1152         aCurve.SetCurve(aCT3D);
1153         if (typl == IntPatch_Parabola) {
1154           Standard_Real aTolF1, aTolF2, aTolBase;
1155
1156           aTolF1 = BRep_Tool::Tolerance(myFace1);
1157           aTolF2 = BRep_Tool::Tolerance(myFace2);
1158           aTolBase=aTolF1+aTolF2;
1159           myTolReached3d=IntTools_Tools::CurveTolerance(aCT3D, aTolBase);
1160         }
1161         //
1162         aCurve.SetCurve(new Geom_TrimmedCurve(newc, fprm, lprm));
1163         if(myApprox1) {
1164           Handle (Geom2d_Curve) C2d;
1165           BuildPCurves(fprm, lprm, Tolpc, myHS1->ChangeSurface().Surface(), newc, C2d);
1166           if(Tolpc>myTolReached2d || myTolReached2d==0.) {
1167             myTolReached2d=Tolpc;
1168           }
1169             //
1170             aCurve.SetFirstCurve2d(new Geom2d_TrimmedCurve(C2d,fprm,lprm));
1171           }
1172           else {
1173             Handle(Geom2d_BSplineCurve) H1;
1174             //
1175             aCurve.SetFirstCurve2d(H1);
1176           }
1177         //
1178         if(myApprox2) {
1179           Handle (Geom2d_Curve) C2d;
1180           BuildPCurves(fprm,lprm,Tolpc,myHS2->ChangeSurface().Surface(),newc,C2d);
1181           if(Tolpc>myTolReached2d || myTolReached2d==0.) {
1182             myTolReached2d=Tolpc;
1183           }
1184           //
1185           aCurve.SetSecondCurve2d(new Geom2d_TrimmedCurve(C2d,fprm,lprm));
1186           }
1187         else {
1188           Handle(Geom2d_BSplineCurve) H1;
1189           //
1190           aCurve.SetSecondCurve2d(H1);
1191         }
1192         mySeqOfCurve.Append(aCurve);
1193       } //if (!bFNIt && !bLPIt) {
1194       else {
1195         //  on regarde si on garde
1196         //
1197         Standard_Real aTestPrm, dT=100.;
1198         //
1199         aTestPrm=0.;
1200         if (bFNIt && !bLPIt) {
1201           aTestPrm=lprm-dT;
1202         }
1203         else if (!bFNIt && bLPIt) {
1204           aTestPrm=fprm+dT;
1205         }
1206         else {
1207           // i.e, if (bFNIt && bLPIt)
1208           aTestPrm=IntTools_Tools::IntermediatePoint(-dT, dT);
1209         }
1210         //
1211         gp_Pnt ptref(newc->Value(aTestPrm));
1212         //
1213         GeomAbs_SurfaceType typS1 = myHS1->GetType();
1214         GeomAbs_SurfaceType typS2 = myHS2->GetType();
1215         if( typS1 == GeomAbs_SurfaceOfExtrusion ||
1216             typS1 == GeomAbs_OffsetSurface ||
1217             typS1 == GeomAbs_SurfaceOfRevolution ||
1218             typS2 == GeomAbs_SurfaceOfExtrusion ||
1219             typS2 == GeomAbs_OffsetSurface ||
1220             typS2 == GeomAbs_SurfaceOfRevolution) {
1221           Handle(Geom2d_BSplineCurve) H1;
1222           mySeqOfCurve.Append(IntTools_Curve(newc, H1, H1));
1223           continue;
1224         }
1225
1226         Standard_Real u1, v1, u2, v2, Tol;
1227
1228         Tol = Precision::Confusion();
1229         Parameters(myHS1, myHS2, ptref,  u1, v1, u2, v2);
1230         ok = (dom1->Classify(gp_Pnt2d(u1, v1), Tol) != TopAbs_OUT);
1231         if(ok) {
1232           ok = (dom2->Classify(gp_Pnt2d(u2,v2),Tol) != TopAbs_OUT);
1233         }
1234         if (ok) {
1235           Handle(Geom2d_BSplineCurve) H1;
1236           mySeqOfCurve.Append(IntTools_Curve(newc, H1, H1));
1237         }
1238       }
1239     }// for (i=1; i<=aNbParts; i++) {
1240   }// case IntPatch_Lin:  case IntPatch_Parabola:  case IntPatch_Hyperbola:
1241     break;
1242
1243   //########################################
1244   // Circle and Ellipse
1245   //########################################
1246   case IntPatch_Circle:
1247   case IntPatch_Ellipse: {
1248
1249     if (typl == IntPatch_Circle) {
1250       newc = new Geom_Circle
1251         (Handle(IntPatch_GLine)::DownCast(L)->Circle());
1252     }
1253     else { //IntPatch_Ellipse
1254       newc = new Geom_Ellipse
1255         (Handle(IntPatch_GLine)::DownCast(L)->Ellipse());
1256     }
1257     //
1258     // myTolReached3d
1259     TolR3d (myFace1, myFace2, myTolReached3d);
1260     //
1261     aNbParts=myLConstruct.NbParts();
1262     //
1263     Standard_Real aPeriod, aNul;
1264     TColStd_SequenceOfReal aSeqFprm,  aSeqLprm;
1265
1266     aNul=0.;
1267     aPeriod=M_PI+M_PI;
1268
1269     for (i=1; i<=aNbParts; i++) {
1270       myLConstruct.Part(i, fprm, lprm);
1271
1272       if (fprm < aNul && lprm > aNul) {
1273         // interval that goes through 0. is divided on two intervals;
1274         while (fprm<aNul || fprm>aPeriod) fprm=fprm+aPeriod;
1275         while (lprm<aNul || lprm>aPeriod) lprm=lprm+aPeriod;
1276         //
1277         if((aPeriod - fprm) > Tolpc) {
1278           aSeqFprm.Append(fprm);
1279           aSeqLprm.Append(aPeriod);
1280         }
1281         else {
1282           gp_Pnt P1 = newc->Value(fprm);
1283           gp_Pnt P2 = newc->Value(aPeriod);
1284           Standard_Real aTolDist = BRep_Tool::Tolerance(myFace1) + BRep_Tool::Tolerance(myFace2);
1285           aTolDist = (myTolReached3d > aTolDist) ? myTolReached3d : aTolDist;
1286
1287           if(P1.Distance(P2) > aTolDist) {
1288             Standard_Real anewpar = fprm;
1289
1290             if(ParameterOutOfBoundary(fprm, newc, myFace1, myFace2,
1291                                       lprm, Standard_False, anewpar, myContext)) {
1292               fprm = anewpar;
1293             }
1294             aSeqFprm.Append(fprm);
1295             aSeqLprm.Append(aPeriod);
1296           }
1297         }
1298
1299         //
1300         if((lprm - aNul) > Tolpc) {
1301           aSeqFprm.Append(aNul);
1302           aSeqLprm.Append(lprm);
1303         }
1304         else {
1305           gp_Pnt P1 = newc->Value(aNul);
1306           gp_Pnt P2 = newc->Value(lprm);
1307           Standard_Real aTolDist = BRep_Tool::Tolerance(myFace1) + BRep_Tool::Tolerance(myFace2);
1308           aTolDist = (myTolReached3d > aTolDist) ? myTolReached3d : aTolDist;
1309
1310           if(P1.Distance(P2) > aTolDist) {
1311             Standard_Real anewpar = lprm;
1312
1313             if(ParameterOutOfBoundary(lprm, newc, myFace1, myFace2,
1314                                       fprm, Standard_True, anewpar, myContext)) {
1315               lprm = anewpar;
1316             }
1317             aSeqFprm.Append(aNul);
1318             aSeqLprm.Append(lprm);
1319           }
1320         }
1321       }
1322       else {
1323         // usual interval
1324         aSeqFprm.Append(fprm);
1325         aSeqLprm.Append(lprm);
1326       }
1327     }
1328
1329     //
1330     aNbParts=aSeqFprm.Length();
1331     for (i=1; i<=aNbParts; i++) {
1332       fprm=aSeqFprm(i);
1333       lprm=aSeqLprm(i);
1334       //
1335       Standard_Real aRealEpsilon=RealEpsilon();
1336       if (Abs(fprm) > aRealEpsilon || Abs(lprm-2.*M_PI) > aRealEpsilon) {
1337         //==============================================
1338         ////
1339         IntTools_Curve aCurve;
1340         Handle(Geom_TrimmedCurve) aTC3D=new Geom_TrimmedCurve(newc,fprm,lprm);
1341         aCurve.SetCurve(aTC3D);
1342         fprm=aTC3D->FirstParameter();
1343         lprm=aTC3D->LastParameter ();
1344         ////
1345         if (typl == IntPatch_Circle || typl == IntPatch_Ellipse) {////
1346           if(myApprox1) {
1347             Handle (Geom2d_Curve) C2d;
1348             BuildPCurves(fprm,lprm,Tolpc,myHS1->ChangeSurface().Surface(),newc,C2d);
1349             if(Tolpc>myTolReached2d || myTolReached2d==0) {
1350               myTolReached2d=Tolpc;
1351             }
1352             //
1353             aCurve.SetFirstCurve2d(C2d);
1354           }
1355           else { ////
1356             Handle(Geom2d_BSplineCurve) H1;
1357             aCurve.SetFirstCurve2d(H1);
1358           }
1359
1360
1361           if(myApprox2) {
1362             Handle (Geom2d_Curve) C2d;
1363             BuildPCurves(fprm,lprm,Tolpc,myHS2->ChangeSurface().Surface(),newc,C2d);
1364             if(Tolpc>myTolReached2d || myTolReached2d==0) {
1365               myTolReached2d=Tolpc;
1366             }
1367             //
1368             aCurve.SetSecondCurve2d(C2d);
1369           }
1370           else {
1371             Handle(Geom2d_BSplineCurve) H1;
1372             aCurve.SetSecondCurve2d(H1);
1373           }
1374         }
1375
1376         else {
1377           Handle(Geom2d_BSplineCurve) H1;
1378           aCurve.SetFirstCurve2d(H1);
1379           aCurve.SetSecondCurve2d(H1);
1380         }
1381         mySeqOfCurve.Append(aCurve);
1382           //==============================================
1383       } //if (Abs(fprm) > RealEpsilon() || Abs(lprm-2.*M_PI) > RealEpsilon())
1384
1385       else {
1386         //  on regarde si on garde
1387         //
1388         if (aNbParts==1) {
1389 //           if (Abs(fprm) < RealEpsilon() &&  Abs(lprm-2.*M_PI) < RealEpsilon()) {
1390           if (Abs(fprm) <= aRealEpsilon && Abs(lprm-2.*M_PI) <= aRealEpsilon) {
1391             IntTools_Curve aCurve;
1392             Handle(Geom_TrimmedCurve) aTC3D=new Geom_TrimmedCurve(newc,fprm,lprm);
1393             aCurve.SetCurve(aTC3D);
1394             fprm=aTC3D->FirstParameter();
1395             lprm=aTC3D->LastParameter ();
1396
1397             if(myApprox1) {
1398               Handle (Geom2d_Curve) C2d;
1399               BuildPCurves(fprm,lprm,Tolpc,myHS1->ChangeSurface().Surface(),newc,C2d);
1400               if(Tolpc>myTolReached2d || myTolReached2d==0) {
1401                 myTolReached2d=Tolpc;
1402               }
1403               //
1404               aCurve.SetFirstCurve2d(C2d);
1405             }
1406             else { ////
1407               Handle(Geom2d_BSplineCurve) H1;
1408               aCurve.SetFirstCurve2d(H1);
1409             }
1410
1411             if(myApprox2) {
1412               Handle (Geom2d_Curve) C2d;
1413               BuildPCurves(fprm,lprm,Tolpc,myHS2->ChangeSurface().Surface(),newc,C2d);
1414               if(Tolpc>myTolReached2d || myTolReached2d==0) {
1415                 myTolReached2d=Tolpc;
1416               }
1417               //
1418               aCurve.SetSecondCurve2d(C2d);
1419             }
1420             else {
1421               Handle(Geom2d_BSplineCurve) H1;
1422               aCurve.SetSecondCurve2d(H1);
1423             }
1424             mySeqOfCurve.Append(aCurve);
1425             break;
1426           }
1427         }
1428         //
1429         Standard_Real aTwoPIdiv17, u1, v1, u2, v2, Tol;
1430
1431         aTwoPIdiv17=2.*M_PI/17.;
1432
1433         for (j=0; j<=17; j++) {
1434           gp_Pnt ptref (newc->Value (j*aTwoPIdiv17));
1435           Tol = Precision::Confusion();
1436
1437           Parameters(myHS1, myHS2, ptref, u1, v1, u2, v2);
1438           ok = (dom1->Classify(gp_Pnt2d(u1,v1),Tol) != TopAbs_OUT);
1439           if(ok) {
1440             ok = (dom2->Classify(gp_Pnt2d(u2,v2),Tol) != TopAbs_OUT);
1441           }
1442           if (ok) {
1443             IntTools_Curve aCurve;
1444             aCurve.SetCurve(newc);
1445             //==============================================
1446             if (typl == IntPatch_Circle || typl == IntPatch_Ellipse) {
1447
1448               if(myApprox1) {
1449                 Handle (Geom2d_Curve) C2d;
1450                 BuildPCurves(fprm, lprm, Tolpc, myHS1->ChangeSurface().Surface(), newc, C2d);
1451                 if(Tolpc>myTolReached2d || myTolReached2d==0) {
1452                   myTolReached2d=Tolpc;
1453                 }
1454                 //
1455                 aCurve.SetFirstCurve2d(C2d);
1456               }
1457               else {
1458                 Handle(Geom2d_BSplineCurve) H1;
1459                 aCurve.SetFirstCurve2d(H1);
1460               }
1461
1462               if(myApprox2) {
1463                 Handle (Geom2d_Curve) C2d;
1464                 BuildPCurves(fprm, lprm, Tolpc,myHS2->ChangeSurface().Surface(), newc, C2d);
1465                 if(Tolpc>myTolReached2d || myTolReached2d==0) {
1466                   myTolReached2d=Tolpc;
1467                 }
1468                 //
1469                 aCurve.SetSecondCurve2d(C2d);
1470               }
1471
1472               else {
1473                 Handle(Geom2d_BSplineCurve) H1;
1474                 aCurve.SetSecondCurve2d(H1);
1475               }
1476             }//  end of if (typl == IntPatch_Circle || typl == IntPatch_Ellipse)
1477
1478             else {
1479               Handle(Geom2d_BSplineCurve) H1;
1480               //
1481               aCurve.SetFirstCurve2d(H1);
1482               aCurve.SetSecondCurve2d(H1);
1483             }
1484             //==============================================
1485             //
1486             mySeqOfCurve.Append(aCurve);
1487             break;
1488
1489             }//  end of if (ok) {
1490           }//  end of for (Standard_Integer j=0; j<=17; j++)
1491         }//  end of else { on regarde si on garde
1492       }// for (i=1; i<=myLConstruct.NbParts(); i++)
1493     }// IntPatch_Circle: IntPatch_Ellipse:
1494     break;
1495
1496   case IntPatch_Analytic: {
1497     IntSurf_Quadric quad1,quad2;
1498     GeomAbs_SurfaceType typs = myHS1->Surface().GetType();
1499
1500     switch (typs) {
1501       case GeomAbs_Plane:
1502         quad1.SetValue(myHS1->Surface().Plane());
1503         break;
1504       case GeomAbs_Cylinder:
1505         quad1.SetValue(myHS1->Surface().Cylinder());
1506         break;
1507       case GeomAbs_Cone:
1508         quad1.SetValue(myHS1->Surface().Cone());
1509         break;
1510       case GeomAbs_Sphere:
1511         quad1.SetValue(myHS1->Surface().Sphere());
1512         break;
1513     case GeomAbs_Torus:
1514       quad1.SetValue(myHS1->Surface().Torus());
1515       break;
1516       default:
1517         Standard_ConstructionError::Raise("GeomInt_IntSS::MakeCurve 1");
1518       }
1519
1520     typs = myHS2->Surface().GetType();
1521
1522     switch (typs) {
1523       case GeomAbs_Plane:
1524         quad2.SetValue(myHS2->Surface().Plane());
1525         break;
1526       case GeomAbs_Cylinder:
1527         quad2.SetValue(myHS2->Surface().Cylinder());
1528         break;
1529       case GeomAbs_Cone:
1530         quad2.SetValue(myHS2->Surface().Cone());
1531         break;
1532       case GeomAbs_Sphere:
1533         quad2.SetValue(myHS2->Surface().Sphere());
1534         break;
1535     case GeomAbs_Torus:
1536       quad2.SetValue(myHS2->Surface().Torus());
1537       break;
1538       default:
1539         Standard_ConstructionError::Raise("GeomInt_IntSS::MakeCurve 2");
1540       }
1541     //
1542     //=========
1543     IntPatch_ALineToWLine convert (quad1, quad2);
1544
1545     if (!myApprox) {
1546       aNbParts=myLConstruct.NbParts();
1547       for (i=1; i<=aNbParts; i++) {
1548         myLConstruct.Part(i, fprm, lprm);
1549         Handle(IntPatch_WLine) WL =
1550           convert.MakeWLine(Handle(IntPatch_ALine)::DownCast(L), fprm, lprm);
1551         //
1552         Handle(Geom2d_BSplineCurve) H1;
1553         Handle(Geom2d_BSplineCurve) H2;
1554
1555         if(myApprox1) {
1556           H1 = MakeBSpline2d(WL, 1, WL->NbPnts(), Standard_True);
1557         }
1558
1559         if(myApprox2) {
1560           H2 = MakeBSpline2d(WL, 1, WL->NbPnts(), Standard_False);
1561         }
1562         //
1563         mySeqOfCurve.Append(IntTools_Curve(MakeBSpline(WL,1,WL->NbPnts()), H1, H2));
1564       }
1565     } // if (!myApprox)
1566
1567     else { // myApprox=TRUE
1568       GeomInt_WLApprox theapp3d;
1569       //
1570       Standard_Real tol2d = myTolApprox;
1571       //
1572       theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, Standard_True);
1573
1574       aNbParts=myLConstruct.NbParts();
1575       for (i=1; i<=aNbParts; i++) {
1576         myLConstruct.Part(i, fprm, lprm);
1577         Handle(IntPatch_WLine) WL =
1578           convert.MakeWLine(Handle(IntPatch_ALine):: DownCast(L),fprm,lprm);
1579
1580         theapp3d.Perform(myHS1,myHS2,WL,Standard_True,myApprox1,myApprox2, 1, WL->NbPnts());
1581
1582         if (!theapp3d.IsDone()) {
1583           //
1584           Handle(Geom2d_BSplineCurve) H1;
1585           Handle(Geom2d_BSplineCurve) H2;
1586
1587           if(myApprox1) {
1588             H1 = MakeBSpline2d(WL, 1, WL->NbPnts(), Standard_True);
1589           }
1590
1591           if(myApprox2) {
1592             H2 = MakeBSpline2d(WL, 1, WL->NbPnts(), Standard_False);
1593           }
1594           //
1595           mySeqOfCurve.Append(IntTools_Curve(MakeBSpline(WL,1,WL->NbPnts()), H1, H2));
1596         }
1597
1598         else {
1599           if(myApprox1 || myApprox2) {
1600             if( theapp3d.TolReached2d()>myTolReached2d || myTolReached2d==0) {
1601               myTolReached2d = theapp3d.TolReached2d();
1602             }
1603           }
1604
1605           if( theapp3d.TolReached3d()>myTolReached3d || myTolReached3d==0) {
1606             myTolReached3d = theapp3d.TolReached3d();
1607           }
1608
1609           Standard_Integer aNbMultiCurves, nbpoles;
1610           aNbMultiCurves=theapp3d.NbMultiCurves();
1611           for (j=1; j<=aNbMultiCurves; j++) {
1612             const AppParCurves_MultiBSpCurve& mbspc = theapp3d.Value(j);
1613             nbpoles = mbspc.NbPoles();
1614
1615             TColgp_Array1OfPnt tpoles(1, nbpoles);
1616             mbspc.Curve(1, tpoles);
1617             Handle(Geom_BSplineCurve) BS=new Geom_BSplineCurve(tpoles,
1618                                                                mbspc.Knots(),
1619                                                                mbspc.Multiplicities(),
1620                                                                mbspc.Degree());
1621
1622             GeomLib_CheckBSplineCurve Check(BS,TOLCHECK,TOLANGCHECK);
1623             Check.FixTangent(Standard_True,Standard_True);
1624             //
1625             IntTools_Curve aCurve;
1626             aCurve.SetCurve(BS);
1627
1628             if(myApprox1) {
1629               TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
1630               mbspc.Curve(2,tpoles2d);
1631               Handle(Geom2d_BSplineCurve) BS2=new Geom2d_BSplineCurve(tpoles2d,
1632                                                                       mbspc.Knots(),
1633                                                                       mbspc.Multiplicities(),
1634                                                                       mbspc.Degree());
1635
1636               GeomLib_Check2dBSplineCurve newCheck(BS2,TOLCHECK,TOLANGCHECK);
1637               newCheck.FixTangent(Standard_True,Standard_True);
1638               //
1639               aCurve.SetFirstCurve2d(BS2);
1640             }
1641             else {
1642               Handle(Geom2d_BSplineCurve) H1;
1643               aCurve.SetFirstCurve2d(H1);
1644             }
1645
1646             if(myApprox2) {
1647               TColgp_Array1OfPnt2d tpoles2d(1, nbpoles);
1648               Standard_Integer TwoOrThree;
1649               TwoOrThree=myApprox1 ? 3 : 2;
1650               mbspc.Curve(TwoOrThree, 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.SetSecondCurve2d(BS2);
1660             }
1661             else {
1662               Handle(Geom2d_BSplineCurve) H2;
1663               aCurve.SetSecondCurve2d(H2);
1664             }
1665             //
1666             mySeqOfCurve.Append(aCurve);
1667
1668           }// for (j=1; j<=aNbMultiCurves; j++) {
1669         }// else from if (!theapp3d.IsDone())
1670       }// for (i=1; i<=aNbParts; i++) {
1671     }// else { // myApprox=TRUE
1672   }// case IntPatch_Analytic:
1673     break;
1674
1675   case IntPatch_Walking:{
1676     Handle(IntPatch_WLine) WL =
1677       Handle(IntPatch_WLine)::DownCast(L);
1678     //
1679     Standard_Integer ifprm, ilprm;
1680     //
1681     if (!myApprox) {
1682       aNbParts = 1;
1683       if(!bAvoidLineConstructor){
1684         aNbParts=myLConstruct.NbParts();
1685       }
1686       for (i=1; i<=aNbParts; ++i) {
1687         Handle(Geom2d_BSplineCurve) H1, H2;
1688         Handle(Geom_Curve) aBSp;
1689         //
1690         if(bAvoidLineConstructor) {
1691           ifprm = 1;
1692           ilprm = WL->NbPnts();
1693         }
1694         else {
1695           myLConstruct.Part(i, fprm, lprm);
1696           ifprm=(Standard_Integer)fprm;
1697           ilprm=(Standard_Integer)lprm;
1698         }
1699         //
1700         if(myApprox1) {
1701           H1 = MakeBSpline2d(WL, ifprm, ilprm, Standard_True);
1702         }
1703         //
1704         if(myApprox2) {
1705           H2 = MakeBSpline2d(WL, ifprm, ilprm, Standard_False);
1706         }
1707         //
1708         aBSp=MakeBSpline(WL, ifprm, ilprm);
1709         IntTools_Curve aIC(aBSp, H1, H2);
1710         mySeqOfCurve.Append(aIC);
1711       }// for (i=1; i<=aNbParts; ++i) {
1712     }// if (!myApprox) {
1713     //
1714     else { // X
1715       Standard_Boolean bIsDecomposited;
1716       Standard_Integer nbiter, aNbSeqOfL;
1717       Standard_Real tol2d, aTolApproxImp;
1718       IntPatch_SequenceOfLine aSeqOfL;
1719       GeomInt_WLApprox theapp3d;
1720       Approx_ParametrizationType aParType = Approx_ChordLength;
1721       //
1722       Standard_Boolean anApprox1 = myApprox1;
1723       Standard_Boolean anApprox2 = myApprox2;
1724       //
1725       aTolApproxImp=1.e-5;
1726       tol2d = myTolApprox;
1727
1728       GeomAbs_SurfaceType typs1, typs2;
1729       typs1 = myHS1->Surface().GetType();
1730       typs2 = myHS2->Surface().GetType();
1731       Standard_Boolean anWithPC = Standard_True;
1732
1733       if(typs1 == GeomAbs_Cylinder && typs2 == GeomAbs_Sphere) {
1734         anWithPC =
1735           ApproxWithPCurves(myHS1->Surface().Cylinder(), myHS2->Surface().Sphere());
1736       }
1737       else if (typs1 == GeomAbs_Sphere && typs2 == GeomAbs_Cylinder) {
1738         anWithPC =
1739           ApproxWithPCurves(myHS2->Surface().Cylinder(), myHS1->Surface().Sphere());
1740       }
1741       //
1742       if(!anWithPC) {
1743         myTolApprox = aTolApproxImp;//1.e-5;
1744         anApprox1 = Standard_False;
1745         anApprox2 = Standard_False;
1746         //
1747         tol2d = myTolApprox;
1748       }
1749
1750       if(myHS1 == myHS2) {
1751         theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, Standard_False, aParType);
1752         rejectSurface = Standard_True;
1753       }
1754       else {
1755         if(reApprox && !rejectSurface)
1756           theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, Standard_False, aParType);
1757         else {
1758           Standard_Integer iDegMax, iDegMin, iNbIter;
1759           //
1760           ApproxParameters(myHS1, myHS2, iDegMin, iDegMax, iNbIter);
1761           theapp3d.SetParameters(myTolApprox, tol2d, iDegMin, iDegMax, iNbIter, Standard_True, aParType);
1762         }
1763       }
1764       //
1765       Standard_Real aReachedTol = Precision::Confusion();
1766       bIsDecomposited=DecompositionOfWLine(WL,
1767                                            myHS1,
1768                                            myHS2,
1769                                            myFace1,
1770                                            myFace2,
1771                                            myLConstruct,
1772                                            bAvoidLineConstructor,
1773                                            aSeqOfL,
1774                                            aReachedTol,
1775                                            myContext);
1776       if ( bIsDecomposited && ( myTolReached3d < aReachedTol ) ) {
1777         myTolReached3d = aReachedTol;
1778       }
1779       //
1780       aNbSeqOfL=aSeqOfL.Length();
1781       //
1782       if (bIsDecomposited) {
1783         nbiter=aNbSeqOfL;
1784       }
1785       else {
1786         nbiter=1;
1787         aNbParts=1;
1788         if (!bAvoidLineConstructor) {
1789           aNbParts=myLConstruct.NbParts();
1790           nbiter=aNbParts;
1791         }
1792       }
1793       //
1794       for(i = 1; i <= nbiter; ++i) {
1795         if(bIsDecomposited) {
1796           WL = Handle(IntPatch_WLine)::DownCast(aSeqOfL.Value(i));
1797           ifprm = 1;
1798           ilprm = WL->NbPnts();
1799         }
1800         else {
1801           if(bAvoidLineConstructor) {
1802             ifprm = 1;
1803             ilprm = WL->NbPnts();
1804           }
1805           else {
1806             myLConstruct.Part(i, fprm, lprm);
1807             ifprm = (Standard_Integer)fprm;
1808             ilprm = (Standard_Integer)lprm;
1809           }
1810         }
1811         //-- lbr :
1812         //-- Si une des surfaces est un plan , on approxime en 2d
1813         //-- sur cette surface et on remonte les points 2d en 3d.
1814         if(typs1 == GeomAbs_Plane) {
1815           theapp3d.Perform(myHS1, myHS2, WL, Standard_False,Standard_True, myApprox2,ifprm,ilprm);
1816         }
1817         else if(typs2 == GeomAbs_Plane) {
1818           theapp3d.Perform(myHS1,myHS2,WL,Standard_False,myApprox1,Standard_True,ifprm,ilprm);
1819         }
1820         else {
1821           //
1822           if (myHS1 != myHS2){
1823             if ((typs1==GeomAbs_BezierSurface || typs1==GeomAbs_BSplineSurface) &&
1824                 (typs2==GeomAbs_BezierSurface || typs2==GeomAbs_BSplineSurface)) {
1825
1826               theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, Standard_True, aParType);
1827
1828               Standard_Boolean bUseSurfaces;
1829               bUseSurfaces=NotUseSurfacesForApprox(myFace1, myFace2, WL, ifprm,  ilprm);
1830               if (bUseSurfaces) {
1831                 // ######
1832                 rejectSurface = Standard_True;
1833                 // ######
1834                 theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, Standard_False, aParType);
1835               }
1836             }
1837           }
1838           //
1839           theapp3d.Perform(myHS1,myHS2,WL,Standard_True,anApprox1,anApprox2,ifprm,ilprm);
1840         }
1841           //
1842         if (!theapp3d.IsDone()) {
1843           Handle(Geom2d_BSplineCurve) H1;
1844           Handle(Geom2d_BSplineCurve) H2;
1845           //
1846           Handle(Geom_Curve) aBSp=MakeBSpline(WL,ifprm, ilprm);
1847           //
1848           if(myApprox1) {
1849             H1 = MakeBSpline2d(WL, ifprm, ilprm, Standard_True);
1850           }
1851           //
1852           if(myApprox2) {
1853             H2 = MakeBSpline2d(WL, ifprm, ilprm, Standard_False);
1854           }
1855           //
1856           IntTools_Curve aIC(aBSp, H1, H2);
1857           mySeqOfCurve.Append(aIC);
1858         }
1859
1860         else {
1861           if(myApprox1 || myApprox2 || (typs1==GeomAbs_Plane || typs2==GeomAbs_Plane)) {
1862             if( theapp3d.TolReached2d()>myTolReached2d || myTolReached2d==0.) {
1863               myTolReached2d = theapp3d.TolReached2d();
1864             }
1865           }
1866           if(typs1==GeomAbs_Plane || typs2==GeomAbs_Plane) {
1867             myTolReached3d = myTolReached2d;
1868             //
1869             if (typs1==GeomAbs_Torus || typs2==GeomAbs_Torus) {
1870               if (myTolReached3d<1.e-6) {
1871                 myTolReached3d = theapp3d.TolReached3d();
1872                 myTolReached3d=1.e-6;
1873               }
1874             }
1875           }
1876           else  if( theapp3d.TolReached3d()>myTolReached3d || myTolReached3d==0.) {
1877             myTolReached3d = theapp3d.TolReached3d();
1878           }
1879
1880           Standard_Integer aNbMultiCurves, nbpoles;
1881           aNbMultiCurves=theapp3d.NbMultiCurves();
1882           for (j=1; j<=aNbMultiCurves; j++) {
1883             if(typs1 == GeomAbs_Plane) {
1884               const AppParCurves_MultiBSpCurve& mbspc = theapp3d.Value(j);
1885               nbpoles = mbspc.NbPoles();
1886
1887               TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
1888               TColgp_Array1OfPnt   tpoles(1,nbpoles);
1889
1890               mbspc.Curve(1,tpoles2d);
1891               const gp_Pln&  Pln = myHS1->Surface().Plane();
1892               //
1893               Standard_Integer ik;
1894               for(ik = 1; ik<= nbpoles; ik++) {
1895                 tpoles.SetValue(ik,
1896                                 ElSLib::Value(tpoles2d.Value(ik).X(),
1897                                               tpoles2d.Value(ik).Y(),
1898                                               Pln));
1899               }
1900               //
1901               Handle(Geom_BSplineCurve) BS =
1902                 new Geom_BSplineCurve(tpoles,
1903                                       mbspc.Knots(),
1904                                       mbspc.Multiplicities(),
1905                                       mbspc.Degree());
1906               GeomLib_CheckBSplineCurve Check(BS,TOLCHECK,TOLANGCHECK);
1907               Check.FixTangent(Standard_True, Standard_True);
1908               //
1909               IntTools_Curve aCurve;
1910               aCurve.SetCurve(BS);
1911
1912               if(myApprox1) {
1913                 Handle(Geom2d_BSplineCurve) BS1 =
1914                   new Geom2d_BSplineCurve(tpoles2d,
1915                                           mbspc.Knots(),
1916                                           mbspc.Multiplicities(),
1917                                           mbspc.Degree());
1918                 GeomLib_Check2dBSplineCurve Check1(BS1,TOLCHECK,TOLANGCHECK);
1919                 Check1.FixTangent(Standard_True,Standard_True);
1920                 //
1921                 // ############################################
1922                 if(!rejectSurface && !reApprox) {
1923                   Standard_Boolean isValid = IsCurveValid(BS1);
1924                   if(!isValid) {
1925                     reApprox = Standard_True;
1926                     goto reapprox;
1927                   }
1928                 }
1929                 // ############################################
1930                 aCurve.SetFirstCurve2d(BS1);
1931               }
1932               else {
1933                 Handle(Geom2d_BSplineCurve) H1;
1934                 aCurve.SetFirstCurve2d(H1);
1935               }
1936
1937               if(myApprox2) {
1938                 mbspc.Curve(2, tpoles2d);
1939
1940                 Handle(Geom2d_BSplineCurve) BS2 = new Geom2d_BSplineCurve(tpoles2d,
1941                                                                           mbspc.Knots(),
1942                                                                           mbspc.Multiplicities(),
1943                                                                           mbspc.Degree());
1944                 GeomLib_Check2dBSplineCurve newCheck(BS2,TOLCHECK,TOLANGCHECK);
1945                 newCheck.FixTangent(Standard_True,Standard_True);
1946
1947                 // ###########################################
1948                 if(!rejectSurface && !reApprox) {
1949                   Standard_Boolean isValid = IsCurveValid(BS2);
1950                   if(!isValid) {
1951                     reApprox = Standard_True;
1952                     goto reapprox;
1953                   }
1954                 }
1955                 // ###########################################
1956                 //
1957                 aCurve.SetSecondCurve2d(BS2);
1958               }
1959               else {
1960                 Handle(Geom2d_BSplineCurve) H2;
1961                 //
1962                   aCurve.SetSecondCurve2d(H2);
1963               }
1964               //
1965               mySeqOfCurve.Append(aCurve);
1966             }//if(typs1 == GeomAbs_Plane) {
1967
1968             else if(typs2 == GeomAbs_Plane) {
1969               const AppParCurves_MultiBSpCurve& mbspc = theapp3d.Value(j);
1970               nbpoles = mbspc.NbPoles();
1971
1972               TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
1973               TColgp_Array1OfPnt   tpoles(1,nbpoles);
1974               mbspc.Curve((myApprox1==Standard_True)? 2 : 1,tpoles2d);
1975               const gp_Pln&  Pln = myHS2->Surface().Plane();
1976               //
1977               Standard_Integer ik;
1978               for(ik = 1; ik<= nbpoles; ik++) {
1979                 tpoles.SetValue(ik,
1980                                 ElSLib::Value(tpoles2d.Value(ik).X(),
1981                                               tpoles2d.Value(ik).Y(),
1982                                               Pln));
1983
1984               }
1985               //
1986               Handle(Geom_BSplineCurve) BS=new Geom_BSplineCurve(tpoles,
1987                                                                  mbspc.Knots(),
1988                                                                  mbspc.Multiplicities(),
1989                                                                  mbspc.Degree());
1990               GeomLib_CheckBSplineCurve Check(BS,TOLCHECK,TOLANGCHECK);
1991               Check.FixTangent(Standard_True,Standard_True);
1992               //
1993               IntTools_Curve aCurve;
1994               aCurve.SetCurve(BS);
1995
1996               if(myApprox2) {
1997                 Handle(Geom2d_BSplineCurve) BS1=new Geom2d_BSplineCurve(tpoles2d,
1998                                                                         mbspc.Knots(),
1999                                                                         mbspc.Multiplicities(),
2000                                                                         mbspc.Degree());
2001                 GeomLib_Check2dBSplineCurve Check1(BS1,TOLCHECK,TOLANGCHECK);
2002                 Check1.FixTangent(Standard_True,Standard_True);
2003                 //
2004                 // ###########################################
2005                 if(!rejectSurface && !reApprox) {
2006                   Standard_Boolean isValid = IsCurveValid(BS1);
2007                   if(!isValid) {
2008                     reApprox = Standard_True;
2009                     goto reapprox;
2010                   }
2011                 }
2012                 // ###########################################
2013                 bPCurvesOk = CheckPCurve(BS1, myFace2);
2014                 aCurve.SetSecondCurve2d(BS1);
2015               }
2016               else {
2017                 Handle(Geom2d_BSplineCurve) H2;
2018                 aCurve.SetSecondCurve2d(H2);
2019               }
2020
2021               if(myApprox1) {
2022                 mbspc.Curve(1,tpoles2d);
2023                 Handle(Geom2d_BSplineCurve) BS2=new Geom2d_BSplineCurve(tpoles2d,
2024                                                                         mbspc.Knots(),
2025                                                                         mbspc.Multiplicities(),
2026                                                                         mbspc.Degree());
2027                 GeomLib_Check2dBSplineCurve Check2(BS2,TOLCHECK,TOLANGCHECK);
2028                 Check2.FixTangent(Standard_True,Standard_True);
2029                 //
2030                 // ###########################################
2031                 if(!rejectSurface && !reApprox) {
2032                   Standard_Boolean isValid = IsCurveValid(BS2);
2033                   if(!isValid) {
2034                     reApprox = Standard_True;
2035                     goto reapprox;
2036                   }
2037                 }
2038                 // ###########################################
2039                 bPCurvesOk = bPCurvesOk && CheckPCurve(BS2, myFace1);
2040                 aCurve.SetFirstCurve2d(BS2);
2041               }
2042               else {
2043                 Handle(Geom2d_BSplineCurve) H1;
2044                 //
2045                 aCurve.SetFirstCurve2d(H1);
2046               }
2047               //
2048               //if points of the pcurves are out of the faces bounds
2049               //create 3d and 2d curves without approximation
2050               if (!bPCurvesOk) {
2051                 Handle(Geom2d_BSplineCurve) H1, H2;
2052                 bPCurvesOk = Standard_True;
2053                 //
2054                 Handle(Geom_Curve) aBSp=MakeBSpline(WL,ifprm, ilprm);
2055
2056                 if(myApprox1) {
2057                   H1 = MakeBSpline2d(WL, ifprm, ilprm, Standard_True);
2058                   bPCurvesOk = CheckPCurve(H1, myFace1);
2059                 }
2060
2061                 if(myApprox2) {
2062                   H2 = MakeBSpline2d(WL, ifprm, ilprm, Standard_False);
2063                   bPCurvesOk = bPCurvesOk && CheckPCurve(H2, myFace2);
2064                 }
2065                 //
2066                 //if pcurves created without approximation are out of the
2067                 //faces bounds, use approximated 3d and 2d curves
2068                 if (bPCurvesOk) {
2069                   IntTools_Curve aIC(aBSp, H1, H2);
2070                   mySeqOfCurve.Append(aIC);
2071                 } else {
2072                   mySeqOfCurve.Append(aCurve);
2073                 }
2074               } else {
2075                 mySeqOfCurve.Append(aCurve);
2076               }
2077             }// else if(typs2 == GeomAbs_Plane)
2078             //
2079             else { //typs2 != GeomAbs_Plane && typs1 != GeomAbs_Plane
2080               Standard_Boolean bIsValid1, bIsValid2;
2081               Handle(Geom_BSplineCurve) BS;
2082               Handle(Geom2d_BSplineCurve) aH2D;
2083               IntTools_Curve aCurve;
2084               //
2085               bIsValid1=Standard_True;
2086               bIsValid2=Standard_True;
2087               //
2088               const AppParCurves_MultiBSpCurve& mbspc = theapp3d.Value(j);
2089               nbpoles = mbspc.NbPoles();
2090               TColgp_Array1OfPnt tpoles(1,nbpoles);
2091               mbspc.Curve(1,tpoles);
2092               BS=new Geom_BSplineCurve(tpoles,
2093                                                                  mbspc.Knots(),
2094                                                                  mbspc.Multiplicities(),
2095                                                                  mbspc.Degree());
2096               GeomLib_CheckBSplineCurve Check(BS,TOLCHECK,TOLANGCHECK);
2097               Check.FixTangent(Standard_True,Standard_True);
2098               //
2099               aCurve.SetCurve(BS);
2100               aCurve.SetFirstCurve2d(aH2D);
2101               aCurve.SetSecondCurve2d(aH2D);
2102               //
2103               if(myApprox1) {
2104                 if(anApprox1) {
2105                   Handle(Geom2d_BSplineCurve) BS1;
2106                   TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
2107                   mbspc.Curve(2,tpoles2d);
2108                   //
2109                   BS1=new Geom2d_BSplineCurve(tpoles2d,
2110                                                                         mbspc.Knots(),
2111                                                                         mbspc.Multiplicities(),
2112                                                                         mbspc.Degree());
2113                   GeomLib_Check2dBSplineCurve newCheck(BS1,TOLCHECK,TOLANGCHECK);
2114                   newCheck.FixTangent(Standard_True,Standard_True);
2115                   //
2116                   if (!reApprox) {
2117                     bIsValid1=CheckPCurve(BS1, myFace1);
2118                   }
2119                   //
2120                   aCurve.SetFirstCurve2d(BS1);
2121                 }
2122                 else {
2123                   Handle(Geom2d_BSplineCurve) BS1;
2124                   fprm = BS->FirstParameter();
2125                   lprm = BS->LastParameter();
2126
2127                   Handle(Geom2d_Curve) C2d;
2128                   Standard_Real aTol = myTolApprox;
2129                   BuildPCurves(fprm, lprm, aTol, myHS1->ChangeSurface().Surface(), BS, C2d);
2130                   BS1 = Handle(Geom2d_BSplineCurve)::DownCast(C2d);
2131                   aCurve.SetFirstCurve2d(BS1);
2132                 }
2133               } // if(myApprox1) {
2134                 //
2135               if(myApprox2) {
2136                 if(anApprox2) {
2137                   Handle(Geom2d_BSplineCurve) BS2;
2138                   TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
2139                   mbspc.Curve((myApprox1==Standard_True)? 3 : 2,tpoles2d);
2140                   BS2=new Geom2d_BSplineCurve(tpoles2d,
2141                                                                         mbspc.Knots(),
2142                                                                         mbspc.Multiplicities(),
2143                                                                         mbspc.Degree());
2144                   GeomLib_Check2dBSplineCurve newCheck(BS2,TOLCHECK,TOLANGCHECK);
2145                   newCheck.FixTangent(Standard_True,Standard_True);
2146                 //
2147                   if (!reApprox) {
2148                     bIsValid2=CheckPCurve(BS2, myFace2);
2149                   }
2150                   aCurve.SetSecondCurve2d(BS2);
2151                 }
2152                 else {
2153                   Handle(Geom2d_BSplineCurve) BS2;
2154                   fprm = BS->FirstParameter();
2155                   lprm = BS->LastParameter();
2156
2157                   Handle(Geom2d_Curve) C2d;
2158                   Standard_Real aTol = myTolApprox;
2159                   BuildPCurves(fprm, lprm, aTol, myHS2->ChangeSurface().Surface(), BS, C2d);
2160                   BS2 = Handle(Geom2d_BSplineCurve)::DownCast(C2d);
2161                   aCurve.SetSecondCurve2d(BS2);
2162                 }
2163               } //if(myApprox2) {
2164               if (!bIsValid1 || !bIsValid2) {
2165                 myTolApprox=aTolApproxImp;//1.e-5;
2166                 tol2d = myTolApprox;
2167                 reApprox = Standard_True;
2168                 goto reapprox;
2169               }
2170                 //
2171               mySeqOfCurve.Append(aCurve);
2172             }
2173           }
2174         }
2175       }
2176     }// else { // X
2177   }// case IntPatch_Walking:{
2178     break;
2179
2180   case IntPatch_Restriction:
2181     break;
2182   default:
2183     break;
2184
2185   }
2186 }
2187
2188 //=======================================================================
2189 //function : BuildPCurves
2190 //purpose  :
2191 //=======================================================================
2192 void BuildPCurves (Standard_Real f,
2193                    Standard_Real l,
2194                    Standard_Real& Tol,
2195                    const Handle (Geom_Surface)& S,
2196                    const Handle (Geom_Curve)&   C,
2197                    Handle (Geom2d_Curve)& C2d)
2198 {
2199   if (!C2d.IsNull()) {
2200     return;
2201   }
2202   //
2203   Standard_Real umin,umax,vmin,vmax;
2204   //
2205   S->Bounds(umin, umax, vmin, vmax);
2206   // in class ProjLib_Function the range of parameters is shrank by 1.e-09
2207   if((l - f) > 2.e-09) {
2208     C2d = GeomProjLib::Curve2d(C,f,l,S,umin,umax,vmin,vmax,Tol);
2209     //
2210     if (C2d.IsNull()) {
2211       // proj. a circle that goes through the pole on a sphere to the sphere
2212       Tol += Precision::Confusion();
2213       C2d = GeomProjLib::Curve2d(C,f,l,S,Tol);
2214     }
2215   }
2216   else {
2217     if((l - f) > Epsilon(Abs(f))) {
2218       GeomAPI_ProjectPointOnSurf aProjector1, aProjector2;
2219       gp_Pnt P1 = C->Value(f);
2220       gp_Pnt P2 = C->Value(l);
2221       aProjector1.Init(P1, S);
2222       aProjector2.Init(P2, S);
2223
2224       if(aProjector1.IsDone() && aProjector2.IsDone()) {
2225         Standard_Real U=0., V=0.;
2226         aProjector1.LowerDistanceParameters(U, V);
2227         gp_Pnt2d p1(U, V);
2228
2229         aProjector2.LowerDistanceParameters(U, V);
2230         gp_Pnt2d p2(U, V);
2231
2232         if(p1.Distance(p2) > gp::Resolution()) {
2233           TColgp_Array1OfPnt2d poles(1,2);
2234           TColStd_Array1OfReal knots(1,2);
2235           TColStd_Array1OfInteger mults(1,2);
2236           poles(1) = p1;
2237           poles(2) = p2;
2238           knots(1) = f;
2239           knots(2) = l;
2240           mults(1) = mults(2) = 2;
2241
2242           C2d = new Geom2d_BSplineCurve(poles,knots,mults,1);
2243
2244           // compute reached tolerance.begin
2245           gp_Pnt PMid = C->Value((f + l) * 0.5);
2246           aProjector1.Perform(PMid);
2247
2248           if(aProjector1.IsDone()) {
2249             aProjector1.LowerDistanceParameters(U, V);
2250             gp_Pnt2d pmidproj(U, V);
2251             gp_Pnt2d pmidcurve2d = C2d->Value((f + l) * 0.5);
2252             Standard_Real adist = pmidcurve2d.Distance(pmidproj);
2253             Tol = (adist > Tol) ? adist : Tol;
2254           }
2255           // compute reached tolerance.end
2256         }
2257       }
2258     }
2259   }
2260   //
2261   if (S->IsUPeriodic() && !C2d.IsNull()) {
2262     // Recadre dans le domaine UV de la face
2263     Standard_Real aTm, U0, aEps, period, du, U0x;
2264     Standard_Boolean bAdjust;
2265     //
2266     aEps = Precision::PConfusion();
2267     period = S->UPeriod();
2268     //
2269     aTm = .5*(f + l);
2270     gp_Pnt2d pm = C2d->Value(aTm);
2271     U0 = pm.X();
2272     //
2273     bAdjust =
2274       GeomInt::AdjustPeriodic(U0, umin, umax, period, U0x, du, aEps);
2275     if (bAdjust) {
2276       gp_Vec2d T1(du, 0.);
2277       C2d->Translate(T1);
2278     }
2279   }
2280 }
2281
2282 //=======================================================================
2283 //function : Parameters
2284 //purpose  :
2285 //=======================================================================
2286  void Parameters(const Handle(GeomAdaptor_HSurface)& HS1,
2287                  const Handle(GeomAdaptor_HSurface)& HS2,
2288                  const gp_Pnt& Ptref,
2289                  Standard_Real& U1,
2290                  Standard_Real& V1,
2291                  Standard_Real& U2,
2292                  Standard_Real& V2)
2293 {
2294
2295   IntSurf_Quadric quad1,quad2;
2296   GeomAbs_SurfaceType typs = HS1->Surface().GetType();
2297
2298   switch (typs) {
2299   case GeomAbs_Plane:
2300     quad1.SetValue(HS1->Surface().Plane());
2301     break;
2302   case GeomAbs_Cylinder:
2303     quad1.SetValue(HS1->Surface().Cylinder());
2304     break;
2305   case GeomAbs_Cone:
2306     quad1.SetValue(HS1->Surface().Cone());
2307     break;
2308   case GeomAbs_Sphere:
2309     quad1.SetValue(HS1->Surface().Sphere());
2310     break;
2311   case GeomAbs_Torus:
2312     quad1.SetValue(HS1->Surface().Torus());
2313     break;
2314   default:
2315     Standard_ConstructionError::Raise("GeomInt_IntSS::MakeCurve");
2316   }
2317
2318   typs = HS2->Surface().GetType();
2319   switch (typs) {
2320   case GeomAbs_Plane:
2321     quad2.SetValue(HS2->Surface().Plane());
2322     break;
2323   case GeomAbs_Cylinder:
2324     quad2.SetValue(HS2->Surface().Cylinder());
2325     break;
2326   case GeomAbs_Cone:
2327     quad2.SetValue(HS2->Surface().Cone());
2328     break;
2329   case GeomAbs_Sphere:
2330     quad2.SetValue(HS2->Surface().Sphere());
2331     break;
2332   case GeomAbs_Torus:
2333     quad2.SetValue(HS2->Surface().Torus());
2334     break;
2335   default:
2336     Standard_ConstructionError::Raise("GeomInt_IntSS::MakeCurve");
2337   }
2338
2339   quad1.Parameters(Ptref,U1,V1);
2340   quad2.Parameters(Ptref,U2,V2);
2341 }
2342
2343 //=======================================================================
2344 //function : MakeBSpline
2345 //purpose  :
2346 //=======================================================================
2347 Handle(Geom_Curve) MakeBSpline  (const Handle(IntPatch_WLine)& WL,
2348                                  const Standard_Integer ideb,
2349                                  const Standard_Integer ifin)
2350 {
2351   Standard_Integer i,nbpnt = ifin-ideb+1;
2352   TColgp_Array1OfPnt poles(1,nbpnt);
2353   TColStd_Array1OfReal knots(1,nbpnt);
2354   TColStd_Array1OfInteger mults(1,nbpnt);
2355   Standard_Integer ipidebm1;
2356   for(i=1,ipidebm1=i+ideb-1; i<=nbpnt;ipidebm1++, i++) {
2357     poles(i) = WL->Point(ipidebm1).Value();
2358     mults(i) = 1;
2359     knots(i) = i-1;
2360   }
2361   mults(1) = mults(nbpnt) = 2;
2362   return
2363     new Geom_BSplineCurve(poles,knots,mults,1);
2364 }
2365 //
2366
2367 //=======================================================================
2368 //function : MakeBSpline2d
2369 //purpose  :
2370 //=======================================================================
2371 Handle(Geom2d_BSplineCurve) MakeBSpline2d(const Handle(IntPatch_WLine)& theWLine,
2372                                           const Standard_Integer ideb,
2373                                           const Standard_Integer ifin,
2374                                           const Standard_Boolean onFirst)
2375 {
2376   Standard_Integer i, nbpnt = ifin-ideb+1;
2377   TColgp_Array1OfPnt2d poles(1,nbpnt);
2378   TColStd_Array1OfReal knots(1,nbpnt);
2379   TColStd_Array1OfInteger mults(1,nbpnt);
2380   Standard_Integer ipidebm1;
2381
2382   for(i = 1, ipidebm1 = i+ideb-1; i <= nbpnt; ipidebm1++, i++) {
2383       Standard_Real U, V;
2384       if(onFirst)
2385         theWLine->Point(ipidebm1).ParametersOnS1(U, V);
2386       else
2387         theWLine->Point(ipidebm1).ParametersOnS2(U, V);
2388       poles(i).SetCoord(U, V);
2389       mults(i) = 1;
2390       knots(i) = i-1;
2391     }
2392     mults(1) = mults(nbpnt) = 2;
2393
2394   return new Geom2d_BSplineCurve(poles,knots,mults,1);
2395 }
2396 //=======================================================================
2397 //function : PrepareLines3D
2398 //purpose  :
2399 //=======================================================================
2400   void IntTools_FaceFace::PrepareLines3D(const Standard_Boolean bToSplit)
2401 {
2402   Standard_Integer i, aNbCurves;
2403   GeomAbs_SurfaceType aType1, aType2;
2404   IntTools_SequenceOfCurves aNewCvs;
2405   //
2406   // 1. Treatment closed  curves
2407   aNbCurves=mySeqOfCurve.Length();
2408   for (i=1; i<=aNbCurves; ++i) {
2409     const IntTools_Curve& aIC=mySeqOfCurve(i);
2410     //
2411     if (bToSplit) {
2412       Standard_Integer j, aNbC;
2413       IntTools_SequenceOfCurves aSeqCvs;
2414       //
2415       aNbC=IntTools_Tools::SplitCurve(aIC, aSeqCvs);
2416       if (aNbC) {
2417         for (j=1; j<=aNbC; ++j) {
2418           const IntTools_Curve& aICNew=aSeqCvs(j);
2419           aNewCvs.Append(aICNew);
2420         }
2421       }
2422       else {
2423         aNewCvs.Append(aIC);
2424       }
2425     }
2426     else {
2427       aNewCvs.Append(aIC);
2428     }
2429   }
2430   //
2431   // 2. Plane\Cone intersection when we had 4 curves
2432   aType1=myHS1->GetType();
2433   aType2=myHS2->GetType();
2434   aNbCurves=aNewCvs.Length();
2435   //
2436   if ((aType1==GeomAbs_Plane && aType2==GeomAbs_Cone) ||
2437       (aType2==GeomAbs_Plane && aType1==GeomAbs_Cone)) {
2438     if (aNbCurves==4) {
2439       GeomAbs_CurveType aCType1;
2440       //
2441       aCType1=aNewCvs(1).Type();
2442       if (aCType1==GeomAbs_Line) {
2443         IntTools_SequenceOfCurves aSeqIn, aSeqOut;
2444         //
2445         for (i=1; i<=aNbCurves; ++i) {
2446           const IntTools_Curve& aIC=aNewCvs(i);
2447           aSeqIn.Append(aIC);
2448         }
2449         //
2450         IntTools_Tools::RejectLines(aSeqIn, aSeqOut);
2451         //
2452         aNewCvs.Clear();
2453         aNbCurves=aSeqOut.Length();
2454         for (i=1; i<=aNbCurves; ++i) {
2455           const IntTools_Curve& aIC=aSeqOut(i);
2456           aNewCvs.Append(aIC);
2457         }
2458       }
2459     }
2460   }// if ((aType1==GeomAbs_Plane && aType2==GeomAbs_Cone)...
2461   //
2462   // 3. Fill  mySeqOfCurve
2463   mySeqOfCurve.Clear();
2464   aNbCurves=aNewCvs.Length();
2465   for (i=1; i<=aNbCurves; ++i) {
2466     const IntTools_Curve& aIC=aNewCvs(i);
2467     mySeqOfCurve.Append(aIC);
2468   }
2469 }
2470 //=======================================================================
2471 //function : CorrectSurfaceBoundaries
2472 //purpose  :
2473 //=======================================================================
2474  void CorrectSurfaceBoundaries(const TopoDS_Face&  theFace,
2475                               const Standard_Real theTolerance,
2476                               Standard_Real&      theumin,
2477                               Standard_Real&      theumax,
2478                               Standard_Real&      thevmin,
2479                               Standard_Real&      thevmax)
2480 {
2481   Standard_Boolean enlarge, isuperiodic, isvperiodic;
2482   Standard_Real uinf, usup, vinf, vsup, delta;
2483   GeomAbs_SurfaceType aType;
2484   Handle(Geom_Surface) aSurface;
2485   //
2486   aSurface = BRep_Tool::Surface(theFace);
2487   aSurface->Bounds(uinf, usup, vinf, vsup);
2488   delta = theTolerance;
2489   enlarge = Standard_False;
2490   //
2491   GeomAdaptor_Surface anAdaptorSurface(aSurface);
2492   //
2493   if(aSurface->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface))) {
2494     Handle(Geom_Surface) aBasisSurface =
2495       (Handle(Geom_RectangularTrimmedSurface)::DownCast(aSurface))->BasisSurface();
2496
2497     if(aBasisSurface->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface)) ||
2498        aBasisSurface->IsKind(STANDARD_TYPE(Geom_OffsetSurface))) {
2499       return;
2500     }
2501   }
2502   //
2503   if(aSurface->IsKind(STANDARD_TYPE(Geom_OffsetSurface))) {
2504     Handle(Geom_Surface) aBasisSurface =
2505       (Handle(Geom_OffsetSurface)::DownCast(aSurface))->BasisSurface();
2506
2507     if(aBasisSurface->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface)) ||
2508        aBasisSurface->IsKind(STANDARD_TYPE(Geom_OffsetSurface))) {
2509       return;
2510     }
2511   }
2512   //
2513   isuperiodic = anAdaptorSurface.IsUPeriodic();
2514   isvperiodic = anAdaptorSurface.IsVPeriodic();
2515   //
2516   aType=anAdaptorSurface.GetType();
2517   if((aType==GeomAbs_BezierSurface) ||
2518      (aType==GeomAbs_BSplineSurface) ||
2519      (aType==GeomAbs_SurfaceOfExtrusion) ||
2520      (aType==GeomAbs_SurfaceOfRevolution) ||
2521      (aType==GeomAbs_Cylinder)) {
2522     enlarge=Standard_True;
2523   }
2524   //
2525   if(!isuperiodic && enlarge) {
2526
2527     if((theumin - uinf) > delta )
2528       theumin -= delta;
2529     else {
2530       theumin = uinf;
2531     }
2532
2533     if((usup - theumax) > delta )
2534       theumax += delta;
2535     else
2536       theumax = usup;
2537   }
2538   //
2539   if(!isvperiodic && enlarge) {
2540     if((thevmin - vinf) > delta ) {
2541       thevmin -= delta;
2542     }
2543     else {
2544       thevmin = vinf;
2545     }
2546     if((vsup - thevmax) > delta ) {
2547       thevmax += delta;
2548     }
2549     else {
2550       thevmax = vsup;
2551     }
2552   }
2553   //
2554   {
2555     Standard_Integer aNbP;
2556     Standard_Real aXP, dXfact, aXmid, aX1, aX2, aTolPA;
2557     //
2558     aTolPA=Precision::Angular();
2559     // U
2560     if (isuperiodic) {
2561       aXP=anAdaptorSurface.UPeriod();
2562       dXfact=theumax-theumin;
2563       if (dXfact-aTolPA>aXP) {
2564         aXmid=0.5*(theumax+theumin);
2565         aNbP=RealToInt(aXmid/aXP);
2566         if (aXmid<0.) {
2567           aNbP=aNbP-1;
2568         }
2569         aX1=aNbP*aXP;
2570         if (theumin>aTolPA) {
2571           aX1=theumin+aNbP*aXP;
2572         }
2573         aX2=aX1+aXP;
2574         if (theumin<aX1) {
2575           theumin=aX1;
2576         }
2577         if (theumax>aX2) {
2578           theumax=aX2;
2579         }
2580       }
2581     }
2582     // V
2583     if (isvperiodic) {
2584       aXP=anAdaptorSurface.VPeriod();
2585       dXfact=thevmax-thevmin;
2586       if (dXfact-aTolPA>aXP) {
2587         aXmid=0.5*(thevmax+thevmin);
2588         aNbP=RealToInt(aXmid/aXP);
2589         if (aXmid<0.) {
2590           aNbP=aNbP-1;
2591         }
2592         aX1=aNbP*aXP;
2593         if (thevmin>aTolPA) {
2594           aX1=thevmin+aNbP*aXP;
2595         }
2596         aX2=aX1+aXP;
2597         if (thevmin<aX1) {
2598           thevmin=aX1;
2599         }
2600         if (thevmax>aX2) {
2601           thevmax=aX2;
2602         }
2603       }
2604     }
2605   }
2606   //
2607   if(isuperiodic || isvperiodic) {
2608     Standard_Boolean correct = Standard_False;
2609     Standard_Boolean correctU = Standard_False;
2610     Standard_Boolean correctV = Standard_False;
2611     Bnd_Box2d aBox;
2612     TopExp_Explorer anExp;
2613
2614     for(anExp.Init(theFace, TopAbs_EDGE); anExp.More(); anExp.Next()) {
2615       if(BRep_Tool::IsClosed(TopoDS::Edge(anExp.Current()), theFace)) {
2616         correct = Standard_True;
2617         Standard_Real f, l;
2618         TopoDS_Edge anEdge = TopoDS::Edge(anExp.Current());
2619
2620         for(Standard_Integer i = 0; i < 2; i++) {
2621           if(i==0) {
2622             anEdge.Orientation(TopAbs_FORWARD);
2623           }
2624           else {
2625             anEdge.Orientation(TopAbs_REVERSED);
2626           }
2627           Handle(Geom2d_Curve) aCurve = BRep_Tool::CurveOnSurface(anEdge, theFace, f, l);
2628
2629           if(aCurve.IsNull()) {
2630             correct = Standard_False;
2631             break;
2632           }
2633           Handle(Geom2d_Line) aLine = Handle(Geom2d_Line)::DownCast(aCurve);
2634
2635           if(aLine.IsNull()) {
2636             correct = Standard_False;
2637             break;
2638           }
2639           gp_Dir2d anUDir(1., 0.);
2640           gp_Dir2d aVDir(0., 1.);
2641           Standard_Real anAngularTolerance = Precision::Angular();
2642
2643           correctU = correctU || aLine->Position().Direction().IsParallel(aVDir, anAngularTolerance);
2644           correctV = correctV || aLine->Position().Direction().IsParallel(anUDir, anAngularTolerance);
2645
2646           gp_Pnt2d pp1 = aCurve->Value(f);
2647           aBox.Add(pp1);
2648           gp_Pnt2d pp2 = aCurve->Value(l);
2649           aBox.Add(pp2);
2650         }
2651         if(!correct)
2652           break;
2653       }
2654     }
2655
2656     if(correct) {
2657       Standard_Real umin, vmin, umax, vmax;
2658       aBox.Get(umin, vmin, umax, vmax);
2659
2660       if(isuperiodic && correctU) {
2661
2662         if(theumin < umin)
2663           theumin = umin;
2664
2665         if(theumax > umax) {
2666           theumax = umax;
2667         }
2668       }
2669       if(isvperiodic && correctV) {
2670
2671         if(thevmin < vmin)
2672           thevmin = vmin;
2673         if(thevmax > vmax)
2674           thevmax = vmax;
2675       }
2676     }
2677   }
2678 }
2679 //
2680 //
2681 // The block is dedicated to determine whether WLine [ifprm, ilprm]
2682 // crosses the degenerated zone on each given surface or not.
2683 // If Yes -> We will not use info about surfaces during approximation
2684 // because inside degenerated zone of the surface the approx. algo.
2685 // uses wrong values of normal, etc., and resulting curve will have
2686 // oscillations that we would not like to have.
2687
2688
2689
2690 static
2691   Standard_Boolean IsDegeneratedZone(const gp_Pnt2d& aP2d,
2692                                      const Handle(Geom_Surface)& aS,
2693                                      const Standard_Integer iDir);
2694 static
2695   Standard_Boolean IsPointInDegeneratedZone(const IntSurf_PntOn2S& aP2S,
2696                                             const TopoDS_Face& aF1,
2697                                             const TopoDS_Face& aF2);
2698 //=======================================================================
2699 //function :  NotUseSurfacesForApprox
2700 //purpose  :
2701 //=======================================================================
2702 Standard_Boolean NotUseSurfacesForApprox(const TopoDS_Face& aF1,
2703                                          const TopoDS_Face& aF2,
2704                                          const Handle(IntPatch_WLine)& WL,
2705                                          const Standard_Integer ifprm,
2706                                          const Standard_Integer ilprm)
2707 {
2708   Standard_Boolean bPInDZ;
2709
2710   Handle(IntSurf_LineOn2S) aLineOn2S=WL->Curve();
2711
2712   const IntSurf_PntOn2S& aP2Sfprm=aLineOn2S->Value(ifprm);
2713   bPInDZ=IsPointInDegeneratedZone(aP2Sfprm, aF1, aF2);
2714   if (bPInDZ) {
2715     return bPInDZ;
2716   }
2717
2718   const IntSurf_PntOn2S& aP2Slprm=aLineOn2S->Value(ilprm);
2719   bPInDZ=IsPointInDegeneratedZone(aP2Slprm, aF1, aF2);
2720
2721   return bPInDZ;
2722 }
2723 //=======================================================================
2724 //function : IsPointInDegeneratedZone
2725 //purpose  :
2726 //=======================================================================
2727 Standard_Boolean IsPointInDegeneratedZone(const IntSurf_PntOn2S& aP2S,
2728                                           const TopoDS_Face& aF1,
2729                                           const TopoDS_Face& aF2)
2730
2731 {
2732   Standard_Boolean bFlag=Standard_True;
2733   Standard_Real US11, US12, VS11, VS12, US21, US22, VS21, VS22;
2734   Standard_Real U1, V1, U2, V2, aDelta, aD;
2735   gp_Pnt2d aP2d;
2736
2737   Handle(Geom_Surface)aS1 = BRep_Tool::Surface(aF1);
2738   aS1->Bounds(US11, US12, VS11, VS12);
2739   GeomAdaptor_Surface aGAS1(aS1);
2740
2741   Handle(Geom_Surface)aS2 = BRep_Tool::Surface(aF2);
2742   aS1->Bounds(US21, US22, VS21, VS22);
2743   GeomAdaptor_Surface aGAS2(aS2);
2744   //
2745   //const gp_Pnt& aP=aP2S.Value();
2746   aP2S.Parameters(U1, V1, U2, V2);
2747   //
2748   aDelta=1.e-7;
2749   // Check on Surf 1
2750   aD=aGAS1.UResolution(aDelta);
2751   aP2d.SetCoord(U1, V1);
2752   if (fabs(U1-US11) < aD) {
2753     bFlag=IsDegeneratedZone(aP2d, aS1, 1);
2754     if (bFlag) {
2755       return bFlag;
2756     }
2757   }
2758   if (fabs(U1-US12) < aD) {
2759     bFlag=IsDegeneratedZone(aP2d, aS1, 1);
2760     if (bFlag) {
2761       return bFlag;
2762     }
2763   }
2764   aD=aGAS1.VResolution(aDelta);
2765   if (fabs(V1-VS11) < aDelta) {
2766     bFlag=IsDegeneratedZone(aP2d, aS1, 2);
2767     if (bFlag) {
2768       return bFlag;
2769     }
2770   }
2771   if (fabs(V1-VS12) < aDelta) {
2772     bFlag=IsDegeneratedZone(aP2d, aS1, 2);
2773     if (bFlag) {
2774       return bFlag;
2775     }
2776   }
2777   // Check on Surf 2
2778   aD=aGAS2.UResolution(aDelta);
2779   aP2d.SetCoord(U2, V2);
2780   if (fabs(U2-US21) < aDelta) {
2781     bFlag=IsDegeneratedZone(aP2d, aS2, 1);
2782     if (bFlag) {
2783       return bFlag;
2784     }
2785   }
2786   if (fabs(U2-US22) < aDelta) {
2787     bFlag=IsDegeneratedZone(aP2d, aS2, 1);
2788     if (bFlag) {
2789       return bFlag;
2790     }
2791   }
2792   aD=aGAS2.VResolution(aDelta);
2793   if (fabs(V2-VS21) < aDelta) {
2794     bFlag=IsDegeneratedZone(aP2d, aS2, 2);
2795     if (bFlag) {
2796       return bFlag;
2797     }
2798   }
2799   if (fabs(V2-VS22) < aDelta) {
2800     bFlag=IsDegeneratedZone(aP2d, aS2, 2);
2801     if (bFlag) {
2802       return bFlag;
2803     }
2804   }
2805   return !bFlag;
2806 }
2807
2808 //=======================================================================
2809 //function : IsDegeneratedZone
2810 //purpose  :
2811 //=======================================================================
2812 Standard_Boolean IsDegeneratedZone(const gp_Pnt2d& aP2d,
2813                                    const Handle(Geom_Surface)& aS,
2814                                    const Standard_Integer iDir)
2815 {
2816   Standard_Boolean bFlag=Standard_True;
2817   Standard_Real US1, US2, VS1, VS2, dY, dX, d1, d2, dD;
2818   Standard_Real aXm, aYm, aXb, aYb, aXe, aYe;
2819   aS->Bounds(US1, US2, VS1, VS2);
2820
2821   gp_Pnt aPm, aPb, aPe;
2822
2823   aXm=aP2d.X();
2824   aYm=aP2d.Y();
2825
2826   aS->D0(aXm, aYm, aPm);
2827
2828   dX=1.e-5;
2829   dY=1.e-5;
2830   dD=1.e-12;
2831
2832   if (iDir==1) {
2833     aXb=aXm;
2834     aXe=aXm;
2835     aYb=aYm-dY;
2836     if (aYb < VS1) {
2837       aYb=VS1;
2838     }
2839     aYe=aYm+dY;
2840     if (aYe > VS2) {
2841       aYe=VS2;
2842     }
2843     aS->D0(aXb, aYb, aPb);
2844     aS->D0(aXe, aYe, aPe);
2845
2846     d1=aPm.Distance(aPb);
2847     d2=aPm.Distance(aPe);
2848     if (d1 < dD && d2 < dD) {
2849       return bFlag;
2850     }
2851     return !bFlag;
2852   }
2853   //
2854   else if (iDir==2) {
2855     aYb=aYm;
2856     aYe=aYm;
2857     aXb=aXm-dX;
2858     if (aXb < US1) {
2859       aXb=US1;
2860     }
2861     aXe=aXm+dX;
2862     if (aXe > US2) {
2863       aXe=US2;
2864     }
2865     aS->D0(aXb, aYb, aPb);
2866     aS->D0(aXe, aYe, aPe);
2867
2868     d1=aPm.Distance(aPb);
2869     d2=aPm.Distance(aPe);
2870     if (d1 < dD && d2 < dD) {
2871       return bFlag;
2872     }
2873     return !bFlag;
2874   }
2875   return !bFlag;
2876 }
2877
2878 //=========================================================================
2879 // static function : ComputePurgedWLine
2880 // purpose : Removes equal points (leave one of equal points) from theWLine
2881 //           and recompute vertex parameters.
2882 //           Returns new WLine or null WLine if the number
2883 //           of the points is less than 2.
2884 //=========================================================================
2885 Handle(IntPatch_WLine) ComputePurgedWLine(const Handle(IntPatch_WLine)& theWLine) {
2886
2887   Standard_Integer i, k, v, nb, nbvtx;
2888   Handle(IntPatch_WLine) aResult;
2889   nbvtx = theWLine->NbVertex();
2890   nb = theWLine->NbPnts();
2891   if (nb==2) {
2892     const IntSurf_PntOn2S& p1 = theWLine->Point(1);
2893     const IntSurf_PntOn2S& p2 = theWLine->Point(2);
2894     if(p1.Value().IsEqual(p2.Value(), gp::Resolution())) {
2895       return aResult;
2896     }
2897   }
2898   //
2899   Handle(IntPatch_WLine) aLocalWLine;
2900   Handle(IntPatch_WLine) aTmpWLine = theWLine;
2901   Handle(IntSurf_LineOn2S) aLineOn2S = new IntSurf_LineOn2S();
2902   aLocalWLine = new IntPatch_WLine(aLineOn2S, Standard_False);
2903   for(i = 1; i <= nb; i++) {
2904     aLineOn2S->Add(theWLine->Point(i));
2905   }
2906
2907   for(v = 1; v <= nbvtx; v++) {
2908     aLocalWLine->AddVertex(theWLine->Vertex(v));
2909   }
2910
2911   for(i = 1; i <= aLineOn2S->NbPoints(); i++) {
2912     Standard_Integer aStartIndex = i + 1;
2913     Standard_Integer anEndIndex = i + 5;
2914     nb = aLineOn2S->NbPoints();
2915     anEndIndex = (anEndIndex > nb) ? nb : anEndIndex;
2916
2917     if((aStartIndex > nb) || (anEndIndex <= 1)) {
2918       continue;
2919     }
2920     k = aStartIndex;
2921
2922     while(k <= anEndIndex) {
2923
2924       if(i != k) {
2925         IntSurf_PntOn2S p1 = aLineOn2S->Value(i);
2926         IntSurf_PntOn2S p2 = aLineOn2S->Value(k);
2927
2928         if(p1.Value().IsEqual(p2.Value(), gp::Resolution())) {
2929           aTmpWLine = aLocalWLine;
2930           aLocalWLine = new IntPatch_WLine(aLineOn2S, Standard_False);
2931
2932           for(v = 1; v <= aTmpWLine->NbVertex(); v++) {
2933             IntPatch_Point aVertex = aTmpWLine->Vertex(v);
2934             Standard_Integer avertexindex = (Standard_Integer)aVertex.ParameterOnLine();
2935
2936             if(avertexindex >= k) {
2937               aVertex.SetParameter(aVertex.ParameterOnLine() - 1.);
2938             }
2939             aLocalWLine->AddVertex(aVertex);
2940           }
2941           aLineOn2S->RemovePoint(k);
2942           anEndIndex--;
2943           continue;
2944         }
2945       }
2946       k++;
2947     }
2948   }
2949
2950   if(aLineOn2S->NbPoints() > 1) {
2951     aResult = aLocalWLine;
2952   }
2953   return aResult;
2954 }
2955
2956 //=======================================================================
2957 //function : TolR3d
2958 //purpose  :
2959 //=======================================================================
2960 void TolR3d(const TopoDS_Face& aF1,
2961             const TopoDS_Face& aF2,
2962             Standard_Real& myTolReached3d)
2963 {
2964   Standard_Real aTolF1, aTolF2, aTolFMax, aTolTresh;
2965
2966   aTolTresh=2.999999e-3;
2967   aTolF1 = BRep_Tool::Tolerance(aF1);
2968   aTolF2 = BRep_Tool::Tolerance(aF2);
2969   aTolFMax=Max(aTolF1, aTolF2);
2970
2971   if (aTolFMax>aTolTresh) {
2972     myTolReached3d=aTolFMax;
2973   }
2974 }
2975 //=======================================================================
2976 //function : IsPointOnBoundary
2977 //purpose  :
2978 //=======================================================================
2979 Standard_Boolean IsPointOnBoundary(const Standard_Real theParameter,
2980                                    const Standard_Real theFirstBoundary,
2981                                    const Standard_Real theSecondBoundary,
2982                                    const Standard_Real theResolution,
2983                                    Standard_Boolean&   IsOnFirstBoundary)
2984 {
2985   Standard_Boolean bRet;
2986   Standard_Integer i;
2987   Standard_Real adist;
2988   //
2989   bRet=Standard_False;
2990   for(i = 0; i < 2; ++i) {
2991     IsOnFirstBoundary = (i == 0);
2992     if (IsOnFirstBoundary) {
2993       adist = fabs(theParameter - theFirstBoundary);
2994     }
2995     else {
2996       adist = fabs(theParameter - theSecondBoundary);
2997     }
2998     if(adist < theResolution) {
2999       return !bRet;
3000     }
3001   }
3002   return bRet;
3003 }
3004 // ------------------------------------------------------------------------------------------------
3005 // static function: FindPoint
3006 // purpose:
3007 // ------------------------------------------------------------------------------------------------
3008 Standard_Boolean FindPoint(const gp_Pnt2d&     theFirstPoint,
3009                            const gp_Pnt2d&     theLastPoint,
3010                            const Standard_Real theUmin,
3011                            const Standard_Real theUmax,
3012                            const Standard_Real theVmin,
3013                            const Standard_Real theVmax,
3014                            gp_Pnt2d&           theNewPoint) {
3015
3016   gp_Vec2d aVec(theFirstPoint, theLastPoint);
3017   Standard_Integer i = 0, j = 0;
3018
3019   for(i = 0; i < 4; i++) {
3020     gp_Vec2d anOtherVec;
3021     gp_Vec2d anOtherVecNormal;
3022     gp_Pnt2d aprojpoint = theLastPoint;
3023
3024     if((i % 2) == 0) {
3025       anOtherVec.SetX(0.);
3026       anOtherVec.SetY(1.);
3027       anOtherVecNormal.SetX(1.);
3028       anOtherVecNormal.SetY(0.);
3029
3030       if(i < 2)
3031         aprojpoint.SetX(theUmin);
3032       else
3033         aprojpoint.SetX(theUmax);
3034     }
3035     else {
3036       anOtherVec.SetX(1.);
3037       anOtherVec.SetY(0.);
3038       anOtherVecNormal.SetX(0.);
3039       anOtherVecNormal.SetY(1.);
3040
3041       if(i < 2)
3042         aprojpoint.SetY(theVmin);
3043       else
3044         aprojpoint.SetY(theVmax);
3045     }
3046     gp_Vec2d anormvec = aVec;
3047     anormvec.Normalize();
3048     RefineVector(anormvec);
3049     Standard_Real adot1 = anormvec.Dot(anOtherVecNormal);
3050
3051     if(fabs(adot1) < Precision::Angular())
3052       continue;
3053     Standard_Real adist = 0.;
3054     Standard_Boolean bIsOut = Standard_False;
3055
3056     if((i % 2) == 0) {
3057       adist = (i < 2) ? fabs(theLastPoint.X() - theUmin) : fabs(theLastPoint.X() - theUmax);
3058       bIsOut = (i < 2) ? (theLastPoint.X() < theUmin) : (theLastPoint.X() > theUmax);
3059     }
3060     else {
3061       adist = (i < 2) ? fabs(theLastPoint.Y() - theVmin) : fabs(theLastPoint.Y() - theVmax);
3062       bIsOut = (i < 2) ? (theLastPoint.Y() < theVmin) : (theLastPoint.Y() > theVmax);
3063     }
3064     Standard_Real anoffset = adist * anOtherVec.Dot(anormvec) / adot1;
3065
3066     for(j = 0; j < 2; j++) {
3067       anoffset = (j == 0) ? anoffset : -anoffset;
3068       gp_Pnt2d acurpoint(aprojpoint.XY() + (anOtherVec.XY()*anoffset));
3069       gp_Vec2d acurvec(theLastPoint, acurpoint);
3070       if ( bIsOut )
3071         acurvec.Reverse();
3072
3073       Standard_Real aDotX, anAngleX;
3074       //
3075       aDotX = aVec.Dot(acurvec);
3076       anAngleX = aVec.Angle(acurvec);
3077       //
3078       if(aDotX > 0. && fabs(anAngleX) < Precision::PConfusion()) {
3079         if((i % 2) == 0) {
3080           if((acurpoint.Y() >= theVmin) &&
3081              (acurpoint.Y() <= theVmax)) {
3082             theNewPoint = acurpoint;
3083             return Standard_True;
3084           }
3085         }
3086         else {
3087           if((acurpoint.X() >= theUmin) &&
3088              (acurpoint.X() <= theUmax)) {
3089             theNewPoint = acurpoint;
3090             return Standard_True;
3091           }
3092         }
3093       }
3094     }
3095   }
3096   return Standard_False;
3097 }
3098
3099
3100 // ------------------------------------------------------------------------------------------------
3101 // static function: FindPoint
3102 // purpose: Find point on the boundary of radial tangent zone
3103 // ------------------------------------------------------------------------------------------------
3104 Standard_Boolean FindPoint(const gp_Pnt2d&     theFirstPoint,
3105                            const gp_Pnt2d&     theLastPoint,
3106                            const Standard_Real theUmin,
3107                            const Standard_Real theUmax,
3108                            const Standard_Real theVmin,
3109                            const Standard_Real theVmax,
3110                            const gp_Pnt2d&     theTanZoneCenter,
3111                            const Standard_Real theZoneRadius,
3112                            Handle(GeomAdaptor_HSurface) theGASurface,
3113                            gp_Pnt2d&           theNewPoint) {
3114   theNewPoint = theLastPoint;
3115
3116   if ( !IsInsideTanZone( theLastPoint, theTanZoneCenter, theZoneRadius, theGASurface) )
3117     return Standard_False;
3118
3119   Standard_Real aUResolution = theGASurface->UResolution( theZoneRadius );
3120   Standard_Real aVResolution = theGASurface->VResolution( theZoneRadius );
3121
3122   Standard_Real aRadius = ( aUResolution < aVResolution ) ? aUResolution : aVResolution;
3123   gp_Ax22d anAxis( theTanZoneCenter, gp_Dir2d(1, 0), gp_Dir2d(0, 1) );
3124   gp_Circ2d aCircle( anAxis, aRadius );
3125
3126   //
3127   gp_Vec2d aDir( theLastPoint.XY() - theFirstPoint.XY() );
3128   Standard_Real aLength = aDir.Magnitude();
3129   if ( aLength <= gp::Resolution() )
3130     return Standard_False;
3131   gp_Lin2d aLine( theFirstPoint, aDir );
3132
3133   //
3134   Handle(Geom2d_Line) aCLine = new Geom2d_Line( aLine );
3135   Handle(Geom2d_TrimmedCurve) aC1 = new Geom2d_TrimmedCurve( aCLine, 0, aLength );
3136   Handle(Geom2d_Circle) aC2 = new Geom2d_Circle( aCircle );
3137
3138   Standard_Real aTol = aRadius * 0.001;
3139   aTol = ( aTol < Precision::PConfusion() ) ? Precision::PConfusion() : aTol;
3140
3141   Geom2dAPI_InterCurveCurve anIntersector;
3142   anIntersector.Init( aC1, aC2, aTol );
3143
3144   if ( anIntersector.NbPoints() == 0 )
3145     return Standard_False;
3146
3147   Standard_Boolean aFound = Standard_False;
3148   Standard_Real aMinDist = aLength * aLength;
3149   Standard_Integer i = 0;
3150   for ( i = 1; i <= anIntersector.NbPoints(); i++ ) {
3151     gp_Pnt2d aPInt = anIntersector.Point( i );
3152     if ( aPInt.SquareDistance( theFirstPoint ) < aMinDist ) {
3153       if ( ( aPInt.X() >= theUmin ) && ( aPInt.X() <= theUmax ) &&
3154            ( aPInt.Y() >= theVmin ) && ( aPInt.Y() <= theVmax ) ) {
3155         theNewPoint = aPInt;
3156         aFound = Standard_True;
3157       }
3158     }
3159   }
3160
3161   return aFound;
3162 }
3163
3164 // ------------------------------------------------------------------------------------------------
3165 // static function: IsInsideTanZone
3166 // purpose: Check if point is inside a radial tangent zone
3167 // ------------------------------------------------------------------------------------------------
3168 Standard_Boolean IsInsideTanZone(const gp_Pnt2d&     thePoint,
3169                                  const gp_Pnt2d&     theTanZoneCenter,
3170                                  const Standard_Real theZoneRadius,
3171                                  Handle(GeomAdaptor_HSurface) theGASurface) {
3172
3173   Standard_Real aUResolution = theGASurface->UResolution( theZoneRadius );
3174   Standard_Real aVResolution = theGASurface->VResolution( theZoneRadius );
3175   Standard_Real aRadiusSQR = ( aUResolution < aVResolution ) ? aUResolution : aVResolution;
3176   aRadiusSQR *= aRadiusSQR;
3177   if ( thePoint.SquareDistance( theTanZoneCenter ) <= aRadiusSQR )
3178     return Standard_True;
3179   return Standard_False;
3180 }
3181
3182 // ------------------------------------------------------------------------------------------------
3183 // static function: CheckTangentZonesExist
3184 // purpose: Check if tangent zone exists
3185 // ------------------------------------------------------------------------------------------------
3186 Standard_Boolean CheckTangentZonesExist( const Handle(GeomAdaptor_HSurface)& theSurface1,
3187                                         const Handle(GeomAdaptor_HSurface)&  theSurface2 )
3188 {
3189   if ( ( theSurface1->GetType() != GeomAbs_Torus ) ||
3190       ( theSurface2->GetType() != GeomAbs_Torus ) )
3191     return Standard_False;
3192
3193   gp_Torus aTor1 = theSurface1->Torus();
3194   gp_Torus aTor2 = theSurface2->Torus();
3195
3196   if ( aTor1.Location().Distance( aTor2.Location() ) > Precision::Confusion() )
3197     return Standard_False;
3198
3199   if ( ( fabs( aTor1.MajorRadius() - aTor2.MajorRadius() ) > Precision::Confusion() ) ||
3200        ( fabs( aTor1.MinorRadius() - aTor2.MinorRadius() ) > Precision::Confusion() ) )
3201     return Standard_False;
3202
3203   if ( ( aTor1.MajorRadius() < aTor1.MinorRadius() ) ||
3204        ( aTor2.MajorRadius() < aTor2.MinorRadius() ) )
3205     return Standard_False;
3206   return Standard_True;
3207 }
3208
3209 // ------------------------------------------------------------------------------------------------
3210 // static function: ComputeTangentZones
3211 // purpose:
3212 // ------------------------------------------------------------------------------------------------
3213 Standard_Integer ComputeTangentZones( const Handle(GeomAdaptor_HSurface)& theSurface1,
3214                                      const Handle(GeomAdaptor_HSurface)&  theSurface2,
3215                                      const TopoDS_Face&                   theFace1,
3216                                      const TopoDS_Face&                   theFace2,
3217                                      Handle(TColgp_HArray1OfPnt2d)&       theResultOnS1,
3218                                      Handle(TColgp_HArray1OfPnt2d)&       theResultOnS2,
3219                                      Handle(TColStd_HArray1OfReal)&       theResultRadius,
3220                                      const Handle(IntTools_Context)& aContext)
3221 {
3222   Standard_Integer aResult = 0;
3223   if ( !CheckTangentZonesExist( theSurface1, theSurface2 ) )
3224     return aResult;
3225
3226
3227   TColgp_SequenceOfPnt2d aSeqResultS1, aSeqResultS2;
3228   TColStd_SequenceOfReal aSeqResultRad;
3229
3230   gp_Torus aTor1 = theSurface1->Torus();
3231   gp_Torus aTor2 = theSurface2->Torus();
3232
3233   gp_Ax2 anax1( aTor1.Location(), aTor1.Axis().Direction() );
3234   gp_Ax2 anax2( aTor2.Location(), aTor2.Axis().Direction() );
3235   Standard_Integer j = 0;
3236
3237   for ( j = 0; j < 2; j++ ) {
3238     Standard_Real aCoef = ( j == 0 ) ? -1 : 1;
3239     Standard_Real aRadius1 = fabs(aTor1.MajorRadius() + aCoef * aTor1.MinorRadius());
3240     Standard_Real aRadius2 = fabs(aTor2.MajorRadius() + aCoef * aTor2.MinorRadius());
3241
3242     gp_Circ aCircle1( anax1, aRadius1 );
3243     gp_Circ aCircle2( anax2, aRadius2 );
3244
3245     // roughly compute radius of tangent zone for perpendicular case
3246     Standard_Real aCriteria = Precision::Confusion() * 0.5;
3247
3248     Standard_Real aT1 = aCriteria;
3249     Standard_Real aT2 = aCriteria;
3250     if ( j == 0 ) {
3251       // internal tangency
3252       Standard_Real aR = ( aRadius1 > aTor2.MinorRadius() ) ? aRadius1 : aTor2.MinorRadius();
3253       //aT1 = aCriteria * aCriteria + aR * aR - ( aR - aCriteria ) * ( aR - aCriteria );
3254       aT1 = 2. * aR * aCriteria;
3255       aT2 = aT1;
3256     }
3257     else {
3258       // external tangency
3259       Standard_Real aRb = ( aRadius1 > aTor2.MinorRadius() ) ? aRadius1 : aTor2.MinorRadius();
3260       Standard_Real aRm = ( aRadius1 < aTor2.MinorRadius() ) ? aRadius1 : aTor2.MinorRadius();
3261       Standard_Real aDelta = aRb - aCriteria;
3262       aDelta *= aDelta;
3263       aDelta -= aRm * aRm;
3264       aDelta /= 2. * (aRb - aRm);
3265       aDelta -= 0.5 * (aRb - aRm);
3266
3267       aT1 = 2. * aRm * (aRm - aDelta);
3268       aT2 = aT1;
3269     }
3270     aCriteria = ( aT1 > aT2) ? aT1 : aT2;
3271     if ( aCriteria > 0 )
3272       aCriteria = sqrt( aCriteria );
3273
3274     if ( aCriteria > 0.5 * aTor1.MinorRadius() ) {
3275       // too big zone -> drop to minimum
3276       aCriteria = Precision::Confusion();
3277     }
3278
3279     GeomAdaptor_Curve aC1( new Geom_Circle(aCircle1) );
3280     GeomAdaptor_Curve aC2( new Geom_Circle(aCircle2) );
3281     Extrema_ExtCC anExtrema(aC1, aC2, 0, 2. * M_PI, 0, 2. * M_PI,
3282                             Precision::PConfusion(), Precision::PConfusion());
3283
3284     if ( anExtrema.IsDone() ) {
3285
3286       Standard_Integer i = 0;
3287       for ( i = 1; i <= anExtrema.NbExt(); i++ ) {
3288         if ( anExtrema.SquareDistance(i) > aCriteria * aCriteria )
3289           continue;
3290
3291         Extrema_POnCurv P1, P2;
3292         anExtrema.Points( i, P1, P2 );
3293
3294         Standard_Boolean bFoundResult = Standard_True;
3295         gp_Pnt2d pr1, pr2;
3296
3297         Standard_Integer surfit = 0;
3298         for ( surfit = 0; surfit < 2; surfit++ ) {
3299           GeomAPI_ProjectPointOnSurf& aProjector =
3300             (surfit == 0) ? aContext->ProjPS(theFace1) : aContext->ProjPS(theFace2);
3301
3302           gp_Pnt aP3d = (surfit == 0) ? P1.Value() : P2.Value();
3303           aProjector.Perform(aP3d);
3304
3305           if(!aProjector.IsDone())
3306             bFoundResult = Standard_False;
3307           else {
3308             if(aProjector.LowerDistance() > aCriteria) {
3309               bFoundResult = Standard_False;
3310             }
3311             else {
3312               Standard_Real foundU = 0, foundV = 0;
3313               aProjector.LowerDistanceParameters(foundU, foundV);
3314               if ( surfit == 0 )
3315                 pr1 = gp_Pnt2d( foundU, foundV );
3316               else
3317                 pr2 = gp_Pnt2d( foundU, foundV );
3318             }
3319           }
3320         }
3321         if ( bFoundResult ) {
3322           aSeqResultS1.Append( pr1 );
3323           aSeqResultS2.Append( pr2 );
3324           aSeqResultRad.Append( aCriteria );
3325
3326           // torus is u and v periodic
3327           const Standard_Real twoPI = M_PI + M_PI;
3328           Standard_Real arr1tmp[2] = {pr1.X(), pr1.Y()};
3329           Standard_Real arr2tmp[2] = {pr2.X(), pr2.Y()};
3330
3331           // iteration on period bounds
3332           for ( Standard_Integer k1 = 0; k1 < 2; k1++ ) {
3333             Standard_Real aBound = ( k1 == 0 ) ? 0 : twoPI;
3334             Standard_Real aShift = ( k1 == 0 ) ? twoPI : -twoPI;
3335
3336             // iteration on surfaces
3337             for ( Standard_Integer k2 = 0; k2 < 2; k2++ ) {
3338               Standard_Real* arr1 = ( k2 == 0 ) ? arr1tmp : arr2tmp;
3339               Standard_Real* arr2 = ( k2 != 0 ) ? arr1tmp : arr2tmp;
3340               TColgp_SequenceOfPnt2d& aSeqS1 = ( k2 == 0 ) ? aSeqResultS1 : aSeqResultS2;
3341               TColgp_SequenceOfPnt2d& aSeqS2 = ( k2 != 0 ) ? aSeqResultS1 : aSeqResultS2;
3342
3343               if (fabs(arr1[0] - aBound) < Precision::PConfusion()) {
3344                 aSeqS1.Append( gp_Pnt2d( arr1[0] + aShift, arr1[1] ) );
3345                 aSeqS2.Append( gp_Pnt2d( arr2[0], arr2[1] ) );
3346                 aSeqResultRad.Append( aCriteria );
3347               }
3348               if (fabs(arr1[1] - aBound) < Precision::PConfusion()) {
3349                 aSeqS1.Append( gp_Pnt2d( arr1[0], arr1[1] + aShift) );
3350                 aSeqS2.Append( gp_Pnt2d( arr2[0], arr2[1] ) );
3351                 aSeqResultRad.Append( aCriteria );
3352               }
3353             }
3354           } //
3355         }
3356       }
3357     }
3358   }
3359   aResult = aSeqResultRad.Length();
3360
3361   if ( aResult > 0 ) {
3362     theResultOnS1 = new TColgp_HArray1OfPnt2d( 1, aResult );
3363     theResultOnS2 = new TColgp_HArray1OfPnt2d( 1, aResult );
3364     theResultRadius = new TColStd_HArray1OfReal( 1, aResult );
3365
3366     for ( Standard_Integer i = 1 ; i <= aResult; i++ ) {
3367       theResultOnS1->SetValue( i, aSeqResultS1.Value(i) );
3368       theResultOnS2->SetValue( i, aSeqResultS2.Value(i) );
3369       theResultRadius->SetValue( i, aSeqResultRad.Value(i) );
3370     }
3371   }
3372   return aResult;
3373 }
3374
3375 // ------------------------------------------------------------------------------------------------
3376 // static function: AdjustByNeighbour
3377 // purpose:
3378 // ------------------------------------------------------------------------------------------------
3379 gp_Pnt2d AdjustByNeighbour(const gp_Pnt2d&     theaNeighbourPoint,
3380                            const gp_Pnt2d&     theOriginalPoint,
3381                            Handle(GeomAdaptor_HSurface) theGASurface) {
3382
3383   gp_Pnt2d ap1 = theaNeighbourPoint;
3384   gp_Pnt2d ap2 = theOriginalPoint;
3385
3386   if ( theGASurface->IsUPeriodic() ) {
3387     Standard_Real aPeriod     = theGASurface->UPeriod();
3388     gp_Pnt2d aPTest = ap2;
3389     Standard_Real aSqDistMin = 1.e+100;
3390
3391     for ( Standard_Integer pIt = -1; pIt <= 1; pIt++) {
3392       aPTest.SetX( theOriginalPoint.X() + aPeriod * pIt );
3393       Standard_Real dd = ap1.SquareDistance( aPTest );
3394
3395       if ( dd < aSqDistMin ) {
3396         ap2 = aPTest;
3397         aSqDistMin = dd;
3398       }
3399     }
3400   }
3401   if ( theGASurface->IsVPeriodic() ) {
3402     Standard_Real aPeriod     = theGASurface->VPeriod();
3403     gp_Pnt2d aPTest = ap2;
3404     Standard_Real aSqDistMin = 1.e+100;
3405
3406     for ( Standard_Integer pIt = -1; pIt <= 1; pIt++) {
3407       aPTest.SetY( theOriginalPoint.Y() + aPeriod * pIt );
3408       Standard_Real dd = ap1.SquareDistance( aPTest );
3409
3410       if ( dd < aSqDistMin ) {
3411         ap2 = aPTest;
3412         aSqDistMin = dd;
3413       }
3414     }
3415   }
3416   return ap2;
3417 }
3418
3419 // ------------------------------------------------------------------------------------------------
3420 //function: DecompositionOfWLine
3421 // purpose:
3422 // ------------------------------------------------------------------------------------------------
3423 Standard_Boolean DecompositionOfWLine(const Handle(IntPatch_WLine)& theWLine,
3424                                       const Handle(GeomAdaptor_HSurface)&            theSurface1,
3425                                       const Handle(GeomAdaptor_HSurface)&            theSurface2,
3426                                       const TopoDS_Face&                             theFace1,
3427                                       const TopoDS_Face&                             theFace2,
3428                                       const GeomInt_LineConstructor&                 theLConstructor,
3429                                       const Standard_Boolean                         theAvoidLConstructor,
3430                                       IntPatch_SequenceOfLine&                       theNewLines,
3431                                       Standard_Real&                                 theReachedTol3d,
3432                                       const Handle(IntTools_Context)& aContext)
3433 {
3434
3435   Standard_Boolean bRet, bAvoidLineConstructor;
3436   Standard_Integer aNbPnts, aNbParts;
3437   //
3438   bRet=Standard_False;
3439   aNbPnts=theWLine->NbPnts();
3440   bAvoidLineConstructor=theAvoidLConstructor;
3441   //
3442   if(!aNbPnts){
3443     return bRet;
3444   }
3445   if (!bAvoidLineConstructor) {
3446     aNbParts=theLConstructor.NbParts();
3447     if (!aNbParts) {
3448       return bRet;
3449     }
3450   }
3451   //
3452   Standard_Boolean bIsPrevPointOnBoundary, bIsPointOnBoundary, bIsCurrentPointOnBoundary;
3453   Standard_Integer nblines, pit, i, j;
3454   Standard_Real aTol;
3455   TColStd_Array1OfListOfInteger anArrayOfLines(1, aNbPnts);
3456   TColStd_Array1OfInteger       anArrayOfLineType(1, aNbPnts);
3457   TColStd_ListOfInteger aListOfPointIndex;
3458
3459   Handle(TColgp_HArray1OfPnt2d) aTanZoneS1;
3460   Handle(TColgp_HArray1OfPnt2d) aTanZoneS2;
3461   Handle(TColStd_HArray1OfReal) aTanZoneRadius;
3462   Standard_Integer aNbZone = ComputeTangentZones( theSurface1, theSurface2, theFace1, theFace2,
3463                                                  aTanZoneS1, aTanZoneS2, aTanZoneRadius, aContext);
3464
3465   //
3466   nblines=0;
3467   aTol=Precision::Confusion();
3468   aTol=0.5*aTol;
3469   bIsPrevPointOnBoundary=Standard_False;
3470   bIsPointOnBoundary=Standard_False;
3471   //
3472   // 1. ...
3473   //
3474   // Points
3475   for(pit = 1; pit <= aNbPnts; ++pit) {
3476     Standard_Boolean bIsOnFirstBoundary, isperiodic;
3477     Standard_Real aResolution, aPeriod, alowerboundary, aupperboundary, U, V;
3478     Standard_Real aParameter, anoffset, anAdjustPar;
3479     Standard_Real umin, umax, vmin, vmax;
3480     //
3481     bIsCurrentPointOnBoundary = Standard_False;
3482     const IntSurf_PntOn2S& aPoint = theWLine->Point(pit);
3483     //
3484     // Surface
3485     for(i = 0; i < 2; ++i) {
3486       Handle(GeomAdaptor_HSurface) aGASurface = (!i) ? theSurface1 : theSurface2;
3487       aGASurface->ChangeSurface().Surface()->Bounds(umin, umax, vmin, vmax);
3488       if(!i) {
3489         aPoint.ParametersOnS1(U, V);
3490       }
3491       else {
3492         aPoint.ParametersOnS2(U, V);
3493       }
3494       // U, V
3495       for(j = 0; j < 2; j++) {
3496         isperiodic = (!j) ? aGASurface->IsUPeriodic() : aGASurface->IsVPeriodic();
3497         if(!isperiodic){
3498           continue;
3499         }
3500         //
3501         if (!j) {
3502           aResolution=aGASurface->UResolution(aTol);
3503           aPeriod=aGASurface->UPeriod();
3504           alowerboundary=umin;
3505           aupperboundary=umax;
3506           aParameter=U;
3507         }
3508         else {
3509           aResolution=aGASurface->VResolution(aTol);
3510           aPeriod=aGASurface->VPeriod();
3511           alowerboundary=vmin;
3512           aupperboundary=vmax;
3513           aParameter=V;
3514         }
3515
3516         GeomInt::AdjustPeriodic(aParameter,
3517                                        alowerboundary,
3518                                        aupperboundary,
3519                                        aPeriod,
3520                                        anAdjustPar,
3521                                        anoffset);
3522         //
3523         bIsOnFirstBoundary = Standard_True;// ?
3524         bIsPointOnBoundary=
3525           IsPointOnBoundary(anAdjustPar,
3526                             alowerboundary,
3527                             aupperboundary,
3528                             aResolution,
3529                             bIsOnFirstBoundary);
3530         //
3531         if(bIsPointOnBoundary) {
3532           bIsCurrentPointOnBoundary = Standard_True;
3533           break;
3534         }
3535         else {
3536           // check if a point belong to a tangent zone. Begin
3537           Standard_Integer zIt = 0;
3538           for ( zIt = 1; zIt <= aNbZone; zIt++ ) {
3539             gp_Pnt2d aPZone = (i == 0) ? aTanZoneS1->Value(zIt) : aTanZoneS2->Value(zIt);
3540             Standard_Real aZoneRadius = aTanZoneRadius->Value(zIt);
3541
3542             if ( IsInsideTanZone(gp_Pnt2d( U, V ), aPZone, aZoneRadius, aGASurface ) ) {
3543               // set boundary flag to split the curve by a tangent zone
3544               bIsPointOnBoundary = Standard_True;
3545               bIsCurrentPointOnBoundary = Standard_True;
3546               if ( theReachedTol3d < aZoneRadius ) {
3547                 theReachedTol3d = aZoneRadius;
3548               }
3549               break;
3550             }
3551           }
3552         }
3553       }//for(j = 0; j < 2; j++) {
3554
3555       if(bIsCurrentPointOnBoundary){
3556         break;
3557       }
3558     }//for(i = 0; i < 2; ++i) {
3559     //
3560     if((bIsCurrentPointOnBoundary != bIsPrevPointOnBoundary)) {
3561       if(!aListOfPointIndex.IsEmpty()) {
3562         nblines++;
3563         anArrayOfLines.SetValue(nblines, aListOfPointIndex);
3564         anArrayOfLineType.SetValue(nblines, bIsPrevPointOnBoundary);
3565         aListOfPointIndex.Clear();
3566       }
3567       bIsPrevPointOnBoundary = bIsCurrentPointOnBoundary;
3568     }
3569     aListOfPointIndex.Append(pit);
3570   } //for(pit = 1; pit <= aNbPnts; ++pit) {
3571   //
3572   if(!aListOfPointIndex.IsEmpty()) {
3573     nblines++;
3574     anArrayOfLines.SetValue(nblines, aListOfPointIndex);
3575     anArrayOfLineType.SetValue(nblines, bIsPrevPointOnBoundary);
3576     aListOfPointIndex.Clear();
3577   }
3578   //
3579   if(nblines<=1) {
3580     return bRet; //Standard_False;
3581   }
3582   //
3583   //
3584   // 2. Correct wlines.begin
3585   TColStd_Array1OfListOfInteger anArrayOfLineEnds(1, nblines);
3586   Handle(IntSurf_LineOn2S) aSeqOfPntOn2S = new IntSurf_LineOn2S();
3587   //
3588   for(i = 1; i <= nblines; i++) {
3589     if(anArrayOfLineType.Value(i) != 0) {
3590       continue;
3591     }
3592     const TColStd_ListOfInteger& aListOfIndex = anArrayOfLines.Value(i);
3593     TColStd_ListOfInteger aListOfFLIndex;
3594
3595     for(j = 0; j < 2; j++) {
3596       Standard_Integer aneighbourindex = (j == 0) ? (i - 1) : (i + 1);
3597
3598       if((aneighbourindex < 1) || (aneighbourindex > nblines))
3599         continue;
3600
3601       if(anArrayOfLineType.Value(aneighbourindex) == 0)
3602         continue;
3603       const TColStd_ListOfInteger& aNeighbour = anArrayOfLines.Value(aneighbourindex);
3604       Standard_Integer anIndex = (j == 0) ? aNeighbour.Last() : aNeighbour.First();
3605       const IntSurf_PntOn2S& aPoint = theWLine->Point(anIndex);
3606
3607       IntSurf_PntOn2S aNewP = aPoint;
3608       if(aListOfIndex.Extent() < 2) {
3609         aSeqOfPntOn2S->Add(aNewP);
3610         aListOfFLIndex.Append(aSeqOfPntOn2S->NbPoints());
3611         continue;
3612       }
3613       //
3614       Standard_Integer iFirst = aListOfIndex.First();
3615       Standard_Integer iLast  = aListOfIndex.Last();
3616       //
3617       for(Standard_Integer surfit = 0; surfit < 2; surfit++) {
3618
3619         Handle(GeomAdaptor_HSurface) aGASurface = (surfit == 0) ? theSurface1 : theSurface2;
3620         Standard_Real umin=0., umax=0., vmin=0., vmax=0.;
3621         aGASurface->ChangeSurface().Surface()->Bounds(umin, umax, vmin, vmax);
3622         Standard_Real U=0., V=0.;
3623
3624         if(surfit == 0)
3625           aNewP.ParametersOnS1(U, V);
3626         else
3627           aNewP.ParametersOnS2(U, V);
3628         Standard_Integer nbboundaries = 0;
3629
3630         Standard_Boolean bIsNearBoundary = Standard_False;
3631         Standard_Integer aZoneIndex = 0;
3632         Standard_Integer bIsUBoundary = Standard_False; // use if nbboundaries == 1
3633         Standard_Integer bIsFirstBoundary = Standard_False; // use if nbboundaries == 1
3634
3635
3636         for(Standard_Integer parit = 0; parit < 2; parit++) {
3637           Standard_Boolean isperiodic = (parit == 0) ? aGASurface->IsUPeriodic() : aGASurface->IsVPeriodic();
3638
3639           Standard_Real aResolution = (parit == 0) ? aGASurface->UResolution(aTol) : aGASurface->VResolution(aTol);
3640           Standard_Real alowerboundary = (parit == 0) ? umin : vmin;
3641           Standard_Real aupperboundary = (parit == 0) ? umax : vmax;
3642
3643           Standard_Real aParameter = (parit == 0) ? U : V;
3644           Standard_Boolean bIsOnFirstBoundary = Standard_True;
3645
3646           if(!isperiodic) {
3647             bIsPointOnBoundary=
3648               IsPointOnBoundary(aParameter, alowerboundary, aupperboundary, aResolution, bIsOnFirstBoundary);
3649             if(bIsPointOnBoundary) {
3650               bIsUBoundary = (parit == 0);
3651               bIsFirstBoundary = bIsOnFirstBoundary;
3652               nbboundaries++;
3653             }
3654           }
3655           else {
3656             Standard_Real aPeriod     = (parit == 0) ? aGASurface->UPeriod() : aGASurface->VPeriod();
3657             Standard_Real anoffset, anAdjustPar;
3658             GeomInt::AdjustPeriodic(aParameter, alowerboundary, aupperboundary,
3659                                            aPeriod, anAdjustPar, anoffset);
3660
3661             bIsPointOnBoundary=
3662               IsPointOnBoundary(anAdjustPar, alowerboundary, aupperboundary, aResolution, bIsOnFirstBoundary);
3663             if(bIsPointOnBoundary) {
3664               bIsUBoundary = (parit == 0);
3665               bIsFirstBoundary = bIsOnFirstBoundary;
3666               nbboundaries++;
3667             }
3668             else {
3669             //check neighbourhood of boundary
3670             Standard_Real anEpsilon = aResolution * 100.;
3671             Standard_Real aPart = ( aupperboundary - alowerboundary ) * 0.1;
3672             anEpsilon = ( anEpsilon > aPart ) ? aPart : anEpsilon;
3673
3674               bIsNearBoundary = IsPointOnBoundary(anAdjustPar, alowerboundary, aupperboundary,
3675                                                   anEpsilon, bIsOnFirstBoundary);
3676
3677             }
3678           }
3679         }
3680
3681         // check if a point belong to a tangent zone. Begin
3682         for ( Standard_Integer zIt = 1; zIt <= aNbZone; zIt++ ) {
3683           gp_Pnt2d aPZone = (surfit == 0) ? aTanZoneS1->Value(zIt) : aTanZoneS2->Value(zIt);
3684           Standard_Real aZoneRadius = aTanZoneRadius->Value(zIt);
3685
3686           Standard_Integer aneighbourpointindex1 = (j == 0) ? iFirst : iLast;
3687           const IntSurf_PntOn2S& aNeighbourPoint = theWLine->Point(aneighbourpointindex1);
3688           Standard_Real nU1, nV1;
3689
3690           if(surfit == 0)
3691             aNeighbourPoint.ParametersOnS1(nU1, nV1);
3692           else
3693             aNeighbourPoint.ParametersOnS2(nU1, nV1);
3694           gp_Pnt2d ap1(nU1, nV1);
3695           gp_Pnt2d ap2 = AdjustByNeighbour( ap1, gp_Pnt2d( U, V ), aGASurface );
3696
3697
3698           if ( IsInsideTanZone( ap2, aPZone, aZoneRadius, aGASurface ) ) {
3699             aZoneIndex = zIt;
3700             bIsNearBoundary = Standard_True;
3701             if ( theReachedTol3d < aZoneRadius ) {
3702               theReachedTol3d = aZoneRadius;
3703             }
3704           }
3705         }
3706         // check if a point belong to a tangent zone. End
3707         Standard_Boolean bComputeLineEnd = Standard_False;
3708
3709         if(nbboundaries == 2) {
3710           //xf
3711           bComputeLineEnd = Standard_True;
3712           //xt
3713         }
3714         else if(nbboundaries == 1) {
3715           Standard_Boolean isperiodic = (bIsUBoundary) ? aGASurface->IsUPeriodic() : aGASurface->IsVPeriodic();
3716
3717           if(isperiodic) {
3718             Standard_Real alowerboundary = (bIsUBoundary) ? umin : vmin;
3719             Standard_Real aupperboundary = (bIsUBoundary) ? umax : vmax;
3720             Standard_Real aPeriod     = (bIsUBoundary) ? aGASurface->UPeriod() : aGASurface->VPeriod();
3721             Standard_Real aParameter = (bIsUBoundary) ? U : V;
3722             Standard_Real anoffset, anAdjustPar;
3723             GeomInt::AdjustPeriodic(aParameter, alowerboundary, aupperboundary,
3724                                            aPeriod, anAdjustPar, anoffset);
3725
3726             Standard_Real adist = (bIsFirstBoundary) ? fabs(anAdjustPar - alowerboundary) : fabs(anAdjustPar - aupperboundary);
3727             Standard_Real anotherPar = (bIsFirstBoundary) ? (aupperboundary - adist) : (alowerboundary + adist);
3728             anotherPar += anoffset;
3729             Standard_Integer aneighbourpointindex = (j == 0) ? aListOfIndex.First() : aListOfIndex.Last();
3730             const IntSurf_PntOn2S& aNeighbourPoint = theWLine->Point(aneighbourpointindex);
3731             Standard_Real nU1, nV1;
3732
3733             if(surfit == 0)
3734               aNeighbourPoint.ParametersOnS1(nU1, nV1);
3735             else
3736               aNeighbourPoint.ParametersOnS2(nU1, nV1);
3737
3738             Standard_Real adist1 = (bIsUBoundary) ? fabs(nU1 - U) : fabs(nV1 - V);
3739             Standard_Real adist2 = (bIsUBoundary) ? fabs(nU1 - anotherPar) : fabs(nV1 - anotherPar);
3740             bComputeLineEnd = Standard_True;
3741             Standard_Boolean bCheckAngle1 = Standard_False;
3742             Standard_Boolean bCheckAngle2 = Standard_False;
3743             gp_Vec2d aNewVec;
3744             Standard_Real anewU = (bIsUBoundary) ? anotherPar : U;
3745             Standard_Real anewV = (bIsUBoundary) ? V : anotherPar;
3746
3747             if(((adist1 - adist2) > Precision::PConfusion()) &&
3748                (adist2 < (aPeriod / 4.))) {
3749               bCheckAngle1 = Standard_True;
3750               aNewVec = gp_Vec2d(gp_Pnt2d(nU1, nV1), gp_Pnt2d(anewU, anewV));
3751
3752               if(aNewVec.SquareMagnitude() < gp::Resolution()) {
3753                 aNewP.SetValue((surfit == 0), anewU, anewV);
3754                 bCheckAngle1 = Standard_False;
3755               }
3756             }
3757             else if(adist1 < (aPeriod / 4.)) {
3758               bCheckAngle2 = Standard_True;
3759               aNewVec = gp_Vec2d(gp_Pnt2d(nU1, nV1), gp_Pnt2d(U, V));
3760
3761               if(aNewVec.SquareMagnitude() < gp::Resolution()) {
3762                 bCheckAngle2 = Standard_False;
3763               }
3764             }
3765
3766             if(bCheckAngle1 || bCheckAngle2) {
3767               // assume there are at least two points in line (see "if" above)
3768               Standard_Integer anindexother = aneighbourpointindex;
3769
3770               while((anindexother <= iLast) && (anindexother >= iFirst)) {
3771                 anindexother = (j == 0) ? (anindexother + 1) : (anindexother - 1);
3772                 const IntSurf_PntOn2S& aPrevNeighbourPoint = theWLine->Point(anindexother);
3773                 Standard_Real nU2, nV2;
3774
3775                 if(surfit == 0)
3776                   aPrevNeighbourPoint.ParametersOnS1(nU2, nV2);
3777                 else
3778                   aPrevNeighbourPoint.ParametersOnS2(nU2, nV2);
3779                 gp_Vec2d aVecOld(gp_Pnt2d(nU2, nV2), gp_Pnt2d(nU1, nV1));
3780
3781                 if(aVecOld.SquareMagnitude() <= gp::Resolution()) {
3782                   continue;
3783                 }
3784                 else {
3785                   Standard_Real anAngle = aNewVec.Angle(aVecOld);
3786
3787                   if((fabs(anAngle) < (M_PI * 0.25)) && (aNewVec.Dot(aVecOld) > 0.)) {
3788
3789                     if(bCheckAngle1) {
3790                       Standard_Real U1, U2, V1, V2;
3791                       IntSurf_PntOn2S atmppoint = aNewP;
3792                       atmppoint.SetValue((surfit == 0), anewU, anewV);
3793                       atmppoint.Parameters(U1, V1, U2, V2);
3794                       gp_Pnt P1 = theSurface1->Value(U1, V1);
3795                       gp_Pnt P2 = theSurface2->Value(U2, V2);
3796                       gp_Pnt P0 = aPoint.Value();
3797
3798                       if(P0.IsEqual(P1, aTol) &&
3799                          P0.IsEqual(P2, aTol) &&
3800                          P1.IsEqual(P2, aTol)) {
3801                         bComputeLineEnd = Standard_False;
3802                         aNewP.SetValue((surfit == 0), anewU, anewV);
3803                       }
3804                     }
3805
3806                     if(bCheckAngle2) {
3807                       bComputeLineEnd = Standard_False;
3808                     }
3809                   }
3810                   break;
3811                 }
3812               } // end while(anindexother...)
3813             }
3814           }
3815         }
3816         else if ( bIsNearBoundary ) {
3817           bComputeLineEnd = Standard_True;
3818         }
3819
3820         if(bComputeLineEnd) {
3821
3822           gp_Pnt2d anewpoint;
3823           Standard_Boolean found = Standard_False;
3824
3825           if ( bIsNearBoundary ) {
3826             // re-compute point near natural boundary or near tangent zone
3827             Standard_Real u1, v1, u2, v2;
3828             aNewP.Parameters( u1, v1, u2, v2 );
3829             if(surfit == 0)
3830               anewpoint = gp_Pnt2d( u1, v1 );
3831             else
3832               anewpoint = gp_Pnt2d( u2, v2 );
3833
3834             Standard_Integer aneighbourpointindex1 = (j == 0) ? iFirst : iLast;
3835             const IntSurf_PntOn2S& aNeighbourPoint = theWLine->Point(aneighbourpointindex1);
3836             Standard_Real nU1, nV1;
3837
3838             if(surfit == 0)
3839               aNeighbourPoint.ParametersOnS1(nU1, nV1);
3840             else
3841               aNeighbourPoint.ParametersOnS2(nU1, nV1);
3842             gp_Pnt2d ap1(nU1, nV1);
3843             gp_Pnt2d ap2;
3844
3845
3846             if ( aZoneIndex ) {
3847               // exclude point from a tangent zone
3848               anewpoint = AdjustByNeighbour( ap1, anewpoint, aGASurface );
3849               gp_Pnt2d aPZone = (surfit == 0) ? aTanZoneS1->Value(aZoneIndex) : aTanZoneS2->Value(aZoneIndex);
3850               Standard_Real aZoneRadius = aTanZoneRadius->Value(aZoneIndex);
3851
3852               if ( FindPoint(ap1, anewpoint, umin, umax, vmin, vmax,
3853                              aPZone, aZoneRadius, aGASurface, ap2) ) {
3854                 anewpoint = ap2;
3855                 found = Standard_True;
3856               }
3857             }
3858             else if ( aGASurface->IsUPeriodic() || aGASurface->IsVPeriodic() ) {
3859               // re-compute point near boundary if shifted on a period
3860               ap2 = AdjustByNeighbour( ap1, anewpoint, aGASurface );
3861
3862               if ( ( ap2.X() < umin ) || ( ap2.X() > umax ) ||
3863                   ( ap2.Y() < vmin ) || ( ap2.Y() > vmax ) ) {
3864                 found = FindPoint(ap1, ap2, umin, umax, vmin, vmax, anewpoint);
3865               }
3866               else {
3867                 anewpoint = ap2;
3868                 aNewP.SetValue( (surfit == 0), anewpoint.X(), anewpoint.Y() );
3869               }
3870             }
3871           }
3872           else {
3873
3874             Standard_Integer aneighbourpointindex1 = (j == 0) ? iFirst : iLast;
3875             const IntSurf_PntOn2S& aNeighbourPoint = theWLine->Point(aneighbourpointindex1);
3876             Standard_Real nU1, nV1;
3877
3878             if(surfit == 0)
3879               aNeighbourPoint.ParametersOnS1(nU1, nV1);
3880             else
3881               aNeighbourPoint.ParametersOnS2(nU1, nV1);
3882             gp_Pnt2d ap1(nU1, nV1);
3883             gp_Pnt2d ap2(nU1, nV1);
3884             Standard_Integer aneighbourpointindex2 = aneighbourpointindex1;
3885
3886             while((aneighbourpointindex2 <= iLast) && (aneighbourpointindex2 >= iFirst)) {
3887               aneighbourpointindex2 = (j == 0) ? (aneighbourpointindex2 + 1) : (aneighbourpointindex2 - 1);
3888               const IntSurf_PntOn2S& aPrevNeighbourPoint = theWLine->Point(aneighbourpointindex2);
3889               Standard_Real nU2, nV2;
3890
3891               if(surfit == 0)
3892                 aPrevNeighbourPoint.ParametersOnS1(nU2, nV2);
3893               else
3894                 aPrevNeighbourPoint.ParametersOnS2(nU2, nV2);
3895               ap2.SetX(nU2);
3896               ap2.SetY(nV2);
3897
3898               if(ap1.SquareDistance(ap2) > gp::Resolution()) {
3899                 break;
3900               }
3901             }
3902             found = FindPoint(ap2, ap1, umin, umax, vmin, vmax, anewpoint);
3903           }
3904
3905           if(found) {
3906             // check point
3907             Standard_Real aCriteria = BRep_Tool::Tolerance(theFace1) + BRep_Tool::Tolerance(theFace2);
3908             GeomAPI_ProjectPointOnSurf& aProjector =
3909               (surfit == 0) ? aContext->ProjPS(theFace2) : aContext->ProjPS(theFace1);
3910             Handle(GeomAdaptor_HSurface) aSurface = (surfit == 0) ? theSurface1 : theSurface2;
3911
3912             Handle(GeomAdaptor_HSurface) aSurfaceOther = (surfit == 0) ? theSurface2 : theSurface1;
3913
3914             gp_Pnt aP3d = aSurface->Value(anewpoint.X(), anewpoint.Y());
3915             aProjector.Perform(aP3d);
3916
3917             if(aProjector.IsDone()) {
3918               if(aProjector.LowerDistance() < aCriteria) {
3919                 Standard_Real foundU = U, foundV = V;
3920                 aProjector.LowerDistanceParameters(foundU, foundV);
3921
3922                 //Correction of projected coordinates. Begin
3923                 //Note, it may be shifted on a period
3924                 Standard_Integer aneindex1 = (j == 0) ? iFirst : iLast;
3925                 const IntSurf_PntOn2S& aNeighbourPoint = theWLine->Point(aneindex1);
3926                 Standard_Real nUn, nVn;
3927
3928                 if(surfit == 0)
3929                   aNeighbourPoint.ParametersOnS2(nUn, nVn);
3930                 else
3931                   aNeighbourPoint.ParametersOnS1(nUn, nVn);
3932                 gp_Pnt2d aNeighbour2d(nUn, nVn);
3933                 gp_Pnt2d anAdjustedPoint = AdjustByNeighbour( aNeighbour2d, gp_Pnt2d(foundU, foundV), aSurfaceOther );
3934                 foundU = anAdjustedPoint.X();
3935                 foundV = anAdjustedPoint.Y();
3936
3937                 if ( ( anAdjustedPoint.X() < umin ) && ( anAdjustedPoint.X() > umax ) &&
3938                     ( anAdjustedPoint.Y() < vmin ) && ( anAdjustedPoint.Y() > vmax ) ) {
3939                   // attempt to roughly re-compute point
3940                   foundU = ( foundU < umin ) ? umin : foundU;
3941                   foundU = ( foundU > umax ) ? umax : foundU;
3942                   foundV = ( foundV < vmin ) ? vmin : foundV;
3943                   foundV = ( foundV > vmax ) ? vmax : foundV;
3944
3945                   GeomAPI_ProjectPointOnSurf& aProjector2 =
3946                     (surfit == 0) ? aContext->ProjPS(theFace1) : aContext->ProjPS(theFace2);
3947
3948                   aP3d = aSurfaceOther->Value(foundU, foundV);
3949                   aProjector2.Perform(aP3d);
3950
3951                   if(aProjector2.IsDone()) {
3952                     if(aProjector2.LowerDistance() < aCriteria) {
3953                       Standard_Real foundU2 = anewpoint.X(), foundV2 = anewpoint.Y();
3954                       aProjector2.LowerDistanceParameters(foundU2, foundV2);
3955                       anewpoint.SetX(foundU2);
3956                       anewpoint.SetY(foundV2);
3957                     }
3958                   }
3959                 }
3960                 //Correction of projected coordinates. End
3961
3962                 if(surfit == 0)
3963                   aNewP.SetValue(aP3d, anewpoint.X(), anewpoint.Y(), foundU, foundV);
3964                 else
3965                   aNewP.SetValue(aP3d, foundU, foundV, anewpoint.X(), anewpoint.Y());
3966               }
3967             }
3968           }
3969         }
3970       }
3971       aSeqOfPntOn2S->Add(aNewP);
3972       aListOfFLIndex.Append(aSeqOfPntOn2S->NbPoints());
3973     }
3974     anArrayOfLineEnds.SetValue(i, aListOfFLIndex);
3975   }
3976   // Correct wlines.end
3977
3978   // Split wlines.begin
3979   Standard_Integer nbiter;
3980   //
3981   nbiter=1;
3982   if (!bAvoidLineConstructor) {
3983     nbiter=theLConstructor.NbParts();
3984   }
3985   //
3986   for(j = 1; j <= nbiter; ++j) {
3987     Standard_Real fprm, lprm;
3988     Standard_Integer ifprm, ilprm;
3989     //
3990     if(bAvoidLineConstructor) {
3991       ifprm = 1;
3992       ilprm = theWLine->NbPnts();
3993     }
3994     else {
3995       theLConstructor.Part(j, fprm, lprm);
3996       ifprm = (Standard_Integer)fprm;
3997       ilprm = (Standard_Integer)lprm;
3998     }
3999
4000     Handle(IntSurf_LineOn2S) aLineOn2S = new IntSurf_LineOn2S();
4001     //
4002     for(i = 1; i <= nblines; i++) {
4003       if(anArrayOfLineType.Value(i) != 0) {
4004         continue;
4005       }
4006       const TColStd_ListOfInteger& aListOfIndex = anArrayOfLines.Value(i);
4007       const TColStd_ListOfInteger& aListOfFLIndex = anArrayOfLineEnds.Value(i);
4008       Standard_Boolean bhasfirstpoint = (aListOfFLIndex.Extent() == 2);
4009       Standard_Boolean bhaslastpoint = (aListOfFLIndex.Extent() == 2);
4010
4011       if(!bhasfirstpoint && !aListOfFLIndex.IsEmpty()) {
4012         bhasfirstpoint = (i != 1);
4013       }
4014
4015       if(!bhaslastpoint && !aListOfFLIndex.IsEmpty()) {
4016         bhaslastpoint = (i != nblines);
4017       }
4018
4019       Standard_Integer iFirst = aListOfIndex.First();
4020       Standard_Integer iLast  = aListOfIndex.Last();
4021       Standard_Boolean bIsFirstInside = ((ifprm >= iFirst) && (ifprm <= iLast));
4022       Standard_Boolean bIsLastInside =  ((ilprm >= iFirst) && (ilprm <= iLast));
4023
4024       if(!bIsFirstInside && !bIsLastInside) {
4025         if((ifprm < iFirst) && (ilprm > iLast)) {
4026           // append whole line, and boundaries if neccesary
4027           if(bhasfirstpoint) {
4028             pit = aListOfFLIndex.First();
4029             const IntSurf_PntOn2S& aP = aSeqOfPntOn2S->Value(pit);
4030             aLineOn2S->Add(aP);
4031           }
4032           TColStd_ListIteratorOfListOfInteger anIt(aListOfIndex);
4033
4034           for(; anIt.More(); anIt.Next()) {
4035             pit = anIt.Value();
4036             const IntSurf_PntOn2S& aP = theWLine->Point(pit);
4037             aLineOn2S->Add(aP);
4038           }
4039
4040           if(bhaslastpoint) {
4041             pit = aListOfFLIndex.Last();
4042             const IntSurf_PntOn2S& aP = aSeqOfPntOn2S->Value(pit);
4043             aLineOn2S->Add(aP);
4044           }
4045
4046           // check end of split line (end is almost always)
4047           Standard_Integer aneighbour = i + 1;
4048           Standard_Boolean bIsEndOfLine = Standard_True;
4049
4050           if(aneighbour <= nblines) {
4051             const TColStd_ListOfInteger& aListOfNeighbourIndex = anArrayOfLines.Value(aneighbour);
4052
4053             if((anArrayOfLineType.Value(aneighbour) != 0) &&
4054                (aListOfNeighbourIndex.IsEmpty())) {
4055               bIsEndOfLine = Standard_False;
4056             }
4057           }
4058
4059           if(bIsEndOfLine) {
4060             if(aLineOn2S->NbPoints() > 1) {
4061               Handle(IntPatch_WLine) aNewWLine =
4062                 new IntPatch_WLine(aLineOn2S, Standard_False);
4063               theNewLines.Append(aNewWLine);
4064             }
4065             aLineOn2S = new IntSurf_LineOn2S();
4066           }
4067         }
4068         continue;
4069       }
4070       // end if(!bIsFirstInside && !bIsLastInside)
4071
4072       if(bIsFirstInside && bIsLastInside) {
4073         // append inside points between ifprm and ilprm
4074         TColStd_ListIteratorOfListOfInteger anIt(aListOfIndex);
4075
4076         for(; anIt.More(); anIt.Next()) {
4077           pit = anIt.Value();
4078           if((pit < ifprm) || (pit > ilprm))
4079             continue;
4080           const IntSurf_PntOn2S& aP = theWLine->Point(pit);
4081           aLineOn2S->Add(aP);
4082         }
4083       }
4084       else {
4085
4086         if(bIsFirstInside) {
4087           // append points from ifprm to last point + boundary point
4088           TColStd_ListIteratorOfListOfInteger anIt(aListOfIndex);
4089
4090           for(; anIt.More(); anIt.Next()) {
4091             pit = anIt.Value();
4092             if(pit < ifprm)
4093               continue;
4094             const IntSurf_PntOn2S& aP = theWLine->Point(pit);
4095             aLineOn2S->Add(aP);
4096           }
4097
4098           if(bhaslastpoint) {
4099             pit = aListOfFLIndex.Last();
4100             const IntSurf_PntOn2S& aP = aSeqOfPntOn2S->Value(pit);
4101             aLineOn2S->Add(aP);
4102           }
4103           // check end of split line (end is almost always)
4104           Standard_Integer aneighbour = i + 1;
4105           Standard_Boolean bIsEndOfLine = Standard_True;
4106
4107           if(aneighbour <= nblines) {
4108             const TColStd_ListOfInteger& aListOfNeighbourIndex = anArrayOfLines.Value(aneighbour);
4109
4110             if((anArrayOfLineType.Value(aneighbour) != 0) &&
4111                (aListOfNeighbourIndex.IsEmpty())) {
4112               bIsEndOfLine = Standard_False;
4113             }
4114           }
4115
4116           if(bIsEndOfLine) {
4117             if(aLineOn2S->NbPoints() > 1) {
4118               Handle(IntPatch_WLine) aNewWLine =
4119                 new IntPatch_WLine(aLineOn2S, Standard_False);
4120               theNewLines.Append(aNewWLine);
4121             }
4122             aLineOn2S = new IntSurf_LineOn2S();
4123           }
4124         }
4125         // end if(bIsFirstInside)
4126
4127         if(bIsLastInside) {
4128           // append points from first boundary point to ilprm
4129           if(bhasfirstpoint) {
4130             pit = aListOfFLIndex.First();
4131             const IntSurf_PntOn2S& aP = aSeqOfPntOn2S->Value(pit);
4132             aLineOn2S->Add(aP);
4133           }
4134           TColStd_ListIteratorOfListOfInteger anIt(aListOfIndex);
4135
4136           for(; anIt.More(); anIt.Next()) {
4137             pit = anIt.Value();
4138             if(pit > ilprm)
4139               continue;
4140             const IntSurf_PntOn2S& aP = theWLine->Point(pit);
4141             aLineOn2S->Add(aP);
4142           }
4143         }
4144         //end if(bIsLastInside)
4145       }
4146     }
4147
4148     if(aLineOn2S->NbPoints() > 1) {
4149       Handle(IntPatch_WLine) aNewWLine =
4150         new IntPatch_WLine(aLineOn2S, Standard_False);
4151       theNewLines.Append(aNewWLine);
4152     }
4153   }
4154   // Split wlines.end
4155
4156   return Standard_True;
4157 }
4158
4159 // ------------------------------------------------------------------------------------------------
4160 // static function: ParameterOutOfBoundary
4161 // purpose:         Computes a new parameter for given curve. The corresponding 2d points
4162 //                  does not lay on any boundary of given faces
4163 // ------------------------------------------------------------------------------------------------
4164 Standard_Boolean ParameterOutOfBoundary(const Standard_Real       theParameter,
4165                                         const Handle(Geom_Curve)& theCurve,
4166                                         const TopoDS_Face&        theFace1,
4167                                         const TopoDS_Face&        theFace2,
4168                                         const Standard_Real       theOtherParameter,
4169                                         const Standard_Boolean    bIncreasePar,
4170                                         Standard_Real&            theNewParameter,
4171                                         const Handle(IntTools_Context)& aContext)
4172 {
4173   Standard_Boolean bIsComputed = Standard_False;
4174   theNewParameter = theParameter;
4175
4176   Standard_Real acurpar = theParameter;
4177   TopAbs_State aState = TopAbs_ON;
4178   Standard_Integer iter = 0;
4179   Standard_Real asumtol = BRep_Tool::Tolerance(theFace1) + BRep_Tool::Tolerance(theFace2);
4180   Standard_Real adelta = asumtol * 0.1;
4181   adelta = (adelta < Precision::Confusion()) ? Precision::Confusion() : adelta;
4182   Handle(Geom_Surface) aSurf1 = BRep_Tool::Surface(theFace1);
4183   Handle(Geom_Surface) aSurf2 = BRep_Tool::Surface(theFace2);
4184
4185   Standard_Real u1, u2, v1, v2;
4186
4187   GeomAPI_ProjectPointOnSurf aPrj1;
4188   aSurf1->Bounds(u1, u2, v1, v2);
4189   aPrj1.Init(aSurf1, u1, u2, v1, v2);
4190
4191   GeomAPI_ProjectPointOnSurf aPrj2;
4192   aSurf2->Bounds(u1, u2, v1, v2);
4193   aPrj2.Init(aSurf2, u1, u2, v1, v2);
4194
4195   while(aState == TopAbs_ON) {
4196     if(bIncreasePar)
4197       acurpar += adelta;
4198     else
4199       acurpar -= adelta;
4200     gp_Pnt aPCurrent = theCurve->Value(acurpar);
4201     aPrj1.Perform(aPCurrent);
4202     Standard_Real U=0., V=0.;
4203
4204     if(aPrj1.IsDone()) {
4205       aPrj1.LowerDistanceParameters(U, V);
4206       aState = aContext->StatePointFace(theFace1, gp_Pnt2d(U, V));
4207     }
4208
4209     if(aState != TopAbs_ON) {
4210       aPrj2.Perform(aPCurrent);
4211
4212       if(aPrj2.IsDone()) {
4213         aPrj2.LowerDistanceParameters(U, V);
4214         aState = aContext->StatePointFace(theFace2, gp_Pnt2d(U, V));
4215       }
4216     }
4217
4218     if(iter > 11) {
4219       break;
4220     }
4221     iter++;
4222   }
4223
4224   if(iter <= 11) {
4225     theNewParameter = acurpar;
4226     bIsComputed = Standard_True;
4227
4228     if(bIncreasePar) {
4229       if(acurpar >= theOtherParameter)
4230         theNewParameter = theOtherParameter;
4231     }
4232     else {
4233       if(acurpar <= theOtherParameter)
4234         theNewParameter = theOtherParameter;
4235     }
4236   }
4237   return bIsComputed;
4238 }
4239
4240 //=======================================================================
4241 //function : IsCurveValid
4242 //purpose  :
4243 //=======================================================================
4244 Standard_Boolean IsCurveValid(Handle(Geom2d_Curve)& thePCurve)
4245 {
4246   if(thePCurve.IsNull())
4247     return Standard_False;
4248
4249   Standard_Real tolint = 1.e-10;
4250   Geom2dAdaptor_Curve PCA;
4251   IntRes2d_Domain PCD;
4252   Geom2dInt_GInter PCI;
4253
4254   Standard_Real pf = 0., pl = 0.;
4255   gp_Pnt2d pntf, pntl;
4256
4257   if(!thePCurve->IsClosed() && !thePCurve->IsPeriodic()) {
4258     pf = thePCurve->FirstParameter();
4259     pl = thePCurve->LastParameter();
4260     pntf = thePCurve->Value(pf);
4261     pntl = thePCurve->Value(pl);
4262     PCA.Load(thePCurve);
4263     if(!PCA.IsPeriodic()) {
4264       if(PCA.FirstParameter() > pf) pf = PCA.FirstParameter();
4265       if(PCA.LastParameter()  < pl) pl = PCA.LastParameter();
4266     }
4267     PCD.SetValues(pntf,pf,tolint,pntl,pl,tolint);
4268     PCI.Perform(PCA,PCD,tolint,tolint);
4269     if(PCI.IsDone())
4270       if(PCI.NbPoints() > 0) {
4271         return Standard_False;
4272       }
4273   }
4274
4275   return Standard_True;
4276 }
4277
4278 //=======================================================================
4279 //static function : ApproxWithPCurves
4280 //purpose  : for bug 20964 only
4281 //=======================================================================
4282 Standard_Boolean ApproxWithPCurves(const gp_Cylinder& theCyl,
4283                                    const gp_Sphere& theSph)
4284 {
4285   Standard_Boolean bRes = Standard_True;
4286   Standard_Real R1 = theCyl.Radius(), R2 = theSph.Radius();
4287   //
4288   {
4289     Standard_Real aD2, aRc2, aEps;
4290     gp_Pnt aApexSph;
4291     //
4292     aEps=1.E-7;
4293     aRc2=R1*R1;
4294     //
4295     const gp_Ax3& aAx3Sph=theSph.Position();
4296     const gp_Pnt& aLocSph=aAx3Sph.Location();
4297     const gp_Dir& aDirSph=aAx3Sph.Direction();
4298     //
4299     const gp_Ax1& aAx1Cyl=theCyl.Axis();
4300     gp_Lin aLinCyl(aAx1Cyl);
4301     //
4302     aApexSph.SetXYZ(aLocSph.XYZ()+R2*aDirSph.XYZ());
4303     aD2=aLinCyl.SquareDistance(aApexSph);
4304     if (fabs(aD2-aRc2)<aEps) {
4305       return !bRes;
4306     }
4307     //
4308     aApexSph.SetXYZ(aLocSph.XYZ()-R2*aDirSph.XYZ());
4309     aD2=aLinCyl.SquareDistance(aApexSph);
4310     if (fabs(aD2-aRc2)<aEps) {
4311       return !bRes;
4312     }
4313   }
4314   //
4315
4316   if(R1 < 2.*R2) {
4317     return bRes;
4318   }
4319   gp_Lin anCylAx(theCyl.Axis());
4320
4321   Standard_Real aDist = anCylAx.Distance(theSph.Location());
4322   Standard_Real aDRel = Abs(aDist - R1)/R2;
4323
4324   if(aDRel > .2) return bRes;
4325
4326   Standard_Real par = ElCLib::Parameter(anCylAx, theSph.Location());
4327   gp_Pnt aP = ElCLib::Value(par, anCylAx);
4328   gp_Vec aV(aP, theSph.Location());
4329
4330   Standard_Real dd = aV.Dot(theSph.Position().XDirection());
4331
4332   if(aDist < R1 && dd > 0.) return Standard_False;
4333   if(aDist > R1 && dd < 0.) return Standard_False;
4334
4335
4336   return bRes;
4337 }
4338 //=======================================================================
4339 //function : PerformPlanes
4340 //purpose  :
4341 //=======================================================================
4342 void  PerformPlanes(const Handle(GeomAdaptor_HSurface)& theS1,
4343                     const Handle(GeomAdaptor_HSurface)& theS2,
4344                     const Standard_Real TolAng,
4345                     const Standard_Real TolTang,
4346                     const Standard_Boolean theApprox1,
4347                     const Standard_Boolean theApprox2,
4348                     IntTools_SequenceOfCurves& theSeqOfCurve,
4349                     Standard_Boolean& theTangentFaces)
4350 {
4351
4352   gp_Pln aPln1 = theS1->Surface().Plane();
4353   gp_Pln aPln2 = theS2->Surface().Plane();
4354
4355   IntAna_QuadQuadGeo aPlnInter(aPln1, aPln2, TolAng, TolTang);
4356
4357   if(!aPlnInter.IsDone()) {
4358     theTangentFaces = Standard_False;
4359     return;
4360   }
4361
4362   IntAna_ResultType aResType = aPlnInter.TypeInter();
4363
4364   if(aResType == IntAna_Same) {
4365     theTangentFaces = Standard_True;
4366     return;
4367   }
4368
4369   theTangentFaces = Standard_False;
4370
4371   if(aResType == IntAna_Empty) {
4372     return;
4373   }
4374
4375   gp_Lin aLin = aPlnInter.Line(1);
4376
4377   ProjLib_Plane aProj;
4378
4379   aProj.Init(aPln1);
4380   aProj.Project(aLin);
4381   gp_Lin2d aLin2d1 = aProj.Line();
4382   //
4383   aProj.Init(aPln2);
4384   aProj.Project(aLin);
4385   gp_Lin2d aLin2d2 = aProj.Line();
4386   //
4387   //classify line2d1 relatively first plane
4388   Standard_Real P11, P12;
4389   Standard_Boolean IsCrossed = ClassifyLin2d(theS1, aLin2d1, TolTang, P11, P12);
4390   if(!IsCrossed) return;
4391   //classify line2d2 relatively second plane
4392   Standard_Real P21, P22;
4393   IsCrossed = ClassifyLin2d(theS2, aLin2d2, TolTang, P21, P22);
4394   if(!IsCrossed) return;
4395
4396   //Analysis of parametric intervals: must have common part
4397
4398   if(P21 >= P12) return;
4399   if(P22 <= P11) return;
4400
4401   Standard_Real pmin, pmax;
4402   pmin = Max(P11, P21);
4403   pmax = Min(P12, P22);
4404
4405   if(pmax - pmin <= TolTang) return;
4406
4407   Handle(Geom_Line) aGLin = new Geom_Line(aLin);
4408
4409   IntTools_Curve aCurve;
4410   Handle(Geom_TrimmedCurve) aGTLin = new Geom_TrimmedCurve(aGLin, pmin, pmax);
4411
4412   aCurve.SetCurve(aGTLin);
4413
4414   if(theApprox1) {
4415     Handle(Geom2d_Line) C2d = new Geom2d_Line(aLin2d1);
4416     aCurve.SetFirstCurve2d(new Geom2d_TrimmedCurve(C2d, pmin, pmax));
4417   }
4418   else {
4419     Handle(Geom2d_Curve) H1;
4420     aCurve.SetFirstCurve2d(H1);
4421   }
4422   if(theApprox2) {
4423     Handle(Geom2d_Line) C2d = new Geom2d_Line(aLin2d2);
4424     aCurve.SetSecondCurve2d(new Geom2d_TrimmedCurve(C2d, pmin, pmax));
4425   }
4426   else {
4427     Handle(Geom2d_Curve) H1;
4428     aCurve.SetFirstCurve2d(H1);
4429   }
4430
4431   theSeqOfCurve.Append(aCurve);
4432
4433 }
4434
4435 //=======================================================================
4436 //function : ClassifyLin2d
4437 //purpose  :
4438 //=======================================================================
4439 static inline Standard_Boolean INTER(const Standard_Real d1,
4440                                      const Standard_Real d2,
4441                                      const Standard_Real tol)
4442 {
4443   return (d1 > tol && d2 < -tol) ||
4444          (d1 < -tol && d2 > tol) ||
4445          ((d1 <= tol && d1 >= -tol) && (d2 > tol || d2 < -tol)) ||
4446          ((d2 <= tol && d2 >= -tol) && (d1 > tol || d1 < -tol));
4447 }
4448 static inline Standard_Boolean COINC(const Standard_Real d1,
4449                                      const Standard_Real d2,
4450                                      const Standard_Real tol)
4451 {
4452   return (d1 <= tol && d1 >= -tol) && (d2 <= tol && d2 >= -tol);
4453 }
4454 Standard_Boolean ClassifyLin2d(const Handle(GeomAdaptor_HSurface)& theS,
4455                                const gp_Lin2d& theLin2d,
4456                                const Standard_Real theTol,
4457                                Standard_Real& theP1,
4458                                Standard_Real& theP2)
4459
4460 {
4461   Standard_Real xmin, xmax, ymin, ymax, d1, d2, A, B, C;
4462   Standard_Real par[2];
4463   Standard_Integer nbi = 0;
4464
4465   xmin = theS->Surface().FirstUParameter();
4466   xmax = theS->Surface().LastUParameter();
4467   ymin = theS->Surface().FirstVParameter();
4468   ymax = theS->Surface().LastVParameter();
4469
4470   theLin2d.Coefficients(A, B, C);
4471
4472   //xmin, ymin <-> xmin, ymax
4473   d1 = A*xmin + B*ymin + C;
4474   d2 = A*xmin + B*ymax + C;
4475
4476   if(INTER(d1, d2, theTol)) {
4477     //Intersection with boundary
4478     Standard_Real y = -(C + A*xmin)/B;
4479     par[nbi] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmin, y));
4480     nbi++;
4481   }
4482   else if (COINC(d1, d2, theTol)) {
4483     //Coincidence with boundary
4484     par[0] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmin, ymin));
4485     par[1] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmin, ymax));
4486     nbi = 2;
4487   }
4488
4489   if(nbi == 2) {
4490
4491     if(fabs(par[0]-par[1]) > theTol) {
4492       theP1 = Min(par[0], par[1]);
4493       theP2 = Max(par[0], par[1]);
4494       return Standard_True;
4495     }
4496     else return Standard_False;
4497
4498   }
4499
4500   //xmin, ymax <-> xmax, ymax
4501   d1 = d2;
4502   d2 = A*xmax + B*ymax + C;
4503
4504   if(d1 > theTol || d1 < -theTol) {//to avoid checking of
4505                                    //coincidence with the same point
4506     if(INTER(d1, d2, theTol)) {
4507       Standard_Real x = -(C + B*ymax)/A;
4508       par[nbi] = ElCLib::Parameter(theLin2d, gp_Pnt2d(x, ymax));
4509       nbi++;
4510     }
4511     else if (COINC(d1, d2, theTol)) {
4512       par[0] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmin, ymax));
4513       par[1] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmax, ymax));
4514       nbi = 2;
4515     }
4516   }
4517
4518   if(nbi == 2) {
4519
4520     if(fabs(par[0]-par[1]) > theTol) {
4521       theP1 = Min(par[0], par[1]);
4522       theP2 = Max(par[0], par[1]);
4523       return Standard_True;
4524     }
4525     else return Standard_False;
4526
4527   }
4528
4529   //xmax, ymax <-> xmax, ymin
4530   d1 = d2;
4531   d2 = A*xmax + B*ymin + C;
4532
4533   if(d1 > theTol || d1 < -theTol) {
4534     if(INTER(d1, d2, theTol)) {
4535       Standard_Real y = -(C + A*xmax)/B;
4536       par[nbi] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmax, y));
4537       nbi++;
4538     }
4539     else if (COINC(d1, d2, theTol)) {
4540       par[0] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmax, ymax));
4541       par[1] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmax, ymin));
4542       nbi = 2;
4543     }
4544   }
4545
4546   if(nbi == 2) {
4547     if(fabs(par[0]-par[1]) > theTol) {
4548       theP1 = Min(par[0], par[1]);
4549       theP2 = Max(par[0], par[1]);
4550       return Standard_True;
4551     }
4552     else return Standard_False;
4553   }
4554
4555   //xmax, ymin <-> xmin, ymin
4556   d1 = d2;
4557   d2 = A*xmin + B*ymin + C;
4558
4559   if(d1 > theTol || d1 < -theTol) {
4560     if(INTER(d1, d2, theTol)) {
4561       Standard_Real x = -(C + B*ymin)/A;
4562       par[nbi] = ElCLib::Parameter(theLin2d, gp_Pnt2d(x, ymin));
4563       nbi++;
4564     }
4565     else if (COINC(d1, d2, theTol)) {
4566       par[0] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmax, ymin));
4567       par[1] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmin, ymin));
4568       nbi = 2;
4569     }
4570   }
4571
4572   if(nbi == 2) {
4573     if(fabs(par[0]-par[1]) > theTol) {
4574       theP1 = Min(par[0], par[1]);
4575       theP2 = Max(par[0], par[1]);
4576       return Standard_True;
4577     }
4578     else return Standard_False;
4579   }
4580
4581   return Standard_False;
4582
4583 }
4584 //
4585 //=======================================================================
4586 //function : ApproxParameters
4587 //purpose  :
4588 //=======================================================================
4589 void ApproxParameters(const Handle(GeomAdaptor_HSurface)& aHS1,
4590                       const Handle(GeomAdaptor_HSurface)& aHS2,
4591                       Standard_Integer& iDegMin,
4592                       Standard_Integer& iDegMax,
4593                       Standard_Integer& iNbIter)
4594
4595 {
4596   GeomAbs_SurfaceType aTS1, aTS2;
4597
4598   //
4599   iNbIter=0;
4600   iDegMin=4;
4601   iDegMax=8;
4602   //
4603   aTS1=aHS1->Surface().GetType();
4604   aTS2=aHS2->Surface().GetType();
4605   //
4606   // Cylinder/Torus
4607   if ((aTS1==GeomAbs_Cylinder && aTS2==GeomAbs_Torus) ||
4608       (aTS2==GeomAbs_Cylinder && aTS1==GeomAbs_Torus)) {
4609     Standard_Real aRC, aRT, dR, aPC;
4610     gp_Cylinder aCylinder;
4611     gp_Torus aTorus;
4612     //
4613     aPC=Precision::Confusion();
4614     //
4615     aCylinder=(aTS1==GeomAbs_Cylinder)? aHS1->Surface().Cylinder() : aHS2->Surface().Cylinder();
4616     aTorus=(aTS1==GeomAbs_Torus)? aHS1->Surface().Torus() : aHS2->Surface().Torus();
4617     //
4618     aRC=aCylinder.Radius();
4619     aRT=aTorus.MinorRadius();
4620     dR=aRC-aRT;
4621     if (dR<0.) {
4622       dR=-dR;
4623     }
4624     //
4625     if (dR<aPC) {
4626       iDegMax=6;
4627     }
4628   }
4629   if (aTS1==GeomAbs_Cylinder && aTS2==GeomAbs_Cylinder) {
4630     iNbIter=1;
4631   }
4632 }
4633 //=======================================================================
4634 //function : Tolerances
4635 //purpose  :
4636 //=======================================================================
4637 void Tolerances(const Handle(GeomAdaptor_HSurface)& aHS1,
4638                 const Handle(GeomAdaptor_HSurface)& aHS2,
4639                 Standard_Real& aTolTang)
4640 {
4641   GeomAbs_SurfaceType aTS1, aTS2;
4642   //
4643   aTS1=aHS1->Surface().GetType();
4644   aTS2=aHS2->Surface().GetType();
4645   //
4646   // Cylinder/Torus
4647   if ((aTS1==GeomAbs_Cylinder && aTS2==GeomAbs_Torus) ||
4648       (aTS2==GeomAbs_Cylinder && aTS1==GeomAbs_Torus)) {
4649     Standard_Real aRC, aRT, dR, aPC;
4650     gp_Cylinder aCylinder;
4651     gp_Torus aTorus;
4652     //
4653     aPC=Precision::Confusion();
4654     //
4655     aCylinder=(aTS1==GeomAbs_Cylinder)? aHS1->Surface().Cylinder() : aHS2->Surface().Cylinder();
4656     aTorus=(aTS1==GeomAbs_Torus)? aHS1->Surface().Torus() : aHS2->Surface().Torus();
4657     //
4658     aRC=aCylinder.Radius();
4659     aRT=aTorus.MinorRadius();
4660     dR=aRC-aRT;
4661     if (dR<0.) {
4662       dR=-dR;
4663     }
4664     //
4665     if (dR<aPC) {
4666       aTolTang=0.1*aTolTang;
4667     }
4668   }
4669 }
4670 //=======================================================================
4671 //function : SortTypes
4672 //purpose  :
4673 //=======================================================================
4674 Standard_Boolean SortTypes(const GeomAbs_SurfaceType aType1,
4675                            const GeomAbs_SurfaceType aType2)
4676 {
4677   Standard_Boolean bRet;
4678   Standard_Integer aI1, aI2;
4679   //
4680   bRet=Standard_False;
4681   //
4682   aI1=IndexType(aType1);
4683   aI2=IndexType(aType2);
4684   if (aI1<aI2){
4685     bRet=!bRet;
4686   }
4687   return bRet;
4688 }
4689 //=======================================================================
4690 //function : IndexType
4691 //purpose  :
4692 //=======================================================================
4693 Standard_Integer IndexType(const GeomAbs_SurfaceType aType)
4694 {
4695   Standard_Integer aIndex;
4696   //
4697   aIndex=11;
4698   //
4699   if (aType==GeomAbs_Plane) {
4700     aIndex=0;
4701   }
4702   else if (aType==GeomAbs_Cylinder) {
4703     aIndex=1;
4704   }
4705   else if (aType==GeomAbs_Cone) {
4706     aIndex=2;
4707   }
4708   else if (aType==GeomAbs_Sphere) {
4709     aIndex=3;
4710   }
4711   else if (aType==GeomAbs_Torus) {
4712     aIndex=4;
4713   }
4714   else if (aType==GeomAbs_BezierSurface) {
4715     aIndex=5;
4716   }
4717   else if (aType==GeomAbs_BSplineSurface) {
4718     aIndex=6;
4719   }
4720   else if (aType==GeomAbs_SurfaceOfRevolution) {
4721     aIndex=7;
4722   }
4723   else if (aType==GeomAbs_SurfaceOfExtrusion) {
4724     aIndex=8;
4725   }
4726   else if (aType==GeomAbs_OffsetSurface) {
4727     aIndex=9;
4728   }
4729   else if (aType==GeomAbs_OtherSurface) {
4730     aIndex=10;
4731   }
4732   return aIndex;
4733 }
4734 #ifdef OCCT_DEBUG_DUMPWLINE
4735 //=======================================================================
4736 //function : DumpWLine
4737 //purpose  :
4738 //=======================================================================
4739 void DumpWLine(const Handle(IntPatch_WLine)& aWLine)
4740 {
4741   Standard_Integer i, aNbPnts;
4742   Standard_Real aX, aY, aZ, aU1, aV1, aU2, aV2;
4743   //
4744   printf(" *WLine\n");
4745   aNbPnts=aWLine->NbPnts();
4746   for (i=1; i<=aNbPnts; ++i) {
4747     const IntSurf_PntOn2S aPntOn2S=aWLine->Point(i);
4748     const gp_Pnt& aP3D=aPntOn2S.Value();
4749     aP3D.Coord(aX, aY, aZ);
4750     aPntOn2S.Parameters(aU1, aV1, aU2, aV2);
4751     //
4752     printf("point p_%d %lf %lf %lf\n", i, aX, aY, aZ);
4753     //printf("point p_%d %20.15lf %20.15lf %20.15lf %20.15lf %20.15lf %20.15lf %20.15lf\n",
4754         //   i, aX, aY, aZ, aU1, aV1, aU2, aV2);
4755   }
4756 }
4757 #endif
4758 //=======================================================================
4759 //function : RefineVector
4760 //purpose  :
4761 //=======================================================================
4762 void RefineVector(gp_Vec2d& aV2D)
4763 {
4764   Standard_Integer k,m;
4765   Standard_Real aC[2], aEps, aR1, aR2, aNum;
4766   //
4767   aEps=RealEpsilon();
4768   aR1=1.-aEps;
4769   aR2=1.+aEps;
4770   //
4771   aV2D.Coord(aC[0], aC[1]);
4772   //
4773   for (k=0; k<2; ++k) {
4774     m=(k+1)%2;
4775     aNum=fabs(aC[k]);
4776     if (aNum>aR1 && aNum<aR2) {
4777       if (aC[k]<0.) {
4778         aC[k]=-1.;
4779       }
4780       else {
4781         aC[k]=1.;
4782       }
4783       aC[m]=0.;
4784       break;
4785     }
4786   }
4787   aV2D.SetCoord(aC[0], aC[1]);
4788 }
4789 //=======================================================================
4790 // Function : FindMaxDistance
4791 // purpose :
4792 //=======================================================================
4793 Standard_Real FindMaxDistance(const Handle(Geom_Curve)& theCurve,
4794                               const Standard_Real theFirst,
4795                               const Standard_Real theLast,
4796                               const TopoDS_Face& theFace,
4797                               const Handle(IntTools_Context)& theContext)
4798 {
4799   Standard_Integer aNbS;
4800   Standard_Real aT1, aT2, aDt, aD, aDMax, anEps;
4801   //
4802   aNbS = 11;
4803   aDt = (theLast - theFirst) / aNbS;
4804   aDMax = 0.;
4805   anEps = 1.e-4 * aDt;
4806   //
4807   GeomAPI_ProjectPointOnSurf& aProjPS = theContext->ProjPS(theFace);
4808   aT2 = theFirst;
4809   for (;;) {
4810     aT1 = aT2;
4811     aT2 += aDt;
4812     //
4813     if (aT2 > theLast) {
4814       break;
4815     }
4816     //
4817     aD = FindMaxDistance(theCurve, aT1, aT2, aProjPS, anEps);
4818     if (aD > aDMax) {
4819       aDMax = aD;
4820     }
4821   }
4822   //
4823   return aDMax;
4824 }
4825 //=======================================================================
4826 // Function : FindMaxDistance
4827 // purpose :
4828 //=======================================================================
4829 Standard_Real FindMaxDistance(const Handle(Geom_Curve)& theC,
4830                               const Standard_Real theFirst,
4831                               const Standard_Real theLast,
4832                               GeomAPI_ProjectPointOnSurf& theProjPS,
4833                               const Standard_Real theEps)
4834 {
4835   Standard_Real aA, aB, aCf, aX, aX1, aX2, aF1, aF2, aF;
4836   //
4837   aCf = 0.61803398874989484820458683436564;//(sqrt(5.)-1)/2.;
4838   aA = theFirst;
4839   aB = theLast;
4840   //
4841   aX1 = aB - aCf * (aB - aA);
4842   aF1 = MaxDistance(theC, aX1, theProjPS);
4843   aX2 = aA + aCf * (aB - aA);
4844   aF2 = MaxDistance(theC, aX2, theProjPS);
4845   //
4846   for (;;) {
4847     if ((aB - aA) < theEps) {
4848       break;
4849     }
4850     //
4851     if (aF1 > aF2) {
4852       aB = aX2;
4853       aX2 = aX1;
4854       aF2 = aF1;
4855       aX1 = aB - aCf * (aB - aA);
4856       aF1 = MaxDistance(theC, aX1, theProjPS);
4857     }
4858     else {
4859       aA = aX1;
4860       aX1 = aX2;
4861       aF1 = aF2;
4862       aX2 = aA + aCf * (aB - aA);
4863       aF2 = MaxDistance(theC, aX2, theProjPS);
4864     }
4865   }
4866   //
4867   aX = 0.5 * (aA + aB);
4868   aF = MaxDistance(theC, aX, theProjPS);
4869   //
4870   if (aF1 > aF) {
4871     aF = aF1;
4872   }
4873   //
4874   if (aF2 > aF) {
4875     aF = aF2;
4876   }
4877   //
4878   return aF;
4879 }
4880 //=======================================================================
4881 // Function : MaxDistance
4882 // purpose :
4883 //=======================================================================
4884 Standard_Real MaxDistance(const Handle(Geom_Curve)& theC,
4885                           const Standard_Real aT,
4886                           GeomAPI_ProjectPointOnSurf& theProjPS)
4887 {
4888   Standard_Real aD;
4889   gp_Pnt aP;
4890   //
4891   theC->D0(aT, aP);
4892   theProjPS.Perform(aP);
4893   aD = theProjPS.NbPoints() ? theProjPS.LowerDistance() : 0.;
4894   //
4895   return aD;
4896 }
4897
4898 //=======================================================================
4899 //function : CheckPCurve
4900 //purpose  : Checks if points of the pcurve are out of the face bounds.
4901 //=======================================================================
4902   Standard_Boolean CheckPCurve(const Handle(Geom2d_Curve)& aPC,
4903                                const TopoDS_Face& aFace)
4904 {
4905   const Standard_Integer NPoints = 23;
4906   Standard_Integer i;
4907   Standard_Real umin,umax,vmin,vmax;
4908
4909   BRepTools::UVBounds(aFace, umin, umax, vmin, vmax);
4910   Standard_Real tolU = Max ((umax-umin)*0.01, Precision::Confusion());
4911   Standard_Real tolV = Max ((vmax-vmin)*0.01, Precision::Confusion());
4912   Standard_Real fp = aPC->FirstParameter();
4913   Standard_Real lp = aPC->LastParameter();
4914
4915
4916   // adjust domain for periodic surfaces
4917   TopLoc_Location aLoc;
4918   Handle(Geom_Surface) aSurf = BRep_Tool::Surface(aFace, aLoc);
4919   if (aSurf->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface))) {
4920     aSurf = (Handle(Geom_RectangularTrimmedSurface)::DownCast(aSurf))->BasisSurface();
4921   }
4922   gp_Pnt2d pnt = aPC->Value((fp+lp)/2);
4923   Standard_Real u,v;
4924   pnt.Coord(u,v);
4925   //
4926   if (aSurf->IsUPeriodic()) {
4927     Standard_Real aPer = aSurf->UPeriod();
4928     Standard_Integer nshift = (Standard_Integer) ((u-umin)/aPer);
4929     if (u < umin+aPer*nshift) nshift--;
4930     umin += aPer*nshift;
4931     umax += aPer*nshift;
4932   }
4933   if (aSurf->IsVPeriodic()) {
4934     Standard_Real aPer = aSurf->VPeriod();
4935     Standard_Integer nshift = (Standard_Integer) ((v-vmin)/aPer);
4936     if (v < vmin+aPer*nshift) nshift--;
4937     vmin += aPer*nshift;
4938     vmax += aPer*nshift;
4939   }
4940   //
4941   //--------------------------------------------------------
4942   Standard_Boolean bRet;
4943   Standard_Integer j, aNbIntervals;
4944   Standard_Real aT, dT;
4945   gp_Pnt2d aP2D;
4946   //
4947   Geom2dAdaptor_Curve aGAC(aPC);
4948   aNbIntervals=aGAC.NbIntervals(GeomAbs_CN);
4949   //
4950   TColStd_Array1OfReal aTI(1, aNbIntervals+1);
4951   aGAC.Intervals(aTI,GeomAbs_CN);
4952   //
4953   bRet=Standard_False;
4954   //
4955   aT=aGAC.FirstParameter();
4956   for (j=1; j<=aNbIntervals; ++j) {
4957     dT=(aTI(j+1)-aTI(j))/NPoints;
4958     //
4959     for (i=1; i<NPoints; i++) {
4960       aT=aT+dT;
4961       aGAC.D0(aT, aP2D);
4962       aP2D.Coord(u,v);
4963     if (umin-u > tolU || u-umax > tolU ||
4964           vmin-v > tolV || v-vmax > tolV) {
4965         return bRet;
4966   }
4967 }
4968   }
4969   return !bRet;
4970 }
4971 //=======================================================================
4972 //function : CorrectPlaneBoundaries
4973 //purpose  :
4974 //=======================================================================
4975  void CorrectPlaneBoundaries(Standard_Real& aUmin,
4976                              Standard_Real& aUmax,
4977                              Standard_Real& aVmin,
4978                              Standard_Real& aVmax)
4979 {
4980   if (!(Precision::IsInfinite(aUmin) ||
4981         Precision::IsInfinite(aUmax))) {
4982     Standard_Real dU;
4983     //
4984     dU=0.1*(aUmax-aUmin);
4985     aUmin=aUmin-dU;
4986     aUmax=aUmax+dU;
4987   }
4988   if (!(Precision::IsInfinite(aVmin) ||
4989         Precision::IsInfinite(aVmax))) {
4990     Standard_Real dV;
4991     //
4992     dV=0.1*(aVmax-aVmin);
4993     aVmin=aVmin-dV;
4994     aVmax=aVmax+dV;
4995   }
4996 }