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
  16 #include <math.hxx>
  17 #include <math_Vector.hxx>
  18 #include <PLib.hxx>
  19 #include <PLib_JacobiPolynomial.hxx>
  20 #include <Standard_ConstructionError.hxx>
  21 #include <Standard_Type.hxx>
  22 #include <TColStd_Array2OfReal.hxx>
  23
  24 IMPLEMENT_STANDARD_RTTIEXT(PLib_JacobiPolynomial,PLib_Base)
  25
  26 #include "PLib_JacobiPolynomial_Data.pxx"
  27
  28 // The possible values for NbGaussPoints
  29 const Standard_Integer NDEG8=8,   NDEG10=10, NDEG15=15, NDEG20=20, NDEG25=25,
  30                        NDEG30=30, NDEG40=40, NDEG50=50, NDEG61=61;
  31
  32 const 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:
  49       throw Standard_ConstructionError("Invalid ConstraintOrder");
  50   }
  51   myDegree = myWorkDegree - 2*(myNivConstr+1);
  52   if (myDegree > 30)
  53     throw Standard_ConstructionError("Invalid Degree");
  54 }
  55
  56 //=======================================================================
  57 //function : Points
  58 //purpose  :
  59 //=======================================================================
  60
  61 void PLib_JacobiPolynomial::Points(const Standard_Integer NbGaussPoints,
  62                                    TColStd_Array1OfReal& TabPoints) const
  63 {
  64   if ((NbGaussPoints != NDEG8  && NbGaussPoints != NDEG10 &&
  65       NbGaussPoints != NDEG15 && NbGaussPoints != NDEG20 &&
  66       NbGaussPoints != NDEG25 && NbGaussPoints != NDEG30 &&
  67       NbGaussPoints != NDEG40 && NbGaussPoints != NDEG50 &&
  68       NbGaussPoints != NDEG61) ||
  69       NbGaussPoints <= myDegree)
  70     throw Standard_ConstructionError("Invalid NbGaussPoints");
  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
  90 void 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
 148 void 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
 166 Standard_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
 194 void 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
 247 Standard_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
 270 void 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
 321 void 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
 462 void 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
 473 void 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
 485 void 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
 498 void 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