projects / occt.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
0024285: Updates of PLib::EvalPolynomial for code acceleration
[occt.git] / src / PLib / PLib_JacobiPolynomial.cxx
CommitLineData
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#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
26const Standard_Integer NDEG8=8, NDEG10=10, NDEG15=15, NDEG20=20, NDEG25=25,
27 NDEG30=30, NDEG40=40, NDEG50=50, NDEG61=61;
28
29const 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
58void 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
87void 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
145void 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
163Standard_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
191void 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
244Standard_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
267void 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
318void 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
459void 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
470void 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
482void 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
495void 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
Open CASCADE Technology repository