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