src/GeomLib/GeomLib.cxx

   1 // File:        GeomLib.cxx
   2 // Created:     Wed Jul  7 15:32:09 1993
   3 // Author:      Jean Claude VAUTHIER
   4
   5 // Copyright:    Matra Datavision
   6 // Version:
   7 // History:     pmn 24/09/96 Ajout du prolongement de courbe.
   8 //              jct 15/04/97 Ajout du prolongement de surface.
   9 //              jct 24/04/97 simplification ou suppression de calculs
  10 //                           inutiles dans ExtendSurfByLength
  11 //                           correction de Tbord et Continuity=0 accepte
  12 //                           correction du calcul de lambda et appel a
  13 //                           TangExtendToConstraint avec lambmin au lieu de 1.
  14 //                           correction du passage Sr rat --> BSp nD
  15 //              xab 26/06/97 treatement partiel anulation des derivees
  16 //                           partiels du denonimateur des Surfaces BSplines Rationnelles
  17 //                           dans le cas de valeurs proportionnelles des denominateurs
  18 //                           en umin umax et/ou vmin vmax.
  19 //              pmn 4/07/97  Gestion de la continuite dans BuildCurve3d (PRO9097)
  20
  21 //              xab 10/07/97 on revient en arriere sur l'ajout du 26/06/97
  22 //              pmn 26/09/97 Ajout des parametres d'approx dans BuildCurve3d
  23 //              xab 29/09/97 on reintegre l'ajout du 26/06/97
  24 //              pmn 31/10/97 Ajoute AdjustExtremity
  25 //              jct 26/11/98 blindage dans ExtendSurf qd NTgte = 0 (CTS21288)
  26 //              jct 19/01/99 traitement de la periodicite dans ExtendSurf
  27 // Design:
  28 // Warning:      None
  29 // References:   None
  30 // Language:     C++2.0
  31 // Purpose:
  32
  33 // Declarations:
  34
  35 #include <GeomLib.ixx>
  36
  37 #include <Precision.hxx>
  38 #include <GeomConvert.hxx>
  39 #include <Hermit.hxx>
  40 #include <Standard_NotImplemented.hxx>
  41 #include <GeomLib_MakeCurvefromApprox.hxx>
  42 #include <GeomLib_DenominatorMultiplier.hxx>
  43 #include <GeomLib_DenominatorMultiplierPtr.hxx>
  44 #include <GeomLib_PolyFunc.hxx>
  45 #include <GeomLib_LogSample.hxx>
  46
  47 #include <AdvApprox_ApproxAFunction.hxx>
  48 #include <AdvApprox_PrefAndRec.hxx>
  49
  50 #include <Adaptor2d_HCurve2d.hxx>
  51 #include <Adaptor3d_HCurve.hxx>
  52 #include <Adaptor3d_HSurface.hxx>
  53 #include <Adaptor3d_CurveOnSurface.hxx>
  54 #include <Geom2dAdaptor_Curve.hxx>
  55 #include <GeomAdaptor_Surface.hxx>
  56 #include <GeomAdaptor_HSurface.hxx>
  57 #include <Geom2dAdaptor_HCurve.hxx>
  58 #include <Geom2dAdaptor_GHCurve.hxx>
  59
  60 #include <Geom2d_BSplineCurve.hxx>
  61 #include <Geom_BSplineCurve.hxx>
  62 #include <Geom2d_BezierCurve.hxx>
  63 #include <Geom_BezierCurve.hxx>
  64 #include <Geom_RectangularTrimmedSurface.hxx>
  65 #include <Geom_Plane.hxx>
  66 #include <Geom_Line.hxx>
  67 #include <Geom2d_Line.hxx>
  68 #include <Geom_Circle.hxx>
  69 #include <Geom2d_Circle.hxx>
  70 #include <Geom_Ellipse.hxx>
  71 #include <Geom2d_Ellipse.hxx>
  72 #include <Geom_Parabola.hxx>
  73 #include <Geom2d_Parabola.hxx>
  74 #include <Geom_Hyperbola.hxx>
  75 #include <Geom2d_Hyperbola.hxx>
  76 #include <Geom_TrimmedCurve.hxx>
  77 #include <Geom2d_TrimmedCurve.hxx>
  78 #include <Geom_OffsetCurve.hxx>
  79 #include <Geom2d_OffsetCurve.hxx>
  80 #include <Geom_BezierSurface.hxx>
  81 #include <Geom_BSplineSurface.hxx>
  82
  83 #include <BSplCLib.hxx>
  84 #include <BSplSLib.hxx>
  85 #include <PLib.hxx>
  86 #include <math_Matrix.hxx>
  87 #include <math_Vector.hxx>
  88 #include <math_Jacobi.hxx>
  89 #include <math.hxx>
  90 #include <math_FunctionAllRoots.hxx>
  91 #include <math_FunctionSample.hxx>
  92
  93 #include <TColStd_HArray1OfReal.hxx>
  94 #include <TColgp_Array1OfPnt.hxx>
  95 #include <TColgp_Array1OfVec.hxx>
  96 #include <TColgp_Array2OfPnt.hxx>
  97 #include <TColgp_HArray2OfPnt.hxx>
  98 #include <TColgp_Array1OfPnt2d.hxx>
  99 #include <TColgp_Array1OfXYZ.hxx>
 100 #include <TColStd_Array1OfReal.hxx>
 101 #include <TColStd_Array2OfReal.hxx>
 102 #include <TColStd_HArray2OfReal.hxx>
 103 #include <TColStd_Array1OfInteger.hxx>
 104
 105 #include <gp_TrsfForm.hxx>
 106 #include <gp_Lin.hxx>
 107 #include <gp_Lin2d.hxx>
 108 #include <gp_Circ.hxx>
 109 #include <gp_Circ2d.hxx>
 110 #include <gp_Elips.hxx>
 111 #include <gp_Elips2d.hxx>
 112 #include <gp_Hypr.hxx>
 113 #include <gp_Hypr2d.hxx>
 114 #include <gp_Parab.hxx>
 115 #include <gp_Parab2d.hxx>
 116 #include <gp_GTrsf2d.hxx>
 117 #include <gp_Trsf2d.hxx>
 118
 119 #include <ElCLib.hxx>
 120 #include <Geom2dConvert.hxx>
 121 #include <GeomConvert_CompCurveToBSplineCurve.hxx>
 122 #include <GeomConvert_ApproxSurface.hxx>
 123
 124
 125 #include <Standard_ConstructionError.hxx>
 126
 127 //=======================================================================
 128 //function : ComputeLambda
 129 //purpose  : Calcul le facteur lambda qui minimise la variation de vittesse
 130 //           sur une interpolation d'hermite d'ordre (i,0)
 131 //=======================================================================
 132 static void ComputeLambda(const math_Matrix& Constraint,
 133                           const math_Matrix& Hermit,
 134                           const Standard_Real Length,
 135                           Standard_Real& Lambda )
 136 {
 137   Standard_Integer size = Hermit.RowNumber();
 138   Standard_Integer Continuity = size-2;
 139   Standard_Integer ii, jj, ip, pp;
 140
 141   //Minimization
 142   math_Matrix HDer(1, size-1, 1, size);
 143   for (jj=1; jj<=size; jj++) {
 144     for (ii=1; ii<size;ii++) {
 145       HDer(ii, jj) = ii*Hermit(jj, ii+1);
 146     }
 147   }
 148
 149   math_Vector V(1, size);
 150   math_Vector Vec1(1, Constraint.RowNumber());
 151   math_Vector Vec2(1, Constraint.RowNumber());
 152   math_Vector Vec3(1, Constraint.RowNumber());
 153   math_Vector Vec4(1, Constraint.RowNumber());
 154
 155   Standard_Real * polynome = &HDer(1,1);
 156   Standard_Real * valhder =  &V(1);
 157   Vec2 =  Constraint.Col(2);
 158   Vec2 /= Length;
 159   Standard_Real t,  squared1 = Vec2.Norm2(), GW;
 160 //  math_Matrix Vec(1, Constraint.RowNumber(), 1, size-1);
 161 //  gp_Vec Vfirst(p0.XYZ()), Vlast(Point.XYZ());
 162 //  TColgp_Array1OfVec Der(2, 4);
 163 //  Der(2) = d1; Der(3) = d2; Der(4) = d3;
 164
 165   Standard_Integer GOrdre = 4 + 4*Continuity,
 166                    DDim=Continuity*(Continuity+2);
 167   math_Vector GaussP(1, GOrdre), GaussW(1, GOrdre),
 168               pol2(1, 2*Continuity+1),
 169               pol4(1, 4*Continuity+1);
 170   math::GaussPoints(GOrdre, GaussP);
 171   math::GaussWeights (GOrdre, GaussW);
 172   pol4.Init(0.);
 173
 174   for (ip=1; ip<=GOrdre; ip++) {
 175     t = (GaussP(ip)+1.)/2;
 176     GW = GaussW(ip);
 177     PLib::NoDerivativeEvalPolynomial(t ,  Continuity, Continuity+2, DDim,
 178                                      polynome[0], valhder[0]);
 179     V /= Length; //Normalisation
 180
 181     //                      i
 182     // C'(t) = SUM Vi*Lambda
 183     Vec1 = Constraint.Col(1);
 184     Vec1 *= V(1);
 185     Vec1 += V(size)*Constraint.Col(size);
 186     Vec2 = Constraint.Col(2);
 187     Vec2 *= V(2);
 188     if (Continuity > 1) {
 189       Vec3 = Constraint.Col(3);
 190       Vec3 *= V(3);
 191       if (Continuity > 2) {
 192         Vec4 = Constraint.Col(4);
 193         Vec4 *= V(4);
 194       }
 195     }
 196
 197
 198     //   2          2
 199     // C'(t) - C'(0)
 200
 201     pol2(1) = Vec1.Norm2();
 202     pol2(2) = 2*(Vec1.Multiplied(Vec2));
 203     pol2(3) = Vec2.Norm2() - squared1;
 204     if (Continuity>1) {
 205       pol2(3) += 2*(Vec1.Multiplied(Vec3));
 206       pol2(4) =  2*(Vec2.Multiplied(Vec3));
 207       pol2(5) =  Vec3.Norm2();
 208       if (Continuity>2) {
 209         pol2(4)+= 2*(Vec1.Multiplied(Vec4));
 210         pol2(5)+= 2*(Vec2.Multiplied(Vec4));
 211         pol2(6) = 2*(Vec3.Multiplied(Vec4));
 212         pol2(7) = Vec4.Norm2();
 213       }
 214     }
 215
 216     //                     2      2  2
 217     // Integrale de ( C'(t) - C'(0) )
 218     for (ii=1; ii<=pol2.Length(); ii++) {
 219       pp = ii;
 220       for(jj=1; jj<ii; jj++, pp++) {
 221         pol4(pp) += 2*GW*pol2(ii)*pol2(jj);
 222       }
 223       pol4(2*ii-1) += GW*Pow(pol2(ii), 2);
 224     }
 225   }
 226
 227   Standard_Real EMin, E;
 228   PLib::NoDerivativeEvalPolynomial(Lambda , pol4.Length()-1, 1,
 229                                    pol4.Length()-1,
 230                                    pol4(1), EMin);
 231
 232   if (EMin > Precision::Confusion()) {
 233     // Recheche des extrema de la fonction
 234     GeomLib_PolyFunc FF(pol4);
 235     GeomLib_LogSample S(Lambda/1000, 50*Lambda, 100);
 236     math_FunctionAllRoots Solve(FF, S, Precision::Confusion(),
 237                                 Precision::Confusion()*(Length+1),
 238                                 1.e-15);
 239     if (Solve.IsDone()) {
 240       for (ii=1; ii<=Solve.NbPoints(); ii++) {
 241         t = Solve.GetPoint(ii);
 242         PLib::NoDerivativeEvalPolynomial(t , pol4.Length()-1, 1,
 243                                          pol4.Length()-1,
 244                                          pol4(1), E);
 245         if (E < EMin) {
 246           Lambda = t;
 247           EMin = E;
 248         }
 249       }
 250     }
 251   }
 252 }
 253
 254 #include <Extrema_LocateExtPC.hxx>
 255 //=======================================================================
 256 //function : RemovePointsFromArray
 257 //purpose  :
 258 //=======================================================================
 259
 260 void GeomLib::RemovePointsFromArray(const Standard_Integer NumPoints,
 261                                     const TColStd_Array1OfReal& InParameters,
 262                                     Handle(TColStd_HArray1OfReal)& OutParameters)
 263 {
 264  Standard_Integer ii,
 265    jj,
 266    add_one_point,
 267    loc_num_points,
 268    num_points,
 269    index ;
 270  Standard_Real delta,
 271    current_parameter ;
 272
 273    loc_num_points = Max(0,NumPoints-2) ;
 274    delta = InParameters(InParameters.Upper()) - InParameters(InParameters.Lower()) ;
 275    delta /= (Standard_Real) (loc_num_points + 1) ;
 276    num_points = 1 ;
 277    current_parameter = InParameters(InParameters.Lower()) + delta * 0.5e0 ;
 278    ii = InParameters.Lower() + 1 ;
 279    for (jj = 0 ; ii < InParameters.Upper() && jj < NumPoints ; jj++) {
 280      add_one_point = 0 ;
 281      while ( ii < InParameters.Upper() && InParameters(ii) < current_parameter) {
 282        ii += 1 ;
 283        add_one_point = 1 ;
 284      }
 285      num_points += add_one_point ;
 286      current_parameter += delta ;
 287    }
 288    if (NumPoints <= 2) {
 289      num_points = 2 ;
 290    }
 291    index = 2 ;
 292    current_parameter = InParameters(InParameters.Lower()) + delta * 0.5e0 ;
 293    OutParameters =
 294      new TColStd_HArray1OfReal(1,num_points) ;
 295    OutParameters->ChangeArray1()(1) = InParameters(InParameters.Lower()) ;
 296    ii = InParameters.Lower() + 1 ;
 297    for (jj = 0 ; ii < InParameters.Upper() && jj < NumPoints ; jj++) {
 298      add_one_point = 0 ;
 299      while (ii < InParameters.Upper() && InParameters(ii) < current_parameter) {
 300        ii += 1 ;
 301        add_one_point = 1 ;
 302      }
 303      if (add_one_point && index <= num_points) {
 304        OutParameters->ChangeArray1()(index) = InParameters(ii-1) ;
 305        index += 1 ;
 306      }
 307      current_parameter += delta ;
 308    }
 309    OutParameters->ChangeArray1()(num_points) = InParameters(InParameters.Upper()) ;
 310 }
 311 //=======================================================================
 312 //function : DensifyArray1OfReal
 313 //purpose  :
 314 //=======================================================================
 315
 316 void GeomLib::DensifyArray1OfReal(const Standard_Integer MinNumPoints,
 317                                   const TColStd_Array1OfReal& InParameters,
 318                                   Handle(TColStd_HArray1OfReal)& OutParameters)
 319 {
 320  Standard_Integer ii,
 321    in_order,
 322    num_points,
 323    num_parameters_to_add,
 324    index ;
 325  Standard_Real delta,
 326    current_parameter ;
 327
 328  in_order = 1 ;
 329  if (MinNumPoints > InParameters.Length()) {
 330
 331    //
 332    // checks the paramaters are in increasing order
 333    //
 334    for (ii = InParameters.Lower() ; ii < InParameters.Upper() ; ii++) {
 335      if (InParameters(ii) > InParameters(ii+1)) {
 336        in_order = 0 ;
 337        break ;
 338      }
 339    }
 340    if (in_order) {
 341      num_parameters_to_add = MinNumPoints - InParameters.Length()  ;
 342      delta = InParameters(InParameters.Upper()) - InParameters(InParameters.Lower()) ;
 343      delta /= (Standard_Real) (num_parameters_to_add + 1) ;
 344      num_points = MinNumPoints ;
 345      OutParameters =
 346        new TColStd_HArray1OfReal(1,num_points) ;
 347      index = 1 ;
 348      current_parameter = InParameters(InParameters.Lower()) ;
 349      OutParameters->ChangeArray1()(index) = current_parameter ;
 350      index += 1 ;
 351      current_parameter += delta ;
 352      for (ii = InParameters.Lower() + 1 ; index <= num_points && ii <= InParameters.Upper() ; ii++) {
 353        while (current_parameter < InParameters(ii) && index <= num_points) {
 354          OutParameters->ChangeArray1()(index) = current_parameter ;
 355          index += 1 ;
 356          current_parameter += delta ;
 357        }
 358        if (index <= num_points) {
 359          OutParameters->ChangeArray1()(index) = InParameters(ii) ;
 360        }
 361        index += 1 ;
 362      }
 363      //
 364      // beware of roundoff !
 365      //
 366      OutParameters->ChangeArray1()(num_points) = InParameters(InParameters.Upper()) ;
 367    }
 368    else {
 369      index = 1 ;
 370      num_points = InParameters.Length() ;
 371      OutParameters =
 372        new TColStd_HArray1OfReal(1,num_points) ;
 373      for (ii = InParameters.Lower()  ; ii <= InParameters.Upper() ; ii++) {
 374        OutParameters->ChangeArray1()(index) = InParameters(ii) ;
 375        index += 1 ;
 376      }
 377    }
 378  }
 379  else {
 380    index = 1 ;
 381    num_points = InParameters.Length() ;
 382    OutParameters =
 383      new TColStd_HArray1OfReal(1,num_points) ;
 384    for (ii = InParameters.Lower()  ; ii <= InParameters.Upper() ; ii++) {
 385      OutParameters->ChangeArray1()(index) = InParameters(ii) ;
 386      index += 1 ;
 387    }
 388  }
 389 }
 390
 391 //=======================================================================
 392 //function : FuseIntervals
 393 //purpose  :
 394 //=======================================================================
 395 void GeomLib::FuseIntervals(const  TColStd_Array1OfReal& I1,
 396                             const  TColStd_Array1OfReal& I2,
 397                             TColStd_SequenceOfReal&  Seq,
 398                             const Standard_Real  Epspar)
 399 {
 400  Standard_Integer ind1=1, ind2=1;
 401  Standard_Real    v1, v2;
 402 // Initialisations : les IND1 et IND2 pointent sur le 1er element
 403 // de chacune des 2 tables a traiter.INDS pointe sur le dernier
 404 // element cree de TABSOR
 405
 406
 407 //--- On remplit TABSOR en parcourant TABLE1 et TABLE2 simultanement ---
 408 //------------------ en eliminant les occurrences multiples ------------
 409
 410  while ((ind1<=I1.Upper()) && (ind2<=I2.Upper())) {
 411       v1 = I1(ind1);
 412       v2 = I2(ind2);
 413       if (Abs(v1-v2)<= Epspar) {
 414 // Ici les elements de I1 et I2 conviennent .
 415          Seq.Append((v1+v2)/2);
 416          ind1++;
 417          ind2++;
 418        }
 419       else if (v1 < v2) {
 420         // Ici l' element de I1 convient.
 421          Seq.Append(v1);
 422          ind1++;
 423        }
 424       else {
 425 // Ici l' element de TABLE2 convient.
 426          Seq.Append(v2);
 427          ind2++;
 428        }
 429     }
 430
 431   if (ind1>I1.Upper()) {
 432 //----- Ici I1 est epuise, on complete avec la fin de TABLE2 -------
 433
 434     for (; ind2<=I2.Upper(); ind2++) {
 435       Seq.Append(I2(ind2));
 436     }
 437   }
 438
 439   if (ind2>I2.Upper()) {
 440 //----- Ici I2 est epuise, on complete avec la fin de I1 -------
 441     for (; ind1<=I1.Upper(); ind1++) {
 442       Seq.Append(I1(ind1));
 443     }
 444   }
 445 }
 446
 447
 448 //=======================================================================
 449 //function : EvalMaxParametricDistance
 450 //purpose  :
 451 //=======================================================================
 452
 453 void GeomLib::EvalMaxParametricDistance(const Adaptor3d_Curve& ACurve,
 454                                const Adaptor3d_Curve& AReferenceCurve,
 455 //                             const Standard_Real  Tolerance,
 456                                const Standard_Real  ,
 457                                const TColStd_Array1OfReal& Parameters,
 458                                Standard_Real& MaxDistance)
 459 {
 460   Standard_Integer ii ;
 461
 462   Standard_Real max_squared = 0.0e0,
 463 //    tolerance_squared,
 464     local_distance_squared ;
 465
 466 //  tolerance_squared = Tolerance * Tolerance ;
 467   gp_Pnt Point1 ;
 468   gp_Pnt Point2 ;
 469   for (ii = Parameters.Lower() ; ii <= Parameters.Upper() ; ii++) {
 470     ACurve.D0(Parameters(ii),
 471               Point1) ;
 472     AReferenceCurve.D0(Parameters(ii),
 473                        Point2) ;
 474     local_distance_squared =
 475       Point1.SquareDistance (Point2) ;
 476     max_squared = Max(max_squared,local_distance_squared) ;
 477   }
 478   if (max_squared > 0.0e0) {
 479     MaxDistance = sqrt(max_squared) ;
 480   }
 481   else {
 482     MaxDistance = 0.0e0 ;
 483   }
 484
 485 }
 486 //=======================================================================
 487 //function : EvalMaxDistanceAlongParameter
 488 //purpose  :
 489 //=======================================================================
 490
 491 void GeomLib::EvalMaxDistanceAlongParameter(const Adaptor3d_Curve& ACurve,
 492                                const Adaptor3d_Curve& AReferenceCurve,
 493                                const Standard_Real  Tolerance,
 494                                const TColStd_Array1OfReal& Parameters,
 495                                Standard_Real& MaxDistance)
 496 {
 497   Standard_Integer ii ;
 498   Standard_Real max_squared = 0.0e0,
 499     tolerance_squared = Tolerance * Tolerance,
 500     other_parameter,
 501     para_tolerance,
 502     local_distance_squared ;
 503   gp_Pnt Point1 ;
 504   gp_Pnt Point2 ;
 505
 506
 507
 508   para_tolerance =
 509     AReferenceCurve.Resolution(Tolerance) ;
 510   other_parameter = Parameters(Parameters.Lower()) ;
 511   ACurve.D0(other_parameter,
 512             Point1) ;
 513   Extrema_LocateExtPC a_projector(Point1,
 514                                   AReferenceCurve,
 515                                   other_parameter,
 516                                   para_tolerance) ;
 517   for (ii = Parameters.Lower() ; ii <= Parameters.Upper() ; ii++) {
 518     ACurve.D0(Parameters(ii),
 519               Point1) ;
 520     AReferenceCurve.D0(Parameters(ii),
 521                        Point2) ;
 522     local_distance_squared =
 523       Point1.SquareDistance (Point2) ;
 524
 525     local_distance_squared =
 526       Point1.SquareDistance (Point2) ;
 527
 528
 529     if (local_distance_squared > tolerance_squared) {
 530
 531
 532       a_projector.Perform(Point1,
 533                           other_parameter) ;
 534       if (a_projector.IsDone()) {
 535         other_parameter =
 536           a_projector.Point().Parameter() ;
 537         AReferenceCurve.D0(other_parameter,
 538                            Point2) ;
 539         local_distance_squared =
 540           Point1.SquareDistance (Point2) ;
 541       }
 542       else {
 543         local_distance_squared = 0.0e0 ;
 544         other_parameter = Parameters(ii) ;
 545       }
 546     }
 547     else {
 548       other_parameter = Parameters(ii) ;
 549     }
 550
 551
 552     max_squared = Max(max_squared,local_distance_squared) ;
 553   }
 554   if (max_squared > tolerance_squared) {
 555     MaxDistance = sqrt(max_squared) ;
 556   }
 557   else {
 558     MaxDistance = Tolerance ;
 559   }
 560 }
 561
 562
 563
 564 // Aliases:
 565
 566 // Global data definitions:
 567
 568 // Methods :
 569
 570
 571 //=======================================================================
 572 //function : To3d
 573 //purpose  :
 574 //=======================================================================
 575
 576 Handle(Geom_Curve) GeomLib::To3d (const gp_Ax2&               Position,
 577                                   const Handle(Geom2d_Curve)& Curve2d  ) {
 578   Handle(Geom_Curve) Curve3d;
 579   Handle(Standard_Type) KindOfCurve = Curve2d->DynamicType();
 580
 581   if (KindOfCurve == STANDARD_TYPE (Geom2d_TrimmedCurve)) {
 582     Handle(Geom2d_TrimmedCurve) Ct =
 583       Handle(Geom2d_TrimmedCurve)::DownCast(Curve2d);
 584     Standard_Real U1 = Ct->FirstParameter ();
 585     Standard_Real U2 = Ct->LastParameter  ();
 586     Handle(Geom2d_Curve) CBasis2d = Ct->BasisCurve();
 587     Handle(Geom_Curve) CC = GeomLib::To3d(Position, CBasis2d);
 588     Curve3d = new Geom_TrimmedCurve (CC, U1, U2);
 589   }
 590   else if (KindOfCurve == STANDARD_TYPE (Geom2d_OffsetCurve)) {
 591     Handle(Geom2d_OffsetCurve) Co =
 592       Handle(Geom2d_OffsetCurve)::DownCast(Curve2d);
 593     Standard_Real Offset = Co->Offset();
 594     Handle(Geom2d_Curve) CBasis2d = Co->BasisCurve();
 595     Handle(Geom_Curve) CC = GeomLib::To3d(Position, CBasis2d);
 596     Curve3d = new Geom_OffsetCurve (CC, Offset, Position.Direction());
 597   }
 598   else if (KindOfCurve == STANDARD_TYPE (Geom2d_BezierCurve)) {
 599     Handle(Geom2d_BezierCurve) CBez2d =
 600       Handle(Geom2d_BezierCurve)::DownCast (Curve2d);
 601     Standard_Integer Nbpoles = CBez2d->NbPoles ();
 602     TColgp_Array1OfPnt2d Poles2d (1, Nbpoles);
 603     CBez2d->Poles (Poles2d);
 604     TColgp_Array1OfPnt Poles3d (1, Nbpoles);
 605     for (Standard_Integer i = 1; i <= Nbpoles; i++) {
 606       Poles3d (i) = ElCLib::To3d (Position, Poles2d (i));
 607     }
 608     Handle(Geom_BezierCurve) CBez3d;
 609     if (CBez2d->IsRational()) {
 610       TColStd_Array1OfReal TheWeights (1, Nbpoles);
 611       CBez2d->Weights (TheWeights);
 612       CBez3d = new Geom_BezierCurve (Poles3d, TheWeights);
 613     }
 614     else {
 615       CBez3d = new Geom_BezierCurve (Poles3d);
 616     }
 617     Curve3d = CBez3d;
 618   }
 619   else if (KindOfCurve == STANDARD_TYPE (Geom2d_BSplineCurve)) {
 620     Handle(Geom2d_BSplineCurve) CBSpl2d =
 621       Handle(Geom2d_BSplineCurve)::DownCast (Curve2d);
 622     Standard_Integer Nbpoles   = CBSpl2d->NbPoles ();
 623     Standard_Integer Nbknots   = CBSpl2d->NbKnots ();
 624     Standard_Integer TheDegree = CBSpl2d->Degree ();
 625     Standard_Boolean IsPeriodic = CBSpl2d->IsPeriodic();
 626     TColgp_Array1OfPnt2d Poles2d (1, Nbpoles);
 627     CBSpl2d->Poles (Poles2d);
 628     TColgp_Array1OfPnt Poles3d (1, Nbpoles);
 629     for (Standard_Integer i = 1; i <= Nbpoles; i++) {
 630       Poles3d (i) = ElCLib::To3d (Position, Poles2d (i));
 631     }
 632     TColStd_Array1OfReal    TheKnots (1, Nbknots);
 633     TColStd_Array1OfInteger TheMults (1, Nbknots);
 634     CBSpl2d->Knots (TheKnots);
 635     CBSpl2d->Multiplicities (TheMults);
 636     Handle(Geom_BSplineCurve) CBSpl3d;
 637     if (CBSpl2d->IsRational()) {
 638       TColStd_Array1OfReal TheWeights (1, Nbpoles);
 639       CBSpl2d->Weights (TheWeights);
 640       CBSpl3d = new Geom_BSplineCurve (Poles3d, TheWeights, TheKnots, TheMults, TheDegree, IsPeriodic);
 641     }
 642     else {
 643       CBSpl3d = new Geom_BSplineCurve (Poles3d, TheKnots, TheMults, TheDegree, IsPeriodic);
 644     }
 645     Curve3d = CBSpl3d;
 646   }
 647   else if (KindOfCurve == STANDARD_TYPE (Geom2d_Line)) {
 648     Handle(Geom2d_Line) Line2d = Handle(Geom2d_Line)::DownCast (Curve2d);
 649     gp_Lin2d L2d = Line2d->Lin2d();
 650     gp_Lin   L3d = ElCLib::To3d (Position, L2d);
 651     Handle(Geom_Line) GeomL3d = new Geom_Line (L3d);
 652     Curve3d = GeomL3d;
 653   }
 654   else if (KindOfCurve == STANDARD_TYPE (Geom2d_Circle)) {
 655     Handle(Geom2d_Circle) Circle2d =
 656       Handle(Geom2d_Circle)::DownCast (Curve2d);
 657     gp_Circ2d C2d = Circle2d->Circ2d();
 658     gp_Circ   C3d = ElCLib::To3d (Position, C2d);
 659     Handle(Geom_Circle) GeomC3d = new Geom_Circle (C3d);
 660     Curve3d = GeomC3d;
 661   }
 662   else if (KindOfCurve == STANDARD_TYPE (Geom2d_Ellipse)) {
 663     Handle(Geom2d_Ellipse) Ellipse2d =
 664       Handle(Geom2d_Ellipse)::DownCast (Curve2d);
 665     gp_Elips2d E2d = Ellipse2d->Elips2d ();
 666     gp_Elips   E3d = ElCLib::To3d (Position, E2d);
 667     Handle(Geom_Ellipse) GeomE3d = new Geom_Ellipse (E3d);
 668     Curve3d = GeomE3d;
 669   }
 670   else if (KindOfCurve == STANDARD_TYPE (Geom2d_Parabola)) {
 671     Handle(Geom2d_Parabola) Parabola2d =
 672       Handle(Geom2d_Parabola)::DownCast (Curve2d);
 673     gp_Parab2d Prb2d = Parabola2d->Parab2d ();
 674     gp_Parab   Prb3d = ElCLib::To3d (Position, Prb2d);
 675     Handle(Geom_Parabola) GeomPrb3d = new Geom_Parabola (Prb3d);
 676     Curve3d = GeomPrb3d;
 677   }
 678   else if (KindOfCurve == STANDARD_TYPE (Geom2d_Hyperbola)) {
 679     Handle(Geom2d_Hyperbola) Hyperbola2d =
 680       Handle(Geom2d_Hyperbola)::DownCast (Curve2d);
 681     gp_Hypr2d H2d = Hyperbola2d->Hypr2d ();
 682     gp_Hypr   H3d = ElCLib::To3d (Position, H2d);
 683     Handle(Geom_Hyperbola) GeomH3d = new Geom_Hyperbola (H3d);
 684     Curve3d = GeomH3d;
 685   }
 686   else {
 687     Standard_NotImplemented::Raise();
 688   }
 689
 690   return Curve3d;
 691 }
 692
 693
 694
 695 //=======================================================================
 696 //function : GTransform
 697 //purpose  :
 698 //=======================================================================
 699
 700 Handle(Geom2d_Curve) GeomLib::GTransform(const Handle(Geom2d_Curve)& Curve,
 701                                          const gp_GTrsf2d&           GTrsf)
 702 {
 703   gp_TrsfForm Form = GTrsf.Form();
 704
 705   if ( Form != gp_Other) {
 706
 707     // Alors, la GTrsf est en fait une Trsf.
 708     // La geometrie des courbes sera alors inchangee.
 709
 710     Handle(Geom2d_Curve) C =
 711       Handle(Geom2d_Curve)::DownCast(Curve->Transformed(GTrsf.Trsf2d()));
 712     return C;
 713   }
 714   else {
 715
 716     // Alors, la GTrsf est une other Transformation.
 717     // La geometrie des courbes est alors changee, et les conics devront
 718     // etre converties en BSplines.
 719
 720     Handle(Standard_Type) TheType = Curve->DynamicType();
 721
 722     if ( TheType == STANDARD_TYPE(Geom2d_TrimmedCurve)) {
 723
 724       // On va recurer sur la BasisCurve
 725
 726       Handle(Geom2d_TrimmedCurve) C =
 727         Handle(Geom2d_TrimmedCurve)::DownCast(Curve->Copy());
 728
 729       Handle(Standard_Type) TheBasisType = (C->BasisCurve())->DynamicType();
 730
 731       if (TheBasisType == STANDARD_TYPE(Geom2d_BSplineCurve) ||
 732           TheBasisType == STANDARD_TYPE(Geom2d_BezierCurve)    ) {
 733
 734         // Dans ces cas le parametrage est conserve sur la courbe transformee
 735         // on peut donc la trimmer avec les parametres de la courbe de base.
 736
 737         Standard_Real U1 = C->FirstParameter();
 738         Standard_Real U2 = C->LastParameter();
 739
 740         Handle(Geom2d_TrimmedCurve) result =
 741           new Geom2d_TrimmedCurve(GTransform(C->BasisCurve(), GTrsf), U1,U2);
 742         return result;
 743       }
 744       else if ( TheBasisType == STANDARD_TYPE(Geom2d_Line)) {
 745
 746         // Dans ce cas, le parametrage n`est plus conserve.
 747         // Il faut recalculer les parametres de Trimming sur la courbe
 748         // resultante. ( Calcul par projection ( ElCLib) des points debut
 749         // et fin transformes)
 750
 751         Handle(Geom2d_Line) L =
 752           Handle(Geom2d_Line)::DownCast(GTransform(C->BasisCurve(), GTrsf));
 753         gp_Lin2d Lin = L->Lin2d();
 754
 755         gp_Pnt2d P1 = C->StartPoint();
 756         gp_Pnt2d P2 = C->EndPoint();
 757         P1.SetXY(GTrsf.Transformed(P1.XY()));
 758         P2.SetXY(GTrsf.Transformed(P2.XY()));
 759         Standard_Real U1 = ElCLib::Parameter(Lin,P1);
 760         Standard_Real U2 = ElCLib::Parameter(Lin,P2);
 761
 762         Handle(Geom2d_TrimmedCurve) result =
 763           new Geom2d_TrimmedCurve(L,U1,U2);
 764         return result;
 765       }
 766       else if (TheBasisType == STANDARD_TYPE(Geom2d_Circle)   ||
 767                TheBasisType == STANDARD_TYPE(Geom2d_Ellipse)  ||
 768                TheBasisType == STANDARD_TYPE(Geom2d_Parabola) ||
 769                TheBasisType == STANDARD_TYPE(Geom2d_Hyperbola)  ) {
 770
 771         // Dans ces cas, la geometrie de la courbe n`est pas conservee
 772         // on la convertir en BSpline avant de lui appliquer la Trsf.
 773
 774         Handle(Geom2d_BSplineCurve) BS =
 775           Geom2dConvert::CurveToBSplineCurve(C);
 776         return GTransform(BS,GTrsf);
 777       }
 778       else {
 779
 780         // La transformee d`une OffsetCurve vaut ????? Sais pas faire !!
 781
 782         Handle(Geom2d_Curve) dummy;
 783         return dummy;
 784       }
 785     }
 786     else if ( TheType == STANDARD_TYPE(Geom2d_Line)) {
 787
 788       Handle(Geom2d_Line) L =
 789         Handle(Geom2d_Line)::DownCast(Curve->Copy());
 790       gp_Lin2d Lin = L->Lin2d();
 791       gp_Pnt2d P  = Lin.Location();
 792       gp_Pnt2d PP = L->Value(10.); // pourquoi pas !!
 793       P.SetXY(GTrsf.Transformed(P.XY()));
 794       PP.SetXY(GTrsf.Transformed(PP.XY()));
 795       L->SetLocation(P);
 796       gp_Vec2d V(P,PP);
 797       L->SetDirection(gp_Dir2d(V));
 798       return L;
 799     }
 800     else if ( TheType == STANDARD_TYPE(Geom2d_BezierCurve)) {
 801
 802       // Les GTrsf etant des operation lineaires, la transformee d`une courbe
 803       // a poles est la courbe dont les poles sont la transformee des poles
 804       // de la courbe de base.
 805
 806       Handle(Geom2d_BezierCurve) C =
 807         Handle(Geom2d_BezierCurve)::DownCast(Curve->Copy());
 808       Standard_Integer NbPoles = C->NbPoles();
 809       TColgp_Array1OfPnt2d Poles(1,NbPoles);
 810       C->Poles(Poles);
 811       for ( Standard_Integer i = 1; i <= NbPoles; i++) {
 812         Poles(i).SetXY(GTrsf.Transformed(Poles(i).XY()));
 813         C->SetPole(i,Poles(i));
 814       }
 815       return C;
 816     }
 817     else if ( TheType == STANDARD_TYPE(Geom2d_BSplineCurve)) {
 818
 819       // Voir commentaire pour les Bezier.
 820
 821       Handle(Geom2d_BSplineCurve) C =
 822         Handle(Geom2d_BSplineCurve)::DownCast(Curve->Copy());
 823       Standard_Integer NbPoles = C->NbPoles();
 824       TColgp_Array1OfPnt2d Poles(1,NbPoles);
 825       C->Poles(Poles);
 826       for ( Standard_Integer i = 1; i <= NbPoles; i++) {
 827         Poles(i).SetXY(GTrsf.Transformed(Poles(i).XY()));
 828         C->SetPole(i,Poles(i));
 829       }
 830       return C;
 831     }
 832     else if ( TheType == STANDARD_TYPE(Geom2d_Circle) ||
 833               TheType == STANDARD_TYPE(Geom2d_Ellipse)  ) {
 834
 835       // Dans ces cas, la geometrie de la courbe n`est pas conservee
 836       // on la convertir en BSpline avant de lui appliquer la Trsf.
 837
 838       Handle(Geom2d_BSplineCurve) C =
 839         Geom2dConvert::CurveToBSplineCurve(Curve);
 840       return GTransform(C, GTrsf);
 841     }
 842     else if ( TheType == STANDARD_TYPE(Geom2d_Parabola)   ||
 843               TheType == STANDARD_TYPE(Geom2d_Hyperbola)  ||
 844               TheType == STANDARD_TYPE(Geom2d_OffsetCurve)  ) {
 845
 846       // On ne sait pas faire : return a null Handle;
 847
 848       Handle(Geom2d_Curve) dummy;
 849       return dummy;
 850     }
 851   }
 852
 853   Handle(Geom2d_Curve) WNT__; // portage Windows.
 854   return WNT__;
 855 }
 856
 857
 858 //=======================================================================
 859 //function : SameRange
 860 //purpose  :
 861 //=======================================================================
 862 void GeomLib::SameRange(const Standard_Real         Tolerance,
 863                         const Handle(Geom2d_Curve)& CurvePtr,
 864                         const Standard_Real         FirstOnCurve,
 865                         const Standard_Real         LastOnCurve,
 866                         const Standard_Real         RequestedFirst,
 867                         const Standard_Real         RequestedLast,
 868                               Handle(Geom2d_Curve)& NewCurvePtr)
 869 {
 870   if(CurvePtr.IsNull()) Standard_Failure::Raise();
 871   if (Abs(LastOnCurve - RequestedLast) <= Tolerance &&
 872       Abs(FirstOnCurve - RequestedFirst) <= Tolerance) {
 873     NewCurvePtr = CurvePtr;
 874     return;
 875   }
 876
 877   // the parametrisation lentgh  must at least be the same.
 878   if (Abs(LastOnCurve - FirstOnCurve - RequestedLast + RequestedFirst)
 879       <= Tolerance) {
 880     if (CurvePtr->IsKind(STANDARD_TYPE(Geom2d_Line))) {
 881       Handle(Geom2d_Line) Line =
 882         Handle(Geom2d_Line)::DownCast(CurvePtr->Copy());
 883       Standard_Real dU = FirstOnCurve - RequestedFirst;
 884       gp_Dir2d D = Line->Direction() ;
 885       Line->Translate(dU * gp_Vec2d(D));
 886       NewCurvePtr = Line;
 887     }
 888     else if (CurvePtr->IsKind(STANDARD_TYPE(Geom2d_Circle))) {
 889       gp_Trsf2d Trsf;
 890       NewCurvePtr = Handle(Geom2d_Curve)::DownCast(CurvePtr->Copy());
 891       Handle(Geom2d_Circle) Circ =
 892         Handle(Geom2d_Circle)::DownCast(NewCurvePtr);
 893       gp_Pnt2d P = Circ->Location();
 894       Standard_Real dU;
 895       if (Circ->Circ2d().IsDirect()) {
 896         dU = FirstOnCurve - RequestedFirst;
 897       }
 898       else {
 899         dU = RequestedFirst - FirstOnCurve;
 900       }
 901       Trsf.SetRotation(P,dU);
 902       NewCurvePtr->Transform(Trsf) ;
 903     }
 904     else if (CurvePtr->IsKind(STANDARD_TYPE(Geom2d_TrimmedCurve))) {
 905       Handle(Geom2d_TrimmedCurve) TC =
 906         Handle(Geom2d_TrimmedCurve)::DownCast(CurvePtr);
 907       GeomLib::SameRange(Tolerance,
 908                          TC->BasisCurve(),
 909                          FirstOnCurve  , LastOnCurve,
 910                          RequestedFirst, RequestedLast,
 911                          NewCurvePtr);
 912       NewCurvePtr = new Geom2d_TrimmedCurve( NewCurvePtr, RequestedFirst, RequestedLast );
 913     }
 914 //
 915 //  attention a des problemes de limitation : utiliser le MEME test que dans
 916 //  Geom2d_TrimmedCurve::SetTrim car sinon comme on risque de relimite sur
 917 //  RequestedFirst et RequestedLast on aura un probleme
 918 //
 919 //
 920     else if (Abs(LastOnCurve - FirstOnCurve) > Precision::PConfusion() ||
 921              Abs(RequestedLast + RequestedFirst) > Precision::PConfusion()) {
 922
 923       Handle(Geom2d_TrimmedCurve) TC =
 924         new Geom2d_TrimmedCurve(CurvePtr,FirstOnCurve,LastOnCurve);
 925
 926       Handle(Geom2d_BSplineCurve) BS =
 927         Geom2dConvert::CurveToBSplineCurve(TC);
 928       TColStd_Array1OfReal Knots(1,BS->NbKnots());
 929       BS->Knots(Knots);
 930
 931       BSplCLib::Reparametrize(RequestedFirst,RequestedLast,Knots);
 932
 933       BS->SetKnots(Knots);
 934       NewCurvePtr = BS;
 935     }
 936
 937   }
 938   else { // On segmente le resultat
 939     Handle(Geom2d_TrimmedCurve) TC =
 940       new Geom2d_TrimmedCurve( CurvePtr, FirstOnCurve, LastOnCurve );
 941
 942     Standard_Real newFirstOnCurve = TC->FirstParameter(), newLastOnCurve = TC->LastParameter();
 943
 944     Handle(Geom2d_BSplineCurve) BS =
 945       Geom2dConvert::CurveToBSplineCurve(TC);
 946
 947     if (BS->IsPeriodic())
 948       BS->Segment( newFirstOnCurve, newLastOnCurve) ;
 949     else
 950       BS->Segment( Max(newFirstOnCurve, BS->FirstParameter()),
 951                    Min(newLastOnCurve,  BS->LastParameter()) );
 952
 953     TColStd_Array1OfReal Knots(1,BS->NbKnots());
 954     BS->Knots(Knots);
 955
 956     BSplCLib::Reparametrize(RequestedFirst,RequestedLast,Knots);
 957
 958     BS->SetKnots(Knots);
 959     NewCurvePtr = BS;
 960   }
 961 }
 962
 963 //=======================================================================
 964 //class : GeomLib_CurveOnSurfaceEvaluator
 965 //purpose: The evaluator for the Curve 3D building
 966 //=======================================================================
 967
 968 class GeomLib_CurveOnSurfaceEvaluator : public AdvApprox_EvaluatorFunction
 969 {
 970  public:
 971   GeomLib_CurveOnSurfaceEvaluator (Adaptor3d_CurveOnSurface& theCurveOnSurface,
 972                                    Standard_Real theFirst, Standard_Real theLast)
 973     : CurveOnSurface(theCurveOnSurface), FirstParam(theFirst), LastParam(theLast) {}
 974
 975   virtual void Evaluate (Standard_Integer *Dimension,
 976                          Standard_Real     StartEnd[2],
 977                          Standard_Real    *Parameter,
 978                          Standard_Integer *DerivativeRequest,
 979                          Standard_Real    *Result, // [Dimension]
 980                          Standard_Integer *ErrorCode);
 981
 982  private:
 983   Adaptor3d_CurveOnSurface& CurveOnSurface;
 984   Standard_Real FirstParam;
 985   Standard_Real LastParam;
 986
 987   Handle(Adaptor3d_HCurve) TrimCurve;
 988 };
 989
 990 void GeomLib_CurveOnSurfaceEvaluator::Evaluate (Standard_Integer *,/*Dimension*/
 991                                                 Standard_Real     DebutFin[2],
 992                                                 Standard_Real    *Parameter,
 993                                                 Standard_Integer *DerivativeRequest,
 994                                                 Standard_Real    *Result,// [Dimension]
 995                                                 Standard_Integer *ReturnCode)
 996 {
 997   register Standard_Integer ii ;
 998   gp_Pnt Point ;
 999
1000   //Gestion des positionnements gauche / droite
1001   if ((DebutFin[0] != FirstParam) || (DebutFin[1] != LastParam))
1002     {
1003       TrimCurve = CurveOnSurface.Trim(DebutFin[0], DebutFin[1], Precision::PConfusion());
1004       FirstParam = DebutFin[0];
1005       LastParam  = DebutFin[1];
1006     }
1007
1008   //Positionemment
1009   if (*DerivativeRequest == 0)
1010     {
1011      TrimCurve->D0((*Parameter), Point) ;
1012
1013      for (ii = 0 ; ii < 3 ; ii++)
1014        Result[ii] = Point.Coord(ii + 1);
1015    }
1016   if (*DerivativeRequest == 1)
1017     {
1018       gp_Vec Vector;
1019       TrimCurve->D1((*Parameter), Point, Vector);
1020       for (ii = 0 ; ii < 3 ; ii++)
1021         Result[ii] = Vector.Coord(ii + 1) ;
1022     }
1023   if (*DerivativeRequest == 2)
1024     {
1025       gp_Vec Vector, VecBis;
1026       TrimCurve->D2((*Parameter), Point, VecBis, Vector);
1027       for (ii = 0 ; ii < 3 ; ii++)
1028         Result[ii] = Vector.Coord(ii + 1) ;
1029     }
1030   ReturnCode[0] = 0;
1031 }
1032
1033 //=======================================================================
1034 //function : BuildCurve3d
1035 //purpose  :
1036 //=======================================================================
1037
1038 void GeomLib::BuildCurve3d(const Standard_Real           Tolerance,
1039                            Adaptor3d_CurveOnSurface&       Curve,
1040                            const Standard_Real           FirstParameter,
1041                            const Standard_Real           LastParameter,
1042                            Handle_Geom_Curve&            NewCurvePtr,
1043                            Standard_Real&                MaxDeviation,
1044                            Standard_Real&                AverageDeviation,
1045                            const GeomAbs_Shape           Continuity,
1046                            const Standard_Integer        MaxDegree,
1047                            const Standard_Integer        MaxSegment)
1048
1049 {
1050
1051
1052   Standard_Integer curve_not_computed = 1 ;
1053   MaxDeviation     = 0.0e0 ;
1054   AverageDeviation = 0.0e0 ;
1055   const Handle(GeomAdaptor_HSurface) &     geom_adaptor_surface_ptr =
1056   Handle(GeomAdaptor_HSurface)::DownCast(Curve.GetSurface()) ;
1057   const Handle(Geom2dAdaptor_HCurve) &     geom_adaptor_curve_ptr =
1058   Handle(Geom2dAdaptor_HCurve)::DownCast(Curve.GetCurve()) ;
1059
1060   if (! geom_adaptor_curve_ptr.IsNull() &&
1061       ! geom_adaptor_surface_ptr.IsNull()) {
1062      Handle(Geom_Plane) P ;
1063      const GeomAdaptor_Surface  &   geom_surface =
1064        * (GeomAdaptor_Surface *) &geom_adaptor_surface_ptr->Surface() ;
1065
1066     Handle(Geom_RectangularTrimmedSurface) RT =
1067       Handle(Geom_RectangularTrimmedSurface)::
1068         DownCast(geom_surface.Surface());
1069     if ( RT.IsNull()) {
1070       P = Handle(Geom_Plane)::DownCast(geom_surface.Surface());
1071     }
1072     else {
1073       P = Handle(Geom_Plane)::DownCast(RT->BasisSurface());
1074     }
1075
1076
1077     if (! P.IsNull()) {
1078       // compute the 3d curve
1079       gp_Ax2 axes = P->Position().Ax2();
1080       const Geom2dAdaptor_Curve & geom2d_curve =
1081         * (Geom2dAdaptor_Curve *) & geom_adaptor_curve_ptr->Curve2d() ;
1082       NewCurvePtr =
1083         GeomLib::To3d(axes,
1084                       geom2d_curve.Curve());
1085      curve_not_computed = 0 ;
1086
1087     }
1088   }
1089   if (curve_not_computed) {
1090
1091       //
1092       // Entree
1093       //
1094     Handle(TColStd_HArray1OfReal)   Tolerance1DPtr,Tolerance2DPtr;
1095     Handle(TColStd_HArray1OfReal) Tolerance3DPtr =
1096       new TColStd_HArray1OfReal(1,1) ;
1097     Tolerance3DPtr->SetValue(1,Tolerance);
1098
1099      // Recherche des discontinuitees
1100      Standard_Integer NbIntervalC2 = Curve.NbIntervals(GeomAbs_C2);
1101      TColStd_Array1OfReal Param_de_decoupeC2 (1, NbIntervalC2+1);
1102      Curve.Intervals(Param_de_decoupeC2, GeomAbs_C2);
1103
1104      Standard_Integer NbIntervalC3 = Curve.NbIntervals(GeomAbs_C3);
1105      TColStd_Array1OfReal Param_de_decoupeC3 (1, NbIntervalC3+1);
1106      Curve.Intervals(Param_de_decoupeC3, GeomAbs_C3);
1107
1108      // Note extension of the parameteric range
1109      // Pour forcer le Trim au premier appel de l'evaluateur
1110      GeomLib_CurveOnSurfaceEvaluator ev (Curve, FirstParameter - 1., LastParameter  + 1.);
1111
1112      // Approximation avec decoupe preferentiel
1113      AdvApprox_PrefAndRec Preferentiel(Param_de_decoupeC2,
1114                                        Param_de_decoupeC3);
1115      AdvApprox_ApproxAFunction  anApproximator(0,
1116                                               0,
1117                                               1,
1118                                               Tolerance1DPtr,
1119                                               Tolerance2DPtr,
1120                                               Tolerance3DPtr,
1121                                               FirstParameter,
1122                                               LastParameter,
1123                                               Continuity,
1124                                               MaxDegree,
1125                                               MaxSegment,
1126                                               ev,
1127 //                                            CurveOnSurfaceEvaluator,
1128                                               Preferentiel) ;
1129
1130     if (anApproximator.HasResult()) {
1131       GeomLib_MakeCurvefromApprox
1132         aCurveBuilder(anApproximator) ;
1133
1134       Handle(Geom_BSplineCurve) aCurvePtr =
1135         aCurveBuilder.Curve(1) ;
1136       // On rend les resultats de l'approx
1137       MaxDeviation = anApproximator.MaxError(3,1) ;
1138       AverageDeviation = anApproximator.AverageError(3,1) ;
1139       NewCurvePtr = aCurvePtr ;
1140     }
1141   }
1142  }
1143
1144 //=======================================================================
1145 //function :  AdjustExtremity
1146 //purpose  :
1147 //=======================================================================
1148
1149 void GeomLib::AdjustExtremity(Handle(Geom_BoundedCurve)& Curve,
1150                               const gp_Pnt& P1,
1151                               const gp_Pnt& P2,
1152                               const gp_Vec& T1,
1153                               const gp_Vec& T2)
1154 {
1155 // il faut Convertir l'entree (en preservant si possible le parametrage)
1156   Handle(Geom_BSplineCurve) aIn, aDef;
1157   aIn = GeomConvert::CurveToBSplineCurve(Curve, Convert_QuasiAngular);
1158
1159   Standard_Integer ii, jj;
1160   gp_Pnt P;
1161   gp_Vec V, Vtan, DV;
1162   TColgp_Array1OfPnt PolesDef(1,4), Coeffs(1,4);
1163   TColStd_Array1OfReal FK(1, 8);
1164   TColStd_Array1OfReal Ti(1, 4);
1165   TColStd_Array1OfInteger Contact(1, 4);
1166
1167   Ti(1) = Ti(2) = aIn->FirstParameter();
1168   Ti(3) = Ti(4) = aIn->LastParameter();
1169   Contact(1) =  Contact(3) = 0;
1170   Contact(2) =  Contact(4) = 1;
1171   for (ii=1; ii<=4; ii++) {
1172     FK(ii) = aIn->FirstParameter();
1173     FK(ii) = aIn->LastParameter();
1174   }
1175
1176   // Calculs des contraintes de deformations
1177   aIn->D1(Ti(1), P, V);
1178   PolesDef(1).ChangeCoord() = P1.XYZ()-P.XYZ();
1179   Vtan = T1;
1180   Vtan.Normalize();
1181   DV = Vtan * (Vtan * V) - V;
1182   PolesDef(2).ChangeCoord() = (Ti(4)-Ti(1))*DV.XYZ();
1183
1184   aIn->D1(Ti(4), P, V);
1185   PolesDef(3).ChangeCoord() = P2.XYZ()-P.XYZ();
1186   Vtan = T2;
1187   Vtan.Normalize();
1188   DV = Vtan * (Vtan * V) - V;
1189   PolesDef(4).ChangeCoord() = (Ti(4)-Ti(1))* DV.XYZ();
1190
1191   // Interpolation des contraintes
1192   math_Matrix Mat(1, 4, 1, 4);
1193   if (!PLib::HermiteCoefficients(0., 1., 1, 1, Mat))
1194     Standard_ConstructionError::Raise();
1195
1196   for (jj=1; jj<=4; jj++) {
1197     gp_XYZ aux(0.,0.,0.);
1198     for (ii=1; ii<=4; ii++) {
1199       aux.SetLinearForm(Mat(ii,jj), PolesDef(ii).XYZ(), aux);
1200     }
1201     Coeffs(jj).SetXYZ(aux);
1202   }
1203
1204   PLib::CoefficientsPoles(Coeffs, PLib::NoWeights(),
1205                           PolesDef,  PLib::NoWeights());
1206
1207   // Ajout de la deformation
1208   TColStd_Array1OfReal K(1, 2);
1209   TColStd_Array1OfInteger M(1, 2);
1210   K(1) = Ti(1);
1211   K(2) = Ti(4);
1212   M.Init(4);
1213
1214   aDef = new (Geom_BSplineCurve) (PolesDef, K, M, 3);
1215   if (aIn->Degree() < 3) aIn->IncreaseDegree(3);
1216   else aDef->IncreaseDegree(aIn->Degree());
1217
1218   for (ii=2; ii<aIn->NbKnots(); ii++) {
1219     aDef->InsertKnot(aIn->Knot(ii), aIn->Multiplicity(ii));
1220   }
1221
1222   if (aDef->NbPoles() != aIn->NbPoles())
1223     Standard_ConstructionError::Raise("Inconsistent poles's number");
1224
1225   for (ii=1; ii<=aDef->NbPoles(); ii++) {
1226     P = aIn->Pole(ii);
1227     P.ChangeCoord() += aDef->Pole(ii).XYZ();
1228     aIn->SetPole(ii, P);
1229   }
1230   Curve = aIn;
1231 }
1232 //=======================================================================
1233 //function : ExtendCurveToPoint
1234 //purpose  :
1235 //=======================================================================
1236
1237 void GeomLib::ExtendCurveToPoint(Handle(Geom_BoundedCurve)& Curve,
1238                                  const gp_Pnt& Point,
1239                                  const Standard_Integer Continuity,
1240                                  const Standard_Boolean After)
1241 {
1242   if(Continuity < 1 || Continuity > 3) return;
1243   Standard_Integer size = Continuity + 2;
1244   Standard_Real Ubord, Tol=1.e-6;
1245   math_Matrix  MatCoefs(1,size, 1,size);
1246   Standard_Real Lambda, L1;
1247   Standard_Integer ii, jj;
1248   gp_Vec d1, d2, d3;
1249   gp_Pnt p0;
1250 // il faut Convertir l'entree (en preservant si possible le parametrage)
1251   GeomConvert_CompCurveToBSplineCurve Concat(Curve, Convert_QuasiAngular);
1252
1253 // Les contraintes de constructions
1254   TColgp_Array1OfXYZ Cont(1,size);
1255   if (After) {
1256      Ubord = Curve->LastParameter();
1257
1258    }
1259   else {
1260      Ubord = Curve->FirstParameter();
1261    }
1262   PLib::HermiteCoefficients(0, 1,           // Les Bornes
1263                             Continuity, 0,  // Les Ordres de contraintes
1264                             MatCoefs);
1265
1266   Curve->D3(Ubord, p0, d1, d2, d3);
1267   if (!After) { // Inversion du parametrage
1268     d1 *= -1;
1269     d3 *= -1;
1270   }
1271
1272   L1 = p0.Distance(Point);
1273   if (L1 > Tol) {
1274     // Lambda est le ratio qu'il faut appliquer a la derive de la courbe
1275     // pour obtenir la derive du prolongement (fixe arbitrairement a la
1276     // longueur du segment bout de la courbe - point cible.
1277     // On essai d'avoir sur le prolongement la vitesse moyenne que l'on
1278     // a sur la courbe.
1279     gp_Vec daux;
1280     gp_Pnt pp;
1281     Standard_Real f= Curve->FirstParameter(), t, dt, norm;
1282     dt = (Curve->LastParameter()-f)/9;
1283     norm = d1.Magnitude();
1284     for (ii=1, t=f+dt; ii<=8; ii++, t+=dt) {
1285       Curve->D1(t, pp, daux);
1286       norm += daux.Magnitude();
1287     }
1288     norm /= 9;
1289     dt = d1.Magnitude() / norm;
1290     if ((dt<1.5) && (dt>0.75)) { // Le bord est dans la moyenne on le garde
1291       Lambda = ((Standard_Real)1) / Max (d1.Magnitude() / L1, Tol);
1292     }
1293     else {
1294       Lambda = ((Standard_Real)1) / Max (norm / L1, Tol);
1295     }
1296   }
1297   else {
1298     return; // Pas d'extension
1299   }
1300
1301   // Optimisation du Lambda
1302   math_Matrix Cons(1, 3, 1, size);
1303   Cons(1,1) = p0.X();  Cons(2,1) = p0.Y(); Cons(3,1) = p0.Z();
1304   Cons(1,2) = d1.X();  Cons(2,2) = d1.Y(); Cons(3,2) = d1.Z();
1305   Cons(1,size) = Point.X();  Cons(2,size) = Point.Y(); Cons(3,size) = Point.Z();
1306   if (Continuity >= 2) {
1307      Cons(1,3) = d2.X();  Cons(2,3) = d2.Y(); Cons(3,3) = d2.Z();
1308   }
1309   if (Continuity >= 3) {
1310      Cons(1,4) = d3.X();  Cons(2,4) = d3.Y(); Cons(3,4) = d3.Z();
1311   }
1312   ComputeLambda(Cons, MatCoefs, L1, Lambda);
1313
1314   // Construction dans la Base Polynomiale
1315   Cont(1) = p0.XYZ();
1316   Cont(2) = d1.XYZ() * Lambda;
1317   if(Continuity >= 2) Cont(3) = d2.XYZ() * Pow(Lambda,2);
1318   if(Continuity >= 3) Cont(4) = d3.XYZ() * Pow(Lambda,3);
1319   Cont(size) = Point.XYZ();
1320
1321
1322   TColgp_Array1OfPnt ExtrapPoles(1, size);
1323   TColgp_Array1OfPnt ExtraCoeffs(1, size);
1324
1325   gp_Pnt PNull(0.,0.,0.);
1326   ExtraCoeffs.Init(PNull);
1327   for (ii=1; ii<=size; ii++) {
1328     for (jj=1; jj<=size; jj++) {
1329       ExtraCoeffs(jj).ChangeCoord() += MatCoefs(ii,jj)*Cont(ii);
1330     }
1331   }
1332
1333   // Convertion Dans la Base de Bernstein
1334   PLib::CoefficientsPoles(ExtraCoeffs,  PLib::NoWeights(),
1335                           ExtrapPoles,  PLib::NoWeights());
1336
1337   Handle(Geom_BezierCurve) Bezier = new (Geom_BezierCurve) (ExtrapPoles);
1338
1339   Standard_Real dist = ExtrapPoles(1).Distance(p0);
1340   Standard_Boolean Ok;
1341   Tol += dist;
1342
1343   // Concatenation
1344   Ok = Concat.Add(Bezier, Tol, After);
1345   if (!Ok) Standard_ConstructionError::Raise("ExtendCurveToPoint");
1346
1347   Curve =  Concat.BSplineCurve();
1348 }
1349
1350
1351 //=======================================================================
1352 //function : ExtendKPart
1353 //purpose  : Extension par longueur des surfaces cannonique
1354 //=======================================================================
1355 static Standard_Boolean
1356 ExtendKPart(Handle(Geom_RectangularTrimmedSurface)& Surface,
1357             const Standard_Real Length,
1358             const Standard_Boolean InU,
1359             const Standard_Boolean After)
1360 {
1361
1362   if  (Surface.IsNull()) return Standard_False;
1363
1364   Standard_Boolean Ok=Standard_True;
1365   Standard_Real Uf, Ul, Vf, Vl;
1366   Handle(Geom_Surface) Support = Surface->BasisSurface();
1367   GeomAbs_SurfaceType Type;
1368
1369   Surface->Bounds(Uf, Ul, Vf, Vl);
1370   GeomAdaptor_Surface AS(Surface);
1371   Type = AS.GetType();
1372
1373   if (InU) {
1374     switch(Type) {
1375     case GeomAbs_Plane :
1376       {
1377         if (After) Ul+=Length;
1378         else       Uf-=Length;
1379         Surface = new (Geom_RectangularTrimmedSurface)
1380           (Support, Uf, Ul, Vf, Vl);
1381         break;
1382       }
1383
1384     default:
1385       Ok = Standard_False;
1386     }
1387   }
1388   else {
1389     switch(Type) {
1390     case GeomAbs_Plane :
1391     case GeomAbs_Cylinder :
1392     case GeomAbs_SurfaceOfExtrusion :
1393       {
1394         if (After) Vl+=Length;
1395         else       Vf-=Length;
1396         Surface = new (Geom_RectangularTrimmedSurface)
1397           (Support, Uf, Ul, Vf, Vl);
1398         break;
1399       }
1400     default:
1401       Ok = Standard_False;
1402     }
1403   }
1404
1405   return Ok;
1406 }
1407
1408 //=======================================================================
1409 //function : ExtendSurfByLength
1410 //purpose  :
1411 //=======================================================================
1412 void GeomLib::ExtendSurfByLength(Handle(Geom_BoundedSurface)& Surface,
1413                                  const Standard_Real Length,
1414                                  const Standard_Integer Continuity,
1415                                  const Standard_Boolean InU,
1416                                  const Standard_Boolean After)
1417 {
1418   if(Continuity < 0 || Continuity > 3) return;
1419   Standard_Integer Cont = Continuity;
1420
1421   // Kpart ?
1422   Handle(Geom_RectangularTrimmedSurface) TS =
1423     Handle(Geom_RectangularTrimmedSurface)::DownCast (Surface);
1424   if (ExtendKPart(TS,Length, InU, After) ) {
1425     Surface = TS;
1426     return;
1427   }
1428
1429 //  format BSplineSurface avec un degre suffisant pour la continuite voulue
1430   Handle(Geom_BSplineSurface) BS =
1431     Handle(Geom_BSplineSurface)::DownCast (Surface);
1432   if (BS.IsNull()) {
1433     //BS = GeomConvert::SurfaceToBSplineSurface(Surface);
1434     Standard_Real Tol = Precision::Confusion(); //1.e-4;
1435     GeomAbs_Shape UCont = GeomAbs_C1, VCont = GeomAbs_C1;
1436     Standard_Integer degU = 14, degV = 14;
1437     Standard_Integer nmax = 16;
1438     Standard_Integer thePrec = 1;
1439     GeomConvert_ApproxSurface theApprox(Surface,Tol,UCont,VCont,degU,degV,nmax,thePrec);
1440     if (theApprox.HasResult())
1441       BS = theApprox.Surface();
1442     else
1443       BS = GeomConvert::SurfaceToBSplineSurface(Surface);
1444   }
1445   if (InU&&(BS->UDegree()<Continuity+1))
1446     BS->IncreaseDegree(Continuity+1,BS->VDegree());
1447   if (!InU&&(BS->VDegree()<Continuity+1))
1448     BS->IncreaseDegree(BS->UDegree(),Continuity+1);
1449
1450   // si BS etait periodique dans le sens de l'extension, elle ne le sera plus
1451   if ( (InU&&(BS->IsUPeriodic())) || (!InU&&(BS->IsVPeriodic())) ) {
1452     Standard_Real U0,U1,V0,V1;
1453     BS->Bounds(U0,U1,V0,V1);
1454     BS->Segment(U0,U1,V0,V1);
1455   }
1456
1457
1458   Standard_Boolean rational = ( InU && BS->IsURational() )
1459                                   || ( !InU && BS->IsVRational() ) ;
1460   Standard_Boolean NullWeight;
1461    Standard_Real EpsW = 10*Precision::PConfusion();
1462   Standard_Integer gap = 3;
1463   if ( rational ) gap++;
1464
1465
1466
1467   Standard_Integer Cdeg, Cdim, NbP, Ksize, Psize ;
1468   Standard_Integer ii, jj, ipole, Kount;
1469   Standard_Real Tbord, lambmin=Length;
1470   Standard_Real * Padr;
1471   Standard_Boolean Ok;
1472   Handle(TColStd_HArray1OfReal)  FKnots, Point, lambda, Tgte, Poles;
1473
1474
1475
1476
1477   for (Kount=0, Ok=Standard_False; Kount<=2 && !Ok; Kount++) {
1478     //  transformation de la surface en une BSpline non rationnelle a une variable
1479     //  de degre UDegree ou VDegree et de dimension 3 ou 4 x NbVpoles ou NbUpoles
1480     //  le nombre de poles egal a NbUpoles ou NbVpoles
1481     //  ATTENTION : dans le cas rationnel, un point de coordonnees (x,y,z)
1482     //              et de poids w devient un point de coordonnees (wx, wy, wz, w )
1483
1484
1485     if (InU) {
1486       Cdeg = BS->UDegree();
1487       NbP = BS->NbUPoles();
1488       Cdim = BS->NbVPoles() * gap;
1489     }
1490     else {
1491       Cdeg = BS->VDegree();
1492       NbP = BS->NbVPoles();
1493       Cdim = BS->NbUPoles() * gap;
1494     }
1495
1496     //  les noeuds plats
1497     Ksize = NbP + Cdeg + 1;
1498     FKnots = new (TColStd_HArray1OfReal) (1,Ksize);
1499     if (InU)
1500       BS->UKnotSequence(FKnots->ChangeArray1());
1501     else
1502       BS->VKnotSequence(FKnots->ChangeArray1());
1503
1504     //  le parametre du noeud de raccord
1505     if (After)
1506       Tbord = FKnots->Value(FKnots->Upper()-Cdeg);
1507     else
1508       Tbord = FKnots->Value(FKnots->Lower()+Cdeg);
1509
1510     //  les poles
1511     Psize = Cdim * NbP;
1512     Poles = new (TColStd_HArray1OfReal) (1,Psize);
1513
1514     if (InU) {
1515       for (ii=1,ipole=1; ii<=NbP; ii++) {
1516         for (jj=1;jj<=BS->NbVPoles();jj++) {
1517           Poles->SetValue(ipole,   BS->Pole(ii,jj).X());
1518           Poles->SetValue(ipole+1, BS->Pole(ii,jj).Y());
1519           Poles->SetValue(ipole+2, BS->Pole(ii,jj).Z());
1520           if (rational) Poles->SetValue(ipole+3, BS->Weight(ii,jj));
1521           ipole+=gap;
1522         }
1523       }
1524     }
1525     else {
1526       for (jj=1,ipole=1; jj<=NbP; jj++) {
1527         for (ii=1;ii<=BS->NbUPoles();ii++) {
1528           Poles->SetValue(ipole,   BS->Pole(ii,jj).X());
1529           Poles->SetValue(ipole+1, BS->Pole(ii,jj).Y());
1530           Poles->SetValue(ipole+2, BS->Pole(ii,jj).Z());
1531           if (rational) Poles->SetValue(ipole+3, BS->Weight(ii,jj));
1532           ipole+=gap;
1533         }
1534       }
1535     }
1536     Padr = (Standard_Real *) &Poles->ChangeValue(1);
1537
1538     //  calcul du point de raccord et de la tangente
1539     Point = new (TColStd_HArray1OfReal)(1,Cdim);
1540     Tgte  = new (TColStd_HArray1OfReal)(1,Cdim);
1541     lambda = new (TColStd_HArray1OfReal)(1,Cdim);
1542
1543     Standard_Boolean  periodic_flag = Standard_False ;
1544     Standard_Integer extrap_mode[2], derivative_request = Max(Continuity,1);
1545     extrap_mode[0] = extrap_mode[1] = Cdeg;
1546     TColStd_Array1OfReal  Result(1, Cdim * (derivative_request+1)) ;
1547
1548     TColStd_Array1OfReal& tgte = Tgte->ChangeArray1();
1549     TColStd_Array1OfReal& point = Point->ChangeArray1();
1550     TColStd_Array1OfReal& lamb = lambda->ChangeArray1();
1551
1552     Standard_Real * Radr = (Standard_Real *) &Result(1) ;
1553
1554     BSplCLib::Eval(Tbord,periodic_flag,derivative_request,extrap_mode[0],
1555                    Cdeg,FKnots->Array1(),Cdim,*Padr,*Radr);
1556     Ok = Standard_True;
1557     for (ii=1;ii<=Cdim;ii++) {
1558       point(ii) = Result(ii);
1559       tgte(ii) = Result(ii+Cdim);
1560     }
1561
1562     //  calcul de la contrainte a atteindre
1563
1564     gp_Vec CurT, OldT;
1565
1566     Standard_Real NTgte, val, Tgtol = 1.e-12, OldN = 0.0;
1567     if (rational) {
1568       for (ii=gap;ii<=Cdim;ii+=gap) {
1569         tgte(ii) = 0.;
1570       }
1571       for (ii=gap;ii<=Cdim;ii+=gap) {
1572         CurT.SetCoord(tgte(ii-3),tgte(ii-2), tgte(ii-1));
1573         NTgte=CurT.Magnitude();
1574         if (NTgte>Tgtol) {
1575           val =  Length/NTgte;
1576           // Attentions aux Cas ou le segment donne par les poles
1577           // est oppose au sens de la derive
1578           // Exemple: Certaine portions de tore.
1579           if ( (OldN > Tgtol) && (CurT.Angle(OldT) > 2)) {
1580             Ok = Standard_False;
1581           }
1582
1583           lamb(ii-1) = lamb(ii-2) = lamb(ii-3) = val;
1584           lamb(ii) = 0.;
1585           lambmin = Min(lambmin, val);
1586         }
1587         else {
1588           lamb(ii-1) = lamb(ii-2) = lamb(ii-3) = 0.;
1589           lamb(ii) = 0.;
1590         }
1591         OldT = CurT;
1592         OldN = NTgte;
1593       }
1594     }
1595     else {
1596       for (ii=gap;ii<=Cdim;ii+=gap) {
1597         CurT.SetCoord(tgte(ii-2),tgte(ii-1), tgte(ii));
1598         NTgte=CurT.Magnitude();
1599         if (NTgte>Tgtol) {
1600           val =  Length/NTgte;
1601           // Attentions aux Cas ou le segment donne par les poles
1602           // est oppose au sens de la derive
1603           // Exemple: Certaine portion de tore.
1604           if ( (OldN > Tgtol) && (CurT.Angle(OldT) > 2)) {
1605              Ok = Standard_False;
1606           }
1607           lamb(ii) = lamb(ii-1) = lamb(ii-2) = val;
1608           lambmin = Min(lambmin, val);
1609         }
1610         else {
1611           lamb(ii) =lamb(ii-1) = lamb(ii-2) = 0.;
1612         }
1613         OldT = CurT;
1614         OldN = NTgte;
1615       }
1616     }
1617     if (!Ok && Kount<2) {
1618       // On augmente le degre de l'iso bord afin de rapprocher les poles de la surface
1619       // Et on ressaye
1620       if (InU) BS->IncreaseDegree(BS->UDegree(), BS->VDegree()+2);
1621       else     BS->IncreaseDegree(BS->UDegree()+2, BS->VDegree());
1622     }
1623   }
1624
1625
1626   TColStd_Array1OfReal ConstraintPoint(1,Cdim);
1627   if (After) {
1628     for (ii=1;ii<=Cdim;ii++) {
1629       ConstraintPoint(ii) = Point->Value(ii) + lambda->Value(ii)*Tgte->Value(ii);
1630     }
1631   }
1632   else {
1633     for (ii=1;ii<=Cdim;ii++) {
1634       ConstraintPoint(ii) = Point->Value(ii) - lambda->Value(ii)*Tgte->Value(ii);
1635     }
1636   }
1637
1638 //  cas particulier du rationnel
1639   if (rational) {
1640     for (ipole=1;ipole<=Psize;ipole+=gap) {
1641       Poles->ChangeValue(ipole) *= Poles->Value(ipole+3);
1642       Poles->ChangeValue(ipole+1) *= Poles->Value(ipole+3);
1643       Poles->ChangeValue(ipole+2) *= Poles->Value(ipole+3);
1644     }
1645     for (ii=1;ii<=Cdim;ii+=gap) {
1646       ConstraintPoint(ii) *= ConstraintPoint(ii+3);
1647       ConstraintPoint(ii+1) *= ConstraintPoint(ii+3);
1648       ConstraintPoint(ii+2) *= ConstraintPoint(ii+3);
1649     }
1650   }
1651
1652 //  tableaux necessaires pour l'extension
1653   Standard_Integer Ksize2 = Ksize+Cdeg, NbPoles, NbKnots;
1654   TColStd_Array1OfReal  FK(1, Ksize2) ;
1655   Standard_Real * FKRadr = &FK(1);
1656
1657   Standard_Integer Psize2 = Psize+Cdeg*Cdim;
1658   TColStd_Array1OfReal  PRes(1, Psize2) ;
1659   Standard_Real * PRadr = &PRes(1);
1660   Standard_Real ww;
1661   Standard_Boolean ExtOk = Standard_False;
1662   Handle(TColgp_HArray2OfPnt) NewPoles;
1663   Handle(TColStd_HArray2OfReal) NewWeights;
1664
1665
1666   for (Kount=1; Kount<=5 && !ExtOk; Kount++) {
1667     //  extension
1668     BSplCLib::TangExtendToConstraint(FKnots->Array1(),
1669                                      lambmin,NbP,*Padr,
1670                                      Cdim,Cdeg,
1671                                      ConstraintPoint, Cont, After,
1672                                      NbPoles, NbKnots,*FKRadr, *PRadr);
1673
1674     //  recopie des poles du resultat sous forme de points 3D et de poids
1675     Standard_Integer NU, NV, indice ;
1676     if (InU) {
1677       NU = NbPoles;
1678       NV = BS->NbVPoles();
1679     }
1680     else {
1681       NU = BS->NbUPoles();
1682       NV = NbPoles;
1683     }
1684
1685     NewPoles = new (TColgp_HArray2OfPnt)(1,NU,1,NV);
1686     TColgp_Array2OfPnt& NewP = NewPoles->ChangeArray2();
1687     NewWeights = new (TColStd_HArray2OfReal) (1,NU,1,NV);
1688     TColStd_Array2OfReal& NewW = NewWeights->ChangeArray2();
1689
1690     if (!rational) NewW.Init(1.);
1691     NullWeight= Standard_False;
1692
1693     if (InU) {
1694       for (ii=1; ii<=NU && !NullWeight; ii++) {
1695         for (jj=1; jj<=NV && !NullWeight; jj++) {
1696           indice = 1+(ii-1)*Cdim+(jj-1)*gap;
1697           NewP(ii,jj).SetCoord(1,PRes(indice));
1698           NewP(ii,jj).SetCoord(2,PRes(indice+1));
1699           NewP(ii,jj).SetCoord(3,PRes(indice+2));
1700           if (rational) {
1701             ww =  PRes(indice+3);
1702             if (ww < EpsW) {
1703               NullWeight = Standard_True;
1704             }
1705             else {
1706               NewW(ii,jj) = ww;
1707               NewP(ii,jj).ChangeCoord() /= ww;
1708             }
1709           }
1710         }
1711       }
1712     }
1713     else {
1714       for (jj=1; jj<=NV && !NullWeight; jj++) {
1715         for (ii=1; ii<=NU && !NullWeight; ii++) {
1716           indice = 1+(ii-1)*gap+(jj-1)*Cdim;
1717           NewP(ii,jj).SetCoord(1,PRes(indice));
1718           NewP(ii,jj).SetCoord(2,PRes(indice+1));
1719           NewP(ii,jj).SetCoord(3,PRes(indice+2));
1720           if (rational) {
1721             ww =  PRes(indice+3);
1722             if (ww < EpsW) {
1723               NullWeight = Standard_True;
1724             }
1725             else {
1726               NewW(ii,jj) = ww;
1727               NewP(ii,jj).ChangeCoord() /= ww;
1728             }
1729           }
1730         }
1731       }
1732     }
1733
1734     if (NullWeight) {
1735 #if DEB
1736       cout << "Echec de l'Extension rationnelle" << endl;
1737 #endif
1738       lambmin /= 3.;
1739       NullWeight = Standard_False;
1740     }
1741     else {
1742       ExtOk = Standard_True;
1743     }
1744   }
1745
1746
1747 // recopie des noeuds plats sous forme de noeuds avec leurs multiplicites
1748 // calcul des degres du resultat
1749   Standard_Integer Usize = BS->NbUKnots(), Vsize = BS->NbVKnots(), UDeg, VDeg;
1750   if (InU)
1751     Usize++;
1752   else
1753     Vsize++;
1754   TColStd_Array1OfReal UKnots(1,Usize);
1755   TColStd_Array1OfReal VKnots(1,Vsize);
1756   TColStd_Array1OfInteger UMults(1,Usize);
1757   TColStd_Array1OfInteger VMults(1,Vsize);
1758   TColStd_Array1OfReal FKRes(1, NbKnots);
1759
1760   for (ii=1; ii<=NbKnots; ii++)
1761      FKRes(ii) = FK(ii);
1762
1763   if (InU) {
1764     BSplCLib::Knots(FKRes, UKnots, UMults);
1765     UDeg = Cdeg;
1766     UMults(Usize) = UDeg+1; // Petite verrue utile quand la continuite
1767                              // n'est pas ok.
1768     BS->VKnots(VKnots);
1769     BS->VMultiplicities(VMults);
1770     VDeg = BS->VDegree();
1771   }
1772   else {
1773     BSplCLib::Knots(FKRes, VKnots, VMults);
1774     VDeg = Cdeg;
1775     VMults(Vsize) = VDeg+1;
1776     BS->UKnots(UKnots);
1777     BS->UMultiplicities(UMults);
1778     UDeg = BS->UDegree();
1779   }
1780
1781 //  construction de la surface BSpline resultat
1782   Handle(Geom_BSplineSurface) Res =
1783     new (Geom_BSplineSurface) (NewPoles->Array2(),
1784                                NewWeights->Array2(),
1785                                UKnots,VKnots,
1786                                UMults,VMults,
1787                                UDeg,VDeg,
1788                                BS->IsUPeriodic(),
1789                                BS->IsVPeriodic());
1790   Surface = Res;
1791 }
1792
1793 //=======================================================================
1794 //function : Inertia
1795 //purpose  :
1796 //=======================================================================
1797 void GeomLib::Inertia(const TColgp_Array1OfPnt& Points,
1798                       gp_Pnt& Bary,
1799                       gp_Dir& XDir,
1800                       gp_Dir& YDir,
1801                       Standard_Real& Xgap,
1802                       Standard_Real& Ygap,
1803                       Standard_Real& Zgap)
1804 {
1805   gp_XYZ GB(0., 0., 0.), Diff;
1806 //  gp_Vec A,B,C,D;
1807
1808   Standard_Integer i,nb=Points.Length();
1809   GB.SetCoord(0.,0.,0.);
1810   for (i=1; i<=nb; i++)
1811       GB += Points(i).XYZ();
1812
1813   GB /= nb;
1814
1815   math_Matrix M (1, 3, 1, 3);
1816   M.Init(0.);
1817   for (i=1; i<=nb; i++) {
1818     Diff.SetLinearForm(-1, Points(i).XYZ(), GB);
1819     M(1,1) += Diff.X() *  Diff.X();
1820     M(2,2) += Diff.Y() *  Diff.Y();
1821     M(3,3) += Diff.Z() *  Diff.Z();
1822     M(1,2) += Diff.X() *  Diff.Y();
1823     M(1,3) += Diff.X() *  Diff.Z();
1824     M(2,3) += Diff.Y() *  Diff.Z();
1825   }
1826
1827   M(2,1)=M(1,2) ;
1828   M(3,1)=M(1,3) ;
1829   M(3,2)=M(2,3) ;
1830
1831   M /= nb;
1832
1833   math_Jacobi J(M);
1834   if (!J.IsDone()) {
1835 #if DEB
1836     cout << "Erreur dans Jacobbi" << endl;
1837     M.Dump(cout);
1838 #endif
1839   }
1840
1841   Standard_Real n1,n2,n3;
1842
1843   n1=J.Value(1);
1844   n2=J.Value(2);
1845   n3=J.Value(3);
1846
1847   Standard_Real r1 = Min(Min(n1,n2),n3), r2;
1848   Standard_Integer m1, m2, m3;
1849   if (r1==n1) {
1850     m1 = 1;
1851     r2 = Min(n2,n3);
1852     if (r2==n2) {
1853       m2 = 2;
1854       m3 = 3;
1855     }
1856     else {
1857       m2 = 3;
1858       m3 = 2;
1859     }
1860   }
1861   else {
1862     if (r1==n2) {
1863       m1 = 2 ;
1864       r2 = Min(n1,n3);
1865       if (r2==n1) {
1866         m2 = 1;
1867         m3 = 3;
1868       }
1869       else {
1870         m2 = 3;
1871         m3 = 1;
1872       }
1873     }
1874     else {
1875       m1 = 3 ;
1876       r2 = Min(n1,n2);
1877       if (r2==n1) {
1878         m2 = 1;
1879         m3 = 2;
1880       }
1881       else {
1882         m2 = 2;
1883         m3 = 1;
1884       }
1885     }
1886   }
1887
1888   math_Vector V2(1,3),V3(1,3);
1889   J.Vector(m2,V2);
1890   J.Vector(m3,V3);
1891
1892   Bary.SetXYZ(GB);
1893   XDir.SetCoord(V3(1),V3(2),V3(3));
1894   YDir.SetCoord(V2(1),V2(2),V2(3));
1895
1896   Zgap = sqrt(Abs(J.Value(m1)));
1897   Ygap = sqrt(Abs(J.Value(m2)));
1898   Xgap = sqrt(Abs(J.Value(m3)));
1899 }
1900 //=======================================================================
1901 //function : AxeOfInertia
1902 //purpose  :
1903 //=======================================================================
1904 void GeomLib::AxeOfInertia(const TColgp_Array1OfPnt& Points,
1905                            gp_Ax2& Axe,
1906                            Standard_Boolean& IsSingular,
1907                            const Standard_Real Tol)
1908 {
1909   gp_Pnt Bary;
1910   gp_Dir OX,OY,OZ;
1911   Standard_Real gx, gy, gz;
1912
1913   GeomLib::Inertia(Points, Bary, OX, OY, gx, gy, gz);
1914
1915   if (gy*Points.Length()<=Tol) {
1916     gp_Ax2 axe (Bary, OX);
1917     OY = axe.XDirection();
1918     IsSingular = Standard_True;
1919   }
1920   else {
1921     IsSingular = Standard_False;
1922   }
1923
1924   OZ = OX^OY;
1925   gp_Ax2 TheAxe(Bary, OZ, OX);
1926   Axe = TheAxe;
1927 }
1928
1929 //=======================================================================
1930 //function : CanBeTreated
1931 //purpose  : indicates if the surface can be treated(if the conditions are
1932 //           filled) and need to be treated(if the surface hasn't been yet
1933 //           treated or if the surface is rationnal and non periodic)
1934 //=======================================================================
1935
1936 static Standard_Boolean CanBeTreated(Handle(Geom_BSplineSurface)& BSurf)
1937
1938 {Standard_Integer i;
1939  Standard_Real    lambda;                                    //proportionnality coefficient
1940  Standard_Boolean AlreadyTreated=Standard_True;
1941
1942  if (!BSurf->IsURational()||(BSurf->IsUPeriodic()))
1943    return Standard_False;
1944  else {
1945    lambda=(BSurf->Weight(1,1)/BSurf->Weight(BSurf->NbUPoles(),1));
1946    for (i=1;i<=BSurf->NbVPoles();i++)      //test of the proportionnality of the denominator on the boundaries
1947      if ((BSurf->Weight(1,i)/(lambda*BSurf->Weight(BSurf->NbUPoles(),i))<(1-Precision::Confusion()))||
1948          (BSurf->Weight(1,i)/(lambda*BSurf->Weight(BSurf->NbUPoles(),i))>(1+Precision::Confusion())))
1949        return Standard_False;
1950    i=1;
1951    while ((AlreadyTreated) && (i<=BSurf->NbVPoles())){        //tests if the surface has already been treated
1952      if (((BSurf->Weight(1,i)/(BSurf->Weight(2,i)))<(1-Precision::Confusion()))||
1953          ((BSurf->Weight(1,i)/(BSurf->Weight(2,i)))>(1+Precision::Confusion()))||
1954          ((BSurf->Weight(BSurf->NbUPoles()-1,i)/(BSurf->Weight(BSurf->NbUPoles(),i)))<(1-Precision::Confusion()))||
1955          ((BSurf->Weight(BSurf->NbUPoles()-1,i)/(BSurf->Weight(BSurf->NbUPoles(),i)))>(1+Precision::Confusion())))
1956        AlreadyTreated=Standard_False;
1957      i++;
1958    }
1959    if (AlreadyTreated)
1960      return Standard_False;
1961  }
1962  return Standard_True;
1963 }
1964
1965 //=======================================================================
1966 //function : law_evaluator
1967 //purpose  : usefull to estimate the value of a function of 2 variables
1968 //=======================================================================
1969
1970 static GeomLib_DenominatorMultiplierPtr MyPtr = NULL ;
1971
1972
1973 static void law_evaluator(const Standard_Integer  DerivativeRequest,
1974                           const Standard_Real     UParameter,
1975                           const Standard_Real     VParameter,
1976                           Standard_Real &         Result,
1977                           Standard_Integer &      ErrorCode) {
1978
1979   ErrorCode = 0 ;
1980
1981   if ((!(MyPtr == NULL)) &&
1982       (DerivativeRequest == 0)) {
1983     Result=MyPtr->Value(UParameter,VParameter);
1984   }
1985   else {
1986     ErrorCode = 1 ;
1987   }
1988 }
1989 //=======================================================================
1990 //function : CheckIfKnotExists
1991 //purpose  : true if the knot already exists in the knot sequence
1992 //=======================================================================
1993
1994 static Standard_Boolean CheckIfKnotExists(const TColStd_Array1OfReal&           surface_knots,
1995                                           const Standard_Real                   knot)
1996
1997 {Standard_Integer    i;
1998  for (i=1;i<=surface_knots.Length();i++)
1999    if ((surface_knots(i)-Precision::Confusion()<=knot)&&(surface_knots(i)+Precision::Confusion()>=knot))
2000      return Standard_True;
2001  return Standard_False;
2002 }
2003
2004 //=======================================================================
2005 //function : AddAKnot
2006 //purpose  : add a knot and its multiplicity to the knot sequence. This knot
2007 //           will be C2 and the degree is increased of deltasurface_degree
2008 //=======================================================================
2009
2010 static void AddAKnot(const TColStd_Array1OfReal&           knots,
2011                      const TColStd_Array1OfInteger&        mults,
2012                      const Standard_Real                   knotinserted,
2013                      const Standard_Integer                deltasurface_degree,
2014                      const Standard_Integer                finalsurfacedegree,
2015                      Handle(TColStd_HArray1OfReal) &       newknots,
2016                      Handle(TColStd_HArray1OfInteger) &    newmults)
2017
2018 {Standard_Integer      i;
2019
2020  newknots=new TColStd_HArray1OfReal(1,knots.Length()+1);
2021  newmults=new TColStd_HArray1OfInteger(1,knots.Length()+1);
2022  i=1;
2023  while (knots(i)<knotinserted){
2024    newknots->SetValue(i,knots(i));
2025    newmults->SetValue(i,mults(i)+deltasurface_degree);
2026    i++;
2027  }
2028  newknots->SetValue(i,knotinserted);                        //insertion of the new knot
2029  newmults->SetValue(i,finalsurfacedegree-2);
2030  i++;
2031  while (i<=newknots->Length()){
2032    newknots->SetValue(i,knots(i-1));
2033    newmults->SetValue(i,mults(i-1)+deltasurface_degree);
2034    i++;
2035  }
2036 }
2037
2038 //=======================================================================
2039 //function : Sort
2040 //purpose  : give the new flat knots(u or v) of the surface
2041 //=======================================================================
2042
2043 static void BuildFlatKnot(const TColStd_Array1OfReal&           surface_knots,
2044                  const TColStd_Array1OfInteger&        surface_mults,
2045                  const Standard_Integer                deltasurface_degree,
2046                  const Standard_Integer                finalsurface_degree,
2047                  const Standard_Real                   knotmin,
2048                  const Standard_Real                   knotmax,
2049                  Handle(TColStd_HArray1OfReal)&        ResultKnots,
2050                  Handle(TColStd_HArray1OfInteger)&     ResultMults)
2051
2052 {
2053   Standard_Integer  i;
2054
2055  if (CheckIfKnotExists(surface_knots,knotmin) &&
2056      CheckIfKnotExists(surface_knots,knotmax)){
2057    ResultKnots=new TColStd_HArray1OfReal(1,surface_knots.Length());
2058    ResultMults=new TColStd_HArray1OfInteger(1,surface_knots.Length());
2059    for (i=1;i<=surface_knots.Length();i++){
2060      ResultKnots->SetValue(i,surface_knots(i));
2061      ResultMults->SetValue(i,surface_mults(i)+deltasurface_degree);
2062    }
2063  }
2064  else{
2065    if ((CheckIfKnotExists(surface_knots,knotmin))&&(!CheckIfKnotExists(surface_knots,knotmax)))
2066      AddAKnot(surface_knots,surface_mults,knotmax,deltasurface_degree,finalsurface_degree,ResultKnots,ResultMults);
2067    else{
2068      if ((!CheckIfKnotExists(surface_knots,knotmin))&&(CheckIfKnotExists(surface_knots,knotmax)))
2069        AddAKnot(surface_knots,surface_mults,knotmin,deltasurface_degree,finalsurface_degree,ResultKnots,ResultMults);
2070      else{
2071        if ((!CheckIfKnotExists(surface_knots,knotmin))&&(!CheckIfKnotExists(surface_knots,knotmax))&&
2072            (knotmin==knotmax)){
2073          AddAKnot(surface_knots,surface_mults,knotmin,deltasurface_degree,finalsurface_degree,ResultKnots,ResultMults);
2074        }
2075        else{
2076          Handle(TColStd_HArray1OfReal)      IntermedKnots;
2077          Handle(TColStd_HArray1OfInteger)   IntermedMults;
2078          AddAKnot(surface_knots,surface_mults,knotmin,deltasurface_degree,finalsurface_degree,IntermedKnots,IntermedMults);
2079          AddAKnot(IntermedKnots->ChangeArray1(),IntermedMults->ChangeArray1(),knotmax,0,finalsurface_degree,ResultKnots,ResultMults);
2080        }
2081      }
2082    }
2083  }
2084 }
2085
2086 //=======================================================================
2087 //function : FunctionMultiply
2088 //purpose  : multiply the surface BSurf by a(u,v) (law_evaluator) on its
2089 //           numerator and denominator
2090 //=======================================================================
2091
2092 static void FunctionMultiply(Handle(Geom_BSplineSurface)&          BSurf,
2093                              const Standard_Real                   knotmin,
2094                              const Standard_Real                   knotmax)
2095
2096 {TColStd_Array1OfReal      surface_u_knots(1,BSurf->NbUKnots()) ;
2097  TColStd_Array1OfInteger   surface_u_mults(1,BSurf->NbUKnots()) ;
2098  TColStd_Array1OfReal      surface_v_knots(1,BSurf->NbVKnots()) ;
2099  TColStd_Array1OfInteger   surface_v_mults(1,BSurf->NbVKnots()) ;
2100  TColgp_Array2OfPnt        surface_poles(1,BSurf->NbUPoles(),
2101                                          1,BSurf->NbVPoles()) ;
2102  TColStd_Array2OfReal      surface_weights(1,BSurf->NbUPoles(),
2103                                            1,BSurf->NbVPoles()) ;
2104  Standard_Integer          i,j,k,status,new_num_u_poles,new_num_v_poles,length=0;
2105  Handle(TColStd_HArray1OfReal)     newuknots,newvknots;
2106  Handle(TColStd_HArray1OfInteger)  newumults,newvmults;
2107
2108  BSurf->UKnots(surface_u_knots) ;
2109  BSurf->UMultiplicities(surface_u_mults) ;
2110  BSurf->VKnots(surface_v_knots) ;
2111  BSurf->VMultiplicities(surface_v_mults) ;
2112  BSurf->Poles(surface_poles) ;
2113  BSurf->Weights(surface_weights) ;
2114
2115  TColStd_Array1OfReal    Knots(1,2);
2116  TColStd_Array1OfInteger Mults(1,2);
2117  Handle(TColStd_HArray1OfReal)      NewKnots;
2118  Handle(TColStd_HArray1OfInteger)   NewMults;
2119
2120  Knots(1)=0;
2121  Knots(2)=1;
2122  Mults(1)=4;
2123  Mults(2)=4;
2124  BuildFlatKnot(Knots,Mults,0,3,knotmin,knotmax,NewKnots,NewMults);
2125
2126  for (i=1;i<=NewMults->Length();i++)
2127    length+=NewMults->Value(i);
2128  TColStd_Array1OfReal       FlatKnots(1,length);
2129  BSplCLib::KnotSequence(NewKnots->ChangeArray1(),NewMults->ChangeArray1(),FlatKnots);
2130
2131  GeomLib_DenominatorMultiplier          local_denominator(BSurf,FlatKnots) ;
2132  MyPtr = &local_denominator ;                 //definition of a(u,v)
2133
2134  BuildFlatKnot(surface_u_knots,
2135                surface_u_mults,
2136                3,
2137                BSurf->UDegree()+3,
2138                knotmin,
2139                knotmax,
2140                newuknots,
2141                newumults);
2142  BuildFlatKnot(surface_v_knots,
2143                surface_v_mults,
2144                BSurf->VDegree(),
2145                2*(BSurf->VDegree()),
2146                1.0,
2147                0.0,
2148                newvknots,
2149                newvmults);
2150  length=0;
2151  for (i=1;i<=newumults->Length();i++)
2152    length+=newumults->Value(i);
2153  new_num_u_poles=(length-BSurf->UDegree()-3-1);
2154  TColStd_Array1OfReal       newuflatknots(1,length);
2155  length=0;
2156  for (i=1;i<=newvmults->Length();i++)
2157    length+=newvmults->Value(i);
2158  new_num_v_poles=(length-2*BSurf->VDegree()-1);
2159  TColStd_Array1OfReal       newvflatknots(1,length);
2160
2161  TColgp_Array2OfPnt        NewNumerator(1,new_num_u_poles,1,new_num_v_poles);
2162  TColStd_Array2OfReal      NewDenominator(1,new_num_u_poles,1,new_num_v_poles);
2163
2164  BSplCLib::KnotSequence(newuknots->ChangeArray1(),newumults->ChangeArray1(),newuflatknots);
2165  BSplCLib::KnotSequence(newvknots->ChangeArray1(),newvmults->ChangeArray1(),newvflatknots);
2166 //POP pour WNT
2167  BSplSLib_EvaluatorFunction ev = law_evaluator;
2168 // BSplSLib::FunctionMultiply(law_evaluator,               //multiplication
2169  BSplSLib::FunctionMultiply(ev,               //multiplication
2170                             BSurf->UDegree(),
2171                             BSurf->VDegree(),
2172                             surface_u_knots,
2173                             surface_v_knots,
2174                             surface_u_mults,
2175                             surface_v_mults,
2176                             surface_poles,
2177                             surface_weights,
2178                             newuflatknots,
2179                             newvflatknots,
2180                             BSurf->UDegree()+3,
2181                             2*(BSurf->VDegree()),
2182                             NewNumerator,
2183                             NewDenominator,
2184                             status);
2185  if (status!=0)
2186    Standard_ConstructionError::Raise("GeomLib Multiplication Error") ;
2187  for (i = 1 ; i <= new_num_u_poles ; i++) {
2188       for (j = 1 ; j <= new_num_v_poles ; j++) {
2189         for (k = 1 ; k <= 3 ; k++) {
2190           NewNumerator(i,j).SetCoord(k,NewNumerator(i,j).Coord(k)/NewDenominator(i,j)) ;
2191         }
2192       }
2193     }
2194  BSurf= new Geom_BSplineSurface(NewNumerator,
2195                                 NewDenominator,
2196                                 newuknots->ChangeArray1(),
2197                                 newvknots->ChangeArray1(),
2198                                 newumults->ChangeArray1(),
2199                                 newvmults->ChangeArray1(),
2200                                 BSurf->UDegree()+3,
2201                                 2*(BSurf->VDegree()) );
2202 }
2203
2204 //=======================================================================
2205 //function : CancelDenominatorDerivative1D
2206 //purpose  : cancel the denominator derivative in one direction
2207 //=======================================================================
2208
2209 static void CancelDenominatorDerivative1D(Handle(Geom_BSplineSurface) & BSurf)
2210
2211 {Standard_Integer            i,j;
2212  Standard_Real               uknotmin=1.0,uknotmax=0.0,
2213                              x,y,
2214                              startu_value,
2215                              endu_value;
2216  TColStd_Array1OfReal        BSurf_u_knots(1,BSurf->NbUKnots()) ;
2217
2218  startu_value=BSurf->UKnot(1);
2219  endu_value=BSurf->UKnot(BSurf->NbUKnots());
2220  BSurf->UKnots(BSurf_u_knots) ;
2221  BSplCLib::Reparametrize(0.0,1.0,BSurf_u_knots);
2222  BSurf->SetUKnots(BSurf_u_knots);                             //reparametrisation of the surface
2223  Handle(Geom_BSplineCurve) BCurve;
2224  TColStd_Array1OfReal      BCurveWeights(1,BSurf->NbUPoles());
2225  TColgp_Array1OfPnt        BCurvePoles(1,BSurf->NbUPoles());
2226  TColStd_Array1OfReal      BCurveKnots(1,BSurf->NbUKnots());
2227  TColStd_Array1OfInteger   BCurveMults(1,BSurf->NbUKnots());
2228
2229  if (CanBeTreated(BSurf)){
2230    for (i=1;i<=BSurf->NbVPoles();i++){  //loop on each pole function
2231      x=1.0;y=0.0;
2232      for (j=1;j<=BSurf->NbUPoles();j++){
2233        BCurveWeights(j)=BSurf->Weight(j,i);
2234        BCurvePoles(j)=BSurf->Pole(j,i);
2235      }
2236      BSurf->UKnots(BCurveKnots);
2237      BSurf->UMultiplicities(BCurveMults);
2238      BCurve = new Geom_BSplineCurve(BCurvePoles, //building of a pole function
2239                                     BCurveWeights,
2240                                     BCurveKnots,
2241                                     BCurveMults,
2242                                     BSurf->UDegree());
2243      Hermit::Solutionbis(BCurve,x,y,Precision::Confusion(),Precision::Confusion());
2244      if (x<uknotmin)
2245        uknotmin=x;    //uknotmin,uknotmax:extremal knots
2246      if ((x!=1.0)&&(x>uknotmax))
2247        uknotmax=x;
2248      if ((y!=0.0)&&(y<uknotmin))
2249        uknotmin=y;
2250      if (y>uknotmax)
2251        uknotmax=y;
2252    }
2253
2254    FunctionMultiply(BSurf,uknotmin,uknotmax);                 //multiplication
2255
2256    BSurf->UKnots(BSurf_u_knots) ;
2257    BSplCLib::Reparametrize(startu_value,endu_value,BSurf_u_knots);
2258    BSurf->SetUKnots(BSurf_u_knots);
2259  }
2260 }
2261
2262 //=======================================================================
2263 //function : CancelDenominatorDerivative
2264 //purpose  :
2265 //=======================================================================
2266
2267 void GeomLib::CancelDenominatorDerivative(Handle(Geom_BSplineSurface)         & BSurf,
2268                                           const Standard_Boolean              udirection,
2269                                           const Standard_Boolean              vdirection)
2270
2271 {if (udirection && !vdirection)
2272    CancelDenominatorDerivative1D(BSurf);
2273  else{
2274    if (!udirection && vdirection) {
2275      BSurf->ExchangeUV();
2276      CancelDenominatorDerivative1D(BSurf);
2277      BSurf->ExchangeUV();
2278    }
2279    else{
2280      if (udirection && vdirection){                            //optimize the treatment
2281        if (BSurf->UDegree()<=BSurf->VDegree()){
2282          CancelDenominatorDerivative1D(BSurf);
2283          BSurf->ExchangeUV();
2284          CancelDenominatorDerivative1D(BSurf);
2285          BSurf->ExchangeUV();
2286        }
2287        else{
2288          BSurf->ExchangeUV();
2289          CancelDenominatorDerivative1D(BSurf);
2290          BSurf->ExchangeUV();
2291          CancelDenominatorDerivative1D(BSurf);
2292        }
2293      }
2294    }
2295  }
2296 }
2297
2298 //=======================================================================
2299 //function : NormEstim
2300 //purpose  :
2301 //=======================================================================
2302
2303 Standard_Integer GeomLib::NormEstim(const Handle(Geom_Surface)& S,
2304                                     const gp_Pnt2d& UV,
2305                                     const Standard_Real Tol, gp_Dir& N)
2306 {
2307   gp_Vec DU, DV;
2308   gp_Pnt DummyPnt;
2309   Standard_Real aTol2 = Square(Tol);
2310
2311   S->D1(UV.X(), UV.Y(), DummyPnt, DU, DV);
2312
2313   Standard_Real MDU = DU.SquareMagnitude(), MDV = DV.SquareMagnitude();
2314
2315   Standard_Real h, sign;
2316   Standard_Boolean AlongV;
2317   Handle(Geom_Curve) Iso;
2318   Standard_Real t, first, last, bid1, bid2;
2319   gp_Vec Tang;
2320
2321   if(MDU >= aTol2 && MDV >= aTol2) {
2322     gp_Vec Norm = DU^DV;
2323     Standard_Real Magn = Norm.SquareMagnitude();
2324     if(Magn < aTol2) return 3;
2325
2326     //Magn = sqrt(Magn);
2327     N.SetXYZ(Norm.XYZ());
2328
2329     return 0;
2330   }
2331   else if(MDU < aTol2 && MDV >= aTol2) {
2332     AlongV = Standard_True;
2333     Iso = S->UIso(UV.X());
2334     t = UV.Y();
2335     S->Bounds(bid1, bid2, first, last);
2336   }
2337   else if(MDU >= aTol2 && MDV < aTol2) {
2338     AlongV = Standard_False;
2339     Iso = S->VIso(UV.Y());
2340     t = UV.X();
2341     S->Bounds(first, last, bid1, bid2);
2342   }
2343   else {
2344     return 3;
2345   }
2346
2347   Standard_Real L = .001;
2348
2349   if(Precision::IsInfinite(Abs(first))) first = t - 1.;
2350   if(Precision::IsInfinite(Abs(last)))  last  = t + 1.;
2351
2352   if(last - t >= t - first) {
2353     sign = 1.;
2354   }
2355   else {
2356     sign = -1.;
2357   }
2358
2359   Standard_Real hmax = .01*(last - first);
2360   if(AlongV) {
2361     h = Min(L/sqrt(MDV), hmax);
2362     S->D1(UV.X(), UV.Y() + sign*h, DummyPnt, DU, DV);
2363   }
2364   else {
2365     h = Min(L/sqrt(MDU), hmax);
2366     S->D1(UV.X() + sign*h, UV.Y(), DummyPnt, DU, DV);
2367   }
2368
2369   gp_Vec DD;
2370
2371   gp_Vec NAux = DU^DV;
2372   Standard_Real h1 = h;
2373   while(NAux.SquareMagnitude() < aTol2) {
2374     h1 += h;
2375     if(AlongV) {
2376       Standard_Real v = UV.Y() + sign*h1;
2377
2378       if(v < first || v > last) return 3;
2379
2380       S->D1(UV.X(), v, DummyPnt, DU, DV);
2381     }
2382     else {
2383       Standard_Real v = UV.X() + sign*h1;
2384
2385       if(v < first || v > last) return 3;
2386
2387       S->D1(v, UV.Y(), DummyPnt, DU, DV);
2388
2389     }
2390     NAux = DU^DV;
2391   }
2392
2393
2394   Iso->D2(t, DummyPnt, Tang, DD);
2395
2396   if(DD.SquareMagnitude() >= aTol2) {
2397     gp_Vec NV = DD * (Tang * Tang) - Tang * (Tang * DD);
2398     Standard_Real MagnNV = NV.SquareMagnitude();
2399     if(MagnNV < aTol2) return 3;
2400
2401     MagnNV = sqrt(MagnNV);
2402     N.SetXYZ(NV.XYZ()/MagnNV);
2403
2404     Standard_Real par = .5*(bid2-bid1);
2405
2406     if(AlongV) {
2407       Iso = S->UIso(par);
2408     }
2409     else {
2410       Iso = S->VIso(par);
2411     }
2412
2413     Iso->D2(t, DummyPnt, Tang, DD);
2414
2415     gp_Vec N1V = DD * (Tang * Tang) - Tang * (Tang * DD);
2416     Standard_Real MagnN1V = N1V.SquareMagnitude();
2417     if(MagnN1V < aTol2) return 3;
2418
2419     MagnN1V = sqrt(MagnN1V);
2420     gp_Dir N1(N1V.XYZ()/MagnN1V);
2421
2422     Standard_Integer res = 1;
2423
2424     if(N*N1 < 1. - Tol) res = 2;
2425
2426     if(N*NAux <= 0.) N.Reverse();
2427
2428     return res;
2429   }
2430   else {
2431     //Seems to be conical singular point
2432
2433     if(AlongV) {
2434       NAux = DU^Tang;
2435     }
2436     else {
2437       NAux = Tang^DV;
2438     }
2439
2440     sign = NAux.Magnitude();
2441
2442     if(sign < Tol) return 3;
2443
2444     N = NAux;
2445
2446     return 2;
2447
2448   }
2449
2450 }
2451
2452