LobattoPolynomial.cpp

Go to the documentation of this file.

/** \file LobattoPolynomial.cpp
 * \brief implementation of multi-grid solver for p- adaptivity
 */
 
 
 
namespace MoFEM {
 
MoFEMErrorCode
LobattoPolynomialCtx::query_interface(boost::typeindex::type_index type_index,
                                      UnknownInterface **iface) const {
  *iface = const_cast<LobattoPolynomialCtx *>(this);
  return 0;
}
 
MoFEMErrorCode
LobattoPolynomial::query_interface(boost::typeindex::type_index type_index,
                                   UnknownInterface **iface) const {
  *iface = const_cast<LobattoPolynomial *>(this);
  return 0;
}
 
MoFEMErrorCode
LobattoPolynomial::getValue(MatrixDouble &pts,
                            boost::shared_ptr<BaseFunctionCtx> ctx_ptr) {
 
  MoFEMFunctionBeginHot;
  auto ctx = ctx_ptr->getInterface<LobattoPolynomialCtx>();
  // Polynomial order start from 2nd order
  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);
}
 
MoFEMErrorCode KernelLobattoPolynomialCtx::query_interface(
    boost::typeindex::type_index type_index, UnknownInterface **iface) const {
  *iface = const_cast<KernelLobattoPolynomialCtx *>(this);
  return 0;
}
 
MoFEMErrorCode KernelLobattoPolynomial::query_interface(
    boost::typeindex::type_index type_index, UnknownInterface **iface) const {
  *iface = const_cast<KernelLobattoPolynomial *>(this);
  return 0;
}
 
MoFEMErrorCode
KernelLobattoPolynomial::getValue(MatrixDouble &pts,
                                  boost::shared_ptr<BaseFunctionCtx> ctx_ptr) {
 
  MoFEMFunctionBeginHot;
  auto ctx = ctx_ptr->getInterface<KernelLobattoPolynomialCtx>();
  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);
}
 
PetscErrorCode Lobatto_polynomials(int p, double s, double *diff_s, double *L,
                                   double *diffL, const 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 < 2)
    SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "p < 2");
#endif // NDEBUG
  double l[p + 1];
  ierr = Legendre_polynomials(p, s, NULL, l, NULL, 1);
  CHKERRQ(ierr);
 
  L[0] = 1;
  if (diffL != NULL) {
    for (int d = 0; d != dim; ++d) {
      diffL[d * (p + 1) + 0] = 0;
    }
  }
  L[1] = s;
  if (diffL != NULL) {
#ifndef NDEBUG
    if (diff_s == NULL) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "diff_s == NULL");
    }
#endif // NDEBUG
    for (int d = 0; d != dim; ++d) {
      diffL[d * (p + 1) + 1] = diff_s[d];
    }
  }
 
  // Integrated Legendre
  for (int k = 2; k <= p; k++) {
    const double factor = 2 * (2 * k - 1);
    L[k] = 1.0 / factor * (l[k] - l[k - 2]);
  }
 
  if (diffL != NULL) {
    for (int k = 2; k <= p; k++) {
      double a = l[k - 1] / 2.;
      for (int d = 0; d != dim; ++d) {
        diffL[d * (p + 1) + k] = a * diff_s[d];
      }
    }
  }
  MoFEMFunctionReturnHot(0);
}
 
} // namespace MoFEM