#include <language/utils/EmbeddedDiscreteFunctionMathFunctions.hpp>

#include <language/utils/EmbeddedDiscreteFunctionUtils.hpp>
#include <mesh/MeshVariant.hpp>
#include <scheme/DiscreteFunctionP0.hpp>
#include <scheme/DiscreteFunctionP0Vector.hpp>
#include <scheme/DiscreteFunctionUtils.hpp>
#include <scheme/DiscreteFunctionVariant.hpp>
#include <scheme/IDiscreteFunctionDescriptor.hpp>

#include <utils/Demangle.hpp>

#define DISCRETE_VH_TO_VH_REAL_FUNCTION_CALL(FUNCTION, ARG)                                \
  return std::visit(                                                                       \
    [&](auto&& discrete_function) -> std::shared_ptr<DiscreteFunctionVariant> {            \
      using DiscreteFunctionT = std::decay_t<decltype(discrete_function)>;                 \
      if constexpr (std::is_same_v<DiscreteFunctionT, DiscreteFunctionP0<const double>>) { \
        return std::make_shared<DiscreteFunctionVariant>(FUNCTION(discrete_function));     \
      } else {                                                                             \
        throw NormalError(EmbeddedDiscreteFunctionUtils::invalidOperandType(ARG));         \
      }                                                                                    \
    },                                                                                     \
    ARG->discreteFunction());

#define DISCRETE_VH_TO_R_CALL(FUNCTION, ARG)                                               \
  return std::visit(                                                                       \
    [&](auto&& discrete_function) -> double {                                              \
      using DiscreteFunctionT = std::decay_t<decltype(discrete_function)>;                 \
      if constexpr (std::is_same_v<DiscreteFunctionT, DiscreteFunctionP0<const double>>) { \
        return FUNCTION(discrete_function);                                                \
      } else {                                                                             \
        throw NormalError(EmbeddedDiscreteFunctionUtils::invalidOperandType(ARG));         \
      }                                                                                    \
    },                                                                                     \
    ARG->discreteFunction());

#define DISCRETE_VH_VH_TO_VH_REAL_FUNCTION_CALL(FUNCTION, ARG0, ARG1)                               \
  if (not hasSameMesh({ARG0, ARG1})) {                                                              \
    throw NormalError("operands are defined on different meshes");                                  \
  }                                                                                                 \
  return std::visit(                                                                                \
    [&](auto&& arg_f, auto&& arg_g) -> std::shared_ptr<DiscreteFunctionVariant> {                   \
      using TypeOfF = std::decay_t<decltype(arg_f)>;                                                \
      using TypeOfG = std::decay_t<decltype(arg_g)>;                                                \
      if constexpr (std::is_same_v<TypeOfF, DiscreteFunctionP0<const double>>) {                    \
        if constexpr (std::is_same_v<TypeOfF, TypeOfG>) {                                           \
          return std::make_shared<DiscreteFunctionVariant>(FUNCTION(arg_f, arg_g));                 \
        } else {                                                                                    \
          throw NormalError(EmbeddedDiscreteFunctionUtils::incompatibleOperandTypes(arg_f, arg_g)); \
        }                                                                                           \
      } else {                                                                                      \
        throw NormalError(EmbeddedDiscreteFunctionUtils::incompatibleOperandTypes(arg_f, arg_g));   \
      }                                                                                             \
    },                                                                                              \
    ARG0->discreteFunction(), ARG1->discreteFunction());

#define DISCRETE_R_VH_TO_VH_REAL_FUNCTION_CALL(FUNCTION, ARG0, ARG1)                                         \
  return std::visit(                                                                                         \
    [&](auto&& discrete_function) -> std::shared_ptr<DiscreteFunctionVariant> {                              \
      using DiscreteFunctionT = std::decay_t<decltype(discrete_function)>;                                   \
      if constexpr (std::is_same_v<DiscreteFunctionT, DiscreteFunctionP0<const double>>) {                   \
        return std::make_shared<DiscreteFunctionVariant>(FUNCTION(ARG0, discrete_function));                 \
      } else {                                                                                               \
        throw NormalError(EmbeddedDiscreteFunctionUtils::incompatibleOperandTypes(ARG0, discrete_function)); \
      }                                                                                                      \
    },                                                                                                       \
    ARG1->discreteFunction());

#define DISCRETE_VH_R_TO_VH_REAL_FUNCTION_CALL(FUNCTION, ARG0, ARG1)                                         \
  return std::visit(                                                                                         \
    [&](auto&& discrete_function) -> std::shared_ptr<DiscreteFunctionVariant> {                              \
      using DiscreteFunctionT = std::decay_t<decltype(discrete_function)>;                                   \
      if constexpr (std::is_same_v<DiscreteFunctionT, DiscreteFunctionP0<const double>>) {                   \
        return std::make_shared<DiscreteFunctionVariant>(FUNCTION(discrete_function, ARG1));                 \
      } else {                                                                                               \
        throw NormalError(EmbeddedDiscreteFunctionUtils::incompatibleOperandTypes(discrete_function, ARG1)); \
      }                                                                                                      \
    },                                                                                                       \
    ARG0->discreteFunction());

std::shared_ptr<const DiscreteFunctionVariant>
sqrt(const std::shared_ptr<const DiscreteFunctionVariant>& f)
{
  DISCRETE_VH_TO_VH_REAL_FUNCTION_CALL(sqrt, f);
}

std::shared_ptr<const DiscreteFunctionVariant>
abs(const std::shared_ptr<const DiscreteFunctionVariant>& f)
{
  DISCRETE_VH_TO_VH_REAL_FUNCTION_CALL(abs, f);
}

std::shared_ptr<const DiscreteFunctionVariant>
sin(const std::shared_ptr<const DiscreteFunctionVariant>& f)
{
  DISCRETE_VH_TO_VH_REAL_FUNCTION_CALL(sin, f);
}

std::shared_ptr<const DiscreteFunctionVariant>
cos(const std::shared_ptr<const DiscreteFunctionVariant>& f)
{
  DISCRETE_VH_TO_VH_REAL_FUNCTION_CALL(cos, f);
}

std::shared_ptr<const DiscreteFunctionVariant>
tan(const std::shared_ptr<const DiscreteFunctionVariant>& f)
{
  DISCRETE_VH_TO_VH_REAL_FUNCTION_CALL(tan, f);
}

std::shared_ptr<const DiscreteFunctionVariant>
asin(const std::shared_ptr<const DiscreteFunctionVariant>& f)
{
  DISCRETE_VH_TO_VH_REAL_FUNCTION_CALL(asin, f);
}

std::shared_ptr<const DiscreteFunctionVariant>
acos(const std::shared_ptr<const DiscreteFunctionVariant>& f)
{
  DISCRETE_VH_TO_VH_REAL_FUNCTION_CALL(acos, f);
}

std::shared_ptr<const DiscreteFunctionVariant>
atan(const std::shared_ptr<const DiscreteFunctionVariant>& f)
{
  DISCRETE_VH_TO_VH_REAL_FUNCTION_CALL(atan, f);
}

std::shared_ptr<const DiscreteFunctionVariant>
atan2(const std::shared_ptr<const DiscreteFunctionVariant>& f_v,
      const std::shared_ptr<const DiscreteFunctionVariant>& g_v)
{
  DISCRETE_VH_VH_TO_VH_REAL_FUNCTION_CALL(atan2, f_v, g_v);
}

std::shared_ptr<const DiscreteFunctionVariant>
atan2(const double a, const std::shared_ptr<const DiscreteFunctionVariant>& f)
{
  DISCRETE_R_VH_TO_VH_REAL_FUNCTION_CALL(atan2, a, f);
}

std::shared_ptr<const DiscreteFunctionVariant>
atan2(const std::shared_ptr<const DiscreteFunctionVariant>& f, const double a)
{
  DISCRETE_VH_R_TO_VH_REAL_FUNCTION_CALL(atan2, f, a);
}

std::shared_ptr<const DiscreteFunctionVariant>
sinh(const std::shared_ptr<const DiscreteFunctionVariant>& f)
{
  DISCRETE_VH_TO_VH_REAL_FUNCTION_CALL(sinh, f);
}

std::shared_ptr<const DiscreteFunctionVariant>
cosh(const std::shared_ptr<const DiscreteFunctionVariant>& f)
{
  DISCRETE_VH_TO_VH_REAL_FUNCTION_CALL(cosh, f);
}

std::shared_ptr<const DiscreteFunctionVariant>
tanh(const std::shared_ptr<const DiscreteFunctionVariant>& f)
{
  DISCRETE_VH_TO_VH_REAL_FUNCTION_CALL(tanh, f);
}

std::shared_ptr<const DiscreteFunctionVariant>
asinh(const std::shared_ptr<const DiscreteFunctionVariant>& f)
{
  DISCRETE_VH_TO_VH_REAL_FUNCTION_CALL(asinh, f);
}

std::shared_ptr<const DiscreteFunctionVariant>
acosh(const std::shared_ptr<const DiscreteFunctionVariant>& f)
{
  DISCRETE_VH_TO_VH_REAL_FUNCTION_CALL(acosh, f);
}

std::shared_ptr<const DiscreteFunctionVariant>
atanh(const std::shared_ptr<const DiscreteFunctionVariant>& f)
{
  DISCRETE_VH_TO_VH_REAL_FUNCTION_CALL(atanh, f);
}

std::shared_ptr<const DiscreteFunctionVariant>
exp(const std::shared_ptr<const DiscreteFunctionVariant>& f)
{
  DISCRETE_VH_TO_VH_REAL_FUNCTION_CALL(exp, f);
}

std::shared_ptr<const DiscreteFunctionVariant>
log(const std::shared_ptr<const DiscreteFunctionVariant>& f)
{
  DISCRETE_VH_TO_VH_REAL_FUNCTION_CALL(log, f);
}

std::shared_ptr<const DiscreteFunctionVariant>
pow(const std::shared_ptr<const DiscreteFunctionVariant>& f, const std::shared_ptr<const DiscreteFunctionVariant>& g)
{
  DISCRETE_VH_VH_TO_VH_REAL_FUNCTION_CALL(pow, f, g);
}

std::shared_ptr<const DiscreteFunctionVariant>
pow(const double a, const std::shared_ptr<const DiscreteFunctionVariant>& f)
{
  DISCRETE_R_VH_TO_VH_REAL_FUNCTION_CALL(pow, a, f);
}

std::shared_ptr<const DiscreteFunctionVariant>
pow(const std::shared_ptr<const DiscreteFunctionVariant>& f, const double a)
{
  DISCRETE_VH_R_TO_VH_REAL_FUNCTION_CALL(pow, f, a);
}

std::shared_ptr<const DiscreteFunctionVariant>
dot(const std::shared_ptr<const DiscreteFunctionVariant>& f_v,
    const std::shared_ptr<const DiscreteFunctionVariant>& g_v)
{
  if (not hasSameMesh({f_v, g_v})) {
    throw NormalError("operands are defined on different meshes");
  }

  return std::visit(
    [&](auto&& f, auto&& g) -> std::shared_ptr<DiscreteFunctionVariant> {
      using TypeOfF = std::decay_t<decltype(f)>;
      using TypeOfG = std::decay_t<decltype(g)>;
      if constexpr (not std::is_same_v<TypeOfF, TypeOfG>) {
        throw NormalError(EmbeddedDiscreteFunctionUtils::incompatibleOperandTypes(f, g));
      } else {
        using DataType = std::decay_t<typename TypeOfF::data_type>;
        if constexpr (is_discrete_function_P0_v<TypeOfF>) {
          if constexpr (is_tiny_vector_v<DataType>) {
            return std::make_shared<DiscreteFunctionVariant>(dot(f, g));
          } else {
            throw NormalError(EmbeddedDiscreteFunctionUtils::invalidOperandType(f));
          }
        } else if constexpr (is_discrete_function_P0_vector_v<TypeOfF>) {
          if (f.size() == g.size()) {
            return std::make_shared<DiscreteFunctionVariant>(dot(f, g));
          } else {
            throw NormalError("operands have different dimension");
          }
        } else {
          throw NormalError(EmbeddedDiscreteFunctionUtils::invalidOperandType(f));
        }
      }
    },
    f_v->discreteFunction(), g_v->discreteFunction());
}

template <size_t VectorDimension>
std::shared_ptr<const DiscreteFunctionVariant>
dot(const std::shared_ptr<const DiscreteFunctionVariant>& f, const TinyVector<VectorDimension>& a)
{
  return std::visit(
    [&](auto&& discrete_function0) -> std::shared_ptr<DiscreteFunctionVariant> {
      using DiscreteFunction0Type = std::decay_t<decltype(discrete_function0)>;
      if constexpr (is_discrete_function_P0_v<DiscreteFunction0Type>) {
        using DataType = std::decay_t<typename DiscreteFunction0Type::data_type>;
        if constexpr (is_tiny_vector_v<DataType>) {
          if constexpr (std::is_same_v<DataType, TinyVector<VectorDimension>>) {
            return std::make_shared<DiscreteFunctionVariant>(dot(discrete_function0, a));
          } else {
            throw NormalError(EmbeddedDiscreteFunctionUtils::incompatibleOperandTypes(f, a));
          }
        } else {
          throw NormalError(EmbeddedDiscreteFunctionUtils::invalidOperandType(f));
        }
      } else {
        throw NormalError(EmbeddedDiscreteFunctionUtils::incompatibleOperandTypes(f, a));
      }
    },
    f->discreteFunction());
}

template <size_t VectorDimension>
std::shared_ptr<const DiscreteFunctionVariant>
dot(const TinyVector<VectorDimension>& a, const std::shared_ptr<const DiscreteFunctionVariant>& f)
{
  return std::visit(
    [&](auto&& discrete_function0) -> std::shared_ptr<DiscreteFunctionVariant> {
      using DiscreteFunction0Type = std::decay_t<decltype(discrete_function0)>;
      if constexpr (is_discrete_function_P0_v<DiscreteFunction0Type>) {
        using DataType = std::decay_t<typename DiscreteFunction0Type::data_type>;
        if constexpr (is_tiny_vector_v<DataType>) {
          if constexpr (std::is_same_v<DataType, TinyVector<VectorDimension>>) {
            return std::make_shared<DiscreteFunctionVariant>(dot(a, discrete_function0));
          } else {
            throw NormalError(EmbeddedDiscreteFunctionUtils::incompatibleOperandTypes(a, f));
          }
        } else {
          throw NormalError(EmbeddedDiscreteFunctionUtils::invalidOperandType(f));
        }
      } else {
        throw NormalError(EmbeddedDiscreteFunctionUtils::incompatibleOperandTypes(a, f));
      }
    },
    f->discreteFunction());
}

template std::shared_ptr<const DiscreteFunctionVariant> dot(const std::shared_ptr<const DiscreteFunctionVariant>&,
                                                            const TinyVector<1>&);

template std::shared_ptr<const DiscreteFunctionVariant> dot(const std::shared_ptr<const DiscreteFunctionVariant>&,
                                                            const TinyVector<2>&);

template std::shared_ptr<const DiscreteFunctionVariant> dot(const std::shared_ptr<const DiscreteFunctionVariant>&,
                                                            const TinyVector<3>&);

template std::shared_ptr<const DiscreteFunctionVariant> dot(const TinyVector<1>&,
                                                            const std::shared_ptr<const DiscreteFunctionVariant>&);

template std::shared_ptr<const DiscreteFunctionVariant> dot(const TinyVector<2>&,
                                                            const std::shared_ptr<const DiscreteFunctionVariant>&);

template std::shared_ptr<const DiscreteFunctionVariant> dot(const TinyVector<3>&,
                                                            const std::shared_ptr<const DiscreteFunctionVariant>&);

std::shared_ptr<const DiscreteFunctionVariant>
doubleDot(const std::shared_ptr<const DiscreteFunctionVariant>& f_v,
          const std::shared_ptr<const DiscreteFunctionVariant>& g_v)
{
  if (not hasSameMesh({f_v, g_v})) {
    throw NormalError("operands are defined on different meshes");
  }

  return std::visit(
    [&](auto&& f, auto&& g) -> std::shared_ptr<DiscreteFunctionVariant> {
      using TypeOfF = std::decay_t<decltype(f)>;
      using TypeOfG = std::decay_t<decltype(g)>;
      if constexpr (not std::is_same_v<TypeOfF, TypeOfG>) {
        throw NormalError(EmbeddedDiscreteFunctionUtils::incompatibleOperandTypes(f, g));
      } else {
        using DataType = std::decay_t<typename TypeOfF::data_type>;
        if constexpr (is_discrete_function_P0_v<TypeOfF>) {
          if constexpr (is_tiny_matrix_v<DataType>) {
            return std::make_shared<DiscreteFunctionVariant>(doubleDot(f, g));
          } else {
            throw NormalError(EmbeddedDiscreteFunctionUtils::invalidOperandType(f));
          }
        } else {
          throw NormalError(EmbeddedDiscreteFunctionUtils::invalidOperandType(f));
        }
      }
    },
    f_v->discreteFunction(), g_v->discreteFunction());
}

template <size_t MatrixDimension>
std::shared_ptr<const DiscreteFunctionVariant>
doubleDot(const std::shared_ptr<const DiscreteFunctionVariant>& f, const TinyMatrix<MatrixDimension>& a)
{
  return std::visit(
    [&](auto&& discrete_function0) -> std::shared_ptr<DiscreteFunctionVariant> {
      using DiscreteFunction0Type = std::decay_t<decltype(discrete_function0)>;
      if constexpr (is_discrete_function_P0_v<DiscreteFunction0Type>) {
        using DataType = std::decay_t<typename DiscreteFunction0Type::data_type>;
        if constexpr (is_tiny_matrix_v<DataType>) {
          if constexpr (std::is_same_v<DataType, TinyMatrix<MatrixDimension>>) {
            return std::make_shared<DiscreteFunctionVariant>(doubleDot(discrete_function0, a));
          } else {
            throw NormalError(EmbeddedDiscreteFunctionUtils::incompatibleOperandTypes(f, a));
          }
        } else {
          throw NormalError(EmbeddedDiscreteFunctionUtils::invalidOperandType(f));
        }
      } else {
        throw NormalError(EmbeddedDiscreteFunctionUtils::incompatibleOperandTypes(f, a));
      }
    },
    f->discreteFunction());
}

template <size_t MatrixDimension>
std::shared_ptr<const DiscreteFunctionVariant>
doubleDot(const TinyMatrix<MatrixDimension>& a, const std::shared_ptr<const DiscreteFunctionVariant>& f)
{
  return std::visit(
    [&](auto&& discrete_function0) -> std::shared_ptr<DiscreteFunctionVariant> {
      using DiscreteFunction0Type = std::decay_t<decltype(discrete_function0)>;
      if constexpr (is_discrete_function_P0_v<DiscreteFunction0Type>) {
        using DataType = std::decay_t<typename DiscreteFunction0Type::data_type>;
        if constexpr (is_tiny_matrix_v<DataType>) {
          if constexpr (std::is_same_v<DataType, TinyMatrix<MatrixDimension>>) {
            return std::make_shared<DiscreteFunctionVariant>(doubleDot(a, discrete_function0));
          } else {
            throw NormalError(EmbeddedDiscreteFunctionUtils::incompatibleOperandTypes(a, f));
          }
        } else {
          throw NormalError(EmbeddedDiscreteFunctionUtils::invalidOperandType(f));
        }
      } else {
        throw NormalError(EmbeddedDiscreteFunctionUtils::incompatibleOperandTypes(a, f));
      }
    },
    f->discreteFunction());
}

template std::shared_ptr<const DiscreteFunctionVariant> doubleDot(const std::shared_ptr<const DiscreteFunctionVariant>&,
                                                                  const TinyMatrix<1>&);

template std::shared_ptr<const DiscreteFunctionVariant> doubleDot(const std::shared_ptr<const DiscreteFunctionVariant>&,
                                                                  const TinyMatrix<2>&);

template std::shared_ptr<const DiscreteFunctionVariant> doubleDot(const std::shared_ptr<const DiscreteFunctionVariant>&,
                                                                  const TinyMatrix<3>&);

template std::shared_ptr<const DiscreteFunctionVariant> doubleDot(
  const TinyMatrix<1>&,
  const std::shared_ptr<const DiscreteFunctionVariant>&);

template std::shared_ptr<const DiscreteFunctionVariant> doubleDot(
  const TinyMatrix<2>&,
  const std::shared_ptr<const DiscreteFunctionVariant>&);

template std::shared_ptr<const DiscreteFunctionVariant> doubleDot(
  const TinyMatrix<3>&,
  const std::shared_ptr<const DiscreteFunctionVariant>&);

std::shared_ptr<const DiscreteFunctionVariant>
tensorProduct(const std::shared_ptr<const DiscreteFunctionVariant>& f_v,
              const std::shared_ptr<const DiscreteFunctionVariant>& g_v)
{
  if (not hasSameMesh({f_v, g_v})) {
    throw NormalError("operands are defined on different meshes");
  }

  return std::visit(
    [&](auto&& f, auto&& g) -> std::shared_ptr<DiscreteFunctionVariant> {
      using TypeOfF = std::decay_t<decltype(f)>;
      using TypeOfG = std::decay_t<decltype(g)>;
      if constexpr (not std::is_same_v<TypeOfF, TypeOfG>) {
        throw NormalError(EmbeddedDiscreteFunctionUtils::incompatibleOperandTypes(f, g));
      } else {
        using DataType = std::decay_t<typename TypeOfF::data_type>;
        if constexpr (is_discrete_function_P0_v<TypeOfF>) {
          if constexpr (is_tiny_vector_v<DataType>) {
            return std::make_shared<DiscreteFunctionVariant>(tensorProduct(f, g));
          } else {
            throw NormalError(EmbeddedDiscreteFunctionUtils::invalidOperandType(f));
          }
        } else {
          throw NormalError(EmbeddedDiscreteFunctionUtils::invalidOperandType(f));
        }
      }
    },
    f_v->discreteFunction(), g_v->discreteFunction());
}

template <size_t VectorDimension>
std::shared_ptr<const DiscreteFunctionVariant>
tensorProduct(const std::shared_ptr<const DiscreteFunctionVariant>& f, const TinyVector<VectorDimension>& a)
{
  return std::visit(
    [&](auto&& discrete_function0) -> std::shared_ptr<DiscreteFunctionVariant> {
      using DiscreteFunction0Type = std::decay_t<decltype(discrete_function0)>;
      if constexpr (is_discrete_function_P0_v<DiscreteFunction0Type>) {
        using DataType = std::decay_t<typename DiscreteFunction0Type::data_type>;
        if constexpr (is_tiny_vector_v<DataType>) {
          if constexpr (std::is_same_v<DataType, TinyVector<VectorDimension>>) {
            return std::make_shared<DiscreteFunctionVariant>(tensorProduct<DataType>(discrete_function0, a));
          } else {
            throw NormalError(EmbeddedDiscreteFunctionUtils::incompatibleOperandTypes(f, a));
          }
        } else {
          throw NormalError(EmbeddedDiscreteFunctionUtils::invalidOperandType(f));
        }
      } else {
        throw NormalError(EmbeddedDiscreteFunctionUtils::incompatibleOperandTypes(f, a));
      }
    },
    f->discreteFunction());
}

template <size_t VectorDimension>
std::shared_ptr<const DiscreteFunctionVariant>
tensorProduct(const TinyVector<VectorDimension>& a, const std::shared_ptr<const DiscreteFunctionVariant>& f)
{
  return std::visit(
    [&](auto&& discrete_function0) -> std::shared_ptr<DiscreteFunctionVariant> {
      using DiscreteFunction0Type = std::decay_t<decltype(discrete_function0)>;
      if constexpr (is_discrete_function_P0_v<DiscreteFunction0Type>) {
        using DataType = std::decay_t<typename DiscreteFunction0Type::data_type>;
        if constexpr (is_tiny_vector_v<DataType>) {
          if constexpr (std::is_same_v<DataType, TinyVector<VectorDimension>>) {
            return std::make_shared<DiscreteFunctionVariant>(tensorProduct<DataType>(a, discrete_function0));
          } else {
            throw NormalError(EmbeddedDiscreteFunctionUtils::incompatibleOperandTypes(a, f));
          }
        } else {
          throw NormalError(EmbeddedDiscreteFunctionUtils::invalidOperandType(f));
        }
      } else {
        throw NormalError(EmbeddedDiscreteFunctionUtils::incompatibleOperandTypes(a, f));
      }
    },
    f->discreteFunction());
}

template std::shared_ptr<const DiscreteFunctionVariant> tensorProduct(
  const std::shared_ptr<const DiscreteFunctionVariant>&,
  const TinyVector<1>&);

template std::shared_ptr<const DiscreteFunctionVariant> tensorProduct(
  const std::shared_ptr<const DiscreteFunctionVariant>&,
  const TinyVector<2>&);

template std::shared_ptr<const DiscreteFunctionVariant> tensorProduct(
  const std::shared_ptr<const DiscreteFunctionVariant>&,
  const TinyVector<3>&);

template std::shared_ptr<const DiscreteFunctionVariant> tensorProduct(
  const TinyVector<1>&,
  const std::shared_ptr<const DiscreteFunctionVariant>&);

template std::shared_ptr<const DiscreteFunctionVariant> tensorProduct(
  const TinyVector<2>&,
  const std::shared_ptr<const DiscreteFunctionVariant>&);

template std::shared_ptr<const DiscreteFunctionVariant> tensorProduct(
  const TinyVector<3>&,
  const std::shared_ptr<const DiscreteFunctionVariant>&);

std::shared_ptr<const DiscreteFunctionVariant>
det(const std::shared_ptr<const DiscreteFunctionVariant>& A)
{
  return std::visit(
    [&](auto&& discrete_function) -> std::shared_ptr<DiscreteFunctionVariant> {
      using DiscreteFunctionT = std::decay_t<decltype(discrete_function)>;
      if constexpr (is_discrete_function_P0_v<DiscreteFunctionT>) {
        if constexpr (is_tiny_matrix_v<std::decay_t<typename DiscreteFunctionT::data_type>>) {
          if constexpr (DiscreteFunctionT::data_type::NumberOfRows == DiscreteFunctionT::data_type::NumberOfColumns) {
            return std::make_shared<DiscreteFunctionVariant>(det(discrete_function));
          } else {
            throw NormalError(EmbeddedDiscreteFunctionUtils::invalidOperandType(A));
          }
        } else {
          throw NormalError(EmbeddedDiscreteFunctionUtils::invalidOperandType(A));
        }
      } else {
        throw NormalError(EmbeddedDiscreteFunctionUtils::invalidOperandType(A));
      }
    },
    A->discreteFunction());
}

std::shared_ptr<const DiscreteFunctionVariant>
trace(const std::shared_ptr<const DiscreteFunctionVariant>& A)
{
  return std::visit(
    [&](auto&& discrete_function) -> std::shared_ptr<DiscreteFunctionVariant> {
      using DiscreteFunctionT = std::decay_t<decltype(discrete_function)>;
      if constexpr (is_discrete_function_P0_v<DiscreteFunctionT>) {
        if constexpr (is_tiny_matrix_v<std::decay_t<typename DiscreteFunctionT::data_type>>) {
          if constexpr (DiscreteFunctionT::data_type::NumberOfRows == DiscreteFunctionT::data_type::NumberOfColumns) {
            return std::make_shared<DiscreteFunctionVariant>(trace(discrete_function));
          } else {
            throw NormalError(EmbeddedDiscreteFunctionUtils::invalidOperandType(A));
          }
        } else {
          throw NormalError(EmbeddedDiscreteFunctionUtils::invalidOperandType(A));
        }
      } else {
        throw NormalError(EmbeddedDiscreteFunctionUtils::invalidOperandType(A));
      }
    },
    A->discreteFunction());
}

std::shared_ptr<const DiscreteFunctionVariant>
inverse(const std::shared_ptr<const DiscreteFunctionVariant>& A)
{
  return std::visit(
    [&](auto&& discrete_function) -> std::shared_ptr<DiscreteFunctionVariant> {
      using DiscreteFunctionT = std::decay_t<decltype(discrete_function)>;
      if constexpr (is_discrete_function_P0_v<DiscreteFunctionT>) {
        if constexpr (is_tiny_matrix_v<std::decay_t<typename DiscreteFunctionT::data_type>>) {
          if constexpr (DiscreteFunctionT::data_type::NumberOfRows == DiscreteFunctionT::data_type::NumberOfColumns) {
            return std::make_shared<DiscreteFunctionVariant>(inverse(discrete_function));
          } else {
            throw NormalError(EmbeddedDiscreteFunctionUtils::invalidOperandType(A));
          }
        } else {
          throw NormalError(EmbeddedDiscreteFunctionUtils::invalidOperandType(A));
        }
      } else {
        throw NormalError(EmbeddedDiscreteFunctionUtils::invalidOperandType(A));
      }
    },
    A->discreteFunction());
}

std::shared_ptr<const DiscreteFunctionVariant>
transpose(const std::shared_ptr<const DiscreteFunctionVariant>& A)
{
  return std::visit(
    [&](auto&& discrete_function) -> std::shared_ptr<DiscreteFunctionVariant> {
      using DiscreteFunctionT = std::decay_t<decltype(discrete_function)>;
      if constexpr (is_discrete_function_P0_v<DiscreteFunctionT>) {
        if constexpr (is_tiny_matrix_v<std::decay_t<typename DiscreteFunctionT::data_type>>) {
          if constexpr (DiscreteFunctionT::data_type::NumberOfRows == DiscreteFunctionT::data_type::NumberOfColumns) {
            return std::make_shared<DiscreteFunctionVariant>(transpose(discrete_function));
          } else {
            throw NormalError(EmbeddedDiscreteFunctionUtils::invalidOperandType(A));
          }
        } else {
          throw NormalError(EmbeddedDiscreteFunctionUtils::invalidOperandType(A));
        }
      } else {
        throw NormalError(EmbeddedDiscreteFunctionUtils::invalidOperandType(A));
      }
    },
    A->discreteFunction());
}

double
min(const std::shared_ptr<const DiscreteFunctionVariant>& f)
{
  DISCRETE_VH_TO_R_CALL(min, f);
}

std::shared_ptr<const DiscreteFunctionVariant>
min(const std::shared_ptr<const DiscreteFunctionVariant>& f, const std::shared_ptr<const DiscreteFunctionVariant>& g)
{
  DISCRETE_VH_VH_TO_VH_REAL_FUNCTION_CALL(min, f, g);
}

std::shared_ptr<const DiscreteFunctionVariant>
min(const double a, const std::shared_ptr<const DiscreteFunctionVariant>& f)
{
  DISCRETE_R_VH_TO_VH_REAL_FUNCTION_CALL(min, a, f);
}

std::shared_ptr<const DiscreteFunctionVariant>
min(const std::shared_ptr<const DiscreteFunctionVariant>& f, const double a)
{
  DISCRETE_VH_R_TO_VH_REAL_FUNCTION_CALL(min, f, a);
}

double
max(const std::shared_ptr<const DiscreteFunctionVariant>& f)
{
  DISCRETE_VH_TO_R_CALL(max, f);
}

std::shared_ptr<const DiscreteFunctionVariant>
max(const std::shared_ptr<const DiscreteFunctionVariant>& f, const std::shared_ptr<const DiscreteFunctionVariant>& g)
{
  DISCRETE_VH_VH_TO_VH_REAL_FUNCTION_CALL(max, f, g);
}

std::shared_ptr<const DiscreteFunctionVariant>
max(const double a, const std::shared_ptr<const DiscreteFunctionVariant>& f)
{
  DISCRETE_R_VH_TO_VH_REAL_FUNCTION_CALL(max, a, f);
}

std::shared_ptr<const DiscreteFunctionVariant>
max(const std::shared_ptr<const DiscreteFunctionVariant>& f, const double a)
{
  DISCRETE_VH_R_TO_VH_REAL_FUNCTION_CALL(max, f, a);
}

template <typename ValueT>
ValueT
sum_of(const std::shared_ptr<const DiscreteFunctionVariant>& f)
{
  ValueT value;
  std::visit(
    [&](auto&& discrete_function) {
      using DiscreteFunctionT = std::decay_t<decltype(discrete_function)>;
      if constexpr (is_discrete_function_P0_v<DiscreteFunctionT>) {
        using DataType = std::decay_t<typename DiscreteFunctionT::data_type>;
        if constexpr (std::is_same_v<ValueT, DataType>) {
          value = sum(discrete_function);
        } else {
          throw NormalError(EmbeddedDiscreteFunctionUtils::invalidOperandType(f));
        }
      } else {
        throw NormalError(EmbeddedDiscreteFunctionUtils::invalidOperandType(f));
      }
    },
    f->discreteFunction());

  return value;
}

std::shared_ptr<const DiscreteFunctionVariant>
sum_of_Vh_components(const std::shared_ptr<const DiscreteFunctionVariant>& f)
{
  return std::visit(
    [&](auto&& discrete_function) -> std::shared_ptr<const DiscreteFunctionVariant> {
      using DiscreteFunctionT = std::decay_t<decltype(discrete_function)>;
      if constexpr (is_discrete_function_P0_vector_v<DiscreteFunctionT>) {
        using DataType = std::decay_t<typename DiscreteFunctionT::data_type>;
        static_assert(std::is_same_v<DataType, double>);
        return std::make_shared<DiscreteFunctionVariant>(sumOfComponents(discrete_function));
      } else {
        throw NormalError(EmbeddedDiscreteFunctionUtils::invalidOperandType(f));
      }
    },
    f->discreteFunction());
}

std::shared_ptr<const DiscreteFunctionVariant>
vectorize(const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_list)
{
  if (hasSameMesh(discrete_function_list)) {
    std::shared_ptr mesh_v = getCommonMesh(discrete_function_list);
    Assert(mesh_v.use_count() > 0);

    DiscreteFunctionP0Vector<double> discrete_vector_function(mesh_v, discrete_function_list.size());

    for (size_t i_discrete_function = 0; i_discrete_function < discrete_function_list.size(); ++i_discrete_function) {
      std::visit(
        [&](auto&& discrete_function) {
          using DiscreteFunctionT = std::decay_t<decltype(discrete_function)>;
          if constexpr (is_discrete_function_P0_v<DiscreteFunctionT>) {
            using DataType = std::remove_const_t<typename DiscreteFunctionT::data_type>;
            if constexpr (std::is_same_v<DataType, double>) {
              parallel_for(
                mesh_v->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
                  discrete_vector_function[cell_id][i_discrete_function] = discrete_function[cell_id];
                });
            } else {
              throw NormalError(EmbeddedDiscreteFunctionUtils::invalidOperandType(discrete_function));
            }
          } else {
            throw NormalError(EmbeddedDiscreteFunctionUtils::invalidOperandType(discrete_function));
          }
        },
        discrete_function_list[i_discrete_function]->discreteFunction());
    }

    return std::make_shared<DiscreteFunctionVariant>(discrete_vector_function);
  } else {
    throw NormalError("discrete functions are not defined on the same mesh");
  }
}

template double sum_of<double>(const std::shared_ptr<const DiscreteFunctionVariant>&);

template TinyVector<1> sum_of<TinyVector<1>>(const std::shared_ptr<const DiscreteFunctionVariant>&);

template TinyVector<2> sum_of<TinyVector<2>>(const std::shared_ptr<const DiscreteFunctionVariant>&);

template TinyVector<3> sum_of<TinyVector<3>>(const std::shared_ptr<const DiscreteFunctionVariant>&);

template TinyMatrix<1> sum_of<TinyMatrix<1>>(const std::shared_ptr<const DiscreteFunctionVariant>&);

template TinyMatrix<2> sum_of<TinyMatrix<2>>(const std::shared_ptr<const DiscreteFunctionVariant>&);

template TinyMatrix<3> sum_of<TinyMatrix<3>>(const std::shared_ptr<const DiscreteFunctionVariant>&);

template <typename ValueT>
ValueT
integral_of(const std::shared_ptr<const DiscreteFunctionVariant>& f)
{
  return std::visit(
    [&](auto&& discrete_function) -> ValueT {
      using DiscreteFunctionT = std::decay_t<decltype(discrete_function)>;
      if constexpr (is_discrete_function_P0_v<DiscreteFunctionT>) {
        using DataType = std::decay_t<typename DiscreteFunctionT::data_type>;
        if constexpr (std::is_same_v<ValueT, DataType>) {
          return integrate(discrete_function);
        } else {
          throw NormalError(EmbeddedDiscreteFunctionUtils::invalidOperandType(f));
        }
      } else {
        throw NormalError(EmbeddedDiscreteFunctionUtils::invalidOperandType(f));
      }
    },
    f->discreteFunction());
}

template double integral_of<double>(const std::shared_ptr<const DiscreteFunctionVariant>&);

template TinyVector<1> integral_of<TinyVector<1>>(const std::shared_ptr<const DiscreteFunctionVariant>&);

template TinyVector<2> integral_of<TinyVector<2>>(const std::shared_ptr<const DiscreteFunctionVariant>&);

template TinyVector<3> integral_of<TinyVector<3>>(const std::shared_ptr<const DiscreteFunctionVariant>&);

template TinyMatrix<1> integral_of<TinyMatrix<1>>(const std::shared_ptr<const DiscreteFunctionVariant>&);

template TinyMatrix<2> integral_of<TinyMatrix<2>>(const std::shared_ptr<const DiscreteFunctionVariant>&);

template TinyMatrix<3> integral_of<TinyMatrix<3>>(const std::shared_ptr<const DiscreteFunctionVariant>&);