1 // Created on: 1997-07-18
2 // Created by: Philippe MANGIN
3 // Copyright (c) 1997-1999 Matra Datavision
4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
6 // This file is part of Open CASCADE Technology software library.
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
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.
14 // Alternatively, this file may be used under the terms of Open CASCADE
15 // commercial license or contractual agreement.
18 #include <Convert_CompPolynomialToPoles.hxx>
19 #include <GeomFill_PolynomialConvertor.hxx>
24 #include <StdFail_NotDone.hxx>
25 #include <TColStd_Array1OfReal.hxx>
26 #include <TColStd_HArray1OfInteger.hxx>
27 #include <TColStd_HArray1OfReal.hxx>
28 #include <TColStd_HArray2OfReal.hxx>
30 GeomFill_PolynomialConvertor::GeomFill_PolynomialConvertor():
32 myinit(Standard_False),
33 BH(1, Ordre, 1, Ordre)
37 Standard_Boolean GeomFill_PolynomialConvertor::Initialized() const
42 void GeomFill_PolynomialConvertor::Init()
44 if (myinit) return; //On n'initialise qu'une fois
45 Standard_Integer ii, jj;
47 math_Matrix H(1,Ordre, 1,Ordre), B(1,Ordre, 1,Ordre);
48 Handle(TColStd_HArray1OfReal)
49 Coeffs = new (TColStd_HArray1OfReal) (1, Ordre*Ordre),
50 TrueInter = new (TColStd_HArray1OfReal) (1,2);
52 Handle(TColStd_HArray2OfReal)
53 Poles1d = new (TColStd_HArray2OfReal) (1, Ordre, 1, Ordre),
54 Inter = new (TColStd_HArray2OfReal) (1,1,1,2);
57 Inter->SetValue(1, 1, -1);
58 Inter->SetValue(1, 2, 1);
59 TrueInter->SetValue(1, -1);
60 TrueInter->SetValue(2, 1);
63 for (ii=1; ii<=Ordre; ii++) { Coeffs->SetValue(ii+(ii-1)*Ordre, 1); }
65 //Convertion ancienne formules
66 Handle(TColStd_HArray1OfInteger) Ncf
67 = new (TColStd_HArray1OfInteger)(1,1);
70 Convert_CompPolynomialToPoles
71 AConverter(1, 1, 8, 8,
76 /* Convert_CompPolynomialToPoles
77 AConverter(8, Ordre-1, Ordre-1,
80 TrueInter); En attente du bon Geomlite*/
81 AConverter.Poles(Poles1d);
83 for (jj=1; jj<=Ordre; jj++) {
84 for (ii=1; ii<=Ordre; ii++) {
85 terme = Poles1d->Value(ii,jj);
86 if (Abs(terme-1) < 1.e-9) terme = 1 ; //petite retouche
87 if (Abs(terme+1) < 1.e-9) terme = -1;
92 myinit = PLib::HermiteCoefficients(-1, 1, Ordre/2-1, Ordre/2-1, H);
101 void GeomFill_PolynomialConvertor::Section(const gp_Pnt& FirstPnt,
102 const gp_Pnt& Center,
104 const Standard_Real Angle,
105 TColgp_Array1OfPnt& Poles) const
107 math_Vector Vx(1, Ordre), Vy(1, Ordre);
108 math_Vector Px(1, Ordre), Py(1, Ordre);
110 Standard_Real Cos_b = Cos(Angle), Sin_b = Sin(Angle);
111 Standard_Real beta, beta2, beta3;
112 gp_Vec V1(Center, FirstPnt), V2;
116 beta3 = beta * beta2;
118 // Calcul de la transformation
119 gp_Mat M(V1.X(), V2.X(), 0,
123 // Calcul des contraintes -----------
124 Vx(1) = 1; Vy(1) = 0;
125 Vx(2) = 0; Vy(2) = beta;
126 Vx(3) = -beta2; Vy(3) = 0;
127 Vx(4) = 0; Vy(4) = -beta3;
128 Vx(5) = Cos_b; Vy(5) = Sin_b;
129 Vx(6) = -beta*Sin_b; Vy(6) = beta*Cos_b;
130 Vx(7) = -beta2*Cos_b; Vy(7) = -beta2*Sin_b;
131 Vx(8) = beta3*Sin_b; Vy(8) = -beta3*Cos_b;
137 for (ii=1; ii<=Ordre; ii++) {
138 pnt.SetCoord(Px(ii), Py(ii), 0);
141 Poles(ii).ChangeCoord() = pnt;
145 void GeomFill_PolynomialConvertor::Section(const gp_Pnt& FirstPnt,
146 const gp_Vec& DFirstPnt,
147 const gp_Pnt& Center,
148 const gp_Vec& DCenter,
151 const Standard_Real Angle,
152 const Standard_Real DAngle,
153 TColgp_Array1OfPnt& Poles,
154 TColgp_Array1OfVec& DPoles) const
156 math_Vector Vx(1, Ordre), Vy(1, Ordre),
157 DVx(1, Ordre), DVy(1, Ordre);
158 math_Vector Px(1, Ordre), Py(1, Ordre),
159 DPx(1, Ordre), DPy(1, Ordre);
161 Standard_Real Cos_b = Cos(Angle), Sin_b = Sin(Angle);
162 Standard_Real beta, beta2, beta3, bprim;
163 gp_Vec V1(Center, FirstPnt), V1Prim, V2;
168 beta3 = beta * beta2;
170 // Calcul des transformations
171 gp_Mat M (V1.X(), V2.X(), 0,
174 V1Prim = DFirstPnt - DCenter;
175 V2 = (DDir^V1) + (Dir^V1Prim);
176 gp_Mat MPrim (V1Prim.X(), V2.X(), 0,
177 V1Prim.Y(), V2.Y(), 0,
178 V1Prim.Z(), V2.Z(), 0);
182 // Calcul des contraintes -----------
183 Vx(1) = 1; Vy(1) = 0;
184 Vx(2) = 0; Vy(2) = beta;
185 Vx(3) = -beta2; Vy(3) = 0;
186 Vx(4) = 0; Vy(4) = -beta3;
187 Vx(5) = Cos_b; Vy(5) = Sin_b;
188 Vx(6) = -beta*Sin_b; Vy(6) = beta*Cos_b;
189 Vx(7) = -beta2*Cos_b; Vy(7) = -beta2*Sin_b;
190 Vx(8) = beta3*Sin_b; Vy(8) = -beta3*Cos_b;
192 Standard_Real b_bprim = bprim*beta,
193 b2_bprim = bprim*beta2;
194 DVx(1) = 0; DVy(1) = 0;
195 DVx(2) = 0; DVy(2) = bprim;
196 DVx(3) = -2*b_bprim; DVy(3) = 0;
197 DVx(4) = 0; DVy(4) = -3*b2_bprim;
198 DVx(5) = -2*bprim*Sin_b; DVy(5) = 2*bprim*Cos_b;
199 DVx(6) = -bprim*Sin_b - 2*b_bprim*Cos_b;
200 DVy(6) = bprim*Cos_b - 2*b_bprim*Sin_b;
201 DVx(7) = 2*b_bprim*(-Cos_b + beta*Sin_b);
202 DVy(7) = -2*b_bprim*(Sin_b+beta*Cos_b);
203 DVx(8) = b2_bprim*(3*Sin_b + 2*beta*Cos_b);
204 DVy(8) = b2_bprim*(2*beta*Sin_b - 3*Cos_b);
213 for (ii=1; ii<=Ordre; ii++) {
214 P.SetCoord(Px(ii), Py(ii), 0);
215 Poles(ii).ChangeCoord() = M*P + Center.XYZ();
217 DP.SetCoord(DPx(ii), DPy(ii), 0);
219 aux.SetLinearForm(1, P, 1, DP, DCenter.XYZ());
220 DPoles(ii).SetXYZ(aux);
224 void GeomFill_PolynomialConvertor::Section(const gp_Pnt& FirstPnt,
225 const gp_Vec& DFirstPnt,
226 const gp_Vec& D2FirstPnt,
227 const gp_Pnt& Center,
228 const gp_Vec& DCenter,
229 const gp_Vec& D2Center,
233 const Standard_Real Angle,
234 const Standard_Real DAngle,
235 const Standard_Real D2Angle,
236 TColgp_Array1OfPnt& Poles,
237 TColgp_Array1OfVec& DPoles,
238 TColgp_Array1OfVec& D2Poles) const
240 math_Vector Vx(1, Ordre), Vy(1, Ordre),
241 DVx(1, Ordre), DVy(1, Ordre),
242 D2Vx(1, Ordre), D2Vy(1, Ordre);
243 math_Vector Px(1, Ordre), Py(1, Ordre),
244 DPx(1, Ordre), DPy(1, Ordre),
245 D2Px(1, Ordre), D2Py(1, Ordre);
248 Standard_Real aux, Cos_b = Cos(Angle), Sin_b = Sin(Angle);
249 Standard_Real beta, beta2, beta3, bprim, bprim2, bsecn;
250 gp_Vec V1(Center, FirstPnt), V1Prim, V1Secn, V2;
257 beta3 = beta * beta2;
258 bprim2 = bprim*bprim;
260 // Calcul des transformations
261 gp_Mat M (V1.X(), V2.X(), 0,
264 V1Prim = DFirstPnt - DCenter;
265 V2 = (DDir^V1) + (Dir^V1Prim);
267 gp_Mat MPrim (V1Prim.X(), V2.X(), 0,
268 V1Prim.Y(), V2.Y(), 0,
269 V1Prim.Z(), V2.Z(), 0);
270 V1Secn = D2FirstPnt - D2Center;
273 V2 += (D2Dir^V1) + (Dir^V1Secn);
276 gp_Mat MSecn (V1Secn.X(), V2.X(), 0,
277 V1Secn.Y(), V2.Y(), 0,
278 V1Secn.Z(), V2.Z(), 0);
281 // Calcul des contraintes -----------
282 Vx(1) = 1; Vy(1) = 0;
283 Vx(2) = 0; Vy(2) = beta;
284 Vx(3) = -beta2; Vy(3) = 0;
285 Vx(4) = 0; Vy(4) = -beta3;
286 Vx(5) = Cos_b; Vy(5) = Sin_b;
287 Vx(6) = -beta*Sin_b; Vy(6) = beta*Cos_b;
288 Vx(7) = -beta2*Cos_b; Vy(7) = -beta2*Sin_b;
289 Vx(8) = beta3*Sin_b; Vy(8) = -beta3*Cos_b;
291 Standard_Real b_bprim = bprim*beta,
292 b2_bprim = bprim*beta2,
293 b_bsecn = bsecn*beta;
294 DVx(1) = 0; DVy(1) = 0;
295 DVx(2) = 0; DVy(2) = bprim;
296 DVx(3) = -2*b_bprim; DVy(3) = 0;
297 DVx(4) = 0; DVy(4) = -3*b2_bprim;
298 DVx(5) = -2*bprim*Sin_b; DVy(5) = 2*bprim*Cos_b;
299 DVx(6) = -bprim*Sin_b - 2*b_bprim*Cos_b;
300 DVy(6) = bprim*Cos_b - 2*b_bprim*Sin_b;
301 DVx(7) = 2*b_bprim*(-Cos_b + beta*Sin_b);
302 DVy(7) = -2*b_bprim*(Sin_b + beta*Cos_b);
303 DVx(8) = b2_bprim*(3*Sin_b + 2*beta*Cos_b);
304 DVy(8) = b2_bprim*(2*beta*Sin_b - 3*Cos_b);
307 D2Vx(1) = 0; D2Vy(1) = 0;
308 D2Vx(2) = 0; D2Vy(2) = bsecn;
309 D2Vx(3) = -2*(bprim2+b_bsecn); D2Vy(3) = 0;
310 D2Vx(4) = 0; D2Vy(4) = -3*beta*(2*bprim2+b_bsecn);
311 D2Vx(5) = -2*(bsecn*Sin_b + 2*bprim2*Cos_b);
312 D2Vy(5) = 2*(bsecn*Cos_b - 2*bprim2*Sin_b);
313 D2Vx(6) = (4*beta*bprim2-bsecn)*Sin_b - 2*(2*bprim2+b_bsecn)*Cos_b;
315 D2Vy(6) = (bsecn - 4*beta*bprim2)*Cos_b
316 - 2*(b_bsecn + 2*bprim2)*Sin_b;
318 aux = 2*(bprim2+b_bsecn);
319 D2Vx(7) = aux*(-Cos_b + beta*Sin_b)
320 + 2*beta*bprim2*(2*beta*Cos_b + 3*Sin_b);
322 D2Vy(7) = -aux*(Sin_b + beta*Cos_b)
323 - 2*beta*bprim2*(3*Cos_b - 2*beta*Sin_b);
325 aux = beta*(2*bprim2+b_bsecn);
326 D2Vx(8) = aux * (3*Sin_b + 2*beta*Cos_b)
327 + 4*beta2*bprim2 * (2*Cos_b - beta*Sin_b);
329 D2Vy(8)= aux * (2*beta*Sin_b - 3*Cos_b)
330 + 4*beta2*bprim2 * (2*Sin_b + beta*Cos_b);
341 gp_XYZ P, DP, D2P, auxyz;
342 for (ii=1; ii<=Ordre; ii++) {
343 P.SetCoord(Px(ii), Py(ii), 0);
344 DP.SetCoord(DPx(ii), DPy(ii), 0);
345 D2P.SetCoord(D2Px(ii), D2Py(ii), 0);
347 Poles(ii).ChangeCoord() = M*P + Center.XYZ();
348 auxyz.SetLinearForm(1, MPrim*P,
351 DPoles(ii).SetXYZ(auxyz);
355 auxyz.SetLinearForm(1, P,
359 D2Poles(ii).SetXYZ(auxyz);