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