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