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