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 <PLib_JacobiPolynomial_0.hxx>
  21 #include <Standard_ConstructionError.hxx>
  22 #include <Standard_Type.hxx>
  23 #include <TColStd_Array2OfReal.hxx>
  24
  25 IMPLEMENT_STANDARD_RTTIEXT(PLib_JacobiPolynomial,PLib_Base)
  26
  27 // The possible values for NbGaussPoints
  28 const Standard_Integer NDEG8=8,   NDEG10=10, NDEG15=15, NDEG20=20, NDEG25=25,
  29                        NDEG30=30, NDEG40=40, NDEG50=50, NDEG61=61;
  30
  31 const 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:
  48       Standard_ConstructionError::Raise("Invalid ConstraintOrder");
  49   }
  50   myDegree = myWorkDegree - 2*(myNivConstr+1);
  51   if (myDegree > 30)
  52     Standard_ConstructionError::Raise("Invalid Degree");
  53 }
  54
  55 //=======================================================================
  56 //function : Points
  57 //purpose  :
  58 //=======================================================================
  59
  60 void PLib_JacobiPolynomial::Points(const Standard_Integer NbGaussPoints,
  61                                    TColStd_Array1OfReal& TabPoints) const
  62 {
  63   if ((NbGaussPoints != NDEG8  && NbGaussPoints != NDEG10 &&
  64       NbGaussPoints != NDEG15 && NbGaussPoints != NDEG20 &&
  65       NbGaussPoints != NDEG25 && NbGaussPoints != NDEG30 &&
  66       NbGaussPoints != NDEG40 && NbGaussPoints != NDEG50 &&
  67       NbGaussPoints != NDEG61) ||
  68       NbGaussPoints <= myDegree)
  69     Standard_ConstructionError::Raise("Invalid NbGaussPoints");
  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
  89 void 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
 147 void 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
 165 Standard_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
 193 void 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
 246 Standard_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
 269 void 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
 320 void 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
 461 void 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
 472 void 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
 484 void 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
 497 void 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