src/IntTools/IntTools_FaceFace.cxx

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