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