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