1 // Created on: 1997-05-28
2 // Created by: Sergey SOKOLOV
3 // Copyright (c) 1997-1999 Matra Datavision
4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
5 //
6 // This file is part of Open CASCADE Technology software library.
7 //
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.
13 //
14 // Alternatively, this file may be used under the terms of Open CASCADE
15 // commercial license or contractual agreement.
16
17
18 #include <math_Vector.hxx>
19 #include <PLib_DoubleJacobiPolynomial.hxx>
20 #include <PLib_JacobiPolynomial.hxx>
21
22 //=======================================================================
23 //function : PLib_DoubleJacobiPolynomial
24 //purpose :
25 //=======================================================================
PLib_DoubleJacobiPolynomial()26 PLib_DoubleJacobiPolynomial::PLib_DoubleJacobiPolynomial()
27
28 {
29 }
30
31 //=======================================================================
32 //function : PLib_DoubleJacobiPolynomial
33 //purpose :
34 //=======================================================================
35
PLib_DoubleJacobiPolynomial(const Handle (PLib_JacobiPolynomial)& JacPolU,const Handle (PLib_JacobiPolynomial)& JacPolV)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
MaxErrorU(const Standard_Integer Dimension,const Standard_Integer DegreeU,const Standard_Integer DegreeV,const Standard_Integer dJacCoeff,const TColStd_Array1OfReal & JacCoeff) const58 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
MaxErrorV(const Standard_Integer Dimension,const Standard_Integer DegreeU,const Standard_Integer DegreeV,const Standard_Integer dJacCoeff,const TColStd_Array1OfReal & JacCoeff) const91 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
MaxError(const Standard_Integer Dimension,const Standard_Integer MinDegreeU,const Standard_Integer MaxDegreeU,const Standard_Integer MinDegreeV,const Standard_Integer MaxDegreeV,const Standard_Integer dJacCoeff,const TColStd_Array1OfReal & JacCoeff,const Standard_Real Error) const124 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
ReduceDegree(const Standard_Integer Dimension,const Standard_Integer MinDegreeU,const Standard_Integer MaxDegreeU,const Standard_Integer MinDegreeV,const Standard_Integer MaxDegreeV,const Standard_Integer dJacCoeff,const TColStd_Array1OfReal & JacCoeff,const Standard_Real EpmsCut,Standard_Real & MaxError,Standard_Integer & NewDegreeU,Standard_Integer & NewDegreeV) const173 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 }
236 while ((ErrU > ErrV && ErrV <= EpmsCut) || (ErrV >= ErrU && ErrU <= EpmsCut));
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
AverageError(const Standard_Integer Dimension,const Standard_Integer DegreeU,const Standard_Integer DegreeV,const Standard_Integer dJacCoeff,const TColStd_Array1OfReal & JacCoeff) const250 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
WDoubleJacobiToCoefficients(const Standard_Integer Dimension,const Standard_Integer DegreeU,const Standard_Integer DegreeV,const TColStd_Array1OfReal & JacCoeff,TColStd_Array1OfReal & Coefficients) const297 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
347