b311480e |
1 | // Created on: 1997-05-28 |
2 | // Created by: Sergey SOKOLOV |
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 | |
42cf5bc1 |
17 | |
7fd59977 |
18 | #include <math_Vector.hxx> |
42cf5bc1 |
19 | #include <PLib_DoubleJacobiPolynomial.hxx> |
20 | #include <PLib_JacobiPolynomial.hxx> |
7fd59977 |
21 | |
22 | //======================================================================= |
23 | //function : PLib_DoubleJacobiPolynomial |
24 | //purpose : |
25 | //======================================================================= |
b311480e |
26 | PLib_DoubleJacobiPolynomial::PLib_DoubleJacobiPolynomial() |
7fd59977 |
27 | |
28 | { |
29 | } |
30 | |
31 | //======================================================================= |
32 | //function : PLib_DoubleJacobiPolynomial |
33 | //purpose : |
34 | //======================================================================= |
35 | |
36 | PLib_DoubleJacobiPolynomial::PLib_DoubleJacobiPolynomial(const Handle(PLib_JacobiPolynomial)& JacPolU, |
37 | const Handle(PLib_JacobiPolynomial)& JacPolV) : |
38 | myJacPolU(JacPolU), |
39 | myJacPolV(JacPolV) |
40 | { |
41 | Handle (TColStd_HArray1OfReal) TabMaxU = |
42 | new TColStd_HArray1OfReal (0,JacPolU->WorkDegree()-2*(JacPolU->NivConstr()+1)); |
43 | JacPolU->MaxValue(TabMaxU->ChangeArray1()); |
44 | myTabMaxU = TabMaxU; |
45 | |
46 | Handle (TColStd_HArray1OfReal) TabMaxV = |
47 | new TColStd_HArray1OfReal (0,JacPolV->WorkDegree()-2*(JacPolV->NivConstr()+1)); |
48 | JacPolV->MaxValue(TabMaxV->ChangeArray1()); |
49 | myTabMaxV = TabMaxV; |
50 | } |
51 | |
52 | //======================================================================= |
53 | //function : MaxErrorU |
54 | //purpose : |
55 | //======================================================================= |
56 | |
57 | Standard_Real |
58 | PLib_DoubleJacobiPolynomial::MaxErrorU(const Standard_Integer Dimension, |
59 | const Standard_Integer DegreeU, |
60 | const Standard_Integer DegreeV, |
61 | const Standard_Integer dJacCoeff, |
62 | const TColStd_Array1OfReal& JacCoeff) const |
63 | { |
64 | Standard_Integer ii,idim,dJac,MinU,MinV,WorkDegreeU,WorkDegreeV; |
65 | Standard_Real Bid0; |
66 | |
67 | math_Vector MaxErrDim(1,Dimension,0.); |
68 | |
69 | MinU = 2*(myJacPolU->NivConstr()+1); |
70 | MinV = 2*(myJacPolV->NivConstr()+1); |
71 | WorkDegreeU = myJacPolU->WorkDegree(); |
72 | WorkDegreeV = myJacPolV->WorkDegree(); |
73 | |
74 | Bid0 = myTabMaxV->Value(DegreeV-MinV); |
75 | for (idim=1; idim<=Dimension; idim++) { |
76 | dJac = dJacCoeff + (idim-1)*(WorkDegreeU+1)*(WorkDegreeV+1); |
77 | for (ii=MinU; ii<=DegreeU; ii++) { |
78 | MaxErrDim(idim) += (Abs(JacCoeff(ii + DegreeV*(WorkDegreeU+1) + dJac)) * |
79 | myTabMaxU->Value(ii-MinU) * Bid0); |
80 | } |
81 | } |
82 | return (MaxErrDim.Norm()); |
83 | } |
84 | |
85 | //======================================================================= |
86 | //function : MaxErrorV |
87 | //purpose : |
88 | //======================================================================= |
89 | |
90 | Standard_Real |
91 | PLib_DoubleJacobiPolynomial::MaxErrorV(const Standard_Integer Dimension, |
92 | const Standard_Integer DegreeU, |
93 | const Standard_Integer DegreeV, |
94 | const Standard_Integer dJacCoeff, |
95 | const TColStd_Array1OfReal& JacCoeff) const |
96 | { |
97 | Standard_Integer jj,idim,dJac,MinU,MinV,WorkDegreeU,WorkDegreeV; |
98 | Standard_Real Bid0; |
99 | |
100 | math_Vector MaxErrDim(1,Dimension,0.); |
101 | |
102 | MinU = 2*(myJacPolU->NivConstr()+1); |
103 | MinV = 2*(myJacPolV->NivConstr()+1); |
104 | WorkDegreeU = myJacPolU->WorkDegree(); |
105 | WorkDegreeV = myJacPolV->WorkDegree(); |
106 | |
107 | Bid0 = myTabMaxU->Value(DegreeU-MinU); |
108 | for (idim=1; idim<=Dimension; idim++) { |
109 | dJac = dJacCoeff + (idim-1)*(WorkDegreeU+1)*(WorkDegreeV+1); |
110 | for (jj=MinV; jj<=DegreeV; jj++) { |
111 | MaxErrDim(idim) += (Abs(JacCoeff(DegreeU + jj*(WorkDegreeU+1) + dJac)) * |
112 | myTabMaxV->Value(jj-MinV) * Bid0); |
113 | } |
114 | } |
115 | return (MaxErrDim.Norm()); |
116 | } |
117 | |
118 | //======================================================================= |
119 | //function : MaxError |
120 | //purpose : |
121 | //======================================================================= |
122 | |
123 | Standard_Real |
124 | PLib_DoubleJacobiPolynomial::MaxError(const Standard_Integer Dimension, |
125 | const Standard_Integer MinDegreeU, |
126 | const Standard_Integer MaxDegreeU, |
127 | const Standard_Integer MinDegreeV, |
128 | const Standard_Integer MaxDegreeV, |
129 | const Standard_Integer dJacCoeff, |
130 | const TColStd_Array1OfReal& JacCoeff, |
131 | const Standard_Real Error) const |
132 | { |
133 | Standard_Integer ii,jj,idim,dJac,MinU,MinV,WorkDegreeU,WorkDegreeV; |
134 | Standard_Real Bid0,Bid1; |
135 | |
136 | math_Vector MaxErrDim(1,Dimension,0.); |
137 | |
138 | MinU = 2*(myJacPolU->NivConstr()+1); |
139 | MinV = 2*(myJacPolV->NivConstr()+1); |
140 | WorkDegreeU = myJacPolU->WorkDegree(); |
141 | WorkDegreeV = myJacPolV->WorkDegree(); |
142 | |
143 | //------------------- Calcul du majorant de l'erreur max --------------- |
144 | //----- lorsque sont enleves les coeff. d'indices MinDegreeU a MaxDegreeU ------ |
145 | //---------------- en U et d'indices MinDegreeV a MaxDegreeV en V -------------- |
146 | |
147 | for (idim=1; idim<=Dimension; idim++) { |
148 | dJac = dJacCoeff + (idim-1)*(WorkDegreeU+1)*(WorkDegreeV+1); |
149 | Bid1 = 0.; |
150 | for (jj=MinDegreeV; jj<=MaxDegreeV; jj++) { |
151 | Bid0 = 0.; |
152 | for (ii=MinDegreeU; ii<=MaxDegreeU; ii++) { |
153 | Bid0 += fabs(JacCoeff(ii + jj*(WorkDegreeU+1) + dJac)) * myTabMaxU->Value(ii-MinU); |
154 | } |
155 | Bid1 += Bid0 * myTabMaxV->Value(jj-MinV); |
156 | } |
157 | MaxErrDim(idim) = Bid1; |
158 | } |
159 | |
160 | //----------------------- Calcul de l' erreur max ---------------------- |
161 | |
162 | math_Vector MaxErr2(1,2); |
163 | MaxErr2(1) = Error; |
164 | MaxErr2(2) = MaxErrDim.Norm(); |
165 | return (MaxErr2.Norm()); |
166 | } |
167 | |
168 | //======================================================================= |
169 | //function : ReduceDegree |
170 | //purpose : |
171 | //======================================================================= |
172 | |
173 | void PLib_DoubleJacobiPolynomial::ReduceDegree(const Standard_Integer Dimension, |
174 | const Standard_Integer MinDegreeU, |
175 | const Standard_Integer MaxDegreeU, |
176 | const Standard_Integer MinDegreeV, |
177 | const Standard_Integer MaxDegreeV, |
178 | const Standard_Integer dJacCoeff, |
179 | const TColStd_Array1OfReal& JacCoeff, |
180 | const Standard_Real EpmsCut, |
181 | Standard_Real& MaxError, |
182 | Standard_Integer& NewDegreeU, |
183 | Standard_Integer& NewDegreeV) const |
184 | { |
185 | Standard_Integer NewU,NewV; |
186 | Standard_Real ErrU,ErrV; |
187 | |
188 | NewU = MaxDegreeU; |
189 | NewV = MaxDegreeV; |
190 | math_Vector MaxErr2(1,2); |
191 | |
192 | //********************************************************************** |
193 | //-------------------- Coupure des coefficients ------------------------ |
194 | //********************************************************************** |
195 | |
196 | do { |
197 | |
198 | //------------------- Calcul du majorant de l'erreur max --------------- |
199 | //----- lorsque sont enleves les coeff. d'indices MinU a NewU ------ |
200 | //---------------- en U, le degre en V etant fixe a NewV ----------------- |
201 | if (NewV > MinDegreeV) |
202 | ErrV = MaxErrorU(Dimension,NewU,NewV,dJacCoeff,JacCoeff); |
203 | else { |
204 | ErrV = 2*EpmsCut; |
205 | } |
206 | |
207 | //------------------- Calcul du majorant de l'erreur max --------------- |
208 | //----- lorsque sont enleves les coeff. d'indices MinV a NewV ------ |
209 | //---------------- en V, le degre en U etant fixe a NewU ----------------- |
210 | if (NewU > MinDegreeU) |
211 | ErrU = MaxErrorV(Dimension,NewU,NewV,dJacCoeff,JacCoeff); |
212 | else { |
213 | ErrU = 2*EpmsCut; |
214 | } |
215 | |
216 | //----------------------- Calcul de l' erreur max ---------------------- |
217 | MaxErr2(1) = MaxError; |
218 | MaxErr2(2) = ErrU; |
219 | ErrU = MaxErr2.Norm(); |
220 | MaxErr2(2) = ErrV; |
221 | ErrV = MaxErr2.Norm(); |
222 | |
223 | if (ErrU > ErrV) { |
224 | if (ErrV < EpmsCut) { |
225 | MaxError = ErrV; |
226 | NewV--; |
227 | } |
228 | } |
229 | else { |
230 | if (ErrU < EpmsCut) { |
231 | MaxError = ErrU; |
232 | NewU--; |
233 | } |
234 | } |
235 | } |
0ebaa4db |
236 | while ((ErrU > ErrV && ErrV <= EpmsCut) || (ErrV >= ErrU && ErrU <= EpmsCut)); |
7fd59977 |
237 | |
238 | //-------------------------- Recuperation des degres ------------------- |
239 | |
240 | NewDegreeU = Max(NewU,1); |
241 | NewDegreeV = Max(NewV,1); |
242 | } |
243 | |
244 | //======================================================================= |
245 | //function : AverageError |
246 | //purpose : |
247 | //======================================================================= |
248 | |
249 | Standard_Real |
250 | PLib_DoubleJacobiPolynomial::AverageError(const Standard_Integer Dimension, |
251 | const Standard_Integer DegreeU, |
252 | const Standard_Integer DegreeV, |
253 | const Standard_Integer dJacCoeff, |
254 | const TColStd_Array1OfReal& JacCoeff) const |
255 | { |
256 | Standard_Integer ii,jj,idim,dJac,IDebU,IDebV,MinU,MinV,WorkDegreeU,WorkDegreeV; |
257 | Standard_Real Bid0,Bid1,AverageErr; |
258 | |
259 | //----------------------------- Initialisations ------------------------ |
260 | |
261 | IDebU = 2*(myJacPolU->NivConstr()+1); |
262 | IDebV = 2*(myJacPolV->NivConstr()+1); |
263 | MinU = Max(IDebU,DegreeU); |
264 | MinV = Max(IDebV,DegreeV); |
265 | WorkDegreeU = myJacPolU->WorkDegree(); |
266 | WorkDegreeV = myJacPolV->WorkDegree(); |
267 | Bid0 = 0.; |
268 | |
269 | //------------------ Calcul du majorant de l'erreur moyenne ------------ |
270 | //----- lorsque sont enleves les coeff. d'indices DegreeU a WorkDegreeU ------ |
271 | //---------------- en U et d'indices DegreeV a WorkDegreeV en V -------------- |
272 | |
273 | for (idim=1; idim<=Dimension; idim++) { |
274 | dJac = dJacCoeff + (idim-1)*(WorkDegreeU+1)*(WorkDegreeV+1); |
275 | for (jj=MinV; jj<=WorkDegreeV; jj++) { |
276 | for (ii=IDebU; ii<=WorkDegreeU; ii++) { |
277 | Bid1 = JacCoeff(ii + jj*(WorkDegreeU+1) + dJac); |
278 | Bid0 += Bid1*Bid1; |
279 | } |
280 | } |
281 | for (jj=IDebV; jj<=MinV-1; jj++) { |
282 | for (ii=MinU; ii<=WorkDegreeU; ii++) { |
283 | Bid1 = JacCoeff(ii + jj*(WorkDegreeU+1) + dJac); |
284 | Bid0 += Bid1*Bid1; |
285 | } |
286 | } |
287 | } |
288 | AverageErr = sqrt(Bid0/4); |
289 | return (AverageErr); |
290 | } |
291 | |
292 | //======================================================================= |
293 | //function : WDoubleJacobiToCoefficients |
294 | //purpose : |
295 | //======================================================================= |
296 | |
297 | void PLib_DoubleJacobiPolynomial::WDoubleJacobiToCoefficients(const Standard_Integer Dimension, |
298 | const Standard_Integer DegreeU, |
299 | const Standard_Integer DegreeV, |
300 | const TColStd_Array1OfReal& JacCoeff, |
301 | TColStd_Array1OfReal& Coefficients) const |
302 | { |
303 | Standard_Integer iu,iv,idim,WorkDegreeU,WorkDegreeV; |
304 | |
305 | Coefficients.Init(0.); |
306 | |
307 | WorkDegreeU = myJacPolU->WorkDegree(); |
308 | WorkDegreeV = myJacPolV->WorkDegree(); |
309 | |
310 | TColStd_Array1OfReal Aux1(0, (DegreeU+1)*(DegreeV+1)*Dimension-1); |
311 | TColStd_Array1OfReal Aux2(0, (DegreeU+1)*(DegreeV+1)*Dimension-1); |
312 | |
313 | for (iu=0; iu<=DegreeU; iu++) { |
314 | for (iv=0; iv<=DegreeV; iv++) { |
315 | for (idim=1; idim<=Dimension; idim++) { |
316 | Aux1(idim-1 + iv*Dimension + iu*Dimension*(DegreeV+1)) = |
317 | JacCoeff(iu + iv*(WorkDegreeU+1) + (idim-1)*(WorkDegreeU+1)*(WorkDegreeV+1)); |
318 | } |
319 | } |
320 | } |
321 | // Passage dans canonique en u. |
322 | myJacPolU->ToCoefficients(Dimension*(DegreeV+1),DegreeU,Aux1,Aux2); |
323 | |
324 | // Permutation des u et des v. |
325 | for (iu=0; iu<=DegreeU; iu++) { |
326 | for (iv=0; iv<=DegreeV; iv++) { |
327 | for (idim=1; idim<=Dimension; idim++) { |
328 | Aux1(idim-1 + iu*Dimension + iv*Dimension*(DegreeU+1)) = |
329 | Aux2(idim-1 + iv*Dimension + iu*Dimension*(DegreeV+1)); |
330 | } |
331 | } |
332 | } |
333 | // Passage dans canonique en v. |
334 | myJacPolV->ToCoefficients(Dimension*(DegreeU+1),DegreeV,Aux1,Aux2); |
335 | |
336 | // Permutation des u et des v. |
337 | for (iu=0; iu<=DegreeU; iu++) { |
338 | for (iv=0; iv<=DegreeV; iv++) { |
339 | for (idim=1; idim<=Dimension; idim++) { |
340 | Coefficients(iu + iv*(DegreeU+1) + (idim-1)*(DegreeU+1)*(DegreeV+1)) = |
341 | Aux2(idim-1 + iu*Dimension + iv*Dimension*(DegreeU+1)); |
342 | } |
343 | } |
344 | } |
345 | } |
346 | |