src/IntCurveSurface/IntCurveSurface_Inter.gxx

   1 // Created on: 1993-04-09
   2 // Created by: Laurent BUCHARD
   3 // Copyright (c) 1993-1999 Matra Datavision
   4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
   5 //
   6 // This file is part of Open CASCADE Technology software library.
   7 //
   8 // This library is free software; you can redistribute it and/or modify it under
   9 // the terms of the GNU Lesser General Public License version 2.1 as published
  10 // by the Free Software Foundation, with special exception defined in the file
  11 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
  12 // distribution for complete text of the license and disclaimer of any warranty.
  13 //
  14 // Alternatively, this file may be used under the terms of Open CASCADE
  15 // commercial license or contractual agreement.
  16
  17 //#ifndef DEB
  18 //#define No_Standard_RangeError
  19 //#define No_Standard_OutOfRange
  20 //#endif
  21
  22
  23 #define  TOLTANGENCY         0.00000001
  24 #define  TOLERANCE_ANGULAIRE 1.e-12//0.00000001
  25 #define  TOLERANCE           0.00000001
  26
  27 #define NBSAMPLESONCIRCLE  32
  28 #define NBSAMPLESONELLIPSE 32
  29 #define NBSAMPLESONPARAB   16
  30 #define NBSAMPLESONHYPR    32
  31
  32 #include <ElSLib.hxx>
  33 #include <Intf_SectionPoint.hxx>
  34 #include <Intf_TangentZone.hxx>
  35 #include <Intf_Tool.hxx>
  36 #include <math_FunctionSetRoot.hxx>
  37 #include <IntCurveSurface_IntersectionPoint.hxx>
  38 #include <IntCurveSurface_IntersectionSegment.hxx>
  39 #include <IntAna_IntConicQuad.hxx>
  40 #include <IntAna_Quadric.hxx>
  41 #include <gp_Vec.hxx>
  42 #include <gp_Dir.hxx>
  43 #include <Precision.hxx>
  44 #include <IntAna_IntLinTorus.hxx>
  45
  46 #include <GeomAbs_Shape.hxx>
  47 #include <GeomAbs_CurveType.hxx>
  48 #include <TColStd_Array1OfReal.hxx>
  49 #include <TColgp_Array1OfPnt.hxx>
  50 #include <gp_Pnt.hxx>
  51 #include <gp_Pln.hxx>
  52 #include <gp_Lin.hxx>
  53 #include <GProp_PEquation.hxx>
  54 #include <GProp_PGProps.hxx>
  55 #include <GProp_PrincipalProps.hxx>
  56 #include <Extrema_ExtElC.hxx>
  57 #include <Extrema_POnCurv.hxx>
  58
  59
  60
  61 #include <ProjLib_Plane.hxx>
  62 #include <IntAna2d_AnaIntersection.hxx>
  63 #include <gp_Lin2d.hxx>
  64 #include <IntAna2d_IntPoint.hxx>
  65 #include <gp_Parab2d.hxx>
  66 #include <IntAna2d_Conic.hxx>
  67 #include <IntAna2d_AnaIntersection.hxx>
  68 #include <gp_Hypr2d.hxx>
  69 #include <Adaptor3d_HCurve.hxx>
  70 #include <Adaptor3d_HSurface.hxx>
  71
  72
  73 #include <TColgp_Array2OfPnt.hxx>
  74 #include <TColStd_HArray1OfReal.hxx>
  75 #include <Adaptor3d_TopolTool.hxx>
  76 #include <ElCLib.hxx>
  77 //#define ICS_DEB
  78 static
  79   void EstLimForInfExtr(const gp_Lin&   Line,
  80                         const TheSurface& surface,
  81                         const Standard_Boolean IsOffSurf,
  82                         const Standard_Integer nbsu,
  83                         const Standard_Boolean U1inf,
  84                         const Standard_Boolean U2inf,
  85                         const Standard_Boolean V1inf,
  86                         const Standard_Boolean V2inf,
  87                         Standard_Real& U1new,
  88                         Standard_Real& U2new,
  89                         Standard_Real& V1new,
  90                         Standard_Real& V2new,
  91                         Standard_Boolean& NoIntersection);
  92
  93 static
  94   void EstLimForInfRevl(const gp_Lin&   Line,
  95                         const TheSurface& surface,
  96                         const Standard_Boolean U1inf,
  97                         const Standard_Boolean U2inf,
  98                         const Standard_Boolean V1inf,
  99                         const Standard_Boolean V2inf,
 100                         Standard_Real& U1new,
 101                         Standard_Real& U2new,
 102                         Standard_Real& V1new,
 103                         Standard_Real& V2new,
 104                         Standard_Boolean& NoIntersection);
 105
 106 static
 107   void EstLimForInfOffs(const gp_Lin&   Line,
 108                         const TheSurface& surface,
 109                         const Standard_Integer nbsu,
 110                         const Standard_Boolean U1inf,
 111                         const Standard_Boolean U2inf,
 112                         const Standard_Boolean V1inf,
 113                         const Standard_Boolean V2inf,
 114                         Standard_Real& U1new,
 115                         Standard_Real& U2new,
 116                         Standard_Real& V1new,
 117                         Standard_Real& V2new,
 118                         Standard_Boolean& NoIntersection);
 119
 120 static
 121   void EstLimForInfSurf(Standard_Real& U1new,
 122                         Standard_Real& U2new,
 123                         Standard_Real& V1new,
 124                         Standard_Real& V2new);
 125
 126 static
 127   void SectionPointToParameters(const Intf_SectionPoint& Sp,
 128                                 const IntCurveSurface_ThePolyhedron& Surf,
 129                                 const IntCurveSurface_ThePolygon&    Curv,
 130                                 Standard_Real& u,
 131                                 Standard_Real& v,
 132                                 Standard_Real& w);
 133
 134 static
 135   void IntCurveSurface_ComputeTransitions(const TheCurve& curve,
 136                                           const Standard_Real w,
 137                                           IntCurveSurface_TransitionOnCurve&   TransOnCurve,
 138                                           const TheSurface& surface,
 139                                           const Standard_Real u,
 140                                           const Standard_Real v);
 141
 142 static
 143   void IntCurveSurface_ComputeParamsOnQuadric(const TheSurface& surface,
 144                                               const gp_Pnt& P,
 145                                               Standard_Real& u,
 146                                               Standard_Real& v);
 147
 148 static
 149   void ProjectIntersectAndEstLim(const gp_Lin&        theLine,
 150                                  const gp_Pln&        thePln,
 151                                  const ProjLib_Plane& theBasCurvProj,
 152                                  Standard_Real&       theVmin,
 153                                  Standard_Real&       theVmax,
 154                                  Standard_Boolean&    theNoIntersection);
 155
 156 //=======================================================================
 157 //function : IntCurveSurface_Inter
 158 //purpose  :
 159 //=======================================================================
 160 IntCurveSurface_Inter::IntCurveSurface_Inter()
 161 {
 162 }
 163 //=======================================================================
 164 //function : DoSurface
 165 //purpose  :
 166 //=======================================================================
 167 void IntCurveSurface_Inter::DoSurface(const TheSurface&   surface,
 168                                       const Standard_Real u0,
 169                                       const Standard_Real u1,
 170                                       const Standard_Real v0,
 171                                       const Standard_Real v1,
 172                                       TColgp_Array2OfPnt& pntsOnSurface,
 173                                       Bnd_Box&            boxSurface,
 174                                       Standard_Real&      gap)
 175 {
 176   Standard_Integer iU = 0, iV = 0;
 177   Standard_Real U = 0., V = 0;
 178 // modified by NIZHNY-MKK  Mon Oct  3 17:38:45 2005
 179 //   Standard_Real dU = fabs(u1-u0)/50., dV = fabs(v1-v0)/50.;
 180   Standard_Real dU = (u1-u0)/50., dV = (v1-v0)/50.;
 181   gp_Pnt aPnt;
 182
 183   for(iU = 0; iU < 50; iU++) {
 184
 185     if(iU == 0)
 186       U = u0;
 187     else if(iU == 49)
 188       U = u1;
 189     else
 190       U = u0 + dU * ((Standard_Real)iU);
 191
 192     for(iV = 0; iV < 50; iV++) {
 193
 194       if(iV == 0)
 195         V = v0;
 196       else if(iV == 49)
 197         V = v1;
 198       else
 199         V = v0 + dV * ((Standard_Real)iV);
 200
 201       TheSurfaceTool::D0(surface,U,V,aPnt);
 202       boxSurface.Add(aPnt);
 203       pntsOnSurface.SetValue(iU+1,iV+1,aPnt);
 204     }
 205   }
 206   Standard_Real Ures = TheSurfaceTool::UResolution(surface,dU);
 207   Standard_Real Vres = TheSurfaceTool::VResolution(surface,dV);
 208   gap = Max(Ures,Vres);
 209 }
 210
 211 //=======================================================================
 212 //function : DoNewBounds
 213 //purpose  :
 214 //=======================================================================
 215 void IntCurveSurface_Inter::DoNewBounds(
 216                                         const TheSurface&           surface,
 217                                         const Standard_Real         u0,
 218                                         const Standard_Real         u1,
 219                                         const Standard_Real         v0,
 220                                         const Standard_Real         v1,
 221                                         const TColgp_Array2OfPnt&   pntsOnSurface,
 222                                         const TColStd_Array1OfReal& X,
 223                                         const TColStd_Array1OfReal& Y,
 224                                         const TColStd_Array1OfReal& Z,
 225                                         TColStd_Array1OfReal&       Bounds)
 226 {
 227   Bounds.SetValue(1,u0);
 228   Bounds.SetValue(2,u1);
 229   Bounds.SetValue(3,v0);
 230   Bounds.SetValue(4,v1);
 231
 232   Standard_Boolean isUClosed = (TheSurfaceTool::IsUClosed(surface) || TheSurfaceTool::IsUPeriodic(surface));
 233   Standard_Boolean isVClosed = (TheSurfaceTool::IsVClosed(surface) || TheSurfaceTool::IsVPeriodic(surface));
 234   Standard_Boolean checkU = (isUClosed) ? Standard_False : Standard_True;
 235   Standard_Boolean checkV = (isVClosed) ? Standard_False : Standard_True;
 236
 237   Standard_Integer i = 0, j = 0, k = 0, iU = 0, iV = 0;
 238   Standard_Integer iUmin = 50, iVmin = 50, iUmax = 1, iVmax = 1;
 239
 240   for(i = 1; i <= 2; i++) {
 241     for(j = 1; j <= 2; j++) {
 242       for(k = 1; k <= 2; k++) {
 243         gp_Pnt aPoint(X(i),Y(j),Z(k));
 244         Standard_Real DistMin = 1.e+100;
 245         Standard_Integer diU = 0, diV = 0;
 246         for(iU = 1; iU <= 50; iU++) {
 247           for(iV = 1; iV <= 50; iV++) {
 248             const gp_Pnt aP = pntsOnSurface.Value(iU,iV);
 249             Standard_Real dist = aP.SquareDistance(aPoint);
 250             if(dist < DistMin) {
 251               DistMin = dist;
 252               diU = iU;
 253               diV = iV;
 254             }
 255           }
 256         }
 257         if(diU > 0 && diU < iUmin) iUmin = diU;
 258         if(diU > 0 && diU > iUmax) iUmax = diU;
 259         if(diV > 0 && diV < iVmin) iVmin = diV;
 260         if(diV > 0 && diV > iVmax) iVmax = diV;
 261       }
 262     }
 263   }
 264
 265 // modified by NIZHNY-MKK  Mon Oct  3 17:38:43 2005
 266 //   Standard_Real dU = fabs(u1-u0)/50., dV = fabs(v1-v0)/50.;
 267   Standard_Real dU = (u1-u0)/50., dV = (v1-v0)/50.;
 268
 269   Standard_Real USmin = u0 + dU * ((Standard_Real)(iUmin - 1));
 270   Standard_Real USmax = u0 + dU * ((Standard_Real)(iUmax - 1));
 271   Standard_Real VSmin = v0 + dV * ((Standard_Real)(iVmin - 1));
 272   Standard_Real VSmax = v0 + dV * ((Standard_Real)(iVmax - 1));
 273
 274   if(USmin > USmax) {
 275     Standard_Real tmp = USmax;
 276     USmax = USmin;
 277     USmin = tmp;
 278   }
 279   if(VSmin > VSmax) {
 280     Standard_Real tmp = VSmax;
 281     VSmax = VSmin;
 282     VSmin = tmp;
 283   }
 284
 285   USmin -= 1.5*dU;
 286   if(USmin < u0) USmin = u0;
 287   USmax += 1.5*dU;
 288   if(USmax > u1) USmax = u1;
 289   VSmin -= 1.5*dV;
 290   if(VSmin < v0) VSmin = v0;
 291   VSmax += 1.5*dV;
 292   if(VSmax > v1) VSmax = v1;
 293
 294   if(checkU) {
 295     Bounds.SetValue(1,USmin);
 296     Bounds.SetValue(2,USmax);
 297   }
 298   if(checkV) {
 299     Bounds.SetValue(3,VSmin);
 300     Bounds.SetValue(4,VSmax);
 301   }
 302 }
 303
 304 //=======================================================================
 305 //function : Perform
 306 //purpose  : Decompose la surface si besoin est
 307 //=======================================================================
 308 void IntCurveSurface_Inter::Perform(const TheCurve&   curve,
 309                                     const TheSurface& surface) {
 310   ResetFields();
 311   done = Standard_True;
 312   Standard_Integer NbUOnS = TheSurfaceTool::NbUIntervals(surface,GeomAbs_C2);
 313   Standard_Integer NbVOnS = TheSurfaceTool::NbVIntervals(surface,GeomAbs_C2);
 314   Standard_Real U0,U1,V0,V1;
 315
 316   if(NbUOnS > 1) {
 317     TColStd_Array1OfReal TabU(1,NbUOnS+1);
 318     TheSurfaceTool::UIntervals(surface,TabU,GeomAbs_C2);
 319     for(Standard_Integer iu = 1;iu <= NbUOnS; iu++) {
 320       U0 = TabU.Value(iu);
 321       U1 = TabU.Value(iu+1);
 322       if(NbVOnS > 1) {
 323         TColStd_Array1OfReal TabV(1,NbVOnS+1);
 324         TheSurfaceTool::VIntervals(surface,TabV,GeomAbs_C2);
 325         for(Standard_Integer iv = 1;iv <= NbVOnS; iv++) {
 326           V0 = TabV.Value(iv);
 327           V1 = TabV.Value(iv+1);
 328           Perform(curve,surface,U0,V0,U1,V1);
 329         }
 330       }
 331       else { //------ composite en U     non composite en V
 332         V0 = TheSurfaceTool::FirstVParameter(surface);
 333         V1 = TheSurfaceTool::LastVParameter(surface);
 334         Perform(curve,surface,U0,V0,U1,V1);
 335       }
 336     }
 337   }
 338   else if(NbVOnS > 1) { //---------- non composite en U  composite en V
 339     U0 = TheSurfaceTool::FirstUParameter(surface);
 340     U1 = TheSurfaceTool::LastUParameter(surface);
 341     TColStd_Array1OfReal TabV(1,NbVOnS+1);
 342     TheSurfaceTool::VIntervals(surface,TabV,GeomAbs_C2);
 343     for(Standard_Integer iv = 1;iv <= NbVOnS; iv++) {
 344       V0 = TabV.Value(iv);
 345       V1 = TabV.Value(iv+1);
 346       Perform(curve,surface,U0,V0,U1,V1);
 347     }
 348   }
 349   else {
 350     V0 = TheSurfaceTool::FirstVParameter(surface);
 351     V1 = TheSurfaceTool::LastVParameter(surface);
 352     U0 = TheSurfaceTool::FirstUParameter(surface);
 353     U1 = TheSurfaceTool::LastUParameter(surface);
 354
 355
 356 #if 0
 357     //-- jgv patch (from)
 358     //Computing of local bounds
 359     Standard_Real LocalU0 = U0, LocalU1 = U1, LocalV0 = V0, LocalV1 = V1;
 360     Standard_Real Umin = RealLast(), Vmin = RealLast(), Umax = RealFirst(), Vmax = RealFirst();
 361     Bnd_Box CurveBox;
 362     Standard_Integer i, j, k;
 363     //Making GeomAdaptor_Curve
 364     Standard_Real f = TheCurveTool::FirstParameter(curve);
 365     Standard_Real l = TheCurveTool::LastParameter(curve);
 366     GeomAbs_CurveType CurveType = TheCurveTool::GetType(curve);
 367     Handle(Geom_Curve) theCurve;
 368     switch (CurveType)
 369       {
 370       case GeomAbs_Line:
 371         theCurve = new Geom_Line( TheCurveTool::Line(curve) );
 372         break;
 373       case GeomAbs_Circle:
 374         theCurve = new Geom_Circle( TheCurveTool::Circle(curve) );
 375         break;
 376       case GeomAbs_Ellipse:
 377         theCurve = new Geom_Ellipse( TheCurveTool::Ellipse(curve) );
 378         break;
 379       case GeomAbs_Hyperbola:
 380         theCurve = new Geom_Hyperbola( TheCurveTool::Hyperbola(curve) );
 381         break;
 382       case GeomAbs_Parabola:
 383         theCurve = new Geom_Parabola( TheCurveTool::Parabola(curve) );
 384         break;
 385       case GeomAbs_BezierCurve:
 386         theCurve = TheCurveTool::Bezier(curve);
 387         break;
 388       case GeomAbs_BSplineCurve:
 389         theCurve = TheCurveTool::BSpline(curve);
 390         break;
 391       }
 392     if (!theCurve.IsNull())
 393       {
 394         GeomAdaptor_Curve GACurve( theCurve, f, l );
 395         BndLib_Add3dCurve::Add( GACurve, Precision::Confusion(), CurveBox );
 396       }
 397     else
 398       {
 399         Standard_Integer nbp = TheCurveTool::NbSamples(curve,f,l);
 400         Standard_Real delta = (f-l)/nbp;
 401         for (i = 0; i <= nbp; i++)
 402           {
 403             gp_Pnt aPnt = TheCurveTool::Value(curve, f + i*delta);
 404             CurveBox.Add(aPnt);
 405           }
 406       }
 407     Standard_Real X[2], Y[2], Z[2];
 408     CurveBox.Get( X[0], Y[0], Z[0], X[1], Y[1], Z[1] );
 409     Standard_Real Ures = TheSurfaceTool::UResolution(surface, Precision::Confusion());
 410     Standard_Real Vres = TheSurfaceTool::VResolution(surface, Precision::Confusion());
 411     //Making GeomAdaptor_Surface
 412     GeomAbs_SurfaceType SurfaceType = TheSurfaceTool::GetType(surface);
 413     Handle(Geom_Surface) theSurface;
 414     switch (SurfaceType)
 415       {
 416       case GeomAbs_Plane:
 417         theSurface = new Geom_Plane( TheSurfaceTool::Plane(surface) );
 418         break;
 419       case GeomAbs_Cylinder:
 420         theSurface = new Geom_CylindricalSurface( TheSurfaceTool::Cylinder(surface) );
 421         break;
 422       case GeomAbs_Cone:
 423         theSurface = new Geom_ConicalSurface( TheSurfaceTool::Cone(surface) );
 424         break;
 425       case GeomAbs_Torus:
 426         theSurface = new Geom_ToroidalSurface( TheSurfaceTool::Torus(surface) );
 427         break;
 428       case GeomAbs_Sphere:
 429         theSurface = new Geom_SphericalSurface( TheSurfaceTool::Sphere(surface) );
 430         break;
 431       case GeomAbs_BezierSurface:
 432         theSurface = TheSurfaceTool::Bezier(surface);
 433         break;
 434       case GeomAbs_BSplineSurface:
 435         theSurface = TheSurfaceTool::BSpline(surface);
 436         break;
 437       case GeomAbs_SurfaceOfRevolution:
 438         {
 439           gp_Ax1 Axis = TheSurfaceTool::AxeOfRevolution(surface);
 440           Handle(Adaptor3d_HCurve) AdBC = TheSurfaceTool::BasisCurve(surface);
 441           Handle(Geom_Curve) BC = GetCurve(AdBC);
 442           if (!BC.IsNull())
 443             theSurface = new Geom_SurfaceOfRevolution( BC, Axis );
 444           break;
 445         }
 446       case GeomAbs_SurfaceOfExtrusion:
 447         {
 448           gp_Dir Direction = TheSurfaceTool::Direction(surface);
 449           Handle(Adaptor3d_HCurve) AdBC = TheSurfaceTool::BasisCurve(surface);
 450           Handle(Geom_Curve) BC = GetCurve(AdBC);
 451           if (!BC.IsNull())
 452             theSurface = new Geom_SurfaceOfLinearExtrusion( BC, Direction );
 453           break;
 454         }
 455       case GeomAbs_OffsetSurface:
 456         {
 457           Standard_Real OffsetValue = TheSurfaceTool::OffsetValue(surface);
 458           Handle(Adaptor3d_HSurface) AdBS = TheSurfaceTool::BasisSurface(surface);
 459           Handle(Geom_Surface) BS = GetSurface(AdBS);
 460           if (!BS.IsNull())
 461             theSurface = new Geom_OffsetSurface( BS, OffsetValue );
 462           break;
 463         }
 464       }
 465     if (!theSurface.IsNull())
 466       {
 467         GeomAdaptor_Surface GASurface( theSurface );
 468         Extrema_ExtPS Projector;
 469         Projector.Initialize( GASurface, U0, U1, V0, V1, Ures, Vres );
 470         for (i = 0; i <= 1; i++)
 471           for (j = 0; j <= 1; j++)
 472             for (k = 0; k <= 1; k++)
 473               {
 474                 gp_Pnt aPoint( X[i], Y[j], Z[k] );
 475                 Projector.Perform( aPoint );
 476                 if (Projector.IsDone() && Projector.NbExt() > 0)
 477                   {
 478                     Standard_Real mindist = RealLast();
 479                     Standard_Integer minind, ind;
 480                     for (ind = 1; ind <= Projector.NbExt(); ind++)
 481                       {
 482                         Standard_Real dist = Projector.Value(ind);
 483                         if (dist < mindist)
 484                           {
 485                             mindist = dist;
 486                             minind  = ind;
 487                           }
 488                       }
 489                     Extrema_POnSurf pons = Projector.Point(minind);
 490                     Standard_Real U, V;
 491                     pons.Parameter( U, V );
 492                     if (U < Umin) Umin = U;
 493                     if (U > Umax) Umax = U;
 494                     if (V < Vmin) Vmin = V;
 495                     if (V > Vmax) Vmax = V;
 496                   }
 497               }
 498         Umin -= Ures; Umax += Ures; Vmin -= Vres; Vmax += Vres;
 499         if (Umin > U0 && Umin <= U1) LocalU0 = Umin;
 500         if (Umax < U1 && Umax >= U0) LocalU1 = Umax;
 501         if (Vmin > V0 && Vmin <= V1) LocalV0 = Vmin;
 502         if (Vmax < V1 && Vmax >= V0) LocalV1 = Vmax;
 503         U0 = LocalU0; U1 = LocalU1; V0 = LocalV0; V1 = LocalV1;
 504       }
 505     //-- jgv patch (to)
 506 #endif
 507     Perform(curve,surface,U0,V0,U1,V1);
 508   }
 509 }
 510 //=======================================================================
 511 //function : Perform
 512 //purpose  :
 513 //=======================================================================
 514 void IntCurveSurface_Inter::Perform(const TheCurve&   curve,
 515                                     const TheSurface& surface,
 516                                     const Standard_Real U1,const Standard_Real V1,
 517                                     const Standard_Real U2,const Standard_Real V2) {
 518
 519   GeomAbs_CurveType CurveType = TheCurveTool::GetType(curve);
 520
 521   switch(CurveType) {
 522   case GeomAbs_Line:
 523     {
 524       PerformConicSurf(TheCurveTool::Line(curve),curve,surface,U1,V1,U2,V2);
 525       break;
 526     }
 527   case GeomAbs_Circle:
 528     {
 529       PerformConicSurf(TheCurveTool::Circle(curve),curve,surface,U1,V1,U2,V2);
 530       break;
 531     }
 532   case GeomAbs_Ellipse:
 533     {
 534       PerformConicSurf(TheCurveTool::Ellipse(curve),curve,surface,U1,V1,U2,V2);
 535       break;
 536     }
 537   case GeomAbs_Parabola:
 538     {
 539       PerformConicSurf(TheCurveTool::Parabola(curve),curve,surface,U1,V1,U2,V2);
 540       break;
 541     }
 542   case GeomAbs_Hyperbola:
 543     {
 544       PerformConicSurf(TheCurveTool::Hyperbola(curve),curve,surface,U1,V1,U2,V2);
 545       break;
 546     }
 547   default:
 548     {
 549       Standard_Integer nbIntervalsOnCurve = TheCurveTool::NbIntervals(curve,GeomAbs_C2);
 550       GeomAbs_SurfaceType SurfaceType = TheSurfaceTool::GetType(surface);
 551       if(   (SurfaceType != GeomAbs_Plane)
 552          && (SurfaceType != GeomAbs_Cylinder)
 553          && (SurfaceType != GeomAbs_Cone)
 554          && (SurfaceType != GeomAbs_Sphere)) {
 555
 556         if(nbIntervalsOnCurve > 1) {
 557           TColStd_Array1OfReal TabW(1,nbIntervalsOnCurve+1);
 558           TheCurveTool::Intervals(curve,TabW,GeomAbs_C2);
 559           for(Standard_Integer i = 1; i<=nbIntervalsOnCurve; i++) {
 560             Standard_Real u1,u2;
 561             u1 = TabW.Value(i);
 562             u2 = TabW.Value(i+1);
 563
 564             Handle(TColStd_HArray1OfReal) aPars;
 565             Standard_Real defl = 0.1;
 566             Standard_Integer NbMin = 10;
 567             TheCurveTool::SamplePars(curve, u1, u2, defl, NbMin, aPars);
 568
 569 //          IntCurveSurface_ThePolygon  polygon(curve,u1,u2,TheCurveTool::NbSamples(curve,u1,u2));
 570             IntCurveSurface_ThePolygon  polygon(curve, aPars->Array1());
 571             InternalPerform(curve,polygon,surface,U1,V1,U2,V2);
 572           }
 573         }
 574         else {
 575           Standard_Real u1,u2;
 576           u1 = TheCurveTool::FirstParameter(curve);
 577           u2 = TheCurveTool::LastParameter(curve);
 578
 579           Handle(TColStd_HArray1OfReal) aPars;
 580           Standard_Real defl = 0.1;
 581           Standard_Integer NbMin = 10;
 582           TheCurveTool::SamplePars(curve, u1, u2, defl, NbMin, aPars);
 583
 584 //        IntCurveSurface_ThePolygon       polygon(curve,TheCurveTool::NbSamples(curve,u1,u2));
 585           IntCurveSurface_ThePolygon  polygon(curve, aPars->Array1());
 586           InternalPerform(curve,polygon,surface,U1,V1,U2,V2);
 587         }
 588       }
 589       else {   //-- la surface est une quadrique
 590         InternalPerformCurveQuadric(curve,surface);
 591       }
 592     }
 593   }
 594 }
 595 //=======================================================================
 596 //function : Perform
 597 //purpose  :
 598 //=======================================================================
 599 void IntCurveSurface_Inter::Perform(const TheCurve&   curve,
 600                                     const IntCurveSurface_ThePolygon& polygon,
 601                                     const TheSurface& surface) {
 602   ResetFields();
 603   done = Standard_True;
 604   Standard_Real u1,v1,u2,v2;
 605   u1 = TheSurfaceTool::FirstUParameter(surface);
 606   v1 = TheSurfaceTool::FirstVParameter(surface);
 607   u2 = TheSurfaceTool::LastUParameter(surface);
 608   v2 = TheSurfaceTool::LastVParameter(surface);
 609   Standard_Integer nbsu,nbsv;
 610   nbsu = TheSurfaceTool::NbSamplesU(surface,u1,u2);
 611   nbsv = TheSurfaceTool::NbSamplesV(surface,v1,v2);
 612   if(nbsu>40) nbsu=40;
 613   if(nbsv>40) nbsv=40;
 614   IntCurveSurface_ThePolyhedron    polyhedron(surface,nbsu,nbsv,u1,v1,u2,v2);
 615   Perform(curve,polygon,surface,polyhedron);
 616 }
 617 //=======================================================================
 618 //function : Perform
 619 //purpose  :
 620 //=======================================================================
 621 void IntCurveSurface_Inter::Perform(const TheCurve&   curve,
 622                                     const TheSurface& surface,
 623                                     const IntCurveSurface_ThePolyhedron& polyhedron) {
 624   ResetFields();
 625   done = Standard_True;
 626   Standard_Real u1 = TheCurveTool::FirstParameter(curve);
 627   Standard_Real u2 = TheCurveTool::LastParameter(curve);
 628   IntCurveSurface_ThePolygon polygon(curve,TheCurveTool::NbSamples(curve,u1,u2));
 629   Perform(curve,polygon,surface,polyhedron);
 630 }
 631 //=======================================================================
 632 //function : Perform
 633 //purpose  :
 634 //=======================================================================
 635 void IntCurveSurface_Inter::Perform(const TheCurve&       curve,
 636                                     const IntCurveSurface_ThePolygon&     polygon,
 637                                     const TheSurface&     surface,
 638                                     const IntCurveSurface_ThePolyhedron&  polyhedron) {
 639   ResetFields();
 640   done = Standard_True;
 641   Standard_Real u1,v1,u2,v2;
 642   u1 = TheSurfaceTool::FirstUParameter(surface);
 643   v1 = TheSurfaceTool::FirstVParameter(surface);
 644   u2 = TheSurfaceTool::LastUParameter(surface);
 645   v2 = TheSurfaceTool::LastVParameter(surface);
 646   InternalPerform(curve,polygon,surface,polyhedron,u1,v1,u2,v2);
 647 }
 648
 649 //=======================================================================
 650 //function : Perform
 651 //purpose  :
 652 //=======================================================================
 653 void IntCurveSurface_Inter::Perform(const TheCurve&       curve,
 654                                     const IntCurveSurface_ThePolygon&     polygon,
 655                                     const TheSurface&     surface,
 656                                     const IntCurveSurface_ThePolyhedron&  polyhedron,
 657                                     Bnd_BoundSortBox& BndBSB) {
 658   ResetFields();
 659   done = Standard_True;
 660   Standard_Real u1,v1,u2,v2;
 661   u1 = TheSurfaceTool::FirstUParameter(surface);
 662   v1 = TheSurfaceTool::FirstVParameter(surface);
 663   u2 = TheSurfaceTool::LastUParameter(surface);
 664   v2 = TheSurfaceTool::LastVParameter(surface);
 665   InternalPerform(curve,polygon,surface,polyhedron,u1,v1,u2,v2,BndBSB);
 666 }
 667 //=======================================================================
 668 //function : InternalPerform
 669 //purpose  : C a l c u l   d u   p o i n t   a p p r o c h e
 670 //==              p u i s   d u   p o i n t   E x a c t
 671 //=======================================================================
 672 void IntCurveSurface_Inter::InternalPerform(const TheCurve&       curve,
 673                                             const IntCurveSurface_ThePolygon&     polygon,
 674                                             const TheSurface&     surface,
 675                                             const IntCurveSurface_ThePolyhedron&  polyhedron,
 676                                             const Standard_Real u0,
 677                                             const Standard_Real v0,
 678                                             const Standard_Real u1,
 679                                             const Standard_Real v1,
 680                                             Bnd_BoundSortBox& BSB) {
 681   IntCurveSurface_TheInterference  interference(polygon,polyhedron,BSB);
 682   IntCurveSurface_TheCSFunction    theicsfunction(surface,curve);
 683   IntCurveSurface_TheExactInter    intersectionExacte(theicsfunction,TOLTANGENCY);
 684   math_FunctionSetRoot rsnld(intersectionExacte.Function());
 685
 686 //  Standard_Real    u,v,w,winit;
 687   Standard_Real    u,v,w;
 688   gp_Pnt P;
 689   Standard_Real winf = polygon.InfParameter();
 690   Standard_Real wsup = polygon.SupParameter();
 691   Standard_Integer NbSectionPoints = interference.NbSectionPoints();
 692   Standard_Integer NbTangentZones  = interference.NbTangentZones();
 693
 694
 695   //-- Les interferences renvoient parfois de nombreuses fois (>20) les memes points
 696
 697   Standard_Integer i,NbStartPoints=NbSectionPoints;
 698   for(i=1; i<= NbTangentZones; i++) {
 699     const Intf_TangentZone& TZ = interference.ZoneValue(i);
 700     Standard_Integer nbpnts = TZ.NumberOfPoints();
 701     NbStartPoints+=nbpnts;
 702   }
 703
 704   if(NbStartPoints) {
 705     Standard_Real *TabU = new Standard_Real [NbStartPoints+1];
 706     Standard_Real *TabV = new Standard_Real [NbStartPoints+1];
 707     Standard_Real *TabW = new Standard_Real [NbStartPoints+1];
 708     Standard_Integer IndexPoint=0;
 709
 710     for(i=1; i<= NbSectionPoints; i++) {
 711       const Intf_SectionPoint& SP = interference.PntValue(i);
 712       SectionPointToParameters(SP,polyhedron,polygon,u,v,w);
 713       TabU[IndexPoint]=u;
 714       TabV[IndexPoint]=v;
 715       TabW[IndexPoint]=w;
 716       IndexPoint++;
 717     }
 718     for(i=1; i<= NbTangentZones; i++) {
 719       const Intf_TangentZone& TZ = interference.ZoneValue(i);
 720       Standard_Integer nbpnts = TZ.NumberOfPoints();
 721       for(Standard_Integer j=1; j<=nbpnts; j++) {
 722         const Intf_SectionPoint& SP = TZ.GetPoint(j);
 723         SectionPointToParameters(SP,polyhedron,polygon,u,v,w);
 724         TabU[IndexPoint]=u;
 725         TabV[IndexPoint]=v;
 726         TabW[IndexPoint]=w;
 727         IndexPoint++;
 728       }
 729     }
 730
 731     //-- Tri
 732     Standard_Real su=0,sv=0,sw=0,ptol;
 733     ptol = 10*Precision::PConfusion();
 734
 735     //-- Tri suivant la variable W
 736     Standard_Boolean Triok;
 737     do {
 738       Triok=Standard_True;
 739       Standard_Integer im1;
 740       for(i=1,im1=0;i<NbStartPoints;im1++,i++) {
 741         if(TabW[i] < TabW[im1]) {
 742           Standard_Real t=TabW[i]; TabW[i]=TabW[im1]; TabW[im1]=t;
 743           t=TabU[i]; TabU[i]=TabU[im1]; TabU[im1]=t;
 744           t=TabV[i]; TabV[i]=TabV[im1]; TabV[im1]=t;
 745           Triok=Standard_False;
 746         }
 747       }
 748     }
 749     while(Triok==Standard_False);
 750
 751     //-- On trie pour des meme W suivant U
 752     do {
 753       Triok=Standard_True;
 754       Standard_Integer im1;
 755       for(i=1,im1=0;i<NbStartPoints;im1++,i++) {
 756 // modified by NIZHNY-MKK  Mon Oct  3 17:38:49 2005
 757 //      if(Abs(TabW[i]-TabW[im1])<ptol) {
 758         if((TabW[i]-TabW[im1])<ptol) {
 759           TabW[i]=TabW[im1];
 760           if(TabU[i] < TabU[im1]) {
 761             Standard_Real t=TabU[i]; TabU[i]=TabU[im1]; TabU[im1]=t;
 762             t=TabV[i]; TabV[i]=TabV[im1]; TabV[im1]=t;
 763             Triok=Standard_False;
 764           }
 765         }
 766       }
 767     }
 768     while(Triok==Standard_False);
 769
 770     //-- On trie pour des meme U et W suivant V
 771     do {
 772       Triok=Standard_True;
 773       Standard_Integer im1;
 774       for(i=1,im1=0;i<NbStartPoints;im1++,i++) {
 775 // modified by NIZHNY-MKK  Mon Oct  3 17:38:52 2005
 776 //      if((Abs(TabW[i]-TabW[im1])<ptol) && (Abs(TabU[i]-TabU[im1])<ptol)) {
 777         if(((TabW[i]-TabW[im1])<ptol) && ((TabU[i]-TabU[im1])<ptol)) {
 778           TabU[i]=TabU[im1];
 779           if(TabV[i] < TabV[im1]) {
 780             Standard_Real t=TabV[i]; TabV[i]=TabV[im1]; TabV[im1]=t;
 781             Triok=Standard_False;
 782           }
 783         }
 784       }
 785     }
 786     while(Triok==Standard_False);
 787
 788
 789     for(i=0;i<NbStartPoints; i++) {
 790       u=TabU[i]; v=TabV[i]; w=TabW[i];
 791       if(i==0) {
 792         su=u-1;
 793       }
 794       if(Abs(u-su)>ptol || Abs(v-sv)>ptol || Abs(w-sw)>ptol) {
 795         intersectionExacte.Perform(u,v,w,rsnld,u0,u1,v0,v1,winf,wsup);
 796         if(intersectionExacte.IsDone()) {
 797           if(!intersectionExacte.IsEmpty()) {
 798             P=intersectionExacte.Point();
 799             w=intersectionExacte.ParameterOnCurve();
 800             intersectionExacte.ParameterOnSurface(u,v);
 801             AppendPoint(curve,w,surface,u,v);
 802           }
 803         }
 804       }
 805       su=TabU[i]; sv=TabV[i]; sw=TabW[i];
 806     }
 807     delete [] TabW;
 808     delete [] TabV;
 809     delete [] TabU;
 810   }
 811 }
 812
 813 //=======================================================================
 814 //function : InternalPerform
 815 //purpose  :
 816 //=======================================================================
 817 void IntCurveSurface_Inter::InternalPerform(const TheCurve&       curve,
 818                                             const IntCurveSurface_ThePolygon&     polygon,
 819                                             const TheSurface&     surface,
 820                                             const IntCurveSurface_ThePolyhedron&  polyhedron,
 821                                             const Standard_Real u0,
 822                                             const Standard_Real v0,
 823                                             const Standard_Real u1,
 824                                             const Standard_Real v1) {
 825   IntCurveSurface_TheInterference  interference(polygon,polyhedron);
 826   IntCurveSurface_TheCSFunction    theicsfunction(surface,curve);
 827   IntCurveSurface_TheExactInter    intersectionExacte(theicsfunction,TOLTANGENCY);
 828   math_FunctionSetRoot rsnld(intersectionExacte.Function());
 829
 830 //  Standard_Real    u,v,w,winit;
 831   Standard_Real    u,v,w;
 832   gp_Pnt P;
 833   Standard_Real winf = polygon.InfParameter();
 834   Standard_Real wsup = polygon.SupParameter();
 835   Standard_Integer NbSectionPoints = interference.NbSectionPoints();
 836   Standard_Integer NbTangentZones  = interference.NbTangentZones();
 837
 838
 839   //-- Les interferences renvoient parfois de nombreuses fois (>20) les memes points
 840
 841   Standard_Integer i,NbStartPoints=NbSectionPoints;
 842   for(i=1; i<= NbTangentZones; i++) {
 843     const Intf_TangentZone& TZ = interference.ZoneValue(i);
 844     Standard_Integer nbpnts = TZ.NumberOfPoints();
 845     NbStartPoints+=nbpnts;
 846   }
 847
 848   if(NbStartPoints) {
 849     Standard_Real *TabU = new Standard_Real [NbStartPoints+1];
 850     Standard_Real *TabV = new Standard_Real [NbStartPoints+1];
 851     Standard_Real *TabW = new Standard_Real [NbStartPoints+1];
 852     Standard_Integer IndexPoint=0;
 853
 854     for(i=1; i<= NbSectionPoints; i++) {
 855       const Intf_SectionPoint& SP = interference.PntValue(i);
 856       SectionPointToParameters(SP,polyhedron,polygon,u,v,w);
 857       TabU[IndexPoint]=u;
 858       TabV[IndexPoint]=v;
 859       TabW[IndexPoint]=w;
 860       IndexPoint++;
 861     }
 862     for(i=1; i<= NbTangentZones; i++) {
 863       const Intf_TangentZone& TZ = interference.ZoneValue(i);
 864       Standard_Integer nbpnts = TZ.NumberOfPoints();
 865       for(Standard_Integer j=1; j<=nbpnts; j++) {
 866         const Intf_SectionPoint& SP = TZ.GetPoint(j);
 867         SectionPointToParameters(SP,polyhedron,polygon,u,v,w);
 868         TabU[IndexPoint]=u;
 869         TabV[IndexPoint]=v;
 870         TabW[IndexPoint]=w;
 871         IndexPoint++;
 872       }
 873     }
 874
 875     //-- Tri
 876     Standard_Real su=0,sv=0,sw=0,ptol;
 877     ptol = 10*Precision::PConfusion();
 878
 879     //-- Tri suivant la variable W
 880     Standard_Boolean Triok;
 881     do {
 882       Triok=Standard_True;
 883       Standard_Integer im1;
 884       for(i=1,im1=0;i<NbStartPoints;im1++,i++) {
 885         if(TabW[i] < TabW[im1]) {
 886           Standard_Real t=TabW[i]; TabW[i]=TabW[im1]; TabW[im1]=t;
 887           t=TabU[i]; TabU[i]=TabU[im1]; TabU[im1]=t;
 888           t=TabV[i]; TabV[i]=TabV[im1]; TabV[im1]=t;
 889           Triok=Standard_False;
 890         }
 891       }
 892     }
 893     while(Triok==Standard_False);
 894
 895     //-- On trie pour des meme W suivant U
 896     do {
 897       Triok=Standard_True;
 898       Standard_Integer im1;
 899       for(i=1,im1=0;i<NbStartPoints;im1++,i++) {
 900 // modified by NIZHNY-MKK  Mon Oct  3 17:38:56 2005
 901 //      if(Abs(TabW[i]-TabW[im1])<ptol) {
 902         if((TabW[i]-TabW[im1])<ptol) {
 903           TabW[i]=TabW[im1];
 904           if(TabU[i] < TabU[im1]) {
 905             Standard_Real t=TabU[i]; TabU[i]=TabU[im1]; TabU[im1]=t;
 906             t=TabV[i]; TabV[i]=TabV[im1]; TabV[im1]=t;
 907             Triok=Standard_False;
 908           }
 909         }
 910       }
 911     }
 912     while(Triok==Standard_False);
 913
 914     //-- On trie pour des meme U et W suivant V
 915     do {
 916       Triok=Standard_True;
 917       Standard_Integer im1;
 918       for(i=1,im1=0;i<NbStartPoints;im1++,i++) {
 919 // modified by NIZHNY-MKK  Mon Oct  3 17:38:58 2005
 920 //      if((Abs(TabW[i]-TabW[im1])<ptol) && (Abs(TabU[i]-TabU[im1])<ptol)) {
 921         if(((TabW[i]-TabW[im1])<ptol) && ((TabU[i]-TabU[im1])<ptol)) {
 922           TabU[i]=TabU[im1];
 923           if(TabV[i] < TabV[im1]) {
 924             Standard_Real t=TabV[i]; TabV[i]=TabV[im1]; TabV[im1]=t;
 925             Triok=Standard_False;
 926           }
 927         }
 928       }
 929     }
 930     while(Triok==Standard_False);
 931
 932
 933     for(i=0;i<NbStartPoints; i++) {
 934       u=TabU[i]; v=TabV[i]; w=TabW[i];
 935       if(i==0) {
 936         su=u-1;
 937       }
 938       if(Abs(u-su)>ptol || Abs(v-sv)>ptol || Abs(w-sw)>ptol) {
 939         intersectionExacte.Perform(u,v,w,rsnld,u0,u1,v0,v1,winf,wsup);
 940         if(intersectionExacte.IsDone()) {
 941           if(!intersectionExacte.IsEmpty()) {
 942             P=intersectionExacte.Point();
 943             w=intersectionExacte.ParameterOnCurve();
 944             intersectionExacte.ParameterOnSurface(u,v);
 945             AppendPoint(curve,w,surface,u,v);
 946           }
 947         }
 948       }
 949       su=TabU[i]; sv=TabV[i]; sw=TabW[i];
 950     }
 951     delete [] TabW;
 952     delete [] TabV;
 953     delete [] TabU;
 954   }
 955 }
 956 //=======================================================================
 957 //function : InternalPerformCurveQuadric
 958 //purpose  :
 959 //=======================================================================
 960 void IntCurveSurface_Inter::InternalPerformCurveQuadric(const TheCurve&                 curve,
 961                                                         const TheSurface&               surface) {
 962   IntCurveSurface_TheQuadCurvExactInter QuadCurv(surface,curve);
 963   if(QuadCurv.IsDone()) {
 964     Standard_Integer NbRoots = QuadCurv.NbRoots();
 965     Standard_Real u,v,w;
 966     for(Standard_Integer i = 1; i<= NbRoots; i++) {
 967       w = QuadCurv.Root(i);
 968       IntCurveSurface_ComputeParamsOnQuadric(surface,TheCurveTool::Value(curve,w),u,v);
 969       AppendPoint(curve,w,surface,u,v);
 970     }
 971     //-- Intervals non traites .............................................
 972   }
 973 }
 974
 975 //=======================================================================
 976 //function : InternalPerform
 977 //purpose  :
 978 //=======================================================================
 979 void IntCurveSurface_Inter::InternalPerform(const TheCurve&       curve,
 980                                             const IntCurveSurface_ThePolygon&  polygon,
 981                                             const TheSurface&     surface,
 982                                             const Standard_Real U1,
 983                                             const Standard_Real V1,
 984                                             const Standard_Real U2,
 985                                             const Standard_Real V2) {
 986   GeomAbs_SurfaceType SurfaceType = TheSurfaceTool::GetType(surface);
 987   if(   (SurfaceType != GeomAbs_Plane)
 988      && (SurfaceType != GeomAbs_Cylinder)
 989      && (SurfaceType != GeomAbs_Cone)
 990      && (SurfaceType != GeomAbs_Sphere) ) {
 991     if(SurfaceType != GeomAbs_BSplineSurface) {
 992       Standard_Integer nbsu,nbsv;
 993       nbsu = TheSurfaceTool::NbSamplesU(surface,U1,U2);
 994       nbsv = TheSurfaceTool::NbSamplesV(surface,V1,V2);
 995       if(nbsu>40) nbsu=40;
 996       if(nbsv>40) nbsv=40;
 997       IntCurveSurface_ThePolyhedron polyhedron(surface,nbsu,nbsv,U1,V1,U2,V2);
 998       InternalPerform(curve,polygon,surface,polyhedron,U1,V1,U2,V2);
 999     }
1000     else {
1001       Handle(Adaptor3d_HSurface) aS = TheSurfaceTool::UTrim(surface, U1, U2, 1.e-9);
1002       aS = aS->VTrim(V1, V2, 1.e-9);
1003       Handle(Adaptor3d_TopolTool) aTopTool = new Adaptor3d_TopolTool(aS);
1004       Standard_Real defl = 0.1;
1005       aTopTool->SamplePnts(defl, 10, 10);
1006
1007       Standard_Integer nbpu = aTopTool->NbSamplesU();
1008       Standard_Integer nbpv = aTopTool->NbSamplesV();
1009       TColStd_Array1OfReal Upars(1, nbpu), Vpars(1, nbpv);
1010       aTopTool->UParameters(Upars);
1011       aTopTool->VParameters(Vpars);
1012
1013       IntCurveSurface_ThePolyhedron polyhedron(surface,Upars, Vpars);
1014       InternalPerform(curve,polygon,surface,polyhedron,U1,V1,U2,V2);
1015     }
1016   }
1017   else {
1018     IntCurveSurface_TheQuadCurvExactInter QuadCurv(surface,curve);
1019     if(QuadCurv.IsDone()) {
1020       Standard_Integer NbRoots = QuadCurv.NbRoots();
1021       Standard_Real u,v,w;
1022       for(Standard_Integer i = 1; i<= NbRoots; i++) {
1023         w = QuadCurv.Root(i);
1024         IntCurveSurface_ComputeParamsOnQuadric(surface,TheCurveTool::Value(curve,w),u,v);
1025         AppendPoint(curve,w,surface,u,v);
1026       }
1027       //-- Intervalles non traites .............................................
1028     }
1029   } //-- Fin : la Surface  est une quadrique
1030 }
1031 //=======================================================================
1032 //function : PerformConicSurf
1033 //purpose  :
1034 //=======================================================================
1035 void IntCurveSurface_Inter::PerformConicSurf(const gp_Lin&   Line,
1036                                              const TheCurve&       curve,
1037                                              const TheSurface&     surface,
1038                                              const Standard_Real U1,
1039                                              const Standard_Real V1,
1040                                              const Standard_Real U2,
1041                                              const Standard_Real V2) {
1042
1043
1044   GeomAbs_SurfaceType SurfaceType = TheSurfaceTool::GetType(surface);
1045   switch(SurfaceType) {
1046   case GeomAbs_Plane:
1047     {
1048       IntAna_IntConicQuad LinPlane(Line,TheSurfaceTool::Plane(surface),TOLERANCE_ANGULAIRE);
1049       AppendIntAna(curve,surface,LinPlane);
1050       break;
1051     }
1052   case GeomAbs_Cylinder:
1053     {
1054       IntAna_IntConicQuad LinCylinder(Line,TheSurfaceTool::Cylinder(surface));
1055       AppendIntAna(curve,surface,LinCylinder);
1056       break;
1057     }
1058   case GeomAbs_Sphere:
1059     {
1060       IntAna_IntConicQuad LinSphere(Line,TheSurfaceTool::Sphere(surface));
1061       AppendIntAna(curve,surface,LinSphere);
1062       break;
1063     }
1064   case GeomAbs_Torus:
1065     {
1066       IntAna_IntLinTorus intlintorus(Line,TheSurfaceTool::Torus(surface));
1067       if(intlintorus.IsDone()) {
1068         Standard_Integer nbp = intlintorus.NbPoints();
1069         Standard_Real fi,theta,w;
1070         for(Standard_Integer i = 1; i<= nbp; i++) {
1071 #ifndef DEB
1072           gp_Pnt P;
1073           P = intlintorus.Value(i);
1074 #else
1075           gp_Pnt P(intlintorus.Value(i));
1076 #endif
1077           w = intlintorus.ParamOnLine(i);
1078           intlintorus.ParamOnTorus(i,fi,theta);
1079           AppendPoint(curve,w,surface,fi,theta);
1080         }
1081         break;
1082       }
1083     }  //-- Si Done retourne False, On passe dans Default !!
1084   case GeomAbs_Cone:
1085     {
1086       //OCC516(apo)->
1087       const Standard_Real correction = 1.E+5*Precision::Angular();
1088       gp_Cone cn = TheSurfaceTool::Cone(surface);
1089       if(Abs(cn.SemiAngle()) < M_PI/2.0 - correction){
1090         IntAna_IntConicQuad LinCone(Line,cn);
1091         AppendIntAna(curve,surface,LinCone);
1092         break;
1093       }//<-OCC516(apo)
1094     }
1095   default:
1096     {
1097       Standard_Integer nbsu,nbsv;
1098       nbsu = TheSurfaceTool::NbSamplesU(surface,U1,U2);
1099       nbsv = TheSurfaceTool::NbSamplesV(surface,V1,V2);
1100
1101       Standard_Boolean U1inf = Precision::IsInfinite(U1);
1102       Standard_Boolean U2inf = Precision::IsInfinite(U2);
1103       Standard_Boolean V1inf = Precision::IsInfinite(V1);
1104       Standard_Boolean V2inf = Precision::IsInfinite(V2);
1105
1106       Standard_Real U1new=U1, U2new=U2, V1new=V1, V2new=V2;
1107
1108       Standard_Boolean NoIntersection = Standard_False;
1109
1110       if(U1inf || U2inf || V1inf || V2inf ) {
1111
1112         if(SurfaceType == GeomAbs_SurfaceOfExtrusion) {
1113
1114           EstLimForInfExtr(Line, surface, Standard_False, nbsu,
1115                            U1inf, U2inf, V1inf, V2inf,
1116                            U1new, U2new, V1new, V2new, NoIntersection);
1117
1118         }
1119         else if(SurfaceType == GeomAbs_SurfaceOfRevolution) {
1120
1121           EstLimForInfRevl(Line, surface,
1122                            U1inf, U2inf, V1inf, V2inf,
1123                            U1new, U2new, V1new, V2new, NoIntersection);
1124
1125         }
1126         else if(SurfaceType == GeomAbs_OffsetSurface) {
1127
1128           EstLimForInfOffs(Line, surface, nbsu,
1129                            U1inf, U2inf, V1inf, V2inf,
1130                            U1new, U2new, V1new, V2new, NoIntersection);
1131         }
1132         else {
1133
1134           EstLimForInfSurf(U1new, U2new, V1new, V2new);
1135         }
1136
1137       }
1138
1139
1140       if(NoIntersection) return;
1141
1142
1143       // modified by NIZHNY-OFV  Mon Aug 20 14:56:47 2001 (60963 begin)
1144       if(nbsu<20) nbsu=20;
1145       if(nbsv<20) nbsv=20;
1146       // modified by NIZHNY-OFV  Mon Aug 20 14:57:06 2001 (60963 end)
1147       IntCurveSurface_ThePolyhedron polyhedron(surface,nbsu,nbsv,U1new,V1new,U2new,V2new);
1148       Intf_Tool                         bndTool;
1149       Bnd_Box                          boxLine;
1150       bndTool.LinBox(Line,polyhedron.Bounding(),boxLine);
1151       for(Standard_Integer nbseg=1; nbseg<= bndTool.NbSegments(); nbseg++) {
1152         Standard_Real pinf = bndTool.BeginParam(nbseg);
1153         Standard_Real psup = bndTool.EndParam(nbseg);
1154         if((psup - pinf)<1e-10) { pinf-=1e-10; psup+=1e-10; }
1155         IntCurveSurface_ThePolygon polygon(curve, pinf,psup,2);
1156         InternalPerform(curve,polygon,surface,polyhedron,U1new,V1new,U2new,V2new);
1157       }
1158     }
1159   }
1160 }
1161 //=======================================================================
1162 //function : PerformConicSurf
1163 //purpose  :
1164 //=======================================================================
1165 void IntCurveSurface_Inter::PerformConicSurf(const gp_Circ&        Circle,
1166                                              const TheCurve&       curve,
1167                                              const TheSurface&     surface,
1168                                              const Standard_Real U1,
1169                                              const Standard_Real V1,
1170                                              const Standard_Real U2,
1171                                              const Standard_Real V2) {
1172
1173   GeomAbs_SurfaceType SurfaceType = TheSurfaceTool::GetType(surface);
1174   switch(SurfaceType) {
1175   case GeomAbs_Plane:
1176     {
1177       IntAna_IntConicQuad CircPlane(Circle,TheSurfaceTool::Plane(surface),TOLERANCE_ANGULAIRE,TOLERANCE);
1178       AppendIntAna(curve,surface,CircPlane);
1179       break;
1180     }
1181   case GeomAbs_Cylinder:
1182     {
1183       IntAna_IntConicQuad CircCylinder(Circle,TheSurfaceTool::Cylinder(surface));
1184       AppendIntAna(curve,surface,CircCylinder);
1185       break;
1186     }
1187   case GeomAbs_Cone:
1188     {
1189       IntAna_IntConicQuad CircCone(Circle,TheSurfaceTool::Cone(surface));
1190       AppendIntAna(curve,surface,CircCone);
1191       break;
1192     }
1193   case GeomAbs_Sphere:
1194     {
1195       IntAna_IntConicQuad CircSphere(Circle,TheSurfaceTool::Sphere(surface));
1196       AppendIntAna(curve,surface,CircSphere);
1197       break;
1198     }
1199   default:
1200     {
1201       IntCurveSurface_ThePolygon polygon(curve,NBSAMPLESONCIRCLE);
1202       InternalPerform(curve,polygon,surface,U1,V1,U2,V2);
1203     }
1204   }
1205 }
1206 //=======================================================================
1207 //function : PerformConicSurf
1208 //purpose  :
1209 //=======================================================================
1210 void IntCurveSurface_Inter::PerformConicSurf(const gp_Elips&       Ellipse,
1211                                              const TheCurve&       curve,
1212                                              const TheSurface&     surface,
1213                                              const Standard_Real U1,
1214                                              const Standard_Real V1,
1215                                              const Standard_Real U2,
1216                                              const Standard_Real V2) {
1217
1218   GeomAbs_SurfaceType SurfaceType = TheSurfaceTool::GetType(surface);
1219   switch(SurfaceType) {
1220   case GeomAbs_Plane:
1221     {
1222       IntAna_IntConicQuad EllipsePlane(Ellipse,TheSurfaceTool::Plane(surface),TOLERANCE_ANGULAIRE,TOLERANCE);
1223       AppendIntAna(curve,surface,EllipsePlane);
1224       break;
1225     }
1226   case GeomAbs_Cylinder:
1227     {
1228       IntAna_IntConicQuad EllipseCylinder(Ellipse,TheSurfaceTool::Cylinder(surface));
1229       AppendIntAna(curve,surface,EllipseCylinder);
1230       break;
1231     }
1232   case GeomAbs_Cone:
1233      {
1234       IntAna_IntConicQuad EllipseCone(Ellipse,TheSurfaceTool::Cone(surface));
1235       AppendIntAna(curve,surface,EllipseCone);
1236       break;
1237     }
1238   case GeomAbs_Sphere:
1239     {
1240       IntAna_IntConicQuad EllipseSphere(Ellipse,TheSurfaceTool::Sphere(surface));
1241       AppendIntAna(curve,surface,EllipseSphere);
1242       break;
1243     }
1244   default:
1245     {
1246       IntCurveSurface_ThePolygon polygon(curve,NBSAMPLESONELLIPSE);
1247       InternalPerform(curve,polygon,surface,U1,V1,U2,V2);
1248     }
1249   }
1250 }
1251 //=======================================================================
1252 //function : PerformConicSurf
1253 //purpose  :
1254 //=======================================================================
1255 void IntCurveSurface_Inter::PerformConicSurf(const gp_Parab&       Parab,
1256                                              const TheCurve&       curve,
1257                                              const TheSurface&     surface,
1258                                              const Standard_Real U1,
1259                                              const Standard_Real V1,
1260                                              const Standard_Real U2,
1261                                              const Standard_Real V2) {
1262
1263   GeomAbs_SurfaceType SurfaceType = TheSurfaceTool::GetType(surface);
1264   switch(SurfaceType) {
1265   case GeomAbs_Plane:
1266     {
1267       IntAna_IntConicQuad ParabPlane(Parab,TheSurfaceTool::Plane(surface),TOLERANCE_ANGULAIRE);
1268       AppendIntAna(curve,surface,ParabPlane);
1269       break;
1270     }
1271   case GeomAbs_Cylinder:
1272     {
1273       IntAna_IntConicQuad ParabCylinder(Parab,TheSurfaceTool::Cylinder(surface));
1274       AppendIntAna(curve,surface,ParabCylinder);
1275       break;
1276     }
1277   case GeomAbs_Cone:
1278      {
1279       IntAna_IntConicQuad ParabCone(Parab,TheSurfaceTool::Cone(surface));
1280       AppendIntAna(curve,surface,ParabCone);
1281       break;
1282     }
1283   case GeomAbs_Sphere:
1284     {
1285       IntAna_IntConicQuad ParabSphere(Parab,TheSurfaceTool::Sphere(surface));
1286       AppendIntAna(curve,surface,ParabSphere);
1287       break;
1288     }
1289   default:
1290     {
1291       Standard_Integer nbsu,nbsv;
1292       nbsu = TheSurfaceTool::NbSamplesU(surface,U1,U2);
1293       nbsv = TheSurfaceTool::NbSamplesV(surface,V1,V2);
1294       if(nbsu>40) nbsu=40;
1295       if(nbsv>40) nbsv=40;
1296       IntCurveSurface_ThePolyhedron polyhedron(surface,nbsu,nbsv,U1,V1,U2,V2);
1297       Intf_Tool                         bndTool;
1298       Bnd_Box                          boxParab;
1299       bndTool.ParabBox(Parab,polyhedron.Bounding(),boxParab);
1300       for(Standard_Integer nbseg=1; nbseg<= bndTool.NbSegments(); nbseg++) {
1301         IntCurveSurface_ThePolygon polygon(curve,
1302                                            bndTool.BeginParam(nbseg),
1303                                            bndTool.EndParam(nbseg),
1304                                            NBSAMPLESONPARAB);
1305         InternalPerform(curve,polygon,surface,polyhedron,U1,V1,U2,V2);
1306       }
1307     }
1308   }
1309 }
1310 //=======================================================================
1311 //function : PerformConicSurf
1312 //purpose  :
1313 //=======================================================================
1314 void IntCurveSurface_Inter::PerformConicSurf(const gp_Hypr&        Hypr,
1315                                              const TheCurve&       curve,
1316                                              const TheSurface&     surface,
1317                                              const Standard_Real U1,
1318                                              const Standard_Real V1,
1319                                              const Standard_Real U2,
1320                                              const Standard_Real V2) {
1321
1322   GeomAbs_SurfaceType SurfaceType = TheSurfaceTool::GetType(surface);
1323   switch(SurfaceType) {
1324   case GeomAbs_Plane:
1325     {
1326       IntAna_IntConicQuad HyprPlane(Hypr,TheSurfaceTool::Plane(surface),TOLERANCE_ANGULAIRE);
1327       AppendIntAna(curve,surface,HyprPlane);
1328       break;
1329     }
1330   case GeomAbs_Cylinder:
1331     {
1332       IntAna_IntConicQuad HyprCylinder(Hypr,TheSurfaceTool::Cylinder(surface));
1333       AppendIntAna(curve,surface,HyprCylinder);
1334       break;
1335     }
1336   case GeomAbs_Cone:
1337      {
1338       IntAna_IntConicQuad HyprCone(Hypr,TheSurfaceTool::Cone(surface));
1339       AppendIntAna(curve,surface,HyprCone);
1340       break;
1341     }
1342   case GeomAbs_Sphere:
1343     {
1344       IntAna_IntConicQuad HyprSphere(Hypr,TheSurfaceTool::Sphere(surface));
1345       AppendIntAna(curve,surface,HyprSphere);
1346       break;
1347     }
1348   default:
1349     {
1350       Standard_Integer nbsu,nbsv;
1351       nbsu = TheSurfaceTool::NbSamplesU(surface,U1,U2);
1352       nbsv = TheSurfaceTool::NbSamplesV(surface,V1,V2);
1353       if(nbsu>40) nbsu=40;
1354       if(nbsv>40) nbsv=40;
1355       IntCurveSurface_ThePolyhedron polyhedron(surface,nbsu,nbsv,U1,V1,U2,V2);
1356       Intf_Tool                         bndTool;
1357       Bnd_Box                          boxHypr;
1358       bndTool.HyprBox(Hypr,polyhedron.Bounding(),boxHypr);
1359       for(Standard_Integer nbseg=1; nbseg<= bndTool.NbSegments(); nbseg++) {
1360         IntCurveSurface_ThePolygon polygon(curve,
1361                                            bndTool.BeginParam(nbseg),
1362                                            bndTool.EndParam(nbseg),
1363                                            NBSAMPLESONHYPR);
1364         InternalPerform(curve,polygon,surface,polyhedron,U1,V1,U2,V2);
1365       }
1366     }
1367   }
1368 }
1369 //=======================================================================
1370 //function : AppendIntAna
1371 //purpose  :
1372 //=======================================================================
1373 void IntCurveSurface_Inter::AppendIntAna(const TheCurve& curve,
1374                                          const TheSurface& surface,
1375                                          const IntAna_IntConicQuad& intana_ConicQuad) {
1376   if(intana_ConicQuad.IsDone()) {
1377     if(intana_ConicQuad.IsInQuadric()) {
1378       //-- cout<<" Courbe Dans la Quadrique !!! Non Traite !!!"<<endl;
1379     }
1380     else if(intana_ConicQuad.IsParallel()) {
1381       //-- cout<<" Courbe // a la Quadrique !!! Non Traite !!!"<<endl;
1382     }
1383     else {
1384       Standard_Integer nbp = intana_ConicQuad.NbPoints();
1385       Standard_Real u,v,w;
1386       for(Standard_Integer i = 1; i<= nbp; i++) {
1387         gp_Pnt P(intana_ConicQuad.Point(i));
1388         w = intana_ConicQuad.ParamOnConic(i);
1389         IntCurveSurface_ComputeParamsOnQuadric(surface,P,u,v);
1390         AppendPoint(curve,w,surface,u,v);
1391       }
1392     }
1393   }
1394   else {
1395     //-- cout<<" IntAna Conic Quad Not Done  "<<endl;
1396   }
1397 }
1398 //=======================================================================
1399 //function : AppendPoint
1400 //purpose  :
1401 //=======================================================================
1402 void IntCurveSurface_Inter::AppendPoint(const TheCurve& curve,
1403                                         const Standard_Real lw,
1404                                         const TheSurface& surface,
1405                                         const Standard_Real su,
1406                                         const Standard_Real sv) {
1407
1408   Standard_Real W0 = TheCurveTool::FirstParameter(curve);
1409   Standard_Real W1 = TheCurveTool::LastParameter(curve);
1410   Standard_Real U0 = TheSurfaceTool::FirstUParameter(surface);
1411   Standard_Real U1 = TheSurfaceTool::LastUParameter(surface);
1412   Standard_Real V0 = TheSurfaceTool::FirstVParameter(surface);
1413   Standard_Real V1 = TheSurfaceTool::LastVParameter(surface);
1414   //-- Test si la courbe est periodique
1415   Standard_Real w = lw, u = su, v = sv;
1416
1417   GeomAbs_CurveType aCType = TheCurveTool::GetType(curve);
1418
1419   if(TheCurveTool::IsPeriodic(curve)
1420      || aCType == GeomAbs_Circle
1421      || aCType == GeomAbs_Ellipse) {
1422     w = ElCLib::InPeriod(w, W0, W0 + TheCurveTool::Period(curve));
1423   }
1424
1425   if((W0 - w) >= TOLTANGENCY || (w - W1) >= TOLTANGENCY) return;
1426
1427   GeomAbs_SurfaceType aSType = TheSurfaceTool::GetType(surface);
1428   if (TheSurfaceTool::IsUPeriodic(surface)
1429       || aSType == GeomAbs_Cylinder
1430       || aSType == GeomAbs_Cone
1431       || aSType == GeomAbs_Sphere) {
1432     u = ElCLib::InPeriod(u, U0, U0 + TheSurfaceTool::UPeriod(surface));
1433   }
1434
1435   if (TheSurfaceTool::IsVPeriodic(surface)) {
1436     v = ElCLib::InPeriod(v, V0, V0 + TheSurfaceTool::VPeriod(surface));
1437   }
1438
1439   if((U0 - u) >= TOLTANGENCY || (u - U1) >= TOLTANGENCY) return;
1440   if((V0 - v) >= TOLTANGENCY || (v - V1) >= TOLTANGENCY) return;
1441
1442   IntCurveSurface_TransitionOnCurve   TransOnCurve;
1443   IntCurveSurface_ComputeTransitions(curve,w,TransOnCurve,
1444                                      surface,u,v);
1445   gp_Pnt P(TheCurveTool::Value(curve,w));
1446   IntCurveSurface_IntersectionPoint IP(P,u,v,w,TransOnCurve);
1447   Append(IP); //-- invoque la methode de IntCurveSurface_Intersection.
1448
1449 }
1450
1451 //=======================================================================
1452 //function : AppendSegment
1453 //purpose  :
1454 //=======================================================================
1455 void IntCurveSurface_Inter::AppendSegment(const TheCurve& ,
1456                                           const Standard_Real ,
1457                                           const Standard_Real ,
1458                                           const TheSurface& ) {
1459   //cout<<" !!! Not Yet Implemented
1460   //IntCurveSurface_Inter::Append(const IntCurveSurf ...)"<<endl;
1461 }
1462
1463 //=======================================================================
1464 //function : SectionPointToParameters
1465 //purpose  : P o i n t   d   i n t e r f e r e n c e   - - >
1466 //    U , V       e t   W
1467 //=======================================================================
1468 void SectionPointToParameters(const Intf_SectionPoint&              Sp,
1469                               const IntCurveSurface_ThePolyhedron&  Polyhedron,
1470                               const IntCurveSurface_ThePolygon&     Polygon,
1471                               Standard_Real& U,
1472                               Standard_Real& V,
1473                               Standard_Real& W) {
1474
1475   Intf_PIType       typ;
1476   Standard_Integer  Adr1,Adr2;
1477   Standard_Real     Param,u,v;
1478   gp_Pnt P(Sp.Pnt());
1479
1480   Standard_Integer Pt1,Pt2,Pt3;
1481   Standard_Real u1 = 0.,v1 = 0.,param;
1482   //----------------------------------------------------------------------
1483   //--          Calcul des parametres approches sur la surface          --
1484   //----------------------------------------------------------------------
1485
1486   Sp.InfoSecond(typ,Adr1,Adr2,Param);
1487   switch(typ) {
1488   case Intf_VERTEX:   //-- Adr1 est le numero du vertex
1489     {
1490       Polyhedron.Parameters(Adr1,u1,v1);
1491       break;
1492     }
1493   case Intf_EDGE:
1494     {
1495       Polyhedron.Parameters(Adr1,u1,v1);
1496       Polyhedron.Parameters(Adr2,u,v);
1497       u1+= Param * (u-u1);
1498       v1+= Param * (v-v1);
1499       break;
1500     }
1501   case Intf_FACE:
1502     {
1503       Standard_Real ua,va,ub,vb,uc,vc,ca,cb,cc,cabc;
1504       Polyhedron.Triangle(Adr1,Pt1,Pt2,Pt3);
1505       gp_Pnt PA(Polyhedron.Point(Pt1));
1506       gp_Pnt PB(Polyhedron.Point(Pt2));
1507       gp_Pnt PC(Polyhedron.Point(Pt3));
1508       Polyhedron.Parameters(Pt1,ua,va);
1509       Polyhedron.Parameters(Pt2,ub,vb);
1510       Polyhedron.Parameters(Pt3,uc,vc);
1511       gp_Vec Normale(gp_Vec(PA,PB).Crossed(gp_Vec(PA,PC)));
1512       cc = (gp_Vec(PA,PB).Crossed(gp_Vec(PA,P))).Dot(Normale);
1513       ca = (gp_Vec(PB,PC).Crossed(gp_Vec(PB,P))).Dot(Normale);
1514       cb = (gp_Vec(PC,PA).Crossed(gp_Vec(PC,P))).Dot(Normale);
1515       cabc = ca + cb + cc;
1516
1517       ca/=cabc;     cb/=cabc;       cc/=cabc;
1518
1519       u1 = ca * ua + cb * ub + cc * uc;
1520       v1 = ca * va + cb * vb + cc * vc;
1521       break;
1522     }
1523   default:
1524     {
1525       cout<<" Default dans SectionPointToParameters "<<endl;
1526       break;
1527     }
1528   }
1529   //----------------------------------------------------------------------
1530   //--                Calcul du point approche sur la Curve             --
1531   //----------------------------------------------------------------------
1532   Standard_Integer SegIndex;
1533
1534   Sp.InfoFirst(typ,SegIndex,param);
1535   W = Polygon.ApproxParamOnCurve(SegIndex,param);
1536   if(param>1.0 || param<0.0) {
1537     //-- IntCurveSurface_ThePolyhedronTool::Dump(Polyhedron);
1538     //-- IntCurveSurface_ThePolygonTool::Dump(Polygon);
1539     }
1540   U = u1;
1541   V = v1;
1542 }
1543 //=======================================================================
1544 //function : IntCurveSurface_ComputeTransitions
1545 //purpose  :
1546 //=======================================================================
1547 void IntCurveSurface_ComputeTransitions(const TheCurve& curve,
1548                                         const Standard_Real w,
1549                                         IntCurveSurface_TransitionOnCurve&   TransOnCurve,
1550                                         const TheSurface& surface,
1551                                         const Standard_Real u,
1552                                         const Standard_Real v) {
1553
1554   gp_Vec NSurf,D1U,D1V;//TgCurv;
1555   gp_Pnt Psurf;
1556   Standard_Real CosDir;
1557
1558
1559   TheSurfaceTool::D1(surface,u,v,Psurf,D1U,D1V);
1560   NSurf = D1U.Crossed(D1V);
1561   TheCurveTool::D1(curve,w,Psurf,D1U);
1562   Standard_Real Norm = NSurf.Magnitude();
1563   if(Norm>TOLERANCE_ANGULAIRE) {
1564     D1U.Normalize();
1565     CosDir = NSurf.Dot(D1U);
1566     CosDir/=Norm;
1567     if( -CosDir > TOLERANCE_ANGULAIRE) {
1568       //--  --Curve--->    <----Surface----
1569       TransOnCurve   = IntCurveSurface_In;
1570     }
1571     else if(CosDir > TOLERANCE_ANGULAIRE) {
1572       //--  --Curve--->  ----Surface-->
1573       TransOnCurve   = IntCurveSurface_Out;
1574     }
1575     else {
1576       TransOnCurve   = IntCurveSurface_Tangent;
1577     }
1578   }
1579   else {
1580     TransOnCurve   = IntCurveSurface_Tangent;
1581   }
1582 }
1583 //=======================================================================
1584 //function : IntCurveSurface_ComputeParamsOnQuadric
1585 //purpose  :
1586 //=======================================================================
1587 void IntCurveSurface_ComputeParamsOnQuadric(const TheSurface& surface,
1588                                             const gp_Pnt& P,
1589                                             Standard_Real& u,
1590                                             Standard_Real& v) {
1591   GeomAbs_SurfaceType SurfaceType = TheSurfaceTool::GetType(surface);
1592   switch(SurfaceType) {
1593   case GeomAbs_Plane:
1594     {
1595       ElSLib::Parameters(TheSurfaceTool::Plane(surface),P,u,v);
1596         break;
1597       }
1598   case GeomAbs_Cylinder:
1599     {
1600       ElSLib::Parameters(TheSurfaceTool::Cylinder(surface),P,u,v);
1601         break;
1602       }
1603   case GeomAbs_Cone:
1604     {
1605       ElSLib::Parameters(TheSurfaceTool::Cone(surface),P,u,v);
1606         break;
1607       }
1608   case GeomAbs_Sphere:
1609     {
1610       ElSLib::Parameters(TheSurfaceTool::Sphere(surface),P,u,v);
1611         break;
1612       }
1613 #ifndef DEB
1614   default:      break;
1615 #endif
1616   }
1617 }
1618 //=======================================================================
1619 //function : EstLimForInfExtr
1620 //purpose  : Estimation of limits for infinite surfaces
1621 //=======================================================================
1622 void EstLimForInfExtr(const gp_Lin&   Line,
1623                       const TheSurface& surface,
1624                       const Standard_Boolean IsOffSurf,
1625                       const Standard_Integer nbsu,
1626                       const Standard_Boolean U1inf,
1627                       const Standard_Boolean U2inf,
1628                       const Standard_Boolean V1inf,
1629                       const Standard_Boolean V2inf,
1630                       Standard_Real& U1new,
1631                       Standard_Real& U2new,
1632                       Standard_Real& V1new,
1633                       Standard_Real& V2new,
1634                       Standard_Boolean& NoIntersection)
1635
1636 {
1637
1638   NoIntersection = Standard_False;
1639
1640   Handle(Adaptor3d_HSurface) aBasSurf;
1641
1642   if(IsOffSurf) aBasSurf = TheSurfaceTool::BasisSurface(surface);
1643
1644   gp_Dir aDirOfExt;
1645
1646   if(IsOffSurf) aDirOfExt = aBasSurf->Direction();
1647   else aDirOfExt = TheSurfaceTool::Direction(surface);
1648
1649   Standard_Real tolang = TOLERANCE_ANGULAIRE;
1650
1651   if(aDirOfExt.IsParallel(Line.Direction(), tolang)) {
1652     NoIntersection = Standard_True;
1653     return;
1654   }
1655
1656   if((V1inf || V2inf) && !(U1inf || U2inf)) {
1657
1658     Standard_Real vmin = RealLast(), vmax = -vmin;
1659     gp_Lin aL;
1660     Standard_Real step = (U2new - U1new) / nbsu;
1661     Standard_Real u = U1new, v;
1662     gp_Pnt aP;
1663     Extrema_POnCurv aP1, aP2;
1664     Standard_Integer i;
1665
1666     for(i = 0; i <= nbsu; i++) {
1667
1668       TheSurfaceTool::D0(surface, u, 0., aP);
1669       aL.SetLocation(aP);
1670       aL.SetDirection(aDirOfExt);
1671
1672       Extrema_ExtElC aExtr(aL, Line, tolang);
1673
1674       if(!aExtr.IsDone()) return;
1675
1676       if(aExtr.IsParallel()) {
1677         NoIntersection = Standard_True;
1678         return;
1679       }
1680
1681       aExtr.Points(1, aP1, aP2);
1682       v = aP1.Parameter();
1683       vmin = Min(vmin, v);
1684       vmax = Max(vmax, v);
1685
1686       u += step;
1687
1688     }
1689
1690     vmin = vmin - Abs(vmin) - 10.;
1691     vmax = vmax + Abs(vmax) + 10.;
1692
1693     V1new = Max(V1new, vmin);
1694     V2new = Min(V2new, vmax);
1695
1696   }
1697   else if(U1inf || U2inf) {
1698
1699     Standard_Real umin = RealLast(), umax = -umin;
1700     Standard_Real u0 = Min(Max(0., U1new), U2new);
1701     Standard_Real v0 = Min(Max(0., V1new), V2new);
1702     gp_Pnt aP;
1703     TheSurfaceTool::D0(surface, u0, v0, aP);
1704     gp_Pln aRefPln(aP, aDirOfExt);
1705
1706     Handle(Adaptor3d_HCurve) aBasCurv;
1707
1708     if(IsOffSurf) aBasCurv = aBasSurf->BasisCurve();
1709     else aBasCurv = TheSurfaceTool::BasisCurve(surface);
1710
1711     ProjLib_Plane Projector(aRefPln);
1712
1713     Projector.Project(Line);
1714
1715     if(!Projector.IsDone()) return;
1716
1717     gp_Lin2d Line2d = Projector.Line();
1718
1719     GeomAbs_CurveType aCurvTyp = aBasCurv->GetType();
1720
1721     if(aCurvTyp == GeomAbs_Line) {
1722
1723       Projector.Project(aBasCurv->Line());
1724
1725       if(!Projector.IsDone()) return;
1726
1727       gp_Lin2d aL2d = Projector.Line();
1728
1729       IntAna2d_AnaIntersection anInter(Line2d, aL2d);
1730
1731       if(!anInter.IsDone()) return;
1732
1733       if(anInter.IsEmpty() || anInter.IdenticalElements() ||
1734                               anInter.ParallelElements()     ) {
1735         NoIntersection = Standard_True;
1736         return;
1737       }
1738
1739       const IntAna2d_IntPoint& anIntPnt = anInter.Point(1);
1740       umin = umax = anIntPnt.ParamOnSecond();
1741     }
1742     else if(aCurvTyp == GeomAbs_Parabola || aCurvTyp == GeomAbs_Hyperbola) {
1743
1744       IntAna2d_Conic aCon(Line2d);
1745       IntAna2d_AnaIntersection anInter;
1746
1747       if(aCurvTyp == GeomAbs_Parabola) {
1748         Projector.Project(aBasCurv->Parabola());
1749         if(!Projector.IsDone()) return;
1750
1751         const gp_Parab2d& aP2d = Projector.Parabola();
1752
1753         anInter.Perform(aP2d, aCon);
1754       }
1755       else {
1756         Projector.Project(aBasCurv->Hyperbola());
1757         if(!Projector.IsDone()) return;
1758
1759         const gp_Hypr2d& aH2d = Projector.Hyperbola();
1760         anInter.Perform(aH2d, aCon);
1761       }
1762
1763       if(!anInter.IsDone()) return;
1764
1765       if(anInter.IsEmpty()) {
1766         NoIntersection = Standard_True;
1767         return;
1768       }
1769
1770       Standard_Integer i, nbint = anInter.NbPoints();
1771       for(i = 1; i <= nbint; i++) {
1772
1773         const IntAna2d_IntPoint& anIntPnt = anInter.Point(i);
1774
1775         umin = Min(anIntPnt.ParamOnFirst(), umin);
1776         umax = Max(anIntPnt.ParamOnFirst(), umax);
1777
1778       }
1779
1780     }
1781     else {
1782       return;
1783     }
1784
1785     umin = umin - Abs(umin) - 10;
1786     umax = umax + Abs(umax) + 10;
1787
1788     U1new = Max(U1new, umin);
1789     U2new = Min(U2new, umax);
1790
1791     if(V1inf || V2inf) {
1792       EstLimForInfExtr(Line, surface, IsOffSurf, nbsu,
1793                        Standard_False, Standard_False, V1inf, V2inf,
1794                        U1new, U2new, V1new, V2new, NoIntersection);
1795     }
1796   }
1797
1798   return;
1799
1800 }
1801 //=======================================================================
1802 //function : ProjectIntersectAndEstLim
1803 //purpose  : project <theLine> and it's X-axe symmetric line to <thePln> and
1804 //           intersect resulting curve with <theBasCurvProj>.
1805 //           Then estimate max and min parameters of intersection on
1806 //           <theBasCurvProj>.
1807 //           Is called from EstLimForInfRevl()
1808 //=======================================================================
1809 void ProjectIntersectAndEstLim(const gp_Lin&        theLine,
1810                                const gp_Pln&        thePln,
1811                                const ProjLib_Plane& theBasCurvProj,
1812                                Standard_Real&       theVmin,
1813                                Standard_Real&       theVmax,
1814                                Standard_Boolean&    theNoIntersection)
1815 {
1816   ProjLib_Plane aLineProj( thePln, theLine );
1817   if (!aLineProj.IsDone()) {
1818 #ifdef DEB
1819   cout
1820   << "Info: IntCurveSurface_Inter::ProjectIntersectAndEstLim(), !aLineProj.IsDone()"
1821   << endl;
1822 #endif
1823     return;
1824   }
1825   gp_Lin2d aLin2d = aLineProj.Line();
1826
1827   // make a second line X-axe symmetric to the first one
1828   gp_Pnt2d aP1 = aLin2d.Location();
1829   gp_Pnt2d aP2 (aP1.XY() + aLin2d.Direction().XY());
1830   gp_Pnt2d aP1sym (aP1.X(), -aP1.Y());
1831   gp_Pnt2d aP2sym (aP2.X(), -aP2.Y());
1832   gp_Lin2d aLin2dsym (aP1sym, gp_Vec2d(aP1sym,aP2sym));
1833
1834   // intersect projections
1835   IntAna2d_Conic aCon    (aLin2d);
1836   IntAna2d_Conic aConSym (aLin2dsym);
1837   IntAna2d_AnaIntersection anIntersect;
1838   IntAna2d_AnaIntersection anIntersectSym;
1839
1840   switch (theBasCurvProj.GetType()) {
1841   case GeomAbs_Line:
1842     anIntersectSym.Perform(theBasCurvProj.Line(), aConSym);
1843     anIntersect.Perform(theBasCurvProj.Line(), aCon); break;
1844   case GeomAbs_Hyperbola:
1845     anIntersectSym.Perform(theBasCurvProj.Hyperbola(), aConSym);
1846     anIntersect.Perform(theBasCurvProj.Hyperbola(), aCon); break;
1847   case GeomAbs_Parabola:
1848     anIntersectSym.Perform(theBasCurvProj.Parabola(), aConSym);
1849     anIntersect.Perform(theBasCurvProj.Parabola(), aCon); break;
1850   default:
1851     return; // not infinite curve
1852   }
1853
1854   // retrieve params of intersections
1855   Standard_Integer aNbIntPnt = anIntersect.IsDone() ? anIntersect.NbPoints() : 0;
1856   Standard_Integer aNbIntPntSym = anIntersectSym.IsDone() ? anIntersectSym.NbPoints() : 0;
1857   Standard_Integer iPnt, aNbPnt = Max (aNbIntPnt, aNbIntPntSym);
1858
1859   if (aNbPnt == 0) {
1860     theNoIntersection = Standard_True;
1861     return;
1862   }
1863   Standard_Real aParam;
1864   for (iPnt = 1; iPnt <= aNbPnt; iPnt++)
1865   {
1866     if (iPnt <= aNbIntPnt) {
1867       const IntAna2d_IntPoint& aIntPnt = anIntersect.Point(iPnt);
1868       aParam = aIntPnt.ParamOnFirst();
1869       theVmin = Min (theVmin, aParam );
1870       theVmax = Max (theVmax, aParam );
1871     }
1872     if (iPnt <= aNbIntPntSym) {
1873       const IntAna2d_IntPoint& aIntPnt = anIntersectSym.Point(iPnt);
1874       aParam = aIntPnt.ParamOnFirst();
1875       theVmin = Min (theVmin, aParam );
1876       theVmax = Max (theVmax, aParam );
1877     }
1878   }
1879
1880   return;
1881 }
1882 //=======================================================================
1883 //function : EstLimForInfRevl
1884 //purpose  : Estimate V1 and V2 to pass to InternalPerform() if they are
1885 //           infinite for Surface of Revolution
1886 //           Algo: intersect projections of Line and basis curve on the
1887 //           plane passing through revolution axe
1888 //=======================================================================
1889 void EstLimForInfRevl(const gp_Lin&   Line,
1890                       const TheSurface& surface,
1891                       const Standard_Boolean U1inf,
1892                       const Standard_Boolean U2inf,
1893                       const Standard_Boolean V1inf,
1894                       const Standard_Boolean V2inf,
1895                       Standard_Real& U1new,
1896                       Standard_Real& U2new,
1897                       Standard_Real& V1new,
1898                       Standard_Real& V2new,
1899                       Standard_Boolean& NoIntersection)
1900 {
1901
1902   NoIntersection = Standard_False;
1903
1904   if (U1inf || U2inf) {
1905     if (U1inf)
1906       U1new = Max (0., U1new);
1907     else
1908       U2new = Min (2 * M_PI, U2new);
1909     if (! V1inf && !V2inf) return;
1910   }
1911
1912   Handle(Adaptor3d_HCurve) aBasisCurve = TheSurfaceTool::BasisCurve(surface);
1913   gp_Ax1 aRevAx = TheSurfaceTool::AxeOfRevolution(surface);
1914   gp_Vec aXVec = aRevAx.Direction();
1915   Standard_Real aTolAng = Precision::Angular();
1916
1917   // make plane to project a basis curve
1918   gp_Pnt O = aRevAx.Location();
1919   Standard_Real aU = 0.;
1920   gp_Pnt P = aBasisCurve->Value(aU);
1921   while (O.SquareDistance(P) <= Precision::PConfusion() ||
1922          aXVec.IsParallel( gp_Vec(O,P), aTolAng)) {
1923     aU += 1.;
1924     P = aBasisCurve->Value(aU);
1925     if (aU > 3)
1926       // basis curve is a line coinciding with aXVec, P is any not on aXVec
1927       P = gp_Pnt(aU, aU+1, aU+2);
1928   }
1929   gp_Vec aNVec = aXVec ^ gp_Vec(O,P);
1930   gp_Pln aPln (gp_Ax3 (O, aNVec ,aXVec));
1931
1932   // project basic curve
1933   ProjLib_Plane aBasCurvProj(aPln);
1934   switch (aBasisCurve->GetType()) {
1935   case GeomAbs_Line:
1936     aBasCurvProj.Project(aBasisCurve->Line     ()); break;
1937   case GeomAbs_Hyperbola:
1938     aBasCurvProj.Project(aBasisCurve->Hyperbola()); break;
1939   case GeomAbs_Parabola:
1940     aBasCurvProj.Project(aBasisCurve->Parabola ()); break;
1941   default:
1942     return; // not infinite curve
1943   }
1944   if (!aBasCurvProj.IsDone()) {
1945 #ifdef DEB
1946     cout << "Info: IntCurveSurface_Inter::EstLimForInfRevl(), !aBasCurvProj.IsDone()" << endl;
1947 #endif
1948     return;
1949   }
1950   // make plane to project Line
1951   if (aXVec.IsParallel( Line.Direction(), aTolAng)) {
1952     P = Line.Location();
1953     while (O.SquareDistance(P) <= Precision::PConfusion()) {
1954       aU += 1.;
1955       P = gp_Pnt(aU, aU+1, aU+2); // any not on aXVec
1956     }
1957     aNVec = aXVec ^ gp_Vec(O,P);
1958   } else
1959     aNVec = aXVec.Crossed( Line.Direction() );
1960
1961   aPln  = gp_Pln (gp_Ax3 (O, aNVec ,aXVec));
1962
1963   // make a second plane perpendicular to the first one, rotated around aXVec
1964   gp_Pln aPlnPrp = aPln.Rotated (gp_Ax1 (O,aXVec), M_PI/2.);
1965
1966   // project Line and it's X-axe symmetric one to plane and intersect
1967   // resulting curve with projection of Basic Curev
1968   Standard_Real aVmin = RealLast(), aVmax = -aVmin;
1969   Standard_Boolean aNoInt1 = Standard_False, aNoInt2 = Standard_False;
1970   ProjectIntersectAndEstLim (Line, aPln,    aBasCurvProj, aVmin, aVmax, aNoInt1);
1971   ProjectIntersectAndEstLim (Line, aPlnPrp, aBasCurvProj, aVmin, aVmax, aNoInt2);
1972
1973   if (aNoInt1 && aNoInt2) {
1974     NoIntersection = Standard_True;
1975     return;
1976   }
1977
1978   aVmin = aVmin - Abs(aVmin) - 10;
1979   aVmax = aVmax + Abs(aVmax) + 10;
1980
1981   if (V1inf) V1new = aVmin;
1982   if (V2inf) V2new = aVmax;
1983
1984   //cout << "EstLimForInfRevl: Vmin " << V1new << " Vmax " << V2new << endl;
1985
1986   return;
1987 }
1988
1989 //=======================================================================
1990 //function : EstLimForInfOffs
1991 //purpose  :
1992 //=======================================================================
1993 void EstLimForInfOffs(const gp_Lin&   Line,
1994                       const TheSurface& surface,
1995                       const Standard_Integer nbsu,
1996                       const Standard_Boolean U1inf,
1997                       const Standard_Boolean U2inf,
1998                       const Standard_Boolean V1inf,
1999                       const Standard_Boolean V2inf,
2000                       Standard_Real& U1new,
2001                       Standard_Real& U2new,
2002                       Standard_Real& V1new,
2003                       Standard_Real& V2new,
2004                       Standard_Boolean& NoIntersection)
2005 {
2006
2007   NoIntersection = Standard_False;
2008
2009   const Handle(Adaptor3d_HSurface)& aBasSurf = TheSurfaceTool::BasisSurface(surface);
2010   Standard_Real anOffVal = TheSurfaceTool::OffsetValue(surface);
2011
2012   GeomAbs_SurfaceType aTypeOfBasSurf = aBasSurf->GetType();
2013
2014   //  case for plane, cylinder and cone - make equivalent surface;
2015   if(aTypeOfBasSurf == GeomAbs_Plane) {
2016
2017     gp_Pln aPln = aBasSurf->Plane();
2018     gp_Vec aT = aPln.Position().XDirection()^aPln.Position().YDirection();
2019     aT *= anOffVal;
2020     aPln.Translate(aT);
2021     IntAna_IntConicQuad LinPlane(Line,aPln,TOLERANCE_ANGULAIRE);
2022
2023     if(!LinPlane.IsDone()) return;
2024
2025     if(LinPlane.IsParallel() || LinPlane.IsInQuadric()) {
2026
2027       NoIntersection = Standard_True;
2028       return;
2029
2030     }
2031
2032     Standard_Real u, v;
2033     ElSLib::Parameters(aPln, LinPlane.Point(1), u, v);
2034     U1new = Max(U1new, u - 10.);
2035     U2new = Min(U2new, u + 10.);
2036     V1new = Max(V1new, v - 10.);
2037     V2new = Min(V2new, v + 10.);
2038
2039   }
2040   else if(aTypeOfBasSurf == GeomAbs_Cylinder) {
2041
2042     gp_Cylinder aCyl = aBasSurf->Cylinder();
2043
2044     Standard_Real aR = aCyl.Radius();
2045     gp_Ax3 anA = aCyl.Position();
2046
2047     if (anA.Direct())
2048       aR += anOffVal;
2049     else
2050       aR -= anOffVal;
2051
2052     if ( aR >= TOLTANGENCY ) {
2053       aCyl.SetRadius(aR);
2054     }
2055     else if ( aR <= -TOLTANGENCY ){
2056       anA.Rotate(gp_Ax1(anA.Location(), anA.Direction()), M_PI);
2057       aCyl.SetPosition(anA);
2058 // modified by NIZHNY-MKK  Mon Oct  3 17:37:54 2005
2059 //       aCyl.SetRadius(Abs(aR));
2060       aCyl.SetRadius(-aR);
2061     }
2062     else {
2063
2064       NoIntersection = Standard_True;
2065       return;
2066
2067     }
2068
2069     IntAna_IntConicQuad LinCylinder(Line, aCyl);
2070
2071     if(!LinCylinder.IsDone()) return;
2072
2073     if(LinCylinder.IsParallel() || LinCylinder.IsInQuadric()) {
2074
2075       NoIntersection = Standard_True;
2076       return;
2077
2078     }
2079
2080     Standard_Integer i, nbp = LinCylinder.NbPoints();
2081     Standard_Real vmin = RealLast(), vmax = -vmin, u, v;
2082
2083     for(i = 1; i <= nbp; i++) {
2084
2085       ElSLib::Parameters(aCyl, LinCylinder.Point(i), u, v);
2086       vmin = Min(vmin, v);
2087       vmax = Max(vmax, v);
2088
2089     }
2090
2091     V1new = Max(V1new, vmin - Abs(vmin) - 10.);
2092     V2new = Min(V2new, vmax + Abs(vmax) + 10.);
2093
2094   }
2095   else if(aTypeOfBasSurf == GeomAbs_Cone) {
2096
2097     gp_Cone aCon = aBasSurf->Cone();
2098     Standard_Real anAng = aCon.SemiAngle();
2099     Standard_Real aR = aCon.RefRadius() + anOffVal * Cos(anAng);
2100     gp_Ax3 anA = aCon.Position();
2101     if ( aR >= 0.) {
2102
2103       gp_Vec aZ( anA.Direction());
2104       aZ *= - anOffVal * Sin(anAng);
2105       anA.Translate(aZ);
2106       aCon.SetPosition(anA);
2107       aCon.SetRadius(aR);
2108       aCon.SetSemiAngle(anAng);
2109
2110     }
2111     else {
2112
2113       return;
2114
2115     }
2116
2117     IntAna_IntConicQuad LinCone(Line, aCon);
2118
2119     if(!LinCone.IsDone()) return;
2120
2121     if(LinCone.IsParallel() || LinCone.IsInQuadric()) {
2122
2123       NoIntersection = Standard_True;
2124       return;
2125
2126     }
2127
2128     Standard_Integer i, nbp = LinCone.NbPoints();
2129     Standard_Real vmin = RealLast(), vmax = -vmin, u, v;
2130
2131     for(i = 1; i <= nbp; i++) {
2132
2133       ElSLib::Parameters(aCon, LinCone.Point(i), u, v);
2134       vmin = Min(vmin, v);
2135       vmax = Max(vmax, v);
2136
2137     }
2138
2139     V1new = Max(V1new, vmin - Abs(vmin) - 10.);
2140     V2new = Min(V2new, vmax + Abs(vmax) + 10.);
2141
2142   }
2143   else if(aTypeOfBasSurf == GeomAbs_SurfaceOfExtrusion) {
2144
2145     Standard_Real anU1 = U1new, anU2 = U2new;
2146
2147     EstLimForInfExtr(Line, surface, Standard_True, nbsu,
2148                      U1inf, U2inf, V1inf, V2inf,
2149                      U1new, U2new, V1new, V2new, NoIntersection);
2150
2151     if(NoIntersection) return;
2152
2153     if(U1inf || U2inf) {
2154
2155       GeomAbs_CurveType aBasCurvType = aBasSurf->BasisCurve()->GetType();
2156       if(aBasCurvType == GeomAbs_Line) {
2157         U1new = Max(anU1, -1.e10);
2158         U2new = Min(anU2,  1.e10);
2159       }
2160       else if(aBasCurvType == GeomAbs_Parabola) {
2161         gp_Parab aPrb  = aBasSurf->BasisCurve()->Parabola();
2162         Standard_Real aF = aPrb.Focal();
2163         Standard_Real dU = 2.e5 * Sqrt(aF);
2164         U1new = Max(anU1, -dU);
2165         U2new = Min(anU2,  dU);
2166       }
2167       else if(aBasCurvType == GeomAbs_Hyperbola) {
2168         U1new = Max(anU1, -30.);
2169         U2new = Min(anU2,  30.);
2170       }
2171       else {
2172         U1new = Max(anU1, -1.e10);
2173         U2new = Min(anU2,  1.e10);
2174       }
2175
2176     }
2177
2178   }
2179   else if(aTypeOfBasSurf == GeomAbs_SurfaceOfRevolution) {
2180
2181     GeomAbs_CurveType aBasCurvType = aBasSurf->BasisCurve()->GetType();
2182     if(aBasCurvType == GeomAbs_Line) {
2183       V1new = Max(V1new, -1.e10);
2184       V2new = Min(V2new,  1.e10);
2185     }
2186     else if(aBasCurvType == GeomAbs_Parabola) {
2187       gp_Parab aPrb  = aBasSurf->BasisCurve()->Parabola();
2188       Standard_Real aF = aPrb.Focal();
2189       Standard_Real dV = 2.e5 * Sqrt(aF);
2190       V1new = Max(V1new, -dV);
2191       V2new = Min(V2new,  dV);
2192     }
2193     else if(aBasCurvType == GeomAbs_Hyperbola) {
2194      V1new  = Max(V1new, -30.);
2195      V2new  = Min(V2new,  30.);
2196     }
2197     else {
2198      V1new  = Max(V1new, -1.e10);
2199      V2new  = Min(V2new,  1.e10);
2200     }
2201
2202   }
2203   else {
2204
2205     V1new = Max(V1new, -1.e10);
2206     V2new = Min(V2new,  1.e10);
2207   }
2208 }
2209 //=======================================================================
2210 //function : EstLimForInfSurf
2211 //purpose  :
2212 //=======================================================================
2213 void EstLimForInfSurf(Standard_Real& U1new,
2214                       Standard_Real& U2new,
2215                       Standard_Real& V1new,
2216                       Standard_Real& V2new)
2217 {
2218     U1new = Max(U1new, -1.e10);
2219     U2new = Min(U2new,  1.e10);
2220     V1new = Max(V1new, -1.e10);
2221     V2new = Min(V2new,  1.e10);
2222 }