0028417: Using PRECOMPILED HEADER to speed up compilation time
[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
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 24IMPLEMENT_STANDARD_RTTIEXT(PLib_JacobiPolynomial,PLib_Base)
25
896faa72 26#include "PLib_JacobiPolynomial_Data.pxx"
27
7fd59977 28// The possible values for NbGaussPoints
29const Standard_Integer NDEG8=8, NDEG10=10, NDEG15=15, NDEG20=20, NDEG25=25,
30 NDEG30=30, NDEG40=40, NDEG50=50, NDEG61=61;
31
32const 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
61void 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
90void 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
148void 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
166Standard_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
194void 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
247Standard_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
270void 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
321void 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
462void 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
473void 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
485void 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
498void 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