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     //
 892     //if(!myListOfPnts.IsEmpty()) {
 893     //  bAvoidLineConstructor = Standard_True;
 894     //}
 895
 896     Standard_Integer nbp = aWLine->NbPnts();
 897     const IntSurf_PntOn2S& p1 = aWLine->Point(1);
 898     const IntSurf_PntOn2S& p2 = aWLine->Point(nbp);
 899
 900     const gp_Pnt& P1 = p1.Value();
 901     const gp_Pnt& P2 = p2.Value();
 902
 903     if(P1.SquareDistance(P2) < 1.e-14) {
 904       bAvoidLineConstructor = Standard_False;
 905     }
 906   }
 907
 908   typl=L->ArcType();
 909
 910   if(typl == IntPatch_Restriction)
 911     bAvoidLineConstructor = Standard_True;
 912
 913   //
 914   // Line Constructor
 915   if(!bAvoidLineConstructor) {
 916     myLConstruct.Perform(L);
 917     //
 918     bDone=myLConstruct.IsDone();
 919     if(!bDone)
 920     {
 921       return;
 922     }
 923
 924     if(typl != IntPatch_Restriction)
 925     {
 926       aNbParts=myLConstruct.NbParts();
 927       if (aNbParts <= 0)
 928       {
 929         return;
 930       }
 931     }
 932   }
 933   // Do the Curve
 934
 935
 936   switch (typl) {
 937   //########################################
 938   // Line, Parabola, Hyperbola
 939   //########################################
 940   case IntPatch_Lin:
 941   case IntPatch_Parabola:
 942   case IntPatch_Hyperbola: {
 943     if (typl == IntPatch_Lin) {
 944       newc =
 945         new Geom_Line (Handle(IntPatch_GLine)::DownCast(L)->Line());
 946     }
 947
 948     else if (typl == IntPatch_Parabola) {
 949       newc =
 950         new Geom_Parabola(Handle(IntPatch_GLine)::DownCast(L)->Parabola());
 951     }
 952
 953     else if (typl == IntPatch_Hyperbola) {
 954       newc =
 955         new Geom_Hyperbola (Handle(IntPatch_GLine)::DownCast(L)->Hyperbola());
 956     }
 957     //
 958     // myTolReached3d
 959     if (typl == IntPatch_Lin) {
 960       TolR3d (myFace1, myFace2, myTolReached3d);
 961     }
 962     //
 963     aNbParts=myLConstruct.NbParts();
 964     for (i=1; i<=aNbParts; i++) {
 965       Standard_Boolean bFNIt, bLPIt;
 966       //
 967       myLConstruct.Part(i, fprm, lprm);
 968         //
 969       bFNIt=Precision::IsNegativeInfinite(fprm);
 970       bLPIt=Precision::IsPositiveInfinite(lprm);
 971       //
 972       if (!bFNIt && !bLPIt) {
 973         //
 974         IntTools_Curve aCurve;
 975         //
 976         Handle(Geom_TrimmedCurve) aCT3D=new Geom_TrimmedCurve(newc, fprm, lprm);
 977         aCurve.SetCurve(aCT3D);
 978         if (typl == IntPatch_Parabola) {
 979           Standard_Real aTolF1, aTolF2, aTolBase;
 980
 981           aTolF1 = BRep_Tool::Tolerance(myFace1);
 982           aTolF2 = BRep_Tool::Tolerance(myFace2);
 983           aTolBase=aTolF1+aTolF2;
 984           myTolReached3d=IntTools_Tools::CurveTolerance(aCT3D, aTolBase);
 985         }
 986         //
 987         aCurve.SetCurve(new Geom_TrimmedCurve(newc, fprm, lprm));
 988         if(myApprox1) {
 989           Handle (Geom2d_Curve) C2d;
 990           GeomInt_IntSS::BuildPCurves(fprm, lprm, Tolpc,
 991                 myHS1->ChangeSurface().Surface(), newc, C2d);
 992           if(Tolpc>myTolReached2d || myTolReached2d==0.) {
 993             myTolReached2d=Tolpc;
 994           }
 995             //
 996             aCurve.SetFirstCurve2d(new Geom2d_TrimmedCurve(C2d,fprm,lprm));
 997           }
 998           else {
 999             Handle(Geom2d_BSplineCurve) H1;
1000             //
1001             aCurve.SetFirstCurve2d(H1);
1002           }
1003         //
1004         if(myApprox2) {
1005           Handle (Geom2d_Curve) C2d;
1006           GeomInt_IntSS::BuildPCurves(fprm, lprm, Tolpc,
1007                     myHS2->ChangeSurface().Surface(), newc, C2d);
1008           if(Tolpc>myTolReached2d || myTolReached2d==0.) {
1009             myTolReached2d=Tolpc;
1010           }
1011           //
1012           aCurve.SetSecondCurve2d(new Geom2d_TrimmedCurve(C2d,fprm,lprm));
1013           }
1014         else {
1015           Handle(Geom2d_BSplineCurve) H1;
1016           //
1017           aCurve.SetSecondCurve2d(H1);
1018         }
1019         mySeqOfCurve.Append(aCurve);
1020       } //if (!bFNIt && !bLPIt) {
1021       else {
1022         //  on regarde si on garde
1023         //
1024         Standard_Real aTestPrm, dT=100.;
1025         //
1026         aTestPrm=0.;
1027         if (bFNIt && !bLPIt) {
1028           aTestPrm=lprm-dT;
1029         }
1030         else if (!bFNIt && bLPIt) {
1031           aTestPrm=fprm+dT;
1032         }
1033         else {
1034           // i.e, if (bFNIt && bLPIt)
1035           aTestPrm=IntTools_Tools::IntermediatePoint(-dT, dT);
1036         }
1037         //
1038         gp_Pnt ptref(newc->Value(aTestPrm));
1039         //
1040         GeomAbs_SurfaceType typS1 = myHS1->GetType();
1041         GeomAbs_SurfaceType typS2 = myHS2->GetType();
1042         if( typS1 == GeomAbs_SurfaceOfExtrusion ||
1043             typS1 == GeomAbs_OffsetSurface ||
1044             typS1 == GeomAbs_SurfaceOfRevolution ||
1045             typS2 == GeomAbs_SurfaceOfExtrusion ||
1046             typS2 == GeomAbs_OffsetSurface ||
1047             typS2 == GeomAbs_SurfaceOfRevolution) {
1048           Handle(Geom2d_BSplineCurve) H1;
1049           mySeqOfCurve.Append(IntTools_Curve(newc, H1, H1));
1050           continue;
1051         }
1052
1053         Standard_Real u1, v1, u2, v2, Tol;
1054
1055         Tol = Precision::Confusion();
1056         Parameters(myHS1, myHS2, ptref,  u1, v1, u2, v2);
1057         ok = (dom1->Classify(gp_Pnt2d(u1, v1), Tol) != TopAbs_OUT);
1058         if(ok) {
1059           ok = (dom2->Classify(gp_Pnt2d(u2,v2),Tol) != TopAbs_OUT);
1060         }
1061         if (ok) {
1062           Handle(Geom2d_BSplineCurve) H1;
1063           mySeqOfCurve.Append(IntTools_Curve(newc, H1, H1));
1064         }
1065       }
1066     }// for (i=1; i<=aNbParts; i++) {
1067   }// case IntPatch_Lin:  case IntPatch_Parabola:  case IntPatch_Hyperbola:
1068     break;
1069
1070   //########################################
1071   // Circle and Ellipse
1072   //########################################
1073   case IntPatch_Circle:
1074   case IntPatch_Ellipse: {
1075
1076     if (typl == IntPatch_Circle) {
1077       newc = new Geom_Circle
1078         (Handle(IntPatch_GLine)::DownCast(L)->Circle());
1079     }
1080     else { //IntPatch_Ellipse
1081       newc = new Geom_Ellipse
1082         (Handle(IntPatch_GLine)::DownCast(L)->Ellipse());
1083     }
1084     //
1085     // myTolReached3d
1086     TolR3d (myFace1, myFace2, myTolReached3d);
1087     //
1088     aNbParts=myLConstruct.NbParts();
1089     //
1090     Standard_Real aPeriod, aNul;
1091     TColStd_SequenceOfReal aSeqFprm,  aSeqLprm;
1092
1093     aNul=0.;
1094     aPeriod=M_PI+M_PI;
1095
1096     for (i=1; i<=aNbParts; i++) {
1097       myLConstruct.Part(i, fprm, lprm);
1098
1099       if (fprm < aNul && lprm > aNul) {
1100         // interval that goes through 0. is divided on two intervals;
1101         while (fprm<aNul || fprm>aPeriod) fprm=fprm+aPeriod;
1102         while (lprm<aNul || lprm>aPeriod) lprm=lprm+aPeriod;
1103         //
1104         if((aPeriod - fprm) > Tolpc) {
1105           aSeqFprm.Append(fprm);
1106           aSeqLprm.Append(aPeriod);
1107         }
1108         else {
1109           gp_Pnt P1 = newc->Value(fprm);
1110           gp_Pnt P2 = newc->Value(aPeriod);
1111           Standard_Real aTolDist = BRep_Tool::Tolerance(myFace1) + BRep_Tool::Tolerance(myFace2);
1112           aTolDist = (myTolReached3d > aTolDist) ? myTolReached3d : aTolDist;
1113
1114           if(P1.Distance(P2) > aTolDist) {
1115             Standard_Real anewpar = fprm;
1116
1117             if(ParameterOutOfBoundary(fprm, newc, myFace1, myFace2,
1118                                       lprm, Standard_False, anewpar, myContext)) {
1119               fprm = anewpar;
1120             }
1121             aSeqFprm.Append(fprm);
1122             aSeqLprm.Append(aPeriod);
1123           }
1124         }
1125
1126         //
1127         if((lprm - aNul) > Tolpc) {
1128           aSeqFprm.Append(aNul);
1129           aSeqLprm.Append(lprm);
1130         }
1131         else {
1132           gp_Pnt P1 = newc->Value(aNul);
1133           gp_Pnt P2 = newc->Value(lprm);
1134           Standard_Real aTolDist = BRep_Tool::Tolerance(myFace1) + BRep_Tool::Tolerance(myFace2);
1135           aTolDist = (myTolReached3d > aTolDist) ? myTolReached3d : aTolDist;
1136
1137           if(P1.Distance(P2) > aTolDist) {
1138             Standard_Real anewpar = lprm;
1139
1140             if(ParameterOutOfBoundary(lprm, newc, myFace1, myFace2,
1141                                       fprm, Standard_True, anewpar, myContext)) {
1142               lprm = anewpar;
1143             }
1144             aSeqFprm.Append(aNul);
1145             aSeqLprm.Append(lprm);
1146           }
1147         }
1148       }
1149       else {
1150         // usual interval
1151         aSeqFprm.Append(fprm);
1152         aSeqLprm.Append(lprm);
1153       }
1154     }
1155
1156     //
1157     aNbParts=aSeqFprm.Length();
1158     for (i=1; i<=aNbParts; i++) {
1159       fprm=aSeqFprm(i);
1160       lprm=aSeqLprm(i);
1161       //
1162       Standard_Real aRealEpsilon=RealEpsilon();
1163       if (Abs(fprm) > aRealEpsilon || Abs(lprm-2.*M_PI) > aRealEpsilon) {
1164         //==============================================
1165         ////
1166         IntTools_Curve aCurve;
1167         Handle(Geom_TrimmedCurve) aTC3D=new Geom_TrimmedCurve(newc,fprm,lprm);
1168         aCurve.SetCurve(aTC3D);
1169         fprm=aTC3D->FirstParameter();
1170         lprm=aTC3D->LastParameter ();
1171         ////
1172         if (typl == IntPatch_Circle || typl == IntPatch_Ellipse) {////
1173           if(myApprox1) {
1174             Handle (Geom2d_Curve) C2d;
1175             GeomInt_IntSS::BuildPCurves(fprm, lprm, Tolpc,
1176                         myHS1->ChangeSurface().Surface(), newc, C2d);
1177             if(Tolpc>myTolReached2d || myTolReached2d==0) {
1178               myTolReached2d=Tolpc;
1179             }
1180             //
1181             aCurve.SetFirstCurve2d(C2d);
1182           }
1183           else { ////
1184             Handle(Geom2d_BSplineCurve) H1;
1185             aCurve.SetFirstCurve2d(H1);
1186           }
1187
1188
1189           if(myApprox2) {
1190             Handle (Geom2d_Curve) C2d;
1191             GeomInt_IntSS::BuildPCurves(fprm,lprm,Tolpc,
1192                     myHS2->ChangeSurface().Surface(),newc,C2d);
1193             if(Tolpc>myTolReached2d || myTolReached2d==0) {
1194               myTolReached2d=Tolpc;
1195             }
1196             //
1197             aCurve.SetSecondCurve2d(C2d);
1198           }
1199           else {
1200             Handle(Geom2d_BSplineCurve) H1;
1201             aCurve.SetSecondCurve2d(H1);
1202           }
1203         }
1204
1205         else {
1206           Handle(Geom2d_BSplineCurve) H1;
1207           aCurve.SetFirstCurve2d(H1);
1208           aCurve.SetSecondCurve2d(H1);
1209         }
1210         mySeqOfCurve.Append(aCurve);
1211           //==============================================
1212       } //if (Abs(fprm) > RealEpsilon() || Abs(lprm-2.*M_PI) > RealEpsilon())
1213
1214       else {
1215         //  on regarde si on garde
1216         //
1217         if (aNbParts==1) {
1218 //           if (Abs(fprm) < RealEpsilon() &&  Abs(lprm-2.*M_PI) < RealEpsilon()) {
1219           if (Abs(fprm) <= aRealEpsilon && Abs(lprm-2.*M_PI) <= aRealEpsilon) {
1220             IntTools_Curve aCurve;
1221             Handle(Geom_TrimmedCurve) aTC3D=new Geom_TrimmedCurve(newc,fprm,lprm);
1222             aCurve.SetCurve(aTC3D);
1223             fprm=aTC3D->FirstParameter();
1224             lprm=aTC3D->LastParameter ();
1225
1226             if(myApprox1) {
1227               Handle (Geom2d_Curve) C2d;
1228               GeomInt_IntSS::BuildPCurves(fprm,lprm,Tolpc,
1229                     myHS1->ChangeSurface().Surface(),newc,C2d);
1230               if(Tolpc>myTolReached2d || myTolReached2d==0) {
1231                 myTolReached2d=Tolpc;
1232               }
1233               //
1234               aCurve.SetFirstCurve2d(C2d);
1235             }
1236             else { ////
1237               Handle(Geom2d_BSplineCurve) H1;
1238               aCurve.SetFirstCurve2d(H1);
1239             }
1240
1241             if(myApprox2) {
1242               Handle (Geom2d_Curve) C2d;
1243               GeomInt_IntSS::BuildPCurves(fprm,lprm,Tolpc,
1244                     myHS2->ChangeSurface().Surface(),newc,C2d);
1245               if(Tolpc>myTolReached2d || myTolReached2d==0) {
1246                 myTolReached2d=Tolpc;
1247               }
1248               //
1249               aCurve.SetSecondCurve2d(C2d);
1250             }
1251             else {
1252               Handle(Geom2d_BSplineCurve) H1;
1253               aCurve.SetSecondCurve2d(H1);
1254             }
1255             mySeqOfCurve.Append(aCurve);
1256             break;
1257           }
1258         }
1259         //
1260         Standard_Real aTwoPIdiv17, u1, v1, u2, v2, Tol;
1261
1262         aTwoPIdiv17=2.*M_PI/17.;
1263
1264         for (j=0; j<=17; j++) {
1265           gp_Pnt ptref (newc->Value (j*aTwoPIdiv17));
1266           Tol = Precision::Confusion();
1267
1268           Parameters(myHS1, myHS2, ptref, u1, v1, u2, v2);
1269           ok = (dom1->Classify(gp_Pnt2d(u1,v1),Tol) != TopAbs_OUT);
1270           if(ok) {
1271             ok = (dom2->Classify(gp_Pnt2d(u2,v2),Tol) != TopAbs_OUT);
1272           }
1273           if (ok) {
1274             IntTools_Curve aCurve;
1275             aCurve.SetCurve(newc);
1276             //==============================================
1277             if (typl == IntPatch_Circle || typl == IntPatch_Ellipse) {
1278
1279               if(myApprox1) {
1280                 Handle (Geom2d_Curve) C2d;
1281                 GeomInt_IntSS::BuildPCurves(fprm, lprm, Tolpc,
1282                         myHS1->ChangeSurface().Surface(), newc, C2d);
1283                 if(Tolpc>myTolReached2d || myTolReached2d==0) {
1284                   myTolReached2d=Tolpc;
1285                 }
1286                 //
1287                 aCurve.SetFirstCurve2d(C2d);
1288               }
1289               else {
1290                 Handle(Geom2d_BSplineCurve) H1;
1291                 aCurve.SetFirstCurve2d(H1);
1292               }
1293
1294               if(myApprox2) {
1295                 Handle (Geom2d_Curve) C2d;
1296                 GeomInt_IntSS::BuildPCurves(fprm, lprm, Tolpc,
1297                         myHS2->ChangeSurface().Surface(), newc, C2d);
1298                 if(Tolpc>myTolReached2d || myTolReached2d==0) {
1299                   myTolReached2d=Tolpc;
1300                 }
1301                 //
1302                 aCurve.SetSecondCurve2d(C2d);
1303               }
1304
1305               else {
1306                 Handle(Geom2d_BSplineCurve) H1;
1307                 aCurve.SetSecondCurve2d(H1);
1308               }
1309             }//  end of if (typl == IntPatch_Circle || typl == IntPatch_Ellipse)
1310
1311             else {
1312               Handle(Geom2d_BSplineCurve) H1;
1313               //
1314               aCurve.SetFirstCurve2d(H1);
1315               aCurve.SetSecondCurve2d(H1);
1316             }
1317             //==============================================
1318             //
1319             mySeqOfCurve.Append(aCurve);
1320             break;
1321
1322             }//  end of if (ok) {
1323           }//  end of for (Standard_Integer j=0; j<=17; j++)
1324         }//  end of else { on regarde si on garde
1325       }// for (i=1; i<=myLConstruct.NbParts(); i++)
1326     }// IntPatch_Circle: IntPatch_Ellipse:
1327     break;
1328
1329   case IntPatch_Analytic: {
1330     IntSurf_Quadric quad1,quad2;
1331     GeomAbs_SurfaceType typs = myHS1->Surface().GetType();
1332
1333     switch (typs) {
1334       case GeomAbs_Plane:
1335         quad1.SetValue(myHS1->Surface().Plane());
1336         break;
1337       case GeomAbs_Cylinder:
1338         quad1.SetValue(myHS1->Surface().Cylinder());
1339         break;
1340       case GeomAbs_Cone:
1341         quad1.SetValue(myHS1->Surface().Cone());
1342         break;
1343       case GeomAbs_Sphere:
1344         quad1.SetValue(myHS1->Surface().Sphere());
1345         break;
1346     case GeomAbs_Torus:
1347       quad1.SetValue(myHS1->Surface().Torus());
1348       break;
1349       default:
1350         Standard_ConstructionError::Raise("GeomInt_IntSS::MakeCurve 1");
1351       }
1352
1353     typs = myHS2->Surface().GetType();
1354
1355     switch (typs) {
1356       case GeomAbs_Plane:
1357         quad2.SetValue(myHS2->Surface().Plane());
1358         break;
1359       case GeomAbs_Cylinder:
1360         quad2.SetValue(myHS2->Surface().Cylinder());
1361         break;
1362       case GeomAbs_Cone:
1363         quad2.SetValue(myHS2->Surface().Cone());
1364         break;
1365       case GeomAbs_Sphere:
1366         quad2.SetValue(myHS2->Surface().Sphere());
1367         break;
1368     case GeomAbs_Torus:
1369       quad2.SetValue(myHS2->Surface().Torus());
1370       break;
1371       default:
1372         Standard_ConstructionError::Raise("GeomInt_IntSS::MakeCurve 2");
1373       }
1374     //
1375     //=========
1376     IntPatch_ALineToWLine convert (quad1, quad2);
1377
1378     if (!myApprox) {
1379       aNbParts=myLConstruct.NbParts();
1380       for (i=1; i<=aNbParts; i++) {
1381         myLConstruct.Part(i, fprm, lprm);
1382         Handle(IntPatch_WLine) WL =
1383           convert.MakeWLine(Handle(IntPatch_ALine)::DownCast(L), fprm, lprm);
1384         //
1385         Handle(Geom2d_BSplineCurve) H1;
1386         Handle(Geom2d_BSplineCurve) H2;
1387
1388         if(myApprox1) {
1389           H1 = GeomInt_IntSS::MakeBSpline2d(WL, 1, WL->NbPnts(), Standard_True);
1390         }
1391
1392         if(myApprox2) {
1393           H2 = GeomInt_IntSS::MakeBSpline2d(WL, 1, WL->NbPnts(), Standard_False);
1394         }
1395         //
1396         mySeqOfCurve.Append(IntTools_Curve(GeomInt_IntSS::MakeBSpline(WL,1,WL->NbPnts()), H1, H2));
1397       }
1398     } // if (!myApprox)
1399
1400     else { // myApprox=TRUE
1401       GeomInt_WLApprox theapp3d;
1402       //
1403       Standard_Real tol2d = myTolApprox;
1404       //
1405       theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, 30, Standard_True);
1406
1407       aNbParts=myLConstruct.NbParts();
1408       for (i=1; i<=aNbParts; i++) {
1409         myLConstruct.Part(i, fprm, lprm);
1410         Handle(IntPatch_WLine) WL =
1411           convert.MakeWLine(Handle(IntPatch_ALine):: DownCast(L),fprm,lprm);
1412
1413         theapp3d.Perform(myHS1,myHS2,WL,Standard_True,myApprox1,myApprox2, 1, WL->NbPnts());
1414
1415         if (!theapp3d.IsDone()) {
1416           //
1417           Handle(Geom2d_BSplineCurve) H1;
1418           Handle(Geom2d_BSplineCurve) H2;
1419
1420           if(myApprox1) {
1421             H1 = GeomInt_IntSS::MakeBSpline2d(WL, 1, WL->NbPnts(), Standard_True);
1422           }
1423
1424           if(myApprox2) {
1425             H2 = GeomInt_IntSS::MakeBSpline2d(WL, 1, WL->NbPnts(), Standard_False);
1426           }
1427           //
1428           mySeqOfCurve.Append(IntTools_Curve(GeomInt_IntSS::MakeBSpline(WL,1,WL->NbPnts()), H1, H2));
1429         }
1430
1431         else {
1432           if(myApprox1 || myApprox2) {
1433             if( theapp3d.TolReached2d()>myTolReached2d || myTolReached2d==0) {
1434               myTolReached2d = theapp3d.TolReached2d();
1435             }
1436           }
1437
1438           if( theapp3d.TolReached3d()>myTolReached3d || myTolReached3d==0) {
1439             myTolReached3d = theapp3d.TolReached3d();
1440           }
1441
1442           Standard_Integer aNbMultiCurves, nbpoles;
1443           aNbMultiCurves=theapp3d.NbMultiCurves();
1444           for (j=1; j<=aNbMultiCurves; j++) {
1445             const AppParCurves_MultiBSpCurve& mbspc = theapp3d.Value(j);
1446             nbpoles = mbspc.NbPoles();
1447
1448             TColgp_Array1OfPnt tpoles(1, nbpoles);
1449             mbspc.Curve(1, tpoles);
1450             Handle(Geom_BSplineCurve) BS=new Geom_BSplineCurve(tpoles,
1451                                                                mbspc.Knots(),
1452                                                                mbspc.Multiplicities(),
1453                                                                mbspc.Degree());
1454
1455             GeomLib_CheckBSplineCurve Check(BS,TOLCHECK,TOLANGCHECK);
1456             Check.FixTangent(Standard_True,Standard_True);
1457             //
1458             IntTools_Curve aCurve;
1459             aCurve.SetCurve(BS);
1460
1461             if(myApprox1) {
1462               TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
1463               mbspc.Curve(2,tpoles2d);
1464               Handle(Geom2d_BSplineCurve) BS2=new Geom2d_BSplineCurve(tpoles2d,
1465                                                                       mbspc.Knots(),
1466                                                                       mbspc.Multiplicities(),
1467                                                                       mbspc.Degree());
1468
1469               GeomLib_Check2dBSplineCurve newCheck(BS2,TOLCHECK,TOLANGCHECK);
1470               newCheck.FixTangent(Standard_True,Standard_True);
1471               //
1472               aCurve.SetFirstCurve2d(BS2);
1473             }
1474             else {
1475               Handle(Geom2d_BSplineCurve) H1;
1476               aCurve.SetFirstCurve2d(H1);
1477             }
1478
1479             if(myApprox2) {
1480               TColgp_Array1OfPnt2d tpoles2d(1, nbpoles);
1481               Standard_Integer TwoOrThree;
1482               TwoOrThree=myApprox1 ? 3 : 2;
1483               mbspc.Curve(TwoOrThree, tpoles2d);
1484               Handle(Geom2d_BSplineCurve) BS2 =new Geom2d_BSplineCurve(tpoles2d,
1485                                                                        mbspc.Knots(),
1486                                                                        mbspc.Multiplicities(),
1487                                                                        mbspc.Degree());
1488
1489               GeomLib_Check2dBSplineCurve newCheck(BS2,TOLCHECK,TOLANGCHECK);
1490               newCheck.FixTangent(Standard_True,Standard_True);
1491               //
1492               aCurve.SetSecondCurve2d(BS2);
1493             }
1494             else {
1495               Handle(Geom2d_BSplineCurve) H2;
1496               aCurve.SetSecondCurve2d(H2);
1497             }
1498             //
1499             mySeqOfCurve.Append(aCurve);
1500
1501           }// for (j=1; j<=aNbMultiCurves; j++) {
1502         }// else from if (!theapp3d.IsDone())
1503       }// for (i=1; i<=aNbParts; i++) {
1504     }// else { // myApprox=TRUE
1505   }// case IntPatch_Analytic:
1506     break;
1507
1508   case IntPatch_Walking:{
1509     Handle(IntPatch_WLine) WL =
1510       Handle(IntPatch_WLine)::DownCast(L);
1511
1512 #ifdef INTTOOLS_FACEFACE_DEBUG
1513     WL->Dump(0);
1514 #endif
1515
1516     //
1517     Standard_Integer ifprm, ilprm;
1518     //
1519     if (!myApprox) {
1520       aNbParts = 1;
1521       if(!bAvoidLineConstructor){
1522         aNbParts=myLConstruct.NbParts();
1523       }
1524       for (i=1; i<=aNbParts; ++i) {
1525         Handle(Geom2d_BSplineCurve) H1, H2;
1526         Handle(Geom_Curve) aBSp;
1527         //
1528         if(bAvoidLineConstructor) {
1529           ifprm = 1;
1530           ilprm = WL->NbPnts();
1531         }
1532         else {
1533           myLConstruct.Part(i, fprm, lprm);
1534           ifprm=(Standard_Integer)fprm;
1535           ilprm=(Standard_Integer)lprm;
1536         }
1537         //
1538         if(myApprox1) {
1539           H1 = GeomInt_IntSS::MakeBSpline2d(WL, ifprm, ilprm, Standard_True);
1540         }
1541         //
1542         if(myApprox2) {
1543           H2 = GeomInt_IntSS::MakeBSpline2d(WL, ifprm, ilprm, Standard_False);
1544         }
1545         //
1546         aBSp=GeomInt_IntSS::MakeBSpline(WL, ifprm, ilprm);
1547         IntTools_Curve aIC(aBSp, H1, H2);
1548         mySeqOfCurve.Append(aIC);
1549       }// for (i=1; i<=aNbParts; ++i) {
1550     }// if (!myApprox) {
1551     //
1552     else { // X
1553       Standard_Boolean bIsDecomposited;
1554       Standard_Integer nbiter, aNbSeqOfL;
1555       Standard_Real tol2d, aTolApproxImp;
1556       IntPatch_SequenceOfLine aSeqOfL;
1557       GeomInt_WLApprox theapp3d;
1558       Approx_ParametrizationType aParType = Approx_ChordLength;
1559       //
1560       Standard_Boolean anApprox1 = myApprox1;
1561       Standard_Boolean anApprox2 = myApprox2;
1562       //
1563       aTolApproxImp=1.e-5;
1564       tol2d = myTolApprox;
1565
1566       GeomAbs_SurfaceType typs1, typs2;
1567       typs1 = myHS1->Surface().GetType();
1568       typs2 = myHS2->Surface().GetType();
1569       Standard_Boolean anWithPC = Standard_True;
1570
1571       if(typs1 == GeomAbs_Cylinder && typs2 == GeomAbs_Sphere) {
1572         anWithPC =
1573           ApproxWithPCurves(myHS1->Surface().Cylinder(), myHS2->Surface().Sphere());
1574       }
1575       else if (typs1 == GeomAbs_Sphere && typs2 == GeomAbs_Cylinder) {
1576         anWithPC =
1577           ApproxWithPCurves(myHS2->Surface().Cylinder(), myHS1->Surface().Sphere());
1578       }
1579       //
1580       if(!anWithPC) {
1581         myTolApprox = aTolApproxImp;//1.e-5;
1582         anApprox1 = Standard_False;
1583         anApprox2 = Standard_False;
1584         //
1585         tol2d = myTolApprox;
1586       }
1587
1588       if(myHS1 == myHS2) {
1589         theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, 30, Standard_False, aParType);
1590         rejectSurface = Standard_True;
1591       }
1592       else {
1593         if(reApprox && !rejectSurface)
1594           theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, 30, Standard_False, aParType);
1595         else {
1596           Standard_Integer iDegMax, iDegMin, iNbIter;
1597           //
1598           ApproxParameters(myHS1, myHS2, iDegMin, iDegMax, iNbIter);
1599           theapp3d.SetParameters(myTolApprox, tol2d, iDegMin, iDegMax,
1600                                                   iNbIter, 30, Standard_True, aParType);
1601         }
1602       }
1603       //
1604       Standard_Real aReachedTol = Precision::Confusion();
1605       bIsDecomposited = IntTools_WLineTool::
1606         DecompositionOfWLine(WL,
1607                              myHS1,
1608                              myHS2,
1609                              myFace1,
1610                              myFace2,
1611                              myLConstruct,
1612                              bAvoidLineConstructor,
1613                              aSeqOfL,
1614                              aReachedTol,
1615                              myContext);
1616       if ( bIsDecomposited && ( myTolReached3d < aReachedTol ) ) {
1617         myTolReached3d = aReachedTol;
1618       }
1619       //
1620       aNbSeqOfL=aSeqOfL.Length();
1621       //
1622       if (bIsDecomposited) {
1623         nbiter=aNbSeqOfL;
1624       }
1625       else {
1626         nbiter=1;
1627         aNbParts=1;
1628         if (!bAvoidLineConstructor) {
1629           aNbParts=myLConstruct.NbParts();
1630           nbiter=aNbParts;
1631         }
1632       }
1633       //
1634       for(i = 1; i <= nbiter; ++i) {
1635         if(bIsDecomposited) {
1636           WL = Handle(IntPatch_WLine)::DownCast(aSeqOfL.Value(i));
1637           ifprm = 1;
1638           ilprm = WL->NbPnts();
1639         }
1640         else {
1641           if(bAvoidLineConstructor) {
1642             ifprm = 1;
1643             ilprm = WL->NbPnts();
1644           }
1645           else {
1646             myLConstruct.Part(i, fprm, lprm);
1647             ifprm = (Standard_Integer)fprm;
1648             ilprm = (Standard_Integer)lprm;
1649           }
1650         }
1651         //-- lbr :
1652         //-- Si une des surfaces est un plan , on approxime en 2d
1653         //-- sur cette surface et on remonte les points 2d en 3d.
1654         if(typs1 == GeomAbs_Plane) {
1655           theapp3d.Perform(myHS1, myHS2, WL, Standard_False,Standard_True, myApprox2,ifprm,ilprm);
1656         }
1657         else if(typs2 == GeomAbs_Plane) {
1658           theapp3d.Perform(myHS1,myHS2,WL,Standard_False,myApprox1,Standard_True,ifprm,ilprm);
1659         }
1660         else {
1661           //
1662           if (myHS1 != myHS2){
1663             if ((typs1==GeomAbs_BezierSurface || typs1==GeomAbs_BSplineSurface) &&
1664                 (typs2==GeomAbs_BezierSurface || typs2==GeomAbs_BSplineSurface)) {
1665
1666               theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, 30,
1667                                                                 Standard_True, aParType);
1668
1669               Standard_Boolean bUseSurfaces;
1670               bUseSurfaces = IntTools_WLineTool::NotUseSurfacesForApprox(myFace1, myFace2, WL, ifprm,  ilprm);
1671               if (bUseSurfaces) {
1672                 // ######
1673                 rejectSurface = Standard_True;
1674                 // ######
1675                 theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, 30,
1676                                                                 Standard_False, aParType);
1677               }
1678             }
1679           }
1680           //
1681           theapp3d.Perform(myHS1,myHS2,WL,Standard_True,anApprox1,anApprox2,ifprm,ilprm);
1682         }
1683           //
1684         if (!theapp3d.IsDone()) {
1685           Handle(Geom2d_BSplineCurve) H1;
1686           Handle(Geom2d_BSplineCurve) H2;
1687           //
1688           Handle(Geom_Curve) aBSp=GeomInt_IntSS::MakeBSpline(WL,ifprm, ilprm);
1689           //
1690           if(myApprox1) {
1691             H1 = GeomInt_IntSS::MakeBSpline2d(WL, ifprm, ilprm, Standard_True);
1692           }
1693           //
1694           if(myApprox2) {
1695             H2 = GeomInt_IntSS::MakeBSpline2d(WL, ifprm, ilprm, Standard_False);
1696           }
1697           //
1698           IntTools_Curve aIC(aBSp, H1, H2);
1699           mySeqOfCurve.Append(aIC);
1700         }
1701
1702         else {
1703           if(myApprox1 || myApprox2 || (typs1==GeomAbs_Plane || typs2==GeomAbs_Plane)) {
1704             if( theapp3d.TolReached2d()>myTolReached2d || myTolReached2d==0.) {
1705               myTolReached2d = theapp3d.TolReached2d();
1706             }
1707           }
1708           if(typs1==GeomAbs_Plane || typs2==GeomAbs_Plane) {
1709             myTolReached3d = myTolReached2d;
1710             //
1711             if (typs1==GeomAbs_Torus || typs2==GeomAbs_Torus) {
1712               if (myTolReached3d<1.e-6) {
1713                 myTolReached3d = theapp3d.TolReached3d();
1714                 myTolReached3d=1.e-6;
1715               }
1716             }
1717           }
1718           else  if( theapp3d.TolReached3d()>myTolReached3d || myTolReached3d==0.) {
1719             myTolReached3d = theapp3d.TolReached3d();
1720           }
1721
1722           Standard_Integer aNbMultiCurves, nbpoles;
1723           aNbMultiCurves=theapp3d.NbMultiCurves();
1724           for (j=1; j<=aNbMultiCurves; j++) {
1725             if(typs1 == GeomAbs_Plane) {
1726               const AppParCurves_MultiBSpCurve& mbspc = theapp3d.Value(j);
1727               nbpoles = mbspc.NbPoles();
1728
1729               TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
1730               TColgp_Array1OfPnt   tpoles(1,nbpoles);
1731
1732               mbspc.Curve(1,tpoles2d);
1733               const gp_Pln&  Pln = myHS1->Surface().Plane();
1734               //
1735               Standard_Integer ik;
1736               for(ik = 1; ik<= nbpoles; ik++) {
1737                 tpoles.SetValue(ik,
1738                                 ElSLib::Value(tpoles2d.Value(ik).X(),
1739                                               tpoles2d.Value(ik).Y(),
1740                                               Pln));
1741               }
1742               //
1743               Handle(Geom_BSplineCurve) BS =
1744                 new Geom_BSplineCurve(tpoles,
1745                                       mbspc.Knots(),
1746                                       mbspc.Multiplicities(),
1747                                       mbspc.Degree());
1748               GeomLib_CheckBSplineCurve Check(BS,TOLCHECK,TOLANGCHECK);
1749               Check.FixTangent(Standard_True, Standard_True);
1750               //
1751               IntTools_Curve aCurve;
1752               aCurve.SetCurve(BS);
1753
1754               if(myApprox1) {
1755                 Handle(Geom2d_BSplineCurve) BS1 =
1756                   new Geom2d_BSplineCurve(tpoles2d,
1757                                           mbspc.Knots(),
1758                                           mbspc.Multiplicities(),
1759                                           mbspc.Degree());
1760                 GeomLib_Check2dBSplineCurve Check1(BS1,TOLCHECK,TOLANGCHECK);
1761                 Check1.FixTangent(Standard_True,Standard_True);
1762                 //
1763                 // ############################################
1764                 if(!rejectSurface && !reApprox) {
1765                   Standard_Boolean isValid = IsCurveValid(BS1);
1766                   if(!isValid) {
1767                     reApprox = Standard_True;
1768                     goto reapprox;
1769                   }
1770                 }
1771                 // ############################################
1772                 aCurve.SetFirstCurve2d(BS1);
1773               }
1774               else {
1775                 Handle(Geom2d_BSplineCurve) H1;
1776                 aCurve.SetFirstCurve2d(H1);
1777               }
1778
1779               if(myApprox2) {
1780                 mbspc.Curve(2, tpoles2d);
1781
1782                 Handle(Geom2d_BSplineCurve) BS2 = new Geom2d_BSplineCurve(tpoles2d,
1783                                                                           mbspc.Knots(),
1784                                                                           mbspc.Multiplicities(),
1785                                                                           mbspc.Degree());
1786                 GeomLib_Check2dBSplineCurve newCheck(BS2,TOLCHECK,TOLANGCHECK);
1787                 newCheck.FixTangent(Standard_True,Standard_True);
1788
1789                 // ###########################################
1790                 if(!rejectSurface && !reApprox) {
1791                   Standard_Boolean isValid = IsCurveValid(BS2);
1792                   if(!isValid) {
1793                     reApprox = Standard_True;
1794                     goto reapprox;
1795                   }
1796                 }
1797                 // ###########################################
1798                 //
1799                 aCurve.SetSecondCurve2d(BS2);
1800               }
1801               else {
1802                 Handle(Geom2d_BSplineCurve) H2;
1803                 //
1804                   aCurve.SetSecondCurve2d(H2);
1805               }
1806               //
1807               mySeqOfCurve.Append(aCurve);
1808
1809             }//if(typs1 == GeomAbs_Plane) {
1810
1811             else if(typs2 == GeomAbs_Plane)
1812             {
1813               const AppParCurves_MultiBSpCurve& mbspc = theapp3d.Value(j);
1814               nbpoles = mbspc.NbPoles();
1815
1816               TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
1817               TColgp_Array1OfPnt   tpoles(1,nbpoles);
1818               mbspc.Curve((myApprox1==Standard_True)? 2 : 1,tpoles2d);
1819               const gp_Pln&  Pln = myHS2->Surface().Plane();
1820               //
1821               Standard_Integer ik;
1822               for(ik = 1; ik<= nbpoles; ik++) {
1823                 tpoles.SetValue(ik,
1824                                 ElSLib::Value(tpoles2d.Value(ik).X(),
1825                                               tpoles2d.Value(ik).Y(),
1826                                               Pln));
1827
1828               }
1829               //
1830               Handle(Geom_BSplineCurve) BS=new Geom_BSplineCurve(tpoles,
1831                                                                  mbspc.Knots(),
1832                                                                  mbspc.Multiplicities(),
1833                                                                  mbspc.Degree());
1834               GeomLib_CheckBSplineCurve Check(BS,TOLCHECK,TOLANGCHECK);
1835               Check.FixTangent(Standard_True,Standard_True);
1836               //
1837               IntTools_Curve aCurve;
1838               aCurve.SetCurve(BS);
1839
1840               if(myApprox2) {
1841                 Handle(Geom2d_BSplineCurve) BS1=new Geom2d_BSplineCurve(tpoles2d,
1842                                                                         mbspc.Knots(),
1843                                                                         mbspc.Multiplicities(),
1844                                                                         mbspc.Degree());
1845                 GeomLib_Check2dBSplineCurve Check1(BS1,TOLCHECK,TOLANGCHECK);
1846                 Check1.FixTangent(Standard_True,Standard_True);
1847                 //
1848                 // ###########################################
1849                 if(!rejectSurface && !reApprox) {
1850                   Standard_Boolean isValid = IsCurveValid(BS1);
1851                   if(!isValid) {
1852                     reApprox = Standard_True;
1853                     goto reapprox;
1854                   }
1855                 }
1856                 // ###########################################
1857                 bPCurvesOk = CheckPCurve(BS1, myFace2);
1858                 aCurve.SetSecondCurve2d(BS1);
1859               }
1860               else {
1861                 Handle(Geom2d_BSplineCurve) H2;
1862                 aCurve.SetSecondCurve2d(H2);
1863               }
1864
1865               if(myApprox1) {
1866                 mbspc.Curve(1,tpoles2d);
1867                 Handle(Geom2d_BSplineCurve) BS2=new Geom2d_BSplineCurve(tpoles2d,
1868                                                                         mbspc.Knots(),
1869                                                                         mbspc.Multiplicities(),
1870                                                                         mbspc.Degree());
1871                 GeomLib_Check2dBSplineCurve Check2(BS2,TOLCHECK,TOLANGCHECK);
1872                 Check2.FixTangent(Standard_True,Standard_True);
1873                 //
1874                 // ###########################################
1875                 if(!rejectSurface && !reApprox) {
1876                   Standard_Boolean isValid = IsCurveValid(BS2);
1877                   if(!isValid) {
1878                     reApprox = Standard_True;
1879                     goto reapprox;
1880                   }
1881                 }
1882                 // ###########################################
1883                 bPCurvesOk = bPCurvesOk && CheckPCurve(BS2, myFace1);
1884                 aCurve.SetFirstCurve2d(BS2);
1885               }
1886               else {
1887                 Handle(Geom2d_BSplineCurve) H1;
1888                 //
1889                 aCurve.SetFirstCurve2d(H1);
1890               }
1891               //
1892               //if points of the pcurves are out of the faces bounds
1893               //create 3d and 2d curves without approximation
1894               if (!bPCurvesOk) {
1895                 Handle(Geom2d_BSplineCurve) H1, H2;
1896                 bPCurvesOk = Standard_True;
1897                 //
1898                 Handle(Geom_Curve) aBSp=GeomInt_IntSS::MakeBSpline(WL,ifprm, ilprm);
1899
1900                 if(myApprox1) {
1901                   H1 = GeomInt_IntSS::MakeBSpline2d(WL, ifprm, ilprm, Standard_True);
1902                   bPCurvesOk = CheckPCurve(H1, myFace1);
1903                 }
1904
1905                 if(myApprox2) {
1906                   H2 = GeomInt_IntSS::MakeBSpline2d(WL, ifprm, ilprm, Standard_False);
1907                   bPCurvesOk = bPCurvesOk && CheckPCurve(H2, myFace2);
1908                 }
1909                 //
1910                 //if pcurves created without approximation are out of the
1911                 //faces bounds, use approximated 3d and 2d curves
1912                 if (bPCurvesOk) {
1913                   IntTools_Curve aIC(aBSp, H1, H2);
1914                   mySeqOfCurve.Append(aIC);
1915                 } else {
1916                   mySeqOfCurve.Append(aCurve);
1917                 }
1918               } else {
1919                 mySeqOfCurve.Append(aCurve);
1920               }
1921
1922             }// else if(typs2 == GeomAbs_Plane)
1923             //
1924             else { //typs2 != GeomAbs_Plane && typs1 != GeomAbs_Plane
1925               Standard_Boolean bIsValid1, bIsValid2;
1926               Handle(Geom_BSplineCurve) BS;
1927               Handle(Geom2d_BSplineCurve) aH2D;
1928               IntTools_Curve aCurve;
1929               //
1930               bIsValid1=Standard_True;
1931               bIsValid2=Standard_True;
1932               //
1933               const AppParCurves_MultiBSpCurve& mbspc = theapp3d.Value(j);
1934               nbpoles = mbspc.NbPoles();
1935               TColgp_Array1OfPnt tpoles(1,nbpoles);
1936               mbspc.Curve(1,tpoles);
1937               BS=new Geom_BSplineCurve(tpoles,
1938                                                                  mbspc.Knots(),
1939                                                                  mbspc.Multiplicities(),
1940                                                                  mbspc.Degree());
1941               GeomLib_CheckBSplineCurve Check(BS,TOLCHECK,TOLANGCHECK);
1942               Check.FixTangent(Standard_True,Standard_True);
1943               //
1944               aCurve.SetCurve(BS);
1945               aCurve.SetFirstCurve2d(aH2D);
1946               aCurve.SetSecondCurve2d(aH2D);
1947               //
1948               if(myApprox1) {
1949                 if(anApprox1) {
1950                   Handle(Geom2d_BSplineCurve) BS1;
1951                   TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
1952                   mbspc.Curve(2,tpoles2d);
1953                   //
1954                   BS1=new Geom2d_BSplineCurve(tpoles2d,
1955                                                                         mbspc.Knots(),
1956                                                                         mbspc.Multiplicities(),
1957                                                                         mbspc.Degree());
1958                   GeomLib_Check2dBSplineCurve newCheck(BS1,TOLCHECK,TOLANGCHECK);
1959                   newCheck.FixTangent(Standard_True,Standard_True);
1960                   //
1961                   if (!reApprox) {
1962                     bIsValid1=CheckPCurve(BS1, myFace1);
1963                   }
1964                   //
1965                   aCurve.SetFirstCurve2d(BS1);
1966                 }
1967                 else {
1968                   Handle(Geom2d_BSplineCurve) BS1;
1969                   fprm = BS->FirstParameter();
1970                   lprm = BS->LastParameter();
1971
1972                   Handle(Geom2d_Curve) C2d;
1973                   Standard_Real aTol = myTolApprox;
1974                   GeomInt_IntSS::BuildPCurves(fprm, lprm, aTol,
1975                             myHS1->ChangeSurface().Surface(), BS, C2d);
1976                   BS1 = Handle(Geom2d_BSplineCurve)::DownCast(C2d);
1977                   aCurve.SetFirstCurve2d(BS1);
1978                 }
1979               } // if(myApprox1) {
1980                 //
1981               if(myApprox2) {
1982                 if(anApprox2) {
1983                   Handle(Geom2d_BSplineCurve) BS2;
1984                   TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
1985                   mbspc.Curve((myApprox1==Standard_True)? 3 : 2,tpoles2d);
1986                   BS2=new Geom2d_BSplineCurve(tpoles2d,
1987                                                                         mbspc.Knots(),
1988                                                                         mbspc.Multiplicities(),
1989                                                                         mbspc.Degree());
1990                   GeomLib_Check2dBSplineCurve newCheck(BS2,TOLCHECK,TOLANGCHECK);
1991                   newCheck.FixTangent(Standard_True,Standard_True);
1992                 //
1993                   if (!reApprox) {
1994                     bIsValid2=CheckPCurve(BS2, myFace2);
1995                   }
1996                   aCurve.SetSecondCurve2d(BS2);
1997                 }
1998                 else {
1999                   Handle(Geom2d_BSplineCurve) BS2;
2000                   fprm = BS->FirstParameter();
2001                   lprm = BS->LastParameter();
2002
2003                   Handle(Geom2d_Curve) C2d;
2004                   Standard_Real aTol = myTolApprox;
2005                   GeomInt_IntSS::BuildPCurves(fprm, lprm, aTol,
2006                             myHS2->ChangeSurface().Surface(), BS, C2d);
2007                   BS2 = Handle(Geom2d_BSplineCurve)::DownCast(C2d);
2008                   aCurve.SetSecondCurve2d(BS2);
2009                 }
2010               } //if(myApprox2) {
2011               if (!bIsValid1 || !bIsValid2) {
2012                 myTolApprox=aTolApproxImp;//1.e-5;
2013                 tol2d = myTolApprox;
2014                 reApprox = Standard_True;
2015                 goto reapprox;
2016               }
2017                 //
2018               mySeqOfCurve.Append(aCurve);
2019             }
2020           }
2021         }
2022       }
2023     }// else { // X
2024   }// case IntPatch_Walking:{
2025     break;
2026
2027   case IntPatch_Restriction:
2028     {
2029       Handle(IntPatch_RLine) RL =
2030         Handle(IntPatch_RLine)::DownCast(L);
2031
2032 #ifdef INTTOOLS_FACEFACE_DEBUG
2033     RL->Dump(0);
2034 #endif
2035
2036       Handle(Geom_Curve) aC3d;
2037       Handle(Geom2d_Curve) aC2d1, aC2d2;
2038       Standard_Real aTolReached;
2039       GeomInt_IntSS::TreatRLine(RL, myHS1, myHS2, aC3d,
2040                                   aC2d1, aC2d2, aTolReached);
2041
2042       if(aC3d.IsNull())
2043         break;
2044
2045       Bnd_Box2d aBox1, aBox2;
2046
2047       const Standard_Real aU1f = myHS1->FirstUParameter(),
2048                           aV1f = myHS1->FirstVParameter(),
2049                           aU1l = myHS1->LastUParameter(),
2050                           aV1l = myHS1->LastVParameter();
2051       const Standard_Real aU2f = myHS2->FirstUParameter(),
2052                           aV2f = myHS2->FirstVParameter(),
2053                           aU2l = myHS2->LastUParameter(),
2054                           aV2l = myHS2->LastVParameter();
2055
2056       aBox1.Add(gp_Pnt2d(aU1f, aV1f));
2057       aBox1.Add(gp_Pnt2d(aU1l, aV1l));
2058       aBox2.Add(gp_Pnt2d(aU2f, aV2f));
2059       aBox2.Add(gp_Pnt2d(aU2l, aV2l));
2060
2061       GeomInt_VectorOfReal anArrayOfParameters;
2062
2063       //We consider here that the intersection line is same-parameter-line
2064       anArrayOfParameters.Append(aC3d->FirstParameter());
2065       anArrayOfParameters.Append(aC3d->LastParameter());
2066
2067       GeomInt_IntSS::
2068         TrimILineOnSurfBoundaries(aC2d1, aC2d2, aBox1, aBox2, anArrayOfParameters);
2069
2070       //Intersect with true boundaries. After that, enlarge bounding-boxes in order to
2071       //correct definition, if point on curve is inscribed in the box.
2072       aBox1.Enlarge(theToler);
2073       aBox2.Enlarge(theToler);
2074
2075       const Standard_Integer aNbIntersSolutionsm1 = anArrayOfParameters.Length() - 1;
2076
2077       //Trim RLine found.
2078       for(Standard_Integer anInd = 0; anInd < aNbIntersSolutionsm1; anInd++)
2079       {
2080         Standard_Real &aParF = anArrayOfParameters(anInd),
2081                       &aParL = anArrayOfParameters(anInd+1);
2082
2083         if((aParL - aParF) <= Precision::PConfusion())
2084         {
2085           //In order to more precise extending to the boundaries of source curves.
2086           if(anInd < aNbIntersSolutionsm1-1)
2087             aParL = aParF;
2088
2089           continue;
2090         }
2091
2092         const Standard_Real aPar = 0.5*(aParF + aParL);
2093         gp_Pnt2d aPt;
2094
2095         Handle(Geom2d_Curve) aCurv2d1, aCurv2d2;
2096         if(!aC2d1.IsNull())
2097         {
2098           aC2d1->D0(aPar, aPt);
2099
2100           if(aBox1.IsOut(aPt))
2101             continue;
2102
2103           if(myApprox1)
2104             aCurv2d1 = new Geom2d_TrimmedCurve(aC2d1, aParF, aParL);
2105         }
2106
2107         if(!aC2d2.IsNull())
2108         {
2109           aC2d2->D0(aPar, aPt);
2110
2111           if(aBox2.IsOut(aPt))
2112             continue;
2113
2114           if(myApprox2)
2115             aCurv2d2 = new Geom2d_TrimmedCurve(aC2d2, aParF, aParL);
2116         }
2117
2118         Handle(Geom_Curve) aCurv3d = new Geom_TrimmedCurve(aC3d, aParF, aParL);
2119
2120         IntTools_Curve aIC(aCurv3d, aCurv2d1, aCurv2d2);
2121         mySeqOfCurve.Append(aIC);
2122       }
2123     }
2124     break;
2125   default:
2126     break;
2127
2128   }
2129 }
2130
2131 //=======================================================================
2132 //function : Parameters
2133 //purpose  :
2134 //=======================================================================
2135  void Parameters(const Handle(GeomAdaptor_HSurface)& HS1,
2136                  const Handle(GeomAdaptor_HSurface)& HS2,
2137                  const gp_Pnt& Ptref,
2138                  Standard_Real& U1,
2139                  Standard_Real& V1,
2140                  Standard_Real& U2,
2141                  Standard_Real& V2)
2142 {
2143
2144   IntSurf_Quadric quad1,quad2;
2145   GeomAbs_SurfaceType typs = HS1->Surface().GetType();
2146
2147   switch (typs) {
2148   case GeomAbs_Plane:
2149     quad1.SetValue(HS1->Surface().Plane());
2150     break;
2151   case GeomAbs_Cylinder:
2152     quad1.SetValue(HS1->Surface().Cylinder());
2153     break;
2154   case GeomAbs_Cone:
2155     quad1.SetValue(HS1->Surface().Cone());
2156     break;
2157   case GeomAbs_Sphere:
2158     quad1.SetValue(HS1->Surface().Sphere());
2159     break;
2160   case GeomAbs_Torus:
2161     quad1.SetValue(HS1->Surface().Torus());
2162     break;
2163   default:
2164     Standard_ConstructionError::Raise("GeomInt_IntSS::MakeCurve");
2165   }
2166
2167   typs = HS2->Surface().GetType();
2168   switch (typs) {
2169   case GeomAbs_Plane:
2170     quad2.SetValue(HS2->Surface().Plane());
2171     break;
2172   case GeomAbs_Cylinder:
2173     quad2.SetValue(HS2->Surface().Cylinder());
2174     break;
2175   case GeomAbs_Cone:
2176     quad2.SetValue(HS2->Surface().Cone());
2177     break;
2178   case GeomAbs_Sphere:
2179     quad2.SetValue(HS2->Surface().Sphere());
2180     break;
2181   case GeomAbs_Torus:
2182     quad2.SetValue(HS2->Surface().Torus());
2183     break;
2184   default:
2185     Standard_ConstructionError::Raise("GeomInt_IntSS::MakeCurve");
2186   }
2187
2188   quad1.Parameters(Ptref,U1,V1);
2189   quad2.Parameters(Ptref,U2,V2);
2190 }
2191
2192 //=======================================================================
2193 //function : MakeBSpline
2194 //purpose  :
2195 //=======================================================================
2196 Handle(Geom_Curve) MakeBSpline  (const Handle(IntPatch_WLine)& WL,
2197                                  const Standard_Integer ideb,
2198                                  const Standard_Integer ifin)
2199 {
2200   Standard_Integer i,nbpnt = ifin-ideb+1;
2201   TColgp_Array1OfPnt poles(1,nbpnt);
2202   TColStd_Array1OfReal knots(1,nbpnt);
2203   TColStd_Array1OfInteger mults(1,nbpnt);
2204   Standard_Integer ipidebm1;
2205   for(i=1,ipidebm1=i+ideb-1; i<=nbpnt;ipidebm1++, i++) {
2206     poles(i) = WL->Point(ipidebm1).Value();
2207     mults(i) = 1;
2208     knots(i) = i-1;
2209   }
2210   mults(1) = mults(nbpnt) = 2;
2211   return
2212     new Geom_BSplineCurve(poles,knots,mults,1);
2213 }
2214
2215 //=======================================================================
2216 //function : PrepareLines3D
2217 //purpose  :
2218 //=======================================================================
2219   void IntTools_FaceFace::PrepareLines3D(const Standard_Boolean bToSplit)
2220 {
2221   Standard_Integer i, aNbCurves;
2222   GeomAbs_SurfaceType aType1, aType2;
2223   IntTools_SequenceOfCurves aNewCvs;
2224   //
2225   // 1. Treatment closed  curves
2226   aNbCurves=mySeqOfCurve.Length();
2227   for (i=1; i<=aNbCurves; ++i) {
2228     const IntTools_Curve& aIC=mySeqOfCurve(i);
2229     //
2230     if (bToSplit) {
2231       Standard_Integer j, aNbC;
2232       IntTools_SequenceOfCurves aSeqCvs;
2233       //
2234       aNbC=IntTools_Tools::SplitCurve(aIC, aSeqCvs);
2235       if (aNbC) {
2236         for (j=1; j<=aNbC; ++j) {
2237           const IntTools_Curve& aICNew=aSeqCvs(j);
2238           aNewCvs.Append(aICNew);
2239         }
2240       }
2241       else {
2242         aNewCvs.Append(aIC);
2243       }
2244     }
2245     else {
2246       aNewCvs.Append(aIC);
2247     }
2248   }
2249   //
2250   // 2. Plane\Cone intersection when we had 4 curves
2251   aType1=myHS1->GetType();
2252   aType2=myHS2->GetType();
2253   aNbCurves=aNewCvs.Length();
2254   //
2255   if ((aType1==GeomAbs_Plane && aType2==GeomAbs_Cone) ||
2256       (aType2==GeomAbs_Plane && aType1==GeomAbs_Cone)) {
2257     if (aNbCurves==4) {
2258       GeomAbs_CurveType aCType1;
2259       //
2260       aCType1=aNewCvs(1).Type();
2261       if (aCType1==GeomAbs_Line) {
2262         IntTools_SequenceOfCurves aSeqIn, aSeqOut;
2263         //
2264         for (i=1; i<=aNbCurves; ++i) {
2265           const IntTools_Curve& aIC=aNewCvs(i);
2266           aSeqIn.Append(aIC);
2267         }
2268         //
2269         IntTools_Tools::RejectLines(aSeqIn, aSeqOut);
2270         //
2271         aNewCvs.Clear();
2272         aNbCurves=aSeqOut.Length();
2273         for (i=1; i<=aNbCurves; ++i) {
2274           const IntTools_Curve& aIC=aSeqOut(i);
2275           aNewCvs.Append(aIC);
2276         }
2277       }
2278     }
2279   }// if ((aType1==GeomAbs_Plane && aType2==GeomAbs_Cone)...
2280   //
2281   // 3. Fill  mySeqOfCurve
2282   mySeqOfCurve.Clear();
2283   aNbCurves=aNewCvs.Length();
2284   for (i=1; i<=aNbCurves; ++i) {
2285     const IntTools_Curve& aIC=aNewCvs(i);
2286     mySeqOfCurve.Append(aIC);
2287   }
2288 }
2289 //=======================================================================
2290 //function : CorrectSurfaceBoundaries
2291 //purpose  :
2292 //=======================================================================
2293  void CorrectSurfaceBoundaries(const TopoDS_Face&  theFace,
2294                               const Standard_Real theTolerance,
2295                               Standard_Real&      theumin,
2296                               Standard_Real&      theumax,
2297                               Standard_Real&      thevmin,
2298                               Standard_Real&      thevmax)
2299 {
2300   Standard_Boolean enlarge, isuperiodic, isvperiodic;
2301   Standard_Real uinf, usup, vinf, vsup, delta;
2302   GeomAbs_SurfaceType aType;
2303   Handle(Geom_Surface) aSurface;
2304   //
2305   aSurface = BRep_Tool::Surface(theFace);
2306   aSurface->Bounds(uinf, usup, vinf, vsup);
2307   delta = theTolerance;
2308   enlarge = Standard_False;
2309   //
2310   GeomAdaptor_Surface anAdaptorSurface(aSurface);
2311   //
2312   if(aSurface->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface))) {
2313     Handle(Geom_Surface) aBasisSurface =
2314       (Handle(Geom_RectangularTrimmedSurface)::DownCast(aSurface))->BasisSurface();
2315
2316     if(aBasisSurface->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface)) ||
2317        aBasisSurface->IsKind(STANDARD_TYPE(Geom_OffsetSurface))) {
2318       return;
2319     }
2320   }
2321   //
2322   if(aSurface->IsKind(STANDARD_TYPE(Geom_OffsetSurface))) {
2323     Handle(Geom_Surface) aBasisSurface =
2324       (Handle(Geom_OffsetSurface)::DownCast(aSurface))->BasisSurface();
2325
2326     if(aBasisSurface->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface)) ||
2327        aBasisSurface->IsKind(STANDARD_TYPE(Geom_OffsetSurface))) {
2328       return;
2329     }
2330   }
2331   //
2332   isuperiodic = anAdaptorSurface.IsUPeriodic();
2333   isvperiodic = anAdaptorSurface.IsVPeriodic();
2334   //
2335   aType=anAdaptorSurface.GetType();
2336   if((aType==GeomAbs_BezierSurface) ||
2337      (aType==GeomAbs_BSplineSurface) ||
2338      (aType==GeomAbs_SurfaceOfExtrusion) ||
2339      (aType==GeomAbs_SurfaceOfRevolution) ||
2340      (aType==GeomAbs_Cylinder)) {
2341     enlarge=Standard_True;
2342   }
2343   //
2344   if(!isuperiodic && enlarge) {
2345
2346     if(!Precision::IsInfinite(theumin) &&
2347         ((theumin - uinf) > delta))
2348       theumin -= delta;
2349     else {
2350       theumin = uinf;
2351     }
2352
2353     if(!Precision::IsInfinite(theumax) &&
2354         ((usup - theumax) > delta))
2355       theumax += delta;
2356     else
2357       theumax = usup;
2358   }
2359   //
2360   if(!isvperiodic && enlarge) {
2361     if(!Precision::IsInfinite(thevmin) &&
2362         ((thevmin - vinf) > delta)) {
2363       thevmin -= delta;
2364     }
2365     else {
2366       thevmin = vinf;
2367     }
2368     if(!Precision::IsInfinite(thevmax) &&
2369         ((vsup - thevmax) > delta)) {
2370       thevmax += delta;
2371     }
2372     else {
2373       thevmax = vsup;
2374     }
2375   }
2376   //
2377   {
2378     Standard_Integer aNbP;
2379     Standard_Real aXP, dXfact, aXmid, aX1, aX2, aTolPA;
2380     //
2381     aTolPA=Precision::Angular();
2382     // U
2383     if (isuperiodic) {
2384       aXP=anAdaptorSurface.UPeriod();
2385       dXfact=theumax-theumin;
2386       if (dXfact-aTolPA>aXP) {
2387         aXmid=0.5*(theumax+theumin);
2388         aNbP=RealToInt(aXmid/aXP);
2389         if (aXmid<0.) {
2390           aNbP=aNbP-1;
2391         }
2392         aX1=aNbP*aXP;
2393         if (theumin>aTolPA) {
2394           aX1=theumin+aNbP*aXP;
2395         }
2396         aX2=aX1+aXP;
2397         if (theumin<aX1) {
2398           theumin=aX1;
2399         }
2400         if (theumax>aX2) {
2401           theumax=aX2;
2402         }
2403       }
2404     }
2405     // V
2406     if (isvperiodic) {
2407       aXP=anAdaptorSurface.VPeriod();
2408       dXfact=thevmax-thevmin;
2409       if (dXfact-aTolPA>aXP) {
2410         aXmid=0.5*(thevmax+thevmin);
2411         aNbP=RealToInt(aXmid/aXP);
2412         if (aXmid<0.) {
2413           aNbP=aNbP-1;
2414         }
2415         aX1=aNbP*aXP;
2416         if (thevmin>aTolPA) {
2417           aX1=thevmin+aNbP*aXP;
2418         }
2419         aX2=aX1+aXP;
2420         if (thevmin<aX1) {
2421           thevmin=aX1;
2422         }
2423         if (thevmax>aX2) {
2424           thevmax=aX2;
2425         }
2426       }
2427     }
2428   }
2429   //
2430   if(isuperiodic || isvperiodic) {
2431     Standard_Boolean correct = Standard_False;
2432     Standard_Boolean correctU = Standard_False;
2433     Standard_Boolean correctV = Standard_False;
2434     Bnd_Box2d aBox;
2435     TopExp_Explorer anExp;
2436
2437     for(anExp.Init(theFace, TopAbs_EDGE); anExp.More(); anExp.Next()) {
2438       if(BRep_Tool::IsClosed(TopoDS::Edge(anExp.Current()), theFace)) {
2439         correct = Standard_True;
2440         Standard_Real f, l;
2441         TopoDS_Edge anEdge = TopoDS::Edge(anExp.Current());
2442
2443         for(Standard_Integer i = 0; i < 2; i++) {
2444           if(i==0) {
2445             anEdge.Orientation(TopAbs_FORWARD);
2446           }
2447           else {
2448             anEdge.Orientation(TopAbs_REVERSED);
2449           }
2450           Handle(Geom2d_Curve) aCurve = BRep_Tool::CurveOnSurface(anEdge, theFace, f, l);
2451
2452           if(aCurve.IsNull()) {
2453             correct = Standard_False;
2454             break;
2455           }
2456           Handle(Geom2d_Line) aLine = Handle(Geom2d_Line)::DownCast(aCurve);
2457
2458           if(aLine.IsNull()) {
2459             correct = Standard_False;
2460             break;
2461           }
2462           gp_Dir2d anUDir(1., 0.);
2463           gp_Dir2d aVDir(0., 1.);
2464           Standard_Real anAngularTolerance = Precision::Angular();
2465
2466           correctU = correctU || aLine->Position().Direction().IsParallel(aVDir, anAngularTolerance);
2467           correctV = correctV || aLine->Position().Direction().IsParallel(anUDir, anAngularTolerance);
2468
2469           gp_Pnt2d pp1 = aCurve->Value(f);
2470           aBox.Add(pp1);
2471           gp_Pnt2d pp2 = aCurve->Value(l);
2472           aBox.Add(pp2);
2473         }
2474         if(!correct)
2475           break;
2476       }
2477     }
2478
2479     if(correct) {
2480       Standard_Real umin, vmin, umax, vmax;
2481       aBox.Get(umin, vmin, umax, vmax);
2482
2483       if(isuperiodic && correctU) {
2484
2485         if(theumin < umin)
2486           theumin = umin;
2487
2488         if(theumax > umax) {
2489           theumax = umax;
2490         }
2491       }
2492       if(isvperiodic && correctV) {
2493
2494         if(thevmin < vmin)
2495           thevmin = vmin;
2496         if(thevmax > vmax)
2497           thevmax = vmax;
2498       }
2499     }
2500   }
2501 }
2502
2503 //=======================================================================
2504 //function : TolR3d
2505 //purpose  :
2506 //=======================================================================
2507 void TolR3d(const TopoDS_Face& aF1,
2508             const TopoDS_Face& aF2,
2509             Standard_Real& myTolReached3d)
2510 {
2511   Standard_Real aTolF1, aTolF2, aTolFMax, aTolTresh;
2512
2513   aTolTresh=2.999999e-3;
2514   aTolF1 = BRep_Tool::Tolerance(aF1);
2515   aTolF2 = BRep_Tool::Tolerance(aF2);
2516   aTolFMax=Max(aTolF1, aTolF2);
2517
2518   if (aTolFMax>aTolTresh) {
2519     myTolReached3d=aTolFMax;
2520   }
2521 }
2522
2523 // ------------------------------------------------------------------------------------------------
2524 // static function: ParameterOutOfBoundary
2525 // purpose:         Computes a new parameter for given curve. The corresponding 2d points
2526 //                  does not lay on any boundary of given faces
2527 // ------------------------------------------------------------------------------------------------
2528 Standard_Boolean ParameterOutOfBoundary(const Standard_Real       theParameter,
2529                                         const Handle(Geom_Curve)& theCurve,
2530                                         const TopoDS_Face&        theFace1,
2531                                         const TopoDS_Face&        theFace2,
2532                                         const Standard_Real       theOtherParameter,
2533                                         const Standard_Boolean    bIncreasePar,
2534                                         Standard_Real&            theNewParameter,
2535                                         const Handle(IntTools_Context)& aContext)
2536 {
2537   Standard_Boolean bIsComputed = Standard_False;
2538   theNewParameter = theParameter;
2539
2540   Standard_Real acurpar = theParameter;
2541   TopAbs_State aState = TopAbs_ON;
2542   Standard_Integer iter = 0;
2543   Standard_Real asumtol = BRep_Tool::Tolerance(theFace1) + BRep_Tool::Tolerance(theFace2);
2544   Standard_Real adelta = asumtol * 0.1;
2545   adelta = (adelta < Precision::Confusion()) ? Precision::Confusion() : adelta;
2546   Handle(Geom_Surface) aSurf1 = BRep_Tool::Surface(theFace1);
2547   Handle(Geom_Surface) aSurf2 = BRep_Tool::Surface(theFace2);
2548
2549   Standard_Real u1, u2, v1, v2;
2550
2551   GeomAPI_ProjectPointOnSurf aPrj1;
2552   aSurf1->Bounds(u1, u2, v1, v2);
2553   aPrj1.Init(aSurf1, u1, u2, v1, v2);
2554
2555   GeomAPI_ProjectPointOnSurf aPrj2;
2556   aSurf2->Bounds(u1, u2, v1, v2);
2557   aPrj2.Init(aSurf2, u1, u2, v1, v2);
2558
2559   while(aState == TopAbs_ON) {
2560     if(bIncreasePar)
2561       acurpar += adelta;
2562     else
2563       acurpar -= adelta;
2564     gp_Pnt aPCurrent = theCurve->Value(acurpar);
2565     aPrj1.Perform(aPCurrent);
2566     Standard_Real U=0., V=0.;
2567
2568     if(aPrj1.IsDone()) {
2569       aPrj1.LowerDistanceParameters(U, V);
2570       aState = aContext->StatePointFace(theFace1, gp_Pnt2d(U, V));
2571     }
2572
2573     if(aState != TopAbs_ON) {
2574       aPrj2.Perform(aPCurrent);
2575
2576       if(aPrj2.IsDone()) {
2577         aPrj2.LowerDistanceParameters(U, V);
2578         aState = aContext->StatePointFace(theFace2, gp_Pnt2d(U, V));
2579       }
2580     }
2581
2582     if(iter > 11) {
2583       break;
2584     }
2585     iter++;
2586   }
2587
2588   if(iter <= 11) {
2589     theNewParameter = acurpar;
2590     bIsComputed = Standard_True;
2591
2592     if(bIncreasePar) {
2593       if(acurpar >= theOtherParameter)
2594         theNewParameter = theOtherParameter;
2595     }
2596     else {
2597       if(acurpar <= theOtherParameter)
2598         theNewParameter = theOtherParameter;
2599     }
2600   }
2601   return bIsComputed;
2602 }
2603
2604 //=======================================================================
2605 //function : IsCurveValid
2606 //purpose  :
2607 //=======================================================================
2608 Standard_Boolean IsCurveValid(Handle(Geom2d_Curve)& thePCurve)
2609 {
2610   if(thePCurve.IsNull())
2611     return Standard_False;
2612
2613   Standard_Real tolint = 1.e-10;
2614   Geom2dAdaptor_Curve PCA;
2615   IntRes2d_Domain PCD;
2616   Geom2dInt_GInter PCI;
2617
2618   Standard_Real pf = 0., pl = 0.;
2619   gp_Pnt2d pntf, pntl;
2620
2621   if(!thePCurve->IsClosed() && !thePCurve->IsPeriodic()) {
2622     pf = thePCurve->FirstParameter();
2623     pl = thePCurve->LastParameter();
2624     pntf = thePCurve->Value(pf);
2625     pntl = thePCurve->Value(pl);
2626     PCA.Load(thePCurve);
2627     if(!PCA.IsPeriodic()) {
2628       if(PCA.FirstParameter() > pf) pf = PCA.FirstParameter();
2629       if(PCA.LastParameter()  < pl) pl = PCA.LastParameter();
2630     }
2631     PCD.SetValues(pntf,pf,tolint,pntl,pl,tolint);
2632     PCI.Perform(PCA,PCD,tolint,tolint);
2633     if(PCI.IsDone())
2634       if(PCI.NbPoints() > 0) {
2635         return Standard_False;
2636       }
2637   }
2638
2639   return Standard_True;
2640 }
2641
2642 //=======================================================================
2643 //static function : ApproxWithPCurves
2644 //purpose  : for bug 20964 only
2645 //=======================================================================
2646 Standard_Boolean ApproxWithPCurves(const gp_Cylinder& theCyl,
2647                                    const gp_Sphere& theSph)
2648 {
2649   Standard_Boolean bRes = Standard_True;
2650   Standard_Real R1 = theCyl.Radius(), R2 = theSph.Radius();
2651   //
2652   {
2653     Standard_Real aD2, aRc2, aEps;
2654     gp_Pnt aApexSph;
2655     //
2656     aEps=1.E-7;
2657     aRc2=R1*R1;
2658     //
2659     const gp_Ax3& aAx3Sph=theSph.Position();
2660     const gp_Pnt& aLocSph=aAx3Sph.Location();
2661     const gp_Dir& aDirSph=aAx3Sph.Direction();
2662     //
2663     const gp_Ax1& aAx1Cyl=theCyl.Axis();
2664     gp_Lin aLinCyl(aAx1Cyl);
2665     //
2666     aApexSph.SetXYZ(aLocSph.XYZ()+R2*aDirSph.XYZ());
2667     aD2=aLinCyl.SquareDistance(aApexSph);
2668     if (fabs(aD2-aRc2)<aEps) {
2669       return !bRes;
2670     }
2671     //
2672     aApexSph.SetXYZ(aLocSph.XYZ()-R2*aDirSph.XYZ());
2673     aD2=aLinCyl.SquareDistance(aApexSph);
2674     if (fabs(aD2-aRc2)<aEps) {
2675       return !bRes;
2676     }
2677   }
2678   //
2679
2680   if(R1 < 2.*R2) {
2681     return bRes;
2682   }
2683   gp_Lin anCylAx(theCyl.Axis());
2684
2685   Standard_Real aDist = anCylAx.Distance(theSph.Location());
2686   Standard_Real aDRel = Abs(aDist - R1)/R2;
2687
2688   if(aDRel > .2) return bRes;
2689
2690   Standard_Real par = ElCLib::Parameter(anCylAx, theSph.Location());
2691   gp_Pnt aP = ElCLib::Value(par, anCylAx);
2692   gp_Vec aV(aP, theSph.Location());
2693
2694   Standard_Real dd = aV.Dot(theSph.Position().XDirection());
2695
2696   if(aDist < R1 && dd > 0.) return Standard_False;
2697   if(aDist > R1 && dd < 0.) return Standard_False;
2698
2699
2700   return bRes;
2701 }
2702 //=======================================================================
2703 //function : PerformPlanes
2704 //purpose  :
2705 //=======================================================================
2706 void  PerformPlanes(const Handle(GeomAdaptor_HSurface)& theS1,
2707                     const Handle(GeomAdaptor_HSurface)& theS2,
2708                     const Standard_Real TolF1,
2709                     const Standard_Real TolF2,
2710                     const Standard_Real TolAng,
2711                     const Standard_Real TolTang,
2712                     const Standard_Boolean theApprox1,
2713                     const Standard_Boolean theApprox2,
2714                     IntTools_SequenceOfCurves& theSeqOfCurve,
2715                     Standard_Boolean& theTangentFaces,
2716                     Standard_Real& TolReached3d)
2717 {
2718
2719   gp_Pln aPln1 = theS1->Surface().Plane();
2720   gp_Pln aPln2 = theS2->Surface().Plane();
2721
2722   IntAna_QuadQuadGeo aPlnInter(aPln1, aPln2, TolAng, TolTang);
2723
2724   if(!aPlnInter.IsDone()) {
2725     theTangentFaces = Standard_False;
2726     return;
2727   }
2728
2729   IntAna_ResultType aResType = aPlnInter.TypeInter();
2730
2731   if(aResType == IntAna_Same) {
2732     theTangentFaces = Standard_True;
2733     return;
2734   }
2735
2736   theTangentFaces = Standard_False;
2737
2738   if(aResType == IntAna_Empty) {
2739     return;
2740   }
2741
2742   gp_Lin aLin = aPlnInter.Line(1);
2743
2744   ProjLib_Plane aProj;
2745
2746   aProj.Init(aPln1);
2747   aProj.Project(aLin);
2748   gp_Lin2d aLin2d1 = aProj.Line();
2749   //
2750   aProj.Init(aPln2);
2751   aProj.Project(aLin);
2752   gp_Lin2d aLin2d2 = aProj.Line();
2753   //
2754   //classify line2d1 relatively first plane
2755   Standard_Real P11, P12;
2756   Standard_Boolean IsCrossed = ClassifyLin2d(theS1, aLin2d1, TolTang, P11, P12);
2757   if(!IsCrossed) return;
2758   //classify line2d2 relatively second plane
2759   Standard_Real P21, P22;
2760   IsCrossed = ClassifyLin2d(theS2, aLin2d2, TolTang, P21, P22);
2761   if(!IsCrossed) return;
2762
2763   //Analysis of parametric intervals: must have common part
2764
2765   if(P21 >= P12) return;
2766   if(P22 <= P11) return;
2767
2768   Standard_Real pmin, pmax;
2769   pmin = Max(P11, P21);
2770   pmax = Min(P12, P22);
2771
2772   if(pmax - pmin <= TolTang) return;
2773
2774   Handle(Geom_Line) aGLin = new Geom_Line(aLin);
2775
2776   IntTools_Curve aCurve;
2777   Handle(Geom_TrimmedCurve) aGTLin = new Geom_TrimmedCurve(aGLin, pmin, pmax);
2778
2779   aCurve.SetCurve(aGTLin);
2780
2781   if(theApprox1) {
2782     Handle(Geom2d_Line) C2d = new Geom2d_Line(aLin2d1);
2783     aCurve.SetFirstCurve2d(new Geom2d_TrimmedCurve(C2d, pmin, pmax));
2784   }
2785   else {
2786     Handle(Geom2d_Curve) H1;
2787     aCurve.SetFirstCurve2d(H1);
2788   }
2789   if(theApprox2) {
2790     Handle(Geom2d_Line) C2d = new Geom2d_Line(aLin2d2);
2791     aCurve.SetSecondCurve2d(new Geom2d_TrimmedCurve(C2d, pmin, pmax));
2792   }
2793   else {
2794     Handle(Geom2d_Curve) H1;
2795     aCurve.SetFirstCurve2d(H1);
2796   }
2797
2798   theSeqOfCurve.Append(aCurve);
2799   //
2800   // computation of the tolerance reached
2801   Standard_Real anAngle, aDt;
2802   gp_Dir aD1, aD2;
2803   //
2804   aD1 = aPln1.Position().Direction();
2805   aD2 = aPln2.Position().Direction();
2806   anAngle = aD1.Angle(aD2);
2807   //
2808   aDt = IntTools_Tools::ComputeIntRange(TolF1, TolF2, anAngle);
2809   TolReached3d = sqrt(aDt*aDt + TolF1*TolF1);
2810 }
2811
2812 //=======================================================================
2813 //function : ClassifyLin2d
2814 //purpose  :
2815 //=======================================================================
2816 static inline Standard_Boolean INTER(const Standard_Real d1,
2817                                      const Standard_Real d2,
2818                                      const Standard_Real tol)
2819 {
2820   return (d1 > tol && d2 < -tol) ||
2821          (d1 < -tol && d2 > tol) ||
2822          ((d1 <= tol && d1 >= -tol) && (d2 > tol || d2 < -tol)) ||
2823          ((d2 <= tol && d2 >= -tol) && (d1 > tol || d1 < -tol));
2824 }
2825 static inline Standard_Boolean COINC(const Standard_Real d1,
2826                                      const Standard_Real d2,
2827                                      const Standard_Real tol)
2828 {
2829   return (d1 <= tol && d1 >= -tol) && (d2 <= tol && d2 >= -tol);
2830 }
2831 Standard_Boolean ClassifyLin2d(const Handle(GeomAdaptor_HSurface)& theS,
2832                                const gp_Lin2d& theLin2d,
2833                                const Standard_Real theTol,
2834                                Standard_Real& theP1,
2835                                Standard_Real& theP2)
2836
2837 {
2838   Standard_Real xmin, xmax, ymin, ymax, d1, d2, A, B, C;
2839   Standard_Real par[2];
2840   Standard_Integer nbi = 0;
2841
2842   xmin = theS->Surface().FirstUParameter();
2843   xmax = theS->Surface().LastUParameter();
2844   ymin = theS->Surface().FirstVParameter();
2845   ymax = theS->Surface().LastVParameter();
2846
2847   theLin2d.Coefficients(A, B, C);
2848
2849   //xmin, ymin <-> xmin, ymax
2850   d1 = A*xmin + B*ymin + C;
2851   d2 = A*xmin + B*ymax + C;
2852
2853   if(INTER(d1, d2, theTol)) {
2854     //Intersection with boundary
2855     Standard_Real y = -(C + A*xmin)/B;
2856     par[nbi] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmin, y));
2857     nbi++;
2858   }
2859   else if (COINC(d1, d2, theTol)) {
2860     //Coincidence with boundary
2861     par[0] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmin, ymin));
2862     par[1] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmin, ymax));
2863     nbi = 2;
2864   }
2865
2866   if(nbi == 2) {
2867
2868     if(fabs(par[0]-par[1]) > theTol) {
2869       theP1 = Min(par[0], par[1]);
2870       theP2 = Max(par[0], par[1]);
2871       return Standard_True;
2872     }
2873     else return Standard_False;
2874
2875   }
2876
2877   //xmin, ymax <-> xmax, ymax
2878   d1 = d2;
2879   d2 = A*xmax + B*ymax + C;
2880
2881   if(d1 > theTol || d1 < -theTol) {//to avoid checking of
2882                                    //coincidence with the same point
2883     if(INTER(d1, d2, theTol)) {
2884       Standard_Real x = -(C + B*ymax)/A;
2885       par[nbi] = ElCLib::Parameter(theLin2d, gp_Pnt2d(x, ymax));
2886       nbi++;
2887     }
2888     else if (COINC(d1, d2, theTol)) {
2889       par[0] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmin, ymax));
2890       par[1] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmax, ymax));
2891       nbi = 2;
2892     }
2893   }
2894
2895   if(nbi == 2) {
2896
2897     if(fabs(par[0]-par[1]) > theTol) {
2898       theP1 = Min(par[0], par[1]);
2899       theP2 = Max(par[0], par[1]);
2900       return Standard_True;
2901     }
2902     else return Standard_False;
2903
2904   }
2905
2906   //xmax, ymax <-> xmax, ymin
2907   d1 = d2;
2908   d2 = A*xmax + B*ymin + C;
2909
2910   if(d1 > theTol || d1 < -theTol) {
2911     if(INTER(d1, d2, theTol)) {
2912       Standard_Real y = -(C + A*xmax)/B;
2913       par[nbi] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmax, y));
2914       nbi++;
2915     }
2916     else if (COINC(d1, d2, theTol)) {
2917       par[0] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmax, ymax));
2918       par[1] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmax, ymin));
2919       nbi = 2;
2920     }
2921   }
2922
2923   if(nbi == 2) {
2924     if(fabs(par[0]-par[1]) > theTol) {
2925       theP1 = Min(par[0], par[1]);
2926       theP2 = Max(par[0], par[1]);
2927       return Standard_True;
2928     }
2929     else return Standard_False;
2930   }
2931
2932   //xmax, ymin <-> xmin, ymin
2933   d1 = d2;
2934   d2 = A*xmin + B*ymin + C;
2935
2936   if(d1 > theTol || d1 < -theTol) {
2937     if(INTER(d1, d2, theTol)) {
2938       Standard_Real x = -(C + B*ymin)/A;
2939       par[nbi] = ElCLib::Parameter(theLin2d, gp_Pnt2d(x, ymin));
2940       nbi++;
2941     }
2942     else if (COINC(d1, d2, theTol)) {
2943       par[0] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmax, ymin));
2944       par[1] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmin, ymin));
2945       nbi = 2;
2946     }
2947   }
2948
2949   if(nbi == 2) {
2950     if(fabs(par[0]-par[1]) > theTol) {
2951       theP1 = Min(par[0], par[1]);
2952       theP2 = Max(par[0], par[1]);
2953       return Standard_True;
2954     }
2955     else return Standard_False;
2956   }
2957
2958   return Standard_False;
2959
2960 }
2961 //
2962 //=======================================================================
2963 //function : ApproxParameters
2964 //purpose  :
2965 //=======================================================================
2966 void ApproxParameters(const Handle(GeomAdaptor_HSurface)& aHS1,
2967                       const Handle(GeomAdaptor_HSurface)& aHS2,
2968                       Standard_Integer& iDegMin,
2969                       Standard_Integer& iDegMax,
2970                       Standard_Integer& iNbIter)
2971
2972 {
2973   GeomAbs_SurfaceType aTS1, aTS2;
2974
2975   //
2976   iNbIter=0;
2977   iDegMin=4;
2978   iDegMax=8;
2979   //
2980   aTS1=aHS1->Surface().GetType();
2981   aTS2=aHS2->Surface().GetType();
2982   //
2983   // Cylinder/Torus
2984   if ((aTS1==GeomAbs_Cylinder && aTS2==GeomAbs_Torus) ||
2985       (aTS2==GeomAbs_Cylinder && aTS1==GeomAbs_Torus)) {
2986     Standard_Real aRC, aRT, dR, aPC;
2987     gp_Cylinder aCylinder;
2988     gp_Torus aTorus;
2989     //
2990     aPC=Precision::Confusion();
2991     //
2992     aCylinder=(aTS1==GeomAbs_Cylinder)? aHS1->Surface().Cylinder() : aHS2->Surface().Cylinder();
2993     aTorus=(aTS1==GeomAbs_Torus)? aHS1->Surface().Torus() : aHS2->Surface().Torus();
2994     //
2995     aRC=aCylinder.Radius();
2996     aRT=aTorus.MinorRadius();
2997     dR=aRC-aRT;
2998     if (dR<0.) {
2999       dR=-dR;
3000     }
3001     //
3002     if (dR<aPC) {
3003       iDegMax=6;
3004     }
3005   }
3006   if (aTS1==GeomAbs_Cylinder && aTS2==GeomAbs_Cylinder) {
3007     iNbIter=1;
3008   }
3009 }
3010 //=======================================================================
3011 //function : Tolerances
3012 //purpose  :
3013 //=======================================================================
3014 void Tolerances(const Handle(GeomAdaptor_HSurface)& aHS1,
3015                 const Handle(GeomAdaptor_HSurface)& aHS2,
3016                 Standard_Real& aTolTang)
3017 {
3018   GeomAbs_SurfaceType aTS1, aTS2;
3019   //
3020   aTS1=aHS1->Surface().GetType();
3021   aTS2=aHS2->Surface().GetType();
3022   //
3023   // Cylinder/Torus
3024   if ((aTS1==GeomAbs_Cylinder && aTS2==GeomAbs_Torus) ||
3025       (aTS2==GeomAbs_Cylinder && aTS1==GeomAbs_Torus)) {
3026     Standard_Real aRC, aRT, dR, aPC;
3027     gp_Cylinder aCylinder;
3028     gp_Torus aTorus;
3029     //
3030     aPC=Precision::Confusion();
3031     //
3032     aCylinder=(aTS1==GeomAbs_Cylinder)? aHS1->Surface().Cylinder() : aHS2->Surface().Cylinder();
3033     aTorus=(aTS1==GeomAbs_Torus)? aHS1->Surface().Torus() : aHS2->Surface().Torus();
3034     //
3035     aRC=aCylinder.Radius();
3036     aRT=aTorus.MinorRadius();
3037     dR=aRC-aRT;
3038     if (dR<0.) {
3039       dR=-dR;
3040     }
3041     //
3042     if (dR<aPC) {
3043       aTolTang=0.1*aTolTang;
3044     }
3045   }
3046 }
3047 //=======================================================================
3048 //function : SortTypes
3049 //purpose  :
3050 //=======================================================================
3051 Standard_Boolean SortTypes(const GeomAbs_SurfaceType aType1,
3052                            const GeomAbs_SurfaceType aType2)
3053 {
3054   Standard_Boolean bRet;
3055   Standard_Integer aI1, aI2;
3056   //
3057   bRet=Standard_False;
3058   //
3059   aI1=IndexType(aType1);
3060   aI2=IndexType(aType2);
3061   if (aI1<aI2){
3062     bRet=!bRet;
3063   }
3064   return bRet;
3065 }
3066 //=======================================================================
3067 //function : IndexType
3068 //purpose  :
3069 //=======================================================================
3070 Standard_Integer IndexType(const GeomAbs_SurfaceType aType)
3071 {
3072   Standard_Integer aIndex;
3073   //
3074   aIndex=11;
3075   //
3076   if (aType==GeomAbs_Plane) {
3077     aIndex=0;
3078   }
3079   else if (aType==GeomAbs_Cylinder) {
3080     aIndex=1;
3081   }
3082   else if (aType==GeomAbs_Cone) {
3083     aIndex=2;
3084   }
3085   else if (aType==GeomAbs_Sphere) {
3086     aIndex=3;
3087   }
3088   else if (aType==GeomAbs_Torus) {
3089     aIndex=4;
3090   }
3091   else if (aType==GeomAbs_BezierSurface) {
3092     aIndex=5;
3093   }
3094   else if (aType==GeomAbs_BSplineSurface) {
3095     aIndex=6;
3096   }
3097   else if (aType==GeomAbs_SurfaceOfRevolution) {
3098     aIndex=7;
3099   }
3100   else if (aType==GeomAbs_SurfaceOfExtrusion) {
3101     aIndex=8;
3102   }
3103   else if (aType==GeomAbs_OffsetSurface) {
3104     aIndex=9;
3105   }
3106   else if (aType==GeomAbs_OtherSurface) {
3107     aIndex=10;
3108   }
3109   return aIndex;
3110 }
3111
3112 //=======================================================================
3113 // Function : FindMaxDistance
3114 // purpose :
3115 //=======================================================================
3116 Standard_Real FindMaxDistance(const Handle(Geom_Curve)& theCurve,
3117                               const Standard_Real theFirst,
3118                               const Standard_Real theLast,
3119                               const TopoDS_Face& theFace,
3120                               const Handle(IntTools_Context)& theContext)
3121 {
3122   Standard_Integer aNbS;
3123   Standard_Real aT1, aT2, aDt, aD, aDMax, anEps;
3124   //
3125   aNbS = 11;
3126   aDt = (theLast - theFirst) / aNbS;
3127   aDMax = 0.;
3128   anEps = 1.e-4 * aDt;
3129   //
3130   GeomAPI_ProjectPointOnSurf& aProjPS = theContext->ProjPS(theFace);
3131   aT2 = theFirst;
3132   for (;;) {
3133     aT1 = aT2;
3134     aT2 += aDt;
3135     //
3136     if (aT2 > theLast) {
3137       break;
3138     }
3139     //
3140     aD = FindMaxDistance(theCurve, aT1, aT2, aProjPS, anEps);
3141     if (aD > aDMax) {
3142       aDMax = aD;
3143     }
3144   }
3145   //
3146   return aDMax;
3147 }
3148
3149 //=======================================================================
3150 // Function : FindMaxDistance
3151 // purpose :
3152 //=======================================================================
3153 Standard_Real FindMaxDistance(const Handle(Geom_Curve)& theC,
3154                               const Standard_Real theFirst,
3155                               const Standard_Real theLast,
3156                               GeomAPI_ProjectPointOnSurf& theProjPS,
3157                               const Standard_Real theEps)
3158 {
3159   Standard_Real aA, aB, aCf, aX, aX1, aX2, aF1, aF2, aF;
3160   //
3161   aCf = 0.61803398874989484820458683436564;//(sqrt(5.)-1)/2.;
3162   aA = theFirst;
3163   aB = theLast;
3164   //
3165   aX1 = aB - aCf * (aB - aA);
3166   aF1 = MaxDistance(theC, aX1, theProjPS);
3167   aX2 = aA + aCf * (aB - aA);
3168   aF2 = MaxDistance(theC, aX2, theProjPS);
3169   //
3170   for (;;) {
3171     if ((aB - aA) < theEps) {
3172       break;
3173     }
3174     //
3175     if (aF1 > aF2) {
3176       aB = aX2;
3177       aX2 = aX1;
3178       aF2 = aF1;
3179       aX1 = aB - aCf * (aB - aA);
3180       aF1 = MaxDistance(theC, aX1, theProjPS);
3181     }
3182     else {
3183       aA = aX1;
3184       aX1 = aX2;
3185       aF1 = aF2;
3186       aX2 = aA + aCf * (aB - aA);
3187       aF2 = MaxDistance(theC, aX2, theProjPS);
3188     }
3189   }
3190   //
3191   aX = 0.5 * (aA + aB);
3192   aF = MaxDistance(theC, aX, theProjPS);
3193   //
3194   if (aF1 > aF) {
3195     aF = aF1;
3196   }
3197   //
3198   if (aF2 > aF) {
3199     aF = aF2;
3200   }
3201   //
3202   return aF;
3203 }
3204
3205 //=======================================================================
3206 // Function : MaxDistance
3207 // purpose :
3208 //=======================================================================
3209 Standard_Real MaxDistance(const Handle(Geom_Curve)& theC,
3210                           const Standard_Real aT,
3211                           GeomAPI_ProjectPointOnSurf& theProjPS)
3212 {
3213   Standard_Real aD;
3214   gp_Pnt aP;
3215   //
3216   theC->D0(aT, aP);
3217   theProjPS.Perform(aP);
3218   aD = theProjPS.NbPoints() ? theProjPS.LowerDistance() : 0.;
3219   //
3220   return aD;
3221 }
3222
3223 //=======================================================================
3224 //function : CheckPCurve
3225 //purpose  : Checks if points of the pcurve are out of the face bounds.
3226 //=======================================================================
3227   Standard_Boolean CheckPCurve(const Handle(Geom2d_Curve)& aPC,
3228                                const TopoDS_Face& aFace)
3229 {
3230   const Standard_Integer NPoints = 23;
3231   Standard_Integer i;
3232   Standard_Real umin,umax,vmin,vmax;
3233
3234   BRepTools::UVBounds(aFace, umin, umax, vmin, vmax);
3235   Standard_Real tolU = Max ((umax-umin)*0.01, Precision::Confusion());
3236   Standard_Real tolV = Max ((vmax-vmin)*0.01, Precision::Confusion());
3237   Standard_Real fp = aPC->FirstParameter();
3238   Standard_Real lp = aPC->LastParameter();
3239
3240
3241   // adjust domain for periodic surfaces
3242   TopLoc_Location aLoc;
3243   Handle(Geom_Surface) aSurf = BRep_Tool::Surface(aFace, aLoc);
3244   if (aSurf->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface))) {
3245     aSurf = (Handle(Geom_RectangularTrimmedSurface)::DownCast(aSurf))->BasisSurface();
3246   }
3247   gp_Pnt2d pnt = aPC->Value((fp+lp)/2);
3248   Standard_Real u,v;
3249   pnt.Coord(u,v);
3250   //
3251   if (aSurf->IsUPeriodic()) {
3252     Standard_Real aPer = aSurf->UPeriod();
3253     Standard_Integer nshift = (Standard_Integer) ((u-umin)/aPer);
3254     if (u < umin+aPer*nshift) nshift--;
3255     umin += aPer*nshift;
3256     umax += aPer*nshift;
3257   }
3258   if (aSurf->IsVPeriodic()) {
3259     Standard_Real aPer = aSurf->VPeriod();
3260     Standard_Integer nshift = (Standard_Integer) ((v-vmin)/aPer);
3261     if (v < vmin+aPer*nshift) nshift--;
3262     vmin += aPer*nshift;
3263     vmax += aPer*nshift;
3264   }
3265   //
3266   //--------------------------------------------------------
3267   Standard_Boolean bRet;
3268   Standard_Integer j, aNbIntervals;
3269   Standard_Real aT, dT;
3270   gp_Pnt2d aP2D;
3271   //
3272   Geom2dAdaptor_Curve aGAC(aPC);
3273   aNbIntervals=aGAC.NbIntervals(GeomAbs_CN);
3274   //
3275   TColStd_Array1OfReal aTI(1, aNbIntervals+1);
3276   aGAC.Intervals(aTI,GeomAbs_CN);
3277   //
3278   bRet=Standard_False;
3279   //
3280   aT=aGAC.FirstParameter();
3281   for (j=1; j<=aNbIntervals; ++j) {
3282     dT=(aTI(j+1)-aTI(j))/NPoints;
3283     //
3284     for (i=1; i<NPoints; i++) {
3285       aT=aT+dT;
3286       aGAC.D0(aT, aP2D);
3287       aP2D.Coord(u,v);
3288     if (umin-u > tolU || u-umax > tolU ||
3289           vmin-v > tolV || v-vmax > tolV) {
3290         return bRet;
3291   }
3292 }
3293   }
3294   return !bRet;
3295 }
3296 //=======================================================================
3297 //function : CorrectPlaneBoundaries
3298 //purpose  :
3299 //=======================================================================
3300  void CorrectPlaneBoundaries(Standard_Real& aUmin,
3301                              Standard_Real& aUmax,
3302                              Standard_Real& aVmin,
3303                              Standard_Real& aVmax)
3304 {
3305   if (!(Precision::IsInfinite(aUmin) ||
3306         Precision::IsInfinite(aUmax))) {
3307     Standard_Real dU;
3308     //
3309     dU=0.1*(aUmax-aUmin);
3310     aUmin=aUmin-dU;
3311     aUmax=aUmax+dU;
3312   }
3313   if (!(Precision::IsInfinite(aVmin) ||
3314         Precision::IsInfinite(aVmax))) {
3315     Standard_Real dV;
3316     //
3317     dV=0.1*(aVmax-aVmin);
3318     aVmin=aVmin-dV;
3319     aVmax=aVmax+dV;
3320   }
3321 }