src/PLib/PLib_JacobiPolynomial.cxx

   1 // Copyright (c) 1997-1999 Matra Datavision
   2 // Copyright (c) 1999-2014 OPEN CASCADE SAS
   3 //
   4 // This file is part of Open CASCADE Technology software library.
   5 //
   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
   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.
  11 //
  12 // Alternatively, this file may be used under the terms of Open CASCADE
  13 // commercial license or contractual agreement.
  14
  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
  26 const Standard_Integer NDEG8=8,   NDEG10=10, NDEG15=15, NDEG20=20, NDEG25=25,
  27                        NDEG30=30, NDEG40=40, NDEG50=50, NDEG61=61;
  28
  29 const 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
  58 void PLib_JacobiPolynomial::Points(const Standard_Integer NbGaussPoints,
  59                                    TColStd_Array1OfReal& TabPoints) const
  60 {
  61   if ((NbGaussPoints != NDEG8  && NbGaussPoints != NDEG10 &&
  62       NbGaussPoints != NDEG15 && NbGaussPoints != NDEG20 &&
  63       NbGaussPoints != NDEG25 && NbGaussPoints != NDEG30 &&
  64       NbGaussPoints != NDEG40 && NbGaussPoints != NDEG50 &&
  65       NbGaussPoints != NDEG61) ||
  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
  87 void 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
 145 void 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
 163 Standard_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
 191 void 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
 244 Standard_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
 267 void 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
 318 void 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
 459 void 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
 470 void 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
 482 void 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
 495 void 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