From fe9878b58dc8792216c0580084cdb93552db9875 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Del=20Pino?= <stephane.delpino44@gmail.com>
Date: Fri, 26 Nov 2021 12:56:32 +0100
Subject: [PATCH] Add DiscreteFunctionIntegrator and
DiscreteFunctionVectorIntegrator
It is now possible to integrate (user) functions within a pugs script
by specifying the mesh and the quadrature formula kind.
The result if a P0 (or P0Vector) discrete function which contains
integrals computed in each cell of the mesh
---
.../GaussLegendreQuadratureDescriptor.hpp | 50 ++
.../GaussLobattoQuadratureDescriptor.hpp | 50 ++
src/analysis/GaussQuadratureDescriptor.hpp | 50 ++
src/analysis/IQuadratureDescriptor.hpp | 25 +
src/analysis/QuadratureManager.hpp | 111 ++++
src/analysis/QuadratureType.hpp | 6 +-
src/language/modules/SchemeModule.cpp | 64 ++
src/language/modules/SchemeModule.hpp | 5 +
src/language/utils/IntegrateCellArray.hpp | 69 ++
src/language/utils/IntegrateCellValue.hpp | 52 ++
src/language/utils/IntegrateOnCells.hpp | 606 ++++++++++++++++++
src/scheme/CMakeLists.txt | 2 +
src/scheme/DiscreteFunctionIntegrator.cpp | 132 ++++
src/scheme/DiscreteFunctionIntegrator.hpp | 39 ++
.../DiscreteFunctionVectorIntegrator.cpp | 70 ++
.../DiscreteFunctionVectorIntegrator.hpp | 47 ++
16 files changed, 1373 insertions(+), 5 deletions(-)
create mode 100644 src/analysis/GaussLegendreQuadratureDescriptor.hpp
create mode 100644 src/analysis/GaussLobattoQuadratureDescriptor.hpp
create mode 100644 src/analysis/GaussQuadratureDescriptor.hpp
create mode 100644 src/analysis/IQuadratureDescriptor.hpp
create mode 100644 src/language/utils/IntegrateCellArray.hpp
create mode 100644 src/language/utils/IntegrateCellValue.hpp
create mode 100644 src/language/utils/IntegrateOnCells.hpp
create mode 100644 src/scheme/DiscreteFunctionIntegrator.cpp
create mode 100644 src/scheme/DiscreteFunctionIntegrator.hpp
create mode 100644 src/scheme/DiscreteFunctionVectorIntegrator.cpp
create mode 100644 src/scheme/DiscreteFunctionVectorIntegrator.hpp
diff --git a/src/analysis/GaussLegendreQuadratureDescriptor.hpp b/src/analysis/GaussLegendreQuadratureDescriptor.hpp
new file mode 100644
index 000000000..4cccc87e0
--- /dev/null
+++ b/src/analysis/GaussLegendreQuadratureDescriptor.hpp
@@ -0,0 +1,50 @@
+#ifndef GAUSS_LEGENDRE_QUADRATURE_DESCRIPTOR_HPP
+#define GAUSS_LEGENDRE_QUADRATURE_DESCRIPTOR_HPP
+
+#include <analysis/IQuadratureDescriptor.hpp>
+
+#include <sstream>
+
+class GaussLegendreQuadratureDescriptor final : public IQuadratureDescriptor
+{
+ private:
+ size_t m_degree;
+
+ public:
+ bool
+ isTensorial() const
+ {
+ return true;
+ }
+
+ QuadratureType
+ type() const
+ {
+ return QuadratureType::GaussLegendre;
+ }
+
+ size_t
+ degree() const
+ {
+ return m_degree;
+ }
+
+ std::string
+ name() const
+ {
+ std::ostringstream os;
+ os << "GaussLegendre(" << m_degree << ")";
+ return os.str();
+ }
+
+ GaussLegendreQuadratureDescriptor(size_t degree) : m_degree{degree} {}
+ GaussLegendreQuadratureDescriptor() noexcept = delete;
+
+ GaussLegendreQuadratureDescriptor(const GaussLegendreQuadratureDescriptor&) = default;
+
+ GaussLegendreQuadratureDescriptor(GaussLegendreQuadratureDescriptor&&) noexcept = default;
+
+ virtual ~GaussLegendreQuadratureDescriptor() noexcept = default;
+};
+
+#endif // GAUSS_LEGENDRE_QUADRATURE_DESCRIPTOR_HPP
diff --git a/src/analysis/GaussLobattoQuadratureDescriptor.hpp b/src/analysis/GaussLobattoQuadratureDescriptor.hpp
new file mode 100644
index 000000000..5b939437b
--- /dev/null
+++ b/src/analysis/GaussLobattoQuadratureDescriptor.hpp
@@ -0,0 +1,50 @@
+#ifndef GAUSS_LOBATTO_QUADRATURE_DESCRIPTOR_HPP
+#define GAUSS_LOBATTO_QUADRATURE_DESCRIPTOR_HPP
+
+#include <analysis/IQuadratureDescriptor.hpp>
+
+#include <sstream>
+
+class GaussLobattoQuadratureDescriptor final : public IQuadratureDescriptor
+{
+ private:
+ size_t m_degree;
+
+ public:
+ bool
+ isTensorial() const
+ {
+ return true;
+ }
+
+ QuadratureType
+ type() const
+ {
+ return QuadratureType::GaussLobatto;
+ }
+
+ size_t
+ degree() const
+ {
+ return m_degree;
+ }
+
+ std::string
+ name() const
+ {
+ std::ostringstream os;
+ os << "GaussLobatto(" << m_degree << ")";
+ return os.str();
+ }
+
+ GaussLobattoQuadratureDescriptor(size_t degree) : m_degree{degree} {}
+ GaussLobattoQuadratureDescriptor() noexcept = delete;
+
+ GaussLobattoQuadratureDescriptor(const GaussLobattoQuadratureDescriptor&) = default;
+
+ GaussLobattoQuadratureDescriptor(GaussLobattoQuadratureDescriptor&&) noexcept = default;
+
+ virtual ~GaussLobattoQuadratureDescriptor() noexcept = default;
+};
+
+#endif // GAUSS_LOBATTO_QUADRATURE_DESCRIPTOR_HPP
diff --git a/src/analysis/GaussQuadratureDescriptor.hpp b/src/analysis/GaussQuadratureDescriptor.hpp
new file mode 100644
index 000000000..13282f65d
--- /dev/null
+++ b/src/analysis/GaussQuadratureDescriptor.hpp
@@ -0,0 +1,50 @@
+#ifndef GAUSS_QUADRATURE_DESCRIPTOR_HPP
+#define GAUSS_QUADRATURE_DESCRIPTOR_HPP
+
+#include <analysis/IQuadratureDescriptor.hpp>
+
+#include <sstream>
+
+class GaussQuadratureDescriptor final : public IQuadratureDescriptor
+{
+ private:
+ size_t m_degree;
+
+ public:
+ bool
+ isTensorial() const
+ {
+ return false;
+ }
+
+ QuadratureType
+ type() const
+ {
+ return QuadratureType::Gauss;
+ }
+
+ size_t
+ degree() const
+ {
+ return m_degree;
+ }
+
+ std::string
+ name() const
+ {
+ std::ostringstream os;
+ os << "Gauss(" << m_degree << ")";
+ return os.str();
+ }
+
+ GaussQuadratureDescriptor(size_t degree) : m_degree{degree} {}
+ GaussQuadratureDescriptor() noexcept = delete;
+
+ GaussQuadratureDescriptor(const GaussQuadratureDescriptor&) = default;
+
+ GaussQuadratureDescriptor(GaussQuadratureDescriptor&&) noexcept = default;
+
+ virtual ~GaussQuadratureDescriptor() noexcept = default;
+};
+
+#endif // GAUSS_QUADRATURE_DESCRIPTOR_HPP
diff --git a/src/analysis/IQuadratureDescriptor.hpp b/src/analysis/IQuadratureDescriptor.hpp
new file mode 100644
index 000000000..5d7e3685a
--- /dev/null
+++ b/src/analysis/IQuadratureDescriptor.hpp
@@ -0,0 +1,25 @@
+#ifndef I_QUADRATURE_DESCRIPTOR_HPP
+#define I_QUADRATURE_DESCRIPTOR_HPP
+
+#include <analysis/QuadratureType.hpp>
+
+#include <string>
+
+class IQuadratureDescriptor
+{
+ public:
+ virtual QuadratureType type() const = 0;
+ virtual bool isTensorial() const = 0;
+ virtual size_t degree() const = 0;
+ virtual std::string name() const = 0;
+
+ IQuadratureDescriptor() noexcept = default;
+
+ IQuadratureDescriptor(const IQuadratureDescriptor&) = default;
+
+ IQuadratureDescriptor(IQuadratureDescriptor&&) noexcept = default;
+
+ virtual ~IQuadratureDescriptor() noexcept = default;
+};
+
+#endif // I_QUADRATURE_DESCRIPTOR_HPP
diff --git a/src/analysis/QuadratureManager.hpp b/src/analysis/QuadratureManager.hpp
index 875a2a3fd..d205c457f 100644
--- a/src/analysis/QuadratureManager.hpp
+++ b/src/analysis/QuadratureManager.hpp
@@ -1,8 +1,10 @@
#ifndef QUADRATURE_MANAGER_HPP
#define QUADRATURE_MANAGER_HPP
+#include <analysis/IQuadratureDescriptor.hpp>
#include <analysis/QuadratureFormula.hpp>
#include <utils/Array.hpp>
+#include <utils/Exceptions.hpp>
#include <utils/PugsAssert.hpp>
#include <utils/PugsMacros.hpp>
@@ -60,6 +62,115 @@ class QuadratureManager
return m_line_gauss_lobatto_formula_list.size() * 2 - 1;
}
+ const QuadratureFormula<1>&
+ getLineFormula(const IQuadratureDescriptor& quadrature_descriptor) const
+ {
+ switch (quadrature_descriptor.type()) {
+ case QuadratureType::Gauss:
+ case QuadratureType::GaussLegendre: {
+ return m_line_gauss_legendre_formula_list[quadrature_descriptor.degree() / 2];
+ break;
+ }
+ case QuadratureType::GaussLobatto: {
+ return m_line_gauss_lobatto_formula_list[quadrature_descriptor.degree() / 2];
+ }
+ default: {
+ throw UnexpectedError("invalid quadrature type");
+ }
+ }
+ }
+
+ const QuadratureFormula<2>&
+ getTriangleFormula(const IQuadratureDescriptor& quadrature_descriptor) const
+ {
+ switch (quadrature_descriptor.type()) {
+ case QuadratureType::Gauss: {
+ return m_triangle_gauss_formula_list[quadrature_descriptor.degree() - 1];
+ }
+ default: {
+ throw UnexpectedError(quadrature_descriptor.name() + " is not defined on triangles");
+ }
+ }
+ }
+
+ const QuadratureFormula<2>&
+ getSquareFormula(const IQuadratureDescriptor& quadrature_descriptor) const
+ {
+ switch (quadrature_descriptor.type()) {
+ case QuadratureType::Gauss: {
+ return m_square_gauss_formula_list[quadrature_descriptor.degree() / 2];
+ }
+ case QuadratureType::GaussLegendre: {
+ return m_square_gauss_legendre_formula_list[quadrature_descriptor.degree() / 2];
+ }
+ case QuadratureType::GaussLobatto: {
+ return m_square_gauss_lobatto_formula_list[quadrature_descriptor.degree() / 2];
+ }
+ default: {
+ throw UnexpectedError("invalid quadrature type");
+ }
+ }
+ }
+
+ const QuadratureFormula<3>&
+ getTetrahedronFormula(const IQuadratureDescriptor& quadrature_descriptor) const
+ {
+ switch (quadrature_descriptor.type()) {
+ case QuadratureType::Gauss: {
+ return m_tetrahedron_gauss_formula_list[quadrature_descriptor.degree() - 1];
+ }
+ default: {
+ throw UnexpectedError(quadrature_descriptor.name() + " is not defined on tetrahedron");
+ }
+ }
+ }
+
+ const QuadratureFormula<3>&
+ getPrismFormula(const IQuadratureDescriptor& quadrature_descriptor) const
+ {
+ switch (quadrature_descriptor.type()) {
+ case QuadratureType::Gauss: {
+ return m_prism_gauss_formula_list[quadrature_descriptor.degree() - 1];
+ }
+ default: {
+ throw UnexpectedError(quadrature_descriptor.name() + " is not defined on prism");
+ }
+ }
+ }
+
+ const QuadratureFormula<3>&
+ getPyramidFormula(const IQuadratureDescriptor& quadrature_descriptor) const
+ {
+ switch (quadrature_descriptor.type()) {
+ case QuadratureType::Gauss: {
+ return m_pyramid_gauss_formula_list[quadrature_descriptor.degree() - 1];
+ }
+ default: {
+ throw UnexpectedError(quadrature_descriptor.name() + " is not defined on pyramid");
+ }
+ }
+ }
+
+ const QuadratureFormula<3>&
+ getCubeFormula(const IQuadratureDescriptor& quadrature_descriptor) const
+ {
+ switch (quadrature_descriptor.type()) {
+ case QuadratureType::Gauss: {
+ return m_cube_gauss_formula_list[quadrature_descriptor.degree() / 2];
+ }
+ case QuadratureType::GaussLegendre: {
+ return m_cube_gauss_legendre_formula_list[quadrature_descriptor.degree() / 2];
+ break;
+ }
+ case QuadratureType::GaussLobatto: {
+ return m_cube_gauss_lobatto_formula_list[quadrature_descriptor.degree() / 2];
+ }
+ default: {
+ throw UnexpectedError("invalid quadrature type");
+ }
+ }
+ }
+
const QuadratureFormula<1>&
getLineGaussLegendreFormula(const size_t degree) const
{
diff --git a/src/analysis/QuadratureType.hpp b/src/analysis/QuadratureType.hpp
index dca48806d..795961591 100644
--- a/src/analysis/QuadratureType.hpp
+++ b/src/analysis/QuadratureType.hpp
@@ -3,13 +3,9 @@
enum class QuadratureType
{
- QT__begin = 0,
- //
- Gauss = QT__begin,
+ Gauss = 0,
GaussLegendre = 1,
GaussLobatto = 2,
- //
- QT__end
};
#endif // QUADRATURE_TYPE_HPP
diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index d989bff45..fcebda3de 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -1,5 +1,8 @@
#include <language/modules/SchemeModule.hpp>
+#include <analysis/GaussLegendreQuadratureDescriptor.hpp>
+#include <analysis/GaussLobattoQuadratureDescriptor.hpp>
+#include <analysis/GaussQuadratureDescriptor.hpp>
#include <language/modules/BinaryOperatorRegisterForVh.hpp>
#include <language/modules/MathFunctionRegisterForVh.hpp>
#include <language/modules/UnaryOperatorRegisterForVh.hpp>
@@ -17,9 +20,11 @@
#include <scheme/DirichletBoundaryConditionDescriptor.hpp>
#include <scheme/DiscreteFunctionDescriptorP0.hpp>
#include <scheme/DiscreteFunctionDescriptorP0Vector.hpp>
+#include <scheme/DiscreteFunctionIntegrator.hpp>
#include <scheme/DiscreteFunctionInterpoler.hpp>
#include <scheme/DiscreteFunctionP0.hpp>
#include <scheme/DiscreteFunctionUtils.hpp>
+#include <scheme/DiscreteFunctionVectorIntegrator.hpp>
#include <scheme/DiscreteFunctionVectorInterpoler.hpp>
#include <scheme/FixedBoundaryConditionDescriptor.hpp>
#include <scheme/FourierBoundaryConditionDescriptor.hpp>
@@ -38,6 +43,7 @@ SchemeModule::SchemeModule()
{
this->_addTypeDescriptor(ast_node_data_type_from<std::shared_ptr<const IDiscreteFunction>>);
this->_addTypeDescriptor(ast_node_data_type_from<std::shared_ptr<const IDiscreteFunctionDescriptor>>);
+ this->_addTypeDescriptor(ast_node_data_type_from<std::shared_ptr<const IQuadratureDescriptor>>);
this->_addTypeDescriptor(ast_node_data_type_from<std::shared_ptr<const IBoundaryDescriptor>>);
this->_addTypeDescriptor(ast_node_data_type_from<std::shared_ptr<const IBoundaryConditionDescriptor>>);
@@ -59,6 +65,64 @@ SchemeModule::SchemeModule()
));
+ this->_addBuiltinFunction("Gauss", std::make_shared<
+ BuiltinFunctionEmbedder<std::shared_ptr<const IQuadratureDescriptor>(uint64_t)>>(
+ [](uint64_t degree) -> std::shared_ptr<const IQuadratureDescriptor> {
+ return std::make_shared<GaussQuadratureDescriptor>(degree);
+ }
+
+ ));
+
+ this->_addBuiltinFunction("GaussLobatto",
+ std::make_shared<
+ BuiltinFunctionEmbedder<std::shared_ptr<const IQuadratureDescriptor>(uint64_t)>>(
+ [](uint64_t degree) -> std::shared_ptr<const IQuadratureDescriptor> {
+ return std::make_shared<GaussLobattoQuadratureDescriptor>(degree);
+ }
+
+ ));
+
+ this->_addBuiltinFunction("GaussLegendre",
+ std::make_shared<
+ BuiltinFunctionEmbedder<std::shared_ptr<const IQuadratureDescriptor>(uint64_t)>>(
+ [](uint64_t degree) -> std::shared_ptr<const IQuadratureDescriptor> {
+ return std::make_shared<GaussLegendreQuadratureDescriptor>(degree);
+ }
+
+ ));
+
+ this
+ ->_addBuiltinFunction("integrate",
+ std::make_shared<BuiltinFunctionEmbedder<
+ std::shared_ptr<const IDiscreteFunction>(std::shared_ptr<const IMesh>,
+ std::shared_ptr<const IQuadratureDescriptor>,
+ std::shared_ptr<const IDiscreteFunctionDescriptor>,
+ const std::vector<FunctionSymbolId>&)>>(
+ [](std::shared_ptr<const IMesh> mesh,
+ std::shared_ptr<const IQuadratureDescriptor> quadrature_descriptor,
+ std::shared_ptr<const IDiscreteFunctionDescriptor> discrete_function_descriptor,
+ const std::vector<FunctionSymbolId>& function_id_list)
+ -> std::shared_ptr<const IDiscreteFunction> {
+ return DiscreteFunctionVectorIntegrator{mesh, quadrature_descriptor,
+ discrete_function_descriptor, function_id_list}
+ .integrate();
+ }
+
+ ));
+
+ this->_addBuiltinFunction("integrate",
+ std::make_shared<BuiltinFunctionEmbedder<
+ std::shared_ptr<const IDiscreteFunction>(std::shared_ptr<const IMesh>,
+ std::shared_ptr<const IQuadratureDescriptor>,
+ const FunctionSymbolId&)>>(
+ [](std::shared_ptr<const IMesh> mesh,
+ std::shared_ptr<const IQuadratureDescriptor> quadrature_descriptor,
+ const FunctionSymbolId& function_id) -> std::shared_ptr<const IDiscreteFunction> {
+ return DiscreteFunctionIntegrator{mesh, quadrature_descriptor, function_id}.integrate();
+ }
+
+ ));
+
this->_addBuiltinFunction(
"interpolate",
std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<
diff --git a/src/language/modules/SchemeModule.hpp b/src/language/modules/SchemeModule.hpp
index 758b5f9c5..52fc163ee 100644
--- a/src/language/modules/SchemeModule.hpp
+++ b/src/language/modules/SchemeModule.hpp
@@ -25,6 +25,11 @@ template <>
inline ASTNodeDataType ast_node_data_type_from<std::shared_ptr<const IDiscreteFunctionDescriptor>> =
ASTNodeDataType::build<ASTNodeDataType::type_id_t>("Vh_type");
+class IQuadratureDescriptor;
+template <>
+inline ASTNodeDataType ast_node_data_type_from<std::shared_ptr<const IQuadratureDescriptor>> =
+ ASTNodeDataType::build<ASTNodeDataType::type_id_t>("quadrature");
+
class SchemeModule : public BuiltinModule
{
friend class MathFunctionRegisterForVh;
diff --git a/src/language/utils/IntegrateCellArray.hpp b/src/language/utils/IntegrateCellArray.hpp
new file mode 100644
index 000000000..41450de3b
--- /dev/null
+++ b/src/language/utils/IntegrateCellArray.hpp
@@ -0,0 +1,69 @@
+#ifndef INTEGRATE_CELL_ARRAY_HPP
+#define INTEGRATE_CELL_ARRAY_HPP
+
+#include <language/utils/IntegrateCellValue.hpp>
+#include <mesh/CellType.hpp>
+#include <mesh/ItemArray.hpp>
+
+template <typename T>
+class IntegrateCellArray;
+template <typename OutputType, typename InputType>
+class IntegrateCellArray<OutputType(InputType)>
+{
+ static constexpr size_t Dimension = OutputType::Dimension;
+
+ public:
+ template <typename MeshType>
+ PUGS_INLINE static CellArray<OutputType>
+ integrate(const std::vector<FunctionSymbolId>& function_symbol_id_list,
+ const IQuadratureDescriptor& quadrature_descriptor,
+ const MeshType& mesh)
+ {
+ CellArray<OutputType> cell_array{mesh.connectivity(), function_symbol_id_list.size()};
+
+ for (size_t i_function_symbol = 0; i_function_symbol < function_symbol_id_list.size(); ++i_function_symbol) {
+ const FunctionSymbolId& function_symbol_id = function_symbol_id_list[i_function_symbol];
+ CellValue<OutputType> cell_value =
+ IntegrateCellValue<OutputType(InputType)>::integrate(function_symbol_id, quadrature_descriptor, mesh);
+ parallel_for(
+ cell_value.numberOfItems(),
+ PUGS_LAMBDA(CellId cell_id) { cell_array[cell_id][i_function_symbol] = cell_value[cell_id]; });
+ }
+
+ return cell_array;
+ }
+
+ template <typename MeshType>
+ PUGS_INLINE static Table<OutputType>
+ integrate(const std::vector<FunctionSymbolId>& function_symbol_id_list,
+ const IQuadratureDescriptor& quadrature_descriptor,
+ const MeshType& mesh,
+ const Array<const CellId>& list_of_cells)
+ {
+ Table<OutputType> table{list_of_cells.size(), function_symbol_id_list.size()};
+
+ for (size_t i_function_symbol = 0; i_function_symbol < function_symbol_id_list.size(); ++i_function_symbol) {
+ const FunctionSymbolId& function_symbol_id = function_symbol_id_list[i_function_symbol];
+ Array<OutputType> array =
+ IntegrateCellValue<OutputType(InputType)>::integrate(function_symbol_id, quadrature_descriptor, mesh,
+ list_of_cells);
+
+ parallel_for(
+ array.size(), PUGS_LAMBDA(size_t i) { table[i][i_function_symbol] = array[i]; });
+ }
+
+ return table;
+ }
+
+ template <typename MeshType>
+ PUGS_INLINE static Table<OutputType>
+ integrate(const std::vector<FunctionSymbolId>& function_symbol_id_list,
+ const IQuadratureDescriptor& quadrature_descriptor,
+ const MeshType& mesh,
+ const Array<CellId>& list_of_cells)
+ {
+ return integrate(function_symbol_id_list, quadrature_descriptor, mesh, Array<const CellId>{list_of_cells});
+ }
+};
+
+#endif // INTEGRATE_CELL_ARRAY_HPP
diff --git a/src/language/utils/IntegrateCellValue.hpp b/src/language/utils/IntegrateCellValue.hpp
new file mode 100644
index 000000000..c4eafc151
--- /dev/null
+++ b/src/language/utils/IntegrateCellValue.hpp
@@ -0,0 +1,52 @@
+#ifndef INTEGRATE_CELL_VALUE_HPP
+#define INTEGRATE_CELL_VALUE_HPP
+
+#include <analysis/IQuadratureDescriptor.hpp>
+#include <language/utils/IntegrateOnCells.hpp>
+#include <mesh/ItemType.hpp>
+#include <mesh/ItemValue.hpp>
+#include <mesh/Mesh.hpp>
+
+template <typename T>
+class IntegrateCellValue;
+template <typename OutputType, typename InputType>
+class IntegrateCellValue<OutputType(InputType)>
+{
+ public:
+ template <typename MeshType>
+ PUGS_INLINE static CellValue<OutputType>
+ integrate(const FunctionSymbolId& function_symbol_id,
+ const IQuadratureDescriptor& quadrature_descriptor,
+ const MeshType& mesh)
+ {
+ CellValue<OutputType> value(mesh.connectivity());
+ IntegrateOnCells<OutputType(const InputType)>::template integrateTo<MeshType>(function_symbol_id,
+ quadrature_descriptor, mesh, value);
+
+ return value;
+ }
+
+ template <typename MeshType>
+ PUGS_INLINE static Array<OutputType>
+ integrate(const FunctionSymbolId& function_symbol_id,
+ const IQuadratureDescriptor& quadrature_descriptor,
+ const MeshType& mesh,
+ const Array<const CellId>& list_of_cells)
+ {
+ return IntegrateOnCells<OutputType(const InputType)>::integrate(function_symbol_id, quadrature_descriptor, mesh,
+ list_of_cells);
+ }
+
+ template <typename MeshType>
+ PUGS_INLINE static Array<OutputType>
+ integrate(const FunctionSymbolId& function_symbol_id,
+ const IQuadratureDescriptor& quadrature_descriptor,
+ const MeshType& mesh,
+ const Array<CellId>& list_of_cells)
+ {
+ return IntegrateOnCells<OutputType(const InputType)>::integrate(function_symbol_id, quadrature_descriptor, mesh,
+ Array<const CellId>{list_of_cells});
+ }
+};
+
+#endif // INTEGRATE_CELL_VALUE_HPP
diff --git a/src/language/utils/IntegrateOnCells.hpp b/src/language/utils/IntegrateOnCells.hpp
new file mode 100644
index 000000000..d4463da96
--- /dev/null
+++ b/src/language/utils/IntegrateOnCells.hpp
@@ -0,0 +1,606 @@
+#ifndef INTEGRATE_ON_CELLS_HPP
+#define INTEGRATE_ON_CELLS_HPP
+
+#include <analysis/IQuadratureDescriptor.hpp>
+#include <analysis/QuadratureManager.hpp>
+#include <geometry/AffineTransformation.hpp>
+#include <geometry/PrismTransformation.hpp>
+#include <geometry/PyramidTransformation.hpp>
+#include <geometry/Q1Transformation.hpp>
+#include <language/utils/PugsFunctionAdapter.hpp>
+#include <mesh/CellType.hpp>
+#include <mesh/Connectivity.hpp>
+#include <mesh/ItemId.hpp>
+#include <utils/Array.hpp>
+
+class FunctionSymbolId;
+
+template <typename T>
+class IntegrateOnCells;
+template <typename OutputType, typename InputType>
+class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<OutputType(InputType)>
+{
+ private:
+ using Adapter = PugsFunctionAdapter<OutputType(InputType)>;
+
+ template <size_t Dimension>
+ class AllCells
+ {
+ private:
+ const Connectivity<Dimension>& m_connectivity;
+
+ public:
+ using index_type = CellId;
+
+ PUGS_INLINE
+ CellId
+ cellId(const CellId cell_id) const
+ {
+ return cell_id;
+ }
+
+ PUGS_INLINE
+ size_t
+ size() const
+ {
+ return m_connectivity.numberOfCells();
+ }
+
+ PUGS_INLINE
+ AllCells(const Connectivity<Dimension>& connectivity) : m_connectivity{connectivity} {}
+ };
+
+ class CellList
+ {
+ private:
+ const Array<CellId>& m_cell_list;
+
+ public:
+ using index_type = Array<CellId>::index_type;
+
+ PUGS_INLINE
+ CellId
+ cellId(const index_type index) const
+ {
+ return m_cell_list[index];
+ }
+
+ PUGS_INLINE
+ size_t
+ size() const
+ {
+ return m_cell_list.size();
+ }
+
+ PUGS_INLINE
+ CellList(const Array<CellId>& cell_list) : m_cell_list{cell_list} {}
+ };
+
+ template <typename MeshType, typename OutputArrayT, typename ListTypeT>
+ static PUGS_INLINE void
+ _tensorialIntegrateTo(const FunctionSymbolId& function_symbol_id,
+ const IQuadratureDescriptor& quadrature_descriptor,
+ const MeshType& mesh,
+ const ListTypeT& cell_list,
+ OutputArrayT& value)
+ {
+ constexpr size_t Dimension = MeshType::Dimension;
+
+ static_assert(std::is_same_v<TinyVector<Dimension>, InputType>, "invalid input data type");
+ static_assert(std::is_same_v<std::remove_const_t<typename OutputArrayT::data_type>, OutputType>,
+ "invalid output data type");
+
+ auto& expression = Adapter::getFunctionExpression(function_symbol_id);
+ auto convert_result = Adapter::getResultConverter(expression.m_data_type);
+
+ auto context_list = Adapter::getContextList(expression);
+
+ using execution_space = typename Kokkos::DefaultExecutionSpace::execution_space;
+ Kokkos::Experimental::UniqueToken<execution_space, Kokkos::Experimental::UniqueTokenScope::Global> tokens;
+
+ if constexpr (std::is_arithmetic_v<OutputType>) {
+ value.fill(0);
+ } else if constexpr (is_tiny_vector_v<OutputType> or is_tiny_matrix_v<OutputType>) {
+ value.fill(zero);
+ } else {
+ static_assert(std::is_same_v<OutputType, double>, "unexpected output type");
+ }
+
+ const auto& connectivity = mesh.connectivity();
+
+ const auto cell_to_node_matrix = connectivity.cellToNodeMatrix();
+ const auto cell_type = connectivity.cellType();
+ const auto xr = mesh.xr();
+
+ auto invalid_cell = std::make_pair(false, CellId{});
+
+ const auto qf = [&] {
+ if constexpr (Dimension == 1) {
+ return QuadratureManager::instance().getLineFormula(quadrature_descriptor);
+ } else if constexpr (Dimension == 2) {
+ return QuadratureManager::instance().getSquareFormula(quadrature_descriptor);
+ } else {
+ static_assert(Dimension == 3);
+ return QuadratureManager::instance().getCubeFormula(quadrature_descriptor);
+ }
+ }();
+
+ parallel_for(cell_list.size(), [=, &expression, &tokens, &invalid_cell, &qf,
+ &cell_list](typename ListTypeT::index_type index) {
+ const int32_t t = tokens.acquire();
+ auto& execution_policy = context_list[t];
+
+ const CellId cell_id = cell_list.cellId(index);
+
+ const auto cell_to_node = cell_to_node_matrix[cell_id];
+
+ if constexpr (Dimension == 1) {
+ switch (cell_type[cell_id]) {
+ case CellType::Line: {
+ const AffineTransformation<1> T({xr[cell_to_node[0]], xr[cell_to_node[1]]});
+
+ for (size_t i_point = 0; i_point < qf.numberOfPoints(); ++i_point) {
+ const auto xi = qf.point(i_point);
+ Adapter::convertArgs(execution_policy.currentContext(), T(xi));
+ auto result = expression.execute(execution_policy);
+ value[index] += qf.weight(i_point) * T.jacobianDeterminant() * convert_result(std::move(result));
+ }
+ break;
+ }
+ default: {
+ invalid_cell = std::make_pair(true, cell_id);
+ break;
+ }
+ }
+ } else if constexpr (Dimension == 2) {
+ switch (cell_type[cell_id]) {
+ case CellType::Triangle: {
+ const Q1Transformation<2> T(
+ {xr[cell_to_node[0]], xr[cell_to_node[1]], xr[cell_to_node[2]], xr[cell_to_node[2]]});
+
+ for (size_t i_point = 0; i_point < qf.numberOfPoints(); ++i_point) {
+ const auto xi = qf.point(i_point);
+ Adapter::convertArgs(execution_policy.currentContext(), T(xi));
+ auto result = expression.execute(execution_policy);
+ value[index] += qf.weight(i_point) * T.jacobianDeterminant(xi) * convert_result(std::move(result));
+ }
+ break;
+ }
+ case CellType::Quadrangle: {
+ const Q1Transformation<2> T(
+ {xr[cell_to_node[0]], xr[cell_to_node[1]], xr[cell_to_node[2]], xr[cell_to_node[3]]});
+
+ for (size_t i_point = 0; i_point < qf.numberOfPoints(); ++i_point) {
+ const auto xi = qf.point(i_point);
+ Adapter::convertArgs(execution_policy.currentContext(), T(xi));
+ auto result = expression.execute(execution_policy);
+ value[index] += qf.weight(i_point) * T.jacobianDeterminant(xi) * convert_result(std::move(result));
+ }
+ break;
+ }
+ default: {
+ invalid_cell = std::make_pair(true, cell_id);
+ break;
+ }
+ }
+ } else {
+ static_assert(Dimension == 3);
+
+ switch (cell_type[cell_id]) {
+ case CellType::Tetrahedron: {
+ const Q1Transformation<3> T({xr[cell_to_node[0]], xr[cell_to_node[1]], xr[cell_to_node[2]],
+ xr[cell_to_node[2]], //
+ xr[cell_to_node[3]], xr[cell_to_node[3]], xr[cell_to_node[3]],
+ xr[cell_to_node[3]]});
+
+ for (size_t i_point = 0; i_point < qf.numberOfPoints(); ++i_point) {
+ const auto xi = qf.point(i_point);
+ Adapter::convertArgs(execution_policy.currentContext(), T(xi));
+ auto result = expression.execute(execution_policy);
+ value[index] += qf.weight(i_point) * T.jacobianDeterminant(xi) * convert_result(std::move(result));
+ }
+ break;
+ }
+ case CellType::Pyramid: {
+ if (cell_to_node.size() == 5) {
+ const Q1Transformation<3> T({xr[cell_to_node[0]], xr[cell_to_node[1]], xr[cell_to_node[2]], //
+ xr[cell_to_node[3]], //
+ xr[cell_to_node[4]], xr[cell_to_node[4]], xr[cell_to_node[4]],
+ xr[cell_to_node[4]]});
+
+ for (size_t i_point = 0; i_point < qf.numberOfPoints(); ++i_point) {
+ const auto xi = qf.point(i_point);
+ Adapter::convertArgs(execution_policy.currentContext(), T(xi));
+ auto result = expression.execute(execution_policy);
+ value[index] += qf.weight(i_point) * T.jacobianDeterminant(xi) * convert_result(std::move(result));
+ }
+ } else {
+ throw NotImplementedError("integration on pyramid with non-quadrangular base (number of nodes " +
+ std::to_string(cell_to_node.size()) + ")");
+ }
+ break;
+ }
+ case CellType::Diamond: {
+ if (cell_to_node.size() == 5) {
+ { // top tetrahedron
+ const Q1Transformation<3> T0({xr[cell_to_node[1]], xr[cell_to_node[2]], xr[cell_to_node[3]],
+ xr[cell_to_node[3]], //
+ xr[cell_to_node[4]], xr[cell_to_node[4]], xr[cell_to_node[4]],
+ xr[cell_to_node[4]]});
+
+ for (size_t i_point = 0; i_point < qf.numberOfPoints(); ++i_point) {
+ const auto xi = qf.point(i_point);
+ Adapter::convertArgs(execution_policy.currentContext(), T0(xi));
+ auto result = expression.execute(execution_policy);
+ value[index] += qf.weight(i_point) * T0.jacobianDeterminant(xi) * convert_result(std::move(result));
+ }
+ }
+ { // bottom tetrahedron
+ const Q1Transformation<3> T1({xr[cell_to_node[3]], xr[cell_to_node[3]], xr[cell_to_node[2]],
+ xr[cell_to_node[1]], //
+ xr[cell_to_node[0]], xr[cell_to_node[0]], xr[cell_to_node[0]],
+ xr[cell_to_node[0]]});
+
+ for (size_t i_point = 0; i_point < qf.numberOfPoints(); ++i_point) {
+ const auto xi = qf.point(i_point);
+ Adapter::convertArgs(execution_policy.currentContext(), T1(xi));
+ auto result = expression.execute(execution_policy);
+ value[index] += qf.weight(i_point) * T1.jacobianDeterminant(xi) * convert_result(std::move(result));
+ }
+ }
+ } else if (cell_to_node.size() == 6) {
+ { // top pyramid
+ const Q1Transformation<3> T0({xr[cell_to_node[1]], xr[cell_to_node[2]], xr[cell_to_node[3]],
+ xr[cell_to_node[4]], //
+ xr[cell_to_node[5]], xr[cell_to_node[5]], xr[cell_to_node[5]],
+ xr[cell_to_node[5]]});
+
+ for (size_t i_point = 0; i_point < qf.numberOfPoints(); ++i_point) {
+ const auto xi = qf.point(i_point);
+ Adapter::convertArgs(execution_policy.currentContext(), T0(xi));
+ auto result = expression.execute(execution_policy);
+ value[index] += qf.weight(i_point) * T0.jacobianDeterminant(xi) * convert_result(std::move(result));
+ }
+ }
+ { // bottom pyramid
+ const Q1Transformation<3> T1({xr[cell_to_node[4]], xr[cell_to_node[3]], xr[cell_to_node[2]],
+ xr[cell_to_node[1]], //
+ xr[cell_to_node[0]], xr[cell_to_node[0]], xr[cell_to_node[0]],
+ xr[cell_to_node[0]]});
+
+ for (size_t i_point = 0; i_point < qf.numberOfPoints(); ++i_point) {
+ const auto xi = qf.point(i_point);
+ Adapter::convertArgs(execution_policy.currentContext(), T1(xi));
+ auto result = expression.execute(execution_policy);
+ value[index] += qf.weight(i_point) * T1.jacobianDeterminant(xi) * convert_result(std::move(result));
+ }
+ }
+ } else {
+ throw NotImplementedError("integration on diamond with non-quadrangular base (" +
+ std::to_string(cell_to_node.size()) + ")");
+ }
+ break;
+ }
+ case CellType::Prism: {
+ const Q1Transformation<3> T({xr[cell_to_node[0]], xr[cell_to_node[1]], //
+ xr[cell_to_node[2]], xr[cell_to_node[2]], //
+ xr[cell_to_node[3]], xr[cell_to_node[4]], //
+ xr[cell_to_node[5]], xr[cell_to_node[5]]});
+
+ for (size_t i_point = 0; i_point < qf.numberOfPoints(); ++i_point) {
+ const auto xi = qf.point(i_point);
+ Adapter::convertArgs(execution_policy.currentContext(), T(xi));
+ auto result = expression.execute(execution_policy);
+ value[index] += qf.weight(i_point) * T.jacobianDeterminant(xi) * convert_result(std::move(result));
+ }
+ break;
+ }
+ case CellType::Hexahedron: {
+ const Q1Transformation<3> T({xr[cell_to_node[0]], xr[cell_to_node[1]], xr[cell_to_node[2]],
+ xr[cell_to_node[3]], xr[cell_to_node[4]], xr[cell_to_node[5]],
+ xr[cell_to_node[6]], xr[cell_to_node[7]]});
+
+ for (size_t i_point = 0; i_point < qf.numberOfPoints(); ++i_point) {
+ const auto xi = qf.point(i_point);
+ Adapter::convertArgs(execution_policy.currentContext(), T(xi));
+ auto result = expression.execute(execution_policy);
+ value[index] += qf.weight(i_point) * T.jacobianDeterminant(xi) * convert_result(std::move(result));
+ }
+ break;
+ }
+ default: {
+ invalid_cell = std::make_pair(true, cell_id);
+ break;
+ }
+ }
+ }
+
+ tokens.release(t);
+ });
+
+ if (invalid_cell.first) {
+ std::ostringstream os;
+ os << "invalid cell type for integration: " << name(cell_type[invalid_cell.second]);
+ throw UnexpectedError(os.str());
+ }
+ }
+
+ template <typename MeshType, typename OutputArrayT, typename ListTypeT>
+ static PUGS_INLINE void
+ _directIntegrateTo(const FunctionSymbolId& function_symbol_id,
+ const IQuadratureDescriptor& quadrature_descriptor,
+ const MeshType& mesh,
+ const ListTypeT& cell_list,
+ OutputArrayT& value)
+ {
+ Assert(not quadrature_descriptor.isTensorial());
+
+ constexpr size_t Dimension = MeshType::Dimension;
+
+ static_assert(std::is_same_v<TinyVector<Dimension>, InputType>, "invalid input data type");
+ static_assert(std::is_same_v<std::remove_const_t<typename OutputArrayT::data_type>, OutputType>,
+ "invalid output data type");
+
+ auto& expression = Adapter::getFunctionExpression(function_symbol_id);
+ auto convert_result = Adapter::getResultConverter(expression.m_data_type);
+
+ auto context_list = Adapter::getContextList(expression);
+
+ using execution_space = typename Kokkos::DefaultExecutionSpace::execution_space;
+ Kokkos::Experimental::UniqueToken<execution_space, Kokkos::Experimental::UniqueTokenScope::Global> tokens;
+
+ if constexpr (std::is_arithmetic_v<OutputType>) {
+ value.fill(0);
+ } else if constexpr (is_tiny_vector_v<OutputType> or is_tiny_matrix_v<OutputType>) {
+ value.fill(zero);
+ } else {
+ static_assert(std::is_same_v<OutputType, double>, "unexpected output type");
+ }
+
+ const auto& connectivity = mesh.connectivity();
+
+ const auto cell_to_node_matrix = connectivity.cellToNodeMatrix();
+ const auto cell_type = connectivity.cellType();
+ const auto xr = mesh.xr();
+
+ auto invalid_cell = std::make_pair(false, CellId{});
+
+ parallel_for(cell_list.size(), [=, &expression, &tokens, &invalid_cell, &quadrature_descriptor,
+ &cell_list](const typename ListTypeT::index_type& index) {
+ const int32_t t = tokens.acquire();
+ auto& execution_policy = context_list[t];
+
+ const CellId cell_id = cell_list.cellId(index);
+
+ const auto cell_to_node = cell_to_node_matrix[cell_id];
+
+ if constexpr (Dimension == 1) {
+ switch (cell_type[cell_id]) {
+ case CellType::Line: {
+ const AffineTransformation<1> T({xr[cell_to_node[0]], xr[cell_to_node[1]]});
+ const auto qf = QuadratureManager::instance().getLineFormula(quadrature_descriptor);
+
+ for (size_t i_point = 0; i_point < qf.numberOfPoints(); ++i_point) {
+ const auto xi = qf.point(i_point);
+ Adapter::convertArgs(execution_policy.currentContext(), T(xi));
+ auto result = expression.execute(execution_policy);
+ value[index] += qf.weight(i_point) * T.jacobianDeterminant() * convert_result(std::move(result));
+ }
+ break;
+ }
+ default: {
+ invalid_cell = std::make_pair(true, cell_id);
+ break;
+ }
+ }
+ } else if constexpr (Dimension == 2) {
+ switch (cell_type[cell_id]) {
+ case CellType::Triangle: {
+ const AffineTransformation<2> T({xr[cell_to_node[0]], xr[cell_to_node[1]], xr[cell_to_node[2]]});
+ const auto qf = QuadratureManager::instance().getTriangleFormula(quadrature_descriptor);
+
+ for (size_t i_point = 0; i_point < qf.numberOfPoints(); ++i_point) {
+ const auto xi = qf.point(i_point);
+ Adapter::convertArgs(execution_policy.currentContext(), T(xi));
+ auto result = expression.execute(execution_policy);
+ value[index] += qf.weight(i_point) * T.jacobianDeterminant() * convert_result(std::move(result));
+ }
+ break;
+ }
+ case CellType::Quadrangle: {
+ const Q1Transformation<2> T(
+ {xr[cell_to_node[0]], xr[cell_to_node[1]], xr[cell_to_node[2]], xr[cell_to_node[3]]});
+ const auto qf = QuadratureManager::instance().getSquareFormula(quadrature_descriptor);
+
+ for (size_t i_point = 0; i_point < qf.numberOfPoints(); ++i_point) {
+ const auto xi = qf.point(i_point);
+ Adapter::convertArgs(execution_policy.currentContext(), T(xi));
+ auto result = expression.execute(execution_policy);
+ value[index] += qf.weight(i_point) * T.jacobianDeterminant(xi) * convert_result(std::move(result));
+ }
+ break;
+ }
+ default: {
+ invalid_cell = std::make_pair(true, cell_id);
+ break;
+ }
+ }
+ } else {
+ static_assert(Dimension == 3);
+
+ switch (cell_type[cell_id]) {
+ case CellType::Tetrahedron: {
+ const AffineTransformation<3> T({xr[cell_to_node[0]], xr[cell_to_node[1]], //
+ xr[cell_to_node[2]], xr[cell_to_node[3]]});
+ const auto qf = QuadratureManager::instance().getTetrahedronFormula(quadrature_descriptor);
+
+ for (size_t i_point = 0; i_point < qf.numberOfPoints(); ++i_point) {
+ const auto xi = qf.point(i_point);
+ Adapter::convertArgs(execution_policy.currentContext(), T(xi));
+ auto result = expression.execute(execution_policy);
+ value[index] += qf.weight(i_point) * T.jacobianDeterminant() * convert_result(std::move(result));
+ }
+ break;
+ }
+ case CellType::Pyramid: {
+ if (cell_to_node.size() == 5) {
+ const PyramidTransformation T({xr[cell_to_node[0]], xr[cell_to_node[1]], xr[cell_to_node[2]], //
+ xr[cell_to_node[3]], xr[cell_to_node[4]]});
+ const auto qf = QuadratureManager::instance().getPyramidFormula(quadrature_descriptor);
+
+ for (size_t i_point = 0; i_point < qf.numberOfPoints(); ++i_point) {
+ const auto xi = qf.point(i_point);
+ Adapter::convertArgs(execution_policy.currentContext(), T(xi));
+ auto result = expression.execute(execution_policy);
+ value[index] += qf.weight(i_point) * T.jacobianDeterminant(xi) * convert_result(std::move(result));
+ }
+ } else {
+ throw NotImplementedError("integration on pyramid with non-quadrangular base");
+ }
+ break;
+ }
+ case CellType::Diamond: {
+ if (cell_to_node.size() == 5) {
+ const auto qf = QuadratureManager::instance().getTetrahedronFormula(quadrature_descriptor);
+ { // top tetrahedron
+ const AffineTransformation<3> T0({xr[cell_to_node[1]], xr[cell_to_node[2]], xr[cell_to_node[3]], //
+ xr[cell_to_node[4]]});
+
+ for (size_t i_point = 0; i_point < qf.numberOfPoints(); ++i_point) {
+ const auto xi = qf.point(i_point);
+ Adapter::convertArgs(execution_policy.currentContext(), T0(xi));
+ auto result = expression.execute(execution_policy);
+ value[index] += qf.weight(i_point) * T0.jacobianDeterminant() * convert_result(std::move(result));
+ }
+ }
+ { // bottom tetrahedron
+ const AffineTransformation<3> T1({xr[cell_to_node[3]], xr[cell_to_node[2]], xr[cell_to_node[1]], //
+ xr[cell_to_node[0]]});
+
+ for (size_t i_point = 0; i_point < qf.numberOfPoints(); ++i_point) {
+ const auto xi = qf.point(i_point);
+ Adapter::convertArgs(execution_policy.currentContext(), T1(xi));
+ auto result = expression.execute(execution_policy);
+ value[index] += qf.weight(i_point) * T1.jacobianDeterminant() * convert_result(std::move(result));
+ }
+ }
+ } else if (cell_to_node.size() == 6) {
+ const auto qf = QuadratureManager::instance().getCubeFormula(quadrature_descriptor);
+ { // top pyramid
+ const Q1Transformation<3> T0({xr[cell_to_node[1]], xr[cell_to_node[2]], xr[cell_to_node[3]],
+ xr[cell_to_node[4]], //
+ xr[cell_to_node[5]], xr[cell_to_node[5]], xr[cell_to_node[5]],
+ xr[cell_to_node[5]]});
+
+ for (size_t i_point = 0; i_point < qf.numberOfPoints(); ++i_point) {
+ const auto xi = qf.point(i_point);
+ Adapter::convertArgs(execution_policy.currentContext(), T0(xi));
+ auto result = expression.execute(execution_policy);
+ value[index] += qf.weight(i_point) * T0.jacobianDeterminant(xi) * convert_result(std::move(result));
+ }
+ }
+ { // bottom pyramid
+ const Q1Transformation<3> T1({xr[cell_to_node[4]], xr[cell_to_node[3]], xr[cell_to_node[2]],
+ xr[cell_to_node[1]], //
+ xr[cell_to_node[0]], xr[cell_to_node[0]], xr[cell_to_node[0]],
+ xr[cell_to_node[0]]});
+
+ for (size_t i_point = 0; i_point < qf.numberOfPoints(); ++i_point) {
+ const auto xi = qf.point(i_point);
+ Adapter::convertArgs(execution_policy.currentContext(), T1(xi));
+ auto result = expression.execute(execution_policy);
+ value[index] += qf.weight(i_point) * T1.jacobianDeterminant(xi) * convert_result(std::move(result));
+ }
+ }
+ } else {
+ throw NotImplementedError("integration on diamond with non-quadrangular base (" +
+ std::to_string(cell_to_node.size()) + ")");
+ }
+ break;
+ }
+ case CellType::Prism: {
+ const PrismTransformation T({xr[cell_to_node[0]], xr[cell_to_node[1]], xr[cell_to_node[2]],
+ xr[cell_to_node[3]], xr[cell_to_node[4]], xr[cell_to_node[5]]});
+ const auto qf = QuadratureManager::instance().getPrismFormula(quadrature_descriptor);
+
+ for (size_t i_point = 0; i_point < qf.numberOfPoints(); ++i_point) {
+ const auto xi = qf.point(i_point);
+ Adapter::convertArgs(execution_policy.currentContext(), T(xi));
+ auto result = expression.execute(execution_policy);
+ value[index] += qf.weight(i_point) * T.jacobianDeterminant(xi) * convert_result(std::move(result));
+ }
+ break;
+ }
+ case CellType::Hexahedron: {
+ const Q1Transformation<3> T({xr[cell_to_node[0]], xr[cell_to_node[1]], xr[cell_to_node[2]],
+ xr[cell_to_node[3]], xr[cell_to_node[4]], xr[cell_to_node[5]],
+ xr[cell_to_node[6]], xr[cell_to_node[7]]});
+ const auto qf = QuadratureManager::instance().getCubeFormula(quadrature_descriptor);
+
+ for (size_t i_point = 0; i_point < qf.numberOfPoints(); ++i_point) {
+ const auto xi = qf.point(i_point);
+ Adapter::convertArgs(execution_policy.currentContext(), T(xi));
+ auto result = expression.execute(execution_policy);
+ value[index] += qf.weight(i_point) * T.jacobianDeterminant(xi) * convert_result(std::move(result));
+ }
+ break;
+ }
+ default: {
+ invalid_cell = std::make_pair(true, cell_id);
+ break;
+ }
+ }
+ }
+
+ tokens.release(t);
+ });
+
+ if (invalid_cell.first) {
+ std::ostringstream os;
+ os << "invalid cell type for integration: " << name(cell_type[invalid_cell.second]);
+ throw UnexpectedError(os.str());
+ }
+ }
+
+ public:
+ template <typename MeshType, typename OutputArrayT>
+ static PUGS_INLINE void
+ integrateTo(const FunctionSymbolId& function_symbol_id,
+ const IQuadratureDescriptor& quadrature_descriptor,
+ const MeshType& mesh,
+ OutputArrayT& value)
+ {
+ constexpr size_t Dimension = MeshType::Dimension;
+
+ if (quadrature_descriptor.isTensorial()) {
+ _tensorialIntegrateTo<MeshType, OutputArrayT>(function_symbol_id, quadrature_descriptor, mesh,
+ AllCells<Dimension>{mesh.connectivity()}, value);
+ } else {
+ _directIntegrateTo<MeshType, OutputArrayT>(function_symbol_id, quadrature_descriptor, mesh,
+ AllCells<Dimension>{mesh.connectivity()}, value);
+ }
+ }
+
+ template <typename MeshType, typename OutputArrayT>
+ static PUGS_INLINE Array<OutputType>
+ integrate(const FunctionSymbolId& function_symbol_id,
+ const IQuadratureDescriptor& quadrature_descriptor,
+ const MeshType& mesh,
+ const Array<CellId>& cell_list)
+ {
+ constexpr size_t Dimension = MeshType::Dimension;
+
+ Array<OutputType> value(size(cell_list));
+ if (quadrature_descriptor.isTensorial()) {
+ _tensorialIntegrateTo<MeshType, OutputArrayT>(function_symbol_id, quadrature_descriptor, mesh,
+ CellList{cell_list}, value);
+ } else {
+ _directIntegrateTo<MeshType, OutputArrayT>(function_symbol_id, quadrature_descriptor, mesh, CellList{cell_list},
+ value);
+ }
+
+ return value;
+ }
+};
+
+#endif // INTEGRATE_ON_CELLS_HPP
diff --git a/src/scheme/CMakeLists.txt b/src/scheme/CMakeLists.txt
index 2131f2fec..115d7bc08 100644
--- a/src/scheme/CMakeLists.txt
+++ b/src/scheme/CMakeLists.txt
@@ -3,6 +3,8 @@
add_library(
PugsScheme
AcousticSolver.cpp
+ DiscreteFunctionIntegrator.cpp
DiscreteFunctionInterpoler.cpp
DiscreteFunctionUtils.cpp
+ DiscreteFunctionVectorIntegrator.cpp
DiscreteFunctionVectorInterpoler.cpp)
diff --git a/src/scheme/DiscreteFunctionIntegrator.cpp b/src/scheme/DiscreteFunctionIntegrator.cpp
new file mode 100644
index 000000000..91cf921df
--- /dev/null
+++ b/src/scheme/DiscreteFunctionIntegrator.cpp
@@ -0,0 +1,132 @@
+#include <scheme/DiscreteFunctionIntegrator.hpp>
+
+#include <language/utils/IntegrateCellValue.hpp>
+#include <scheme/DiscreteFunctionP0.hpp>
+#include <utils/Exceptions.hpp>
+
+template <size_t Dimension, typename DataType, typename ValueType>
+std::shared_ptr<IDiscreteFunction>
+DiscreteFunctionIntegrator::_integrate() const
+{
+ using MeshType = Mesh<Connectivity<Dimension>>;
+ std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(m_mesh);
+
+ if constexpr (std::is_same_v<DataType, ValueType>) {
+ return std::make_shared<
+ DiscreteFunctionP0<Dimension, ValueType>>(mesh,
+ IntegrateCellValue<DataType(TinyVector<Dimension>)>::template integrate<
+ MeshType>(m_function_id, *m_quadrature_descriptor, *mesh));
+ } else {
+ static_assert(std::is_convertible_v<DataType, ValueType>);
+
+ CellValue<DataType> cell_data =
+ IntegrateCellValue<DataType(TinyVector<Dimension>)>::template integrate<MeshType>(m_function_id,
+ *m_quadrature_descriptor,
+ *mesh);
+
+ CellValue<ValueType> cell_value{mesh->connectivity()};
+
+ parallel_for(
+ mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { cell_value[cell_id] = cell_data[cell_id]; });
+
+ return std::make_shared<DiscreteFunctionP0<Dimension, ValueType>>(mesh, cell_value);
+ }
+}
+
+template <size_t Dimension>
+std::shared_ptr<IDiscreteFunction>
+DiscreteFunctionIntegrator::_integrate() const
+{
+ const auto& function_descriptor = m_function_id.descriptor();
+ Assert(function_descriptor.domainMappingNode().children[1]->m_data_type == ASTNodeDataType::typename_t);
+
+ const ASTNodeDataType& data_type = function_descriptor.domainMappingNode().children[1]->m_data_type.contentType();
+
+ switch (data_type) {
+ case ASTNodeDataType::bool_t: {
+ return this->_integrate<Dimension, bool, double>();
+ }
+ case ASTNodeDataType::unsigned_int_t: {
+ return this->_integrate<Dimension, uint64_t, double>();
+ }
+ case ASTNodeDataType::int_t: {
+ return this->_integrate<Dimension, int64_t, double>();
+ }
+ case ASTNodeDataType::double_t: {
+ return this->_integrate<Dimension, double>();
+ }
+ case ASTNodeDataType::vector_t: {
+ switch (data_type.dimension()) {
+ case 1: {
+ return this->_integrate<Dimension, TinyVector<1>>();
+ }
+ case 2: {
+ return this->_integrate<Dimension, TinyVector<2>>();
+ }
+ case 3: {
+ return this->_integrate<Dimension, TinyVector<3>>();
+ }
+ // LCOV_EXCL_START
+ default: {
+ std::ostringstream os;
+ os << "invalid vector dimension " << rang::fgB::red << data_type.dimension() << rang::style::reset;
+
+ throw UnexpectedError(os.str());
+ }
+ // LCOV_EXCL_STOP
+ }
+ }
+ case ASTNodeDataType::matrix_t: {
+ Assert(data_type.numberOfColumns() == data_type.numberOfRows(), "undefined matrix type");
+ switch (data_type.numberOfColumns()) {
+ case 1: {
+ return this->_integrate<Dimension, TinyMatrix<1>>();
+ }
+ case 2: {
+ return this->_integrate<Dimension, TinyMatrix<2>>();
+ }
+ case 3: {
+ return this->_integrate<Dimension, TinyMatrix<3>>();
+ }
+ // LCOV_EXCL_START
+ default: {
+ std::ostringstream os;
+ os << "invalid vector dimension " << rang::fgB::red << data_type.dimension() << rang::style::reset;
+
+ throw UnexpectedError(os.str());
+ }
+ // LCOV_EXCL_STOP
+ }
+ }
+ // LCOV_EXCL_START
+ default: {
+ std::ostringstream os;
+ os << "invalid integrate value type: " << rang::fgB::red << dataTypeName(data_type) << rang::style::reset;
+
+ throw UnexpectedError(os.str());
+ }
+ // LCOV_EXCL_STOP
+ }
+}
+
+std::shared_ptr<IDiscreteFunction>
+DiscreteFunctionIntegrator::integrate() const
+{
+ std::shared_ptr<IDiscreteFunction> discrete_function;
+ switch (m_mesh->dimension()) {
+ case 1: {
+ return this->_integrate<1>();
+ }
+ case 2: {
+ return this->_integrate<2>();
+ }
+ case 3: {
+ return this->_integrate<3>();
+ }
+ // LCOV_EXCL_START
+ default: {
+ throw UnexpectedError("invalid dimension");
+ }
+ // LCOV_EXCL_STOP
+ }
+}
diff --git a/src/scheme/DiscreteFunctionIntegrator.hpp b/src/scheme/DiscreteFunctionIntegrator.hpp
new file mode 100644
index 000000000..1784f891e
--- /dev/null
+++ b/src/scheme/DiscreteFunctionIntegrator.hpp
@@ -0,0 +1,39 @@
+#ifndef DISCRETE_FUNCTION_INTEGRATOR_HPP
+#define DISCRETE_FUNCTION_INTEGRATOR_HPP
+
+#include <analysis/IQuadratureDescriptor.hpp>
+#include <language/utils/FunctionSymbolId.hpp>
+#include <mesh/IMesh.hpp>
+#include <scheme/IDiscreteFunction.hpp>
+
+#include <memory>
+
+class DiscreteFunctionIntegrator
+{
+ private:
+ std::shared_ptr<const IMesh> m_mesh;
+ std::shared_ptr<const IQuadratureDescriptor> m_quadrature_descriptor;
+ const FunctionSymbolId m_function_id;
+
+ template <size_t Dimension, typename DataType, typename ValueType = DataType>
+ std::shared_ptr<IDiscreteFunction> _integrate() const;
+
+ template <size_t Dimension>
+ std::shared_ptr<IDiscreteFunction> _integrate() const;
+
+ public:
+ std::shared_ptr<IDiscreteFunction> integrate() const;
+
+ DiscreteFunctionIntegrator(const std::shared_ptr<const IMesh>& mesh,
+ const std::shared_ptr<const IQuadratureDescriptor>& quadrature_descriptor,
+ const FunctionSymbolId& function_id)
+ : m_mesh{mesh}, m_quadrature_descriptor{quadrature_descriptor}, m_function_id{function_id}
+ {}
+
+ DiscreteFunctionIntegrator(const DiscreteFunctionIntegrator&) = delete;
+ DiscreteFunctionIntegrator(DiscreteFunctionIntegrator&&) = delete;
+
+ ~DiscreteFunctionIntegrator() = default;
+};
+
+#endif // DISCRETE_FUNCTION_INTEGRATOR_HPP
diff --git a/src/scheme/DiscreteFunctionVectorIntegrator.cpp b/src/scheme/DiscreteFunctionVectorIntegrator.cpp
new file mode 100644
index 000000000..6d8937f84
--- /dev/null
+++ b/src/scheme/DiscreteFunctionVectorIntegrator.cpp
@@ -0,0 +1,70 @@
+#include <scheme/DiscreteFunctionVectorIntegrator.hpp>
+
+#include <language/utils/IntegrateCellArray.hpp>
+
+#include <scheme/DiscreteFunctionP0Vector.hpp>
+#include <utils/Exceptions.hpp>
+
+template <size_t Dimension, typename DataType>
+std::shared_ptr<IDiscreteFunction>
+DiscreteFunctionVectorIntegrator::_integrate() const
+{
+ std::shared_ptr mesh = std::dynamic_pointer_cast<const Mesh<Connectivity<Dimension>>>(m_mesh);
+
+ return std::make_shared<
+ DiscreteFunctionP0Vector<Dimension, DataType>>(mesh, IntegrateCellArray<DataType(TinyVector<Dimension>)>::
+ template integrate(m_function_id_list,
+ *m_quadrature_descriptor, *mesh));
+}
+
+template <size_t Dimension>
+std::shared_ptr<IDiscreteFunction>
+DiscreteFunctionVectorIntegrator::_integrate() const
+{
+ for (const auto& function_id : m_function_id_list) {
+ const auto& function_descriptor = function_id.descriptor();
+ Assert(function_descriptor.domainMappingNode().children[1]->m_data_type == ASTNodeDataType::typename_t);
+ const ASTNodeDataType& data_type = function_descriptor.domainMappingNode().children[1]->m_data_type.contentType();
+
+ switch (data_type) {
+ case ASTNodeDataType::bool_t:
+ case ASTNodeDataType::unsigned_int_t:
+ case ASTNodeDataType::int_t:
+ case ASTNodeDataType::double_t: {
+ break;
+ }
+ default: {
+ std::ostringstream os;
+ os << "vector functions require scalar value type.\n"
+ << "Invalid interpolation value type: " << rang::fgB::red << dataTypeName(data_type) << rang::style::reset;
+ throw NormalError(os.str());
+ }
+ }
+ }
+ return this->_integrate<Dimension, double>();
+}
+
+std::shared_ptr<IDiscreteFunction>
+DiscreteFunctionVectorIntegrator::integrate() const
+{
+ if (m_discrete_function_descriptor->type() != DiscreteFunctionType::P0Vector) {
+ throw NormalError("invalid discrete function type for vector interpolation");
+ }
+
+ switch (m_mesh->dimension()) {
+ case 1: {
+ return this->_integrate<1>();
+ }
+ case 2: {
+ return this->_integrate<2>();
+ }
+ case 3: {
+ return this->_integrate<3>();
+ }
+ // LCOV_EXCL_START
+ default: {
+ throw UnexpectedError("invalid dimension");
+ }
+ // LCOV_EXCL_STOP
+ }
+}
diff --git a/src/scheme/DiscreteFunctionVectorIntegrator.hpp b/src/scheme/DiscreteFunctionVectorIntegrator.hpp
new file mode 100644
index 000000000..01e086e7e
--- /dev/null
+++ b/src/scheme/DiscreteFunctionVectorIntegrator.hpp
@@ -0,0 +1,47 @@
+#ifndef DISCRETE_FUNCTION_VECTOR_INTEGRATOR_HPP
+#define DISCRETE_FUNCTION_VECTOR_INTEGRATOR_HPP
+
+#include <analysis/IQuadratureDescriptor.hpp>
+#include <language/utils/FunctionSymbolId.hpp>
+#include <mesh/IMesh.hpp>
+#include <scheme/IDiscreteFunction.hpp>
+#include <scheme/IDiscreteFunctionDescriptor.hpp>
+
+#include <memory>
+#include <vector>
+
+class DiscreteFunctionVectorIntegrator
+{
+ private:
+ std::shared_ptr<const IMesh> m_mesh;
+ std::shared_ptr<const IQuadratureDescriptor> m_quadrature_descriptor;
+ std::shared_ptr<const IDiscreteFunctionDescriptor> m_discrete_function_descriptor;
+ const std::vector<FunctionSymbolId> m_function_id_list;
+
+ template <size_t Dimension, typename DataType>
+ std::shared_ptr<IDiscreteFunction> _integrate() const;
+
+ template <size_t Dimension>
+ std::shared_ptr<IDiscreteFunction> _integrate() const;
+
+ public:
+ std::shared_ptr<IDiscreteFunction> integrate() const;
+
+ DiscreteFunctionVectorIntegrator(
+ const std::shared_ptr<const IMesh>& mesh,
+ const std::shared_ptr<const IQuadratureDescriptor>& quadrature_descriptor,
+ const std::shared_ptr<const IDiscreteFunctionDescriptor>& discrete_function_descriptor,
+ const std::vector<FunctionSymbolId>& function_id_list)
+ : m_mesh{mesh},
+ m_quadrature_descriptor(quadrature_descriptor),
+ m_discrete_function_descriptor{discrete_function_descriptor},
+ m_function_id_list{function_id_list}
+ {}
+
+ DiscreteFunctionVectorIntegrator(const DiscreteFunctionVectorIntegrator&) = delete;
+ DiscreteFunctionVectorIntegrator(DiscreteFunctionVectorIntegrator&&) = delete;
+
+ ~DiscreteFunctionVectorIntegrator() = default;
+};
+
+#endif // DISCRETE_FUNCTION_VECTOR_INTEGRATOR_HPP
--
GitLab