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