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