FieldEvaluator.cpp

Go to the documentation of this file.

/** \file FieldEvaluator.cpp
 * \brief Field evaluator
 *
 */
 
namespace MoFEM {
 
FieldEvaluatorInterface::FieldEvaluatorInterface(const MoFEM::Core &core)
    : cOre(const_cast<MoFEM::Core &>(core)) {
 
  if (!LogManager::checkIfChannelExist("FieldEvaluatorWorld")) {
    auto core_log = logging::core::get();
 
    core_log->add_sink(LogManager::createSink(LogManager::getStrmWorld(),
                                              "FieldEvaluatorWorld"));
    core_log->add_sink(LogManager::createSink(LogManager::getStrmSync(),
                                              "FieldEvaluatorSync"));
    core_log->add_sink(LogManager::createSink(LogManager::getStrmSelf(),
                                              "FieldEvaluatorSelf"));
 
    LogManager::setLog("FieldEvaluatorWorld");
    LogManager::setLog("FieldEvaluatorSync");
    LogManager::setLog("FieldEvaluatorSelf");
 
    MOFEM_LOG_TAG("FieldEvaluatorWorld", "FieldEvaluator");
    MOFEM_LOG_TAG("FieldEvaluatorSync", "FieldEvaluator");
    MOFEM_LOG_TAG("FieldEvaluatorSelf", "FieldEvaluator");
  }
 
  MOFEM_LOG("FieldEvaluatorWorld", Sev::noisy) << "Field evaluator intreface";
}
 
MoFEMErrorCode FieldEvaluatorInterface::query_interface(
    boost::typeindex::type_index type_index, UnknownInterface **iface) const {
  MoFEMFunctionBeginHot;
  *iface = const_cast<FieldEvaluatorInterface *>(this);
  MoFEMFunctionReturnHot(0);
}
 
template <int D>
MoFEMErrorCode FieldEvaluatorInterface::buildTreeImpl(
    boost::shared_ptr<SetIntegrationPtsMethodData> spd_ptr,
    const std::string finite_element, BitRefLevel bit, BitRefLevel mask) {
  CoreInterface &m_field = cOre;
  MoFEMFunctionBegin;
  EntityHandle fe_meshset;
  fe_meshset = m_field.get_finite_element_meshset(finite_element);
  Range entities;
  CHKERR m_field.get_moab().get_entities_by_dimension(fe_meshset, D, entities,
                                                      true);
  if (bit.any()) {
    CHKERR m_field.getInterface<BitRefManager>()->filterEntitiesByRefLevel(
        bit, mask, entities);
  }
  if (spd_ptr->treePtr) {
    spd_ptr->treePtr->reset_tree();
  } else {
    spd_ptr->treePtr = boost::make_shared<AdaptiveKDTree>(&m_field.get_moab());
  }
  CHKERR spd_ptr->treePtr->build_tree(entities, &spd_ptr->rooTreeSet);
  MoFEMFunctionReturn(0);
}
 
template <>
MoFEMErrorCode FieldEvaluatorInterface::buildTree<2>(
    boost::shared_ptr<SetIntegrationPtsMethodData> spd_ptr,
    const std::string finite_element, BitRefLevel bit, BitRefLevel mask) {
  return buildTreeImpl<2>(spd_ptr, finite_element, bit, mask);
}
 
template <>
MoFEMErrorCode FieldEvaluatorInterface::buildTree<3>(
    boost::shared_ptr<SetIntegrationPtsMethodData> spd_ptr,
    const std::string finite_element, BitRefLevel bit, BitRefLevel mask) {
  return buildTreeImpl<3>(spd_ptr, finite_element, bit, mask);
}
 
FieldEvaluatorInterface::SetIntegrationPtsMethod::SetIntegrationPtsMethod(
    boost::shared_ptr<SetIntegrationPtsMethodData> data_ptr)
    : dataPtr(data_ptr) {}
 
MoFEMErrorCode FieldEvaluatorInterface::SetIntegrationPtsMethod::operator()(
    ForcesAndSourcesCore *fe_raw_ptr, int order_row, int order_col,
    int order_data) {
  MoFEMFunctionBegin;
 
  if (auto data_ptr = dataPtr.lock()) {
 
    const auto nb_eval_points = data_ptr->nbEvalPoints;
    const auto &shape_functions = data_ptr->shapeFunctions;
    const auto &eval_pointentity_handle = data_ptr->evalPointEntityHandle;
 
    if (auto fe_ptr = data_ptr->feMethodPtr) {
 
      auto &fe = static_cast<ForcesAndSourcesCore &>(*fe_ptr);
 
#ifndef NDEBUG
      if (fe_ptr.get() != fe_raw_ptr)
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Wrong FE ptr");
#endif
 
      const auto fe_ent = fe.numeredEntFiniteElementPtr->getEnt();
      const auto fe_type = type_from_handle(fe_ent);
      const auto fe_dim = moab::CN::Dimension(fe_type);
 
      auto &gauss_pts = fe.gaussPts;
      int nb_gauss_pts;
 
      if (fe_dim == 3) {
        gauss_pts.resize(4, nb_eval_points, false);
        FTensor::Tensor1<FTensor::PackPtr<const double *, 4>, 4> t_shape{
            &shape_functions(0, 0), &shape_functions(0, 1),
            &shape_functions(0, 2), &shape_functions(0, 3)};
        FTensor::Tensor1<FTensor::PackPtr<double *, 1>, 3> t_gauss_pts{
            &gauss_pts(0, 0), &gauss_pts(1, 0), &gauss_pts(2, 0)};
 
        nb_gauss_pts = 0;
        for (int nn = 0; nn != nb_eval_points; ++nn) {
          if (eval_pointentity_handle[nn] == fe_ent) {
            for (const int i : {0, 1, 2}) {
              t_gauss_pts(i) = t_shape(i + 1);
            }
            gauss_pts(3, nb_gauss_pts) = nn;
            ++t_gauss_pts;
            ++nb_gauss_pts;
          }
          ++t_shape;
        }
        gauss_pts.resize(4, nb_gauss_pts, true);
      } else if (fe_dim == 2) {
        gauss_pts.resize(3, nb_eval_points, false);
        FTensor::Tensor1<FTensor::PackPtr<const double *, 3>, 3> t_shape{
            &shape_functions(0, 0), &shape_functions(0, 1),
            &shape_functions(0, 2)};
        FTensor::Tensor1<FTensor::PackPtr<double *, 1>, 2> t_gauss_pts{
            &gauss_pts(0, 0), &gauss_pts(1, 0)};
        nb_gauss_pts = 0;
        for (int nn = 0; nn != nb_eval_points; ++nn) {
          if (eval_pointentity_handle[nn] == fe_ent) {
            for (const int i : {0, 1}) {
              t_gauss_pts(i) = t_shape(i + 1);
            }
            gauss_pts(2, nb_gauss_pts) = nn;
            ++t_gauss_pts;
            ++nb_gauss_pts;
          }
          ++t_shape;
        }
        gauss_pts.resize(3, nb_gauss_pts, true);
      } else {
        SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
                "Dimension not implemented");
      }
 
#ifndef NDEBUG
      MOFEM_LOG("SELF", Sev::noisy)
          << "nbEvalOPoints / nbGau_sspt_s: " << nb_eval_points << " / "
          << nb_gauss_pts;
      MOFEM_LOG("SELF", Sev::noisy) << "gauss pts: " << gauss_pts;
#endif
 
    } else
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Pointer to element does not exists");
 
  } else
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Pointer to data does not exists");
 
  MoFEMFunctionReturn(0);
}
 
template <int D>
MoFEMErrorCode FieldEvaluatorInterface::evalFEAtThePointImpl(
    const double *const point, const double distance, const std::string problem,
    const std::string finite_element,
    boost::shared_ptr<SetIntegrationPtsMethodData> data_ptr, int lower_rank,
    int upper_rank, BitRefLevel bit, BitRefLevel mask,
    boost::shared_ptr<CacheTuple> cache_ptr, MoFEMTypes bh,
    VERBOSITY_LEVELS verb) {
  CoreInterface &m_field = cOre;
  MoFEMFunctionBegin;
 
  std::vector<EntityHandle> leafs_out;
  CHKERR data_ptr->treePtr->distance_search(point, distance, leafs_out);
  Range tree_ents;
  for (auto lit : leafs_out) //{
    CHKERR m_field.get_moab().get_entities_by_dimension(lit, D, tree_ents,
                                                        true);
 
  if (verb >= NOISY)
    MOFEM_LOG("SELF", Sev::noisy) << "tree entities: " << tree_ents;
 
  data_ptr->evalPointEntityHandle.resize(data_ptr->nbEvalPoints);
  std::fill(data_ptr->evalPointEntityHandle.begin(),
            data_ptr->evalPointEntityHandle.end(), 0);
 
  std::vector<double> local_coords;
  std::vector<double> shape;
 
  auto nb_eval_points = data_ptr->nbEvalPoints;
 
  for (auto tet : tree_ents) {
    if ((type_from_handle(tet) != moab::MBTRI) &&
        (type_from_handle(tet) != moab::MBTET)) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
              "Wrong element type, Field Evaluator not implemented for Quads "
              "and Hexes");
    }
    const EntityHandle *conn;
    int num_nodes;
    CHKERR m_field.get_moab().get_connectivity(tet, conn, num_nodes, true);
 
    if constexpr (D == 3) {
 
      local_coords.resize(3 * nb_eval_points);
      shape.resize(4 * nb_eval_points);
 
      std::array<double, 12> coords;
      CHKERR m_field.get_moab().get_coords(conn, num_nodes, coords.data());
 
      CHKERR Tools::getLocalCoordinatesOnReferenceFourNodeTet(
          coords.data(), &data_ptr->evalPoints[0], nb_eval_points,
          &local_coords[0]);
      CHKERR Tools::shapeFunMBTET<3>(&shape[0], &local_coords[0],
                                     &local_coords[1], &local_coords[2],
                                     nb_eval_points);
 
      FTensor::Index<'i', 4> i4;
      FTensor::Tensor1<FTensor::PackPtr<double *, 4>, 4> t_shape{
          &shape[0], &shape[1], &shape[2], &shape[3]};
      FTensor::Tensor1<FTensor::PackPtr<double *, 4>, 4> t_shape_data{
          &data_ptr->shapeFunctions(0, 0), &data_ptr->shapeFunctions(0, 1),
          &data_ptr->shapeFunctions(0, 2), &data_ptr->shapeFunctions(0, 3)};
 
      FTensor::Index<'j', 3> j3;
      FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_local{
          &local_coords[0], &local_coords[1], &local_coords[2]};
      FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_local_data{
          &(data_ptr->localCoords(0, 0)), &(data_ptr->localCoords(0, 1)),
          &(data_ptr->localCoords(0, 2))};
 
      for (int n = 0; n != nb_eval_points; ++n) {
 
        const double eps = data_ptr->eps;
        if (t_shape(0) >= 0 - eps && t_shape(0) <= 1 + eps &&
 
            t_shape(1) >= 0 - eps && t_shape(1) <= 1 + eps &&
 
            t_shape(2) >= 0 - eps && t_shape(2) <= 1 + eps &&
 
            t_shape(3) >= 0 - eps && t_shape(3) <= 1 + eps) {
 
          data_ptr->evalPointEntityHandle[n] = tet;
          t_shape_data(i4) = t_shape(i4);
          t_local_data(j3) = t_local(j3);
        }
 
        ++t_shape;
        ++t_shape_data;
        ++t_local;
        ++t_local_data;
      }
    }
 
    if constexpr (D == 2) {
 
      local_coords.resize(2 * nb_eval_points);
      shape.resize(3 * nb_eval_points);
 
      std::array<double, 9> coords;
      CHKERR m_field.get_moab().get_coords(conn, num_nodes, coords.data());
 
      CHKERR Tools::getLocalCoordinatesOnReferenceThreeNodeTri(
          coords.data(), &data_ptr->evalPoints[0], nb_eval_points,
          &local_coords[0]);
      CHKERR Tools::shapeFunMBTRI<2>(&shape[0], &local_coords[0],
                                     &local_coords[1], nb_eval_points);
 
      FTensor::Index<'i', 3> i3;
      FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_shape{
          &shape[0], &shape[1], &shape[2]};
      FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_shape_data{
          &data_ptr->shapeFunctions(0, 0), &data_ptr->shapeFunctions(0, 1),
          &data_ptr->shapeFunctions(0, 2)};
 
      FTensor::Index<'j', 2> j2;
      FTensor::Tensor1<FTensor::PackPtr<double *, 2>, 2> t_local{
          &local_coords[0], &local_coords[1]};
      data_ptr->localCoords.resize(nb_eval_points, 2);
      FTensor::Tensor1<FTensor::PackPtr<double *, 2>, 2> t_local_data{
          &(data_ptr->localCoords(0, 0)), &(data_ptr->localCoords(0, 1))};
 
      for (int n = 0; n != nb_eval_points; ++n) {
 
        const double eps = data_ptr->eps;
        if (t_shape(0) >= 0 - eps && t_shape(0) <= 1 + eps &&
 
            t_shape(1) >= 0 - eps && t_shape(1) <= 1 + eps &&
 
            t_shape(2) >= 0 - eps && t_shape(2) <= 1 + eps) {
 
          data_ptr->evalPointEntityHandle[n] = tet;
          t_shape_data(i3) = t_shape(i3);
          t_local_data(j2) = t_local(j2);
        }
 
        ++t_shape;
        ++t_shape_data;
        ++t_local;
        ++t_local_data;
      }
    }
  }
 
  const Problem *prb_ptr;
  CHKERR m_field.get_problem(problem, &prb_ptr);
  boost::shared_ptr<NumeredEntFiniteElement_multiIndex> numered_fes(
      new NumeredEntFiniteElement_multiIndex());
 
  Range in_tets;
  in_tets.insert_list(data_ptr->evalPointEntityHandle.begin(),
                      data_ptr->evalPointEntityHandle.end());
  in_tets = in_tets.subset_by_dimension(D);
 
  if (verb >= NOISY)
    MOFEM_LOG("SELF", Sev::noisy) << "in tets: " << in_tets << endl;
 
  auto fe_miit =
      m_field.get_finite_elements()->get<FiniteElement_name_mi_tag>().find(
          finite_element);
  if (fe_miit !=
      m_field.get_finite_elements()->get<FiniteElement_name_mi_tag>().end()) {
 
    for (auto peit = in_tets.pair_begin(); peit != in_tets.pair_end(); ++peit) {
 
      auto lo =
          prb_ptr->numeredFiniteElementsPtr->get<Unique_mi_tag>().lower_bound(
              EntFiniteElement::getLocalUniqueIdCalculate(
                  peit->first, (*fe_miit)->getFEUId()));
      auto hi =
          prb_ptr->numeredFiniteElementsPtr->get<Unique_mi_tag>().upper_bound(
              EntFiniteElement::getLocalUniqueIdCalculate(
                  peit->second, (*fe_miit)->getFEUId()));
 
      if (bit.any()) {
        if (verb >= VERY_NOISY)
          MOFEM_LOG("FieldEvaluatorSelf", Sev::noisy);
        for (auto it = lo; it != hi; ++it) {
          auto fe_bit = (*it)->getBitRefLevel();
          if ((fe_bit & mask) != fe_bit)
            continue;
          if ((fe_bit & bit).none())
            continue;
          numered_fes->insert(*it);
          if (verb >= VERY_NOISY)
            MOFEM_LOG("FieldEvaluatorSelf", Sev::noisy) << **it << endl;
        }
      } else {
        numered_fes->insert(lo, hi);
        if (verb >= VERY_NOISY) {
          MOFEM_LOG("FieldEvaluatorSelf", Sev::noisy)
              << "numbered elements:" << std::endl;
          for (; lo != hi; ++lo)
            MOFEM_LOG("FieldEvaluatorSelf", Sev::noisy) << **lo << endl;
        }
      }
    }
 
    if (auto fe_ptr = data_ptr->feMethodPtr) {
 
      MOFEM_LOG_C("FieldEvaluatorSelf", Sev::noisy,
                  "Number elements %zu to evaluate at proc %d",
                  numered_fes->size(), m_field.get_comm_rank());
 
      if (!cache_ptr) {
        cache_ptr = boost::make_shared<CacheTuple>();
        CHKERR m_field.cache_problem_entities(prb_ptr->getName(), cache_ptr);
 
        MOFEM_LOG("FieldEvaluatorSelf", Sev::noisy)
            << "If you call that function many times in the loop consider to "
               "set cache_ptr outside of the loop. Otherwise code can be slow.";
      }
 
      CHKERR m_field.loop_finite_elements(prb_ptr, finite_element, *fe_ptr,
                                          lower_rank, upper_rank, numered_fes,
                                          bh, cache_ptr, verb);
 
    } else
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Pointer to element does not exists");
  }
 
  MoFEMFunctionReturn(0);
}
 
template <>
MoFEMErrorCode FieldEvaluatorInterface::evalFEAtThePoint<2>(
    const double *const point, const double distance, const std::string problem,
    const std::string finite_element,
    boost::shared_ptr<SetIntegrationPtsMethodData> data_ptr, int lower_rank,
    int upper_rank, BitRefLevel bit, BitRefLevel mask,
    boost::shared_ptr<CacheTuple> cache_ptr, MoFEMTypes bh,
    VERBOSITY_LEVELS verb) {
  return evalFEAtThePointImpl<2>(point, distance, problem, finite_element,
                                 data_ptr, lower_rank, upper_rank, bit, mask,
                                 cache_ptr, bh, verb);
}
 
template <>
MoFEMErrorCode FieldEvaluatorInterface::evalFEAtThePoint<3>(
    const double *const point, const double distance, const std::string problem,
    const std::string finite_element,
    boost::shared_ptr<SetIntegrationPtsMethodData> data_ptr, int lower_rank,
    int upper_rank, BitRefLevel bit, BitRefLevel mask,
    boost::shared_ptr<CacheTuple> cache_ptr, MoFEMTypes bh,
    VERBOSITY_LEVELS verb) {
  return evalFEAtThePointImpl<3>(point, distance, problem, finite_element,
                                 data_ptr, lower_rank, upper_rank, bit, mask,
                                 cache_ptr, bh, verb);
}
 
template <int D>
ForcesAndSourcesCore::UserDataOperator *
FieldEvaluatorInterface::getDataOperatorImpl(
    FieldEvaluatorInterface::MatDataPts mat_data_pts,
    const std::string finite_element,
    boost::shared_ptr<SetIntegrationPtsMethodData> data_ptr, int lower_rank,
    int upper_rank, BitRefLevel bit, BitRefLevel mask,
    MoFEMTypes bh, VERBOSITY_LEVELS verb, DataLayout data_layout) {
 
  using BaseOP = ForcesAndSourcesCore::UserDataOperator;
 
  struct Op : public BaseOP {
    Op(
 
        FieldEvaluatorInterface::MatDataPts mat_data_pts,
        const std::string finite_element,
        boost::shared_ptr<SetIntegrationPtsMethodData> data_ptr, int lower_rank,
        int upper_rank, BitRefLevel bit, BitRefLevel mask, double eps,
        MoFEMTypes bh, VERBOSITY_LEVELS verb, DataLayout data_layout
 
        )
        : BaseOP(NOSPACE, BaseOP::OPSPACE), 
          matDataPts(mat_data_pts), finiteElement(finite_element),
          dataPtr(data_ptr), lowerRank(lower_rank), upperRank(upper_rank),
          bitRefLevel(bit), maskRefLevel(mask), eps(eps), bhType(bh),
          verbosity(verb), dataLayout(data_layout) {}
 
    MoFEMErrorCode doWork(int side, EntityType type,
                          EntitiesFieldData::EntData &data) {
      MoFEMFunctionBegin;
 
      MOFEM_LOG_CHANNEL("SELF");
      MOFEM_LOG_TAG("SELF", "FieldEvaluator");
 
      MoFEM::Interface &m_field = this->getPtrFE()->mField;
      auto field_eval_ptr = m_field.getInterface<FieldEvaluatorInterface>();
 
      auto nb_integration_pts = getGaussPts().size2();
 
      auto prb_name = getPtrFE()->problemPtr->getName();
      auto cache_ptr = getPtrFE()->getCacheWeakPtr();
 
      auto coords_at_pts = getCoordsAtGaussPts();
      for (int gg = 0; gg != nb_integration_pts; ++gg) {
        double *coords_ptr = &coords_at_pts(gg, 0);
 
        dataPtr->setEvalPoints(coords_ptr, 1);
        CHKERR field_eval_ptr->evalFEAtThePoint<D>(
            coords_ptr, eps, prb_name, finiteElement, dataPtr, lowerRank,
            upperRank, bitRefLevel, maskRefLevel, cache_ptr.lock(), MF_EXIST,
            QUIET);
 
        // Get number of found finite elements, that consist gauss point.
        // E.g. for edge flip expected is that only one element will be found.
        auto loop_size = dataPtr->feMethodPtr->getLoopSize();
        if (loop_size != 1) {
          if (verbosity >= QUIET) {
            MOFEM_LOG_ATTRIBUTES("SELF",
                                 LogManager::BitLineID | LogManager::BitScope);
            MOFEM_LOG("SELF", Sev::warning)
                << "Field evaluator wrong number of elements found: "
                << loop_size;
          }
        }
 
        for (auto &m : matDataPts) {
 
          const auto nb_data_rows = m.first->size1();
          const auto nb_data_cols = m.first->size2();
 
#ifndef NDEBUG
          const auto nb_data_gauss_pts =
              dataLayout == DataLayout::GaussByCoeffs ? nb_data_rows
                                                      : nb_data_cols;
          if (nb_data_gauss_pts != loop_size) {
            SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                    "Inconsistent source data layout for size (%ld, %ld), "
                    "loop size %d, and requested layout %s",
                    nb_data_rows, nb_data_cols, loop_size,
                    dataLayout == DataLayout::GaussByCoeffs ? "GaussByCoeffs"
                                                            : "CoeffsByGauss");
          }
#endif
 
          const auto nb_data_coeffs =
              dataLayout == DataLayout::GaussByCoeffs ? nb_data_cols
                                                      : nb_data_rows;
 
          if (dataLayout == DataLayout::GaussByCoeffs) {
            m.second->resize(nb_integration_pts, nb_data_coeffs, false);
          } else {
            m.second->resize(nb_data_coeffs, nb_integration_pts, false);
          }
          if (gg == 0)
            m.second->clear();
          for (auto ll = 0; ll != loop_size; ++ll) {
            for (int dd = 0; dd != nb_data_coeffs; ++dd) {
              const auto value = dataLayout == DataLayout::GaussByCoeffs
                                     ? (*m.first)(ll, dd)
                                     : (*m.first)(dd, ll);
              if (dataLayout == DataLayout::GaussByCoeffs) {
                (*m.second)(gg, dd) += value / loop_size;
              } else {
                (*m.second)(dd, gg) += value / loop_size;
              }
            }
          }
        }
      }
 
      MOFEM_LOG_CHANNEL("SELF");
 
      MoFEMFunctionReturn(0);
    }
 
  private:
    FieldEvaluatorInterface::MatDataPts matDataPts;
    std::string finiteElement;
    boost::shared_ptr<SetIntegrationPtsMethodData> dataPtr;
    int lowerRank;
    int upperRank;
    BitRefLevel bitRefLevel;
    BitRefLevel maskRefLevel;
    double eps;
    MoFEMTypes bhType;
    VERBOSITY_LEVELS verbosity;
    DataLayout dataLayout;
  };
 
  double min[3];
  double max[3];
  unsigned int dep;
  CHKERR data_ptr->treePtr->get_info(data_ptr->rooTreeSet, min, max, dep);
  double dist = std::sqrt(pow(max[0] - min[0], 2) + pow(max[1] - min[1], 2) +
                          pow(max[2] - min[2], 2));
 
  return new Op(mat_data_pts, finite_element, data_ptr, lower_rank,
                upper_rank, bit, mask, 1e-6 * dist, bh, verb, data_layout);
}
 
FieldEvaluatorInterface::SetIntegrationPtsMethodData::
    SetIntegrationPtsMethodData(
        boost::shared_ptr<MoFEM::ForcesAndSourcesCore> fe_method_ptr,
        const double *eval_points, const int nb_eval_points, const double eps,
        VERBOSITY_LEVELS verb)
    : feMethodPtr(fe_method_ptr), evalPoints(eval_points),
      nbEvalPoints(nb_eval_points), eps(eps), verb(verb) {
  localCoords.resize(nbEvalPoints, 3);
  shapeFunctions.resize(nbEvalPoints, 4);
}
 
void FieldEvaluatorInterface::SetIntegrationPtsMethodData::setEvalPoints(
    const double *ptr, const int nb_eval_points) {
  evalPoints = ptr;
  nbEvalPoints = nb_eval_points;
  localCoords.resize(nbEvalPoints, 3, false);
  shapeFunctions.resize(nbEvalPoints, 4, false);
}
 
template <>
ForcesAndSourcesCore::UserDataOperator *
FieldEvaluatorInterface::getDataOperator<2>(
    FieldEvaluatorInterface::MatDataPts mat_data_pts,
    const std::string finite_element,
    boost::shared_ptr<SetIntegrationPtsMethodData> data_ptr, int lower_rank,
    int upper_rank, BitRefLevel bit, BitRefLevel mask, MoFEMTypes bh,
    VERBOSITY_LEVELS verb, DataLayout data_layout) {
  return getDataOperatorImpl<2>(mat_data_pts, finite_element,
                                data_ptr, lower_rank, upper_rank, bit, mask,
                                bh, verb, data_layout);
}
 
template <>
ForcesAndSourcesCore::UserDataOperator *
FieldEvaluatorInterface::getDataOperator<3>(
    FieldEvaluatorInterface::MatDataPts mat_data_pts,
    const std::string finite_element,
    boost::shared_ptr<SetIntegrationPtsMethodData> data_ptr, int lower_rank,
    int upper_rank, BitRefLevel bit, BitRefLevel mask, MoFEMTypes bh,
    VERBOSITY_LEVELS verb, DataLayout data_layout) {
  return getDataOperatorImpl<3>(mat_data_pts, finite_element,
                                data_ptr, lower_rank, upper_rank, bit, mask,
                                bh, verb, data_layout);
}
 
} // namespace MoFEM