CGGUserPolynomialBase.cpp

Go to the documentation of this file.

/**
 * @file CGGUserPolynomialBase.cpp
 * @brief Implementation of CGGUserPolynomialBase class
 *
 */
 
CGGUserPolynomialBase::CGGUserPolynomialBase(
    boost::shared_ptr<CachePhi> cache_phi_otr)
    : TetPolynomialBase(), cachePhiPtr(cache_phi_otr) {}
 
MoFEMErrorCode CGGUserPolynomialBase::query_interface(
    boost::typeindex::type_index type_index,
    BaseFunctionUnknownInterface **iface) const {
  *iface = const_cast<CGGUserPolynomialBase *>(this);
  return 0;
}
 
MoFEMErrorCode
CGGUserPolynomialBase::getValue(MatrixDouble &pts,
                                boost::shared_ptr<BaseFunctionCtx> ctx_ptr) {
  MoFEMFunctionBeginHot;
 
  this->cTx = ctx_ptr->getInterface<EntPolynomialBaseCtx>();
 
  int nb_gauss_pts = pts.size2();
  if (!nb_gauss_pts) {
    MoFEMFunctionReturnHot(0);
  }
 
  if (pts.size1() < 3) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Wrong dimension of pts, should be at least 3 rows with "
            "coordinates");
  }
 
  const auto base = this->cTx->bAse;
  EntitiesFieldData &data = this->cTx->dAta;
 
  switch (this->cTx->sPace) {
  case HDIV:
    CHKERR getValueHdivForCGGBubble(pts);
    break;
  case L2:
    data.dataOnEntities[MBVERTEX][0].getN(base).resize(nb_gauss_pts, 4, false);
    CHKERR Tools::shapeFunMBTET(
        &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
        &pts(0, 0), &pts(1, 0), &pts(2, 0), nb_gauss_pts);
    data.dataOnEntities[MBVERTEX][0].getDiffN(base).resize(4, 3, false);
    std::copy(Tools::diffShapeFunMBTET.begin(), Tools::diffShapeFunMBTET.end(),
              data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin());
    this->cTx->basePolynomialsType0 = Legendre_polynomials;
    CHKERR getValueL2AinsworthBase(pts);
    break;
  default:
    SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not yet implemented");
  }
 
  MoFEMFunctionReturnHot(0);
}
 
MoFEMErrorCode
CGGUserPolynomialBase::getValueHdivForCGGBubble(MatrixDouble &pts) {
  MoFEMFunctionBegin;
 
  constexpr bool calculate_diff = false;
  constexpr bool debug = false;
 
  // This should be used only in case USER_BASE is selected
  if (this->cTx->bAse != USER_BASE) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Wrong base, should be USER_BASE");
  }
  // get access to data structures on element
  EntitiesFieldData &data = this->cTx->dAta;
  // Get approximation order on element. Note that bubble functions are only
  // on tetrahedron.
  const int order = data.dataOnEntities[MBTET][0].getOrder();
  /// number of integration points
  const int nb_gauss_pts = pts.size2();
 
  auto check_cache = [this](int order, int nb_gauss_pts) -> bool {
    if (cachePhiPtr) {
      return cachePhiPtr->get<0>() == order &&
             cachePhiPtr->get<1>() == nb_gauss_pts;
    } else {
      return false;
    }
  };
 
  if (check_cache(order, nb_gauss_pts)) {
    auto &mat = cachePhiPtr->get<2>();
    auto &phi = data.dataOnEntities[MBTET][0].getN(USER_BASE);
    phi.resize(mat.size1(), mat.size2(), false);
    noalias(phi) = mat;
    if constexpr (calculate_diff) {
      auto &diff_phi = data.dataOnEntities[MBTET][0].getDiffN(USER_BASE);
      auto &diff_mat = cachePhiPtr->get<3>();
      diff_phi.resize(diff_mat.size1(), diff_mat.size2(), false);
      noalias(diff_phi) = diff_mat;
    }
  } else {
    // calculate shape functions, i.e. barycentric coordinates
    shapeFun.resize(nb_gauss_pts, 4, false);
    CHKERR ShapeMBTET(&*shapeFun.data().begin(), &pts(0, 0), &pts(1, 0),
                      &pts(2, 0), nb_gauss_pts);
    // derivatives of shape functions
    double diff_shape_fun[12];
    CHKERR ShapeDiffMBTET(diff_shape_fun);
 
    const int nb_base_functions = NBVOLUMETET_CCG_BUBBLE(order);
    // get base functions and set size
    MatrixDouble &phi = data.dataOnEntities[MBTET][0].getN(USER_BASE);
    phi.resize(nb_gauss_pts, 9 * nb_base_functions, false);
    // finally calculate base functions
    FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3> t_phi(
        &phi(0, 0), &phi(0, 1), &phi(0, 2),
 
        &phi(0, 3), &phi(0, 4), &phi(0, 5),
 
        &phi(0, 6), &phi(0, 7), &phi(0, 8));
    CHKERR CGG_BubbleBase_MBTET(order, &shapeFun(0, 0), diff_shape_fun, t_phi,
                                nb_gauss_pts);
 
    auto imagingray_fd_directive = [&](auto &diff_phi) {
      MoFEMFunctionBegin;
 
      MatrixComplexDouble temp_phi;
      temp_phi.resize(nb_gauss_pts, 9 * nb_base_functions, false);
 
      MatrixComplexDouble diff_shape_fun;
      diff_shape_fun.resize(4, 3, false);
      diff_shape_fun(0, 0) = diffN_MBTET0x;
      diff_shape_fun(0, 1) = diffN_MBTET0y;
      diff_shape_fun(0, 2) = diffN_MBTET0z;
      diff_shape_fun(1, 0) = diffN_MBTET1x;
      diff_shape_fun(1, 1) = diffN_MBTET1y;
      diff_shape_fun(1, 2) = diffN_MBTET1z;
      diff_shape_fun(2, 0) = diffN_MBTET2x;
      diff_shape_fun(2, 1) = diffN_MBTET2y;
      diff_shape_fun(2, 2) = diffN_MBTET2z;
      diff_shape_fun(3, 0) = diffN_MBTET3x;
      diff_shape_fun(3, 1) = diffN_MBTET3y;
      diff_shape_fun(3, 2) = diffN_MBTET3z;
 
      MatrixComplexDouble shape_fun;
      shape_fun.resize(nb_gauss_pts, 4, false);
      for (int dd = 0; dd < 3; ++dd) {
        constexpr double eps = 1e-16;
 
        for (auto gg = 0; gg < nb_gauss_pts; ++gg) {
          std::array<std::complex<double>, 3> ksi{pts(0, gg), pts(1, gg),
                                                  pts(2, gg)};
 
          ksi[dd] += std::complex<double>(0, eps);
 
          shape_fun(gg, 0) = N_MBTET0(ksi[0], ksi[1], ksi[2]);
          shape_fun(gg, 1) = N_MBTET1(ksi[0], ksi[1], ksi[2]);
          shape_fun(gg, 2) = N_MBTET2(ksi[0], ksi[1], ksi[2]);
          shape_fun(gg, 3) = N_MBTET3(ksi[0], ksi[1], ksi[2]);
        }
 
        auto t_complex_phi_in =
            FTensor::Tensor2<FTensor::PackPtr<std::complex<double> *, 9>, 3, 3>{
                &temp_phi(0, 0), &temp_phi(0, 1), &temp_phi(0, 2),
                &temp_phi(0, 3), &temp_phi(0, 4), &temp_phi(0, 5),
                &temp_phi(0, 6), &temp_phi(0, 7), &temp_phi(0, 8)};
 
        CHKERR CGG_BubbleBase_MBTET(order, shape_fun.data().data(),
                                    diff_shape_fun.data().data(),
                                    t_complex_phi_in, nb_gauss_pts);
 
        auto t_complex_phi_out =
            FTensor::Tensor2<FTensor::PackPtr<std::complex<double> *, 9>, 3, 3>{
                &temp_phi(0, 0), &temp_phi(0, 1), &temp_phi(0, 2),
                &temp_phi(0, 3), &temp_phi(0, 4), &temp_phi(0, 5),
                &temp_phi(0, 6), &temp_phi(0, 7), &temp_phi(0, 8)};
        auto t_diff_phi = getFTensor3FromPtr<3, 3, 3>(diff_phi.data().data());
        for (auto gg = 0; gg < nb_gauss_pts; ++gg) {
          for (auto bb = 0; bb << nb_base_functions; ++bb) {
            t_diff_phi(0, 0, dd) = std::imag(t_complex_phi_out(0, 0)) / eps;
            t_diff_phi(0, 1, dd) = std::imag(t_complex_phi_out(0, 1)) / eps;
            t_diff_phi(0, 2, dd) = std::imag(t_complex_phi_out(0, 2)) / eps;
            t_diff_phi(1, 0, dd) = std::imag(t_complex_phi_out(1, 0)) / eps;
            t_diff_phi(1, 1, dd) = std::imag(t_complex_phi_out(1, 1)) / eps;
            t_diff_phi(1, 2, dd) = std::imag(t_complex_phi_out(1, 2)) / eps;
            t_diff_phi(2, 0, dd) = std::imag(t_complex_phi_out(2, 0)) / eps;
            t_diff_phi(2, 1, dd) = std::imag(t_complex_phi_out(2, 1)) / eps;
            t_diff_phi(2, 2, dd) = std::imag(t_complex_phi_out(2, 2)) / eps;
            if (debug) {
              MOFEM_LOG("SELF", Sev::noisy)
                  << "Gauss pt " << gg << " direction " << dd
                  << " diff_phi: " << t_diff_phi(0, 0, dd) << ", "
                  << t_diff_phi(0, 1, dd) << ", " << t_diff_phi(0, 2, dd)
                  << "; " << t_diff_phi(1, 0, dd) << ", "
                  << t_diff_phi(1, 1, dd) << ", " << t_diff_phi(1, 2, dd)
                  << "; " << t_diff_phi(2, 0, dd) << ", "
                  << t_diff_phi(2, 1, dd) << ", " << t_diff_phi(2, 2, dd);
            }
            ++t_complex_phi_out;
            ++t_diff_phi;
          }
        }
 
        if (debug) {
          auto t_phi = FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3>{
              &phi(0, 0), &phi(0, 1), &phi(0, 2), &phi(0, 3), &phi(0, 4),
              &phi(0, 5), &phi(0, 6), &phi(0, 7), &phi(0, 8)};
          auto t_complex_phi_out =
              FTensor::Tensor2<FTensor::PackPtr<std::complex<double> *, 9>, 3,
                               3>{
                  &temp_phi(0, 0), &temp_phi(0, 1), &temp_phi(0, 2),
                  &temp_phi(0, 3), &temp_phi(0, 4), &temp_phi(0, 5),
                  &temp_phi(0, 6), &temp_phi(0, 7), &temp_phi(0, 8)};
 
          for (auto gg = 0; gg < nb_gauss_pts; ++gg) {
            FTensor::Index<'i', 3> i;
            FTensor::Index<'j', 3> j;
            FTensor::Tensor2<std::complex<double>, 3, 3> t_diff;
            for (auto bb = 0; bb << nb_base_functions; ++bb) {
              t_diff(i, j) = t_complex_phi_out(i, j) - t_phi(i, j);
              auto nrm2_complex = t_diff(i, j) * t_diff(i, j);
              if (std::abs(nrm2_complex) > 1e-12) {
                MOFEM_LOG("SELF", Sev::warning)
                    << "Difference detected at gauss pt " << gg
                    << " for direction " << dd
                    << " nrm2 = " << std::abs(nrm2_complex);
                SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                        "Difference detected between complex step "
                        "differentiation and real valued differentiation");
              }
              ++t_phi;
              ++t_complex_phi_out;
            }
          }
        }
 
      
 
      }
 
      MoFEMFunctionReturn(0);
    };
 
    MatrixDouble &diff_phi = data.dataOnEntities[MBTET][0].getDiffN(USER_BASE);
    if constexpr (calculate_diff) {
      diff_phi.resize(nb_gauss_pts, 27 * nb_base_functions, false);
      CHKERR imagingray_fd_directive(diff_phi);
    }
 
    if (cachePhiPtr) {
      cachePhiPtr->get<0>() = order;
      cachePhiPtr->get<1>() = nb_gauss_pts;
      cachePhiPtr->get<2>().resize(phi.size1(), phi.size2(), false);
      noalias(cachePhiPtr->get<2>()) = phi;
      if constexpr (calculate_diff) {
        cachePhiPtr->get<3>().resize(diff_phi.size1(), diff_phi.size2(), false);
        noalias(cachePhiPtr->get<3>()) = diff_phi;
      }
    }
  }
 
  MoFEMFunctionReturn(0);
}