b311480e |
1 | // Created on: 1997-10-22 |
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 | |
7fd59977 |
17 | #include <PLib_HermitJacobi.ixx> |
18 | #include <PLib.hxx> |
f7b4312f |
19 | #include <NCollection_LocalArray.hxx> |
7fd59977 |
20 | #include <TColStd_HArray1OfReal.hxx> |
21 | |
22 | //======================================================================= |
23 | //function : PLib_HermitJacobi |
24 | //purpose : |
25 | //======================================================================= |
26 | |
27 | PLib_HermitJacobi::PLib_HermitJacobi(const Standard_Integer WorkDegree, |
28 | const GeomAbs_Shape ConstraintOrder) : |
29 | myH(1,2*(PLib::NivConstr(ConstraintOrder)+1), |
30 | 1,2*(PLib::NivConstr(ConstraintOrder)+1)), |
31 | myWCoeff(1,2*(PLib::NivConstr(ConstraintOrder)+1)+1) |
32 | { |
33 | Standard_Integer NivConstr = PLib::NivConstr(ConstraintOrder); |
34 | PLib::HermiteCoefficients(-1.,1.,NivConstr,NivConstr,myH); |
35 | |
36 | myJacobi = new PLib_JacobiPolynomial (WorkDegree,ConstraintOrder); |
37 | |
38 | myWCoeff.Init(0.); |
39 | myWCoeff(1) = 1.; |
40 | switch(NivConstr) { |
41 | case 0: myWCoeff(3) = -1.; break; |
42 | case 1: myWCoeff(3) = -2.; myWCoeff(5) = 1.; break; |
43 | case 2: myWCoeff(3) = -3.; myWCoeff(5) = 3.; myWCoeff(7) = -1.; break; |
44 | } |
45 | } |
46 | |
47 | //======================================================================= |
48 | //function : MaxError |
49 | //purpose : |
50 | //======================================================================= |
51 | |
52 | Standard_Real PLib_HermitJacobi::MaxError(const Standard_Integer Dimension, |
53 | Standard_Real& HermJacCoeff, |
54 | const Standard_Integer NewDegree) const |
55 | { |
56 | return myJacobi->MaxError(Dimension,HermJacCoeff,NewDegree); |
57 | } |
58 | |
59 | //======================================================================= |
60 | //function : ReduceDegree |
61 | //purpose : |
62 | //======================================================================= |
63 | |
64 | void PLib_HermitJacobi::ReduceDegree(const Standard_Integer Dimension, |
65 | const Standard_Integer MaxDegree, |
66 | const Standard_Real Tol, |
67 | Standard_Real& HermJacCoeff, |
68 | Standard_Integer& NewDegree, |
69 | Standard_Real& MaxError) const |
70 | { |
71 | myJacobi->ReduceDegree(Dimension,MaxDegree,Tol, |
72 | HermJacCoeff,NewDegree,MaxError); |
73 | } |
74 | |
75 | //======================================================================= |
76 | //function : AverageError |
77 | //purpose : |
78 | //======================================================================= |
79 | |
80 | Standard_Real PLib_HermitJacobi::AverageError(const Standard_Integer Dimension, |
81 | Standard_Real& HermJacCoeff, |
82 | const Standard_Integer NewDegree) const |
83 | { |
84 | return myJacobi->AverageError(Dimension,HermJacCoeff,NewDegree); |
85 | } |
86 | |
87 | //======================================================================= |
88 | //function : ToCoefficients |
89 | //purpose : |
90 | //======================================================================= |
91 | |
92 | void PLib_HermitJacobi::ToCoefficients(const Standard_Integer Dimension, |
93 | const Standard_Integer Degree, |
94 | const TColStd_Array1OfReal& HermJacCoeff, |
95 | TColStd_Array1OfReal& Coefficients) const |
96 | { |
97 | Standard_Integer i,k,idim,i1,i2; |
98 | Standard_Real h1, h2; |
99 | Standard_Integer NivConstr = this->NivConstr(), |
100 | DegreeH = 2*NivConstr+1; |
101 | Standard_Integer ibegHJC = HermJacCoeff.Lower(), kdim; |
102 | |
103 | TColStd_Array1OfReal AuxCoeff(0,(Degree+1)*Dimension-1); |
104 | AuxCoeff.Init(0.); |
105 | |
106 | for (k=0; k<=DegreeH; k++) { |
107 | kdim = k*Dimension; |
108 | for (i=0; i<=NivConstr; i++) { |
109 | h1 = myH(i+1, k+1); |
110 | h2 = myH(i+NivConstr+2,k+1); |
111 | i1 = ibegHJC + i*Dimension; |
112 | i2 = ibegHJC + (i+NivConstr+1)*Dimension; |
113 | |
114 | for (idim=0; idim<Dimension; idim++) { |
115 | AuxCoeff(idim + kdim) += HermJacCoeff(i1 + idim) * h1 + |
116 | HermJacCoeff(i2 + idim) * h2; |
117 | } |
118 | } |
119 | } |
120 | kdim = (Degree+1)*Dimension; |
121 | for (k=(DegreeH+1)*Dimension; k<kdim; k++) { |
122 | AuxCoeff(k) = HermJacCoeff(ibegHJC + k); |
123 | } |
124 | |
125 | if(Degree > DegreeH) |
126 | myJacobi->ToCoefficients(Dimension,Degree,AuxCoeff,Coefficients); |
127 | else { |
128 | Standard_Integer ibegC = Coefficients.Lower(); |
129 | kdim = (Degree+1)*Dimension; |
130 | for(k=0; k < kdim; k++) |
131 | Coefficients(ibegC+k) = AuxCoeff(k); |
132 | } |
133 | } |
134 | |
135 | //======================================================================= |
136 | //function : D0123 |
137 | //purpose : common part of D0,D1,D2,D3 (FORTRAN subroutine MPOBAS) |
138 | //======================================================================= |
139 | |
140 | void PLib_HermitJacobi::D0123(const Standard_Integer NDeriv, |
141 | const Standard_Real U, |
142 | TColStd_Array1OfReal& BasisValue, |
143 | TColStd_Array1OfReal& BasisD1, |
144 | TColStd_Array1OfReal& BasisD2, |
145 | TColStd_Array1OfReal& BasisD3) |
146 | { |
f7b4312f |
147 | NCollection_LocalArray<Standard_Real> jac0 (4 * 20); |
148 | NCollection_LocalArray<Standard_Real> jac1 (4 * 20); |
149 | NCollection_LocalArray<Standard_Real> jac2 (4 * 20); |
150 | NCollection_LocalArray<Standard_Real> jac3 (4 * 20); |
151 | NCollection_LocalArray<Standard_Real> wvalues (4); |
7fd59977 |
152 | |
153 | Standard_Integer i, j; |
154 | Standard_Integer NivConstr = this->NivConstr(), |
155 | WorkDegree = this->WorkDegree(), |
156 | DegreeH = 2*NivConstr+1; |
157 | Standard_Integer ibeg0 = BasisValue.Lower(), |
158 | ibeg1 = BasisD1.Lower(), |
159 | ibeg2 = BasisD2.Lower(), |
160 | ibeg3 = BasisD3.Lower(); |
161 | Standard_Integer JacDegree = WorkDegree-DegreeH-1; |
162 | Standard_Real W0; |
163 | |
164 | TColStd_Array1OfReal JacValue0(jac0[0], 0, Max(0,JacDegree)); |
165 | TColStd_Array1OfReal WValues(wvalues[0],0,NDeriv); |
166 | WValues.Init(0.); |
167 | |
168 | // Evaluation des polynomes d'hermite |
169 | math_Matrix HermitValues(0,DegreeH, 0, NDeriv, 0.); |
170 | if(NDeriv == 0) |
171 | for (i=0; i<=DegreeH; i++) { |
172 | PLib::NoDerivativeEvalPolynomial(U,DegreeH,1, DegreeH, |
173 | myH(i+1,1), HermitValues(i,0)); |
174 | } |
175 | else |
176 | for (i=0; i<=DegreeH; i++) { |
177 | PLib::EvalPolynomial(U,NDeriv,DegreeH,1, |
178 | myH(i+1,1), HermitValues(i,0)); |
179 | } |
180 | |
181 | // Evaluation des polynomes de Jaccobi |
182 | if(JacDegree >= 0) { |
183 | |
184 | switch (NDeriv) { |
185 | case 0 : |
186 | myJacobi->D0(U,JacValue0); |
187 | break; |
188 | case 1 : |
189 | { |
190 | TColStd_Array1OfReal JacValue1(jac1[0], 0, JacDegree); |
191 | myJacobi->D1(U,JacValue0,JacValue1); |
192 | break; |
193 | } |
194 | case 2 : |
195 | { |
196 | TColStd_Array1OfReal JacValue1(jac1[0], 0, JacDegree); |
197 | TColStd_Array1OfReal JacValue2(jac2[0], 0, JacDegree); |
198 | myJacobi->D2(U,JacValue0,JacValue1,JacValue2); |
199 | break; |
200 | } |
201 | case 3 : |
202 | { |
203 | TColStd_Array1OfReal JacValue1(jac1[0], 0, JacDegree); |
204 | TColStd_Array1OfReal JacValue2(jac2[0], 0, JacDegree); |
205 | TColStd_Array1OfReal JacValue3(jac3[0], 0, JacDegree); |
206 | myJacobi->D3(U,JacValue0,JacValue1,JacValue2,JacValue3); |
207 | } |
208 | } |
209 | |
210 | // Evaluation de W(t) |
211 | if(NDeriv == 0) |
212 | PLib::NoDerivativeEvalPolynomial(U,DegreeH+1,1,DegreeH+1,myWCoeff(1), WValues(0)); |
213 | else |
214 | PLib::EvalPolynomial(U,NDeriv,DegreeH+1,1,myWCoeff(1), WValues(0)); |
215 | } |
216 | |
217 | // Evaluation a l'ordre 0 |
218 | for (i=0; i<=DegreeH; i++) { |
219 | BasisValue(ibeg0+i) = HermitValues(i,0); |
220 | } |
221 | W0 = WValues(0); |
222 | for (i=DegreeH+1, j=0; i<=WorkDegree; i++, j++) { |
223 | BasisValue(ibeg0+i) = W0 * jac0[j]; |
224 | } |
225 | |
226 | // Evaluation a l'ordre 1 |
227 | if (NDeriv >= 1) { |
228 | Standard_Real W1=WValues(1); |
229 | for (i=0; i<=DegreeH; i++) { |
230 | BasisD1(ibeg1+i) = HermitValues(i,1); |
231 | } |
232 | for (i=DegreeH+1, j=0; i<=WorkDegree; i++, j++) { |
233 | BasisD1(ibeg1+i) = W0 * jac1[j] + |
234 | W1 * jac0[j]; |
235 | } |
236 | // Evaluation a l'ordre 2 |
237 | if (NDeriv >= 2) { |
238 | Standard_Real W2=WValues(2); |
239 | for (i=0; i<=DegreeH; i++) { |
240 | BasisD2(ibeg2+i) = HermitValues(i,2); |
241 | } |
242 | for (i=DegreeH+1, j=0; i<=WorkDegree; i++, j++) { |
243 | BasisD2(ibeg2+i) = |
244 | W0 * jac2[j] + 2 * W1 * jac1[j] + W2 * jac0[j]; |
245 | } |
246 | |
247 | // Evaluation a l'ordre 3 |
248 | if (NDeriv == 3) { |
249 | Standard_Real W3 = WValues(3); |
250 | for (i=0; i<=DegreeH; i++) { |
251 | BasisD3(ibeg3+i) = HermitValues(i,3); |
252 | } |
253 | for (i=DegreeH+1, j=0; i<=WorkDegree; i++,j++) { |
254 | BasisD3(ibeg3+i) = W0*jac3[j] + W3*jac0[j] |
255 | + 3*(W1*jac2[j] + W2*jac1[j]); |
256 | } |
257 | } |
258 | } |
259 | } |
260 | } |
261 | |
262 | //======================================================================= |
263 | //function : D0 |
264 | //purpose : |
265 | //======================================================================= |
266 | |
267 | void PLib_HermitJacobi::D0(const Standard_Real U, TColStd_Array1OfReal& BasisValue) |
268 | { |
269 | D0123(0,U,BasisValue,BasisValue,BasisValue,BasisValue); |
270 | } |
271 | |
272 | //======================================================================= |
273 | //function : D1 |
274 | //purpose : |
275 | //======================================================================= |
276 | |
277 | void PLib_HermitJacobi::D1(const Standard_Real U, |
278 | TColStd_Array1OfReal& BasisValue, TColStd_Array1OfReal& BasisD1) |
279 | { |
280 | D0123(1,U,BasisValue,BasisD1,BasisD1,BasisD1); |
281 | } |
282 | |
283 | //======================================================================= |
284 | //function : D2 |
285 | //purpose : |
286 | //======================================================================= |
287 | |
288 | void PLib_HermitJacobi::D2(const Standard_Real U, TColStd_Array1OfReal& BasisValue, |
289 | TColStd_Array1OfReal& BasisD1,TColStd_Array1OfReal& BasisD2) |
290 | { |
291 | D0123(2,U,BasisValue,BasisD1,BasisD2,BasisD2); |
292 | } |
293 | |
294 | //======================================================================= |
295 | //function : D3 |
296 | //purpose : |
297 | //======================================================================= |
298 | |
299 | void PLib_HermitJacobi::D3(const Standard_Real U, |
300 | TColStd_Array1OfReal& BasisValue, TColStd_Array1OfReal& BasisD1, |
301 | TColStd_Array1OfReal& BasisD2, TColStd_Array1OfReal& BasisD3) |
302 | { |
303 | D0123(3,U,BasisValue,BasisD1,BasisD2,BasisD3); |
304 | } |