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
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 <math_Vector.hxx>
19 #include <PLib_DoubleJacobiPolynomial.hxx>
20 #include <PLib_JacobiPolynomial.hxx>
22 //=======================================================================
23 //function : PLib_DoubleJacobiPolynomial
25 //=======================================================================
26 PLib_DoubleJacobiPolynomial::PLib_DoubleJacobiPolynomial()
31 //=======================================================================
32 //function : PLib_DoubleJacobiPolynomial
34 //=======================================================================
36 PLib_DoubleJacobiPolynomial::PLib_DoubleJacobiPolynomial(const Handle(PLib_JacobiPolynomial)& JacPolU,
37 const Handle(PLib_JacobiPolynomial)& JacPolV) :
41 Handle (TColStd_HArray1OfReal) TabMaxU =
42 new TColStd_HArray1OfReal (0,JacPolU->WorkDegree()-2*(JacPolU->NivConstr()+1));
43 JacPolU->MaxValue(TabMaxU->ChangeArray1());
46 Handle (TColStd_HArray1OfReal) TabMaxV =
47 new TColStd_HArray1OfReal (0,JacPolV->WorkDegree()-2*(JacPolV->NivConstr()+1));
48 JacPolV->MaxValue(TabMaxV->ChangeArray1());
52 //=======================================================================
53 //function : MaxErrorU
55 //=======================================================================
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
64 Standard_Integer ii,idim,dJac,MinU,MinV,WorkDegreeU,WorkDegreeV;
67 math_Vector MaxErrDim(1,Dimension,0.);
69 MinU = 2*(myJacPolU->NivConstr()+1);
70 MinV = 2*(myJacPolV->NivConstr()+1);
71 WorkDegreeU = myJacPolU->WorkDegree();
72 WorkDegreeV = myJacPolV->WorkDegree();
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);
82 return (MaxErrDim.Norm());
85 //=======================================================================
86 //function : MaxErrorV
88 //=======================================================================
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
97 Standard_Integer jj,idim,dJac,MinU,MinV,WorkDegreeU,WorkDegreeV;
100 math_Vector MaxErrDim(1,Dimension,0.);
102 MinU = 2*(myJacPolU->NivConstr()+1);
103 MinV = 2*(myJacPolV->NivConstr()+1);
104 WorkDegreeU = myJacPolU->WorkDegree();
105 WorkDegreeV = myJacPolV->WorkDegree();
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);
115 return (MaxErrDim.Norm());
118 //=======================================================================
119 //function : MaxError
121 //=======================================================================
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
133 Standard_Integer ii,jj,idim,dJac,MinU,MinV,WorkDegreeU,WorkDegreeV;
134 Standard_Real Bid0,Bid1;
136 math_Vector MaxErrDim(1,Dimension,0.);
138 MinU = 2*(myJacPolU->NivConstr()+1);
139 MinV = 2*(myJacPolV->NivConstr()+1);
140 WorkDegreeU = myJacPolU->WorkDegree();
141 WorkDegreeV = myJacPolV->WorkDegree();
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 --------------
147 for (idim=1; idim<=Dimension; idim++) {
148 dJac = dJacCoeff + (idim-1)*(WorkDegreeU+1)*(WorkDegreeV+1);
150 for (jj=MinDegreeV; jj<=MaxDegreeV; jj++) {
152 for (ii=MinDegreeU; ii<=MaxDegreeU; ii++) {
153 Bid0 += fabs(JacCoeff(ii + jj*(WorkDegreeU+1) + dJac)) * myTabMaxU->Value(ii-MinU);
155 Bid1 += Bid0 * myTabMaxV->Value(jj-MinV);
157 MaxErrDim(idim) = Bid1;
160 //----------------------- Calcul de l' erreur max ----------------------
162 math_Vector MaxErr2(1,2);
164 MaxErr2(2) = MaxErrDim.Norm();
165 return (MaxErr2.Norm());
168 //=======================================================================
169 //function : ReduceDegree
171 //=======================================================================
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
185 Standard_Integer NewU,NewV;
186 Standard_Real ErrU,ErrV;
190 math_Vector MaxErr2(1,2);
192 //**********************************************************************
193 //-------------------- Coupure des coefficients ------------------------
194 //**********************************************************************
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);
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);
216 //----------------------- Calcul de l' erreur max ----------------------
217 MaxErr2(1) = MaxError;
219 ErrU = MaxErr2.Norm();
221 ErrV = MaxErr2.Norm();
224 if (ErrV < EpmsCut) {
230 if (ErrU < EpmsCut) {
236 while ((ErrU > ErrV && ErrV <= EpmsCut) || (ErrV >= ErrU && ErrU <= EpmsCut));
238 //-------------------------- Recuperation des degres -------------------
240 NewDegreeU = Max(NewU,1);
241 NewDegreeV = Max(NewV,1);
244 //=======================================================================
245 //function : AverageError
247 //=======================================================================
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
256 Standard_Integer ii,jj,idim,dJac,IDebU,IDebV,MinU,MinV,WorkDegreeU,WorkDegreeV;
257 Standard_Real Bid0,Bid1,AverageErr;
259 //----------------------------- Initialisations ------------------------
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();
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 --------------
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);
281 for (jj=IDebV; jj<=MinV-1; jj++) {
282 for (ii=MinU; ii<=WorkDegreeU; ii++) {
283 Bid1 = JacCoeff(ii + jj*(WorkDegreeU+1) + dJac);
288 AverageErr = sqrt(Bid0/4);
292 //=======================================================================
293 //function : WDoubleJacobiToCoefficients
295 //=======================================================================
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
303 Standard_Integer iu,iv,idim,WorkDegreeU,WorkDegreeV;
305 Coefficients.Init(0.);
307 WorkDegreeU = myJacPolU->WorkDegree();
308 WorkDegreeV = myJacPolV->WorkDegree();
310 TColStd_Array1OfReal Aux1(0, (DegreeU+1)*(DegreeV+1)*Dimension-1);
311 TColStd_Array1OfReal Aux2(0, (DegreeU+1)*(DegreeV+1)*Dimension-1);
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));
321 // Passage dans canonique en u.
322 myJacPolU->ToCoefficients(Dimension*(DegreeV+1),DegreeU,Aux1,Aux2);
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));
333 // Passage dans canonique en v.
334 myJacPolV->ToCoefficients(Dimension*(DegreeU+1),DegreeV,Aux1,Aux2);
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));