diff --git a/src/language/modules/MathFunctionRegisterForVh.cpp b/src/language/modules/MathFunctionRegisterForVh.cpp
index 7617463648fdee030760ef65ed9dec50bdb21232..22f7866bddd9a0539d9d2cb669c69d69c0be23a8 100644
--- a/src/language/modules/MathFunctionRegisterForVh.cpp
+++ b/src/language/modules/MathFunctionRegisterForVh.cpp
@@ -253,49 +253,97 @@ MathFunctionRegisterForVh::MathFunctionRegisterForVh(SchemeModule& scheme_module
double b) -> std::shared_ptr<const IDiscreteFunction> { return max(a, b); }));
scheme_module
- ._addBuiltinFunction("sum_to_R",
+ ._addBuiltinFunction("sum_of_R",
std::make_shared<BuiltinFunctionEmbedder<double(std::shared_ptr<const IDiscreteFunction>)>>(
- [](std::shared_ptr<const IDiscreteFunction> a) -> double { return sum_to<double>(a); }));
+ [](std::shared_ptr<const IDiscreteFunction> a) -> double { return sum_of<double>(a); }));
- scheme_module._addBuiltinFunction("sum_to_R1",
+ scheme_module._addBuiltinFunction("sum_of_R1",
std::make_shared<
BuiltinFunctionEmbedder<TinyVector<1>(std::shared_ptr<const IDiscreteFunction>)>>(
[](std::shared_ptr<const IDiscreteFunction> a) -> TinyVector<1> {
- return sum_to<TinyVector<1>>(a);
+ return sum_of<TinyVector<1>>(a);
}));
- scheme_module._addBuiltinFunction("sum_to_R2",
+ scheme_module._addBuiltinFunction("sum_of_R2",
std::make_shared<
BuiltinFunctionEmbedder<TinyVector<2>(std::shared_ptr<const IDiscreteFunction>)>>(
[](std::shared_ptr<const IDiscreteFunction> a) -> TinyVector<2> {
- return sum_to<TinyVector<2>>(a);
+ return sum_of<TinyVector<2>>(a);
}));
- scheme_module._addBuiltinFunction("sum_to_R3",
+ scheme_module._addBuiltinFunction("sum_of_R3",
std::make_shared<
BuiltinFunctionEmbedder<TinyVector<3>(std::shared_ptr<const IDiscreteFunction>)>>(
[](std::shared_ptr<const IDiscreteFunction> a) -> TinyVector<3> {
- return sum_to<TinyVector<3>>(a);
+ return sum_of<TinyVector<3>>(a);
}));
- scheme_module._addBuiltinFunction("sum_to_R1x1",
+ scheme_module._addBuiltinFunction("sum_of_R1x1",
std::make_shared<
BuiltinFunctionEmbedder<TinyMatrix<1>(std::shared_ptr<const IDiscreteFunction>)>>(
[](std::shared_ptr<const IDiscreteFunction> a) -> TinyMatrix<1> {
- return sum_to<TinyMatrix<1>>(a);
+ return sum_of<TinyMatrix<1>>(a);
}));
- scheme_module._addBuiltinFunction("sum_to_R2x2",
+ scheme_module._addBuiltinFunction("sum_of_R2x2",
std::make_shared<
BuiltinFunctionEmbedder<TinyMatrix<2>(std::shared_ptr<const IDiscreteFunction>)>>(
[](std::shared_ptr<const IDiscreteFunction> a) -> TinyMatrix<2> {
- return sum_to<TinyMatrix<2>>(a);
+ return sum_of<TinyMatrix<2>>(a);
}));
- scheme_module._addBuiltinFunction("sum_to_R3x3",
+ scheme_module._addBuiltinFunction("sum_of_R3x3",
std::make_shared<
BuiltinFunctionEmbedder<TinyMatrix<3>(std::shared_ptr<const IDiscreteFunction>)>>(
[](std::shared_ptr<const IDiscreteFunction> a) -> TinyMatrix<3> {
- return sum_to<TinyMatrix<3>>(a);
+ return sum_of<TinyMatrix<3>>(a);
+ }));
+ scheme_module
+ ._addBuiltinFunction("integral_of_R",
+ std::make_shared<BuiltinFunctionEmbedder<double(std::shared_ptr<const IDiscreteFunction>)>>(
+ [](std::shared_ptr<const IDiscreteFunction> a) -> double {
+ return integral_of<double>(a);
+ }));
+
+ scheme_module._addBuiltinFunction("integral_of_R1",
+ std::make_shared<
+ BuiltinFunctionEmbedder<TinyVector<1>(std::shared_ptr<const IDiscreteFunction>)>>(
+ [](std::shared_ptr<const IDiscreteFunction> a) -> TinyVector<1> {
+ return integral_of<TinyVector<1>>(a);
+ }));
+
+ scheme_module._addBuiltinFunction("integral_of_R2",
+ std::make_shared<
+ BuiltinFunctionEmbedder<TinyVector<2>(std::shared_ptr<const IDiscreteFunction>)>>(
+ [](std::shared_ptr<const IDiscreteFunction> a) -> TinyVector<2> {
+ return integral_of<TinyVector<2>>(a);
+ }));
+
+ scheme_module._addBuiltinFunction("integral_of_R3",
+ std::make_shared<
+ BuiltinFunctionEmbedder<TinyVector<3>(std::shared_ptr<const IDiscreteFunction>)>>(
+ [](std::shared_ptr<const IDiscreteFunction> a) -> TinyVector<3> {
+ return integral_of<TinyVector<3>>(a);
+ }));
+
+ scheme_module._addBuiltinFunction("integral_of_R1x1",
+ std::make_shared<
+ BuiltinFunctionEmbedder<TinyMatrix<1>(std::shared_ptr<const IDiscreteFunction>)>>(
+ [](std::shared_ptr<const IDiscreteFunction> a) -> TinyMatrix<1> {
+ return integral_of<TinyMatrix<1>>(a);
+ }));
+
+ scheme_module._addBuiltinFunction("integral_of_R2x2",
+ std::make_shared<
+ BuiltinFunctionEmbedder<TinyMatrix<2>(std::shared_ptr<const IDiscreteFunction>)>>(
+ [](std::shared_ptr<const IDiscreteFunction> a) -> TinyMatrix<2> {
+ return integral_of<TinyMatrix<2>>(a);
+ }));
+
+ scheme_module._addBuiltinFunction("integral_of_R3x3",
+ std::make_shared<
+ BuiltinFunctionEmbedder<TinyMatrix<3>(std::shared_ptr<const IDiscreteFunction>)>>(
+ [](std::shared_ptr<const IDiscreteFunction> a) -> TinyMatrix<3> {
+ return integral_of<TinyMatrix<3>>(a);
}));
}
diff --git a/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.cpp b/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.cpp
index 0c5f23fb63210a8f72ddda25688fc758f743c13c..a688acf1372f93b10192fec99abe0bd15de52683 100644
--- a/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.cpp
+++ b/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.cpp
@@ -722,7 +722,7 @@ max(const std::shared_ptr<const IDiscreteFunction>& f, const double a)
template <typename ValueT>
ValueT
-sum_to(const std::shared_ptr<const IDiscreteFunction>& f)
+sum_of(const std::shared_ptr<const IDiscreteFunction>& f)
{
if (f->dataType() == ast_node_data_type_from<ValueT> and f->descriptor().type() == DiscreteFunctionType::P0) {
switch (f->mesh()->dimension()) {
@@ -747,16 +747,57 @@ sum_to(const std::shared_ptr<const IDiscreteFunction>& f)
}
}
-template double sum_to<double>(const std::shared_ptr<const IDiscreteFunction>&);
+template double sum_of<double>(const std::shared_ptr<const IDiscreteFunction>&);
-template TinyVector<1> sum_to<TinyVector<1>>(const std::shared_ptr<const IDiscreteFunction>&);
+template TinyVector<1> sum_of<TinyVector<1>>(const std::shared_ptr<const IDiscreteFunction>&);
-template TinyVector<2> sum_to<TinyVector<2>>(const std::shared_ptr<const IDiscreteFunction>&);
+template TinyVector<2> sum_of<TinyVector<2>>(const std::shared_ptr<const IDiscreteFunction>&);
-template TinyVector<3> sum_to<TinyVector<3>>(const std::shared_ptr<const IDiscreteFunction>&);
+template TinyVector<3> sum_of<TinyVector<3>>(const std::shared_ptr<const IDiscreteFunction>&);
-template TinyMatrix<1> sum_to<TinyMatrix<1>>(const std::shared_ptr<const IDiscreteFunction>&);
+template TinyMatrix<1> sum_of<TinyMatrix<1>>(const std::shared_ptr<const IDiscreteFunction>&);
-template TinyMatrix<2> sum_to<TinyMatrix<2>>(const std::shared_ptr<const IDiscreteFunction>&);
+template TinyMatrix<2> sum_of<TinyMatrix<2>>(const std::shared_ptr<const IDiscreteFunction>&);
-template TinyMatrix<3> sum_to<TinyMatrix<3>>(const std::shared_ptr<const IDiscreteFunction>&);
+template TinyMatrix<3> sum_of<TinyMatrix<3>>(const std::shared_ptr<const IDiscreteFunction>&);
+
+template <typename ValueT>
+ValueT
+integral_of(const std::shared_ptr<const IDiscreteFunction>& f)
+{
+ if (f->dataType() == ast_node_data_type_from<ValueT> and f->descriptor().type() == DiscreteFunctionType::P0) {
+ switch (f->mesh()->dimension()) {
+ case 1: {
+ using DiscreteFunctionType = DiscreteFunctionP0<1, ValueT>;
+ return integrate(dynamic_cast<const DiscreteFunctionType&>(*f));
+ }
+ case 2: {
+ using DiscreteFunctionType = DiscreteFunctionP0<2, ValueT>;
+ return integrate(dynamic_cast<const DiscreteFunctionType&>(*f));
+ }
+ case 3: {
+ using DiscreteFunctionType = DiscreteFunctionP0<3, ValueT>;
+ return integrate(dynamic_cast<const DiscreteFunctionType&>(*f));
+ }
+ default: {
+ throw UnexpectedError("invalid mesh dimension");
+ }
+ }
+ } else {
+ throw NormalError("invalid operand type " + operand_type_name(f));
+ }
+}
+
+template double integral_of<double>(const std::shared_ptr<const IDiscreteFunction>&);
+
+template TinyVector<1> integral_of<TinyVector<1>>(const std::shared_ptr<const IDiscreteFunction>&);
+
+template TinyVector<2> integral_of<TinyVector<2>>(const std::shared_ptr<const IDiscreteFunction>&);
+
+template TinyVector<3> integral_of<TinyVector<3>>(const std::shared_ptr<const IDiscreteFunction>&);
+
+template TinyMatrix<1> integral_of<TinyMatrix<1>>(const std::shared_ptr<const IDiscreteFunction>&);
+
+template TinyMatrix<2> integral_of<TinyMatrix<2>>(const std::shared_ptr<const IDiscreteFunction>&);
+
+template TinyMatrix<3> integral_of<TinyMatrix<3>>(const std::shared_ptr<const IDiscreteFunction>&);
diff --git a/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.hpp b/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.hpp
index 6eef58101016b95411613677cef2d27e5e52eac1..582e63d8551842e9f93e3962a08669ba34baf756 100644
--- a/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.hpp
+++ b/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.hpp
@@ -84,6 +84,9 @@ std::shared_ptr<const IDiscreteFunction> max(const double, const std::shared_ptr
std::shared_ptr<const IDiscreteFunction> max(const std::shared_ptr<const IDiscreteFunction>&, const double);
template <typename ValueT>
-ValueT sum_to(const std::shared_ptr<const IDiscreteFunction>&);
+ValueT sum_of(const std::shared_ptr<const IDiscreteFunction>&);
+
+template <typename ValueT>
+ValueT integral_of(const std::shared_ptr<const IDiscreteFunction>&);
#endif // EMBEDDED_I_DISCRETE_FUNCTION_MATH_FUNCTIONS_HPP
diff --git a/src/mesh/ItemValueUtils.hpp b/src/mesh/ItemValueUtils.hpp
index c307ce4d6010e9ce0d40b4afdcb2ab4466e1d009..ab39d5c5c6423dac3d8c7d3fd2c1cbeffd1a9809 100644
--- a/src/mesh/ItemValueUtils.hpp
+++ b/src/mesh/ItemValueUtils.hpp
@@ -7,6 +7,7 @@
#include <mesh/ItemValue.hpp>
#include <mesh/Synchronizer.hpp>
#include <mesh/SynchronizerManager.hpp>
+#include <utils/PugsTraits.hpp>
#include <iostream>
@@ -244,6 +245,7 @@ sum(const ItemValue<DataType, item_type>& item_value)
if constexpr (std::is_arithmetic_v<data_type>) {
value = 0;
} else {
+ static_assert(is_tiny_vector_v<data_type> or is_tiny_matrix_v<data_type>, "invalid data type");
value = zero;
}
}
diff --git a/src/scheme/DiscreteFunctionP0.hpp b/src/scheme/DiscreteFunctionP0.hpp
index ccedb957d9e1fdb9b0c6643d66618709a9717cad..bb6b484a6bc10bd803e16b8bb54a7cba25a349de 100644
--- a/src/scheme/DiscreteFunctionP0.hpp
+++ b/src/scheme/DiscreteFunctionP0.hpp
@@ -659,6 +659,73 @@ class DiscreteFunctionP0 : public IDiscreteFunction
return sum(f.m_cell_values);
}
+ private:
+ class Integrator
+ {
+ private:
+ const DiscreteFunctionP0 m_function;
+ CellValue<const double> m_cell_volume;
+ CellValue<const bool> m_cell_is_owned;
+
+ public:
+ PUGS_INLINE
+ operator data_type()
+ {
+ data_type reduced_value;
+ parallel_reduce(m_cell_volume.numberOfItems(), *this, reduced_value);
+ return reduced_value;
+ }
+
+ PUGS_INLINE
+ void
+ operator()(const CellId& cell_id, data_type& data) const
+ {
+ if (m_cell_is_owned[cell_id]) {
+ data += m_cell_volume[cell_id] * m_function[cell_id];
+ }
+ }
+
+ PUGS_INLINE
+ void
+ join(volatile data_type& dst, const volatile data_type& src) const
+ {
+ dst += src;
+ }
+
+ PUGS_INLINE
+ void
+ init(data_type& value) const
+ {
+ if constexpr (std::is_arithmetic_v<data_type>) {
+ value = 0;
+ } else {
+ static_assert(is_tiny_vector_v<data_type> or is_tiny_matrix_v<data_type>, "invalid data type");
+ value = zero;
+ }
+ }
+
+ PUGS_INLINE
+ Integrator(const DiscreteFunctionP0& function) : m_function(function)
+ {
+ const Mesh<Connectivity<Dimension>>& mesh = *function.m_mesh;
+ m_cell_volume = MeshDataManager::instance().getMeshData(mesh).Vj();
+ m_cell_is_owned = mesh.connectivity().cellIsOwned();
+ }
+
+ PUGS_INLINE
+ ~Integrator() = default;
+ };
+
+ public:
+ PUGS_INLINE friend DataType
+ integrate(const DiscreteFunctionP0& f)
+ {
+ Assert(f.m_cell_values.isBuilt());
+ const DataType integral = Integrator{f};
+
+ return parallel::allReduceSum(integral);
+ }
+
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)
diff --git a/src/utils/ArrayUtils.hpp b/src/utils/ArrayUtils.hpp
index d85fb162756deb20d19971470ea1f40421dfdfa1..3b3438294291af1deafe733374e7de9ddf4982b6 100644
--- a/src/utils/ArrayUtils.hpp
+++ b/src/utils/ArrayUtils.hpp
@@ -2,6 +2,7 @@
#define ARRAY_UTILS_HPP
#include <utils/PugsMacros.hpp>
+#include <utils/PugsTraits.hpp>
#include <utils/PugsUtils.hpp>
#include <utils/Types.hpp>
@@ -137,7 +138,7 @@ sum(const ArrayT<DataType>& array)
class ArraySum
{
private:
- const ArrayType& m_array;
+ const ArrayType m_array;
public:
PUGS_INLINE
@@ -169,6 +170,7 @@ sum(const ArrayT<DataType>& array)
if constexpr (std::is_arithmetic_v<data_type>) {
value = 0;
} else {
+ static_assert(is_tiny_vector_v<data_type> or is_tiny_matrix_v<data_type>, "invalid data type");
value = zero;
}
}
diff --git a/src/utils/PugsUtils.hpp b/src/utils/PugsUtils.hpp
index 91267a2af92928084ce07e1021baf6336e8a8f98..76d8ffd3feafd53d6fa2d5326fa18e8e65d95ca6 100644
--- a/src/utils/PugsUtils.hpp
+++ b/src/utils/PugsUtils.hpp
@@ -13,11 +13,11 @@ parallel_for(size_t size, const FunctionType& lambda, const std::string& label =
Kokkos::parallel_for(size, lambda, label);
}
-template <typename ArrayType, typename ReturnType>
+template <typename FunctorType, typename ReturnType>
PUGS_FORCEINLINE void
-parallel_reduce(size_t size, const ArrayType& array, ReturnType& value, const std::string& label = "")
+parallel_reduce(size_t size, const FunctorType& functor, ReturnType& value, const std::string& label = "")
{
- Kokkos::parallel_reduce(label, size, array, value);
+ Kokkos::parallel_reduce(label, size, functor, value);
}
void setDefaultOMPEnvironment();