src/IntTools/IntTools_FaceFace.cxx

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