EdgePolynomialBase.cpp

Go to the documentation of this file.

/** \file EdgePolynomialBase.cpp
\brief Implementation of Ainsworth-Cole H1 base on edge
*/
 
using namespace MoFEM;
 
MoFEMErrorCode
EdgePolynomialBase::query_interface(boost::typeindex::type_index type_index,
                                    UnknownInterface **iface) const {
  *iface = const_cast<EdgePolynomialBase *>(this);
  return 0;
}
 
MoFEMErrorCode
EdgePolynomialBase::getValue(MatrixDouble &pts,
                             boost::shared_ptr<BaseFunctionCtx> ctx_ptr) {
  MoFEMFunctionBegin;
 
  cTx = ctx_ptr->getInterface<EntPolynomialBaseCtx>();
 
  int nb_gauss_pts = pts.size2();
  if (!nb_gauss_pts)
    MoFEMFunctionReturnHot(0);
 
  if (pts.size1() < 1)
    SETERRQ(
        PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
        "Wrong dimension of pts, should be at least 1 row with coordinates");
 
  const FieldApproximationBase base = cTx->bAse;
  EntitiesFieldData &data = cTx->dAta;
 
  if (base != AINSWORTH_BERNSTEIN_BEZIER_BASE) {
    if (cTx->copyNodeBase == LASTBASE) {
      data.dataOnEntities[MBVERTEX][0].getN(base).resize(nb_gauss_pts, 2,
                                                         false);
      CHKERR ShapeMBEDGE(
          &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
          &pts(0, 0), nb_gauss_pts);
    } else
      data.dataOnEntities[MBVERTEX][0].getN(base) =
          data.dataOnEntities[MBVERTEX][0].getN(cTx->copyNodeBase);
 
    if (data.dataOnEntities[MBVERTEX][0].getN(base).size1() !=
        (unsigned int)nb_gauss_pts)
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
               "Base functions or nodes has wrong number of integration points "
               "for base %s",
               ApproximationBaseNames[base]);
 
    data.dataOnEntities[MBVERTEX][0].getDiffN(base).resize(nb_gauss_pts, 2 * 1,
                                                           false);
    for (auto gg = 0; gg != nb_gauss_pts; ++gg)
      std::copy(Tools::diffShapeFunMBEDGE.begin(),
                Tools::diffShapeFunMBEDGE.end(),
                &data.dataOnEntities[MBVERTEX][0].getDiffN(base)(gg, 0));
  }
 
  switch (cTx->sPace) {
  case H1:
    MoFEMFunctionReturnHot(getValueH1(pts));
  case HDIV:
    MoFEMFunctionReturnHot(getValueHdiv(pts));
  case HCURL:
    MoFEMFunctionReturnHot(getValueHcurl(pts));
  case L2:
    MoFEMFunctionReturnHot(getValueL2(pts));
  default:
    SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not yet implemented");
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode EdgePolynomialBase::getValueH1(MatrixDouble &pts) {
  MoFEMFunctionBeginHot;
 
  switch (cTx->bAse) {
  case AINSWORTH_LEGENDRE_BASE:
  case AINSWORTH_LOBATTO_BASE:
    MoFEMFunctionReturnHot(getValueH1AinsworthBase(pts));
    break;
  case DEMKOWICZ_JACOBI_BASE:
    MoFEMFunctionReturnHot(getValueH1DemkowiczBase(pts));
    break;
  case AINSWORTH_BERNSTEIN_BEZIER_BASE:
    MoFEMFunctionReturnHot(getValueH1BernsteinBezierBase(pts));
    break;
  default:
    SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
  }
 
  MoFEMFunctionReturnHot(0);
}
 
MoFEMErrorCode
EdgePolynomialBase::getValueH1BernsteinBezierBase(MatrixDouble &pts) {
  MoFEMFunctionBegin;
 
  EntitiesFieldData &data = cTx->dAta;
  const std::string &field_name = cTx->fieldName;
  const int nb_gauss_pts = pts.size2();
 
  auto get_alpha = [field_name](auto &data) -> MatrixInt & {
    auto &ptr = data.getBBAlphaIndicesSharedPtr(field_name);
    if (!ptr)
      ptr.reset(new MatrixInt());
    return *ptr;
  };
 
  auto get_base = [field_name](auto &data) -> MatrixDouble & {
    auto &ptr = data.getBBNSharedPtr(field_name);
    if (!ptr)
      ptr.reset(new MatrixDouble());
    return *ptr;
  };
 
  auto get_diff_base = [field_name](auto &data) -> MatrixDouble & {
    auto &ptr = data.getBBDiffNSharedPtr(field_name);
    if (!ptr)
      ptr.reset(new MatrixDouble());
    return *ptr;
  };
 
  if (data.dataOnEntities[MBVERTEX].size() != 1 ||
      data.dataOnEntities[MBEDGE].size() != 1)
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Wrong size ent of ent data");
  auto &vert_ent_data = data.dataOnEntities[MBVERTEX][0];
  auto &edge_ent_data = data.dataOnEntities[MBEDGE][0];
 
  if (vert_ent_data.getN(NOBASE).size1() != (unsigned int)nb_gauss_pts)
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
             "Base functions or nodes has wrong number of integration points "
             "for base %s",
             ApproximationBaseNames[NOBASE]);
  auto &lambda = vert_ent_data.getN(NOBASE);
 
  // The BB base is not hierarchical, so mixing degrees between the vertices
  // and the edge loses completeness. H1 requires the edge order on the edge
  // (when it carries DOFs) and vertex degrees between one and that order.
  auto &node_order = vert_ent_data.getBBNodeOrder();
  if (node_order.size() != 2)
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Wrong number of vertex orders for BB base: %d",
            static_cast<int>(node_order.size()));
  const bool edge_has_h1 = data.spacesOnEntities[MBEDGE].test(H1);
  const int edge_order = edge_has_h1 ? edge_ent_data.getOrder() : 0;
  const int order = std::max({node_order[0], node_order[1], edge_order});
  if (order < 1)
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Bernstein-Bezier H1 base requires approximation order >= 1");
  if (static_cast<size_t>(order) >=
      EntitiesFieldData::EntData::MaxBernsteinBezierOrder)
    SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
            "Bernstein-Bezier order %d exceeds the supported maximum %d", order,
            static_cast<int>(
                EntitiesFieldData::EntData::MaxBernsteinBezierOrder) -
                1);
  for (int n = 0; n != 2; ++n)
    if (node_order[n] < 1)
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Bernstein-Bezier H1 base requires vertex order between 1 and "
              "the element order %d; vertex %d has order %d",
              order, n, node_order[n]);
  if (NBEDGE_H1(order) && edge_has_h1 && edge_order != order)
    SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
            "Bernstein-Bezier base in H1 space requires uniform approximation "
            "order on the element: edge has order %d, element order is %d. "
            "BB functions are not hierarchical, so mixed degrees do not span "
            "a complete polynomial space. Set one order for the whole H1 "
            "field (L2 fields are not restricted).",
            edge_order, order);
 
  auto &vert_get_n = get_base(vert_ent_data);
  auto &vert_get_diff_n = get_diff_base(vert_ent_data);
  vert_get_n.resize(nb_gauss_pts, 2, false);
  vert_get_diff_n.resize(nb_gauss_pts, 2, false);
  vert_get_n.clear();
  vert_get_diff_n.clear();
 
  auto &vertex_alpha = get_alpha(vert_ent_data);
  vertex_alpha.resize(2, 2, false);
  vertex_alpha(0, 0) = node_order[0];
  vertex_alpha(0, 1) = 0;
  vertex_alpha(1, 0) = 0;
  vertex_alpha(1, 1) = node_order[1];
 
  CHKERR BernsteinBezier::baseFunctionsEdge(
      1, nb_gauss_pts, vertex_alpha.size1(), &vertex_alpha(0, 0), &lambda(0, 0),
      Tools::diffShapeFunMBEDGE.data(), &vert_get_n(0, 0),
      &vert_get_diff_n(0, 0));
  std::array<double, 2> f = {boost::math::factorial<double>(node_order[0]),
                             boost::math::factorial<double>(node_order[1])};
 
  for (int g = 0; g != nb_gauss_pts; ++g)
    for (int n = 0; n != 2; ++n) {
      vert_get_n(g, n) *= f[n];
      vert_get_diff_n(g, n) *= f[n];
    }
 
  // Edge interior functions. Without H1 on the edge, or without edge DOFs,
  // leave empty matrices and no alpha indices for this field.
  const int nb_dofs = edge_has_h1 ? NBEDGE_H1(edge_order) : 0;
  auto &edge_get_n = get_base(edge_ent_data);
  auto &edge_get_diff_n = get_diff_base(edge_ent_data);
  edge_get_n.resize(nb_gauss_pts, nb_dofs, false);
  edge_get_diff_n.resize(nb_gauss_pts, nb_dofs, false);
 
  if (nb_dofs) {
    auto &edge_alpha = get_alpha(edge_ent_data);
    edge_alpha.resize(nb_dofs, 2, false);
    CHKERR BernsteinBezier::generateIndicesEdgeEdge(edge_order,
                                                    &edge_alpha(0, 0));
    CHKERR BernsteinBezier::baseFunctionsEdge(
        edge_order, nb_gauss_pts, edge_alpha.size1(), &edge_alpha(0, 0),
        &lambda(0, 0), Tools::diffShapeFunMBEDGE.data(), &edge_get_n(0, 0),
        &edge_get_diff_n(0, 0));
  } else {
    edge_ent_data.getBBAlphaIndicesSharedPtr(field_name).reset();
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode EdgePolynomialBase::getValueH1AinsworthBase(MatrixDouble &pts) {
  MoFEMFunctionBeginHot;
 
  EntitiesFieldData &data = cTx->dAta;
  const FieldApproximationBase base = cTx->bAse;
  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
                                     double *diffL, const int dim) =
      cTx->basePolynomialsType0;
 
  int nb_gauss_pts = pts.size2();
 
  const int side_number = 0;
  int sense = data.dataOnEntities[MBEDGE][side_number].getSense();
  int order = data.dataOnEntities[MBEDGE][side_number].getOrder();
  data.dataOnEntities[MBEDGE][side_number].getN(base).resize(
      nb_gauss_pts, NBEDGE_H1(order), false);
  data.dataOnEntities[MBEDGE][side_number].getDiffN(base).resize(
      nb_gauss_pts, NBEDGE_H1(order), false);
 
  data.dataOnEntities[MBEDGE][side_number].getN(base).clear();
  data.dataOnEntities[MBEDGE][side_number].getDiffN(base).clear();
 
  L.resize(NBEDGE_H1(order), false);
  diffL.resize(NBEDGE_H1(order), false);
 
  if (data.dataOnEntities[MBEDGE][side_number].getOrder() > 1) {
 
#ifndef NDEBUG
    if (sense != 1 && sense != -1)
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Wrong sense of edge, should be 1 or -1");
#endif
 
    for (int gg = 0; gg != nb_gauss_pts; gg++) {
 
      double s = (2. * pts(0, gg) - 1) * sense; // from [-1..1] * sense
      double diff_s = 2. * sense; // oriented s = s(xi), ds/dxi = 2 * sense
 
      // calculate Legendre polynomials at integration points
      CHKERR base_polynomials(NBEDGE_H1(order) - 1, s, &diff_s,
                              &*L.data().begin(), &*diffL.data().begin(), 1);
 
      for (unsigned int pp = 0;
           pp < data.dataOnEntities[MBEDGE][side_number].getN(base).size2();
           pp++) {
 
        // Calculate edge shape functions N0*N1*L(p), where N0 and N1 are nodal
        // shape functions
        double v = data.dataOnEntities[MBVERTEX][0].getN(base)(gg, 0) *
                   data.dataOnEntities[MBVERTEX][0].getN(base)(gg, 1);
        data.dataOnEntities[MBEDGE][side_number].getN(base)(gg, pp) = v * L(pp);
 
        // Calculate derivative edge shape functions
        // dN/dksi = dN0/dxi*N1*L + N0*dN1/ksi*L + N0*N1*dL/dxi
        data.dataOnEntities[MBEDGE][side_number].getDiffN(base)(gg, pp) =
            ((-1.) * data.dataOnEntities[MBVERTEX][0].getN(base)(gg, 1) +
             data.dataOnEntities[MBVERTEX][0].getN(base)(gg, 0) * (+1.)) *
                L(pp) +
            v * diffL(pp);
      }
    }
  }
 
  MoFEMFunctionReturnHot(0);
}
 
MoFEMErrorCode EdgePolynomialBase::getValueH1DemkowiczBase(MatrixDouble &pts) {
  MoFEMFunctionBeginHot;
 
  EntitiesFieldData &data = cTx->dAta;
  const FieldApproximationBase base = cTx->bAse;
 
  int nb_gauss_pts = pts.size2();
 
  const int side_number = 0;
  int order = data.dataOnEntities[MBEDGE][side_number].getOrder();
 
  data.dataOnEntities[MBEDGE][side_number].getN(base).resize(
      nb_gauss_pts, NBEDGE_H1(order), false);
  data.dataOnEntities[MBEDGE][side_number].getDiffN(base).resize(
      nb_gauss_pts, NBEDGE_H1(order), false);
 
  data.dataOnEntities[MBEDGE][side_number].getN(base).clear();
  data.dataOnEntities[MBEDGE][side_number].getDiffN(base).clear();
 
  double diff_shape[] = {-1, 1};
  MatrixDouble shape(nb_gauss_pts, 2);
 
  if (NBEDGE_L2(order)) {
    for (int gg = 0; gg != nb_gauss_pts; gg++) {
      const double mu0 = 1.0 - pts(0, gg); // pts ranges over [0, 1]
      const double mu1 = pts(0, gg);
      shape(gg, 0) = mu0;
      shape(gg, 1) = mu1;
    }
 
    CHKERR DemkowiczHexAndQuad::H1_BubbleShapeFunctions_ONSEGMENT(
        order, &*shape.data().begin(), diff_shape,
        &*data.dataOnEntities[MBEDGE][side_number].getN(base).data().begin(),
        &*data.dataOnEntities[MBEDGE][side_number]
              .getDiffN(base)
              .data()
              .begin(),
        nb_gauss_pts);
  }
 
  MoFEMFunctionReturnHot(0);
}
 
MoFEMErrorCode EdgePolynomialBase::getValueL2(MatrixDouble &pts) {
  MoFEMFunctionBeginHot;
 
  switch (cTx->bAse) {
  case AINSWORTH_LEGENDRE_BASE:
  case AINSWORTH_LOBATTO_BASE:
    MoFEMFunctionReturnHot(getValueL2AinsworthBase(pts));
  case DEMKOWICZ_JACOBI_BASE:
    MoFEMFunctionReturnHot(getValueL2DemkowiczBase(pts));
  default:
    SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
  }
 
  MoFEMFunctionReturnHot(0);
}
 
MoFEMErrorCode EdgePolynomialBase::getValueL2AinsworthBase(MatrixDouble &pts) {
  MoFEMFunctionBeginHot;
 
  EntitiesFieldData &data = cTx->dAta;
  FieldApproximationBase base = cTx->bAse;
 
  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
                                     double *diffL, const int dim) =
      cTx->basePolynomialsType0;
 
  if(base == AINSWORTH_LOBATTO_BASE)  
    base_polynomials = Lobatto_polynomials;
 
  int nb_gauss_pts = pts.size2();
 
  constexpr int side_number = 0;
  int order = data.dataOnEntities[MBEDGE][side_number].getOrder();
 
  data.dataOnEntities[MBEDGE][side_number].getN(base).resize(
      nb_gauss_pts, NBEDGE_L2(order), false);
  data.dataOnEntities[MBEDGE][side_number].getDiffN(base).resize(
      nb_gauss_pts, NBEDGE_L2(order), false);
 
  auto *fun_n =
      &*data.dataOnEntities[MBEDGE][side_number].getN(base).data().begin();
  auto *diff_fun_n =
      &*data.dataOnEntities[MBEDGE][side_number].getDiffN(base).data().begin();
 
  if (NBEDGE_L2(order)) {
 
    double diff_mu = 2;
    double l[NBEDGE_L2(order)];
    double diff_l[NBEDGE_L2(order)];
    for (int gg = 0; gg != nb_gauss_pts; gg++) {
      double mu = 2 * pts(0, gg) - 1;
      CHKERR base_polynomials(order, mu, &diff_mu, l, diff_l, 1);
      int qd_shift = NBEDGE_L2(order) * gg;
      for (int n = 0; n != NBEDGE_L2(order); n++) {
        fun_n[qd_shift + n] = l[n];
        diff_fun_n[qd_shift + n] = diff_l[n];
      }
    }
 
  }
 
  MoFEMFunctionReturnHot(0);
}
 
MoFEMErrorCode EdgePolynomialBase::getValueL2DemkowiczBase(MatrixDouble &pts) {
  return getValueL2AinsworthBase(pts);
}
 
MoFEMErrorCode EdgePolynomialBase::getValueHdiv(MatrixDouble &pts) {
  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
          "Make no sense, unless problem is 2d (2d not implemented yet)");
}
 
MoFEMErrorCode
EdgePolynomialBase::getValueHcurlAinsworthBase(MatrixDouble &pts) {
  MoFEMFunctionBegin;
 
  EntitiesFieldData &data = cTx->dAta;
  const FieldApproximationBase base = cTx->bAse;
 
  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
                                     double *diffL, const int dim) =
      cTx->basePolynomialsType0;
 
  int nb_gauss_pts = pts.size2();
  if (data.spacesOnEntities[MBEDGE].test(HCURL)) {
 
    if (data.dataOnEntities[MBEDGE].size() != 1)
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
 
    int sense = data.dataOnEntities[MBEDGE][0].getSense();
    int order = data.dataOnEntities[MBEDGE][0].getOrder();
    int nb_dofs =
        NBEDGE_AINSWORTH_HCURL(data.dataOnEntities[MBEDGE][0].getOrder());
    data.dataOnEntities[MBEDGE][0].getN(base).resize(nb_gauss_pts, 3 * nb_dofs,
                                                     false);
    data.dataOnEntities[MBEDGE][0].getDiffN(base).resize(nb_gauss_pts, 0,
                                                         false);
 
    CHKERR Hcurl_Ainsworth_EdgeBaseFunctions_MBTET_ON_EDGE(
        sense, order, &data.dataOnEntities[MBVERTEX][0].getN(base)(0, 0),
        &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
        &*data.dataOnEntities[MBEDGE][0].getN(base).data().begin(), nullptr,
        nb_gauss_pts, base_polynomials);
 
  } else {
    data.dataOnEntities[MBEDGE][0].getN(base).resize(nb_gauss_pts, 0, false);
    data.dataOnEntities[MBEDGE][0].getDiffN(base).resize(nb_gauss_pts, 0,
                                                         false);
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode
EdgePolynomialBase::getValueHcurlDemkowiczBase(MatrixDouble &pts) {
  MoFEMFunctionBegin;
 
  EntitiesFieldData &data = cTx->dAta;
  const FieldApproximationBase base = cTx->bAse;
 
  int nb_gauss_pts = pts.size2();
  if (data.spacesOnEntities[MBEDGE].test(HCURL)) {
 
    if (data.dataOnEntities[MBEDGE].size() != 1)
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "No data structure to store base functions");
 
    int sense = data.dataOnEntities[MBEDGE][0].getSense();
    int order = data.dataOnEntities[MBEDGE][0].getOrder();
    int nb_dofs =
        NBEDGE_DEMKOWICZ_HCURL(data.dataOnEntities[MBEDGE][0].getOrder());
    data.dataOnEntities[MBEDGE][0].getN(base).resize(nb_gauss_pts, 3 * nb_dofs,
                                                     false);
    data.dataOnEntities[MBEDGE][0].getDiffN(base).resize(nb_gauss_pts, 0,
                                                         false);
    if (nb_dofs) {
      CHKERR Hcurl_Demkowicz_EdgeBaseFunctions_MBEDGE(
          sense, order, &data.dataOnEntities[MBVERTEX][0].getN(base)(0, 0),
          &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
          &*data.dataOnEntities[MBEDGE][0].getN(base).data().begin(), NULL,
          nb_gauss_pts);
    }
  } else {
 
    data.dataOnEntities[MBEDGE][0].getN(base).resize(nb_gauss_pts, 0, false);
    data.dataOnEntities[MBEDGE][0].getDiffN(base).resize(nb_gauss_pts, 0,
                                                         false);
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode EdgePolynomialBase::getValueHcurl(MatrixDouble &pts) {
  MoFEMFunctionBeginHot;
 
  switch (cTx->bAse) {
  case AINSWORTH_LEGENDRE_BASE:
  case AINSWORTH_LOBATTO_BASE:
    MoFEMFunctionReturnHot(getValueHcurlAinsworthBase(pts));
  case DEMKOWICZ_JACOBI_BASE:
    MoFEMFunctionReturnHot(getValueHcurlDemkowiczBase(pts));
  default:
    SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
  }
 
  MoFEMFunctionReturnHot(0);
}