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