1 // File:        GeomFill_PolynomialConvertor.cxx
2 // Created:     Fri Jul 18 17:43:05 1997
3 // Author:      Philippe MANGIN
4 //              <pmn@sgi29>
7 #include <GeomFill_PolynomialConvertor.ixx>
10 #include <PLib.hxx>
11 #include <gp_Mat.hxx>
13 #include <Convert_CompPolynomialToPoles.hxx>
15 #include <TColStd_Array1OfReal.hxx>
16 #include <TColStd_HArray2OfReal.hxx>
17 #include <TColStd_HArray1OfInteger.hxx>
18 #include <TColStd_HArray1OfReal.hxx>
21 GeomFill_PolynomialConvertor::GeomFill_PolynomialConvertor():
22                     Ordre(8),
23                     myinit(Standard_False),
24                     BH(1, Ordre, 1, Ordre)
25 {
26 }
28 Standard_Boolean GeomFill_PolynomialConvertor::Initialized() const
29 {
30  return myinit;
31 }
33 void GeomFill_PolynomialConvertor::Init()
34 {
35   if (myinit) return; //On n'initialise qu'une fois
36   Standard_Integer ii, jj;
37   Standard_Real terme;
38   math_Matrix H(1,Ordre, 1,Ordre), B(1,Ordre, 1,Ordre);
39   Handle(TColStd_HArray1OfReal)
40     Coeffs = new (TColStd_HArray1OfReal) (1, Ordre*Ordre),
41     TrueInter  = new  (TColStd_HArray1OfReal) (1,2);
43   Handle(TColStd_HArray2OfReal)
44     Poles1d = new (TColStd_HArray2OfReal) (1, Ordre, 1, Ordre),
45     Inter   = new (TColStd_HArray2OfReal) (1,1,1,2);
47   //Calcul de B
48   Inter->SetValue(1, 1, -1);
49   Inter->SetValue(1, 2,  1);
50   TrueInter->SetValue(1, -1);
51   TrueInter->SetValue(2,  1);
53   Coeffs->Init(0);
54   for (ii=1; ii<=Ordre; ii++) {  Coeffs->SetValue(ii+(ii-1)*Ordre, 1); }
56    //Convertion ancienne formules
57      Handle(TColStd_HArray1OfInteger) Ncf
58        = new (TColStd_HArray1OfInteger)(1,1);
59      Ncf->Init(Ordre);
61      Convert_CompPolynomialToPoles
62           AConverter(1, 1, 8, 8,
63                      Ncf,
64                      Coeffs,
65                      Inter,
66                      TrueInter);
67 /*  Convert_CompPolynomialToPoles
68          AConverter(8, Ordre-1,  Ordre-1,
69                     Coeffs,
70                     Inter,
71                     TrueInter); En attente du bon Geomlite*/
72   AConverter.Poles(Poles1d);
74   for (jj=1; jj<=Ordre; jj++) {
75       for (ii=1; ii<=Ordre; ii++) {
76         terme = Poles1d->Value(ii,jj);
77         if (Abs(terme-1) < 1.e-9) terme = 1 ; //petite retouche
78         if (Abs(terme+1) < 1.e-9) terme = -1;
79         B(ii, jj) =  terme;
80       }
81   }
82   //Calcul de H
83   myinit = PLib::HermiteCoefficients(-1, 1, Ordre/2-1, Ordre/2-1, H);
84   H.Transpose();
86   if (!myinit) return;
88   // reste l'essentiel
89   BH = B * H;
90 }
92 void GeomFill_PolynomialConvertor::Section(const gp_Pnt& FirstPnt,
93                                        const gp_Pnt& Center,
94                                        const gp_Vec& Dir,
95                                        const Standard_Real Angle,
96                                        TColgp_Array1OfPnt& Poles) const
97 {
98   math_Vector Vx(1, Ordre), Vy(1, Ordre);
99   math_Vector Px(1, Ordre), Py(1, Ordre);
100   Standard_Integer ii;
101   Standard_Real Cos_b = Cos(Angle), Sin_b = Sin(Angle);
102   Standard_Real beta, beta2, beta3;
103   gp_Vec V1(Center, FirstPnt), V2;
104   V2 =  Dir^V1;
105   beta = Angle/2;
106   beta2 = beta * beta;
107   beta3 = beta * beta2;
109   // Calcul de la transformation
110   gp_Mat M(V1.X(), V2.X(), 0,
111            V1.Y(), V2.Y(), 0,
112            V1.Z(), V2.Z(), 0);
114   // Calcul des contraintes  -----------
115   Vx(1) = 1;                   Vy(1) = 0;
116   Vx(2) = 0;                   Vy(2) = beta;
117   Vx(3) = -beta2;              Vy(3) = 0;
118   Vx(4) = 0;                   Vy(4) = -beta3;
119   Vx(5) = Cos_b;               Vy(5) = Sin_b;
120   Vx(6) = -beta*Sin_b;         Vy(6) = beta*Cos_b;
121   Vx(7) = -beta2*Cos_b;        Vy(7) = -beta2*Sin_b;
122   Vx(8) = beta3*Sin_b;         Vy(8) = -beta3*Cos_b;
124   // Calcul des poles
125   Px = BH * Vx;
126   Py = BH * Vy;
127   gp_XYZ pnt;
128   for (ii=1; ii<=Ordre; ii++) {
129      pnt.SetCoord(Px(ii), Py(ii), 0);
130      pnt *= M;
131      pnt += Center.XYZ();
132      Poles(ii).ChangeCoord() = pnt;
133   }
134 }
136 void GeomFill_PolynomialConvertor::Section(const gp_Pnt& FirstPnt,
137                                        const gp_Vec& DFirstPnt,
138                                        const gp_Pnt& Center,
139                                        const gp_Vec& DCenter,
140                                        const gp_Vec& Dir,
141                                        const gp_Vec& DDir,
142                                        const Standard_Real Angle,
143                                        const Standard_Real DAngle,
144                                        TColgp_Array1OfPnt& Poles,
145                                        TColgp_Array1OfVec& DPoles) const
146 {
147   math_Vector Vx(1, Ordre), Vy(1, Ordre),
148               DVx(1, Ordre), DVy(1, Ordre);
149   math_Vector Px(1, Ordre), Py(1, Ordre),
150               DPx(1, Ordre), DPy(1, Ordre);
151   Standard_Integer ii;
152   Standard_Real Cos_b = Cos(Angle), Sin_b = Sin(Angle);
153   Standard_Real beta, beta2, beta3, bprim;
154   gp_Vec V1(Center, FirstPnt), V1Prim, V2;
155   V2 =  Dir^V1;
156   beta = Angle/2;
157   bprim = DAngle/2;
158   beta2 = beta * beta;
159   beta3 = beta * beta2;
161   // Calcul des  transformations
162   gp_Mat M (V1.X(), V2.X(), 0,
163             V1.Y(), V2.Y(), 0,
164             V1.Z(), V2.Z(), 0);
165   V1Prim = DFirstPnt - DCenter;
166   V2 = (DDir^V1) + (Dir^V1Prim);
167   gp_Mat MPrim (V1Prim.X(), V2.X(), 0,
168                 V1Prim.Y(), V2.Y(), 0,
169                 V1Prim.Z(), V2.Z(), 0);
173   // Calcul des contraintes  -----------
174   Vx(1) = 1;                   Vy(1) = 0;
175   Vx(2) = 0;                   Vy(2) = beta;
176   Vx(3) = -beta2;              Vy(3) = 0;
177   Vx(4) = 0;                   Vy(4) = -beta3;
178   Vx(5) = Cos_b;               Vy(5) = Sin_b;
179   Vx(6) = -beta*Sin_b;         Vy(6) = beta*Cos_b;
180   Vx(7) = -beta2*Cos_b;        Vy(7) = -beta2*Sin_b;
181   Vx(8) = beta3*Sin_b;         Vy(8) = -beta3*Cos_b;
183   Standard_Real b_bprim = bprim*beta,
184                 b2_bprim = bprim*beta2;
185   DVx(1) = 0;                  DVy(1) = 0;
186   DVx(2) = 0;                  DVy(2) = bprim;
187   DVx(3) = -2*b_bprim;         DVy(3) = 0;
188   DVx(4) = 0;                  DVy(4) = -3*b2_bprim;
189   DVx(5) = -2*bprim*Sin_b;     DVy(5) = 2*bprim*Cos_b;
190   DVx(6) = -bprim*Sin_b - 2*b_bprim*Cos_b;
191                   DVy(6) = bprim*Cos_b - 2*b_bprim*Sin_b;
192   DVx(7) = 2*b_bprim*(-Cos_b + beta*Sin_b);
193                   DVy(7) = -2*b_bprim*(Sin_b+beta*Cos_b);
194   DVx(8) = b2_bprim*(3*Sin_b + 2*beta*Cos_b);
195                   DVy(8) = b2_bprim*(2*beta*Sin_b - 3*Cos_b);
197   // Calcul des poles
198   Px = BH * Vx;
199   Py = BH * Vy;
200   DPx = BH * DVx;
201   DPy = BH * DVy;
202   gp_XYZ P, DP, aux;
204   for (ii=1; ii<=Ordre; ii++) {
205      P.SetCoord(Px(ii), Py(ii), 0);
206      Poles(ii).ChangeCoord() = M*P + Center.XYZ();
207      P *= MPrim;
208      DP.SetCoord(DPx(ii), DPy(ii), 0);
209      DP *= M;
210      aux.SetLinearForm(1, P, 1, DP, DCenter.XYZ());
211      DPoles(ii).SetXYZ(aux);
212   }
213 }
215 void GeomFill_PolynomialConvertor::Section(const gp_Pnt& FirstPnt,
216                                        const gp_Vec& DFirstPnt,
217                                        const gp_Vec& D2FirstPnt,
218                                        const gp_Pnt& Center,
219                                        const gp_Vec& DCenter,
220                                        const gp_Vec& D2Center,
221                                        const gp_Vec& Dir,
222                                        const gp_Vec& DDir,
223                                        const gp_Vec& D2Dir,
224                                        const Standard_Real Angle,
225                                        const Standard_Real DAngle,
226                                        const Standard_Real D2Angle,
227                                        TColgp_Array1OfPnt& Poles,
228                                        TColgp_Array1OfVec& DPoles,
229                                        TColgp_Array1OfVec& D2Poles) const
230 {
231   math_Vector Vx(1, Ordre), Vy(1, Ordre),
232               DVx(1, Ordre), DVy(1, Ordre),
233               D2Vx(1, Ordre), D2Vy(1, Ordre);
234   math_Vector Px(1, Ordre), Py(1, Ordre),
235               DPx(1, Ordre), DPy(1, Ordre),
236               D2Px(1, Ordre), D2Py(1, Ordre);
238   Standard_Integer ii;
239   Standard_Real aux, Cos_b = Cos(Angle), Sin_b = Sin(Angle);
240   Standard_Real beta, beta2, beta3, bprim, bprim2, bsecn;
241   gp_Vec V1(Center, FirstPnt), V1Prim, V1Secn, V2;
242   V2 =  Dir^V1;
243   beta = Angle/2;
244   bprim = DAngle/2;
245   bsecn = D2Angle/2;
246   bsecn = D2Angle/2;
247   beta2 = beta * beta;
248   beta3 = beta * beta2;
249   bprim2 =  bprim*bprim;
251   // Calcul des  transformations
252   gp_Mat M (V1.X(), V2.X(), 0,
253             V1.Y(), V2.Y(), 0,
254             V1.Z(), V2.Z(), 0);
255   V1Prim = DFirstPnt - DCenter;
256   V2 = (DDir^V1) + (Dir^V1Prim);
258   gp_Mat MPrim (V1Prim.X(), V2.X(), 0,
259                 V1Prim.Y(), V2.Y(), 0,
260                 V1Prim.Z(), V2.Z(), 0);
261   V1Secn = D2FirstPnt - D2Center;
262   V2 = DDir^V1Prim;
263   V2 *= 2;
264   V2 +=  (D2Dir^V1) + (Dir^V1Secn);
267   gp_Mat MSecn (V1Secn.X(), V2.X(), 0,
268                 V1Secn.Y(), V2.Y(), 0,
269                 V1Secn.Z(), V2.Z(), 0);
272   // Calcul des contraintes  -----------
273   Vx(1) = 1;                   Vy(1) = 0;
274   Vx(2) = 0;                   Vy(2) = beta;
275   Vx(3) = -beta2;              Vy(3) = 0;
276   Vx(4) = 0;                   Vy(4) = -beta3;
277   Vx(5) = Cos_b;               Vy(5) = Sin_b;
278   Vx(6) = -beta*Sin_b;         Vy(6) = beta*Cos_b;
279   Vx(7) = -beta2*Cos_b;        Vy(7) = -beta2*Sin_b;
280   Vx(8) = beta3*Sin_b;         Vy(8) = -beta3*Cos_b;
282   Standard_Real b_bprim = bprim*beta,
283                 b2_bprim = bprim*beta2,
284                 b_bsecn  = bsecn*beta;
285   DVx(1) = 0;                  DVy(1) = 0;
286   DVx(2) = 0;                  DVy(2) = bprim;
287   DVx(3) = -2*b_bprim;         DVy(3) = 0;
288   DVx(4) = 0;                  DVy(4) = -3*b2_bprim;
289   DVx(5) = -2*bprim*Sin_b;     DVy(5) = 2*bprim*Cos_b;
290   DVx(6) = -bprim*Sin_b - 2*b_bprim*Cos_b;
291                   DVy(6) = bprim*Cos_b - 2*b_bprim*Sin_b;
292   DVx(7) = 2*b_bprim*(-Cos_b + beta*Sin_b);
293                   DVy(7) = -2*b_bprim*(Sin_b + beta*Cos_b);
294   DVx(8) = b2_bprim*(3*Sin_b + 2*beta*Cos_b);
295                   DVy(8) = b2_bprim*(2*beta*Sin_b - 3*Cos_b);
298   D2Vx(1) = 0;                  D2Vy(1) = 0;
299   D2Vx(2) = 0;                  D2Vy(2) = bsecn;
300   D2Vx(3) = -2*(bprim2+b_bsecn); D2Vy(3) = 0;
301   D2Vx(4) = 0;    D2Vy(4) = -3*beta*(2*bprim2+b_bsecn);
302   D2Vx(5) = -2*(bsecn*Sin_b + 2*bprim2*Cos_b);
303                   D2Vy(5) = 2*(bsecn*Cos_b - 2*bprim2*Sin_b);
304   D2Vx(6) = (4*beta*bprim2-bsecn)*Sin_b - 2*(2*bprim2+b_bsecn)*Cos_b;
306                   D2Vy(6) = (bsecn - 4*beta*bprim2)*Cos_b
307                           - 2*(b_bsecn + 2*bprim2)*Sin_b;
309   aux =  2*(bprim2+b_bsecn);
310   D2Vx(7) = aux*(-Cos_b + beta*Sin_b)
311           + 2*beta*bprim2*(2*beta*Cos_b + 3*Sin_b);
313                   D2Vy(7) = -aux*(Sin_b + beta*Cos_b)
314                           - 2*beta*bprim2*(3*Cos_b - 2*beta*Sin_b);
316   aux = beta*(2*bprim2+b_bsecn);
317   D2Vx(8) = aux * (3*Sin_b + 2*beta*Cos_b)
318           + 4*beta2*bprim2 * (2*Cos_b - beta*Sin_b);
320                   D2Vy(8)= aux * (2*beta*Sin_b - 3*Cos_b)
321                          + 4*beta2*bprim2 * (2*Sin_b + beta*Cos_b);
324   // Calcul des poles
325   Px = BH * Vx;
326   Py = BH * Vy;
327   DPx = BH * DVx;
328   DPy = BH * DVy;
329   D2Px = BH * D2Vx;
330   D2Py = BH * D2Vy;
332   gp_XYZ P, DP, D2P, auxyz;
333   for (ii=1; ii<=Ordre; ii++) {
334      P.SetCoord(Px(ii), Py(ii), 0);
335      DP.SetCoord(DPx(ii), DPy(ii), 0);
336      D2P.SetCoord(D2Px(ii), D2Py(ii), 0);
338      Poles(ii).ChangeCoord() = M*P + Center.XYZ();
339      auxyz.SetLinearForm(1, MPrim*P,
340                          1, M*DP,
341                          DCenter.XYZ());
342      DPoles(ii).SetXYZ(auxyz);
343      P  *= MSecn;
344      DP *= MPrim;
345      D2P*= M;
346      auxyz.SetLinearForm(1, P,
347                          2, DP,
348                          1, D2P,
349                          D2Center.XYZ());
350      D2Poles(ii).SetXYZ(auxyz);
351   }
352 }