projects / occt.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
0022792: Globally defined symbol PI conflicts with VTK definition (Intel compiler)
[occt.git] / src / GeomFill / GeomFill.cxx
CommitLineData
7fd59977 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>
7fd59977 29
7fd59977 30
31//=======================================================================
32//function : Surface
33//purpose :
34//=======================================================================
35Handle(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
182void 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 =
c6541a0c 206 (Standard_Integer)(Ceiling(3.*Abs(MaxAng)/2./M_PI));
7fd59977 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
229void 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
267void 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
291void 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
326Standard_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//============================================================================
355void 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;
c6541a0c 383 else Angle = 2.*M_PI - Angle;
7fd59977 384 }
385
386 switch (TConv) {
387 case Convert_QuasiAngular:
f69df442
RK		 388 {
389 GeomFill_QuasiAngularConvertor QConvertor;
390 QConvertor.Init();
391 QConvertor.Section(pts1, Center, nplan, Angle, Poles, Weights);
7fd59977 392 break;
393 }
394 case Convert_Polynomial:
f69df442
RK		 395 {
396 GeomFill_PolynomialConvertor PConvertor;
397 PConvertor.Init();
398 PConvertor.Section(pts1, Center, nplan, Angle, Poles);
7fd59977 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
438Standard_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;
c6541a0c 469 else Angle = 2.*M_PI - Angle;
7fd59977 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 {
f69df442
RK		 495 GeomFill_QuasiAngularConvertor QConvertor;
496 QConvertor.Init();
497 QConvertor.Section(pts1, tang1,
7fd59977 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 {
f69df442
RK		 507 GeomFill_PolynomialConvertor PConvertor;
508 PConvertor.Init();
509 PConvertor.Section(pts1, tang1,
7fd59977 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
569Standard_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;
c6541a0c 609 else Angle = 2.*M_PI - Angle;
7fd59977 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 {
f69df442
RK		 649 GeomFill_QuasiAngularConvertor QConvertor;
650 QConvertor.Init();
651 QConvertor.Section(pts1, tang1, Dtang1,
7fd59977 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 {
f69df442
RK		 661 GeomFill_PolynomialConvertor PConvertor;
662 PConvertor.Init();
663 PConvertor.Section(pts1, tang1, Dtang1,
7fd59977 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}
Open CASCADE Technology repository