LegendrePolynomial.cpp

Go to the documentation of this file.

/** \file LegendrePolynomial.cpp
 * \brief Implementation of base for Legendre polynomials
 */
 
namespace MoFEM {
 
MoFEMErrorCode
LegendrePolynomialCtx::query_interface(boost::typeindex::type_index type_index,
                                       UnknownInterface **iface) const {
  *iface = const_cast<LegendrePolynomialCtx *>(this);
  return 0;
}
 
MoFEMErrorCode
LegendrePolynomial::query_interface(boost::typeindex::type_index type_index,
                                    UnknownInterface **iface) const {
  *iface = const_cast<LegendrePolynomial *>(this);
  return 9;
}
 
MoFEMErrorCode
LegendrePolynomial::getValue(MatrixDouble &pts,
                             boost::shared_ptr<BaseFunctionCtx> ctx_ptr) {
 
  MoFEMFunctionBeginHot;
  auto ctx = ctx_ptr->getInterface<LegendrePolynomialCtx>();
  ctx->baseFunPtr->resize(pts.size2(), ctx->P + 1, false);
  ctx->baseDiffFunPtr->resize(pts.size2(), ctx->dIm * (ctx->P + 1), false);
  double *l = NULL;
  double *diff_l = NULL;
  for (unsigned int gg = 0; gg < pts.size2(); gg++) {
    if (ctx->baseFunPtr)
      l = &((*ctx->baseFunPtr)(gg, 0));
    if (ctx->baseDiffFunPtr)
      diff_l = &((*ctx->baseDiffFunPtr)(gg, 0));
    ierr = (ctx->basePolynomialsType0)(ctx->P, pts(0, gg), ctx->diffS, l,
                                       diff_l, ctx->dIm);
    CHKERRG(ierr);
  }
  MoFEMFunctionReturnHot(0);
}
 
/**
\brief Calculate Legendre approximation basis
 
\ingroup mofem_base_functions
 
Lagrange polynomial is given by
\f[
L_0(s)=1;\quad L_1(s) = s
\f]
and following terms are generated inductively
\f[
L_{l+1}=\frac{2l+1}{l+1}sL_l(s)-\frac{l}{l+1}L_{l-1}(s)
\f]
 
Note that:
\f[
s\in[-1,1] \quad \textrm{and}\; s=s(\xi_0,\xi_1,\xi_2)
\f]
where \f$\xi_i\f$ are barycentric coordinates of element.
 
\param p is approximation order
\param s is position \f$s\in[-1,1]\f$
\param diff_s derivatives of shape functions, i.e. \f$\frac{\partial s}{\partial
\xi_i}\f$ \retval L approximation functions \retval diffL derivatives, i.e.
\f$\frac{\partial L}{\partial \xi_i}\f$ \param dim dimension \return error code
 
*/
template <typename Scalar, typename DsScalar>
inline MoFEMErrorCode Legendre_polynomials_T(
    int p,
    Scalar &s, // value of the coordinate
    DsScalar *
        diff_s, // gradient of s wrt physical coords; may be nullptr iff diffL==nullptr
    Scalar *L,     // size p+1
    Scalar *diffL, // size dim*(p+1); may be nullptr
    int dim) {
  MoFEMFunctionBeginHot;
 
#ifndef NDEBUG
  if (dim < 1)
    SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "dim < 1");
  if (dim > 3)
    SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "dim > 3");
  if (p < 0)
    SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "p < 0");
  if (diffL != nullptr && diff_s == nullptr)
    SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
            "diff_s == NULL while diffL requested");
#endif
 
  // P0
  L[0] = Scalar(1);
  if (diffL) {
    for (int d = 0; d < dim; ++d)
      diffL[d * (p + 1) + 0] = Scalar(0);
  }
  if (p == 0)
    MoFEMFunctionReturnHot(0);
 
  // P1
  L[1] = s;
  if (diffL) {
    for (int d = 0; d < dim; ++d)
      diffL[d * (p + 1) + 1] = static_cast<Scalar>(diff_s[d]);
  }
  if (p == 1)
    MoFEMFunctionReturnHot(0);
 
  // Recurrence: P_{l+1} = A_l * s * P_l - B_l * P_{l-1}
  for (int l = 1; l < p; ++l) {
    const Scalar A = Scalar((2.0 * l + 1.0) / (l + 1.0));
    const Scalar B = Scalar(double(l) / (l + 1.0));
 
    L[l + 1] = A * s * L[l] - B * L[l - 1];
 
    if (diffL) {
      for (int d = 0; d < dim; ++d) {
        const Scalar ds = static_cast<Scalar>(diff_s[d]);
        // Chain rule: dP_{l+1} = A*( ds*P_l + s*dP_l ) - B*dP_{l-1}
        diffL[d * (p + 1) + (l + 1)] =
            A * (ds * L[l] + s * diffL[d * (p + 1) + l]) -
            B * diffL[d * (p + 1) + (l - 1)];
      }
    }
  }
 
  MoFEMFunctionReturnHot(0);
};
 
PetscErrorCode Legendre_polynomials(int p, double s, double *diff_s, double *L,
                                    double *diffL, int dim) {
  return Legendre_polynomials_T(p, s, diff_s, L, diffL, dim);
}
 
PetscErrorCode Legendre_polynomials(int p, std::complex<double> s,
                                    std::complex<double> *diff_s,
                                    std::complex<double> *L,
                                    std::complex<double> *diffL, int dim) {
  return Legendre_polynomials_T<std::complex<double>, std::complex<double>>(
      p, s, diff_s, L, diffL, dim);
}
 
} // namespace MoFEM