src/IntTools/IntTools_FaceFace.cxx

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