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