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