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