src/IntTools/IntTools_FaceFace.cxx

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