src/IntTools/IntTools_FaceFace.cxx

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