MooseFunctor.h

Go to the documentation of this file.

//* This file is part of the MOOSE framework
//* https://mooseframework.inl.gov
//*
//* All rights reserved, see COPYRIGHT for full restrictions
//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
//*
//* Licensed under LGPL 2.1, please see LICENSE for details
//* https://www.gnu.org/licenses/lgpl-2.1.html
 
#pragma once
 
#include <tuple>
 
#include "MooseFunctorForward.h"
#include "MooseFunctorArguments.h"
#include "FaceArgInterface.h"
#include "MooseMesh.h"
#include "MooseTypes.h"
#include "MooseError.h"
#include "MooseUtils.h"
 
#include "libmesh/remote_elem.h"
#include "libmesh/tensor_tools.h"
 
#include "metaphysicl/ct_types.h"
 
#include <unordered_map>
#include <functional>
 
namespace Moose
{
enum class FunctorEvaluationKind
{
  Value,
  Gradient,
  Dot,
  GradDot
};
 
template <typename, FunctorEvaluationKind>
struct FunctorReturnType;
 
template <typename T>
struct FunctorReturnType<T, FunctorEvaluationKind::Value>
{
  typedef T type;
};
 
template <typename T>
struct FunctorReturnType<T, FunctorEvaluationKind::Gradient>
{
  typedef typename MetaPhysicL::ReplaceAlgebraicType<
      T,
      typename libMesh::TensorTools::IncrementRank<
          typename MetaPhysicL::ValueType<T>::type>::type>::type type;
};
 
template <typename T>
struct FunctorReturnType<T, FunctorEvaluationKind::Dot>
{
  typedef T type;
};
 
template <typename T>
struct FunctorReturnType<T, FunctorEvaluationKind::GradDot>
{
  typedef typename FunctorReturnType<T, FunctorEvaluationKind::Gradient>::type type;
};
 
template <FunctorEvaluationKind>
struct FunctorGradientEvaluationKind;
 
template <>
struct FunctorGradientEvaluationKind<FunctorEvaluationKind::Value>
{
  static constexpr FunctorEvaluationKind value = FunctorEvaluationKind::Gradient;
};
 
template <>
struct FunctorGradientEvaluationKind<FunctorEvaluationKind::Dot>
{
  static constexpr FunctorEvaluationKind value = FunctorEvaluationKind::GradDot;
};
 
class FunctorAbstract : public FaceArgInterface
{
public:
  virtual void residualSetup() = 0;
  virtual void jacobianSetup() = 0;
  virtual void timestepSetup() = 0;
  virtual void customSetup(const ExecFlagType & exec_type) = 0;
};
 
template <typename T>
class FunctorBase : public FunctorAbstract
{
public:
  using FunctorType = FunctorBase<T>;
  using ValueType = T;
  using GradientType = typename FunctorReturnType<T, FunctorEvaluationKind::Gradient>::type;
  using DotType = ValueType;
 
  virtual ~FunctorBase() = default;
  FunctorBase(const MooseFunctorName & name,
              const std::set<ExecFlagType> & clearance_schedule = {EXEC_ALWAYS})
    : _always_evaluate(true), _functor_name(name)
  {
    setCacheClearanceSchedule(clearance_schedule);
  }
 
#ifdef MOOSE_KOKKOS_ENABLED
  FunctorBase(const FunctorBase<T> &, const Moose::Kokkos::FunctorCopy &) {}
#endif
 
  template <FunctorEvaluationKind FET, typename Space, typename State>
  typename FunctorReturnType<T, FET>::type genericEvaluate(const Space & r,
                                                           const State & state) const;
 
  const MooseFunctorName & functorName() const { return _functor_name; }
 
 
  ValueType operator()(const ElemArg & elem, const StateArg & state) const;
  ValueType operator()(const FaceArg & face, const StateArg & state) const;
  ValueType operator()(const ElemQpArg & qp, const StateArg & state) const;
  ValueType operator()(const ElemSideQpArg & qp, const StateArg & state) const;
  ValueType operator()(const ElemPointArg & elem_point, const StateArg & state) const;
  ValueType operator()(const NodeArg & node, const StateArg & state) const;
 
 
  GradientType gradient(const ElemArg & elem, const StateArg & state) const;
  GradientType gradient(const FaceArg & face, const StateArg & state) const;
  GradientType gradient(const ElemQpArg & qp, const StateArg & state) const;
  GradientType gradient(const ElemSideQpArg & qp, const StateArg & state) const;
  GradientType gradient(const ElemPointArg & elem_point, const StateArg & state) const;
  GradientType gradient(const NodeArg & node, const StateArg & state) const;
 
 
  DotType dot(const ElemArg & elem, const StateArg & state) const;
  DotType dot(const FaceArg & face, const StateArg & state) const;
  DotType dot(const ElemQpArg & qp, const StateArg & state) const;
  DotType dot(const ElemSideQpArg & qp, const StateArg & state) const;
  DotType dot(const ElemPointArg & elem_point, const StateArg & state) const;
  DotType dot(const NodeArg & node, const StateArg & state) const;
 
 
  GradientType gradDot(const ElemArg & elem, const StateArg & state) const;
  GradientType gradDot(const FaceArg & face, const StateArg & state) const;
  GradientType gradDot(const ElemQpArg & qp, const StateArg & state) const;
  GradientType gradDot(const ElemSideQpArg & qp, const StateArg & state) const;
  GradientType gradDot(const ElemPointArg & elem_point, const StateArg & state) const;
  GradientType gradDot(const NodeArg & node, const StateArg & state) const;
 
  virtual void residualSetup() override;
  virtual void jacobianSetup() override;
  virtual void timestepSetup() override;
  virtual void customSetup(const ExecFlagType & exec_type) override;
 
  void setCacheClearanceSchedule(const std::set<ExecFlagType> & clearance_schedule);
 
  virtual bool hasBlocks(SubdomainID /* id */) const
  {
    mooseError("Block restriction has not been implemented for functor " + functorName());
    return false;
  }
 
  virtual bool isExtrapolatedBoundaryFace(const FaceInfo &, const Elem *, const StateArg &) const
  {
    mooseError("not implemented");
  }
 
  bool isInternalFace(const FaceInfo &) const;
 
  virtual bool isConstant() const { return false; }
 
  virtual bool hasFaceSide(const FaceInfo & fi, const bool fi_elem_side) const override;
 
  void checkFace(const Moose::FaceArg & face) const;
 
  virtual bool supportsFaceArg() const = 0;
 
  virtual bool supportsElemSideQpArg() const = 0;
 
protected:
 
  virtual ValueType evaluate(const ElemArg & elem, const StateArg & state) const = 0;
 
  virtual ValueType evaluate(const FaceArg & face, const StateArg & state) const = 0;
 
  virtual ValueType evaluate(const ElemQpArg & qp, const StateArg & state) const = 0;
 
  virtual ValueType evaluate(const ElemSideQpArg & side_qp, const StateArg & state) const = 0;
 
  virtual ValueType evaluate(const ElemPointArg & elem_point, const StateArg & state) const = 0;
 
  virtual ValueType evaluate(const NodeArg & node, const StateArg & state) const = 0;
 
  virtual GradientType evaluateGradient(const ElemArg &, const StateArg &) const
  {
    mooseError("Element gradient not implemented for functor " + functorName());
  }
 
  virtual GradientType evaluateGradient(const FaceArg &, const StateArg &) const
  {
    mooseError("Face gradient not implemented for functor " + functorName());
  }
 
  virtual GradientType evaluateGradient(const ElemQpArg &, const StateArg &) const
  {
    mooseError("Element quadrature point gradient not implemented for functor " + functorName());
  }
 
  virtual GradientType evaluateGradient(const ElemSideQpArg &, const StateArg &) const
  {
    mooseError("Element side quadrature point gradient not implemented for functor " +
               functorName());
  }
 
  virtual GradientType evaluateGradient(const ElemPointArg &, const StateArg &) const
  {
    mooseError("Element-point gradient not implemented for functor " + functorName());
  }
 
  virtual GradientType evaluateGradient(const NodeArg &, const StateArg &) const
  {
    mooseError("Gradient at node not implemented for functor " + functorName());
  }
 
  virtual DotType evaluateDot(const ElemArg &, const StateArg &) const
  {
    mooseError("Element time derivative not implemented for functor " + functorName());
  }
 
  virtual DotType evaluateDot(const FaceArg &, const StateArg &) const
  {
    mooseError("Face time derivative not implemented for functor " + functorName());
  }
 
  virtual DotType evaluateDot(const ElemQpArg &, const StateArg &) const
  {
    mooseError("Element quadrature point time derivative not implemented for functor " +
               functorName());
  }
 
  virtual DotType evaluateDot(const ElemSideQpArg &, const StateArg &) const
  {
    mooseError("Element side quadrature point time derivative not implemented for functor " +
               functorName());
  }
 
  virtual DotType evaluateDot(const ElemPointArg &, const StateArg &) const
  {
    mooseError("Element-point time derivative not implemented for functor " + functorName());
  }
 
  virtual DotType evaluateDot(const NodeArg &, const StateArg &) const
  {
    mooseError("Time derivative at node not implemented for functor " + functorName());
  }
 
  virtual GradientType evaluateGradDot(const ElemArg &, const StateArg &) const
  {
    mooseError("Element gradient-dot not implemented for functor " + functorName());
  }
 
  virtual GradientType evaluateGradDot(const FaceArg &, const StateArg &) const
  {
    mooseError("Face gradient-dot not implemented for functor " + functorName());
  }
 
  virtual GradientType evaluateGradDot(const ElemQpArg &, const StateArg &) const
  {
    mooseError("Element quadrature point gradient-dot not implemented for functor " +
               functorName());
  }
 
  virtual GradientType evaluateGradDot(const ElemSideQpArg &, const StateArg &) const
  {
    mooseError("Element side quadrature point gradient-dot not implemented for functor " +
               functorName());
  }
 
  virtual GradientType evaluateGradDot(const ElemPointArg &, const StateArg &) const
  {
    mooseError("Element-point gradient-dot not implemented for functor " + functorName());
  }
 
  virtual GradientType evaluateGradDot(const NodeArg &, const StateArg &) const
  {
    mooseError("Gradient-dot at node not implemented for functor " + functorName());
  }
 
private:
  void clearCacheData();
 
  template <typename SpaceArg, typename StateArg>
  ValueType queryQpCache(unsigned int qp,
                         const libMesh::QBase & qrule,
                         std::vector<std::pair<bool, T>> & qp_cache_data,
                         const SpaceArg & space,
                         const StateArg & state) const;
 
  template <typename SpaceArg>
  ValueType queryFVArgCache(std::map<SpaceArg, ValueType> & cache_data,
                            const SpaceArg & space) const;
 
  std::set<ExecFlagType> _clearance_schedule;
 
  bool _always_evaluate;
 
  // Data for traditional element-quadrature point property evaluations which are useful for
  // caching implementation
 
  mutable dof_id_type _current_qp_map_key = libMesh::DofObject::invalid_id;
 
  mutable std::vector<std::pair<bool, ValueType>> * _current_qp_map_value = nullptr;
 
  mutable std::unordered_map<dof_id_type, std::vector<std::pair<bool, ValueType>>> _qp_to_value;
 
  // Data for traditional element-side-quadrature point property evaluations which are useful for
  // caching implementation
 
  mutable dof_id_type _current_side_qp_map_key = libMesh::DofObject::invalid_id;
 
  mutable std::vector<std::vector<std::pair<bool, ValueType>>> * _current_side_qp_map_value =
      nullptr;
 
  mutable std::unordered_map<dof_id_type, std::vector<std::vector<std::pair<bool, ValueType>>>>
      _side_qp_to_value;
 
  mutable std::map<ElemArg, ValueType> _elem_arg_to_value;
 
  mutable std::map<FaceArg, ValueType> _face_arg_to_value;
 
  mutable std::map<NodeArg, ValueType> _node_arg_to_value;
 
  MooseFunctorName _functor_name;
};
 
template <typename T>
bool
FunctorBase<T>::isInternalFace(const FaceInfo & fi) const
{
  if (!fi.neighborPtr())
    return false;
 
  return hasBlocks(fi.elem().subdomain_id()) && hasBlocks(fi.neighborPtr()->subdomain_id());
}
 
template <typename T>
template <typename SpaceArg>
typename FunctorBase<T>::ValueType
FunctorBase<T>::queryFVArgCache(std::map<SpaceArg, ValueType> & cache_data,
                                const SpaceArg & space) const
{
  // We don't want to evaluate if the key already exists, so instead we value initialize
  auto [it, inserted] = cache_data.try_emplace(space, ValueType());
  auto & value = it->second;
 
  if (inserted)
    // value not ready to go
    // this function is only called from functions that assert we are in the current time state
    value = evaluate(space, currentState());
 
  return value;
}
 
template <typename T>
typename FunctorBase<T>::ValueType
FunctorBase<T>::operator()(const ElemArg & elem, const StateArg & state) const
{
  if (_always_evaluate)
    return evaluate(elem, state);
 
  mooseAssert(state.state == 0,
              "Cached evaluations are only currently supported for the current state.");
 
  return queryFVArgCache(_elem_arg_to_value, elem);
}
 
template <typename T>
typename FunctorBase<T>::ValueType
FunctorBase<T>::operator()(const FaceArg & face_in, const StateArg & state) const
{
  checkFace(face_in);
 
  if (_always_evaluate)
    return evaluate(face_in, state);
 
  mooseAssert(state.state == 0,
              "Cached evaluations are only currently supported for the current state.");
 
  return queryFVArgCache(_face_arg_to_value, face_in);
}
 
template <typename T>
template <typename SpaceArg, typename StateArg>
typename FunctorBase<T>::ValueType
FunctorBase<T>::queryQpCache(const unsigned int qp,
                             const libMesh::QBase & qrule,
                             std::vector<std::pair<bool, ValueType>> & qp_cache_data,
                             const SpaceArg & space,
                             const StateArg & state) const
{
  // Check and see whether we even have sized for this quadrature point. If we haven't then we
  // must evaluate
  if (qp >= qp_cache_data.size())
  {
    qp_cache_data.resize(qrule.n_points(), std::make_pair(false, ValueType()));
    auto & pr = qp_cache_data[qp];
    pr.second = evaluate(space, state);
    pr.first = true;
    return pr.second;
  }
 
  // We've already sized for this qp, so let's see whether we have a valid cache value
  auto & pr = qp_cache_data[qp];
  if (pr.first)
    return pr.second;
 
  // No valid cache value so evaluate
  pr.second = evaluate(space, state);
  pr.first = true;
  return pr.second;
}
 
template <typename T>
typename FunctorBase<T>::ValueType
FunctorBase<T>::operator()(const ElemQpArg & elem_qp, const StateArg & state) const
{
  if (_always_evaluate)
    return evaluate(elem_qp, state);
 
  const auto elem_id = elem_qp.elem->id();
  if (elem_id != _current_qp_map_key)
  {
    _current_qp_map_key = elem_id;
    _current_qp_map_value = &_qp_to_value[elem_id];
  }
  auto & qp_data = *_current_qp_map_value;
  const auto qp = elem_qp.qp;
  const auto * const qrule = elem_qp.qrule;
  mooseAssert(qrule, "qrule must be non-null");
 
  return queryQpCache(qp, *qrule, qp_data, elem_qp, state);
}
 
template <typename T>
typename FunctorBase<T>::ValueType
FunctorBase<T>::operator()(const ElemSideQpArg & elem_side_qp, const StateArg & state) const
{
  if (_always_evaluate)
    return evaluate(elem_side_qp, state);
 
  const Elem * const elem = elem_side_qp.elem;
  mooseAssert(elem, "elem must be non-null");
  const auto elem_id = elem->id();
  if (elem_id != _current_side_qp_map_key)
  {
    _current_side_qp_map_key = elem_id;
    _current_side_qp_map_value = &_side_qp_to_value[elem_id];
  }
  auto & side_qp_data = *_current_side_qp_map_value;
  const auto side = elem_side_qp.side;
  const auto qp = elem_side_qp.qp;
  const auto * const qrule = elem_side_qp.qrule;
  mooseAssert(qrule, "qrule must be non-null");
 
  // Check and see whether we even have sized for this side
  if (side >= side_qp_data.size())
    side_qp_data.resize(elem->n_sides());
 
  // Ok we were sized enough for our side
  auto & qp_data = side_qp_data[side];
  return queryQpCache(qp, *qrule, qp_data, elem_side_qp, state);
}
 
template <typename T>
typename FunctorBase<T>::ValueType
FunctorBase<T>::operator()(const ElemPointArg & elem_point, const StateArg & state) const
{
  return evaluate(elem_point, state);
}
 
template <typename T>
void
FunctorBase<T>::setCacheClearanceSchedule(const std::set<ExecFlagType> & clearance_schedule)
{
  if (clearance_schedule.count(EXEC_ALWAYS))
    _always_evaluate = true;
 
  _clearance_schedule = clearance_schedule;
}
 
template <typename T>
typename FunctorBase<T>::ValueType
FunctorBase<T>::operator()(const NodeArg & node, const StateArg & state) const
{
  mooseAssert(node.subdomain_ids, "Subdomain IDs must be supplied to the node argument");
  return evaluate(node, state);
}
 
template <typename T>
void
FunctorBase<T>::checkFace(const Moose::FaceArg &
#if DEBUG
                              face
#endif
) const
{
#if DEBUG
  const Elem * const elem = face.face_side;
  const FaceInfo * const fi = face.fi;
  mooseAssert(fi, "face info should be non-null");
  bool check_elem_def = false;
  bool check_neighbor_def = false;
  // We check if the functor is defined on both sides of the face
  if (!elem)
  {
    if (!hasFaceSide(*fi, true))
      check_neighbor_def = true;
    else if (!hasFaceSide(*fi, false))
      check_elem_def = true;
  }
  else if (elem == fi->elemPtr())
    check_elem_def = true;
  else
  {
    mooseAssert(elem == fi->neighborPtr(), "This has to match something");
    check_neighbor_def = true;
  }
 
  if (check_elem_def && !hasFaceSide(*fi, true))
  {
    std::string additional_message = "It is not defined on the neighbor side either.";
    if (hasFaceSide(*fi, false))
      additional_message = "It is however defined on the neighbor side.";
    additional_message += " Face centroid: " + Moose::stringify(fi->faceCentroid());
    mooseError(_functor_name,
               " is not defined on the element side of the face information, but a face argument "
               "producer "
               "(e.g. residual object, postprocessor, etc.) has requested evaluation there.\n",
               additional_message);
  }
  if (check_neighbor_def && !hasFaceSide(*fi, false))
  {
    std::string additional_message = "It is not defined on the element side either.";
    if (hasFaceSide(*fi, true))
      additional_message = "It is however defined on the element side.";
    additional_message += " Face centroid: " + Moose::stringify(fi->faceCentroid());
    mooseError(
        _functor_name,
        " is not defined on the neighbor side of the face information, but a face argument "
        "producer (e.g. residual object, postprocessor, etc.) has requested evaluation there.\n",
        additional_message);
  }
#endif
}
 
template <typename T>
void
FunctorBase<T>::clearCacheData()
{
  for (auto & map_pr : _qp_to_value)
    for (auto & pr : map_pr.second)
      pr.first = false;
 
  for (auto & map_pr : _side_qp_to_value)
  {
    auto & side_vector = map_pr.second;
    for (auto & qp_vector : side_vector)
      for (auto & pr : qp_vector)
        pr.first = false;
  }
 
  _current_qp_map_key = libMesh::DofObject::invalid_id;
  _current_qp_map_value = nullptr;
  _current_side_qp_map_key = libMesh::DofObject::invalid_id;
  _current_side_qp_map_value = nullptr;
 
  _elem_arg_to_value.clear();
  _face_arg_to_value.clear();
  _node_arg_to_value.clear();
}
 
template <typename T>
void
FunctorBase<T>::timestepSetup()
{
  if (_clearance_schedule.count(EXEC_TIMESTEP_BEGIN))
    clearCacheData();
}
 
template <typename T>
void
FunctorBase<T>::residualSetup()
{
  if (_clearance_schedule.count(EXEC_LINEAR))
    clearCacheData();
}
 
template <typename T>
void
FunctorBase<T>::jacobianSetup()
{
  if (_clearance_schedule.count(EXEC_NONLINEAR))
    clearCacheData();
}
 
template <typename T>
void
FunctorBase<T>::customSetup(const ExecFlagType & exec_type)
{
  if (_clearance_schedule.count(exec_type))
    clearCacheData();
}
 
template <typename T>
typename FunctorBase<T>::GradientType
FunctorBase<T>::gradient(const ElemArg & elem, const StateArg & state) const
{
  return evaluateGradient(elem, state);
}
 
template <typename T>
typename FunctorBase<T>::GradientType
FunctorBase<T>::gradient(const FaceArg & face, const StateArg & state) const
{
  checkFace(face);
  return evaluateGradient(face, state);
}
 
template <typename T>
typename FunctorBase<T>::GradientType
FunctorBase<T>::gradient(const ElemQpArg & elem_qp, const StateArg & state) const
{
  return evaluateGradient(elem_qp, state);
}
 
template <typename T>
typename FunctorBase<T>::GradientType
FunctorBase<T>::gradient(const ElemSideQpArg & elem_side_qp, const StateArg & state) const
{
  return evaluateGradient(elem_side_qp, state);
}
 
template <typename T>
typename FunctorBase<T>::GradientType
FunctorBase<T>::gradient(const ElemPointArg & elem_point, const StateArg & state) const
{
  return evaluateGradient(elem_point, state);
}
 
template <typename T>
typename FunctorBase<T>::GradientType
FunctorBase<T>::gradient(const NodeArg & node, const StateArg & state) const
{
  return evaluateGradient(node, state);
}
 
template <typename T>
typename FunctorBase<T>::DotType
FunctorBase<T>::dot(const ElemArg & elem, const StateArg & state) const
{
  return evaluateDot(elem, state);
}
 
template <typename T>
typename FunctorBase<T>::DotType
FunctorBase<T>::dot(const FaceArg & face, const StateArg & state) const
{
  checkFace(face);
  return evaluateDot(face, state);
}
 
template <typename T>
typename FunctorBase<T>::DotType
FunctorBase<T>::dot(const ElemQpArg & elem_qp, const StateArg & state) const
{
  return evaluateDot(elem_qp, state);
}
 
template <typename T>
typename FunctorBase<T>::DotType
FunctorBase<T>::dot(const ElemSideQpArg & elem_side_qp, const StateArg & state) const
{
  return evaluateDot(elem_side_qp, state);
}
 
template <typename T>
typename FunctorBase<T>::DotType
FunctorBase<T>::dot(const ElemPointArg & elem_point, const StateArg & state) const
{
  return evaluateDot(elem_point, state);
}
 
template <typename T>
typename FunctorBase<T>::DotType
FunctorBase<T>::dot(const NodeArg & node, const StateArg & state) const
{
  return evaluateDot(node, state);
}
 
template <typename T>
typename FunctorBase<T>::GradientType
FunctorBase<T>::gradDot(const ElemArg & elem, const StateArg & state) const
{
  return evaluateGradDot(elem, state);
}
 
template <typename T>
typename FunctorBase<T>::GradientType
FunctorBase<T>::gradDot(const FaceArg & face, const StateArg & state) const
{
  checkFace(face);
  return evaluateGradDot(face, state);
}
 
template <typename T>
typename FunctorBase<T>::GradientType
FunctorBase<T>::gradDot(const ElemQpArg & elem_qp, const StateArg & state) const
{
  return evaluateGradDot(elem_qp, state);
}
 
template <typename T>
typename FunctorBase<T>::GradientType
FunctorBase<T>::gradDot(const ElemSideQpArg & elem_side_qp, const StateArg & state) const
{
  return evaluateGradDot(elem_side_qp, state);
}
 
template <typename T>
typename FunctorBase<T>::GradientType
FunctorBase<T>::gradDot(const ElemPointArg & elem_point, const StateArg & state) const
{
  return evaluateGradDot(elem_point, state);
}
 
template <typename T>
typename FunctorBase<T>::GradientType
FunctorBase<T>::gradDot(const NodeArg & node, const StateArg & state) const
{
  return evaluateGradDot(node, state);
}
 
template <typename T>
bool
FunctorBase<T>::hasFaceSide(const FaceInfo & fi, const bool fi_elem_side) const
{
  if (fi_elem_side)
    return hasBlocks(fi.elem().subdomain_id());
  else
    return fi.neighborPtr() && hasBlocks(fi.neighbor().subdomain_id());
}
 
template <typename T>
template <FunctorEvaluationKind FET, typename Space, typename State>
typename FunctorReturnType<T, FET>::type
FunctorBase<T>::genericEvaluate(const Space & r, const State & state) const
{
  if constexpr (FET == FunctorEvaluationKind::Value)
    return (*this)(r, state);
  else if constexpr (FET == FunctorEvaluationKind::Gradient)
    return gradient(r, state);
  else if constexpr (FET == FunctorEvaluationKind::Dot)
    return dot(r, state);
  else
    return gradDot(r, state);
}
 
class FunctorEnvelopeBase
{
public:
  FunctorEnvelopeBase() = default;
  virtual ~FunctorEnvelopeBase() = default;
 
  virtual bool wrapsNull() const = 0;
 
  virtual std::string returnType() const = 0;
 
  virtual bool isConstant() const = 0;
 
  virtual bool ownsWrappedFunctor() const = 0;
};
 
template <typename T>
class FunctorEnvelope final : public FunctorBase<T>, public FunctorEnvelopeBase
{
public:
  using typename Moose::FunctorBase<T>::ValueType;
  using typename Moose::FunctorBase<T>::GradientType;
  using typename Moose::FunctorBase<T>::DotType;
 
  FunctorEnvelope(const FunctorBase<T> & wrapped)
    : FunctorBase<T>("wraps_" + wrapped.functorName()), FunctorEnvelopeBase(), _wrapped(&wrapped)
  {
  }
 
  FunctorEnvelope(std::unique_ptr<FunctorBase<T>> && wrapped)
    : FunctorBase<T>("wraps_" + wrapped->functorName()),
      FunctorEnvelopeBase(),
      _owned(std::move(wrapped)),
      _wrapped(_owned.get())
  {
  }
 
  FunctorEnvelope(FunctorBase<T> &&) = delete;
 
  void assign(const FunctorBase<T> & wrapped)
  {
    _owned.reset();
    _wrapped = &wrapped;
  }
 
  void assign(std::unique_ptr<FunctorBase<T>> && wrapped)
  {
    _owned = std::move(wrapped);
    _wrapped = _owned.get();
  }
 
  void assign(FunctorBase<T> &&) = delete;
 
  FunctorEnvelope(const FunctorEnvelope &) = delete;
  FunctorEnvelope(FunctorEnvelope &&) = delete;
  FunctorEnvelope & operator=(const FunctorEnvelope &) = delete;
  FunctorEnvelope & operator=(FunctorEnvelope &&) = delete;
 
  virtual ~FunctorEnvelope() = default;
 
  virtual bool wrapsNull() const override { return wrapsType<NullFunctor<T>>(); }
 
  virtual std::string returnType() const override { return libMesh::demangle(typeid(T).name()); }
 
  virtual bool ownsWrappedFunctor() const override { return _owned.get(); }
 
  template <typename T2>
  bool wrapsType() const
  {
    return dynamic_cast<const T2 *>(_wrapped);
  }
 
  virtual bool isExtrapolatedBoundaryFace(const FaceInfo & fi,
                                          const Elem * const elem,
                                          const StateArg & state) const override
  {
    return _wrapped->isExtrapolatedBoundaryFace(fi, elem, state);
  }
  virtual bool isConstant() const override { return _wrapped->isConstant(); }
  virtual bool hasBlocks(const SubdomainID id) const override { return _wrapped->hasBlocks(id); }
  virtual bool hasFaceSide(const FaceInfo & fi, const bool fi_elem_side) const override
  {
    return _wrapped->hasFaceSide(fi, fi_elem_side);
  }
 
  bool supportsFaceArg() const override final { return true; }
  bool supportsElemSideQpArg() const override final { return true; }
 
protected:
 
  virtual ValueType evaluate(const ElemArg & elem, const StateArg & state) const override
  {
    return _wrapped->operator()(elem, state);
  }
  virtual ValueType evaluate(const FaceArg & face, const StateArg & state) const override
  {
    return _wrapped->operator()(face, state);
  }
  virtual ValueType evaluate(const ElemQpArg & qp, const StateArg & state) const override
  {
    return _wrapped->operator()(qp, state);
  }
  virtual ValueType evaluate(const ElemSideQpArg & qp, const StateArg & state) const override
  {
    return _wrapped->operator()(qp, state);
  }
  virtual ValueType evaluate(const ElemPointArg & elem_point, const StateArg & state) const override
  {
    return _wrapped->operator()(elem_point, state);
  }
  virtual ValueType evaluate(const NodeArg & node, const StateArg & state) const override
  {
    return _wrapped->operator()(node, state);
  }
 
  virtual GradientType evaluateGradient(const ElemArg & elem, const StateArg & state) const override
  {
    return _wrapped->gradient(elem, state);
  }
  virtual GradientType evaluateGradient(const FaceArg & face, const StateArg & state) const override
  {
    return _wrapped->gradient(face, state);
  }
  virtual GradientType evaluateGradient(const ElemQpArg & qp, const StateArg & state) const override
  {
    return _wrapped->gradient(qp, state);
  }
  virtual GradientType evaluateGradient(const ElemSideQpArg & qp,
                                        const StateArg & state) const override
  {
    return _wrapped->gradient(qp, state);
  }
  virtual GradientType evaluateGradient(const ElemPointArg & elem_point,
                                        const StateArg & state) const override
  {
    return _wrapped->gradient(elem_point, state);
  }
  virtual GradientType evaluateGradient(const NodeArg & node, const StateArg & state) const override
  {
    return _wrapped->gradient(node, state);
  }
 
  virtual DotType evaluateDot(const ElemArg & elem, const StateArg & state) const override
  {
    return _wrapped->dot(elem, state);
  }
  virtual DotType evaluateDot(const FaceArg & face, const StateArg & state) const override
  {
    return _wrapped->dot(face, state);
  }
  virtual DotType evaluateDot(const ElemQpArg & qp, const StateArg & state) const override
  {
    return _wrapped->dot(qp, state);
  }
  virtual DotType evaluateDot(const ElemSideQpArg & qp, const StateArg & state) const override
  {
    return _wrapped->dot(qp, state);
  }
  virtual DotType evaluateDot(const ElemPointArg & elem_point,
                              const StateArg & state) const override
  {
    return _wrapped->dot(elem_point, state);
  }
  virtual DotType evaluateDot(const NodeArg & node, const StateArg & state) const override
  {
    return _wrapped->dot(node, state);
  }
 
  virtual GradientType evaluateGradDot(const ElemArg & elem, const StateArg & state) const override
  {
    return _wrapped->gradDot(elem, state);
  }
  virtual GradientType evaluateGradDot(const FaceArg & face, const StateArg & state) const override
  {
    return _wrapped->gradDot(face, state);
  }
  virtual GradientType evaluateGradDot(const ElemQpArg & qp, const StateArg & state) const override
  {
    return _wrapped->gradDot(qp, state);
  }
  virtual GradientType evaluateGradDot(const ElemSideQpArg & qp,
                                       const StateArg & state) const override
  {
    return _wrapped->gradDot(qp, state);
  }
  virtual GradientType evaluateGradDot(const ElemPointArg & elem_point,
                                       const StateArg & state) const override
  {
    return _wrapped->gradDot(elem_point, state);
  }
  virtual GradientType evaluateGradDot(const NodeArg & node, const StateArg & state) const override
  {
    return _wrapped->gradDot(node, state);
  }
 
private:
  std::unique_ptr<FunctorBase<T>> _owned;
  const FunctorBase<T> * _wrapped;
 
  friend class ::SubProblem;
};
 
template <typename T>
class ConstantFunctor final : public FunctorBase<T>
{
public:
  using typename FunctorBase<T>::FunctorType;
  using typename FunctorBase<T>::ValueType;
  using typename FunctorBase<T>::GradientType;
  using typename FunctorBase<T>::DotType;
 
  ConstantFunctor(const ValueType & value)
    : FunctorBase<T>("constant_" + std::to_string(value)), _value(value)
  {
  }
  ConstantFunctor(ValueType && value)
    : FunctorBase<T>("constant_" + std::to_string(MetaPhysicL::raw_value(value))), _value(value)
  {
  }
 
  virtual bool isConstant() const override { return true; }
 
  bool hasBlocks(SubdomainID /* id */) const override { return true; }
 
  bool supportsFaceArg() const override final { return true; }
  bool supportsElemSideQpArg() const override final { return true; }
 
private:
  ValueType evaluate(const ElemArg &, const StateArg &) const override { return _value; }
  ValueType evaluate(const FaceArg &, const StateArg &) const override { return _value; }
  ValueType evaluate(const ElemQpArg &, const StateArg &) const override { return _value; }
  ValueType evaluate(const ElemSideQpArg &, const StateArg &) const override { return _value; }
  ValueType evaluate(const ElemPointArg &, const StateArg &) const override { return _value; }
  ValueType evaluate(const NodeArg &, const StateArg &) const override { return _value; }
 
  GradientType evaluateGradient(const ElemArg &, const StateArg &) const override { return 0; }
  GradientType evaluateGradient(const FaceArg &, const StateArg &) const override { return 0; }
  GradientType evaluateGradient(const ElemQpArg &, const StateArg &) const override { return 0; }
  GradientType evaluateGradient(const ElemSideQpArg &, const StateArg &) const override
  {
    return 0;
  }
  GradientType evaluateGradient(const ElemPointArg &, const StateArg &) const override { return 0; }
  GradientType evaluateGradient(const NodeArg &, const StateArg &) const override { return 0; }
 
  DotType evaluateDot(const ElemArg &, const StateArg &) const override { return 0; }
  DotType evaluateDot(const FaceArg &, const StateArg &) const override { return 0; }
  DotType evaluateDot(const ElemQpArg &, const StateArg &) const override { return 0; }
  DotType evaluateDot(const ElemSideQpArg &, const StateArg &) const override { return 0; }
  DotType evaluateDot(const ElemPointArg &, const StateArg &) const override { return 0; }
  DotType evaluateDot(const NodeArg &, const StateArg &) const override { return 0; }
 
  GradientType evaluateGradDot(const ElemArg &, const StateArg &) const override { return 0; }
  GradientType evaluateGradDot(const FaceArg &, const StateArg &) const override { return 0; }
  GradientType evaluateGradDot(const ElemQpArg &, const StateArg &) const override { return 0; }
  GradientType evaluateGradDot(const ElemSideQpArg &, const StateArg &) const override { return 0; }
  GradientType evaluateGradDot(const ElemPointArg &, const StateArg &) const override { return 0; }
  GradientType evaluateGradDot(const NodeArg &, const StateArg &) const override { return 0; }
 
private:
  ValueType _value;
};
 
template <typename T>
class NullFunctor final : public FunctorBase<T>
{
public:
  using typename FunctorBase<T>::FunctorType;
  using typename FunctorBase<T>::ValueType;
  using typename FunctorBase<T>::GradientType;
  using typename FunctorBase<T>::DotType;
 
  NullFunctor() : FunctorBase<T>("null") {}
 
  // For backwards compatiblity of unit testing
  bool hasFaceSide(const FaceInfo & fi, bool) const override;
 
  bool supportsFaceArg() const override final { return false; }
  bool supportsElemSideQpArg() const override final { return false; }
 
private:
  ValueType evaluate(const ElemArg &, const StateArg &) const override
  {
    mooseError("We should never get here. If you have, contact a MOOSE developer and tell them "
               "they've written broken code");
  }
  ValueType evaluate(const FaceArg &, const StateArg &) const override
  {
    mooseError("We should never get here. If you have, contact a MOOSE developer and tell them "
               "they've written broken code");
  }
  ValueType evaluate(const ElemQpArg &, const StateArg &) const override
  {
    mooseError("We should never get here. If you have, contact a MOOSE developer and tell them "
               "they've written broken code");
  }
  ValueType evaluate(const ElemSideQpArg &, const StateArg &) const override
  {
    mooseError("We should never get here. If you have, contact a MOOSE developer and tell them "
               "they've written broken code");
  }
  ValueType evaluate(const ElemPointArg &, const StateArg &) const override
  {
    mooseError("We should never get here. If you have, contact a MOOSE developer and tell them "
               "they've written broken code");
  }
  ValueType evaluate(const NodeArg &, const StateArg &) const override
  {
    mooseError("We should never get here. If you have, contact a MOOSE developer and tell them "
               "they've written broken code");
  }
};
 
template <typename T>
bool
NullFunctor<T>::hasFaceSide(const FaceInfo &, const bool) const
{
  // For backwards compatiblity of unit testing
  return true;
}
}