Skip to content
Snippets Groups Projects
Commit 9ed28691 authored by Stéphane Del Pino's avatar Stéphane Del Pino
Browse files

Remove legacy IntegrationTools and associated tests

parent dbcabee3
No related branches found
No related tags found
1 merge request!124Add files for high order integration with quadratures
Show whitespace changes
Inline Side-by-side
src/algebra/IntegrationTools.hpp deleted 100644 → 0
+ 0
− 197
View file @ dbcabee3
#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
tests/CMakeLists.txt
+ 0
− 1
View file @ 9ed28691
...@@ -73,7 +73,6 @@ add_executable (unit_tests ...@@ -73,7 +73,6 @@ add_executable (unit_tests
test_EmbeddedData.cpp test_EmbeddedData.cpp
test_EmbeddedIDiscreteFunctionUtils.cpp test_EmbeddedIDiscreteFunctionUtils.cpp
test_EscapedString.cpp test_EscapedString.cpp
test_FunctionEvaluation.cpp
test_Exceptions.cpp test_Exceptions.cpp
test_ExecutionPolicy.cpp test_ExecutionPolicy.cpp
test_FakeProcessor.cpp test_FakeProcessor.cpp
......
tests/test_FunctionEvaluation.cpp deleted 100644 → 0
+ 0
− 140
View file @ dbcabee3
#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));
}
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment