src/IntTools/IntTools_FaceFace.cxx

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