FieldEvaluator.hpp

Go to the documentation of this file.

/** \file FieldEvaluator.hpp
 * \brief Field Evaluator
 *
 * Evaluate field at given coordinate
 *
 *
 * \ingroup field_evaluator
 */
 
#ifndef __FIELD_EVALUATOR_HPP__
  #define __FIELD_EVALUATOR_HPP__
 
  #include "UnknownInterface.hpp"
 
namespace MoFEM {
 
/** \brief Field evaluator interface
 
  * \ingroup field_evaluator
  */
struct FieldEvaluatorInterface : public UnknownInterface {
 
  struct SetIntegrationPtsMethodData;
  using SetPtsData = SetIntegrationPtsMethodData;
 
  /**
   * @brief Default evaluator for setting integration points
   *
   */
  struct SetIntegrationPtsMethod;
 
  MoFEMErrorCode query_interface(boost::typeindex::type_index type_index,
                                 UnknownInterface **iface) const;
 
  FieldEvaluatorInterface(const MoFEM::Core &core);
 
  /**
   * @brief Get the Data object
   *
   * Pack pointers with data structures for field evaluator and finite
   * element. Function return shared pointer if returned shared pointer
   * is reset; all data are destroyed. The idea is to pack all data in
   * one structure, create shared pointer to it and return aliased shared
   * pointer to one of the elements of the data structure. It is a bit
   * complicated, but has great flexibility.
   *
   * @tparam VE
   * @tparam SetIntegrationPtsMethodData
   * @tparam SetIntegrationPtsMethod
   * @param ptr
   * @param nb_eval_points
   * @param eps
   * @param verb
   * @return boost::shared_ptr<SPD>
   */
  template <typename VE, typename SPD = SetIntegrationPtsMethodData,
            typename SP = SetIntegrationPtsMethod>
  boost::shared_ptr<SPD>
  getData(const double *ptr = nullptr, const int nb_eval_points = 0,
          const double eps = 1e-12, VERBOSITY_LEVELS verb = QUIET);
 
  /**
   * @brief Build spatial tree
   *
   * @param finite_element finite element name
   * @return MoFEMErrorCode
   */
  template <int D>
  inline MoFEMErrorCode
  buildTree(boost::shared_ptr<SetIntegrationPtsMethodData> spd_ptr,
            const std::string finite_element) {
    return buildTree<D>(spd_ptr, finite_element, BitRefLevel(), BitRefLevel());
  }
 
  /**
   * @brief Build spatial tree
   *
   * @param finite_element finite element name
   * @param bit bit reference level
   * @param mask mask reference level
   * @return MoFEMErrorCode
   */
  template <int D>
  MoFEMErrorCode
  buildTree(boost::shared_ptr<SetIntegrationPtsMethodData> spd_ptr,
            const std::string finite_element, BitRefLevel bit,
            BitRefLevel mask);
 
  /**
   * @brief Evaluate field at arbitrary position
   *
   * \code
 
    std::array<double, 3> point = {0, 0, 0};
    const double dist = 0.3;
    std::array<double, 6> eval_points = {-1.,  -1., -1., 1., 1., 1. };
 
    using VolEle = VolumeElementForcesAndSourcesCore;
    auto data_ptr =
      m_field.getInterface<FieldEvaluatorInterface>()->
      getData<VolEle>(point.data(), point.size/3);
 
    if(auto vol_ele = data_ptr->feMethod.lock()) {
      // push operators to finite element instance, e.g.
      vol_ele->getOpPtrVector().push_back(new MyOp());
      // iterate over elemnts with evaluated points
      auto cache_ptr = boost::make_shared<CacheTuple>();
      CHKERR m_field.cache_problem_entities(prb_ptr->getName(), cache_ptr);
      // SPACE_DIM is 3 in this example
      CHKERR m_field.getInterface<FieldEvaluatorInterface>()
        ->evalFEAtThePoint<3>(point.data(), dist, prb_ptr->getName(),
                             "FINITE_ELEMENT_NAME",
                             data_ptr, m_field.get_comm_rank(),
                             m_field.get_comm_rank(), cache_ptr);
    }
 
   * \endcode
   *
   * @param point point used to find tetrahedrons
   * @param distance distance from the point where tetrahedrons are searched
   * @param problem problem name
   * @param finite_element finite element name
   * @param data_ptr pointer to data abut gauss points
   * @param lower_rank lower processor
   * @param upper_rank upper process
   * @param cache_ptr cache or problem entities
   * @param bh control looping over entities, e.g. throwing error if element not
   found
   * @param verb verbosity level
   * @return MoFEMErrorCode
   */
  template <int D>
  inline MoFEMErrorCode
  evalFEAtThePoint(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,
                   boost::shared_ptr<CacheTuple> cache_ptr,
                   MoFEMTypes bh = MF_EXIST, VERBOSITY_LEVELS verb = QUIET) {
    return evalFEAtThePoint<D>(point, distance, problem, finite_element,
                               data_ptr, lower_rank, upper_rank, BitRefLevel(),
                               BitRefLevel(), cache_ptr, bh, verb);
  }
 
  using MatDataPts = std::vector<
 
      std::pair<boost::shared_ptr<MatrixDouble>,
                boost::shared_ptr<MatrixDouble>>
 
      >;
 
  /**
   * @brief Evaluate field values at integration points of finite elements
   * to which operator is pushed
   * 
   * @param mat_data_pts vector of pairs of pointers to matrices with field
   * values at integration points
   * @param finite_element finite element name on which field is evaluated
   * @param data_ptr data about field evaluator
   * @param lower_rank lower processor
   * @param upper_rank upper processor
   * @param bit bit reference level
   * @param mask mask reference level
   * @param bh control looping over entities, e.g. throwing error if element not
   * found
   * @param verb verbosity level
   * @param data_layout layout of evaluated output matrices
   * @return ForcesAndSourcesCore::UserDataOperator*
   */
  template <int D>
  ForcesAndSourcesCore::UserDataOperator *
  getDataOperator(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 = MF_EXIST,
                  VERBOSITY_LEVELS verb = QUIET,
                  DataLayout data_layout = DataLayout::GaussByCoeffs);
 
  /** @brief Evaluate field at arbitrary position
   * 
   * @param point point used to find tetrahedrons
   * @param distance distance from the point where tetrahedrons are searched
   * @param problem problem name
   * @param finite_element finite element name
   * @param data_ptr pointer to data abut gauss points
   * @param lower_rank lower processor
   * @param upper_rank upper process
   * @param bit bit reference level
   * @param mask mask reference level
   * @param cache_ptr cache or problem entities
   * @param bh control looping over entities, e.g. throwing error if element not
   found
   * @param verb verbosity level
   * @return MoFEMErrorCode
   */
  template <int D>
  MoFEMErrorCode
  evalFEAtThePoint(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 = MF_EXIST, VERBOSITY_LEVELS verb = QUIET);
 
private:
  MoFEM::Core &cOre;
 
  template <int D>
  MoFEMErrorCode
  buildTreeImpl(boost::shared_ptr<SetIntegrationPtsMethodData> spd_ptr,
                const std::string finite_element, BitRefLevel bit,
                BitRefLevel mask);
 
  template <int D>
  MoFEMErrorCode 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 = MF_EXIST,
      VERBOSITY_LEVELS verb = QUIET);
 
  template <int D>
  ForcesAndSourcesCore::UserDataOperator *
  getDataOperatorImpl(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 = MF_EXIST,
                      VERBOSITY_LEVELS verb = QUIET,
                      DataLayout data_layout = DataLayout::GaussByCoeffs);
};
 
template <>
MoFEMErrorCode FieldEvaluatorInterface::buildTree<2>(
    boost::shared_ptr<SetIntegrationPtsMethodData> spd_ptr,
    const std::string finite_element, BitRefLevel bit, BitRefLevel mask);
 
template <>
MoFEMErrorCode FieldEvaluatorInterface::buildTree<3>(
    boost::shared_ptr<SetIntegrationPtsMethodData> spd_ptr,
    const std::string finite_element, BitRefLevel bit, BitRefLevel mask);
 
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);
 
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);
 
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);
 
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);
 
struct FieldEvaluatorInterface::SetIntegrationPtsMethodData {
 
  SetIntegrationPtsMethodData() = delete;
 
  /**
     * @brief Set the Gauss Pts data
     *
     * @param fe_method_ptr pointer to finite element instance
     * @param eval_points pointer to array with evaluation points
     * @param nb_eval_points number of evaluated points
     * @param eps tolerance used to find if point is in the element
     * @param verb
     */
  SetIntegrationPtsMethodData(
      boost::shared_ptr<MoFEM::ForcesAndSourcesCore> fe_method_ptr,
      const double *eval_points, const int nb_eval_points, const double eps,
      VERBOSITY_LEVELS verb = QUIET);
 
  void setEvalPoints(const double *ptr, const int nb_eval_points);
 
  boost::shared_ptr<MoFEM::ForcesAndSourcesCore> feMethodPtr;
  EntityHandle rooTreeSet;
  boost::shared_ptr<AdaptiveKDTree> treePtr;
 
protected:
  const double *evalPoints;
  int nbEvalPoints;
  double eps;
  VERBOSITY_LEVELS verb;
 
  MatrixDouble localCoords;
  MatrixDouble shapeFunctions;
  std::vector<EntityHandle> evalPointEntityHandle;
 
  friend FieldEvaluatorInterface;
  friend struct SetIntegrationPtsMethod;
};
 
struct FieldEvaluatorInterface::SetIntegrationPtsMethod {
  SetIntegrationPtsMethod() = delete;
  SetIntegrationPtsMethod(
      boost::shared_ptr<SetIntegrationPtsMethodData> data_ptr);
  MoFEMErrorCode operator()(ForcesAndSourcesCore *fe_raw_ptr, int order_row,
                            int order_col, int order_data);
 
private:
  boost::weak_ptr<SetIntegrationPtsMethodData> dataPtr;
};
 
/*
  template <typename VE, typename SPD = SetIntegrationPtsMethodData,
            typename SP = SetIntegrationPtsMethod>
*/
template <typename VE, typename SPD, typename SP>
boost::shared_ptr<SPD>
FieldEvaluatorInterface::getData(const double *ptr, const int nb_eval_points,
                                 const double eps, VERBOSITY_LEVELS verb) {
  MoFEM::Interface &m_field = cOre;
  auto ve_method_ptr = boost::make_shared<VE>(m_field);
  auto spd_ptr =
      boost::make_shared<SPD>(ve_method_ptr, ptr, nb_eval_points, eps, verb);
  ve_method_ptr->setRuleHook = SP(spd_ptr);
  return spd_ptr;
}
 
} // namespace MoFEM
 
#endif // __FIELD_EVALUATOR_HPP__
 
/**
 * \defgroup field_evaluator Field Evaluator
 * \brief Evaluate field at the point
 *
 * \ingroup mofem
 */