DiscreteFunctionP0.hpp

#ifndef DISCRETE_FUNCTION_P0_HPP
#define DISCRETE_FUNCTION_P0_HPP

#include <scheme/IDiscreteFunction.hpp>

#include <mesh/Connectivity.hpp>
#include <mesh/Mesh.hpp>
#include <mesh/MeshData.hpp>
#include <mesh/MeshDataManager.hpp>
#include <scheme/DiscreteFunctionDescriptorP0.hpp>

template <size_t Dimension, typename DataType>
class DiscreteFunctionP0 : public IDiscreteFunction
{
 public:
  using data_type = DataType;
  using MeshType  = Mesh<Connectivity<Dimension>>;

  static constexpr HandledItemDataType handled_data_type = HandledItemDataType::value;

  friend class DiscreteFunctionP0<Dimension, std::add_const_t<DataType>>;
  friend class DiscreteFunctionP0<Dimension, std::remove_const_t<DataType>>;

 private:
  std::shared_ptr<const MeshType> m_mesh;
  CellValue<DataType> m_cell_values;

  DiscreteFunctionDescriptorP0 m_discrete_function_descriptor;

 public:
  PUGS_INLINE
  ASTNodeDataType
  dataType() const final
  {
    return ast_node_data_type_from<DataType>;
  }

  PUGS_INLINE
  const CellValue<DataType>&
  cellValues() const
  {
    return m_cell_values;
  }

  PUGS_INLINE
  std::shared_ptr<const IMesh>
  mesh() const
  {
    return m_mesh;
  }

  PUGS_INLINE
  const IDiscreteFunctionDescriptor&
  descriptor() const final
  {
    return m_discrete_function_descriptor;
  }

  PUGS_FORCEINLINE
  operator DiscreteFunctionP0<Dimension, const DataType>() const
  {
    return DiscreteFunctionP0<Dimension, const DataType>(m_mesh, m_cell_values);
  }

  PUGS_INLINE
  void
  fill(const DataType& data) const noexcept
  {
    static_assert(not std::is_const_v<DataType>, "Cannot modify ItemValue of const");
    m_cell_values.fill(data);
  }

  PUGS_INLINE DiscreteFunctionP0
  operator=(const DiscreteFunctionP0& f)
  {
    Assert(m_mesh == f.m_mesh);
    m_cell_values = f.m_cell_values;
    return *this;
  }

  friend PUGS_INLINE DiscreteFunctionP0<Dimension, std::remove_const_t<DataType>>
  copy(const DiscreteFunctionP0& source)
  {
    return DiscreteFunctionP0<Dimension, std::remove_const_t<DataType>>{source.m_mesh, copy(source.cellValues())};
  }

  friend PUGS_INLINE void
  copy_to(const DiscreteFunctionP0<Dimension, DataType>& source,
          DiscreteFunctionP0<Dimension, std::remove_const_t<DataType>>& destination)
  {
    Assert(source.m_mesh == destination.m_mesh);
    copy_to(source.m_cell_values, destination.m_cell_values);
  }

  PUGS_FORCEINLINE DataType&
  operator[](const CellId cell_id) const noexcept(NO_ASSERT)
  {
    return m_cell_values[cell_id];
  }

  PUGS_INLINE DiscreteFunctionP0<Dimension, std::remove_const_t<DataType>>
  operator-() const
  {
    Assert(m_cell_values.isBuilt());
    DiscreteFunctionP0<Dimension, std::remove_const_t<DataType>> opposite = copy(*this);
    parallel_for(
      m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { opposite[cell_id] = -opposite[cell_id]; });

    return opposite;
  }

  template <typename DataType2T>
  PUGS_INLINE DiscreteFunctionP0<Dimension, decltype(DataType{} + DataType2T{})>
  operator+(const DiscreteFunctionP0<Dimension, DataType2T>& g) const
  {
    const DiscreteFunctionP0& f = *this;
    Assert(f.m_cell_values.isBuilt() and g.m_cell_values.isBuilt());
    Assert(f.mesh() == g.mesh(), "functions are not defined on the same mesh");
    DiscreteFunctionP0<Dimension, decltype(DataType{} + DataType2T{})> sum(f.m_mesh);
    parallel_for(
      f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { sum[cell_id] = f[cell_id] + g[cell_id]; });
    return sum;
  }

  template <typename LHSDataType>
  PUGS_INLINE friend DiscreteFunctionP0<Dimension, decltype(LHSDataType{} + DataType{})>
  operator+(const LHSDataType& a, const DiscreteFunctionP0& g)
  {
    using SumDataType = decltype(LHSDataType{} + DataType{});
    Assert(g.m_cell_values.isBuilt());
    DiscreteFunctionP0<Dimension, SumDataType> sum(g.m_mesh);
    parallel_for(
      g.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { sum[cell_id] = a + g[cell_id]; });
    return sum;
  }

  template <typename RHSDataType>
  PUGS_INLINE friend DiscreteFunctionP0<Dimension, decltype(DataType{} + RHSDataType{})>
  operator+(const DiscreteFunctionP0& f, const RHSDataType& b)
  {
    using SumDataType = decltype(DataType{} + RHSDataType{});
    Assert(f.m_cell_values.isBuilt());
    DiscreteFunctionP0<Dimension, SumDataType> sum(f.m_mesh);
    parallel_for(
      f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { sum[cell_id] = f[cell_id] + b; });
    return sum;
  }

  template <typename DataType2T>
  PUGS_INLINE DiscreteFunctionP0<Dimension, decltype(DataType{} - DataType2T{})>
  operator-(const DiscreteFunctionP0<Dimension, DataType2T>& g) const
  {
    const DiscreteFunctionP0& f = *this;
    Assert(f.m_cell_values.isBuilt() and g.m_cell_values.isBuilt());
    Assert(f.mesh() == g.mesh(), "functions are not defined on the same mesh");
    DiscreteFunctionP0<Dimension, decltype(DataType{} - DataType2T{})> difference(f.m_mesh);
    parallel_for(
      f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { difference[cell_id] = f[cell_id] - g[cell_id]; });
    return difference;
  }

  template <typename LHSDataType>
  PUGS_INLINE friend DiscreteFunctionP0<Dimension, decltype(LHSDataType{} - DataType{})>
  operator-(const LHSDataType& a, const DiscreteFunctionP0& g)
  {
    using DifferenceDataType = decltype(LHSDataType{} - DataType{});
    Assert(g.m_cell_values.isBuilt());
    DiscreteFunctionP0<Dimension, DifferenceDataType> difference(g.m_mesh);
    parallel_for(
      g.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { difference[cell_id] = a - g[cell_id]; });
    return difference;
  }

  template <typename RHSDataType>
  PUGS_INLINE friend DiscreteFunctionP0<Dimension, decltype(DataType{} - RHSDataType{})>
  operator-(const DiscreteFunctionP0& f, const RHSDataType& b)
  {
    using DifferenceDataType = decltype(DataType{} - RHSDataType{});
    Assert(f.m_cell_values.isBuilt());
    DiscreteFunctionP0<Dimension, DifferenceDataType> difference(f.m_mesh);
    parallel_for(
      f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { difference[cell_id] = f[cell_id] - b; });
    return difference;
  }

  template <typename DataType2T>
  PUGS_INLINE DiscreteFunctionP0<Dimension, decltype(DataType{} * DataType2T{})>
  operator*(const DiscreteFunctionP0<Dimension, DataType2T>& g) const
  {
    const DiscreteFunctionP0& f = *this;
    Assert(f.m_cell_values.isBuilt());
    Assert(f.mesh() == g.mesh(), "functions are not defined on the same mesh");
    DiscreteFunctionP0<Dimension, decltype(DataType{} * DataType2T{})> product(f.m_mesh);
    parallel_for(
      f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { product[cell_id] = f[cell_id] * g[cell_id]; });
    return product;
  }

  template <typename LHSDataType>
  PUGS_INLINE friend DiscreteFunctionP0<Dimension, decltype(LHSDataType{} * DataType{})>
  operator*(const LHSDataType& a, const DiscreteFunctionP0& f)
  {
    using ProductDataType = decltype(LHSDataType{} * DataType{});
    Assert(f.m_cell_values.isBuilt());
    DiscreteFunctionP0<Dimension, ProductDataType> product(f.m_mesh);
    parallel_for(
      f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { product[cell_id] = a * f[cell_id]; });
    return product;
  }

  template <typename RHSDataType>
  PUGS_INLINE friend DiscreteFunctionP0<Dimension, decltype(DataType{} * RHSDataType{})>
  operator*(const DiscreteFunctionP0& f, const RHSDataType& b)
  {
    using ProductDataType = decltype(DataType{} * RHSDataType{});
    Assert(f.m_cell_values.isBuilt());
    DiscreteFunctionP0<Dimension, ProductDataType> product(f.m_mesh);
    parallel_for(
      f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { product[cell_id] = f[cell_id] * b; });
    return product;
  }

  template <typename DataType2T>
  PUGS_INLINE DiscreteFunctionP0<Dimension, decltype(DataType{} / DataType2T{})>
  operator/(const DiscreteFunctionP0<Dimension, DataType2T>& g) const
  {
    const DiscreteFunctionP0& f = *this;
    Assert(f.m_cell_values.isBuilt() and g.m_cell_values.isBuilt());
    Assert(f.mesh() == g.mesh(), "functions are not defined on the same mesh");
    DiscreteFunctionP0<Dimension, decltype(DataType{} / DataType2T{})> ratio(f.m_mesh);
    parallel_for(
      f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { ratio[cell_id] = f[cell_id] / g[cell_id]; });
    return ratio;
  }

  template <typename LHSDataType>
  PUGS_INLINE friend DiscreteFunctionP0<Dimension, decltype(LHSDataType{} / DataType{})>
  operator/(const LHSDataType& a, const DiscreteFunctionP0& f)
  {
    using RatioDataType = decltype(LHSDataType{} / DataType{});
    Assert(f.m_cell_values.isBuilt());
    DiscreteFunctionP0<Dimension, RatioDataType> ratio(f.m_mesh);
    parallel_for(
      f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { ratio[cell_id] = a / f[cell_id]; });
    return ratio;
  }

  template <typename RHSDataType>
  PUGS_INLINE friend DiscreteFunctionP0<Dimension, decltype(DataType{} / RHSDataType{})>
  operator/(const DiscreteFunctionP0& f, const RHSDataType& b)
  {
    using RatioDataType = decltype(DataType{} / RHSDataType{});
    Assert(f.m_cell_values.isBuilt());
    DiscreteFunctionP0<Dimension, RatioDataType> ratio(f.m_mesh);
    parallel_for(
      f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { ratio[cell_id] = f[cell_id] / b; });
    return ratio;
  }

  PUGS_INLINE friend DiscreteFunctionP0<Dimension, double>
  sqrt(const DiscreteFunctionP0& f)
  {
    static_assert(std::is_arithmetic_v<DataType>);
    Assert(f.m_cell_values.isBuilt());
    DiscreteFunctionP0<Dimension, std::remove_const_t<DataType>> result{f.m_mesh};
    parallel_for(
      f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { result[cell_id] = std::sqrt(f[cell_id]); });

    return result;
  }

  PUGS_INLINE friend DiscreteFunctionP0<Dimension, double>
  abs(const DiscreteFunctionP0& f)
  {
    static_assert(std::is_arithmetic_v<DataType>);
    Assert(f.m_cell_values.isBuilt());
    DiscreteFunctionP0<Dimension, std::remove_const_t<DataType>> result{f.m_mesh};
    parallel_for(
      f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { result[cell_id] = std::abs(f[cell_id]); });

    return result;
  }

  PUGS_INLINE friend DiscreteFunctionP0<Dimension, double>
  sin(const DiscreteFunctionP0& f)
  {
    static_assert(std::is_arithmetic_v<DataType>);
    Assert(f.m_cell_values.isBuilt());
    DiscreteFunctionP0<Dimension, std::remove_const_t<DataType>> result{f.m_mesh};
    parallel_for(
      f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { result[cell_id] = std::sin(f[cell_id]); });

    return result;
  }

  PUGS_INLINE friend DiscreteFunctionP0<Dimension, double>
  cos(const DiscreteFunctionP0& f)
  {
    static_assert(std::is_arithmetic_v<DataType>);
    Assert(f.m_cell_values.isBuilt());
    DiscreteFunctionP0<Dimension, std::remove_const_t<DataType>> result{f.m_mesh};
    parallel_for(
      f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { result[cell_id] = std::cos(f[cell_id]); });

    return result;
  }

  PUGS_INLINE friend DiscreteFunctionP0<Dimension, double>
  tan(const DiscreteFunctionP0& f)
  {
    static_assert(std::is_arithmetic_v<DataType>);
    Assert(f.m_cell_values.isBuilt());
    DiscreteFunctionP0<Dimension, std::remove_const_t<DataType>> result{f.m_mesh};
    parallel_for(
      f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { result[cell_id] = std::tan(f[cell_id]); });

    return result;
  }

  PUGS_INLINE friend DiscreteFunctionP0<Dimension, double>
  asin(const DiscreteFunctionP0& f)
  {
    static_assert(std::is_arithmetic_v<DataType>);
    Assert(f.m_cell_values.isBuilt());
    DiscreteFunctionP0<Dimension, std::remove_const_t<DataType>> result{f.m_mesh};
    parallel_for(
      f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { result[cell_id] = std::asin(f[cell_id]); });

    return result;
  }

  PUGS_INLINE friend DiscreteFunctionP0<Dimension, double>
  acos(const DiscreteFunctionP0& f)
  {
    static_assert(std::is_arithmetic_v<DataType>);
    Assert(f.m_cell_values.isBuilt());
    DiscreteFunctionP0<Dimension, std::remove_const_t<DataType>> result{f.m_mesh};
    parallel_for(
      f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { result[cell_id] = std::acos(f[cell_id]); });

    return result;
  }

  PUGS_INLINE friend DiscreteFunctionP0<Dimension, double>
  atan(const DiscreteFunctionP0& f)
  {
    static_assert(std::is_arithmetic_v<DataType>);
    Assert(f.m_cell_values.isBuilt());
    DiscreteFunctionP0<Dimension, std::remove_const_t<DataType>> result{f.m_mesh};
    parallel_for(
      f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { result[cell_id] = std::atan(f[cell_id]); });

    return result;
  }

  PUGS_INLINE friend DiscreteFunctionP0<Dimension, double>
  sinh(const DiscreteFunctionP0& f)
  {
    static_assert(std::is_arithmetic_v<DataType>);
    Assert(f.m_cell_values.isBuilt());
    DiscreteFunctionP0<Dimension, std::remove_const_t<DataType>> result{f.m_mesh};
    parallel_for(
      f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { result[cell_id] = std::sinh(f[cell_id]); });

    return result;
  }

  PUGS_INLINE friend DiscreteFunctionP0<Dimension, double>
  cosh(const DiscreteFunctionP0& f)
  {
    static_assert(std::is_arithmetic_v<DataType>);
    Assert(f.m_cell_values.isBuilt());
    DiscreteFunctionP0<Dimension, std::remove_const_t<DataType>> result{f.m_mesh};
    parallel_for(
      f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { result[cell_id] = std::cosh(f[cell_id]); });

    return result;
  }

  PUGS_INLINE friend DiscreteFunctionP0<Dimension, double>
  tanh(const DiscreteFunctionP0& f)
  {
    static_assert(std::is_arithmetic_v<DataType>);
    Assert(f.m_cell_values.isBuilt());
    DiscreteFunctionP0<Dimension, std::remove_const_t<DataType>> result{f.m_mesh};
    parallel_for(
      f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { result[cell_id] = std::tanh(f[cell_id]); });

    return result;
  }

  PUGS_INLINE friend DiscreteFunctionP0<Dimension, double>
  asinh(const DiscreteFunctionP0& f)
  {
    static_assert(std::is_arithmetic_v<DataType>);
    Assert(f.m_cell_values.isBuilt());
    DiscreteFunctionP0<Dimension, std::remove_const_t<DataType>> result{f.m_mesh};
    parallel_for(
      f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { result[cell_id] = std::asinh(f[cell_id]); });

    return result;
  }

  PUGS_INLINE friend DiscreteFunctionP0<Dimension, double>
  acosh(const DiscreteFunctionP0& f)
  {
    static_assert(std::is_arithmetic_v<DataType>);
    Assert(f.m_cell_values.isBuilt());
    DiscreteFunctionP0<Dimension, std::remove_const_t<DataType>> result{f.m_mesh};
    parallel_for(
      f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { result[cell_id] = std::acosh(f[cell_id]); });

    return result;
  }

  PUGS_INLINE friend DiscreteFunctionP0<Dimension, double>
  atanh(const DiscreteFunctionP0& f)
  {
    static_assert(std::is_arithmetic_v<DataType>);
    Assert(f.m_cell_values.isBuilt());
    DiscreteFunctionP0<Dimension, std::remove_const_t<DataType>> result{f.m_mesh};
    parallel_for(
      f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { result[cell_id] = std::atanh(f[cell_id]); });

    return result;
  }

  PUGS_INLINE friend DiscreteFunctionP0<Dimension, double>
  exp(const DiscreteFunctionP0& f)
  {
    static_assert(std::is_arithmetic_v<DataType>);
    Assert(f.m_cell_values.isBuilt());
    DiscreteFunctionP0<Dimension, std::remove_const_t<DataType>> result{f.m_mesh};
    parallel_for(
      f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { result[cell_id] = std::exp(f[cell_id]); });

    return result;
  }

  PUGS_INLINE friend DiscreteFunctionP0<Dimension, double>
  log(const DiscreteFunctionP0& f)
  {
    static_assert(std::is_arithmetic_v<DataType>);
    Assert(f.m_cell_values.isBuilt());
    DiscreteFunctionP0<Dimension, std::remove_const_t<DataType>> result{f.m_mesh};
    parallel_for(
      f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { result[cell_id] = std::log(f[cell_id]); });

    return result;
  }

  PUGS_INLINE friend DiscreteFunctionP0
  atan2(const DiscreteFunctionP0& f, const DiscreteFunctionP0& g)
  {
    static_assert(std::is_arithmetic_v<DataType>);
    Assert(f.m_cell_values.isBuilt() and g.m_cell_values.isBuilt());
    Assert(f.m_mesh == g.m_mesh);
    DiscreteFunctionP0<Dimension, std::remove_const_t<DataType>> result{f.m_mesh};
    parallel_for(
      f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { result[cell_id] = std::atan2(f[cell_id], g[cell_id]); });

    return result;
  }

  PUGS_INLINE friend DiscreteFunctionP0
  atan2(const double a, const DiscreteFunctionP0& f)
  {
    static_assert(std::is_arithmetic_v<DataType>);
    Assert(f.m_cell_values.isBuilt());
    DiscreteFunctionP0<Dimension, std::remove_const_t<DataType>> result{f.m_mesh};
    parallel_for(
      f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { result[cell_id] = std::atan2(a, f[cell_id]); });

    return result;
  }

  PUGS_INLINE friend DiscreteFunctionP0
  atan2(const DiscreteFunctionP0& f, const double a)
  {
    static_assert(std::is_arithmetic_v<DataType>);
    Assert(f.m_cell_values.isBuilt());
    DiscreteFunctionP0<Dimension, std::remove_const_t<DataType>> result{f.m_mesh};
    parallel_for(
      f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { result[cell_id] = std::atan2(f[cell_id], a); });

    return result;
  }

  PUGS_INLINE friend DiscreteFunctionP0
  pow(const DiscreteFunctionP0& f, const DiscreteFunctionP0& g)
  {
    static_assert(std::is_arithmetic_v<DataType>);
    Assert(f.m_cell_values.isBuilt() and g.m_cell_values.isBuilt());
    Assert(f.m_mesh == g.m_mesh);
    DiscreteFunctionP0<Dimension, std::remove_const_t<DataType>> result{f.m_mesh};
    parallel_for(
      f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { result[cell_id] = std::pow(f[cell_id], g[cell_id]); });

    return result;
  }

  PUGS_INLINE friend DiscreteFunctionP0
  pow(const double a, const DiscreteFunctionP0& f)
  {
    static_assert(std::is_arithmetic_v<DataType>);
    Assert(f.m_cell_values.isBuilt());
    DiscreteFunctionP0<Dimension, std::remove_const_t<DataType>> result{f.m_mesh};
    parallel_for(
      f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { result[cell_id] = std::pow(a, f[cell_id]); });

    return result;
  }

  PUGS_INLINE friend DiscreteFunctionP0
  pow(const DiscreteFunctionP0& f, const double a)
  {
    static_assert(std::is_arithmetic_v<DataType>);
    Assert(f.m_cell_values.isBuilt());
    DiscreteFunctionP0<Dimension, std::remove_const_t<DataType>> result{f.m_mesh};
    parallel_for(
      f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { result[cell_id] = std::pow(f[cell_id], a); });

    return result;
  }

  PUGS_INLINE friend DiscreteFunctionP0<Dimension, double>
  dot(const DiscreteFunctionP0& f, const DiscreteFunctionP0& g)
  {
    static_assert(is_tiny_vector_v<DataType>);
    Assert(f.m_cell_values.isBuilt() and g.m_cell_values.isBuilt());
    DiscreteFunctionP0<Dimension, double> result{f.m_mesh};
    parallel_for(
      f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { result[cell_id] = dot(f[cell_id], g[cell_id]); });

    return result;
  }

  PUGS_INLINE friend DiscreteFunctionP0<Dimension, double>
  dot(const DiscreteFunctionP0& f, const DataType& a)
  {
    static_assert(is_tiny_vector_v<DataType>);
    Assert(f.m_cell_values.isBuilt());
    DiscreteFunctionP0<Dimension, double> result{f.m_mesh};
    parallel_for(
      f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { result[cell_id] = dot(f[cell_id], a); });

    return result;
  }

  PUGS_INLINE friend DiscreteFunctionP0<Dimension, double>
  dot(const DataType& a, const DiscreteFunctionP0& f)
  {
    static_assert(is_tiny_vector_v<DataType>);
    Assert(f.m_cell_values.isBuilt());
    DiscreteFunctionP0<Dimension, double> result{f.m_mesh};
    parallel_for(
      f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { result[cell_id] = dot(a, f[cell_id]); });

    return result;
  }

  DiscreteFunctionP0(const std::shared_ptr<const MeshType>& mesh) : m_mesh{mesh}, m_cell_values{mesh->connectivity()} {}

  DiscreteFunctionP0(const std::shared_ptr<const MeshType>& mesh, const CellValue<DataType>& cell_value)
    : m_mesh{mesh}, m_cell_values{cell_value}
  {
    Assert(mesh->shared_connectivity() == cell_value.connectivity_ptr());
  }

  DiscreteFunctionP0() noexcept = delete;

  DiscreteFunctionP0(const DiscreteFunctionP0&) noexcept = default;
  DiscreteFunctionP0(DiscreteFunctionP0&&) noexcept      = default;

  ~DiscreteFunctionP0() = default;
};

#endif   // DISCRETE_FUNCTION_P0_HPP