b311480e |
1 | // Copyright (c) 1997-1999 Matra Datavision |
973c2be1 |
2 | // Copyright (c) 1999-2014 OPEN CASCADE SAS |
b311480e |
3 | // |
973c2be1 |
4 | // This file is part of Open CASCADE Technology software library. |
b311480e |
5 | // |
973c2be1 |
6 | // This library is free software; you can redistribute it and / or modify it |
7 | // under the terms of the GNU Lesser General Public version 2.1 as published |
8 | // by the Free Software Foundation, with special exception defined in the file |
9 | // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT |
10 | // distribution for complete text of the license and disclaimer of any warranty. |
b311480e |
11 | // |
973c2be1 |
12 | // Alternatively, this file may be used under the terms of Open CASCADE |
13 | // commercial license or contractual agreement. |
b311480e |
14 | |
7fd59977 |
15 | #include <PLib_JacobiPolynomial.ixx> |
16 | |
17 | #include <math.hxx> |
18 | #include <math_Vector.hxx> |
19 | #include <TColStd_Array2OfReal.hxx> |
20 | #include <PLib.hxx> |
21 | #include <Standard_ConstructionError.hxx> |
22 | |
23 | #include <PLib_JacobiPolynomial_0.hxx> |
24 | |
25 | // The possible values for NbGaussPoints |
26 | const Standard_Integer NDEG8=8, NDEG10=10, NDEG15=15, NDEG20=20, NDEG25=25, |
27 | NDEG30=30, NDEG40=40, NDEG50=50, NDEG61=61; |
28 | |
29 | const Standard_Integer UNDEFINED=-999; |
30 | |
31 | //======================================================================= |
32 | //function : PLib_JacobiPolynomial |
33 | //purpose : |
34 | //======================================================================= |
35 | |
36 | PLib_JacobiPolynomial::PLib_JacobiPolynomial (const Standard_Integer WorkDegree, |
37 | const GeomAbs_Shape ConstraintOrder) |
38 | { |
39 | myWorkDegree = WorkDegree; |
40 | |
41 | switch (ConstraintOrder) { |
42 | case GeomAbs_C0: myNivConstr = 0; break; |
43 | case GeomAbs_C1: myNivConstr = 1; break; |
44 | case GeomAbs_C2: myNivConstr = 2; break; |
45 | default: |
46 | Standard_ConstructionError::Raise("Invalid ConstraintOrder"); |
47 | } |
48 | myDegree = myWorkDegree - 2*(myNivConstr+1); |
49 | if (myDegree > 30) |
50 | Standard_ConstructionError::Raise("Invalid Degree"); |
51 | } |
52 | |
53 | //======================================================================= |
54 | //function : Points |
55 | //purpose : |
56 | //======================================================================= |
57 | |
58 | void PLib_JacobiPolynomial::Points(const Standard_Integer NbGaussPoints, |
59 | TColStd_Array1OfReal& TabPoints) const |
60 | { |
0ebaa4db |
61 | if ((NbGaussPoints != NDEG8 && NbGaussPoints != NDEG10 && |
7fd59977 |
62 | NbGaussPoints != NDEG15 && NbGaussPoints != NDEG20 && |
63 | NbGaussPoints != NDEG25 && NbGaussPoints != NDEG30 && |
64 | NbGaussPoints != NDEG40 && NbGaussPoints != NDEG50 && |
0ebaa4db |
65 | NbGaussPoints != NDEG61) || |
7fd59977 |
66 | NbGaussPoints <= myDegree) |
67 | Standard_ConstructionError::Raise("Invalid NbGaussPoints"); |
68 | |
69 | math_Vector DecreasingPoints(1,NbGaussPoints); |
70 | |
71 | math::GaussPoints(NbGaussPoints,DecreasingPoints); |
72 | |
73 | // TabPoints consist of only positive increasing values |
74 | for (Standard_Integer i=1; i<=NbGaussPoints/2; i++) |
75 | TabPoints(i) = DecreasingPoints(NbGaussPoints/2-i+1); |
76 | if (NbGaussPoints % 2 == 1) |
77 | TabPoints(0) = 0.; |
78 | else |
79 | TabPoints(0) = UNDEFINED; |
80 | } |
81 | |
82 | //======================================================================= |
83 | //function : Weights |
84 | //purpose : |
85 | //======================================================================= |
86 | |
87 | void PLib_JacobiPolynomial::Weights(const Standard_Integer NbGaussPoints, |
88 | TColStd_Array2OfReal& TabWeights) const |
89 | { |
90 | |
91 | Standard_Integer i,j; |
92 | Standard_Real const *pdb=NULL; // the current pointer to WeightsDB |
93 | switch (myNivConstr) { |
94 | case 0: pdb = WeightsDB_C0; break; |
95 | case 1: pdb = WeightsDB_C1; break; |
96 | case 2: pdb = WeightsDB_C2; break; |
97 | } |
98 | Standard_Integer infdg = 2*(myNivConstr+1); |
99 | if (NbGaussPoints > NDEG8) pdb += (NDEG8 *(NDEG8 -infdg)/2); |
100 | if (NbGaussPoints > NDEG10) pdb += (NDEG10*(NDEG10-infdg)/2); |
101 | if (NbGaussPoints > NDEG15) pdb += (((NDEG15-1)/2)*(NDEG15-infdg)); |
102 | if (NbGaussPoints > NDEG20) pdb += (NDEG20*(NDEG20-infdg)/2); |
103 | if (NbGaussPoints > NDEG25) pdb += (((NDEG25-1)/2)*(NDEG25-infdg)); |
104 | if (NbGaussPoints > NDEG30) pdb += (NDEG30*(NDEG30-infdg)/2); |
105 | if (NbGaussPoints > NDEG40) pdb += (NDEG40*(NDEG40-infdg)/2); |
106 | if (NbGaussPoints > NDEG50) pdb += (NDEG50*(NDEG50-infdg)/2); |
107 | |
108 | // the copy of TabWeightsDB into TabWeights |
109 | for (j=0; j<=myDegree; j++) { |
110 | for (i=1; i<=NbGaussPoints/2; i++) { |
111 | TabWeights.SetValue(i,j,*pdb++); |
112 | } |
113 | } |
114 | |
115 | if (NbGaussPoints % 2 == 1) { |
116 | // NbGaussPoints is odd - the values addition for 0. |
117 | Standard_Real const *pdb0=NULL; // the current pointer to WeightsDB0 |
118 | switch (myNivConstr) { |
119 | case 0: pdb0 = WeightsDB0_C0; break; |
120 | case 1: pdb0 = WeightsDB0_C1; break; |
121 | case 2: pdb0 = WeightsDB0_C2; break; |
122 | } |
123 | |
124 | if (NbGaussPoints > NDEG15) pdb0 += ((NDEG15-1-infdg)/2 + 1); |
125 | if (NbGaussPoints > NDEG25) pdb0 += ((NDEG25-1-infdg)/2 + 1); |
126 | |
127 | // the copy of TabWeightsDB0 into TabWeights |
128 | for (j=0; j<=myDegree; j+=2) |
129 | TabWeights.SetValue(0,j,*pdb0++); |
130 | for (j=1; j<=myDegree; j+=2) |
131 | TabWeights.SetValue(0,j,0.); |
132 | } |
133 | else { |
134 | for (j=0; j<=myDegree; j++) { |
135 | TabWeights.SetValue(0,j,UNDEFINED); |
136 | } |
137 | } |
138 | } |
139 | |
140 | //======================================================================= |
141 | //function : MaxValue |
142 | //purpose : |
143 | //======================================================================= |
144 | |
145 | void PLib_JacobiPolynomial::MaxValue(TColStd_Array1OfReal& TabMax) const |
146 | { |
147 | Standard_Real const *pdb=NULL; // the pointer to MaxValues |
148 | switch (myNivConstr) { |
149 | case 0: pdb = MaxValuesDB_C0; break; |
150 | case 1: pdb = MaxValuesDB_C1; break; |
151 | case 2: pdb = MaxValuesDB_C2; break; |
152 | } |
153 | for (Standard_Integer i=TabMax.Lower(); i <= TabMax.Upper(); i++) { |
154 | TabMax.SetValue(i,*pdb++); |
155 | } |
156 | } |
157 | |
158 | //======================================================================= |
159 | //function : MaxError |
160 | //purpose : |
161 | //======================================================================= |
162 | |
163 | Standard_Real PLib_JacobiPolynomial::MaxError(const Standard_Integer Dimension, |
164 | Standard_Real& JacCoeff, |
165 | const Standard_Integer NewDegree) const |
166 | { |
167 | Standard_Integer i,idim,ibeg,icut; |
168 | |
169 | math_Vector MaxErrDim(1,Dimension,0.); |
170 | |
171 | TColStd_Array1OfReal TabMax(0, myDegree+1); |
172 | MaxValue(TabMax); |
173 | |
174 | ibeg = 2*(myNivConstr+1); |
175 | icut = Max (ibeg, NewDegree+1); |
176 | Standard_Real * JacArray = &JacCoeff; |
177 | for (idim=1; idim<=Dimension; idim++) { |
178 | for (i=icut; i<=myWorkDegree; i++) { |
179 | MaxErrDim(idim) += Abs(JacArray[i*Dimension+idim-1]) * TabMax(i-ibeg); |
180 | } |
181 | } |
182 | Standard_Real MaxErr = MaxErrDim.Norm(); |
183 | return (MaxErr); |
184 | } |
185 | |
186 | //======================================================================= |
187 | //function : ReduceDegree |
188 | //purpose : |
189 | //======================================================================= |
190 | |
191 | void PLib_JacobiPolynomial::ReduceDegree(const Standard_Integer Dimension, |
192 | const Standard_Integer MaxDegree, |
193 | const Standard_Real Tol, |
194 | Standard_Real& JacCoeff, |
195 | Standard_Integer& NewDegree, |
196 | Standard_Real& MaxError) const |
197 | { |
198 | Standard_Integer i,idim,icut, ia = 2*(myNivConstr+1)-1; |
199 | Standard_Real Bid,Eps1,Error; |
200 | |
201 | math_Vector MaxErrDim(1,Dimension,0.); |
202 | |
203 | NewDegree = ia; |
204 | MaxError = 0.; |
205 | Error = 0.; |
206 | icut=ia+1; |
207 | |
208 | TColStd_Array1OfReal TabMax(0, myDegree+1); |
209 | MaxValue(TabMax); |
210 | Standard_Real * JacArray = &JacCoeff; |
211 | for (i=myWorkDegree; i>=icut; i--) { |
212 | for (idim=1; idim<=Dimension; idim++) { |
213 | MaxErrDim(idim) += Abs(JacArray[i*Dimension+idim-1]) * TabMax(i-icut); |
214 | } |
215 | Error = MaxErrDim.Norm(); |
216 | if (Error > Tol && i <= MaxDegree) { |
217 | NewDegree = i; |
218 | break; |
219 | } |
220 | else |
221 | MaxError = Error; |
222 | } |
223 | if (NewDegree==ia) { |
224 | Eps1=0.000000001; |
225 | NewDegree = 0; |
226 | for (i=ia; i>=1; i--) { |
227 | Bid = 0.; |
228 | for (idim=1; idim<=Dimension; idim++) { |
229 | Bid += Abs(JacArray[i*Dimension+idim-1]); |
230 | } |
231 | if (Bid > Eps1) { |
232 | NewDegree = i; |
233 | break; |
234 | } |
235 | } |
236 | } |
237 | } |
238 | |
239 | //======================================================================= |
240 | //function : AverageError |
241 | //purpose : |
242 | //======================================================================= |
243 | |
244 | Standard_Real PLib_JacobiPolynomial::AverageError(const Standard_Integer Dimension, |
245 | Standard_Real& JacCoeff, |
246 | const Standard_Integer NewDegree) |
247 | const |
248 | { |
249 | Standard_Integer i,idim, icut = Max (2*(myNivConstr+1)+1, NewDegree+1); |
250 | Standard_Real BidJ, AverageErr = 0.; |
251 | Standard_Real * JacArray = &JacCoeff; |
252 | for (idim=1; idim<=Dimension; idim++) { |
253 | for (i=icut; i<=myDegree; i++) { |
254 | BidJ = JacArray[i*Dimension+idim-1]; |
255 | AverageErr += BidJ*BidJ; |
256 | } |
257 | } |
258 | AverageErr = sqrt(AverageErr/2); |
259 | return (AverageErr); |
260 | } |
261 | |
262 | //======================================================================= |
263 | //function :ToCoefficients |
264 | //purpose : |
265 | //======================================================================= |
266 | |
267 | void PLib_JacobiPolynomial::ToCoefficients(const Standard_Integer Dimension, |
268 | const Standard_Integer Degree, |
269 | const TColStd_Array1OfReal& JacCoeff, |
270 | TColStd_Array1OfReal& Coefficients) const |
271 | { |
272 | const Standard_Integer MAXM=31; |
273 | Standard_Integer i,iptt,j,idim, ii, jj; |
274 | Standard_Real const *pTr=NULL; // the pointer to TransMatrix |
275 | Standard_Real Bid; |
276 | Standard_Integer ibegJC=JacCoeff.Lower(), ibegC=Coefficients.Lower(); |
277 | |
278 | switch (myNivConstr) { |
279 | case 0: pTr = &TransMatrix_C0[0][0]; break; |
280 | case 1: pTr = &TransMatrix_C1[0][0]; break; |
281 | case 2: pTr = &TransMatrix_C2[0][0]; break; |
282 | } |
283 | // the conversation for even elements of JacCoeff |
284 | for (i=0; i<=Degree/2; i++) { |
285 | iptt = i*MAXM-(i+1)*i/2; |
286 | for (idim=1; idim<=Dimension; idim++) { |
287 | Bid = 0.; |
288 | for (j=i; j<=Degree/2; j++) { |
289 | Bid += (*(pTr+iptt+j)) * JacCoeff(2*j*Dimension+idim-1); |
290 | } |
291 | Coefficients.SetValue(2*i*Dimension+idim-1, Bid); |
292 | } |
293 | } |
294 | |
295 | if (Degree == 0) return; |
296 | |
297 | // the conversation for odd elements of JacCoeff |
298 | pTr += MAXM*(MAXM+1)/2; |
299 | for (i=0; i<=(Degree-1)/2; i++) { |
300 | iptt = i*MAXM-(i+1)*i/2; |
301 | ii = ibegC+(2*i+1)*Dimension; |
302 | for (idim=1; idim<=Dimension; idim++, ii++) { |
303 | Bid = 0.; |
304 | jj = ibegJC+(2*i+1)*Dimension+idim-1; |
305 | for (j=i; j<=(Degree-1)/2; j++, jj+=2*Dimension) { |
306 | Bid += (*(pTr+iptt+j)) * JacCoeff(jj); |
307 | } |
308 | Coefficients(ii) = Bid; |
309 | } |
310 | } |
311 | } |
312 | |
313 | //======================================================================= |
314 | //function : D0123 |
315 | //purpose : common part of D0,D1,D2,D3 (FORTRAN subroutine MPOJAC) |
316 | //======================================================================= |
317 | |
318 | void PLib_JacobiPolynomial::D0123(const Standard_Integer NDeriv, |
319 | const Standard_Real U, |
320 | TColStd_Array1OfReal& BasisValue, |
321 | TColStd_Array1OfReal& BasisD1, |
322 | TColStd_Array1OfReal& BasisD2, |
323 | TColStd_Array1OfReal& BasisD3) |
324 | { |
325 | Standard_Integer i,j, HermitNivConstr = 2*(myNivConstr+1); |
326 | Standard_Real Aux1,Aux2; |
327 | |
328 | if (myTNorm.IsNull()) { |
329 | |
330 | // Inizialization of myTNorm,myCofA,myCofB,myDenom |
331 | |
332 | myTNorm = new TColStd_HArray1OfReal(0,myDegree); |
333 | for (i=0; i<=myDegree; i++) { |
334 | Aux2 = 1.; |
335 | for (j=1; j<=HermitNivConstr; j++) { |
336 | Aux2 *= ((Standard_Real)(i+HermitNivConstr+j)/(Standard_Real)(i+j)); |
337 | } |
338 | myTNorm->SetValue(i, Sqrt (Aux2 * (2*i+2*HermitNivConstr+1) / |
339 | (Pow (2,2*HermitNivConstr+1)))); |
340 | } |
341 | |
342 | if(myDegree >= 2) { |
343 | myCofA = new TColStd_HArray1OfReal(0,myDegree); |
344 | myCofB = new TColStd_HArray1OfReal(0,myDegree); |
345 | myDenom = new TColStd_HArray1OfReal(0,myDegree); |
346 | for (i=2; i<=myDegree; i++) { |
347 | Aux1 = HermitNivConstr+i-1; |
348 | Aux2 = 2 * Aux1; |
349 | myCofA ->SetValue(i, Aux2*(Aux2+1)*(Aux2+2)); |
350 | myCofB ->SetValue(i, -2. *(Aux2+2) * Aux1* Aux1); |
351 | myDenom->SetValue(i, 1./(2. * i * ( i-1 + 2*HermitNivConstr+1) * Aux2)); |
352 | } |
353 | } |
354 | } |
355 | |
356 | // --- Positionements triviaux ----- |
357 | Standard_Integer ibeg0 = BasisValue.Lower(); |
358 | Standard_Integer ibeg1 = BasisD1.Lower(); |
359 | Standard_Integer ibeg2 = BasisD2.Lower(); |
360 | Standard_Integer ibeg3 = BasisD3.Lower(); |
361 | Standard_Integer i0, i1, i2, i3; |
362 | |
363 | |
364 | if (myDegree == 0) { |
365 | BasisValue(ibeg0+0) = 1.; |
366 | if (NDeriv >= 1) { |
367 | BasisD1(ibeg1+0) = 0.; |
368 | if (NDeriv >= 2) { |
369 | BasisD2(ibeg2+0) = 0.; |
370 | if (NDeriv == 3) |
371 | BasisD3(ibeg3+0) = 0.; |
372 | } |
373 | } |
374 | } |
375 | else { |
376 | BasisValue(ibeg0+0) = 1.; |
377 | Aux1 = HermitNivConstr+1; |
378 | BasisValue(ibeg0+1) = Aux1 * U; |
379 | if (NDeriv >= 1) { |
380 | BasisD1(ibeg1+0) = 0.; |
381 | BasisD1(ibeg1+1) = Aux1; |
382 | if (NDeriv >= 2) { |
383 | BasisD2(ibeg2+0) = 0.; |
384 | BasisD2(ibeg2+1) = 0.; |
385 | if (NDeriv == 3) { |
386 | BasisD3(ibeg3+0) = 0.; |
387 | BasisD3(ibeg3+1) = 0.; |
388 | } |
389 | } |
390 | } |
391 | } |
392 | |
393 | |
394 | // --- Positionement par reccurence |
395 | if (myDegree > 1) { |
396 | if (NDeriv == 0) { |
397 | Standard_Real * BV = &BasisValue(ibeg0); |
398 | Standard_Real * CofA = &myCofA->ChangeValue(0); |
399 | Standard_Real * CofB = &myCofB->ChangeValue(0); |
400 | Standard_Real * Denom = &myDenom->ChangeValue(0); |
401 | for (i=2; i<=myDegree; i++) { |
402 | BV[i] = (CofA[i]*U*BV[i-1] + CofB[i]*BV[i-2])*Denom[i]; |
403 | } |
404 | } |
405 | |
406 | else { |
407 | Standard_Real CofA, CofB, Denom; |
408 | for (i=2; i<=myDegree; i++) { |
409 | i0=i+ibeg0; |
410 | i1=i+ibeg1; |
411 | CofA = myCofA->Value(i); |
412 | CofB = myCofB->Value(i); |
413 | Denom = myDenom->Value(i); |
414 | |
415 | BasisValue(i0) = (CofA * U * BasisValue(i0-1) + |
416 | CofB * BasisValue(i0-2)) * Denom; |
417 | BasisD1(i1) = (CofA * (U * BasisD1(i1-1) + BasisValue(i0-1)) + |
418 | CofB * BasisD1(i1-2)) * Denom; |
419 | if (NDeriv >= 2) { |
420 | i2=i+ibeg2; |
421 | BasisD2(i2) = ( CofA * (U*BasisD2(i2-1) + 2*BasisD1(i1-1)) + |
422 | CofB*BasisD2(i2-2)) * Denom; |
423 | if (NDeriv == 3) { |
424 | i3=i+ibeg3; |
425 | BasisD3(i3) = (CofA * (U*BasisD3(i3-1) + 3*BasisD2(i2-1)) + |
426 | CofB*BasisD3(i3-2)) * Denom; |
427 | } |
428 | } |
429 | } |
430 | } |
431 | } |
432 | |
433 | // Normalization |
434 | if (NDeriv == 0) { |
435 | Standard_Real * BV = &BasisValue(ibeg0); |
436 | Standard_Real * TNorm = &myTNorm->ChangeValue(0); |
437 | for (i=0; i<=myDegree; i++) |
438 | BV[i] *= TNorm[i]; |
439 | } |
440 | else { |
441 | Standard_Real TNorm; |
442 | for (i=0; i<=myDegree; i++) { |
443 | TNorm = myTNorm->Value(i); |
444 | BasisValue(i+ibeg0) *= TNorm; |
445 | BasisD1(i+ibeg1) *= TNorm; |
446 | if (NDeriv >= 2) { |
447 | BasisD2(i+ibeg2) *= TNorm; |
448 | if (NDeriv >= 3) BasisD3(i+ibeg3) *= TNorm; |
449 | } |
450 | } |
451 | } |
452 | } |
453 | |
454 | //======================================================================= |
455 | //function : D0 |
456 | //purpose : |
457 | //======================================================================= |
458 | |
459 | void PLib_JacobiPolynomial::D0(const Standard_Real U, |
460 | TColStd_Array1OfReal& BasisValue) |
461 | { |
462 | D0123(0,U,BasisValue,BasisValue,BasisValue,BasisValue); |
463 | } |
464 | |
465 | //======================================================================= |
466 | //function : D1 |
467 | //purpose : |
468 | //======================================================================= |
469 | |
470 | void PLib_JacobiPolynomial::D1(const Standard_Real U, |
471 | TColStd_Array1OfReal& BasisValue, |
472 | TColStd_Array1OfReal& BasisD1) |
473 | { |
474 | D0123(1,U,BasisValue,BasisD1,BasisD1,BasisD1); |
475 | } |
476 | |
477 | //======================================================================= |
478 | //function : D2 |
479 | //purpose : |
480 | //======================================================================= |
481 | |
482 | void PLib_JacobiPolynomial::D2(const Standard_Real U, |
483 | TColStd_Array1OfReal& BasisValue, |
484 | TColStd_Array1OfReal& BasisD1, |
485 | TColStd_Array1OfReal& BasisD2) |
486 | { |
487 | D0123(2,U,BasisValue,BasisD1,BasisD2,BasisD2); |
488 | } |
489 | |
490 | //======================================================================= |
491 | //function : D3 |
492 | //purpose : |
493 | //======================================================================= |
494 | |
495 | void PLib_JacobiPolynomial::D3(const Standard_Real U, |
496 | TColStd_Array1OfReal& BasisValue, |
497 | TColStd_Array1OfReal& BasisD1, |
498 | TColStd_Array1OfReal& BasisD2, |
499 | TColStd_Array1OfReal& BasisD3) |
500 | { |
501 | D0123(3,U,BasisValue,BasisD1,BasisD2,BasisD3); |
502 | } |
503 | |