projects / occt.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
0026937: Eliminate NO_CXX_EXCEPTION macro support
[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>
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 25IMPLEMENT_STANDARD_RTTIEXT(PLib_JacobiPolynomial,PLib_Base)
26
7fd59977 27// The possible values for NbGaussPoints
28const Standard_Integer NDEG8=8, NDEG10=10, NDEG15=15, NDEG20=20, NDEG25=25,
29 NDEG30=30, NDEG40=40, NDEG50=50, NDEG61=61;
30
31const 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
60void 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
89void 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
147void 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
165Standard_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
193void 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
246Standard_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
269void 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
320void 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
461void 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
472void 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
484void 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
497void 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
Open CASCADE Technology repository