src/GeomLib/GeomLib.cxx

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