karniadakisOp.hh

Go to the documentation of this file.

 
#ifndef karniadakisOp_hh
#define karniadakisOp_hh
 
 
#include "toolbox/array.hh"
 
 
namespace concepts {
 
 inline const concepts::Array<Real> computeKarniadakisValues(uint np, const Real* absc,
    uint npx, uint type) {
 
  //at least first order functions
  conceptsAssert(np > 0, concepts::Assertion());
 
  //only value and derivative - evaluation implemented
  if (type > 1)
    throw conceptsException(concepts::MissingFeature("Only Karniadakis type evaluation not supported yet."));
 
  //values of derivative or values
  concepts::Array<Real> val((np + 1) * npx);
 
 
  //jacobi(alpha=1,beta=1) polynomials evaluations on [0,1] representations
  concepts::Array<Real> jacobi;
  if (np >= 2) {
    jacobi = concepts::Array<Real>((np - 1) * npx);
 
    //values for recursion formula
    Real an, bn, dn;
 
    for (uint i = 0; i < npx; ++i) {
      Real x = absc[i];
      // Q_0
      jacobi[i] = 1.0;
      if (np >= 3) {
        //Q_1
        jacobi[npx + i] = 4 * x - 2.0;
        //recursive formula
        //dn*Q_{n}(x) = (2an * x - an)Q_{n-1}(x) - bn * Q_{n-2}(x)
        for (uint n = 2; n < np - 1; ++n) {
          an = 2 * n * (2 * n + 1) * (2 * n + 2);
          bn = 2 * n * n * (2 * n + 2);
          dn = 2 * n * 2 * n * (n + 2);
          jacobi[n * npx + i] = ((2 * an * x - an) * jacobi[(n - 1)
              * npx + i] - bn * jacobi[(n - 2) * npx + i]) / dn;
        }// for j
      }// if np >= 3
    }//for i
  }// if np >= 2
 
 
  switch (type) {
  case (0): {
 
    //evaluate transformed karniadakis polynomials K_0 and K_1
    for (uint i = 0; i < npx; ++i) {
      val[i] = 1 - absc[i];
      val[npx + i] = absc[i];
    }
    //build karniadakis K_m evaluations
    //loop over jacobi polynomials
    for (uint n = 0; n < np - 1; ++n) {
      uint m = n + 2;
      //loop over abscissa points
      for (uint i = 0; i < npx; ++i) {
        // K_m =  (1-x)*x * Q_{m-2} = K_0*K_1 * Q_{m-2}, m = n + 2
        val[m * npx + i] = (val[i] * val[npx + i])
            * jacobi[n * npx + i];
      }
    }
    break;
  }//case normal values
 
  case (1): {
    //derivative of the polynomial = JacobiPolynomial ( alpha = 2, beta = 2)
    concepts::Array<Real> jacobiD;
    if (np >= 3) {
      jacobiD = concepts::Array<Real>((np - 2) * npx);
      //values for recursion formula
      Real an, bn, dn;
      for (uint i = 0; i < npx; ++i) {
        Real x = absc[i];
        // Q_0
        jacobiD[i] = 1.0;
        if (np >= 4) {
          //Q_1
          jacobiD[npx + i] = 6.0 * x - 3.0;
          //recursive formula
          //dn*Q_{n}(x) = (2an * x - an)Q_{n-1}(x) - bn * Q_{n-2}(x)
          for (uint n = 2; n < np - 2; ++n) {
            an = (2 * n + 2) * (2 * n + 3) * (2 * n + 4);
            bn = 2 * (n + 1) * (n + 1) * (2 * n + 4);
            dn = 2 * n * (n + 4) * (2 * n + 2);
            jacobiD[n * npx + i] = ((2 * an * x - an) * jacobiD[(n
                - 1) * npx + i] - bn * jacobiD[(n - 2) * npx
                + i]) / dn;
          }// for j
        }// if np >= 3
      }//for i
    }
    //evaluate karniadakis polynomials K_0 and K_1
    for (uint i = 0; i < npx; ++i) {
      val[i] = -0.5;
      val[npx + i] = 0.5;
    }
 
    if (np >= 2) {
      // (-0.25z^2+0.25)' ->  -0.5z   -> z = 2x-1 -> -x+0.5 (derivative of p=2 karniadakis)
 
      // loop over all fcts.
      for (uint n = 0; n < np - 1; ++n) {
        uint m = n + 2;
        // loop over all points
        for (uint i = 0; i < npx; ++i) {
          Real x = absc[i];
          // dx [ - 0.25z^2 +0.25 ] * dx/dz , z = 2x-1, as substituted   multiplied with P_n^(1,1)(2x-1)
          Real& v = val[m * npx + i] = (0.5 - x)
              * jacobi[n * npx + i];
 
          if (n > 0)
            // +  (- 0.25z^2 +0.25) *  d/dx P_n(1,1)(z) * dx/dz, z = 2x-1
            // =  (- x^2+x) *  (n+3)/2*P_{n-1}^(2,2)(2x-1) * 1/2
            v += (x - x * x) * (n + 3) / 2 * jacobiD[(n - 1) * npx
                + i];
        } //loop i
      } //loop n
    } //np>=2
    break;
  } //case == derivative
 
  default:
    //already thrown above
    throw conceptsException(concepts::MissingFeature("Karniadakis type evaluation not supported yet."));
    break;
  }
 
  return val;
}
 
 
 
} // namespace concepts
 
#endif // karniadakisOp_hh