src/IntTools/IntTools_FaceFace.cxx

   1 // Created on: 2000-11-23
   2 // Created by: Michael KLOKOV
   3 // Copyright (c) 2000-2014 OPEN CASCADE SAS
   4 //
   5 // This file is part of Open CASCADE Technology software library.
   6 //
   7 // This library is free software; you can redistribute it and/or modify it under
   8 // the terms of the GNU Lesser General Public License version 2.1 as published
   9 // by the Free Software Foundation, with special exception defined in the file
  10 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
  11 // distribution for complete text of the license and disclaimer of any warranty.
  12 //
  13 // Alternatively, this file may be used under the terms of Open CASCADE
  14 // commercial license or contractual agreement.
  15
  16 #include <IntTools_FaceFace.hxx>
  17
  18 #include <Precision.hxx>
  19
  20 #include <TColStd_HArray1OfReal.hxx>
  21 #include <TColStd_Array1OfReal.hxx>
  22 #include <TColStd_Array1OfInteger.hxx>
  23 #include <TColStd_SequenceOfReal.hxx>
  24 #include <TColStd_ListOfInteger.hxx>
  25 #include <TColStd_ListIteratorOfListOfInteger.hxx>
  26 #include <TColStd_Array1OfListOfInteger.hxx>
  27
  28 #include <Bnd_Box.hxx>
  29
  30 #include <TColgp_HArray1OfPnt2d.hxx>
  31 #include <TColgp_SequenceOfPnt2d.hxx>
  32 #include <TColgp_Array1OfPnt.hxx>
  33 #include <TColgp_Array1OfPnt2d.hxx>
  34
  35 #include <IntAna_QuadQuadGeo.hxx>
  36
  37 #include <IntSurf_PntOn2S.hxx>
  38 #include <IntSurf_LineOn2S.hxx>
  39 #include <IntSurf_PntOn2S.hxx>
  40 #include <IntSurf_ListOfPntOn2S.hxx>
  41 #include <IntRes2d_Domain.hxx>
  42 #include <ProjLib_Plane.hxx>
  43
  44 #include <IntPatch_GLine.hxx>
  45 #include <IntPatch_RLine.hxx>
  46 #include <IntPatch_WLine.hxx>
  47 #include <IntPatch_ALine.hxx>
  48 #include <IntPatch_ALineToWLine.hxx>
  49
  50 #include <ElSLib.hxx>
  51 #include <ElCLib.hxx>
  52
  53 #include <BndLib_AddSurface.hxx>
  54
  55 #include <Adaptor3d_SurfacePtr.hxx>
  56 #include <Adaptor2d_HLine2d.hxx>
  57
  58 #include <GeomAbs_SurfaceType.hxx>
  59 #include <GeomAbs_CurveType.hxx>
  60
  61 #include <Geom_Surface.hxx>
  62 #include <Geom_Line.hxx>
  63 #include <Geom_Circle.hxx>
  64 #include <Geom_Ellipse.hxx>
  65 #include <Geom_Parabola.hxx>
  66 #include <Geom_Hyperbola.hxx>
  67 #include <Geom_TrimmedCurve.hxx>
  68 #include <Geom_BSplineCurve.hxx>
  69 #include <Geom_RectangularTrimmedSurface.hxx>
  70 #include <Geom_OffsetSurface.hxx>
  71 #include <Geom_Curve.hxx>
  72 #include <Geom_Conic.hxx>
  73
  74 #include <Geom2d_TrimmedCurve.hxx>
  75 #include <Geom2d_BSplineCurve.hxx>
  76 #include <Geom2d_Line.hxx>
  77 #include <Geom2d_Curve.hxx>
  78
  79 #include <Geom2dInt_GInter.hxx>
  80 #include <Geom2dAdaptor.hxx>
  81 #include <GeomAdaptor_HSurface.hxx>
  82 #include <GeomAdaptor_Surface.hxx>
  83 #include <GeomLib_CheckBSplineCurve.hxx>
  84 #include <GeomLib_Check2dBSplineCurve.hxx>
  85
  86 #include <GeomInt_WLApprox.hxx>
  87 #include <GeomProjLib.hxx>
  88 #include <GeomAPI_ProjectPointOnSurf.hxx>
  89 #include <Geom2dAdaptor_Curve.hxx>
  90 #include <TopoDS.hxx>
  91 #include <TopoDS_Edge.hxx>
  92 #include <TopExp_Explorer.hxx>
  93
  94 #include <BRep_Tool.hxx>
  95 #include <BRepTools.hxx>
  96 #include <BRepAdaptor_Surface.hxx>
  97
  98 #include <IntTools_Curve.hxx>
  99 #include <IntTools_Tools.hxx>
 100 #include <IntTools_Tools.hxx>
 101 #include <IntTools_TopolTool.hxx>
 102 #include <IntTools_PntOnFace.hxx>
 103 #include <IntTools_PntOn2Faces.hxx>
 104 #include <IntTools_Context.hxx>
 105 #include <IntSurf_ListIteratorOfListOfPntOn2S.hxx>
 106
 107 #include <Approx_CurveOnSurface.hxx>
 108 #include <GeomAdaptor.hxx>
 109 #include <GeomInt_IntSS.hxx>
 110 #include <IntTools_WLineTool.hxx>
 111 #include <IntPatch_WLineTool.hxx>
 112
 113 //#ifdef OCCT_DEBUG_DUMPWLINE
 114 //static
 115 //  void DumpWLine(const Handle(IntPatch_WLine)& aWLine);
 116 //#endif
 117 ////
 118 static
 119   void TolR3d(const TopoDS_Face& ,
 120               const TopoDS_Face& ,
 121               Standard_Real& );
 122
 123 static
 124   void Parameters(const Handle(GeomAdaptor_HSurface)&,
 125                   const Handle(GeomAdaptor_HSurface)&,
 126                   const gp_Pnt&,
 127                   Standard_Real&,
 128                   Standard_Real&,
 129                   Standard_Real&,
 130                   Standard_Real&);
 131
 132 static
 133   void CorrectSurfaceBoundaries(const TopoDS_Face&  theFace,
 134                                 const Standard_Real theTolerance,
 135                                 Standard_Real&      theumin,
 136                                 Standard_Real&      theumax,
 137                                 Standard_Real&      thevmin,
 138                                 Standard_Real&      thevmax);
 139
 140 static
 141   Standard_Boolean ParameterOutOfBoundary(const Standard_Real       theParameter,
 142                                           const Handle(Geom_Curve)& theCurve,
 143                                           const TopoDS_Face&        theFace1,
 144                                           const TopoDS_Face&        theFace2,
 145                                           const Standard_Real       theOtherParameter,
 146                                           const Standard_Boolean    bIncreasePar,
 147                                           Standard_Real&            theNewParameter,
 148                                           const Handle(IntTools_Context)& );
 149
 150 static
 151   Standard_Boolean IsCurveValid(Handle(Geom2d_Curve)& thePCurve);
 152
 153 static
 154   Standard_Boolean  ApproxWithPCurves(const gp_Cylinder& theCyl,
 155                                       const gp_Sphere& theSph);
 156
 157 static void  PerformPlanes(const Handle(GeomAdaptor_HSurface)& theS1,
 158                            const Handle(GeomAdaptor_HSurface)& theS2,
 159                            const Standard_Real TolAng,
 160                            const Standard_Real TolTang,
 161                            const Standard_Boolean theApprox1,
 162                            const Standard_Boolean theApprox2,
 163                            IntTools_SequenceOfCurves& theSeqOfCurve,
 164                            Standard_Boolean& theTangentFaces);
 165
 166 static Standard_Boolean ClassifyLin2d(const Handle(GeomAdaptor_HSurface)& theS,
 167                                       const gp_Lin2d& theLin2d,
 168                                       const Standard_Real theTol,
 169                                       Standard_Real& theP1,
 170                                       Standard_Real& theP2);
 171 //
 172 static
 173   void ApproxParameters(const Handle(GeomAdaptor_HSurface)& aHS1,
 174                         const Handle(GeomAdaptor_HSurface)& aHS2,
 175                         Standard_Integer& iDegMin,
 176                         Standard_Integer& iNbIter,
 177                         Standard_Integer& iDegMax);
 178
 179 static
 180   void Tolerances(const Handle(GeomAdaptor_HSurface)& aHS1,
 181                   const Handle(GeomAdaptor_HSurface)& aHS2,
 182                   Standard_Real& aTolTang);
 183
 184 static
 185   Standard_Boolean SortTypes(const GeomAbs_SurfaceType aType1,
 186                              const GeomAbs_SurfaceType aType2);
 187 static
 188   Standard_Integer IndexType(const GeomAbs_SurfaceType aType);
 189
 190 //
 191 static
 192   Standard_Boolean CheckPCurve(const Handle(Geom2d_Curve)& aPC,
 193                                const TopoDS_Face& aFace);
 194
 195 static
 196   Standard_Real MaxDistance(const Handle(Geom_Curve)& theC,
 197                             const Standard_Real aT,
 198                             GeomAPI_ProjectPointOnSurf& theProjPS);
 199
 200 static
 201   Standard_Real FindMaxDistance(const Handle(Geom_Curve)& theC,
 202                                 const Standard_Real theFirst,
 203                                 const Standard_Real theLast,
 204                                 GeomAPI_ProjectPointOnSurf& theProjPS,
 205                                 const Standard_Real theEps);
 206
 207 static
 208   Standard_Real FindMaxDistance(const Handle(Geom_Curve)& theCurve,
 209                                 const Standard_Real theFirst,
 210                                 const Standard_Real theLast,
 211                                 const TopoDS_Face& theFace,
 212                                 const Handle(IntTools_Context)& theContext);
 213
 214 static
 215   void CorrectPlaneBoundaries(Standard_Real& aUmin,
 216                               Standard_Real& aUmax,
 217                               Standard_Real& aVmin,
 218                               Standard_Real& aVmax);
 219
 220 //=======================================================================
 221 //function :
 222 //purpose  :
 223 //=======================================================================
 224 IntTools_FaceFace::IntTools_FaceFace()
 225 {
 226   myIsDone=Standard_False;
 227   myTangentFaces=Standard_False;
 228   //
 229   myHS1 = new GeomAdaptor_HSurface ();
 230   myHS2 = new GeomAdaptor_HSurface ();
 231   myTolReached2d=0.;
 232   myTolReached3d=0.;
 233   SetParameters(Standard_True, Standard_True, Standard_True, 1.e-07);
 234
 235 }
 236 //=======================================================================
 237 //function : SetContext
 238 //purpose  :
 239 //=======================================================================
 240 void IntTools_FaceFace::SetContext(const Handle(IntTools_Context)& aContext)
 241 {
 242   myContext=aContext;
 243 }
 244 //=======================================================================
 245 //function : Context
 246 //purpose  :
 247 //=======================================================================
 248 const Handle(IntTools_Context)& IntTools_FaceFace::Context()const
 249 {
 250   return myContext;
 251 }
 252 //=======================================================================
 253 //function : Face1
 254 //purpose  :
 255 //=======================================================================
 256 const TopoDS_Face& IntTools_FaceFace::Face1() const
 257 {
 258   return myFace1;
 259 }
 260 //=======================================================================
 261 //function : Face2
 262 //purpose  :
 263 //=======================================================================
 264 const TopoDS_Face& IntTools_FaceFace::Face2() const
 265 {
 266   return myFace2;
 267 }
 268 //=======================================================================
 269 //function : TangentFaces
 270 //purpose  :
 271 //=======================================================================
 272 Standard_Boolean IntTools_FaceFace::TangentFaces() const
 273 {
 274   return myTangentFaces;
 275 }
 276 //=======================================================================
 277 //function : Points
 278 //purpose  :
 279 //=======================================================================
 280 const IntTools_SequenceOfPntOn2Faces& IntTools_FaceFace::Points() const
 281 {
 282   return myPnts;
 283 }
 284 //=======================================================================
 285 //function : IsDone
 286 //purpose  :
 287 //=======================================================================
 288 Standard_Boolean IntTools_FaceFace::IsDone() const
 289 {
 290   return myIsDone;
 291 }
 292 //=======================================================================
 293 //function : TolReached3d
 294 //purpose  :
 295 //=======================================================================
 296 Standard_Real IntTools_FaceFace::TolReached3d() const
 297 {
 298   return myTolReached3d;
 299 }
 300 //=======================================================================
 301 //function : Lines
 302 //purpose  : return lines of intersection
 303 //=======================================================================
 304 const IntTools_SequenceOfCurves& IntTools_FaceFace::Lines() const
 305 {
 306   StdFail_NotDone_Raise_if
 307     (!myIsDone,
 308      "IntTools_FaceFace::Lines() => myIntersector NOT DONE");
 309   return mySeqOfCurve;
 310 }
 311 //=======================================================================
 312 //function : TolReached2d
 313 //purpose  :
 314 //=======================================================================
 315 Standard_Real IntTools_FaceFace::TolReached2d() const
 316 {
 317   return myTolReached2d;
 318 }
 319 // =======================================================================
 320 // function: SetParameters
 321 //
 322 // =======================================================================
 323 void IntTools_FaceFace::SetParameters(const Standard_Boolean ToApproxC3d,
 324                                       const Standard_Boolean ToApproxC2dOnS1,
 325                                       const Standard_Boolean ToApproxC2dOnS2,
 326                                       const Standard_Real ApproximationTolerance)
 327 {
 328   myApprox = ToApproxC3d;
 329   myApprox1 = ToApproxC2dOnS1;
 330   myApprox2 = ToApproxC2dOnS2;
 331   myTolApprox = ApproximationTolerance;
 332 }
 333 //=======================================================================
 334 //function : SetList
 335 //purpose  :
 336 //=======================================================================
 337 void IntTools_FaceFace::SetList(IntSurf_ListOfPntOn2S& aListOfPnts)
 338 {
 339   myListOfPnts = aListOfPnts;
 340 }
 341
 342
 343 static Standard_Boolean isTreatAnalityc(const TopoDS_Face& theF1,
 344                                         const TopoDS_Face& theF2)
 345 {
 346   const Standard_Real Tolang = 1.e-8;
 347   const Standard_Real aTolF1=BRep_Tool::Tolerance(theF1);
 348   const Standard_Real aTolF2=BRep_Tool::Tolerance(theF2);
 349   const Standard_Real aTolSum = aTolF1 + aTolF2;
 350   Standard_Real aHigh = 0.0;
 351
 352   const BRepAdaptor_Surface aBAS1(theF1), aBAS2(theF2);
 353   const GeomAbs_SurfaceType aType1=aBAS1.GetType();
 354   const GeomAbs_SurfaceType aType2=aBAS2.GetType();
 355
 356   gp_Pln aS1;
 357   gp_Cylinder aS2;
 358   if(aType1 == GeomAbs_Plane)
 359   {
 360     aS1=aBAS1.Plane();
 361   }
 362   else if(aType2 == GeomAbs_Plane)
 363   {
 364     aS1=aBAS2.Plane();
 365   }
 366   else
 367   {
 368     return Standard_True;
 369   }
 370
 371   if(aType1 == GeomAbs_Cylinder)
 372   {
 373     aS2=aBAS1.Cylinder();
 374     const Standard_Real VMin = aBAS1.FirstVParameter();
 375     const Standard_Real VMax = aBAS1.LastVParameter();
 376
 377     if( Precision::IsNegativeInfinite(VMin) ||
 378         Precision::IsPositiveInfinite(VMax))
 379           return Standard_True;
 380     else
 381       aHigh = VMax - VMin;
 382   }
 383   else if(aType2 == GeomAbs_Cylinder)
 384   {
 385     aS2=aBAS2.Cylinder();
 386
 387     const Standard_Real VMin = aBAS2.FirstVParameter();
 388     const Standard_Real VMax = aBAS2.LastVParameter();
 389
 390     if( Precision::IsNegativeInfinite(VMin) ||
 391         Precision::IsPositiveInfinite(VMax))
 392           return Standard_True;
 393     else
 394       aHigh = VMax - VMin;
 395   }
 396   else
 397   {
 398     return Standard_True;
 399   }
 400
 401   IntAna_QuadQuadGeo inter;
 402   inter.Perform(aS1,aS2,Tolang,aTolSum, aHigh);
 403   if(inter.TypeInter() == IntAna_Ellipse)
 404   {
 405     const gp_Elips anEl = inter.Ellipse(1);
 406     const Standard_Real aMajorR = anEl.MajorRadius();
 407     const Standard_Real aMinorR = anEl.MinorRadius();
 408
 409     return (aMajorR < 100000.0 * aMinorR);
 410   }
 411   else
 412   {
 413     return inter.IsDone();
 414   }
 415 }
 416 //=======================================================================
 417 //function : Perform
 418 //purpose  : intersect surfaces of the faces
 419 //=======================================================================
 420 void IntTools_FaceFace::Perform(const TopoDS_Face& aF1,
 421                                 const TopoDS_Face& aF2)
 422 {
 423   Standard_Boolean RestrictLine = Standard_False, hasCone = Standard_False;
 424
 425   if (myContext.IsNull()) {
 426     myContext=new IntTools_Context;
 427   }
 428
 429   mySeqOfCurve.Clear();
 430   myTolReached2d=0.;
 431   myTolReached3d=0.;
 432   myIsDone = Standard_False;
 433   myNbrestr=0;//?
 434
 435   myFace1=aF1;
 436   myFace2=aF2;
 437
 438   const BRepAdaptor_Surface aBAS1(myFace1, Standard_False);
 439   const BRepAdaptor_Surface aBAS2(myFace2, Standard_False);
 440   GeomAbs_SurfaceType aType1=aBAS1.GetType();
 441   GeomAbs_SurfaceType aType2=aBAS2.GetType();
 442
 443   const Standard_Boolean bReverse=SortTypes(aType1, aType2);
 444   if (bReverse)
 445   {
 446     myFace1=aF2;
 447     myFace2=aF1;
 448     aType1=aBAS2.GetType();
 449     aType2=aBAS1.GetType();
 450
 451     if (myListOfPnts.Extent())
 452     {
 453       Standard_Real aU1,aV1,aU2,aV2;
 454       IntSurf_ListIteratorOfListOfPntOn2S aItP2S;
 455       //
 456       aItP2S.Initialize(myListOfPnts);
 457       for (; aItP2S.More(); aItP2S.Next())
 458       {
 459         IntSurf_PntOn2S& aP2S=aItP2S.Value();
 460         aP2S.Parameters(aU1,aV1,aU2,aV2);
 461         aP2S.SetValue(aU2,aV2,aU1,aV1);
 462       }
 463     }
 464     //
 465     Standard_Boolean anAproxTmp = myApprox1;
 466     myApprox1 = myApprox2;
 467     myApprox2 = anAproxTmp;
 468   }
 469
 470
 471   const Handle(Geom_Surface) S1=BRep_Tool::Surface(myFace1);
 472   const Handle(Geom_Surface) S2=BRep_Tool::Surface(myFace2);
 473
 474   const Standard_Real aTolF1=BRep_Tool::Tolerance(myFace1);
 475   const Standard_Real aTolF2=BRep_Tool::Tolerance(myFace2);
 476
 477   Standard_Real TolArc = aTolF1 + aTolF2;
 478   Standard_Real TolTang = TolArc;
 479
 480   const Standard_Boolean isFace1Quad = (aType1 == GeomAbs_Cylinder ||
 481                                         aType1 == GeomAbs_Cone ||
 482                                         aType1 == GeomAbs_Torus);
 483
 484   const Standard_Boolean isFace2Quad = (aType2 == GeomAbs_Cylinder ||
 485                                         aType2 == GeomAbs_Cone ||
 486                                         aType2 == GeomAbs_Torus);
 487
 488   if(aType1==GeomAbs_Plane && aType2==GeomAbs_Plane)  {
 489     Standard_Real umin, umax, vmin, vmax;
 490     //
 491     BRepTools::UVBounds(myFace1, umin, umax, vmin, vmax);
 492     CorrectPlaneBoundaries(umin, umax, vmin, vmax);
 493     myHS1->ChangeSurface().Load(S1, umin, umax, vmin, vmax);
 494     //
 495     BRepTools::UVBounds(myFace2, umin, umax, vmin, vmax);
 496     CorrectPlaneBoundaries(umin, umax, vmin, vmax);
 497     myHS2->ChangeSurface().Load(S2, umin, umax, vmin, vmax);
 498     //
 499     Standard_Real TolAng = 1.e-8;
 500     //
 501     PerformPlanes(myHS1, myHS2, TolAng, TolTang, myApprox1, myApprox2,
 502                   mySeqOfCurve, myTangentFaces);
 503     //
 504     myIsDone = Standard_True;
 505
 506     if(!myTangentFaces) {
 507       const Standard_Integer NbLinPP = mySeqOfCurve.Length();
 508       if(NbLinPP) {
 509         Standard_Real aTolFMax;
 510         myTolReached3d = 1.e-7;
 511         aTolFMax=Max(aTolF1, aTolF2);
 512         if (aTolFMax>myTolReached3d) {
 513           myTolReached3d=aTolFMax;
 514         }
 515         //
 516         myTolReached2d = myTolReached3d;
 517
 518         if (bReverse) {
 519           Handle(Geom2d_Curve) aC2D1, aC2D2;
 520           const Standard_Integer aNbLin = mySeqOfCurve.Length();
 521           for (Standard_Integer i = 1; i <= aNbLin; ++i) {
 522             IntTools_Curve& aIC=mySeqOfCurve(i);
 523             aC2D1=aIC.FirstCurve2d();
 524             aC2D2=aIC.SecondCurve2d();
 525             aIC.SetFirstCurve2d(aC2D2);
 526             aIC.SetSecondCurve2d(aC2D1);
 527           }
 528         }
 529       }
 530     }
 531     return;
 532   }//if(aType1==GeomAbs_Plane && aType2==GeomAbs_Plane){
 533
 534   if ((aType1==GeomAbs_Plane) && isFace2Quad)
 535   {
 536     Standard_Real umin, umax, vmin, vmax;
 537     // F1
 538     BRepTools::UVBounds(myFace1, umin, umax, vmin, vmax);
 539     CorrectPlaneBoundaries(umin, umax, vmin, vmax);
 540     myHS1->ChangeSurface().Load(S1, umin, umax, vmin, vmax);
 541     // F2
 542     BRepTools::UVBounds(myFace2, umin, umax, vmin, vmax);
 543     CorrectSurfaceBoundaries(myFace2, (aTolF1 + aTolF2) * 2., umin, umax, vmin, vmax);
 544     myHS2->ChangeSurface().Load(S2, umin, umax, vmin, vmax);
 545     //
 546     if( aType2==GeomAbs_Cone ) {
 547       TolArc = 0.0001;
 548       hasCone = Standard_True;
 549     }
 550   }
 551   else if ((aType2==GeomAbs_Plane) && isFace1Quad)
 552   {
 553     Standard_Real umin, umax, vmin, vmax;
 554     //F1
 555     BRepTools::UVBounds(myFace1, umin, umax, vmin, vmax);
 556     CorrectSurfaceBoundaries(myFace1, (aTolF1 + aTolF2) * 2., umin, umax, vmin, vmax);
 557     myHS1->ChangeSurface().Load(S1, umin, umax, vmin, vmax);
 558     // F2
 559     BRepTools::UVBounds(myFace2, umin, umax, vmin, vmax);
 560     CorrectPlaneBoundaries(umin, umax, vmin, vmax);
 561     myHS2->ChangeSurface().Load(S2, umin, umax, vmin, vmax);
 562     //
 563     if( aType1==GeomAbs_Cone ) {
 564       TolArc = 0.0001;
 565       hasCone = Standard_True;
 566     }
 567   }
 568   else
 569   {
 570     Standard_Real umin, umax, vmin, vmax;
 571     BRepTools::UVBounds(myFace1, umin, umax, vmin, vmax);
 572     CorrectSurfaceBoundaries(myFace1, (aTolF1 + aTolF2) * 2., umin, umax, vmin, vmax);
 573     myHS1->ChangeSurface().Load(S1, umin, umax, vmin, vmax);
 574     BRepTools::UVBounds(myFace2, umin, umax, vmin, vmax);
 575     CorrectSurfaceBoundaries(myFace2, (aTolF1 + aTolF2) * 2., umin, umax, vmin, vmax);
 576     myHS2->ChangeSurface().Load(S2, umin, umax, vmin, vmax);
 577   }
 578
 579   const Handle(IntTools_TopolTool) dom1 = new IntTools_TopolTool(myHS1);
 580   const Handle(IntTools_TopolTool) dom2 = new IntTools_TopolTool(myHS2);
 581
 582   myLConstruct.Load(dom1, dom2, myHS1, myHS2);
 583
 584
 585   Tolerances(myHS1, myHS2, TolTang);
 586
 587   {
 588     const Standard_Real UVMaxStep = 0.001;
 589     const Standard_Real Deflection = (hasCone) ? 0.085 : 0.1;
 590     myIntersector.SetTolerances(TolArc, TolTang, UVMaxStep, Deflection);
 591   }
 592
 593   if((myHS1->IsUClosed() && !myHS1->IsUPeriodic()) ||
 594      (myHS1->IsVClosed() && !myHS1->IsVPeriodic()) ||
 595      (myHS2->IsUClosed() && !myHS2->IsUPeriodic()) ||
 596      (myHS2->IsVClosed() && !myHS2->IsVPeriodic()))
 597   {
 598     RestrictLine = Standard_True;
 599   }
 600   //
 601   if((aType1 != GeomAbs_BSplineSurface) &&
 602       (aType1 != GeomAbs_BezierSurface)  &&
 603      (aType1 != GeomAbs_OtherSurface)  &&
 604      (aType2 != GeomAbs_BSplineSurface) &&
 605       (aType2 != GeomAbs_BezierSurface)  &&
 606      (aType2 != GeomAbs_OtherSurface))
 607   {
 608     RestrictLine = Standard_True;
 609
 610     if ((aType1 == GeomAbs_Torus) ||
 611         (aType2 == GeomAbs_Torus))
 612     {
 613       myListOfPnts.Clear();
 614     }
 615   }
 616
 617   //
 618   if(!RestrictLine)
 619   {
 620     TopExp_Explorer aExp;
 621     for(Standard_Integer i = 0; (!RestrictLine) && (i < 2); i++)
 622     {
 623       const TopoDS_Face& aF=(!i) ? myFace1 : myFace2;
 624       aExp.Init(aF, TopAbs_EDGE);
 625       for(; aExp.More(); aExp.Next())
 626       {
 627         const TopoDS_Edge& aE=TopoDS::Edge(aExp.Current());
 628
 629         if(BRep_Tool::Degenerated(aE))
 630         {
 631           RestrictLine = Standard_True;
 632           break;
 633         }
 634       }
 635     }
 636   }
 637
 638   const Standard_Boolean isGeomInt = isTreatAnalityc(aF1, aF2);
 639   myIntersector.Perform(myHS1, dom1, myHS2, dom2, TolArc, TolTang,
 640                                   myListOfPnts, RestrictLine, isGeomInt);
 641
 642   myIsDone = myIntersector.IsDone();
 643
 644   if (myIsDone)
 645   {
 646     myTangentFaces=myIntersector.TangentFaces();
 647     if (myTangentFaces) {
 648       return;
 649     }
 650     //
 651     if(RestrictLine) {
 652       myListOfPnts.Clear(); // to use LineConstructor
 653     }
 654     //
 655     const Standard_Integer aNbLinIntersector = myIntersector.NbLines();
 656     for (Standard_Integer i=1; i <= aNbLinIntersector; ++i) {
 657       MakeCurve(i, dom1, dom2, TolArc);
 658     }
 659     //
 660     ComputeTolReached3d();
 661     //
 662     if (bReverse) {
 663       Handle(Geom2d_Curve) aC2D1, aC2D2;
 664       //
 665       const Standard_Integer aNbLinSeqOfCurve =mySeqOfCurve.Length();
 666       for (Standard_Integer i=1; i<=aNbLinSeqOfCurve; ++i)
 667       {
 668         IntTools_Curve& aIC=mySeqOfCurve(i);
 669         aC2D1=aIC.FirstCurve2d();
 670         aC2D2=aIC.SecondCurve2d();
 671         aIC.SetFirstCurve2d(aC2D2);
 672         aIC.SetSecondCurve2d(aC2D1);
 673       }
 674     }
 675
 676     // Points
 677     Standard_Boolean bValid2D1, bValid2D2;
 678     Standard_Real U1,V1,U2,V2;
 679     IntTools_PntOnFace aPntOnF1, aPntOnF2;
 680     IntTools_PntOn2Faces aPntOn2Faces;
 681     //
 682     const Standard_Integer aNbPnts = myIntersector.NbPnts();
 683     for (Standard_Integer i=1; i <= aNbPnts; ++i)
 684     {
 685       const IntSurf_PntOn2S& aISPnt=myIntersector.Point(i).PntOn2S();
 686       const gp_Pnt& aPnt=aISPnt.Value();
 687       aISPnt.Parameters(U1,V1,U2,V2);
 688       //
 689       // check the validity of the intersection point for the faces
 690       bValid2D1 = myContext->IsPointInOnFace(myFace1, gp_Pnt2d(U1, V1));
 691       if (!bValid2D1) {
 692         continue;
 693       }
 694       //
 695       bValid2D2 = myContext->IsPointInOnFace(myFace2, gp_Pnt2d(U2, V2));
 696       if (!bValid2D2) {
 697         continue;
 698       }
 699       //
 700       // add the intersection point
 701       aPntOnF1.Init(myFace1, aPnt, U1, V1);
 702       aPntOnF2.Init(myFace2, aPnt, U2, V2);
 703       //
 704       if (!bReverse)
 705       {
 706         aPntOn2Faces.SetP1(aPntOnF1);
 707         aPntOn2Faces.SetP2(aPntOnF2);
 708       }
 709       else
 710       {
 711         aPntOn2Faces.SetP2(aPntOnF1);
 712         aPntOn2Faces.SetP1(aPntOnF2);
 713       }
 714
 715       myPnts.Append(aPntOn2Faces);
 716     }
 717   }
 718 }
 719
 720 //=======================================================================
 721 //function : ComputeTolerance
 722 //purpose  :
 723 //=======================================================================
 724 Standard_Real IntTools_FaceFace::ComputeTolerance()
 725 {
 726   Standard_Integer i, j, aNbLin;
 727   Standard_Real aFirst, aLast, aD, aDMax, aT;
 728   Handle(Geom_Surface) aS1, aS2;
 729   //
 730   aDMax = 0;
 731   aNbLin = mySeqOfCurve.Length();
 732   //
 733   aS1 = myHS1->ChangeSurface().Surface();
 734   aS2 = myHS2->ChangeSurface().Surface();
 735   //
 736   for (i = 1; i <= aNbLin; ++i)
 737   {
 738     const IntTools_Curve& aIC = mySeqOfCurve(i);
 739     const Handle(Geom_Curve)& aC3D = aIC.Curve();
 740     if (aC3D.IsNull())
 741     {
 742       continue;
 743     }
 744     //
 745     aFirst = aC3D->FirstParameter();
 746     aLast  = aC3D->LastParameter();
 747     //
 748     const Handle(Geom2d_Curve)& aC2D1 = aIC.FirstCurve2d();
 749     const Handle(Geom2d_Curve)& aC2D2 = aIC.SecondCurve2d();
 750     //
 751     for (j = 0; j < 2; ++j)
 752     {
 753       const Handle(Geom2d_Curve)& aC2D = !j ? aC2D1 : aC2D2;
 754       const Handle(Geom_Surface)& aS = !j ? aS1 : aS2;
 755       //
 756       if (!aC2D.IsNull())
 757       {
 758         if (IntTools_Tools::ComputeTolerance
 759             (aC3D, aC2D, aS, aFirst, aLast, aD, aT))
 760         {
 761           if (aD > aDMax)
 762           {
 763             aDMax = aD;
 764           }
 765         }
 766       }
 767       else
 768       {
 769         const TopoDS_Face& aF = !j ? myFace1 : myFace2;
 770         aD = FindMaxDistance(aC3D, aFirst, aLast, aF, myContext);
 771         if (aD > aDMax)
 772         {
 773           aDMax = aD;
 774         }
 775       }
 776     }
 777   }
 778   //
 779   return aDMax;
 780 }
 781
 782 //=======================================================================
 783 //function :ComputeTolReached3d
 784 //purpose  :
 785 //=======================================================================
 786   void IntTools_FaceFace::ComputeTolReached3d()
 787 {
 788   Standard_Integer aNbLin;
 789   GeomAbs_SurfaceType aType1, aType2;
 790   //
 791   aNbLin=myIntersector.NbLines();
 792   if (!aNbLin) {
 793     return;
 794   }
 795   //
 796   aType1=myHS1->Surface().GetType();
 797   aType2=myHS2->Surface().GetType();
 798   //
 799   if (aType1==GeomAbs_Cylinder && aType2==GeomAbs_Cylinder)
 800   {
 801     if (aNbLin==2)
 802     {
 803       Handle(IntPatch_Line) aIL1, aIL2;
 804       IntPatch_IType aTL1, aTL2;
 805       //
 806       aIL1=myIntersector.Line(1);
 807       aIL2=myIntersector.Line(2);
 808       aTL1=aIL1->ArcType();
 809       aTL2=aIL2->ArcType();
 810       if (aTL1==IntPatch_Lin && aTL2==IntPatch_Lin) {
 811         Standard_Real aD, aDTresh, dTol;
 812         gp_Lin aL1, aL2;
 813         //
 814         dTol=1.e-8;
 815         aDTresh=1.5e-6;
 816         //
 817         aL1=Handle(IntPatch_GLine)::DownCast(aIL1)->Line();
 818         aL2=Handle(IntPatch_GLine)::DownCast(aIL2)->Line();
 819         aD=aL1.Distance(aL2);
 820         aD=0.5*aD;
 821         if (aD<aDTresh)
 822         {//In order to avoid creation too thin face
 823           myTolReached3d=aD+dTol;
 824         }
 825       }
 826     }
 827   }// if (aType1==GeomAbs_Cylinder && aType2==GeomAbs_Cylinder) {
 828   //
 829
 830   Standard_Real aDMax = ComputeTolerance();
 831   if (aDMax > myTolReached3d)
 832   {
 833       myTolReached3d = aDMax;
 834     }
 835   }
 836
 837 //=======================================================================
 838 //function : MakeCurve
 839 //purpose  :
 840 //=======================================================================
 841 void IntTools_FaceFace::MakeCurve(const Standard_Integer Index,
 842                                   const Handle(Adaptor3d_TopolTool)& dom1,
 843                                   const Handle(Adaptor3d_TopolTool)& dom2,
 844                                   const Standard_Real theToler)
 845 {
 846   Standard_Boolean bDone, rejectSurface, reApprox, bAvoidLineConstructor;
 847   Standard_Boolean ok, bPCurvesOk;
 848   Standard_Integer i, j, aNbParts;
 849   Standard_Real fprm, lprm;
 850   Standard_Real Tolpc;
 851   Handle(IntPatch_Line) L;
 852   IntPatch_IType typl;
 853   Handle(Geom_Curve) newc;
 854   //
 855   const Standard_Real TOLCHECK   =0.0000001;
 856   const Standard_Real TOLANGCHECK=0.1;
 857   //
 858   rejectSurface = Standard_False;
 859   reApprox = Standard_False;
 860   //
 861   bPCurvesOk = Standard_True;
 862
 863  reapprox:;
 864
 865   Tolpc = myTolApprox;
 866   bAvoidLineConstructor = Standard_False;
 867   L = myIntersector.Line(Index);
 868   typl = L->ArcType();
 869   //
 870   if(typl==IntPatch_Walking) {
 871     Handle(IntPatch_WLine) aWLine (Handle(IntPatch_WLine)::DownCast(L));
 872     if(aWLine.IsNull()) {
 873       return;
 874     }
 875     L = aWLine;
 876
 877     //
 878     //if(!myListOfPnts.IsEmpty()) {
 879     //  bAvoidLineConstructor = Standard_True;
 880     //}
 881
 882     Standard_Integer nbp = aWLine->NbPnts();
 883     const IntSurf_PntOn2S& p1 = aWLine->Point(1);
 884     const IntSurf_PntOn2S& p2 = aWLine->Point(nbp);
 885
 886     const gp_Pnt& P1 = p1.Value();
 887     const gp_Pnt& P2 = p2.Value();
 888
 889     if(P1.SquareDistance(P2) < 1.e-14) {
 890       bAvoidLineConstructor = Standard_False;
 891     }
 892   }
 893
 894   typl=L->ArcType();
 895
 896   //
 897   // Line Constructor
 898   if(!bAvoidLineConstructor) {
 899     myLConstruct.Perform(L);
 900     //
 901     bDone=myLConstruct.IsDone();
 902     if(!bDone)
 903     {
 904       return;
 905     }
 906
 907     if(typl != IntPatch_Restriction)
 908     {
 909       aNbParts=myLConstruct.NbParts();
 910       if (aNbParts <= 0)
 911       {
 912         return;
 913       }
 914     }
 915   }
 916   // Do the Curve
 917
 918
 919   switch (typl) {
 920   //########################################
 921   // Line, Parabola, Hyperbola
 922   //########################################
 923   case IntPatch_Lin:
 924   case IntPatch_Parabola:
 925   case IntPatch_Hyperbola: {
 926     if (typl == IntPatch_Lin) {
 927       newc =
 928         new Geom_Line (Handle(IntPatch_GLine)::DownCast(L)->Line());
 929     }
 930
 931     else if (typl == IntPatch_Parabola) {
 932       newc =
 933         new Geom_Parabola(Handle(IntPatch_GLine)::DownCast(L)->Parabola());
 934     }
 935
 936     else if (typl == IntPatch_Hyperbola) {
 937       newc =
 938         new Geom_Hyperbola (Handle(IntPatch_GLine)::DownCast(L)->Hyperbola());
 939     }
 940     //
 941     // myTolReached3d
 942     if (typl == IntPatch_Lin) {
 943       TolR3d (myFace1, myFace2, myTolReached3d);
 944     }
 945     //
 946     aNbParts=myLConstruct.NbParts();
 947     for (i=1; i<=aNbParts; i++) {
 948       Standard_Boolean bFNIt, bLPIt;
 949       //
 950       myLConstruct.Part(i, fprm, lprm);
 951         //
 952       bFNIt=Precision::IsNegativeInfinite(fprm);
 953       bLPIt=Precision::IsPositiveInfinite(lprm);
 954       //
 955       if (!bFNIt && !bLPIt) {
 956         //
 957         IntTools_Curve aCurve;
 958         //
 959         Handle(Geom_TrimmedCurve) aCT3D=new Geom_TrimmedCurve(newc, fprm, lprm);
 960         aCurve.SetCurve(aCT3D);
 961         if (typl == IntPatch_Parabola) {
 962           Standard_Real aTolF1, aTolF2, aTolBase;
 963
 964           aTolF1 = BRep_Tool::Tolerance(myFace1);
 965           aTolF2 = BRep_Tool::Tolerance(myFace2);
 966           aTolBase=aTolF1+aTolF2;
 967           myTolReached3d=IntTools_Tools::CurveTolerance(aCT3D, aTolBase);
 968         }
 969         //
 970         aCurve.SetCurve(new Geom_TrimmedCurve(newc, fprm, lprm));
 971         if(myApprox1) {
 972           Handle (Geom2d_Curve) C2d;
 973           GeomInt_IntSS::BuildPCurves(fprm, lprm, Tolpc,
 974                 myHS1->ChangeSurface().Surface(), newc, C2d);
 975           if(Tolpc>myTolReached2d || myTolReached2d==0.) {
 976             myTolReached2d=Tolpc;
 977           }
 978             //
 979             aCurve.SetFirstCurve2d(new Geom2d_TrimmedCurve(C2d,fprm,lprm));
 980           }
 981           else {
 982             Handle(Geom2d_BSplineCurve) H1;
 983             //
 984             aCurve.SetFirstCurve2d(H1);
 985           }
 986         //
 987         if(myApprox2) {
 988           Handle (Geom2d_Curve) C2d;
 989           GeomInt_IntSS::BuildPCurves(fprm, lprm, Tolpc,
 990                     myHS2->ChangeSurface().Surface(), newc, C2d);
 991           if(Tolpc>myTolReached2d || myTolReached2d==0.) {
 992             myTolReached2d=Tolpc;
 993           }
 994           //
 995           aCurve.SetSecondCurve2d(new Geom2d_TrimmedCurve(C2d,fprm,lprm));
 996           }
 997         else {
 998           Handle(Geom2d_BSplineCurve) H1;
 999           //
1000           aCurve.SetSecondCurve2d(H1);
1001         }
1002         mySeqOfCurve.Append(aCurve);
1003       } //if (!bFNIt && !bLPIt) {
1004       else {
1005         //  on regarde si on garde
1006         //
1007         Standard_Real aTestPrm, dT=100.;
1008         //
1009         aTestPrm=0.;
1010         if (bFNIt && !bLPIt) {
1011           aTestPrm=lprm-dT;
1012         }
1013         else if (!bFNIt && bLPIt) {
1014           aTestPrm=fprm+dT;
1015         }
1016         else {
1017           // i.e, if (bFNIt && bLPIt)
1018           aTestPrm=IntTools_Tools::IntermediatePoint(-dT, dT);
1019         }
1020         //
1021         gp_Pnt ptref(newc->Value(aTestPrm));
1022         //
1023         GeomAbs_SurfaceType typS1 = myHS1->GetType();
1024         GeomAbs_SurfaceType typS2 = myHS2->GetType();
1025         if( typS1 == GeomAbs_SurfaceOfExtrusion ||
1026             typS1 == GeomAbs_OffsetSurface ||
1027             typS1 == GeomAbs_SurfaceOfRevolution ||
1028             typS2 == GeomAbs_SurfaceOfExtrusion ||
1029             typS2 == GeomAbs_OffsetSurface ||
1030             typS2 == GeomAbs_SurfaceOfRevolution) {
1031           Handle(Geom2d_BSplineCurve) H1;
1032           mySeqOfCurve.Append(IntTools_Curve(newc, H1, H1));
1033           continue;
1034         }
1035
1036         Standard_Real u1, v1, u2, v2, Tol;
1037
1038         Tol = Precision::Confusion();
1039         Parameters(myHS1, myHS2, ptref,  u1, v1, u2, v2);
1040         ok = (dom1->Classify(gp_Pnt2d(u1, v1), Tol) != TopAbs_OUT);
1041         if(ok) {
1042           ok = (dom2->Classify(gp_Pnt2d(u2,v2),Tol) != TopAbs_OUT);
1043         }
1044         if (ok) {
1045           Handle(Geom2d_BSplineCurve) H1;
1046           mySeqOfCurve.Append(IntTools_Curve(newc, H1, H1));
1047         }
1048       }
1049     }// for (i=1; i<=aNbParts; i++) {
1050   }// case IntPatch_Lin:  case IntPatch_Parabola:  case IntPatch_Hyperbola:
1051     break;
1052
1053   //########################################
1054   // Circle and Ellipse
1055   //########################################
1056   case IntPatch_Circle:
1057   case IntPatch_Ellipse: {
1058
1059     if (typl == IntPatch_Circle) {
1060       newc = new Geom_Circle
1061         (Handle(IntPatch_GLine)::DownCast(L)->Circle());
1062     }
1063     else { //IntPatch_Ellipse
1064       newc = new Geom_Ellipse
1065         (Handle(IntPatch_GLine)::DownCast(L)->Ellipse());
1066     }
1067     //
1068     // myTolReached3d
1069     TolR3d (myFace1, myFace2, myTolReached3d);
1070     //
1071     aNbParts=myLConstruct.NbParts();
1072     //
1073     Standard_Real aPeriod, aNul;
1074     TColStd_SequenceOfReal aSeqFprm,  aSeqLprm;
1075
1076     aNul=0.;
1077     aPeriod=M_PI+M_PI;
1078
1079     for (i=1; i<=aNbParts; i++) {
1080       myLConstruct.Part(i, fprm, lprm);
1081
1082       if (fprm < aNul && lprm > aNul) {
1083         // interval that goes through 0. is divided on two intervals;
1084         while (fprm<aNul || fprm>aPeriod) fprm=fprm+aPeriod;
1085         while (lprm<aNul || lprm>aPeriod) lprm=lprm+aPeriod;
1086         //
1087         if((aPeriod - fprm) > Tolpc) {
1088           aSeqFprm.Append(fprm);
1089           aSeqLprm.Append(aPeriod);
1090         }
1091         else {
1092           gp_Pnt P1 = newc->Value(fprm);
1093           gp_Pnt P2 = newc->Value(aPeriod);
1094           Standard_Real aTolDist = BRep_Tool::Tolerance(myFace1) + BRep_Tool::Tolerance(myFace2);
1095           aTolDist = (myTolReached3d > aTolDist) ? myTolReached3d : aTolDist;
1096
1097           if(P1.Distance(P2) > aTolDist) {
1098             Standard_Real anewpar = fprm;
1099
1100             if(ParameterOutOfBoundary(fprm, newc, myFace1, myFace2,
1101                                       lprm, Standard_False, anewpar, myContext)) {
1102               fprm = anewpar;
1103             }
1104             aSeqFprm.Append(fprm);
1105             aSeqLprm.Append(aPeriod);
1106           }
1107         }
1108
1109         //
1110         if((lprm - aNul) > Tolpc) {
1111           aSeqFprm.Append(aNul);
1112           aSeqLprm.Append(lprm);
1113         }
1114         else {
1115           gp_Pnt P1 = newc->Value(aNul);
1116           gp_Pnt P2 = newc->Value(lprm);
1117           Standard_Real aTolDist = BRep_Tool::Tolerance(myFace1) + BRep_Tool::Tolerance(myFace2);
1118           aTolDist = (myTolReached3d > aTolDist) ? myTolReached3d : aTolDist;
1119
1120           if(P1.Distance(P2) > aTolDist) {
1121             Standard_Real anewpar = lprm;
1122
1123             if(ParameterOutOfBoundary(lprm, newc, myFace1, myFace2,
1124                                       fprm, Standard_True, anewpar, myContext)) {
1125               lprm = anewpar;
1126             }
1127             aSeqFprm.Append(aNul);
1128             aSeqLprm.Append(lprm);
1129           }
1130         }
1131       }
1132       else {
1133         // usual interval
1134         aSeqFprm.Append(fprm);
1135         aSeqLprm.Append(lprm);
1136       }
1137     }
1138
1139     //
1140     aNbParts=aSeqFprm.Length();
1141     for (i=1; i<=aNbParts; i++) {
1142       fprm=aSeqFprm(i);
1143       lprm=aSeqLprm(i);
1144       //
1145       Standard_Real aRealEpsilon=RealEpsilon();
1146       if (Abs(fprm) > aRealEpsilon || Abs(lprm-2.*M_PI) > aRealEpsilon) {
1147         //==============================================
1148         ////
1149         IntTools_Curve aCurve;
1150         Handle(Geom_TrimmedCurve) aTC3D=new Geom_TrimmedCurve(newc,fprm,lprm);
1151         aCurve.SetCurve(aTC3D);
1152         fprm=aTC3D->FirstParameter();
1153         lprm=aTC3D->LastParameter ();
1154         ////
1155         if (typl == IntPatch_Circle || typl == IntPatch_Ellipse) {////
1156           if(myApprox1) {
1157             Handle (Geom2d_Curve) C2d;
1158             GeomInt_IntSS::BuildPCurves(fprm, lprm, Tolpc,
1159                         myHS1->ChangeSurface().Surface(), newc, C2d);
1160             if(Tolpc>myTolReached2d || myTolReached2d==0) {
1161               myTolReached2d=Tolpc;
1162             }
1163             //
1164             aCurve.SetFirstCurve2d(C2d);
1165           }
1166           else { ////
1167             Handle(Geom2d_BSplineCurve) H1;
1168             aCurve.SetFirstCurve2d(H1);
1169           }
1170
1171
1172           if(myApprox2) {
1173             Handle (Geom2d_Curve) C2d;
1174             GeomInt_IntSS::BuildPCurves(fprm,lprm,Tolpc,
1175                     myHS2->ChangeSurface().Surface(),newc,C2d);
1176             if(Tolpc>myTolReached2d || myTolReached2d==0) {
1177               myTolReached2d=Tolpc;
1178             }
1179             //
1180             aCurve.SetSecondCurve2d(C2d);
1181           }
1182           else {
1183             Handle(Geom2d_BSplineCurve) H1;
1184             aCurve.SetSecondCurve2d(H1);
1185           }
1186         }
1187
1188         else {
1189           Handle(Geom2d_BSplineCurve) H1;
1190           aCurve.SetFirstCurve2d(H1);
1191           aCurve.SetSecondCurve2d(H1);
1192         }
1193         mySeqOfCurve.Append(aCurve);
1194           //==============================================
1195       } //if (Abs(fprm) > RealEpsilon() || Abs(lprm-2.*M_PI) > RealEpsilon())
1196
1197       else {
1198         //  on regarde si on garde
1199         //
1200         if (aNbParts==1) {
1201 //           if (Abs(fprm) < RealEpsilon() &&  Abs(lprm-2.*M_PI) < RealEpsilon()) {
1202           if (Abs(fprm) <= aRealEpsilon && Abs(lprm-2.*M_PI) <= aRealEpsilon) {
1203             IntTools_Curve aCurve;
1204             Handle(Geom_TrimmedCurve) aTC3D=new Geom_TrimmedCurve(newc,fprm,lprm);
1205             aCurve.SetCurve(aTC3D);
1206             fprm=aTC3D->FirstParameter();
1207             lprm=aTC3D->LastParameter ();
1208
1209             if(myApprox1) {
1210               Handle (Geom2d_Curve) C2d;
1211               GeomInt_IntSS::BuildPCurves(fprm,lprm,Tolpc,
1212                     myHS1->ChangeSurface().Surface(),newc,C2d);
1213               if(Tolpc>myTolReached2d || myTolReached2d==0) {
1214                 myTolReached2d=Tolpc;
1215               }
1216               //
1217               aCurve.SetFirstCurve2d(C2d);
1218             }
1219             else { ////
1220               Handle(Geom2d_BSplineCurve) H1;
1221               aCurve.SetFirstCurve2d(H1);
1222             }
1223
1224             if(myApprox2) {
1225               Handle (Geom2d_Curve) C2d;
1226               GeomInt_IntSS::BuildPCurves(fprm,lprm,Tolpc,
1227                     myHS2->ChangeSurface().Surface(),newc,C2d);
1228               if(Tolpc>myTolReached2d || myTolReached2d==0) {
1229                 myTolReached2d=Tolpc;
1230               }
1231               //
1232               aCurve.SetSecondCurve2d(C2d);
1233             }
1234             else {
1235               Handle(Geom2d_BSplineCurve) H1;
1236               aCurve.SetSecondCurve2d(H1);
1237             }
1238             mySeqOfCurve.Append(aCurve);
1239             break;
1240           }
1241         }
1242         //
1243         Standard_Real aTwoPIdiv17, u1, v1, u2, v2, Tol;
1244
1245         aTwoPIdiv17=2.*M_PI/17.;
1246
1247         for (j=0; j<=17; j++) {
1248           gp_Pnt ptref (newc->Value (j*aTwoPIdiv17));
1249           Tol = Precision::Confusion();
1250
1251           Parameters(myHS1, myHS2, ptref, u1, v1, u2, v2);
1252           ok = (dom1->Classify(gp_Pnt2d(u1,v1),Tol) != TopAbs_OUT);
1253           if(ok) {
1254             ok = (dom2->Classify(gp_Pnt2d(u2,v2),Tol) != TopAbs_OUT);
1255           }
1256           if (ok) {
1257             IntTools_Curve aCurve;
1258             aCurve.SetCurve(newc);
1259             //==============================================
1260             if (typl == IntPatch_Circle || typl == IntPatch_Ellipse) {
1261
1262               if(myApprox1) {
1263                 Handle (Geom2d_Curve) C2d;
1264                 GeomInt_IntSS::BuildPCurves(fprm, lprm, Tolpc,
1265                         myHS1->ChangeSurface().Surface(), newc, C2d);
1266                 if(Tolpc>myTolReached2d || myTolReached2d==0) {
1267                   myTolReached2d=Tolpc;
1268                 }
1269                 //
1270                 aCurve.SetFirstCurve2d(C2d);
1271               }
1272               else {
1273                 Handle(Geom2d_BSplineCurve) H1;
1274                 aCurve.SetFirstCurve2d(H1);
1275               }
1276
1277               if(myApprox2) {
1278                 Handle (Geom2d_Curve) C2d;
1279                 GeomInt_IntSS::BuildPCurves(fprm, lprm, Tolpc,
1280                         myHS2->ChangeSurface().Surface(), newc, C2d);
1281                 if(Tolpc>myTolReached2d || myTolReached2d==0) {
1282                   myTolReached2d=Tolpc;
1283                 }
1284                 //
1285                 aCurve.SetSecondCurve2d(C2d);
1286               }
1287
1288               else {
1289                 Handle(Geom2d_BSplineCurve) H1;
1290                 aCurve.SetSecondCurve2d(H1);
1291               }
1292             }//  end of if (typl == IntPatch_Circle || typl == IntPatch_Ellipse)
1293
1294             else {
1295               Handle(Geom2d_BSplineCurve) H1;
1296               //
1297               aCurve.SetFirstCurve2d(H1);
1298               aCurve.SetSecondCurve2d(H1);
1299             }
1300             //==============================================
1301             //
1302             mySeqOfCurve.Append(aCurve);
1303             break;
1304
1305             }//  end of if (ok) {
1306           }//  end of for (Standard_Integer j=0; j<=17; j++)
1307         }//  end of else { on regarde si on garde
1308       }// for (i=1; i<=myLConstruct.NbParts(); i++)
1309     }// IntPatch_Circle: IntPatch_Ellipse:
1310     break;
1311
1312   case IntPatch_Analytic: {
1313     IntSurf_Quadric quad1,quad2;
1314     GeomAbs_SurfaceType typs = myHS1->Surface().GetType();
1315
1316     switch (typs) {
1317       case GeomAbs_Plane:
1318         quad1.SetValue(myHS1->Surface().Plane());
1319         break;
1320       case GeomAbs_Cylinder:
1321         quad1.SetValue(myHS1->Surface().Cylinder());
1322         break;
1323       case GeomAbs_Cone:
1324         quad1.SetValue(myHS1->Surface().Cone());
1325         break;
1326       case GeomAbs_Sphere:
1327         quad1.SetValue(myHS1->Surface().Sphere());
1328         break;
1329     case GeomAbs_Torus:
1330       quad1.SetValue(myHS1->Surface().Torus());
1331       break;
1332       default:
1333         Standard_ConstructionError::Raise("GeomInt_IntSS::MakeCurve 1");
1334       }
1335
1336     typs = myHS2->Surface().GetType();
1337
1338     switch (typs) {
1339       case GeomAbs_Plane:
1340         quad2.SetValue(myHS2->Surface().Plane());
1341         break;
1342       case GeomAbs_Cylinder:
1343         quad2.SetValue(myHS2->Surface().Cylinder());
1344         break;
1345       case GeomAbs_Cone:
1346         quad2.SetValue(myHS2->Surface().Cone());
1347         break;
1348       case GeomAbs_Sphere:
1349         quad2.SetValue(myHS2->Surface().Sphere());
1350         break;
1351     case GeomAbs_Torus:
1352       quad2.SetValue(myHS2->Surface().Torus());
1353       break;
1354       default:
1355         Standard_ConstructionError::Raise("GeomInt_IntSS::MakeCurve 2");
1356       }
1357     //
1358     //=========
1359     IntPatch_ALineToWLine convert (quad1, quad2);
1360
1361     if (!myApprox) {
1362       aNbParts=myLConstruct.NbParts();
1363       for (i=1; i<=aNbParts; i++) {
1364         myLConstruct.Part(i, fprm, lprm);
1365         Handle(IntPatch_WLine) WL =
1366           convert.MakeWLine(Handle(IntPatch_ALine)::DownCast(L), fprm, lprm);
1367         //
1368         Handle(Geom2d_BSplineCurve) H1;
1369         Handle(Geom2d_BSplineCurve) H2;
1370
1371         if(myApprox1) {
1372           H1 = GeomInt_IntSS::MakeBSpline2d(WL, 1, WL->NbPnts(), Standard_True);
1373         }
1374
1375         if(myApprox2) {
1376           H2 = GeomInt_IntSS::MakeBSpline2d(WL, 1, WL->NbPnts(), Standard_False);
1377         }
1378         //
1379         mySeqOfCurve.Append(IntTools_Curve(GeomInt_IntSS::MakeBSpline(WL,1,WL->NbPnts()), H1, H2));
1380       }
1381     } // if (!myApprox)
1382
1383     else { // myApprox=TRUE
1384       GeomInt_WLApprox theapp3d;
1385       //
1386       Standard_Real tol2d = myTolApprox;
1387       //
1388       theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, 30, Standard_True);
1389
1390       aNbParts=myLConstruct.NbParts();
1391       for (i=1; i<=aNbParts; i++) {
1392         myLConstruct.Part(i, fprm, lprm);
1393         Handle(IntPatch_WLine) WL =
1394           convert.MakeWLine(Handle(IntPatch_ALine):: DownCast(L),fprm,lprm);
1395
1396         theapp3d.Perform(myHS1,myHS2,WL,Standard_True,myApprox1,myApprox2, 1, WL->NbPnts());
1397
1398         if (!theapp3d.IsDone()) {
1399           //
1400           Handle(Geom2d_BSplineCurve) H1;
1401           Handle(Geom2d_BSplineCurve) H2;
1402
1403           if(myApprox1) {
1404             H1 = GeomInt_IntSS::MakeBSpline2d(WL, 1, WL->NbPnts(), Standard_True);
1405           }
1406
1407           if(myApprox2) {
1408             H2 = GeomInt_IntSS::MakeBSpline2d(WL, 1, WL->NbPnts(), Standard_False);
1409           }
1410           //
1411           mySeqOfCurve.Append(IntTools_Curve(GeomInt_IntSS::MakeBSpline(WL,1,WL->NbPnts()), H1, H2));
1412         }
1413
1414         else {
1415           if(myApprox1 || myApprox2) {
1416             if( theapp3d.TolReached2d()>myTolReached2d || myTolReached2d==0) {
1417               myTolReached2d = theapp3d.TolReached2d();
1418             }
1419           }
1420
1421           if( theapp3d.TolReached3d()>myTolReached3d || myTolReached3d==0) {
1422             myTolReached3d = theapp3d.TolReached3d();
1423           }
1424
1425           Standard_Integer aNbMultiCurves, nbpoles;
1426           aNbMultiCurves=theapp3d.NbMultiCurves();
1427           for (j=1; j<=aNbMultiCurves; j++) {
1428             const AppParCurves_MultiBSpCurve& mbspc = theapp3d.Value(j);
1429             nbpoles = mbspc.NbPoles();
1430
1431             TColgp_Array1OfPnt tpoles(1, nbpoles);
1432             mbspc.Curve(1, tpoles);
1433             Handle(Geom_BSplineCurve) BS=new Geom_BSplineCurve(tpoles,
1434                                                                mbspc.Knots(),
1435                                                                mbspc.Multiplicities(),
1436                                                                mbspc.Degree());
1437
1438             GeomLib_CheckBSplineCurve Check(BS,TOLCHECK,TOLANGCHECK);
1439             Check.FixTangent(Standard_True,Standard_True);
1440             //
1441             IntTools_Curve aCurve;
1442             aCurve.SetCurve(BS);
1443
1444             if(myApprox1) {
1445               TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
1446               mbspc.Curve(2,tpoles2d);
1447               Handle(Geom2d_BSplineCurve) BS2=new Geom2d_BSplineCurve(tpoles2d,
1448                                                                       mbspc.Knots(),
1449                                                                       mbspc.Multiplicities(),
1450                                                                       mbspc.Degree());
1451
1452               GeomLib_Check2dBSplineCurve newCheck(BS2,TOLCHECK,TOLANGCHECK);
1453               newCheck.FixTangent(Standard_True,Standard_True);
1454               //
1455               aCurve.SetFirstCurve2d(BS2);
1456             }
1457             else {
1458               Handle(Geom2d_BSplineCurve) H1;
1459               aCurve.SetFirstCurve2d(H1);
1460             }
1461
1462             if(myApprox2) {
1463               TColgp_Array1OfPnt2d tpoles2d(1, nbpoles);
1464               Standard_Integer TwoOrThree;
1465               TwoOrThree=myApprox1 ? 3 : 2;
1466               mbspc.Curve(TwoOrThree, tpoles2d);
1467               Handle(Geom2d_BSplineCurve) BS2 =new Geom2d_BSplineCurve(tpoles2d,
1468                                                                        mbspc.Knots(),
1469                                                                        mbspc.Multiplicities(),
1470                                                                        mbspc.Degree());
1471
1472               GeomLib_Check2dBSplineCurve newCheck(BS2,TOLCHECK,TOLANGCHECK);
1473               newCheck.FixTangent(Standard_True,Standard_True);
1474               //
1475               aCurve.SetSecondCurve2d(BS2);
1476             }
1477             else {
1478               Handle(Geom2d_BSplineCurve) H2;
1479               aCurve.SetSecondCurve2d(H2);
1480             }
1481             //
1482             mySeqOfCurve.Append(aCurve);
1483
1484           }// for (j=1; j<=aNbMultiCurves; j++) {
1485         }// else from if (!theapp3d.IsDone())
1486       }// for (i=1; i<=aNbParts; i++) {
1487     }// else { // myApprox=TRUE
1488   }// case IntPatch_Analytic:
1489     break;
1490
1491   case IntPatch_Walking:{
1492     Handle(IntPatch_WLine) WL =
1493       Handle(IntPatch_WLine)::DownCast(L);
1494
1495 #ifdef OCCT_DEBUG
1496     //WL->Dump(0);
1497 #endif
1498
1499     //
1500     Standard_Integer ifprm, ilprm;
1501     //
1502     if (!myApprox) {
1503       aNbParts = 1;
1504       if(!bAvoidLineConstructor){
1505         aNbParts=myLConstruct.NbParts();
1506       }
1507       for (i=1; i<=aNbParts; ++i) {
1508         Handle(Geom2d_BSplineCurve) H1, H2;
1509         Handle(Geom_Curve) aBSp;
1510         //
1511         if(bAvoidLineConstructor) {
1512           ifprm = 1;
1513           ilprm = WL->NbPnts();
1514         }
1515         else {
1516           myLConstruct.Part(i, fprm, lprm);
1517           ifprm=(Standard_Integer)fprm;
1518           ilprm=(Standard_Integer)lprm;
1519         }
1520         //
1521         if(myApprox1) {
1522           H1 = GeomInt_IntSS::MakeBSpline2d(WL, ifprm, ilprm, Standard_True);
1523         }
1524         //
1525         if(myApprox2) {
1526           H2 = GeomInt_IntSS::MakeBSpline2d(WL, ifprm, ilprm, Standard_False);
1527         }
1528         //
1529         aBSp=GeomInt_IntSS::MakeBSpline(WL, ifprm, ilprm);
1530         IntTools_Curve aIC(aBSp, H1, H2);
1531         mySeqOfCurve.Append(aIC);
1532       }// for (i=1; i<=aNbParts; ++i) {
1533     }// if (!myApprox) {
1534     //
1535     else { // X
1536       Standard_Boolean bIsDecomposited;
1537       Standard_Integer nbiter, aNbSeqOfL;
1538       Standard_Real tol2d, aTolApproxImp;
1539       IntPatch_SequenceOfLine aSeqOfL;
1540       GeomInt_WLApprox theapp3d;
1541       Approx_ParametrizationType aParType = Approx_ChordLength;
1542       //
1543       Standard_Boolean anApprox1 = myApprox1;
1544       Standard_Boolean anApprox2 = myApprox2;
1545       //
1546       aTolApproxImp=1.e-5;
1547       tol2d = myTolApprox;
1548
1549       GeomAbs_SurfaceType typs1, typs2;
1550       typs1 = myHS1->Surface().GetType();
1551       typs2 = myHS2->Surface().GetType();
1552       Standard_Boolean anWithPC = Standard_True;
1553
1554       if(typs1 == GeomAbs_Cylinder && typs2 == GeomAbs_Sphere) {
1555         anWithPC =
1556           ApproxWithPCurves(myHS1->Surface().Cylinder(), myHS2->Surface().Sphere());
1557       }
1558       else if (typs1 == GeomAbs_Sphere && typs2 == GeomAbs_Cylinder) {
1559         anWithPC =
1560           ApproxWithPCurves(myHS2->Surface().Cylinder(), myHS1->Surface().Sphere());
1561       }
1562       //
1563       if(!anWithPC) {
1564         myTolApprox = aTolApproxImp;//1.e-5;
1565         anApprox1 = Standard_False;
1566         anApprox2 = Standard_False;
1567         //
1568         tol2d = myTolApprox;
1569       }
1570
1571       if(myHS1 == myHS2) {
1572         theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, 30, Standard_False, aParType);
1573         rejectSurface = Standard_True;
1574       }
1575       else {
1576         if(reApprox && !rejectSurface)
1577           theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, 30, Standard_False, aParType);
1578         else {
1579           Standard_Integer iDegMax, iDegMin, iNbIter;
1580           //
1581           ApproxParameters(myHS1, myHS2, iDegMin, iDegMax, iNbIter);
1582           theapp3d.SetParameters(myTolApprox, tol2d, iDegMin, iDegMax,
1583                                                   iNbIter, 30, Standard_True, aParType);
1584         }
1585       }
1586       //
1587       Standard_Real aReachedTol = Precision::Confusion();
1588       bIsDecomposited = IntTools_WLineTool::
1589         DecompositionOfWLine(WL,
1590                              myHS1,
1591                              myHS2,
1592                              myFace1,
1593                              myFace2,
1594                              myLConstruct,
1595                              bAvoidLineConstructor,
1596                              aSeqOfL,
1597                              aReachedTol,
1598                              myContext);
1599       if ( bIsDecomposited && ( myTolReached3d < aReachedTol ) ) {
1600         myTolReached3d = aReachedTol;
1601       }
1602       //
1603       aNbSeqOfL=aSeqOfL.Length();
1604       //
1605       if (bIsDecomposited) {
1606         nbiter=aNbSeqOfL;
1607       }
1608       else {
1609         nbiter=1;
1610         aNbParts=1;
1611         if (!bAvoidLineConstructor) {
1612           aNbParts=myLConstruct.NbParts();
1613           nbiter=aNbParts;
1614         }
1615       }
1616       //
1617       for(i = 1; i <= nbiter; ++i) {
1618         if(bIsDecomposited) {
1619           WL = Handle(IntPatch_WLine)::DownCast(aSeqOfL.Value(i));
1620           ifprm = 1;
1621           ilprm = WL->NbPnts();
1622         }
1623         else {
1624           if(bAvoidLineConstructor) {
1625             ifprm = 1;
1626             ilprm = WL->NbPnts();
1627           }
1628           else {
1629             myLConstruct.Part(i, fprm, lprm);
1630             ifprm = (Standard_Integer)fprm;
1631             ilprm = (Standard_Integer)lprm;
1632           }
1633         }
1634         //-- lbr :
1635         //-- Si une des surfaces est un plan , on approxime en 2d
1636         //-- sur cette surface et on remonte les points 2d en 3d.
1637         if(typs1 == GeomAbs_Plane) {
1638           theapp3d.Perform(myHS1, myHS2, WL, Standard_False,Standard_True, myApprox2,ifprm,ilprm);
1639         }
1640         else if(typs2 == GeomAbs_Plane) {
1641           theapp3d.Perform(myHS1,myHS2,WL,Standard_False,myApprox1,Standard_True,ifprm,ilprm);
1642         }
1643         else {
1644           //
1645           if (myHS1 != myHS2){
1646             if ((typs1==GeomAbs_BezierSurface || typs1==GeomAbs_BSplineSurface) &&
1647                 (typs2==GeomAbs_BezierSurface || typs2==GeomAbs_BSplineSurface)) {
1648
1649               theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, 30,
1650                                                                 Standard_True, aParType);
1651
1652               Standard_Boolean bUseSurfaces;
1653               bUseSurfaces = IntTools_WLineTool::NotUseSurfacesForApprox(myFace1, myFace2, WL, ifprm,  ilprm);
1654               if (bUseSurfaces) {
1655                 // ######
1656                 rejectSurface = Standard_True;
1657                 // ######
1658                 theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, 30,
1659                                                                 Standard_False, aParType);
1660               }
1661             }
1662           }
1663           //
1664           theapp3d.Perform(myHS1,myHS2,WL,Standard_True,anApprox1,anApprox2,ifprm,ilprm);
1665         }
1666           //
1667         if (!theapp3d.IsDone()) {
1668           Handle(Geom2d_BSplineCurve) H1;
1669           Handle(Geom2d_BSplineCurve) H2;
1670           //
1671           Handle(Geom_Curve) aBSp=GeomInt_IntSS::MakeBSpline(WL,ifprm, ilprm);
1672           //
1673           if(myApprox1) {
1674             H1 = GeomInt_IntSS::MakeBSpline2d(WL, ifprm, ilprm, Standard_True);
1675           }
1676           //
1677           if(myApprox2) {
1678             H2 = GeomInt_IntSS::MakeBSpline2d(WL, ifprm, ilprm, Standard_False);
1679           }
1680           //
1681           IntTools_Curve aIC(aBSp, H1, H2);
1682           mySeqOfCurve.Append(aIC);
1683         }
1684
1685         else {
1686           if(myApprox1 || myApprox2 || (typs1==GeomAbs_Plane || typs2==GeomAbs_Plane)) {
1687             if( theapp3d.TolReached2d()>myTolReached2d || myTolReached2d==0.) {
1688               myTolReached2d = theapp3d.TolReached2d();
1689             }
1690           }
1691           if(typs1==GeomAbs_Plane || typs2==GeomAbs_Plane) {
1692             myTolReached3d = myTolReached2d;
1693             //
1694             if (typs1==GeomAbs_Torus || typs2==GeomAbs_Torus) {
1695               if (myTolReached3d<1.e-6) {
1696                 myTolReached3d = theapp3d.TolReached3d();
1697                 myTolReached3d=1.e-6;
1698               }
1699             }
1700           }
1701           else  if( theapp3d.TolReached3d()>myTolReached3d || myTolReached3d==0.) {
1702             myTolReached3d = theapp3d.TolReached3d();
1703           }
1704
1705           Standard_Integer aNbMultiCurves, nbpoles;
1706           aNbMultiCurves=theapp3d.NbMultiCurves();
1707           for (j=1; j<=aNbMultiCurves; j++) {
1708             if(typs1 == GeomAbs_Plane) {
1709               const AppParCurves_MultiBSpCurve& mbspc = theapp3d.Value(j);
1710               nbpoles = mbspc.NbPoles();
1711
1712               TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
1713               TColgp_Array1OfPnt   tpoles(1,nbpoles);
1714
1715               mbspc.Curve(1,tpoles2d);
1716               const gp_Pln&  Pln = myHS1->Surface().Plane();
1717               //
1718               Standard_Integer ik;
1719               for(ik = 1; ik<= nbpoles; ik++) {
1720                 tpoles.SetValue(ik,
1721                                 ElSLib::Value(tpoles2d.Value(ik).X(),
1722                                               tpoles2d.Value(ik).Y(),
1723                                               Pln));
1724               }
1725               //
1726               Handle(Geom_BSplineCurve) BS =
1727                 new Geom_BSplineCurve(tpoles,
1728                                       mbspc.Knots(),
1729                                       mbspc.Multiplicities(),
1730                                       mbspc.Degree());
1731               GeomLib_CheckBSplineCurve Check(BS,TOLCHECK,TOLANGCHECK);
1732               Check.FixTangent(Standard_True, Standard_True);
1733               //
1734               IntTools_Curve aCurve;
1735               aCurve.SetCurve(BS);
1736
1737               if(myApprox1) {
1738                 Handle(Geom2d_BSplineCurve) BS1 =
1739                   new Geom2d_BSplineCurve(tpoles2d,
1740                                           mbspc.Knots(),
1741                                           mbspc.Multiplicities(),
1742                                           mbspc.Degree());
1743                 GeomLib_Check2dBSplineCurve Check1(BS1,TOLCHECK,TOLANGCHECK);
1744                 Check1.FixTangent(Standard_True,Standard_True);
1745                 //
1746                 // ############################################
1747                 if(!rejectSurface && !reApprox) {
1748                   Standard_Boolean isValid = IsCurveValid(BS1);
1749                   if(!isValid) {
1750                     reApprox = Standard_True;
1751                     goto reapprox;
1752                   }
1753                 }
1754                 // ############################################
1755                 aCurve.SetFirstCurve2d(BS1);
1756               }
1757               else {
1758                 Handle(Geom2d_BSplineCurve) H1;
1759                 aCurve.SetFirstCurve2d(H1);
1760               }
1761
1762               if(myApprox2) {
1763                 mbspc.Curve(2, tpoles2d);
1764
1765                 Handle(Geom2d_BSplineCurve) BS2 = new Geom2d_BSplineCurve(tpoles2d,
1766                                                                           mbspc.Knots(),
1767                                                                           mbspc.Multiplicities(),
1768                                                                           mbspc.Degree());
1769                 GeomLib_Check2dBSplineCurve newCheck(BS2,TOLCHECK,TOLANGCHECK);
1770                 newCheck.FixTangent(Standard_True,Standard_True);
1771
1772                 // ###########################################
1773                 if(!rejectSurface && !reApprox) {
1774                   Standard_Boolean isValid = IsCurveValid(BS2);
1775                   if(!isValid) {
1776                     reApprox = Standard_True;
1777                     goto reapprox;
1778                   }
1779                 }
1780                 // ###########################################
1781                 //
1782                 aCurve.SetSecondCurve2d(BS2);
1783               }
1784               else {
1785                 Handle(Geom2d_BSplineCurve) H2;
1786                 //
1787                   aCurve.SetSecondCurve2d(H2);
1788               }
1789               //
1790               mySeqOfCurve.Append(aCurve);
1791
1792             }//if(typs1 == GeomAbs_Plane) {
1793
1794             else if(typs2 == GeomAbs_Plane)
1795             {
1796               const AppParCurves_MultiBSpCurve& mbspc = theapp3d.Value(j);
1797               nbpoles = mbspc.NbPoles();
1798
1799               TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
1800               TColgp_Array1OfPnt   tpoles(1,nbpoles);
1801               mbspc.Curve((myApprox1==Standard_True)? 2 : 1,tpoles2d);
1802               const gp_Pln&  Pln = myHS2->Surface().Plane();
1803               //
1804               Standard_Integer ik;
1805               for(ik = 1; ik<= nbpoles; ik++) {
1806                 tpoles.SetValue(ik,
1807                                 ElSLib::Value(tpoles2d.Value(ik).X(),
1808                                               tpoles2d.Value(ik).Y(),
1809                                               Pln));
1810
1811               }
1812               //
1813               Handle(Geom_BSplineCurve) BS=new Geom_BSplineCurve(tpoles,
1814                                                                  mbspc.Knots(),
1815                                                                  mbspc.Multiplicities(),
1816                                                                  mbspc.Degree());
1817               GeomLib_CheckBSplineCurve Check(BS,TOLCHECK,TOLANGCHECK);
1818               Check.FixTangent(Standard_True,Standard_True);
1819               //
1820               IntTools_Curve aCurve;
1821               aCurve.SetCurve(BS);
1822
1823               if(myApprox2) {
1824                 Handle(Geom2d_BSplineCurve) BS1=new Geom2d_BSplineCurve(tpoles2d,
1825                                                                         mbspc.Knots(),
1826                                                                         mbspc.Multiplicities(),
1827                                                                         mbspc.Degree());
1828                 GeomLib_Check2dBSplineCurve Check1(BS1,TOLCHECK,TOLANGCHECK);
1829                 Check1.FixTangent(Standard_True,Standard_True);
1830                 //
1831                 // ###########################################
1832                 if(!rejectSurface && !reApprox) {
1833                   Standard_Boolean isValid = IsCurveValid(BS1);
1834                   if(!isValid) {
1835                     reApprox = Standard_True;
1836                     goto reapprox;
1837                   }
1838                 }
1839                 // ###########################################
1840                 bPCurvesOk = CheckPCurve(BS1, myFace2);
1841                 aCurve.SetSecondCurve2d(BS1);
1842               }
1843               else {
1844                 Handle(Geom2d_BSplineCurve) H2;
1845                 aCurve.SetSecondCurve2d(H2);
1846               }
1847
1848               if(myApprox1) {
1849                 mbspc.Curve(1,tpoles2d);
1850                 Handle(Geom2d_BSplineCurve) BS2=new Geom2d_BSplineCurve(tpoles2d,
1851                                                                         mbspc.Knots(),
1852                                                                         mbspc.Multiplicities(),
1853                                                                         mbspc.Degree());
1854                 GeomLib_Check2dBSplineCurve Check2(BS2,TOLCHECK,TOLANGCHECK);
1855                 Check2.FixTangent(Standard_True,Standard_True);
1856                 //
1857                 // ###########################################
1858                 if(!rejectSurface && !reApprox) {
1859                   Standard_Boolean isValid = IsCurveValid(BS2);
1860                   if(!isValid) {
1861                     reApprox = Standard_True;
1862                     goto reapprox;
1863                   }
1864                 }
1865                 // ###########################################
1866                 bPCurvesOk = bPCurvesOk && CheckPCurve(BS2, myFace1);
1867                 aCurve.SetFirstCurve2d(BS2);
1868               }
1869               else {
1870                 Handle(Geom2d_BSplineCurve) H1;
1871                 //
1872                 aCurve.SetFirstCurve2d(H1);
1873               }
1874               //
1875               //if points of the pcurves are out of the faces bounds
1876               //create 3d and 2d curves without approximation
1877               if (!bPCurvesOk) {
1878                 Handle(Geom2d_BSplineCurve) H1, H2;
1879                 bPCurvesOk = Standard_True;
1880                 //
1881                 Handle(Geom_Curve) aBSp=GeomInt_IntSS::MakeBSpline(WL,ifprm, ilprm);
1882
1883                 if(myApprox1) {
1884                   H1 = GeomInt_IntSS::MakeBSpline2d(WL, ifprm, ilprm, Standard_True);
1885                   bPCurvesOk = CheckPCurve(H1, myFace1);
1886                 }
1887
1888                 if(myApprox2) {
1889                   H2 = GeomInt_IntSS::MakeBSpline2d(WL, ifprm, ilprm, Standard_False);
1890                   bPCurvesOk = bPCurvesOk && CheckPCurve(H2, myFace2);
1891                 }
1892                 //
1893                 //if pcurves created without approximation are out of the
1894                 //faces bounds, use approximated 3d and 2d curves
1895                 if (bPCurvesOk) {
1896                   IntTools_Curve aIC(aBSp, H1, H2);
1897                   mySeqOfCurve.Append(aIC);
1898                 } else {
1899                   mySeqOfCurve.Append(aCurve);
1900                 }
1901               } else {
1902                 mySeqOfCurve.Append(aCurve);
1903               }
1904
1905             }// else if(typs2 == GeomAbs_Plane)
1906             //
1907             else { //typs2 != GeomAbs_Plane && typs1 != GeomAbs_Plane
1908               Standard_Boolean bIsValid1, bIsValid2;
1909               Handle(Geom_BSplineCurve) BS;
1910               Handle(Geom2d_BSplineCurve) aH2D;
1911               IntTools_Curve aCurve;
1912               //
1913               bIsValid1=Standard_True;
1914               bIsValid2=Standard_True;
1915               //
1916               const AppParCurves_MultiBSpCurve& mbspc = theapp3d.Value(j);
1917               nbpoles = mbspc.NbPoles();
1918               TColgp_Array1OfPnt tpoles(1,nbpoles);
1919               mbspc.Curve(1,tpoles);
1920               BS=new Geom_BSplineCurve(tpoles,
1921                                                                  mbspc.Knots(),
1922                                                                  mbspc.Multiplicities(),
1923                                                                  mbspc.Degree());
1924               GeomLib_CheckBSplineCurve Check(BS,TOLCHECK,TOLANGCHECK);
1925               Check.FixTangent(Standard_True,Standard_True);
1926               //
1927               aCurve.SetCurve(BS);
1928               aCurve.SetFirstCurve2d(aH2D);
1929               aCurve.SetSecondCurve2d(aH2D);
1930               //
1931               if(myApprox1) {
1932                 if(anApprox1) {
1933                   Handle(Geom2d_BSplineCurve) BS1;
1934                   TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
1935                   mbspc.Curve(2,tpoles2d);
1936                   //
1937                   BS1=new Geom2d_BSplineCurve(tpoles2d,
1938                                                                         mbspc.Knots(),
1939                                                                         mbspc.Multiplicities(),
1940                                                                         mbspc.Degree());
1941                   GeomLib_Check2dBSplineCurve newCheck(BS1,TOLCHECK,TOLANGCHECK);
1942                   newCheck.FixTangent(Standard_True,Standard_True);
1943                   //
1944                   if (!reApprox) {
1945                     bIsValid1=CheckPCurve(BS1, myFace1);
1946                   }
1947                   //
1948                   aCurve.SetFirstCurve2d(BS1);
1949                 }
1950                 else {
1951                   Handle(Geom2d_BSplineCurve) BS1;
1952                   fprm = BS->FirstParameter();
1953                   lprm = BS->LastParameter();
1954
1955                   Handle(Geom2d_Curve) C2d;
1956                   Standard_Real aTol = myTolApprox;
1957                   GeomInt_IntSS::BuildPCurves(fprm, lprm, aTol,
1958                             myHS1->ChangeSurface().Surface(), BS, C2d);
1959                   BS1 = Handle(Geom2d_BSplineCurve)::DownCast(C2d);
1960                   aCurve.SetFirstCurve2d(BS1);
1961                 }
1962               } // if(myApprox1) {
1963                 //
1964               if(myApprox2) {
1965                 if(anApprox2) {
1966                   Handle(Geom2d_BSplineCurve) BS2;
1967                   TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
1968                   mbspc.Curve((myApprox1==Standard_True)? 3 : 2,tpoles2d);
1969                   BS2=new Geom2d_BSplineCurve(tpoles2d,
1970                                                                         mbspc.Knots(),
1971                                                                         mbspc.Multiplicities(),
1972                                                                         mbspc.Degree());
1973                   GeomLib_Check2dBSplineCurve newCheck(BS2,TOLCHECK,TOLANGCHECK);
1974                   newCheck.FixTangent(Standard_True,Standard_True);
1975                 //
1976                   if (!reApprox) {
1977                     bIsValid2=CheckPCurve(BS2, myFace2);
1978                   }
1979                   aCurve.SetSecondCurve2d(BS2);
1980                 }
1981                 else {
1982                   Handle(Geom2d_BSplineCurve) BS2;
1983                   fprm = BS->FirstParameter();
1984                   lprm = BS->LastParameter();
1985
1986                   Handle(Geom2d_Curve) C2d;
1987                   Standard_Real aTol = myTolApprox;
1988                   GeomInt_IntSS::BuildPCurves(fprm, lprm, aTol,
1989                             myHS2->ChangeSurface().Surface(), BS, C2d);
1990                   BS2 = Handle(Geom2d_BSplineCurve)::DownCast(C2d);
1991                   aCurve.SetSecondCurve2d(BS2);
1992                 }
1993               } //if(myApprox2) {
1994               if (!bIsValid1 || !bIsValid2) {
1995                 myTolApprox=aTolApproxImp;//1.e-5;
1996                 tol2d = myTolApprox;
1997                 reApprox = Standard_True;
1998                 goto reapprox;
1999               }
2000                 //
2001               mySeqOfCurve.Append(aCurve);
2002             }
2003           }
2004         }
2005       }
2006     }// else { // X
2007   }// case IntPatch_Walking:{
2008     break;
2009
2010   case IntPatch_Restriction:
2011     {
2012       Handle(IntPatch_RLine) RL =
2013         Handle(IntPatch_RLine)::DownCast(L);
2014       Handle(Geom_Curve) aC3d;
2015       Handle(Geom2d_Curve) aC2d1, aC2d2;
2016       Standard_Real aTolReached;
2017       GeomInt_IntSS::TreatRLine(RL, myHS1, myHS2, aC3d,
2018                                   aC2d1, aC2d2, aTolReached);
2019
2020       if(aC3d.IsNull())
2021         break;
2022
2023       Bnd_Box2d aBox1, aBox2;
2024
2025       const Standard_Real aU1f = myHS1->FirstUParameter(),
2026                           aV1f = myHS1->FirstVParameter(),
2027                           aU1l = myHS1->LastUParameter(),
2028                           aV1l = myHS1->LastVParameter();
2029       const Standard_Real aU2f = myHS2->FirstUParameter(),
2030                           aV2f = myHS2->FirstVParameter(),
2031                           aU2l = myHS2->LastUParameter(),
2032                           aV2l = myHS2->LastVParameter();
2033
2034       aBox1.Add(gp_Pnt2d(aU1f, aV1f));
2035       aBox1.Add(gp_Pnt2d(aU1l, aV1l));
2036       aBox2.Add(gp_Pnt2d(aU2f, aV2f));
2037       aBox2.Add(gp_Pnt2d(aU2l, aV2l));
2038
2039       GeomInt_VectorOfReal anArrayOfParameters;
2040
2041       //We consider here that the intersection line is same-parameter-line
2042       anArrayOfParameters.Append(aC3d->FirstParameter());
2043       anArrayOfParameters.Append(aC3d->LastParameter());
2044
2045       GeomInt_IntSS::
2046         TrimILineOnSurfBoundaries(aC2d1, aC2d2, aBox1, aBox2, anArrayOfParameters);
2047
2048       //Intersect with true boundaries. After that, enlarge bounding-boxes in order to
2049       //correct definition, if point on curve is inscribed in the box.
2050       aBox1.Enlarge(theToler);
2051       aBox2.Enlarge(theToler);
2052
2053       const Standard_Integer aNbIntersSolutionsm1 = anArrayOfParameters.Length() - 1;
2054
2055       //Trim RLine found.
2056       for(Standard_Integer anInd = 0; anInd < aNbIntersSolutionsm1; anInd++)
2057       {
2058         Standard_Real &aParF = anArrayOfParameters(anInd),
2059                       &aParL = anArrayOfParameters(anInd+1);
2060
2061         if((aParL - aParF) <= Precision::PConfusion())
2062         {
2063           //In order to more precise extending to the boundaries of source curves.
2064           if(anInd < aNbIntersSolutionsm1-1)
2065             aParL = aParF;
2066
2067           continue;
2068         }
2069
2070         const Standard_Real aPar = 0.5*(aParF + aParL);
2071         gp_Pnt2d aPt;
2072
2073         Handle(Geom2d_Curve) aCurv2d1, aCurv2d2;
2074         if(!aC2d1.IsNull())
2075         {
2076           aC2d1->D0(aPar, aPt);
2077
2078           if(aBox1.IsOut(aPt))
2079             continue;
2080
2081           if(myApprox1)
2082             aCurv2d1 = new Geom2d_TrimmedCurve(aC2d1, aParF, aParL);
2083         }
2084
2085         if(!aC2d2.IsNull())
2086         {
2087           aC2d2->D0(aPar, aPt);
2088
2089           if(aBox2.IsOut(aPt))
2090             continue;
2091
2092           if(myApprox2)
2093             aCurv2d2 = new Geom2d_TrimmedCurve(aC2d2, aParF, aParL);
2094         }
2095
2096         Handle(Geom_Curve) aCurv3d = new Geom_TrimmedCurve(aC3d, aParF, aParL);
2097
2098         IntTools_Curve aIC(aCurv3d, aCurv2d1, aCurv2d2);
2099         mySeqOfCurve.Append(aIC);
2100       }
2101     }
2102     break;
2103   default:
2104     break;
2105
2106   }
2107 }
2108
2109 //=======================================================================
2110 //function : Parameters
2111 //purpose  :
2112 //=======================================================================
2113  void Parameters(const Handle(GeomAdaptor_HSurface)& HS1,
2114                  const Handle(GeomAdaptor_HSurface)& HS2,
2115                  const gp_Pnt& Ptref,
2116                  Standard_Real& U1,
2117                  Standard_Real& V1,
2118                  Standard_Real& U2,
2119                  Standard_Real& V2)
2120 {
2121
2122   IntSurf_Quadric quad1,quad2;
2123   GeomAbs_SurfaceType typs = HS1->Surface().GetType();
2124
2125   switch (typs) {
2126   case GeomAbs_Plane:
2127     quad1.SetValue(HS1->Surface().Plane());
2128     break;
2129   case GeomAbs_Cylinder:
2130     quad1.SetValue(HS1->Surface().Cylinder());
2131     break;
2132   case GeomAbs_Cone:
2133     quad1.SetValue(HS1->Surface().Cone());
2134     break;
2135   case GeomAbs_Sphere:
2136     quad1.SetValue(HS1->Surface().Sphere());
2137     break;
2138   case GeomAbs_Torus:
2139     quad1.SetValue(HS1->Surface().Torus());
2140     break;
2141   default:
2142     Standard_ConstructionError::Raise("GeomInt_IntSS::MakeCurve");
2143   }
2144
2145   typs = HS2->Surface().GetType();
2146   switch (typs) {
2147   case GeomAbs_Plane:
2148     quad2.SetValue(HS2->Surface().Plane());
2149     break;
2150   case GeomAbs_Cylinder:
2151     quad2.SetValue(HS2->Surface().Cylinder());
2152     break;
2153   case GeomAbs_Cone:
2154     quad2.SetValue(HS2->Surface().Cone());
2155     break;
2156   case GeomAbs_Sphere:
2157     quad2.SetValue(HS2->Surface().Sphere());
2158     break;
2159   case GeomAbs_Torus:
2160     quad2.SetValue(HS2->Surface().Torus());
2161     break;
2162   default:
2163     Standard_ConstructionError::Raise("GeomInt_IntSS::MakeCurve");
2164   }
2165
2166   quad1.Parameters(Ptref,U1,V1);
2167   quad2.Parameters(Ptref,U2,V2);
2168 }
2169
2170 //=======================================================================
2171 //function : MakeBSpline
2172 //purpose  :
2173 //=======================================================================
2174 Handle(Geom_Curve) MakeBSpline  (const Handle(IntPatch_WLine)& WL,
2175                                  const Standard_Integer ideb,
2176                                  const Standard_Integer ifin)
2177 {
2178   Standard_Integer i,nbpnt = ifin-ideb+1;
2179   TColgp_Array1OfPnt poles(1,nbpnt);
2180   TColStd_Array1OfReal knots(1,nbpnt);
2181   TColStd_Array1OfInteger mults(1,nbpnt);
2182   Standard_Integer ipidebm1;
2183   for(i=1,ipidebm1=i+ideb-1; i<=nbpnt;ipidebm1++, i++) {
2184     poles(i) = WL->Point(ipidebm1).Value();
2185     mults(i) = 1;
2186     knots(i) = i-1;
2187   }
2188   mults(1) = mults(nbpnt) = 2;
2189   return
2190     new Geom_BSplineCurve(poles,knots,mults,1);
2191 }
2192
2193 //=======================================================================
2194 //function : PrepareLines3D
2195 //purpose  :
2196 //=======================================================================
2197   void IntTools_FaceFace::PrepareLines3D(const Standard_Boolean bToSplit)
2198 {
2199   Standard_Integer i, aNbCurves;
2200   GeomAbs_SurfaceType aType1, aType2;
2201   IntTools_SequenceOfCurves aNewCvs;
2202   //
2203   // 1. Treatment closed  curves
2204   aNbCurves=mySeqOfCurve.Length();
2205   for (i=1; i<=aNbCurves; ++i) {
2206     const IntTools_Curve& aIC=mySeqOfCurve(i);
2207     //
2208     if (bToSplit) {
2209       Standard_Integer j, aNbC;
2210       IntTools_SequenceOfCurves aSeqCvs;
2211       //
2212       aNbC=IntTools_Tools::SplitCurve(aIC, aSeqCvs);
2213       if (aNbC) {
2214         for (j=1; j<=aNbC; ++j) {
2215           const IntTools_Curve& aICNew=aSeqCvs(j);
2216           aNewCvs.Append(aICNew);
2217         }
2218       }
2219       else {
2220         aNewCvs.Append(aIC);
2221       }
2222     }
2223     else {
2224       aNewCvs.Append(aIC);
2225     }
2226   }
2227   //
2228   // 2. Plane\Cone intersection when we had 4 curves
2229   aType1=myHS1->GetType();
2230   aType2=myHS2->GetType();
2231   aNbCurves=aNewCvs.Length();
2232   //
2233   if ((aType1==GeomAbs_Plane && aType2==GeomAbs_Cone) ||
2234       (aType2==GeomAbs_Plane && aType1==GeomAbs_Cone)) {
2235     if (aNbCurves==4) {
2236       GeomAbs_CurveType aCType1;
2237       //
2238       aCType1=aNewCvs(1).Type();
2239       if (aCType1==GeomAbs_Line) {
2240         IntTools_SequenceOfCurves aSeqIn, aSeqOut;
2241         //
2242         for (i=1; i<=aNbCurves; ++i) {
2243           const IntTools_Curve& aIC=aNewCvs(i);
2244           aSeqIn.Append(aIC);
2245         }
2246         //
2247         IntTools_Tools::RejectLines(aSeqIn, aSeqOut);
2248         //
2249         aNewCvs.Clear();
2250         aNbCurves=aSeqOut.Length();
2251         for (i=1; i<=aNbCurves; ++i) {
2252           const IntTools_Curve& aIC=aSeqOut(i);
2253           aNewCvs.Append(aIC);
2254         }
2255       }
2256     }
2257   }// if ((aType1==GeomAbs_Plane && aType2==GeomAbs_Cone)...
2258   //
2259   // 3. Fill  mySeqOfCurve
2260   mySeqOfCurve.Clear();
2261   aNbCurves=aNewCvs.Length();
2262   for (i=1; i<=aNbCurves; ++i) {
2263     const IntTools_Curve& aIC=aNewCvs(i);
2264     mySeqOfCurve.Append(aIC);
2265   }
2266 }
2267 //=======================================================================
2268 //function : CorrectSurfaceBoundaries
2269 //purpose  :
2270 //=======================================================================
2271  void CorrectSurfaceBoundaries(const TopoDS_Face&  theFace,
2272                               const Standard_Real theTolerance,
2273                               Standard_Real&      theumin,
2274                               Standard_Real&      theumax,
2275                               Standard_Real&      thevmin,
2276                               Standard_Real&      thevmax)
2277 {
2278   Standard_Boolean enlarge, isuperiodic, isvperiodic;
2279   Standard_Real uinf, usup, vinf, vsup, delta;
2280   GeomAbs_SurfaceType aType;
2281   Handle(Geom_Surface) aSurface;
2282   //
2283   aSurface = BRep_Tool::Surface(theFace);
2284   aSurface->Bounds(uinf, usup, vinf, vsup);
2285   delta = theTolerance;
2286   enlarge = Standard_False;
2287   //
2288   GeomAdaptor_Surface anAdaptorSurface(aSurface);
2289   //
2290   if(aSurface->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface))) {
2291     Handle(Geom_Surface) aBasisSurface =
2292       (Handle(Geom_RectangularTrimmedSurface)::DownCast(aSurface))->BasisSurface();
2293
2294     if(aBasisSurface->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface)) ||
2295        aBasisSurface->IsKind(STANDARD_TYPE(Geom_OffsetSurface))) {
2296       return;
2297     }
2298   }
2299   //
2300   if(aSurface->IsKind(STANDARD_TYPE(Geom_OffsetSurface))) {
2301     Handle(Geom_Surface) aBasisSurface =
2302       (Handle(Geom_OffsetSurface)::DownCast(aSurface))->BasisSurface();
2303
2304     if(aBasisSurface->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface)) ||
2305        aBasisSurface->IsKind(STANDARD_TYPE(Geom_OffsetSurface))) {
2306       return;
2307     }
2308   }
2309   //
2310   isuperiodic = anAdaptorSurface.IsUPeriodic();
2311   isvperiodic = anAdaptorSurface.IsVPeriodic();
2312   //
2313   aType=anAdaptorSurface.GetType();
2314   if((aType==GeomAbs_BezierSurface) ||
2315      (aType==GeomAbs_BSplineSurface) ||
2316      (aType==GeomAbs_SurfaceOfExtrusion) ||
2317      (aType==GeomAbs_SurfaceOfRevolution) ||
2318      (aType==GeomAbs_Cylinder)) {
2319     enlarge=Standard_True;
2320   }
2321   //
2322   if(!isuperiodic && enlarge) {
2323
2324     if(!Precision::IsInfinite(theumin) &&
2325         ((theumin - uinf) > delta))
2326       theumin -= delta;
2327     else {
2328       theumin = uinf;
2329     }
2330
2331     if(!Precision::IsInfinite(theumax) &&
2332         ((usup - theumax) > delta))
2333       theumax += delta;
2334     else
2335       theumax = usup;
2336   }
2337   //
2338   if(!isvperiodic && enlarge) {
2339     if(!Precision::IsInfinite(thevmin) &&
2340         ((thevmin - vinf) > delta)) {
2341       thevmin -= delta;
2342     }
2343     else {
2344       thevmin = vinf;
2345     }
2346     if(!Precision::IsInfinite(thevmax) &&
2347         ((vsup - thevmax) > delta)) {
2348       thevmax += delta;
2349     }
2350     else {
2351       thevmax = vsup;
2352     }
2353   }
2354   //
2355   {
2356     Standard_Integer aNbP;
2357     Standard_Real aXP, dXfact, aXmid, aX1, aX2, aTolPA;
2358     //
2359     aTolPA=Precision::Angular();
2360     // U
2361     if (isuperiodic) {
2362       aXP=anAdaptorSurface.UPeriod();
2363       dXfact=theumax-theumin;
2364       if (dXfact-aTolPA>aXP) {
2365         aXmid=0.5*(theumax+theumin);
2366         aNbP=RealToInt(aXmid/aXP);
2367         if (aXmid<0.) {
2368           aNbP=aNbP-1;
2369         }
2370         aX1=aNbP*aXP;
2371         if (theumin>aTolPA) {
2372           aX1=theumin+aNbP*aXP;
2373         }
2374         aX2=aX1+aXP;
2375         if (theumin<aX1) {
2376           theumin=aX1;
2377         }
2378         if (theumax>aX2) {
2379           theumax=aX2;
2380         }
2381       }
2382     }
2383     // V
2384     if (isvperiodic) {
2385       aXP=anAdaptorSurface.VPeriod();
2386       dXfact=thevmax-thevmin;
2387       if (dXfact-aTolPA>aXP) {
2388         aXmid=0.5*(thevmax+thevmin);
2389         aNbP=RealToInt(aXmid/aXP);
2390         if (aXmid<0.) {
2391           aNbP=aNbP-1;
2392         }
2393         aX1=aNbP*aXP;
2394         if (thevmin>aTolPA) {
2395           aX1=thevmin+aNbP*aXP;
2396         }
2397         aX2=aX1+aXP;
2398         if (thevmin<aX1) {
2399           thevmin=aX1;
2400         }
2401         if (thevmax>aX2) {
2402           thevmax=aX2;
2403         }
2404       }
2405     }
2406   }
2407   //
2408   if(isuperiodic || isvperiodic) {
2409     Standard_Boolean correct = Standard_False;
2410     Standard_Boolean correctU = Standard_False;
2411     Standard_Boolean correctV = Standard_False;
2412     Bnd_Box2d aBox;
2413     TopExp_Explorer anExp;
2414
2415     for(anExp.Init(theFace, TopAbs_EDGE); anExp.More(); anExp.Next()) {
2416       if(BRep_Tool::IsClosed(TopoDS::Edge(anExp.Current()), theFace)) {
2417         correct = Standard_True;
2418         Standard_Real f, l;
2419         TopoDS_Edge anEdge = TopoDS::Edge(anExp.Current());
2420
2421         for(Standard_Integer i = 0; i < 2; i++) {
2422           if(i==0) {
2423             anEdge.Orientation(TopAbs_FORWARD);
2424           }
2425           else {
2426             anEdge.Orientation(TopAbs_REVERSED);
2427           }
2428           Handle(Geom2d_Curve) aCurve = BRep_Tool::CurveOnSurface(anEdge, theFace, f, l);
2429
2430           if(aCurve.IsNull()) {
2431             correct = Standard_False;
2432             break;
2433           }
2434           Handle(Geom2d_Line) aLine = Handle(Geom2d_Line)::DownCast(aCurve);
2435
2436           if(aLine.IsNull()) {
2437             correct = Standard_False;
2438             break;
2439           }
2440           gp_Dir2d anUDir(1., 0.);
2441           gp_Dir2d aVDir(0., 1.);
2442           Standard_Real anAngularTolerance = Precision::Angular();
2443
2444           correctU = correctU || aLine->Position().Direction().IsParallel(aVDir, anAngularTolerance);
2445           correctV = correctV || aLine->Position().Direction().IsParallel(anUDir, anAngularTolerance);
2446
2447           gp_Pnt2d pp1 = aCurve->Value(f);
2448           aBox.Add(pp1);
2449           gp_Pnt2d pp2 = aCurve->Value(l);
2450           aBox.Add(pp2);
2451         }
2452         if(!correct)
2453           break;
2454       }
2455     }
2456
2457     if(correct) {
2458       Standard_Real umin, vmin, umax, vmax;
2459       aBox.Get(umin, vmin, umax, vmax);
2460
2461       if(isuperiodic && correctU) {
2462
2463         if(theumin < umin)
2464           theumin = umin;
2465
2466         if(theumax > umax) {
2467           theumax = umax;
2468         }
2469       }
2470       if(isvperiodic && correctV) {
2471
2472         if(thevmin < vmin)
2473           thevmin = vmin;
2474         if(thevmax > vmax)
2475           thevmax = vmax;
2476       }
2477     }
2478   }
2479 }
2480
2481 //=======================================================================
2482 //function : TolR3d
2483 //purpose  :
2484 //=======================================================================
2485 void TolR3d(const TopoDS_Face& aF1,
2486             const TopoDS_Face& aF2,
2487             Standard_Real& myTolReached3d)
2488 {
2489   Standard_Real aTolF1, aTolF2, aTolFMax, aTolTresh;
2490
2491   aTolTresh=2.999999e-3;
2492   aTolF1 = BRep_Tool::Tolerance(aF1);
2493   aTolF2 = BRep_Tool::Tolerance(aF2);
2494   aTolFMax=Max(aTolF1, aTolF2);
2495
2496   if (aTolFMax>aTolTresh) {
2497     myTolReached3d=aTolFMax;
2498   }
2499 }
2500
2501 // ------------------------------------------------------------------------------------------------
2502 // static function: ParameterOutOfBoundary
2503 // purpose:         Computes a new parameter for given curve. The corresponding 2d points
2504 //                  does not lay on any boundary of given faces
2505 // ------------------------------------------------------------------------------------------------
2506 Standard_Boolean ParameterOutOfBoundary(const Standard_Real       theParameter,
2507                                         const Handle(Geom_Curve)& theCurve,
2508                                         const TopoDS_Face&        theFace1,
2509                                         const TopoDS_Face&        theFace2,
2510                                         const Standard_Real       theOtherParameter,
2511                                         const Standard_Boolean    bIncreasePar,
2512                                         Standard_Real&            theNewParameter,
2513                                         const Handle(IntTools_Context)& aContext)
2514 {
2515   Standard_Boolean bIsComputed = Standard_False;
2516   theNewParameter = theParameter;
2517
2518   Standard_Real acurpar = theParameter;
2519   TopAbs_State aState = TopAbs_ON;
2520   Standard_Integer iter = 0;
2521   Standard_Real asumtol = BRep_Tool::Tolerance(theFace1) + BRep_Tool::Tolerance(theFace2);
2522   Standard_Real adelta = asumtol * 0.1;
2523   adelta = (adelta < Precision::Confusion()) ? Precision::Confusion() : adelta;
2524   Handle(Geom_Surface) aSurf1 = BRep_Tool::Surface(theFace1);
2525   Handle(Geom_Surface) aSurf2 = BRep_Tool::Surface(theFace2);
2526
2527   Standard_Real u1, u2, v1, v2;
2528
2529   GeomAPI_ProjectPointOnSurf aPrj1;
2530   aSurf1->Bounds(u1, u2, v1, v2);
2531   aPrj1.Init(aSurf1, u1, u2, v1, v2);
2532
2533   GeomAPI_ProjectPointOnSurf aPrj2;
2534   aSurf2->Bounds(u1, u2, v1, v2);
2535   aPrj2.Init(aSurf2, u1, u2, v1, v2);
2536
2537   while(aState == TopAbs_ON) {
2538     if(bIncreasePar)
2539       acurpar += adelta;
2540     else
2541       acurpar -= adelta;
2542     gp_Pnt aPCurrent = theCurve->Value(acurpar);
2543     aPrj1.Perform(aPCurrent);
2544     Standard_Real U=0., V=0.;
2545
2546     if(aPrj1.IsDone()) {
2547       aPrj1.LowerDistanceParameters(U, V);
2548       aState = aContext->StatePointFace(theFace1, gp_Pnt2d(U, V));
2549     }
2550
2551     if(aState != TopAbs_ON) {
2552       aPrj2.Perform(aPCurrent);
2553
2554       if(aPrj2.IsDone()) {
2555         aPrj2.LowerDistanceParameters(U, V);
2556         aState = aContext->StatePointFace(theFace2, gp_Pnt2d(U, V));
2557       }
2558     }
2559
2560     if(iter > 11) {
2561       break;
2562     }
2563     iter++;
2564   }
2565
2566   if(iter <= 11) {
2567     theNewParameter = acurpar;
2568     bIsComputed = Standard_True;
2569
2570     if(bIncreasePar) {
2571       if(acurpar >= theOtherParameter)
2572         theNewParameter = theOtherParameter;
2573     }
2574     else {
2575       if(acurpar <= theOtherParameter)
2576         theNewParameter = theOtherParameter;
2577     }
2578   }
2579   return bIsComputed;
2580 }
2581
2582 //=======================================================================
2583 //function : IsCurveValid
2584 //purpose  :
2585 //=======================================================================
2586 Standard_Boolean IsCurveValid(Handle(Geom2d_Curve)& thePCurve)
2587 {
2588   if(thePCurve.IsNull())
2589     return Standard_False;
2590
2591   Standard_Real tolint = 1.e-10;
2592   Geom2dAdaptor_Curve PCA;
2593   IntRes2d_Domain PCD;
2594   Geom2dInt_GInter PCI;
2595
2596   Standard_Real pf = 0., pl = 0.;
2597   gp_Pnt2d pntf, pntl;
2598
2599   if(!thePCurve->IsClosed() && !thePCurve->IsPeriodic()) {
2600     pf = thePCurve->FirstParameter();
2601     pl = thePCurve->LastParameter();
2602     pntf = thePCurve->Value(pf);
2603     pntl = thePCurve->Value(pl);
2604     PCA.Load(thePCurve);
2605     if(!PCA.IsPeriodic()) {
2606       if(PCA.FirstParameter() > pf) pf = PCA.FirstParameter();
2607       if(PCA.LastParameter()  < pl) pl = PCA.LastParameter();
2608     }
2609     PCD.SetValues(pntf,pf,tolint,pntl,pl,tolint);
2610     PCI.Perform(PCA,PCD,tolint,tolint);
2611     if(PCI.IsDone())
2612       if(PCI.NbPoints() > 0) {
2613         return Standard_False;
2614       }
2615   }
2616
2617   return Standard_True;
2618 }
2619
2620 //=======================================================================
2621 //static function : ApproxWithPCurves
2622 //purpose  : for bug 20964 only
2623 //=======================================================================
2624 Standard_Boolean ApproxWithPCurves(const gp_Cylinder& theCyl,
2625                                    const gp_Sphere& theSph)
2626 {
2627   Standard_Boolean bRes = Standard_True;
2628   Standard_Real R1 = theCyl.Radius(), R2 = theSph.Radius();
2629   //
2630   {
2631     Standard_Real aD2, aRc2, aEps;
2632     gp_Pnt aApexSph;
2633     //
2634     aEps=1.E-7;
2635     aRc2=R1*R1;
2636     //
2637     const gp_Ax3& aAx3Sph=theSph.Position();
2638     const gp_Pnt& aLocSph=aAx3Sph.Location();
2639     const gp_Dir& aDirSph=aAx3Sph.Direction();
2640     //
2641     const gp_Ax1& aAx1Cyl=theCyl.Axis();
2642     gp_Lin aLinCyl(aAx1Cyl);
2643     //
2644     aApexSph.SetXYZ(aLocSph.XYZ()+R2*aDirSph.XYZ());
2645     aD2=aLinCyl.SquareDistance(aApexSph);
2646     if (fabs(aD2-aRc2)<aEps) {
2647       return !bRes;
2648     }
2649     //
2650     aApexSph.SetXYZ(aLocSph.XYZ()-R2*aDirSph.XYZ());
2651     aD2=aLinCyl.SquareDistance(aApexSph);
2652     if (fabs(aD2-aRc2)<aEps) {
2653       return !bRes;
2654     }
2655   }
2656   //
2657
2658   if(R1 < 2.*R2) {
2659     return bRes;
2660   }
2661   gp_Lin anCylAx(theCyl.Axis());
2662
2663   Standard_Real aDist = anCylAx.Distance(theSph.Location());
2664   Standard_Real aDRel = Abs(aDist - R1)/R2;
2665
2666   if(aDRel > .2) return bRes;
2667
2668   Standard_Real par = ElCLib::Parameter(anCylAx, theSph.Location());
2669   gp_Pnt aP = ElCLib::Value(par, anCylAx);
2670   gp_Vec aV(aP, theSph.Location());
2671
2672   Standard_Real dd = aV.Dot(theSph.Position().XDirection());
2673
2674   if(aDist < R1 && dd > 0.) return Standard_False;
2675   if(aDist > R1 && dd < 0.) return Standard_False;
2676
2677
2678   return bRes;
2679 }
2680 //=======================================================================
2681 //function : PerformPlanes
2682 //purpose  :
2683 //=======================================================================
2684 void  PerformPlanes(const Handle(GeomAdaptor_HSurface)& theS1,
2685                     const Handle(GeomAdaptor_HSurface)& theS2,
2686                     const Standard_Real TolAng,
2687                     const Standard_Real TolTang,
2688                     const Standard_Boolean theApprox1,
2689                     const Standard_Boolean theApprox2,
2690                     IntTools_SequenceOfCurves& theSeqOfCurve,
2691                     Standard_Boolean& theTangentFaces)
2692 {
2693
2694   gp_Pln aPln1 = theS1->Surface().Plane();
2695   gp_Pln aPln2 = theS2->Surface().Plane();
2696
2697   IntAna_QuadQuadGeo aPlnInter(aPln1, aPln2, TolAng, TolTang);
2698
2699   if(!aPlnInter.IsDone()) {
2700     theTangentFaces = Standard_False;
2701     return;
2702   }
2703
2704   IntAna_ResultType aResType = aPlnInter.TypeInter();
2705
2706   if(aResType == IntAna_Same) {
2707     theTangentFaces = Standard_True;
2708     return;
2709   }
2710
2711   theTangentFaces = Standard_False;
2712
2713   if(aResType == IntAna_Empty) {
2714     return;
2715   }
2716
2717   gp_Lin aLin = aPlnInter.Line(1);
2718
2719   ProjLib_Plane aProj;
2720
2721   aProj.Init(aPln1);
2722   aProj.Project(aLin);
2723   gp_Lin2d aLin2d1 = aProj.Line();
2724   //
2725   aProj.Init(aPln2);
2726   aProj.Project(aLin);
2727   gp_Lin2d aLin2d2 = aProj.Line();
2728   //
2729   //classify line2d1 relatively first plane
2730   Standard_Real P11, P12;
2731   Standard_Boolean IsCrossed = ClassifyLin2d(theS1, aLin2d1, TolTang, P11, P12);
2732   if(!IsCrossed) return;
2733   //classify line2d2 relatively second plane
2734   Standard_Real P21, P22;
2735   IsCrossed = ClassifyLin2d(theS2, aLin2d2, TolTang, P21, P22);
2736   if(!IsCrossed) return;
2737
2738   //Analysis of parametric intervals: must have common part
2739
2740   if(P21 >= P12) return;
2741   if(P22 <= P11) return;
2742
2743   Standard_Real pmin, pmax;
2744   pmin = Max(P11, P21);
2745   pmax = Min(P12, P22);
2746
2747   if(pmax - pmin <= TolTang) return;
2748
2749   Handle(Geom_Line) aGLin = new Geom_Line(aLin);
2750
2751   IntTools_Curve aCurve;
2752   Handle(Geom_TrimmedCurve) aGTLin = new Geom_TrimmedCurve(aGLin, pmin, pmax);
2753
2754   aCurve.SetCurve(aGTLin);
2755
2756   if(theApprox1) {
2757     Handle(Geom2d_Line) C2d = new Geom2d_Line(aLin2d1);
2758     aCurve.SetFirstCurve2d(new Geom2d_TrimmedCurve(C2d, pmin, pmax));
2759   }
2760   else {
2761     Handle(Geom2d_Curve) H1;
2762     aCurve.SetFirstCurve2d(H1);
2763   }
2764   if(theApprox2) {
2765     Handle(Geom2d_Line) C2d = new Geom2d_Line(aLin2d2);
2766     aCurve.SetSecondCurve2d(new Geom2d_TrimmedCurve(C2d, pmin, pmax));
2767   }
2768   else {
2769     Handle(Geom2d_Curve) H1;
2770     aCurve.SetFirstCurve2d(H1);
2771   }
2772
2773   theSeqOfCurve.Append(aCurve);
2774
2775 }
2776
2777 //=======================================================================
2778 //function : ClassifyLin2d
2779 //purpose  :
2780 //=======================================================================
2781 static inline Standard_Boolean INTER(const Standard_Real d1,
2782                                      const Standard_Real d2,
2783                                      const Standard_Real tol)
2784 {
2785   return (d1 > tol && d2 < -tol) ||
2786          (d1 < -tol && d2 > tol) ||
2787          ((d1 <= tol && d1 >= -tol) && (d2 > tol || d2 < -tol)) ||
2788          ((d2 <= tol && d2 >= -tol) && (d1 > tol || d1 < -tol));
2789 }
2790 static inline Standard_Boolean COINC(const Standard_Real d1,
2791                                      const Standard_Real d2,
2792                                      const Standard_Real tol)
2793 {
2794   return (d1 <= tol && d1 >= -tol) && (d2 <= tol && d2 >= -tol);
2795 }
2796 Standard_Boolean ClassifyLin2d(const Handle(GeomAdaptor_HSurface)& theS,
2797                                const gp_Lin2d& theLin2d,
2798                                const Standard_Real theTol,
2799                                Standard_Real& theP1,
2800                                Standard_Real& theP2)
2801
2802 {
2803   Standard_Real xmin, xmax, ymin, ymax, d1, d2, A, B, C;
2804   Standard_Real par[2];
2805   Standard_Integer nbi = 0;
2806
2807   xmin = theS->Surface().FirstUParameter();
2808   xmax = theS->Surface().LastUParameter();
2809   ymin = theS->Surface().FirstVParameter();
2810   ymax = theS->Surface().LastVParameter();
2811
2812   theLin2d.Coefficients(A, B, C);
2813
2814   //xmin, ymin <-> xmin, ymax
2815   d1 = A*xmin + B*ymin + C;
2816   d2 = A*xmin + B*ymax + C;
2817
2818   if(INTER(d1, d2, theTol)) {
2819     //Intersection with boundary
2820     Standard_Real y = -(C + A*xmin)/B;
2821     par[nbi] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmin, y));
2822     nbi++;
2823   }
2824   else if (COINC(d1, d2, theTol)) {
2825     //Coincidence with boundary
2826     par[0] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmin, ymin));
2827     par[1] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmin, ymax));
2828     nbi = 2;
2829   }
2830
2831   if(nbi == 2) {
2832
2833     if(fabs(par[0]-par[1]) > theTol) {
2834       theP1 = Min(par[0], par[1]);
2835       theP2 = Max(par[0], par[1]);
2836       return Standard_True;
2837     }
2838     else return Standard_False;
2839
2840   }
2841
2842   //xmin, ymax <-> xmax, ymax
2843   d1 = d2;
2844   d2 = A*xmax + B*ymax + C;
2845
2846   if(d1 > theTol || d1 < -theTol) {//to avoid checking of
2847                                    //coincidence with the same point
2848     if(INTER(d1, d2, theTol)) {
2849       Standard_Real x = -(C + B*ymax)/A;
2850       par[nbi] = ElCLib::Parameter(theLin2d, gp_Pnt2d(x, ymax));
2851       nbi++;
2852     }
2853     else if (COINC(d1, d2, theTol)) {
2854       par[0] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmin, ymax));
2855       par[1] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmax, ymax));
2856       nbi = 2;
2857     }
2858   }
2859
2860   if(nbi == 2) {
2861
2862     if(fabs(par[0]-par[1]) > theTol) {
2863       theP1 = Min(par[0], par[1]);
2864       theP2 = Max(par[0], par[1]);
2865       return Standard_True;
2866     }
2867     else return Standard_False;
2868
2869   }
2870
2871   //xmax, ymax <-> xmax, ymin
2872   d1 = d2;
2873   d2 = A*xmax + B*ymin + C;
2874
2875   if(d1 > theTol || d1 < -theTol) {
2876     if(INTER(d1, d2, theTol)) {
2877       Standard_Real y = -(C + A*xmax)/B;
2878       par[nbi] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmax, y));
2879       nbi++;
2880     }
2881     else if (COINC(d1, d2, theTol)) {
2882       par[0] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmax, ymax));
2883       par[1] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmax, ymin));
2884       nbi = 2;
2885     }
2886   }
2887
2888   if(nbi == 2) {
2889     if(fabs(par[0]-par[1]) > theTol) {
2890       theP1 = Min(par[0], par[1]);
2891       theP2 = Max(par[0], par[1]);
2892       return Standard_True;
2893     }
2894     else return Standard_False;
2895   }
2896
2897   //xmax, ymin <-> xmin, ymin
2898   d1 = d2;
2899   d2 = A*xmin + B*ymin + C;
2900
2901   if(d1 > theTol || d1 < -theTol) {
2902     if(INTER(d1, d2, theTol)) {
2903       Standard_Real x = -(C + B*ymin)/A;
2904       par[nbi] = ElCLib::Parameter(theLin2d, gp_Pnt2d(x, ymin));
2905       nbi++;
2906     }
2907     else if (COINC(d1, d2, theTol)) {
2908       par[0] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmax, ymin));
2909       par[1] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmin, ymin));
2910       nbi = 2;
2911     }
2912   }
2913
2914   if(nbi == 2) {
2915     if(fabs(par[0]-par[1]) > theTol) {
2916       theP1 = Min(par[0], par[1]);
2917       theP2 = Max(par[0], par[1]);
2918       return Standard_True;
2919     }
2920     else return Standard_False;
2921   }
2922
2923   return Standard_False;
2924
2925 }
2926 //
2927 //=======================================================================
2928 //function : ApproxParameters
2929 //purpose  :
2930 //=======================================================================
2931 void ApproxParameters(const Handle(GeomAdaptor_HSurface)& aHS1,
2932                       const Handle(GeomAdaptor_HSurface)& aHS2,
2933                       Standard_Integer& iDegMin,
2934                       Standard_Integer& iDegMax,
2935                       Standard_Integer& iNbIter)
2936
2937 {
2938   GeomAbs_SurfaceType aTS1, aTS2;
2939
2940   //
2941   iNbIter=0;
2942   iDegMin=4;
2943   iDegMax=8;
2944   //
2945   aTS1=aHS1->Surface().GetType();
2946   aTS2=aHS2->Surface().GetType();
2947   //
2948   // Cylinder/Torus
2949   if ((aTS1==GeomAbs_Cylinder && aTS2==GeomAbs_Torus) ||
2950       (aTS2==GeomAbs_Cylinder && aTS1==GeomAbs_Torus)) {
2951     Standard_Real aRC, aRT, dR, aPC;
2952     gp_Cylinder aCylinder;
2953     gp_Torus aTorus;
2954     //
2955     aPC=Precision::Confusion();
2956     //
2957     aCylinder=(aTS1==GeomAbs_Cylinder)? aHS1->Surface().Cylinder() : aHS2->Surface().Cylinder();
2958     aTorus=(aTS1==GeomAbs_Torus)? aHS1->Surface().Torus() : aHS2->Surface().Torus();
2959     //
2960     aRC=aCylinder.Radius();
2961     aRT=aTorus.MinorRadius();
2962     dR=aRC-aRT;
2963     if (dR<0.) {
2964       dR=-dR;
2965     }
2966     //
2967     if (dR<aPC) {
2968       iDegMax=6;
2969     }
2970   }
2971   if (aTS1==GeomAbs_Cylinder && aTS2==GeomAbs_Cylinder) {
2972     iNbIter=1;
2973   }
2974 }
2975 //=======================================================================
2976 //function : Tolerances
2977 //purpose  :
2978 //=======================================================================
2979 void Tolerances(const Handle(GeomAdaptor_HSurface)& aHS1,
2980                 const Handle(GeomAdaptor_HSurface)& aHS2,
2981                 Standard_Real& aTolTang)
2982 {
2983   GeomAbs_SurfaceType aTS1, aTS2;
2984   //
2985   aTS1=aHS1->Surface().GetType();
2986   aTS2=aHS2->Surface().GetType();
2987   //
2988   // Cylinder/Torus
2989   if ((aTS1==GeomAbs_Cylinder && aTS2==GeomAbs_Torus) ||
2990       (aTS2==GeomAbs_Cylinder && aTS1==GeomAbs_Torus)) {
2991     Standard_Real aRC, aRT, dR, aPC;
2992     gp_Cylinder aCylinder;
2993     gp_Torus aTorus;
2994     //
2995     aPC=Precision::Confusion();
2996     //
2997     aCylinder=(aTS1==GeomAbs_Cylinder)? aHS1->Surface().Cylinder() : aHS2->Surface().Cylinder();
2998     aTorus=(aTS1==GeomAbs_Torus)? aHS1->Surface().Torus() : aHS2->Surface().Torus();
2999     //
3000     aRC=aCylinder.Radius();
3001     aRT=aTorus.MinorRadius();
3002     dR=aRC-aRT;
3003     if (dR<0.) {
3004       dR=-dR;
3005     }
3006     //
3007     if (dR<aPC) {
3008       aTolTang=0.1*aTolTang;
3009     }
3010   }
3011 }
3012 //=======================================================================
3013 //function : SortTypes
3014 //purpose  :
3015 //=======================================================================
3016 Standard_Boolean SortTypes(const GeomAbs_SurfaceType aType1,
3017                            const GeomAbs_SurfaceType aType2)
3018 {
3019   Standard_Boolean bRet;
3020   Standard_Integer aI1, aI2;
3021   //
3022   bRet=Standard_False;
3023   //
3024   aI1=IndexType(aType1);
3025   aI2=IndexType(aType2);
3026   if (aI1<aI2){
3027     bRet=!bRet;
3028   }
3029   return bRet;
3030 }
3031 //=======================================================================
3032 //function : IndexType
3033 //purpose  :
3034 //=======================================================================
3035 Standard_Integer IndexType(const GeomAbs_SurfaceType aType)
3036 {
3037   Standard_Integer aIndex;
3038   //
3039   aIndex=11;
3040   //
3041   if (aType==GeomAbs_Plane) {
3042     aIndex=0;
3043   }
3044   else if (aType==GeomAbs_Cylinder) {
3045     aIndex=1;
3046   }
3047   else if (aType==GeomAbs_Cone) {
3048     aIndex=2;
3049   }
3050   else if (aType==GeomAbs_Sphere) {
3051     aIndex=3;
3052   }
3053   else if (aType==GeomAbs_Torus) {
3054     aIndex=4;
3055   }
3056   else if (aType==GeomAbs_BezierSurface) {
3057     aIndex=5;
3058   }
3059   else if (aType==GeomAbs_BSplineSurface) {
3060     aIndex=6;
3061   }
3062   else if (aType==GeomAbs_SurfaceOfRevolution) {
3063     aIndex=7;
3064   }
3065   else if (aType==GeomAbs_SurfaceOfExtrusion) {
3066     aIndex=8;
3067   }
3068   else if (aType==GeomAbs_OffsetSurface) {
3069     aIndex=9;
3070   }
3071   else if (aType==GeomAbs_OtherSurface) {
3072     aIndex=10;
3073   }
3074   return aIndex;
3075 }
3076 #ifdef OCCT_DEBUG_DUMPWLINE
3077 //=======================================================================
3078 //function : DumpWLine
3079 //purpose  :
3080 //=======================================================================
3081 void DumpWLine(const Handle(IntPatch_WLine)& aWLine)
3082 {
3083   Standard_Integer i, aNbPnts;
3084   Standard_Real aX, aY, aZ, aU1, aV1, aU2, aV2;
3085   //
3086   printf(" *WLine\n");
3087   aNbPnts=aWLine->NbPnts();
3088   for (i=1; i<=aNbPnts; ++i) {
3089     const IntSurf_PntOn2S aPntOn2S=aWLine->Point(i);
3090     const gp_Pnt& aP3D=aPntOn2S.Value();
3091     aP3D.Coord(aX, aY, aZ);
3092     aPntOn2S.Parameters(aU1, aV1, aU2, aV2);
3093     //
3094     printf("point p_%d %lf %lf %lf\n", i, aX, aY, aZ);
3095     //printf("point p_%d %20.15lf %20.15lf %20.15lf %20.15lf %20.15lf %20.15lf %20.15lf\n",
3096         //   i, aX, aY, aZ, aU1, aV1, aU2, aV2);
3097   }
3098 }
3099 #endif
3100
3101 //=======================================================================
3102 // Function : FindMaxDistance
3103 // purpose :
3104 //=======================================================================
3105 Standard_Real FindMaxDistance(const Handle(Geom_Curve)& theCurve,
3106                               const Standard_Real theFirst,
3107                               const Standard_Real theLast,
3108                               const TopoDS_Face& theFace,
3109                               const Handle(IntTools_Context)& theContext)
3110 {
3111   Standard_Integer aNbS;
3112   Standard_Real aT1, aT2, aDt, aD, aDMax, anEps;
3113   //
3114   aNbS = 11;
3115   aDt = (theLast - theFirst) / aNbS;
3116   aDMax = 0.;
3117   anEps = 1.e-4 * aDt;
3118   //
3119   GeomAPI_ProjectPointOnSurf& aProjPS = theContext->ProjPS(theFace);
3120   aT2 = theFirst;
3121   for (;;) {
3122     aT1 = aT2;
3123     aT2 += aDt;
3124     //
3125     if (aT2 > theLast) {
3126       break;
3127     }
3128     //
3129     aD = FindMaxDistance(theCurve, aT1, aT2, aProjPS, anEps);
3130     if (aD > aDMax) {
3131       aDMax = aD;
3132     }
3133   }
3134   //
3135   return aDMax;
3136 }
3137
3138 //=======================================================================
3139 // Function : FindMaxDistance
3140 // purpose :
3141 //=======================================================================
3142 Standard_Real FindMaxDistance(const Handle(Geom_Curve)& theC,
3143                               const Standard_Real theFirst,
3144                               const Standard_Real theLast,
3145                               GeomAPI_ProjectPointOnSurf& theProjPS,
3146                               const Standard_Real theEps)
3147 {
3148   Standard_Real aA, aB, aCf, aX, aX1, aX2, aF1, aF2, aF;
3149   //
3150   aCf = 0.61803398874989484820458683436564;//(sqrt(5.)-1)/2.;
3151   aA = theFirst;
3152   aB = theLast;
3153   //
3154   aX1 = aB - aCf * (aB - aA);
3155   aF1 = MaxDistance(theC, aX1, theProjPS);
3156   aX2 = aA + aCf * (aB - aA);
3157   aF2 = MaxDistance(theC, aX2, theProjPS);
3158   //
3159   for (;;) {
3160     if ((aB - aA) < theEps) {
3161       break;
3162     }
3163     //
3164     if (aF1 > aF2) {
3165       aB = aX2;
3166       aX2 = aX1;
3167       aF2 = aF1;
3168       aX1 = aB - aCf * (aB - aA);
3169       aF1 = MaxDistance(theC, aX1, theProjPS);
3170     }
3171     else {
3172       aA = aX1;
3173       aX1 = aX2;
3174       aF1 = aF2;
3175       aX2 = aA + aCf * (aB - aA);
3176       aF2 = MaxDistance(theC, aX2, theProjPS);
3177     }
3178   }
3179   //
3180   aX = 0.5 * (aA + aB);
3181   aF = MaxDistance(theC, aX, theProjPS);
3182   //
3183   if (aF1 > aF) {
3184     aF = aF1;
3185   }
3186   //
3187   if (aF2 > aF) {
3188     aF = aF2;
3189   }
3190   //
3191   return aF;
3192 }
3193
3194 //=======================================================================
3195 // Function : MaxDistance
3196 // purpose :
3197 //=======================================================================
3198 Standard_Real MaxDistance(const Handle(Geom_Curve)& theC,
3199                           const Standard_Real aT,
3200                           GeomAPI_ProjectPointOnSurf& theProjPS)
3201 {
3202   Standard_Real aD;
3203   gp_Pnt aP;
3204   //
3205   theC->D0(aT, aP);
3206   theProjPS.Perform(aP);
3207   aD = theProjPS.NbPoints() ? theProjPS.LowerDistance() : 0.;
3208   //
3209   return aD;
3210 }
3211
3212 //=======================================================================
3213 //function : CheckPCurve
3214 //purpose  : Checks if points of the pcurve are out of the face bounds.
3215 //=======================================================================
3216   Standard_Boolean CheckPCurve(const Handle(Geom2d_Curve)& aPC,
3217                                const TopoDS_Face& aFace)
3218 {
3219   const Standard_Integer NPoints = 23;
3220   Standard_Integer i;
3221   Standard_Real umin,umax,vmin,vmax;
3222
3223   BRepTools::UVBounds(aFace, umin, umax, vmin, vmax);
3224   Standard_Real tolU = Max ((umax-umin)*0.01, Precision::Confusion());
3225   Standard_Real tolV = Max ((vmax-vmin)*0.01, Precision::Confusion());
3226   Standard_Real fp = aPC->FirstParameter();
3227   Standard_Real lp = aPC->LastParameter();
3228
3229
3230   // adjust domain for periodic surfaces
3231   TopLoc_Location aLoc;
3232   Handle(Geom_Surface) aSurf = BRep_Tool::Surface(aFace, aLoc);
3233   if (aSurf->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface))) {
3234     aSurf = (Handle(Geom_RectangularTrimmedSurface)::DownCast(aSurf))->BasisSurface();
3235   }
3236   gp_Pnt2d pnt = aPC->Value((fp+lp)/2);
3237   Standard_Real u,v;
3238   pnt.Coord(u,v);
3239   //
3240   if (aSurf->IsUPeriodic()) {
3241     Standard_Real aPer = aSurf->UPeriod();
3242     Standard_Integer nshift = (Standard_Integer) ((u-umin)/aPer);
3243     if (u < umin+aPer*nshift) nshift--;
3244     umin += aPer*nshift;
3245     umax += aPer*nshift;
3246   }
3247   if (aSurf->IsVPeriodic()) {
3248     Standard_Real aPer = aSurf->VPeriod();
3249     Standard_Integer nshift = (Standard_Integer) ((v-vmin)/aPer);
3250     if (v < vmin+aPer*nshift) nshift--;
3251     vmin += aPer*nshift;
3252     vmax += aPer*nshift;
3253   }
3254   //
3255   //--------------------------------------------------------
3256   Standard_Boolean bRet;
3257   Standard_Integer j, aNbIntervals;
3258   Standard_Real aT, dT;
3259   gp_Pnt2d aP2D;
3260   //
3261   Geom2dAdaptor_Curve aGAC(aPC);
3262   aNbIntervals=aGAC.NbIntervals(GeomAbs_CN);
3263   //
3264   TColStd_Array1OfReal aTI(1, aNbIntervals+1);
3265   aGAC.Intervals(aTI,GeomAbs_CN);
3266   //
3267   bRet=Standard_False;
3268   //
3269   aT=aGAC.FirstParameter();
3270   for (j=1; j<=aNbIntervals; ++j) {
3271     dT=(aTI(j+1)-aTI(j))/NPoints;
3272     //
3273     for (i=1; i<NPoints; i++) {
3274       aT=aT+dT;
3275       aGAC.D0(aT, aP2D);
3276       aP2D.Coord(u,v);
3277     if (umin-u > tolU || u-umax > tolU ||
3278           vmin-v > tolV || v-vmax > tolV) {
3279         return bRet;
3280   }
3281 }
3282   }
3283   return !bRet;
3284 }
3285 //=======================================================================
3286 //function : CorrectPlaneBoundaries
3287 //purpose  :
3288 //=======================================================================
3289  void CorrectPlaneBoundaries(Standard_Real& aUmin,
3290                              Standard_Real& aUmax,
3291                              Standard_Real& aVmin,
3292                              Standard_Real& aVmax)
3293 {
3294   if (!(Precision::IsInfinite(aUmin) ||
3295         Precision::IsInfinite(aUmax))) {
3296     Standard_Real dU;
3297     //
3298     dU=0.1*(aUmax-aUmin);
3299     aUmin=aUmin-dU;
3300     aUmax=aUmax+dU;
3301   }
3302   if (!(Precision::IsInfinite(aVmin) ||
3303         Precision::IsInfinite(aVmax))) {
3304     Standard_Real dV;
3305     //
3306     dV=0.1*(aVmax-aVmin);
3307     aVmin=aVmin-dV;
3308     aVmax=aVmax+dV;
3309   }
3310 }