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