OCC22579 Improving thread-safety of GeomFill
[occt.git] / src / GeomFill / GeomFill.cxx
1 // File:        GeomFill.cxx
2 // Created:     Fri Feb 25 09:49:20 1994
3 // Author:      Bruno DUMORTIER
4 //              <dub@fuegox>
5
6
7 #include <GeomFill.ixx>
8 #include <GeomFill_Generator.hxx>
9
10 #include <Geom_RectangularTrimmedSurface.hxx>
11 #include <Geom_CylindricalSurface.hxx>
12 #include <Geom_ConicalSurface.hxx>
13 #include <Geom_Plane.hxx>
14 #include <Geom_TrimmedCurve.hxx>
15 #include <Geom_BSplineCurve.hxx>
16 #include <Geom_Line.hxx>
17 #include <Geom_Circle.hxx>
18
19 #include <gp_Lin.hxx>
20 #include <gp_Circ.hxx>
21 #include <gp_Dir.hxx>
22 #include <gp_Ax3.hxx>
23 #include <gp_Vec.hxx>
24
25 #include <GeomConvert.hxx>
26 #include <GeomFill_PolynomialConvertor.hxx>
27 #include <GeomFill_QuasiAngularConvertor.hxx>
28 #include <Precision.hxx>
29
30
31 //=======================================================================
32 //function : Surface
33 //purpose  : 
34 //=======================================================================
35 Handle(Geom_Surface) GeomFill::Surface
36   (const Handle(Geom_Curve)& Curve1, 
37    const Handle(Geom_Curve)& Curve2)
38
39 {
40   Handle(Geom_Curve) TheCurve1, TheCurve2;
41   Handle(Geom_Surface) Surf;
42   
43   // recherche du type de la surface resultat:
44   // les surfaces reglees particulieres sont :
45   //     - les plans
46   //     - les cylindres
47   //     - les cones
48   // dans ces trois cas les courbes doivent etre de meme type :
49   //     - ou 2 droites
50   //     - ou 2 cercles
51   
52   
53   Standard_Real a1=0, a2=0, b1=0, b2=0;
54   Standard_Boolean Trim1= Standard_False, Trim2 = Standard_False;
55   if ( Curve1->IsKind(STANDARD_TYPE(Geom_TrimmedCurve))) {
56     Handle(Geom_TrimmedCurve) Ctrim 
57       = Handle(Geom_TrimmedCurve)::DownCast(Curve1);
58     TheCurve1 = Ctrim->BasisCurve();
59     a1 = Ctrim->FirstParameter();
60     b1 = Ctrim->LastParameter();
61     Trim1 = Standard_True;
62   }
63   else {
64     TheCurve1 = Handle(Geom_Curve)::DownCast(Curve1->Copy());
65   }
66   if ( Curve2->IsKind(STANDARD_TYPE(Geom_TrimmedCurve))) {
67     Handle(Geom_TrimmedCurve) Ctrim 
68       = Handle(Geom_TrimmedCurve)::DownCast(Curve2);
69     TheCurve2 = Ctrim->BasisCurve();
70     a2 = Ctrim->FirstParameter();
71     b2 = Ctrim->LastParameter();
72     Trim2 = Standard_True;
73   }
74   else {
75     TheCurve2 = Handle(Geom_Curve)::DownCast(Curve2->Copy());
76   }
77
78   Standard_Boolean IsDone = Standard_False;
79   // Les deux courbes sont des droites.
80   if ( TheCurve1->IsKind(STANDARD_TYPE(Geom_Line)) &&
81        TheCurve2->IsKind(STANDARD_TYPE(Geom_Line)) &&
82        Trim1 && Trim2                     ) {
83
84     gp_Lin L1 = (Handle(Geom_Line)::DownCast(TheCurve1))->Lin();
85     gp_Lin L2 = (Handle(Geom_Line)::DownCast(TheCurve2))->Lin();
86     gp_Dir D1 = L1.Direction();
87     gp_Dir D2 = L2.Direction();
88     
89     if ( D1.IsParallel(D2, Precision::Angular())) {
90       gp_Vec P1P2(L1.Location(),L2.Location()); 
91       Standard_Real proj = P1P2.Dot(D1);
92
93       if ( D1.IsEqual(D2, Precision::Angular())) {
94         if ( Abs( a1 - proj - a2 ) <= Precision::Confusion() &&
95              Abs( b1 - proj - b2 ) <= Precision::Confusion()    ) {
96           gp_Ax3 Ax(L1.Location(), gp_Dir(D1.Crossed(P1P2)),D1);
97           Handle(Geom_Plane) P = new Geom_Plane(Ax);
98           Standard_Real V = P1P2.Dot( Ax.YDirection());
99           Surf = new Geom_RectangularTrimmedSurface( P , a1, b1,
100                                                     Min(0.,V),Max(0.,V));
101           IsDone = Standard_True;
102         }
103       }
104       if ( D1.IsOpposite(D2, Precision::Angular())) {
105         if ( Abs( a1 - proj + b2 ) <= Precision::Confusion() &&
106              Abs( b1 - proj + a2 ) <= Precision::Confusion()    ) {
107           gp_Ax3 Ax(L1.Location(), gp_Dir(D1.Crossed(P1P2)),D1);
108           Handle(Geom_Plane) P = new Geom_Plane(Ax);
109           Standard_Real V = P1P2.Dot( Ax.YDirection());
110           Surf = new Geom_RectangularTrimmedSurface( P , a1, b1,
111                                                     Min(0.,V),Max(0.,V));
112           IsDone = Standard_True;
113         }
114       }
115     }
116   }
117   
118   
119   // Les deux courbes sont des cercles.
120   else if ( TheCurve1->IsKind(STANDARD_TYPE(Geom_Circle)) &&
121             TheCurve2->IsKind(STANDARD_TYPE(Geom_Circle))   ) {
122
123     gp_Circ C1 = (Handle(Geom_Circle)::DownCast(TheCurve1))->Circ();
124     gp_Circ C2 = (Handle(Geom_Circle)::DownCast(TheCurve2))->Circ();
125
126     gp_Ax3 A1 = C1.Position();
127     gp_Ax3 A2 = C2.Position();
128
129     // first, A1 & A2 must be coaxials
130     if ( A1.Axis().IsCoaxial(A2.Axis(),Precision::Angular(),
131                                        Precision::Confusion()) ) {
132       Standard_Real V = 
133         gp_Vec( A1.Location(),A2.Location()).Dot(gp_Vec(A1.Direction()));
134       if ( !Trim1 && !Trim2) {
135         if ( Abs( C1.Radius() - C2.Radius()) < Precision::Confusion()) {
136           Handle(Geom_CylindricalSurface) C = 
137             new Geom_CylindricalSurface( A1, C1.Radius());
138           Surf = new Geom_RectangularTrimmedSurface
139             ( C, Min(0.,V), Max(0.,V), Standard_False);
140         }
141         else {
142           Standard_Real Rad = C2.Radius() - C1.Radius();
143           Standard_Real Ang = ATan( Rad / V);
144           if ( Ang < 0.) {
145             A1.ZReverse();
146             V = -V;
147             Ang = -Ang;
148           }
149           Handle(Geom_ConicalSurface) C = 
150             new Geom_ConicalSurface( A1, Ang, C1.Radius());
151           V /= Cos(Ang);
152           Surf = new Geom_RectangularTrimmedSurface
153             ( C, Min(0.,V), Max(0.,V), Standard_False);
154         }
155         IsDone = Standard_True;
156       }
157       else if ( Trim1 && Trim2) {
158
159
160       }
161
162     }
163
164   }
165
166   if ( !IsDone) {
167     GeomFill_Generator Generator;
168     Generator.AddCurve(Curve1);
169     Generator.AddCurve(Curve2);
170     Generator.Perform(Precision::PConfusion());
171     Surf = Generator.Surface();
172   }
173   
174   return Surf;
175 }
176
177 //=======================================================================
178 //function : GetShape
179 //purpose  : 
180 //=======================================================================
181
182 void GeomFill::GetShape (const Standard_Real MaxAng,
183                          Standard_Integer& NbPoles,
184                          Standard_Integer& NbKnots,
185                          Standard_Integer& Degree,
186                          Convert_ParameterisationType& TConv)
187 {
188   switch (TConv) {
189   case Convert_QuasiAngular:
190     {
191       NbPoles = 7 ;
192       NbKnots = 2 ;
193       Degree  = 6 ;
194     }
195     break;
196   case Convert_Polynomial:
197     {
198       NbPoles = 8;
199       NbKnots = 2;
200       Degree = 7;
201     }
202     break;
203   default:
204     {
205       Standard_Integer NbSpan =
206         (Standard_Integer)(Ceiling(3.*Abs(MaxAng)/2./PI));
207       NbPoles = 2*NbSpan+1;
208       NbKnots = NbSpan+1;
209       Degree = 2;
210       if (NbSpan == 1) {
211         TConv = Convert_TgtThetaOver2_1;
212       }
213       else if (NbSpan == 2) {
214         TConv = Convert_TgtThetaOver2_2;
215       }
216       else if (NbSpan == 3) {
217         TConv = Convert_TgtThetaOver2_3;
218       }
219     }
220   }
221 }
222
223 //=======================================================================
224 //function : GetMinimalWeights
225 //purpose  : On suppose les extremum de poids sont obtenus pour les
226 //           extremums d'angles (A verifier eventuelement pour Quasi-Angular)
227 //=======================================================================
228
229 void GeomFill::GetMinimalWeights(const Convert_ParameterisationType  TConv,
230                                  const Standard_Real MinAng,
231                                  const Standard_Real MaxAng,
232                                  TColStd_Array1OfReal& Weights)
233
234 {
235   if (TConv == Convert_Polynomial) Weights.Init(1);
236   else {
237     gp_Ax2 popAx2(gp_Pnt(0, 0, 0), gp_Dir(0,0,1));
238     gp_Circ C (popAx2, 1);
239     Handle(Geom_TrimmedCurve) Sect1 = 
240       new Geom_TrimmedCurve(new Geom_Circle(C), 0., MaxAng);
241     Handle(Geom_BSplineCurve) CtoBspl = 
242       GeomConvert::CurveToBSplineCurve(Sect1, TConv);
243     CtoBspl->Weights(Weights);
244       
245     TColStd_Array1OfReal poids (Weights.Lower(),  Weights.Upper());
246     Standard_Real angle_min = Max(Precision::PConfusion(), MinAng);
247    
248     Handle(Geom_TrimmedCurve) Sect2 = 
249       new Geom_TrimmedCurve(new Geom_Circle(C), 0., angle_min);
250     CtoBspl = GeomConvert::CurveToBSplineCurve(Sect2, TConv);
251     CtoBspl->Weights(poids);
252
253     for (Standard_Integer ii=Weights.Lower(); ii<=Weights.Upper(); ii++) {
254       if (poids(ii) < Weights(ii)) {
255         Weights(ii) = poids(ii);
256       }
257     }
258   }
259 }
260
261
262 //=======================================================================
263 //function : Knots
264 //purpose  : 
265 //=======================================================================
266
267 void GeomFill::Knots(const Convert_ParameterisationType  TConv,
268                      TColStd_Array1OfReal& TKnots)
269 {
270   if ((TConv!=Convert_QuasiAngular) && 
271       (TConv!=Convert_Polynomial) ) {
272     Standard_Integer i;
273     Standard_Real val = 0.;
274     for (i=TKnots.Lower(); i<=TKnots.Upper(); i++) {
275       TKnots(i) = val;
276       val = val+1.;
277     }
278   }
279   else {
280     TKnots(1) = 0.;
281     TKnots(2) = 1.;
282   }
283 }
284
285
286 //=======================================================================
287 //function : Mults
288 //purpose  : 
289 //=======================================================================
290
291 void GeomFill::Mults(const Convert_ParameterisationType  TConv,
292                      TColStd_Array1OfInteger& TMults)
293 {
294   switch (TConv) {
295   case Convert_QuasiAngular :  
296     {
297       TMults(1) = 7;
298       TMults(2) = 7;
299     }
300     break;
301   case Convert_Polynomial :
302     {
303       TMults(1) = 8;
304       TMults(2) = 8;
305     }
306     break;
307
308   default : 
309     {
310       // Cas rational classsique
311       Standard_Integer i;
312       TMults(TMults.Lower())=3;
313       for (i=TMults.Lower()+1; i<=TMults.Upper()-1; i++) {
314         TMults(i) = 2;
315       }
316       TMults(TMults.Upper())=3;
317     }
318   }
319 }
320 //=======================================================================
321 //function : GetTolerance
322 //purpose  : Determiner la Tolerance 3d permetant de respecter la Tolerance
323 //           de continuite G1.
324 //=======================================================================
325
326 Standard_Real GeomFill::GetTolerance(const Convert_ParameterisationType TConv,
327                                      const Standard_Real AngleMin, 
328                                      const Standard_Real Radius, 
329                                      const Standard_Real AngularTol, 
330                                      const Standard_Real SpatialTol)
331
332   gp_Ax2 popAx2(gp_Pnt(0, 0, 0), gp_Dir(0,0,1));
333   gp_Circ C (popAx2, Radius);
334   Handle(Geom_Circle) popCircle = new Geom_Circle(C);
335   Handle(Geom_TrimmedCurve) Sect = 
336     new Geom_TrimmedCurve(popCircle ,
337                           0.,Max(AngleMin, 0.02) );
338   // 0.02 est proche d'1 degree, en desous on ne se preocupe pas de la tngence
339   // afin d'eviter des tolerances d'approximation tendant vers 0 !
340   Handle(Geom_BSplineCurve) CtoBspl = 
341         GeomConvert::CurveToBSplineCurve(Sect, TConv);
342   Standard_Real Dist;
343   Dist = CtoBspl->Pole(1).Distance(CtoBspl->Pole(2)) + SpatialTol;
344   return Dist*AngularTol/2;
345 }
346
347 //===========================================================================
348 //function : GetCircle
349 //purpose  : Calculs les poles et poids d'un cercle definie par ses extremites
350 //           et son rayon.
351 //   On evite (si possible) de passer par les convertions pour
352 //  1) Des problemes de performances.
353 //  2) Assurer la coherance entre cette methode est celle qui donne la derive
354 //============================================================================
355 void GeomFill::GetCircle( const Convert_ParameterisationType  TConv,
356                           const gp_Vec& ns1, // Normal rentrente au premier point
357                           const gp_Vec& ns2, // Normal rentrente au second point
358                           const gp_Vec& nplan, // Normal au plan
359                           const gp_Pnt& pts1, 
360                           const gp_Pnt& pts2,
361                           const Standard_Real Rayon, // Rayon (doit etre positif)
362                           const gp_Pnt& Center, 
363                           TColgp_Array1OfPnt& Poles, 
364                           TColStd_Array1OfReal& Weights)
365 {
366   // La classe de convertion
367
368   Standard_Integer i, jj;
369   Standard_Real Cosa,Sina,Angle,Alpha,Cosas2,lambda;
370   gp_Vec temp, np2;
371   Standard_Integer low = Poles.Lower();
372   Standard_Integer upp = Poles.Upper();
373
374   Cosa = ns1.Dot(ns2);
375   Sina = nplan.Dot(ns1.Crossed(ns2));
376
377   if (Cosa<-1.) {Cosa=-1; Sina = 0;}
378   if (Cosa>1.) {Cosa=1; Sina = 0;}
379   Angle = ACos(Cosa);
380   // Recadrage sur ]-pi/2, 3pi/2]
381   if (Sina <0.) {
382     if (Cosa > 0.) Angle = -Angle;
383     else           Angle =  2.*PI - Angle;
384   }
385
386   switch (TConv) {
387   case Convert_QuasiAngular:
388     {
389       GeomFill_QuasiAngularConvertor QConvertor;
390       QConvertor.Init();
391       QConvertor.Section(pts1, Center, nplan, Angle, Poles, Weights);
392       break;
393     }
394   case Convert_Polynomial:
395     { 
396       GeomFill_PolynomialConvertor PConvertor;
397       PConvertor.Init();
398       PConvertor.Section(pts1, Center, nplan, Angle, Poles);
399       Weights.Init(1);
400       break;
401     }
402   default:
403     {
404       // Cas Rational, on utilise une expression directe beaucoup plus
405       // performente que GeomConvert
406       Standard_Integer NbSpan=(Poles.Length()-1)/2;
407
408       Poles(low) = pts1;
409       Poles(upp) = pts2;
410       Weights(low)    = 1;
411       Weights(upp)    = 1; 
412  
413       np2 = nplan.Crossed(ns1);
414
415       Alpha = Angle/((Standard_Real)(NbSpan));
416       Cosas2 = Cos(Alpha/2);
417   
418       for (i=1, jj=low+2; i<= NbSpan-1; i++, jj+=2) {
419         lambda = ((Standard_Real)(i))*Alpha;
420         Cosa = Cos(lambda);
421         Sina = Sin(lambda);
422         temp.SetLinearForm(Cosa-1, ns1, Sina, np2);
423         Poles(jj).SetXYZ(pts1.XYZ() + Rayon*temp.XYZ());
424         Weights(jj)    = 1; 
425       }
426   
427       lambda = 1./(2.*Cosas2*Cosas2);
428       for (i=1, jj=low+1; i<=NbSpan; i++, jj+=2) {
429         temp.SetXYZ(Poles(jj-1).XYZ() + Poles(jj+1).XYZ()
430                     -2.*Center.XYZ());
431         Poles(jj).SetXYZ(Center.XYZ() + lambda*temp.XYZ());
432         Weights(jj)    = Cosas2;
433       }
434     }
435   }
436 }
437
438 Standard_Boolean GeomFill::GetCircle(const Convert_ParameterisationType  TConv, 
439                                       const gp_Vec& ns1,   const gp_Vec& ns2, 
440                                       const gp_Vec& dn1w,  const gp_Vec& dn2w, 
441                                       const gp_Vec& nplan, const gp_Vec& dnplan, 
442                                       const gp_Pnt& pts1,  const gp_Pnt& pts2, 
443                                       const gp_Vec& tang1, const gp_Vec& tang2, 
444                                       const Standard_Real Rayon, 
445                                       const Standard_Real DRayon, 
446                                       const gp_Pnt& Center, 
447                                       const gp_Vec& DCenter, 
448                                       TColgp_Array1OfPnt& Poles, 
449                                       TColgp_Array1OfVec& DPoles,
450                                       TColStd_Array1OfReal& Weights, 
451                                       TColStd_Array1OfReal& DWeights)
452 {
453   Standard_Real Cosa,Sina,Cosas2,Sinas2,Angle,DAngle,Alpha,lambda,Dlambda;
454   gp_Vec temp, np2, dnp2;
455   Standard_Integer i, jj;
456   Standard_Integer NbSpan=(Poles.Length()-1)/2;
457   Standard_Integer low = Poles.Lower();
458   Standard_Integer upp = Poles.Upper();
459
460   Cosa = ns1.Dot(ns2);
461   Sina = nplan.Dot(ns1.Crossed(ns2));
462
463   if (Cosa<-1.){Cosa=-1; Sina = 0;}
464   if (Cosa>1.) {Cosa=1; Sina = 0;}
465   Angle = ACos(Cosa);
466   // Recadrage sur ]-pi/2, 3pi/2]
467   if (Sina <0.) {
468     if (Cosa > 0.) Angle = -Angle;
469     else           Angle =  2.*PI - Angle;
470   }
471
472   if (Abs(Sina)>Abs(Cosa)) {
473     DAngle = -(dn1w.Dot(ns2) + ns1.Dot(dn2w))/Sina;
474   }
475   else{
476     DAngle = (dnplan.Dot(ns1.Crossed(ns2)) + nplan.Dot(dn1w.Crossed(ns2) 
477                                            + ns1.Crossed(dn2w)))/Cosa;
478   }
479
480 // Aux Extremites.  
481   Poles(low) = pts1;
482   Poles(upp) = pts2;
483   Weights(low)  = 1;
484   Weights(upp) = 1;
485
486   DPoles(low) = tang1;
487   DPoles(upp) = tang2;
488   DWeights(low) = 0;
489   DWeights(upp) = 0;
490
491
492   switch (TConv) {
493   case Convert_QuasiAngular:
494     {
495       GeomFill_QuasiAngularConvertor QConvertor;
496       QConvertor.Init();
497       QConvertor.Section(pts1, tang1, 
498                         Center, DCenter, 
499                         nplan, dnplan,
500                         Angle, DAngle,
501                         Poles, DPoles,
502                         Weights, DWeights);
503       return Standard_True;
504   }
505   case Convert_Polynomial:
506     {
507       GeomFill_PolynomialConvertor PConvertor;
508       PConvertor.Init();
509       PConvertor.Section(pts1, tang1, 
510                         Center, DCenter, 
511                         nplan, dnplan,
512                         Angle, DAngle,
513                         Poles, DPoles);
514       Weights.Init(1);
515       DWeights.Init(0);
516       return Standard_True;
517   }
518
519   default:
520      // Cas rationel classique
521     { 
522       np2 = nplan.Crossed(ns1);
523       dnp2 = dnplan.Crossed(ns1).Added(nplan.Crossed(dn1w));
524  
525       Alpha = Angle/((Standard_Real)(NbSpan));
526       Cosas2 = Cos(Alpha/2);
527       Sinas2 = Sin(Alpha/2);
528
529       for (i=1, jj=low+2; i<= NbSpan-1; i++, jj+=2) {
530         lambda = ((Standard_Real)(i))*Alpha;
531         Cosa = Cos(lambda);
532         Sina = Sin(lambda);
533         temp.SetLinearForm(Cosa-1,ns1,Sina,np2);
534         Poles(jj).SetXYZ(pts1.XYZ() + Rayon*temp.XYZ());
535     
536         DPoles(jj).SetLinearForm(DRayon, temp, tang1);
537         temp.SetLinearForm(-Sina,ns1,Cosa,np2);
538         temp.Multiply(((Standard_Real)(i))/((Standard_Real)(NbSpan))*DAngle);
539         temp.Add(((Cosa-1)*dn1w).Added(Sina*dnp2));
540         DPoles(jj)+= Rayon*temp;
541       }
542   
543       lambda = 1./(2.*Cosas2*Cosas2);
544       Dlambda = (lambda*Sinas2*DAngle)/(Cosas2*NbSpan);
545
546       for (i=1, jj=low; i<=NbSpan; i++, jj+=2) {
547         temp.SetXYZ(Poles(jj).XYZ() + Poles(jj+2).XYZ()
548                     -2.*Center.XYZ());
549         Poles(jj+1).SetXYZ(Center.XYZ()+lambda*temp.XYZ());
550         DPoles(jj+1).SetLinearForm(Dlambda, temp,
551                                    1.-2*lambda, DCenter,
552                                    lambda, (DPoles(jj)+ DPoles(jj+2)));       
553       }
554   
555       // Les poids
556       Dlambda = -Sinas2*DAngle/(2*NbSpan);
557       for (i=low; i<upp; i+=2) {
558         Weights(i)    = 1.;
559         Weights(i+1)  = Cosas2;
560         DWeights(i)   = 0.;
561         DWeights(i+1) = Dlambda;
562       }
563     }
564     return Standard_True;
565   }
566   return Standard_False;
567 }
568
569 Standard_Boolean GeomFill::GetCircle(const Convert_ParameterisationType  TConv, 
570                                      const gp_Vec& ns1,   const gp_Vec& ns2, 
571                                      const gp_Vec& dn1w,  const gp_Vec& dn2w,
572                                      const gp_Vec& d2n1w, const gp_Vec& d2n2w,
573                                      const gp_Vec& nplan, const gp_Vec& dnplan, 
574                                      const gp_Vec& d2nplan, 
575                                      const gp_Pnt& pts1,  const gp_Pnt& pts2, 
576                                      const gp_Vec& tang1, const gp_Vec& tang2,
577                                      const gp_Vec& Dtang1, const gp_Vec& Dtang2, 
578                                      const Standard_Real Rayon, 
579                                      const Standard_Real DRayon,
580                                      const Standard_Real D2Rayon, 
581                                      const gp_Pnt& Center, 
582                                      const gp_Vec& DCenter,
583                                      const gp_Vec& D2Center,
584                                      TColgp_Array1OfPnt& Poles, 
585                                      TColgp_Array1OfVec& DPoles,
586                                      TColgp_Array1OfVec& D2Poles,
587                                      TColStd_Array1OfReal& Weights, 
588                                      TColStd_Array1OfReal& DWeights,
589                                      TColStd_Array1OfReal& D2Weights)
590 {
591   Standard_Real Cosa,Sina,Cosas2,Sinas2;
592   Standard_Real Angle, DAngle, D2Angle, Alpha;
593   Standard_Real lambda, Dlambda, D2lambda, aux;
594   gp_Vec temp, dtemp, np2, dnp2, d2np2;
595   Standard_Integer i, jj;
596   Standard_Integer NbSpan=(Poles.Length()-1)/2;
597   Standard_Integer low = Poles.Lower();
598   Standard_Integer upp = Poles.Upper();
599
600   Cosa = ns1.Dot(ns2);
601   Sina = nplan.Dot(ns1.Crossed(ns2));
602
603   if (Cosa<-1.){Cosa=-1; Sina = 0;}
604   if (Cosa>1.) {Cosa=1; Sina = 0;}
605   Angle = ACos(Cosa);
606   // Recadrage sur ]-pi/2, 3pi/2]
607   if (Sina <0.) {
608     if (Cosa > 0.) Angle = -Angle;
609     else           Angle =  2.*PI - Angle;
610   }
611
612   if (Abs(Sina)>Abs(Cosa)) {
613     aux = dn1w.Dot(ns2) + ns1.Dot(dn2w);
614     DAngle  = -aux/Sina;
615     D2Angle = -(d2n1w.Dot(ns2) + 2*dn1w.Dot(dn2w) + ns1.Dot(d2n2w))/Sina
616               + aux*(dnplan.Dot(ns1.Crossed(ns2)) + nplan.Dot(dn1w.Crossed(ns2) 
617               + ns1.Crossed(dn2w)))/(Sina*Sina); 
618   }
619   else{
620     temp = dn1w.Crossed(ns2) + ns1.Crossed(dn2w);
621     DAngle = (dnplan.Dot(ns1.Crossed(ns2)) + nplan.Dot(temp))/Cosa;
622     D2Angle = ( d2nplan.Dot(ns1.Crossed(ns2)) +2*dnplan.Dot(temp)
623               + nplan.Dot(d2n1w.Crossed(ns2) + 2*dn1w.Crossed(dn2w) 
624               + ns1.Crossed(d2n2w)) )/Cosa
625               - ( dn1w.Dot(ns2) + ns1.Dot(dn2w))
626                 * (dnplan.Dot(ns1.Crossed(ns2)) + nplan.Dot(temp)) /(Cosa*Cosa);
627   }
628
629 // Aux Extremites.  
630   Poles(low) = pts1;
631   Poles(upp) = pts2;
632   Weights(low)  = 1;
633   Weights(upp) = 1;
634
635   DPoles(low) = tang1;
636   DPoles(upp) = tang2;
637   DWeights(low) = 0;
638   DWeights(upp) = 0;
639
640   D2Poles(low) = Dtang1;
641   D2Poles(upp) = Dtang2;
642   D2Weights(low) = 0;
643   D2Weights(upp) = 0;  
644
645
646   switch (TConv) {
647   case Convert_QuasiAngular:
648     {
649       GeomFill_QuasiAngularConvertor QConvertor;
650       QConvertor.Init();
651       QConvertor.Section(pts1, tang1, Dtang1,
652                            Center, DCenter, D2Center,
653                            nplan, dnplan, d2nplan,
654                            Angle, DAngle, D2Angle,
655                            Poles, DPoles, D2Poles,
656                            Weights, DWeights, D2Weights);
657       return Standard_True;  
658   }
659   case Convert_Polynomial:
660     {
661       GeomFill_PolynomialConvertor PConvertor;
662       PConvertor.Init();
663       PConvertor.Section(pts1, tang1, Dtang1,
664                            Center, DCenter, D2Center,
665                            nplan, dnplan,  d2nplan,
666                            Angle, DAngle, D2Angle,
667                            Poles, DPoles, D2Poles);
668       Weights.Init(1);
669       DWeights.Init(0);
670       D2Weights.Init(0);
671       return Standard_True;
672   }
673
674   default:
675     { 
676       np2 = nplan.Crossed(ns1);
677       dnp2 = dnplan.Crossed(ns1).Added(nplan.Crossed(dn1w));
678       d2np2 = d2nplan.Crossed(ns1).Added(nplan.Crossed(dn2w));
679       d2np2 += 2*dnplan.Crossed(dn1w);
680  
681       Alpha = Angle/((Standard_Real)(NbSpan));
682       Cosas2 = Cos(Alpha/2);
683       Sinas2 = Sin(Alpha/2);
684
685       for (i=1, jj=low+2; i<= NbSpan-1; i++, jj+=2) {
686         lambda = ((Standard_Real)(i))*Alpha;
687         Cosa = Cos(lambda);
688         Sina = Sin(lambda);
689         temp.SetLinearForm(Cosa-1,ns1,Sina,np2);
690         Poles(jj).SetXYZ(pts1.XYZ() + Rayon*temp.XYZ());
691     
692         DPoles(jj).SetLinearForm(DRayon, temp, tang1);
693         dtemp.SetLinearForm(-Sina,ns1,Cosa,np2);
694         aux = ((Standard_Real)(i))/((Standard_Real)(NbSpan));
695         dtemp.Multiply(aux*DAngle);
696         dtemp.Add(((Cosa-1)*dn1w).Added(Sina*dnp2));
697         DPoles(jj)+= Rayon*dtemp;
698
699         D2Poles(jj).SetLinearForm(D2Rayon, temp, 
700                                   2*DRayon, dtemp, Dtang1);
701         temp.SetLinearForm(Cosa-1, dn2w, Sina, d2np2);
702         dtemp.SetLinearForm(-Sina,ns1,Cosa,np2);
703         temp+= (aux*aux*D2Angle)*dtemp;
704         dtemp.SetLinearForm(-Sina, dn1w+np2, Cosa, dnp2,
705                             -Cosa, ns1);
706         temp+=(aux*DAngle)*dtemp;
707         D2Poles(jj)+= Rayon*temp;       
708       }
709   
710       lambda = 1./(2.*Cosas2*Cosas2);
711       Dlambda = (lambda*Sinas2*DAngle)/(Cosas2*NbSpan);
712       aux = Sinas2/Cosas2;
713       D2lambda = (  Dlambda * aux*DAngle 
714                   + D2Angle * aux*lambda
715                   + (1+aux*aux)*(DAngle/(2*NbSpan)) * DAngle*lambda ) 
716                / NbSpan;
717       for (i=1, jj=low; i<=NbSpan; i++, jj+=2) {
718         temp.SetXYZ(Poles(jj).XYZ() + Poles(jj+2).XYZ()
719                     -2.*Center.XYZ());
720         Poles(jj+1).SetXYZ(Center.XYZ()+lambda*temp.XYZ());
721
722
723         dtemp.SetXYZ(DPoles(jj).XYZ() + DPoles(jj+2).XYZ()
724                     -2.*DCenter.XYZ());
725         DPoles(jj+1).SetLinearForm(Dlambda, temp,
726                                    lambda, dtemp,
727                                    DCenter);
728         D2Poles(jj+1).SetLinearForm(D2lambda, temp,
729                                     2*Dlambda, dtemp,
730                                     lambda, (D2Poles(jj)+ D2Poles(jj+2)));
731         D2Poles(jj+1)+= (1-2*lambda)*D2Center;
732       }
733   
734       // Les poids
735       Dlambda = -Sinas2*DAngle/(2*NbSpan);
736       D2lambda = -Sinas2*D2Angle/(2*NbSpan)
737                  -Cosas2*Pow(DAngle/(2*NbSpan),2);
738
739       for (i=low; i<upp; i+=2) {
740         Weights(i)    = 1.;
741         Weights(i+1)  = Cosas2;
742         DWeights(i)   = 0.;
743         DWeights(i+1) = Dlambda;
744         D2Weights(i)  = 0.;
745         D2Weights(i+1)  = D2lambda;
746       }
747     }
748     return Standard_True;
749   }
750   return Standard_False;
751 }