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