diff --git a/src/algebra/IntegrationTools.hpp b/src/algebra/IntegrationTools.hpp
deleted file mode 100644
index 77e79851e0a8000807f6901373c8801c26130bb2..0000000000000000000000000000000000000000
--- a/src/algebra/IntegrationTools.hpp
+++ /dev/null
@@ -1,197 +0,0 @@
-#ifndef INTEGRATION_TOOLS_HPP
-#define INTEGRATION_TOOLS_HPP
-
-#include <language/utils/EvaluateAtPoints.hpp>
-#include <utils/Exceptions.hpp>
-#include <utils/PugsAssert.hpp>
-#include <utils/PugsMacros.hpp>
-#include <utils/Types.hpp>
-
-#include <analysis/QuadratureManager.hpp>
-#include <analysis/QuadratureType.hpp>
-#include <analysis/TensorialGaussLegendreQuadrature.hpp>
-#include <analysis/TensorialGaussLobattoQuadrature.hpp>
-
-#include <cmath>
-
-template <size_t Order, typename T = double>
-class IntegrationTools
-{
- public:
- using data_type = T;
- SmallArray<const double> m_weights;
- SmallArray<const TinyVector<1>> m_quadrature_positions;
-
- public:
- PUGS_INLINE T
- integrate(const FunctionSymbolId& function, const double& a, const double& b) const
- {
- T result = 0;
- SmallArray<const TinyVector<1>> positions = quadraturePoints(a, b);
- SmallArray<const double> weights = m_weights;
- SmallArray<T> values(positions.size());
- EvaluateAtPoints<double(const TinyVector<1>)>::evaluateTo(function, positions, values);
- Assert(positions.size() == weights.size(), "Wrong number of quadrature points and/or weights in Gauss quadrature");
- for (size_t i = 0; i < values.size(); ++i) {
- result += weights[i] * values[i];
- }
- return (b - a) / 2 * result;
- }
-
- PUGS_INLINE SmallArray<T>
- integrateFunction(const FunctionSymbolId& function, const SmallArray<TinyVector<1>>& vertices) const
- {
- SmallArray<const TinyVector<1>> positions = quadraturePoints(vertices);
- SmallArray<T> result{vertices.size() - 1};
- SmallArray<double> interval_size{vertices.size() - 1};
- SmallArray<const double> weights = m_weights;
- SmallArray<T> values(positions.size());
- EvaluateAtPoints<double(const TinyVector<1>)>::evaluateTo(function, positions, values);
-
- for (size_t i = 0; i < interval_size.size(); ++i) {
- interval_size[i] = vertices[i + 1][0] - vertices[i][0];
- }
- for (size_t i = 0; i < interval_size.size(); ++i) {
- for (size_t j = 0; j < weights.size(); j++) {
- result[i] += 0.5 * interval_size[i] * weights[j] * values[i * weights.size() + j];
- }
- }
- return result;
- }
-
- PUGS_INLINE T
- testIntegrate(const double& a, const double& b) const
- {
- T result = 0;
- SmallArray<const TinyVector<1>> positions = quadraturePoints(a, b);
-
- SmallArray<const double> weights = m_weights;
-
- SmallArray<T> values(weights.size());
- for (size_t j = 0; j < weights.size(); ++j) {
- values[j] = std::pow(positions[j][0], 3.);
- }
-
- Assert(positions.size() == weights.size(), "Wrong number of quadrature points and/or weights in Gauss quadrature");
- for (size_t i = 0; i < values.size(); ++i) {
- result += weights[i] * values[i];
- }
-
- return 0.5 * (b - a) * result;
- }
-
- PUGS_INLINE SmallArray<T>
- testIntegrateFunction(const SmallArray<TinyVector<1>>& vertices) const
- {
- SmallArray<const TinyVector<1>> positions = quadraturePoints(vertices);
- SmallArray<double> interval_size{vertices.size() - 1};
- SmallArray<const double> weights = this->weights();
- SmallArray<T> values(positions.size());
- for (size_t j = 0; j < positions.size(); ++j) {
- values[j] = std::pow(positions[j][0], 3);
- }
-
- for (size_t i = 0; i < interval_size.size(); ++i) {
- interval_size[i] = vertices[i + 1][0] - vertices[i][0];
- }
- SmallArray<T> result{vertices.size() - 1};
- if constexpr (std::is_arithmetic_v<T>) {
- result.fill(0);
- } else if constexpr (is_tiny_vector_v<T> or is_tiny_matrix_v<T>) {
- result.fill(zero);
- } else {
- static_assert(std::is_arithmetic_v<T>, "incorrect data type");
- }
- for (size_t i = 0; i < interval_size.size(); ++i) {
- for (size_t j = 0; j < weights.size(); j++) {
- result[i] += 0.5 * interval_size[i] * weights[j] * values[i * weights.size() + j];
- }
- }
- return result;
- }
-
- PUGS_INLINE SmallArray<const TinyVector<1>>
- quadraturePoints(const double& a, const double& b) const
- {
- SmallArray<TinyVector<1>> points{m_quadrature_positions.size()};
- double length = b - a;
- for (size_t i = 0; i < m_quadrature_positions.size(); i++) {
- TinyVector<1> qpos{m_quadrature_positions[i]};
- points[i] = 0.5 * (length * qpos + a + b);
- }
- return points;
- }
-
- PUGS_INLINE
- SmallArray<TinyVector<1>>
- quadraturePoints(const SmallArray<TinyVector<1>>& positions) const
- {
- size_t number_of_intervals = positions.size() - 1;
- size_t quadrature_size = m_quadrature_positions.size();
- SmallArray<TinyVector<1>> quadratures{quadrature_size * number_of_intervals};
- for (size_t j = 0; j < number_of_intervals; ++j) {
- double a = positions[j][0];
- double b = positions[j + 1][0];
- auto intermediate_positions = this->quadraturePoints(a, b);
- for (size_t k = 0; k < quadrature_size; k++) {
- quadratures[j * quadrature_size + k] = intermediate_positions[k];
- }
- }
- return quadratures;
- }
-
- PUGS_INLINE SmallArray<const double>
- weights() const
- {
- return m_weights;
- }
-
- PUGS_INLINE
- constexpr size_t
- order() const
- {
- return Order;
- }
-
- PUGS_INLINE
- IntegrationTools(QuadratureType quadrature)
- {
- switch (quadrature) {
- case QuadratureType::GaussLobatto: {
- auto tensorial_gauss_lobatto = QuadratureManager::instance().getLineGaussLobattoFormula(Order);
-
- m_quadrature_positions = tensorial_gauss_lobatto.pointList();
- m_weights = tensorial_gauss_lobatto.weightList();
- break;
- }
- case QuadratureType::GaussLegendre: {
- auto tensorial_gauss_legendre = QuadratureManager::instance().getLineGaussLegendreFormula(Order);
-
- m_quadrature_positions = tensorial_gauss_legendre.pointList();
- m_weights = tensorial_gauss_legendre.weightList();
- break;
- }
- case QuadratureType::QT__end: {
- throw NormalError("QuadratureType is not defined!");
- }
- default: {
- }
- }
- }
-
- // One does not use the '=default' constructor to avoid
- // (zero-initialization) performances issues
- PUGS_INLINE
- constexpr IntegrationTools() noexcept {}
-
- PUGS_INLINE
- constexpr IntegrationTools(const IntegrationTools&) noexcept = default;
-
- PUGS_INLINE
- constexpr IntegrationTools(IntegrationTools&& v) noexcept = default;
-
- PUGS_INLINE
- ~IntegrationTools() noexcept = default;
-};
-
-#endif // INTEGRATION_TOOLS_HPP
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index ba8ceb6c279f71dddff992aa4aebac11a629f79c..edd77f19444d13ec9024241d686eaa66085f6ada 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -73,7 +73,6 @@ add_executable (unit_tests
test_EmbeddedData.cpp
test_EmbeddedIDiscreteFunctionUtils.cpp
test_EscapedString.cpp
- test_FunctionEvaluation.cpp
test_Exceptions.cpp
test_ExecutionPolicy.cpp
test_FakeProcessor.cpp
diff --git a/tests/test_FunctionEvaluation.cpp b/tests/test_FunctionEvaluation.cpp
deleted file mode 100644
index 5702f2a2a1da295b853c086bdb342d0106d62fce..0000000000000000000000000000000000000000
--- a/tests/test_FunctionEvaluation.cpp
+++ /dev/null
@@ -1,140 +0,0 @@
-#include <catch2/catch_approx.hpp>
-#include <catch2/catch_test_macros.hpp>
-
-#include <algebra/IntegrationTools.hpp>
-#include <algebra/TinyVector.hpp>
-#include <language/utils/EvaluateAtPoints.hpp>
-#include <language/utils/FunctionSymbolId.hpp>
-#include <language/utils/SymbolTable.hpp>
-#include <utils/PugsAssert.hpp>
-
-// clazy:excludeall=non-pod-global-static
-
-TEST_CASE("IntegrationTools", "[integration]")
-{
- IntegrationTools<2> integrationtools(QuadratureType::GaussLobatto);
- double a = 1;
- double b = 2;
-
- SmallArray<TinyVector<1>> vertices(4);
- for (size_t i = 0; i < 4; i++) {
- vertices[0] = {a};
- vertices[1] = {b};
- vertices[2] = {b};
- vertices[3] = {a};
- }
- SECTION("Order")
- {
- size_t order = integrationtools.order();
- REQUIRE(order == 2);
- }
- SECTION("QuadraturePoints")
- {
- SmallArray<const TinyVector<1>> quadratures = integrationtools.quadraturePoints(a, b);
- REQUIRE(quadratures[0] == 1);
- REQUIRE(quadratures[1] == 1.5);
- REQUIRE(quadratures[2] == 2);
- }
- SECTION("Weights")
- {
- SmallArray<const double> weights = integrationtools.weights();
- REQUIRE(((weights[0] == 1. / 3.) and (weights[1] == 4. / 3.) and (weights[2] == 1. / 3.)));
- }
- SECTION("TestOrder3")
- {
- double result = integrationtools.testIntegrate(a, b);
- REQUIRE(result == Catch::Approx((0.25 * (std::pow(2, 4) - std::pow(1, 4)))).epsilon(1E-14));
- SmallArray<double> results = integrationtools.testIntegrateFunction(vertices);
- REQUIRE(results.size() == 3);
- REQUIRE(results[0] == Catch::Approx((0.25 * (std::pow(2, 4) - std::pow(1, 4)))).epsilon(1E-14));
- REQUIRE(results[1] == 0);
- REQUIRE(results[2] == Catch::Approx((-0.25 * (std::pow(2, 4) - std::pow(1, 4)))).epsilon(1E-14));
- }
- IntegrationTools<5> integrationtools5(QuadratureType::GaussLobatto);
- SECTION("TestOrder5")
- {
- double result = integrationtools5.testIntegrate(a, b);
- REQUIRE(result == Catch::Approx((0.25 * (std::pow(2, 4) - std::pow(1, 4)))).epsilon(1E-14));
- SmallArray<double> results = integrationtools.testIntegrateFunction(vertices);
- REQUIRE(results[0] == Catch::Approx((0.25 * (std::pow(2, 4) - std::pow(1, 4)))).epsilon(1E-14));
- REQUIRE(results[1] == 0);
- REQUIRE(results[2] == Catch::Approx((-0.25 * (std::pow(2, 4) - std::pow(1, 4)))).epsilon(1E-14));
- }
- IntegrationTools<7> integrationtools7(QuadratureType::GaussLobatto);
- SECTION("TestOrder7")
- {
- double result = integrationtools7.testIntegrate(a, b);
- REQUIRE(result == Catch::Approx((0.25 * (std::pow(2, 4) - std::pow(1, 4)))).epsilon(1E-14));
- SmallArray<double> results = integrationtools.testIntegrateFunction(vertices);
- REQUIRE(results[0] == Catch::Approx((0.25 * (std::pow(2, 4) - std::pow(1, 4)))).epsilon(1E-14));
- REQUIRE(results[1] == 0);
- REQUIRE(results[2] == Catch::Approx((-0.25 * (std::pow(2, 4) - std::pow(1, 4)))).epsilon(1E-14));
- }
- IntegrationTools<9> integrationtools9(QuadratureType::GaussLobatto);
- SECTION("TestOrder9")
- {
- double result = integrationtools9.testIntegrate(a, b);
- REQUIRE(result == Catch::Approx((0.25 * (std::pow(2, 4) - std::pow(1, 4)))).epsilon(1E-14));
- SmallArray<double> results = integrationtools.testIntegrateFunction(vertices);
- REQUIRE(results[0] == Catch::Approx((0.25 * (std::pow(2, 4) - std::pow(1, 4)))).epsilon(1E-14));
- REQUIRE(results[1] == 0);
- REQUIRE(results[2] == Catch::Approx((-0.25 * (std::pow(2, 4) - std::pow(1, 4)))).epsilon(1E-14));
- }
- IntegrationTools<11> integrationtools11(QuadratureType::GaussLobatto);
- SECTION("TestOrder11")
- {
- double result = integrationtools11.testIntegrate(a, b);
- REQUIRE(result == Catch::Approx((0.25 * (std::pow(2, 4) - std::pow(1, 4)))).epsilon(1E-14));
- SmallArray<double> results = integrationtools.testIntegrateFunction(vertices);
- REQUIRE(results[0] == Catch::Approx((0.25 * (std::pow(2, 4) - std::pow(1, 4)))).epsilon(1E-14));
- REQUIRE(results[1] == 0);
- REQUIRE(results[2] == Catch::Approx((-0.25 * (std::pow(2, 4) - std::pow(1, 4)))).epsilon(1E-14));
- }
-
- IntegrationTools<2> integrationtoolslegendre(QuadratureType::GaussLegendre);
- SECTION("Order-Legendre")
- {
- size_t order = integrationtools.order();
- REQUIRE(order == 2);
- }
- SECTION("TestOrder3-Legendre")
- {
- double result = integrationtoolslegendre.testIntegrate(a, b);
- REQUIRE(result == Catch::Approx((0.25 * (std::pow(2, 4) - std::pow(1, 4)))).epsilon(1E-14));
- SmallArray<double> results = integrationtoolslegendre.testIntegrateFunction(vertices);
- REQUIRE(results.size() == 3);
- REQUIRE(results[0] == Catch::Approx((0.25 * (std::pow(2, 4) - std::pow(1, 4)))).epsilon(1E-14));
- REQUIRE(results[1] == 0);
- REQUIRE(results[2] == Catch::Approx((-0.25 * (std::pow(2, 4) - std::pow(1, 4)))).epsilon(1E-14));
- }
- IntegrationTools<5> integrationtoolslegendre5(QuadratureType::GaussLegendre);
- SECTION("TestOrder5-Legendre")
- {
- double result = integrationtoolslegendre5.testIntegrate(a, b);
- REQUIRE(result == Catch::Approx((0.25 * (std::pow(2, 4) - std::pow(1, 4)))).epsilon(1E-14));
- SmallArray<double> results = integrationtoolslegendre5.testIntegrateFunction(vertices);
- REQUIRE(results[0] == Catch::Approx((0.25 * (std::pow(2, 4) - std::pow(1, 4)))).epsilon(1E-14));
- REQUIRE(results[1] == 0);
- REQUIRE(results[2] == Catch::Approx((-0.25 * (std::pow(2, 4) - std::pow(1, 4)))).epsilon(1E-14));
- }
- IntegrationTools<7> integrationtoolslegendre7(QuadratureType::GaussLegendre);
- SECTION("TestOrder7-Legendre")
- {
- double result = integrationtoolslegendre7.testIntegrate(a, b);
- REQUIRE(result == Catch::Approx((0.25 * (std::pow(2, 4) - std::pow(1, 4)))).epsilon(1E-14));
- SmallArray<double> results = integrationtoolslegendre7.testIntegrateFunction(vertices);
- REQUIRE(results[0] == Catch::Approx((0.25 * (std::pow(2, 4) - std::pow(1, 4)))).epsilon(1E-14));
- REQUIRE(results[1] == 0);
- REQUIRE(results[2] == Catch::Approx((-0.25 * (std::pow(2, 4) - std::pow(1, 4)))).epsilon(1E-14));
- }
- IntegrationTools<9> integrationtoolslegendre9(QuadratureType::GaussLegendre);
- SECTION("TestOrder9-Legendre")
- {
- double result = integrationtoolslegendre9.testIntegrate(a, b);
- REQUIRE(result == Catch::Approx((0.25 * (std::pow(2, 4) - std::pow(1, 4)))).epsilon(1E-14));
- SmallArray<double> results = integrationtoolslegendre9.testIntegrateFunction(vertices);
- REQUIRE(results[0] == Catch::Approx((0.25 * (std::pow(2, 4) - std::pow(1, 4)))).epsilon(1E-14));
- REQUIRE(results[1] == 0);
- REQUIRE(results[2] == Catch::Approx((-0.25 * (std::pow(2, 4) - std::pow(1, 4)))).epsilon(1E-14));
- }
-}