From 546807e762af6b2f3d83eb6940e24b34b490606e Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Del=20Pino?= <stephane.delpino44@gmail.com>
Date: Thu, 30 Sep 2021 18:51:17 +0200
Subject: [PATCH] Begin quadrature formula rework
Add bunch of Gauss-Legendre formula and begin testing them
---
src/algebra/IntegrationTools.hpp | 604 +++++++++++++++++++++++-------
tests/test_FunctionEvaluation.cpp | 87 ++++-
2 files changed, 552 insertions(+), 139 deletions(-)
diff --git a/src/algebra/IntegrationTools.hpp b/src/algebra/IntegrationTools.hpp
index bdf6de692..756fe09b9 100644
--- a/src/algebra/IntegrationTools.hpp
+++ b/src/algebra/IntegrationTools.hpp
@@ -1,15 +1,13 @@
#ifndef INTEGRATION_TOOLS_HPP
#define INTEGRATION_TOOLS_HPP
-#include <iostream>
-
+#include <language/utils/EvaluateAtPoints.hpp>
#include <utils/Exceptions.hpp>
#include <utils/PugsAssert.hpp>
#include <utils/PugsMacros.hpp>
+#include <utils/Types.hpp>
#include <cmath>
-#include <language/utils/EvaluateAtPoints.hpp>
-#include <utils/Types.hpp>
enum class QuadratureType : int8_t
{
@@ -21,15 +19,362 @@ enum class QuadratureType : int8_t
QT__end
};
-template <size_t Order>
-class IntegrationMethod
+template <size_t Dimension>
+class QuadratureBase
+{
+ protected:
+ size_t m_order;
+ SmallArray<const TinyVector<Dimension>> m_point_list;
+ SmallArray<const double> m_weight_list;
+
+ virtual void _buildPointAndWeightLists() = 0;
+
+ public:
+ size_t
+ numberOfPoints() const
+ {
+ Assert(m_point_list.size() == m_weight_list.size());
+ return m_point_list.size();
+ }
+
+ const SmallArray<const TinyVector<Dimension>>&
+ pointList() const
+ {
+ return m_point_list;
+ }
+
+ const SmallArray<const double>&
+ weightList() const
+ {
+ return m_weight_list;
+ }
+
+ QuadratureBase(QuadratureBase&&) = default;
+ QuadratureBase(const QuadratureBase&) = default;
+
+ protected:
+ explicit QuadratureBase(const size_t order) : m_order{order} {}
+
+ public:
+ QuadratureBase() = delete;
+ virtual ~QuadratureBase() = default;
+};
+
+template <size_t Dimension>
+class TensorialGaussLegendre final : public QuadratureBase<Dimension>
{
+ private:
+ void _buildPointAndWeightLists();
+
public:
- // PUGS_INLINE virtual T integrate(const FunctionSymbolId& function, const double& a, const double& b) = 0;
+ TensorialGaussLegendre(TensorialGaussLegendre&&) = default;
+ TensorialGaussLegendre(const TensorialGaussLegendre&) = default;
+
+ explicit TensorialGaussLegendre(const size_t order) : QuadratureBase<Dimension>(order)
+ {
+ this->_buildPointAndWeightLists();
+ }
+
+ TensorialGaussLegendre() = delete;
+ ~TensorialGaussLegendre() = default;
+};
+
+template <>
+void
+TensorialGaussLegendre<1>::_buildPointAndWeightLists()
+{
+ const size_t nb_points = m_order / 2 + 1;
+
+ SmallArray<TinyVector<1>> point_list(nb_points);
+ SmallArray<double> weight_list(nb_points);
+
+ switch (m_order) {
+ case 0:
+ case 1: {
+ point_list[0] = 0;
+
+ weight_list[0] = 2;
+ break;
+ }
+ case 2:
+ case 3: {
+ point_list[0] = -0.5773502691896257645091488;
+ point_list[1] = +0.5773502691896257645091488;
+
+ weight_list[0] = 1;
+ weight_list[1] = 1;
+ break;
+ }
+ case 4:
+ case 5: {
+ point_list[0] = -0.7745966692414833770358531;
+ point_list[1] = 0;
+ point_list[2] = +0.7745966692414833770358531;
+
+ weight_list[0] = 0.5555555555555555555555556;
+ weight_list[1] = 0.8888888888888888888888889;
+ weight_list[2] = 0.5555555555555555555555556;
+ break;
+ }
+ case 6:
+ case 7: {
+ point_list[0] = -0.8611363115940525752239465;
+ point_list[1] = -0.3399810435848562648026658;
+ point_list[2] = +0.3399810435848562648026658;
+ point_list[3] = +0.8611363115940525752239465;
+
+ weight_list[0] = 0.3478548451374538573730639;
+ weight_list[1] = 0.6521451548625461426269361;
+ weight_list[2] = 0.6521451548625461426269361;
+ weight_list[3] = 0.3478548451374538573730639;
+ break;
+ }
+ case 8:
+ case 9: {
+ point_list[0] = -0.9061798459386639927976269;
+ point_list[1] = -0.5384693101056830910363144;
+ point_list[2] = 0;
+ point_list[3] = +0.5384693101056830910363144;
+ point_list[4] = +0.9061798459386639927976269;
+
+ weight_list[0] = 0.2369268850561890875142640;
+ weight_list[1] = 0.4786286704993664680412915;
+ weight_list[2] = 0.5688888888888888888888889;
+ weight_list[3] = 0.4786286704993664680412915;
+ weight_list[4] = 0.2369268850561890875142640;
+ break;
+ }
+ case 10:
+ case 11: {
+ point_list[0] = -0.9324695142031520278123016;
+ point_list[1] = -0.6612093864662645136613996;
+ point_list[2] = -0.2386191860831969086305017;
+ point_list[3] = +0.2386191860831969086305017;
+ point_list[4] = +0.6612093864662645136613996;
+ point_list[5] = +0.9324695142031520278123016;
+
+ weight_list[0] = 0.1713244923791703450402961;
+ weight_list[1] = 0.3607615730481386075698335;
+ weight_list[2] = 0.4679139345726910473898703;
+ weight_list[3] = 0.4679139345726910473898703;
+ weight_list[4] = 0.3607615730481386075698335;
+ weight_list[5] = 0.1713244923791703450402961;
+ break;
+ }
+ case 12:
+ case 13: {
+ point_list[0] = -0.9491079123427585245261897;
+ point_list[1] = -0.7415311855993944398638648;
+ point_list[2] = -0.4058451513773971669066064;
+ point_list[3] = 0;
+ point_list[4] = +0.4058451513773971669066064;
+ point_list[5] = +0.7415311855993944398638648;
+ point_list[6] = +0.9491079123427585245261897;
+
+ weight_list[0] = 0.1294849661688696932706114;
+ weight_list[1] = 0.2797053914892766679014678;
+ weight_list[2] = 0.3818300505051189449503698;
+ weight_list[3] = 0.4179591836734693877551020;
+ weight_list[4] = 0.3818300505051189449503698;
+ weight_list[5] = 0.2797053914892766679014678;
+ weight_list[6] = 0.1294849661688696932706114;
+ break;
+ }
+ case 14:
+ case 15: {
+ point_list[0] = -0.9602898564975362316835609;
+ point_list[1] = -0.7966664774136267395915539;
+ point_list[2] = -0.5255324099163289858177390;
+ point_list[3] = -0.1834346424956498049394761;
+ point_list[4] = +0.1834346424956498049394761;
+ point_list[5] = +0.5255324099163289858177390;
+ point_list[6] = +0.7966664774136267395915539;
+ point_list[7] = +0.9602898564975362316835609;
+
+ weight_list[0] = 0.1012285362903762591525314;
+ weight_list[1] = 0.2223810344533744705443560;
+ weight_list[2] = 0.3137066458778872873379622;
+ weight_list[3] = 0.3626837833783619829651504;
+ weight_list[4] = 0.3626837833783619829651504;
+ weight_list[5] = 0.3137066458778872873379622;
+ weight_list[6] = 0.2223810344533744705443560;
+ weight_list[7] = 0.1012285362903762591525314;
+ break;
+ }
+ case 16:
+ case 17: {
+ point_list[0] = -0.9681602395076260898355762;
+ point_list[1] = -0.8360311073266357942994298;
+ point_list[2] = -0.6133714327005903973087020;
+ point_list[3] = -0.3242534234038089290385380;
+ point_list[4] = 0;
+ point_list[5] = +0.3242534234038089290385380;
+ point_list[6] = +0.6133714327005903973087020;
+ point_list[7] = +0.8360311073266357942994298;
+ point_list[8] = +0.9681602395076260898355762;
+
+ weight_list[0] = 0.0812743883615744119718922;
+ weight_list[1] = 0.1806481606948574040584720;
+ weight_list[2] = 0.2606106964029354623187429;
+ weight_list[3] = 0.3123470770400028400686304;
+ weight_list[4] = 0.3302393550012597631645251;
+ weight_list[5] = 0.3123470770400028400686304;
+ weight_list[6] = 0.2606106964029354623187429;
+ weight_list[7] = 0.1806481606948574040584720;
+ weight_list[8] = 0.0812743883615744119718922;
+ break;
+ }
+ case 18:
+ case 19: {
+ point_list[0] = -0.9739065285171717200779640;
+ point_list[1] = -0.8650633666889845107320967;
+ point_list[2] = -0.6794095682990244062343274;
+ point_list[3] = -0.4333953941292471907992659;
+ point_list[4] = -0.1488743389816312108848260;
+ point_list[5] = +0.1488743389816312108848260;
+ point_list[6] = +0.4333953941292471907992659;
+ point_list[7] = +0.6794095682990244062343274;
+ point_list[8] = +0.8650633666889845107320967;
+ point_list[9] = +0.9739065285171717200779640;
+
+ weight_list[0] = 0.0666713443086881375935688;
+ weight_list[1] = 0.1494513491505805931457763;
+ weight_list[2] = 0.2190863625159820439955349;
+ weight_list[3] = 0.2692667193099963550912269;
+ weight_list[4] = 0.2955242247147528701738930;
+ weight_list[5] = 0.2955242247147528701738930;
+ weight_list[6] = 0.2692667193099963550912269;
+ weight_list[7] = 0.2190863625159820439955349;
+ weight_list[8] = 0.1494513491505805931457763;
+ weight_list[9] = 0.0666713443086881375935688;
+ break;
+ }
+ case 20:
+ case 21: {
+ point_list[0] = -0.9782286581460569928039380;
+ point_list[1] = -0.8870625997680952990751578;
+ point_list[2] = -0.7301520055740493240934163;
+ point_list[3] = -0.5190961292068118159257257;
+ point_list[4] = -0.2695431559523449723315320;
+ point_list[5] = 0;
+ point_list[6] = +0.2695431559523449723315320;
+ point_list[7] = +0.5190961292068118159257257;
+ point_list[8] = +0.7301520055740493240934163;
+ point_list[9] = +0.8870625997680952990751578;
+ point_list[10] = +0.9782286581460569928039380;
+
+ weight_list[0] = 0.0556685671161736664827537;
+ weight_list[1] = 0.1255803694649046246346940;
+ weight_list[2] = 0.1862902109277342514260980;
+ weight_list[3] = 0.2331937645919904799185237;
+ weight_list[4] = 0.2628045445102466621806889;
+ weight_list[5] = 0.2729250867779006307144835;
+ weight_list[6] = 0.2628045445102466621806889;
+ weight_list[7] = 0.2331937645919904799185237;
+ weight_list[8] = 0.1862902109277342514260980;
+ weight_list[9] = 0.1255803694649046246346940;
+ weight_list[10] = 0.0556685671161736664827537;
+ break;
+ }
+ case 22:
+ case 23: {
+ point_list[0] = -0.9815606342467192506905491;
+ point_list[1] = -0.9041172563704748566784659;
+ point_list[2] = -0.7699026741943046870368938;
+ point_list[3] = -0.5873179542866174472967024;
+ point_list[4] = -0.3678314989981801937526915;
+ point_list[5] = -0.1252334085114689154724414;
+ point_list[6] = +0.1252334085114689154724414;
+ point_list[7] = +0.3678314989981801937526915;
+ point_list[8] = +0.5873179542866174472967024;
+ point_list[9] = +0.7699026741943046870368938;
+ point_list[10] = +0.9041172563704748566784659;
+ point_list[11] = +0.9815606342467192506905491;
+
+ weight_list[0] = 0.0471753363865118271946160;
+ weight_list[1] = 0.1069393259953184309602547;
+ weight_list[2] = 0.1600783285433462263346525;
+ weight_list[3] = 0.2031674267230659217490645;
+ weight_list[4] = 0.2334925365383548087608499;
+ weight_list[5] = 0.2491470458134027850005624;
+ weight_list[6] = 0.2491470458134027850005624;
+ weight_list[7] = 0.2334925365383548087608499;
+ weight_list[8] = 0.2031674267230659217490645;
+ weight_list[9] = 0.1600783285433462263346525;
+ weight_list[10] = 0.1069393259953184309602547;
+ weight_list[11] = 0.0471753363865118271946160;
+ break;
+ }
+ default: {
+ throw UnexpectedError("Gauss-Legendre quadratures handle orders up to 23");
+ break;
+ }
+ }
+
+ m_point_list = point_list;
+ m_weight_list = weight_list;
+}
- PUGS_INLINE virtual Array<TinyVector<1>> quadraturePoints(const double& a, const double& b) = 0;
- PUGS_INLINE virtual size_t numberOfPoints() = 0;
- PUGS_INLINE virtual Array<double> weights() = 0;
+template <>
+void
+TensorialGaussLegendre<2>::_buildPointAndWeightLists()
+{
+ const size_t nb_points_1d = this->m_order / 2 + 1;
+ const size_t nb_points = nb_points_1d * nb_points_1d;
+
+ SmallArray<TinyVector<2>> point_list(nb_points);
+ SmallArray<double> weight_list(nb_points);
+
+ TensorialGaussLegendre<1> gauss_legendre_1d(m_order);
+ const auto& point_list_1d = gauss_legendre_1d.pointList();
+ const auto& weight_list_1d = gauss_legendre_1d.weightList();
+
+ size_t l = 0;
+ for (size_t i = 0; i < nb_points_1d; ++i) {
+ for (size_t j = 0; j < nb_points_1d; ++j, ++l) {
+ point_list[l] = TinyVector<2>{point_list_1d[i][0], point_list_1d[j][0]};
+ weight_list[l] = weight_list_1d[i] * weight_list_1d[j];
+ }
+ }
+
+ this->m_point_list = point_list;
+ this->m_weight_list = weight_list;
+}
+
+template <>
+void
+TensorialGaussLegendre<3>::_buildPointAndWeightLists()
+{
+ const size_t nb_points_1d = this->m_order / 2 + 1;
+ const size_t nb_points = nb_points_1d * nb_points_1d * nb_points_1d;
+
+ SmallArray<TinyVector<3>> point_list(nb_points);
+ SmallArray<double> weight_list(nb_points);
+
+ TensorialGaussLegendre<1> gauss_legendre_1d(m_order);
+ const auto& point_list_1d = gauss_legendre_1d.pointList();
+ const auto& weight_list_1d = gauss_legendre_1d.weightList();
+
+ size_t l = 0;
+ for (size_t i = 0; i < nb_points_1d; ++i) {
+ for (size_t j = 0; j < nb_points_1d; ++j) {
+ for (size_t k = 0; k < nb_points_1d; ++k, ++l) {
+ point_list[l] = TinyVector<3>{point_list_1d[i][0], point_list_1d[j][0], point_list_1d[k][0]};
+ weight_list[l] = weight_list_1d[i] * weight_list_1d[j] * weight_list_1d[k];
+ }
+ }
+ }
+
+ this->m_point_list = point_list;
+ this->m_weight_list = weight_list;
+}
+
+class IntegrationMethod
+{
+ public:
+ PUGS_INLINE virtual SmallArray<TinyVector<1>> quadraturePoints(const double& a, const double& b) = 0;
+ PUGS_INLINE virtual size_t numberOfPoints() = 0;
+ PUGS_INLINE virtual SmallArray<double> weights() = 0;
// One does not use the '=default' constructor to avoid
// (zero-initialization) performances issues
PUGS_INLINE
@@ -45,16 +390,16 @@ class IntegrationMethod
};
template <size_t Order>
-class IntegrationMethodLobatto : public IntegrationMethod<Order>
+class IntegrationMethodLobatto : public IntegrationMethod
{
private:
static constexpr size_t number_points = (Order < 4 ? 3 : (Order < 6 ? 4 : (Order < 8 ? 5 : Order < 10 ? 6 : 7)));
- Array<double> m_weights;
- Array<TinyVector<1>> m_quadrature_positions;
- PUGS_INLINE Array<TinyVector<1>>
+ SmallArray<double> m_weights;
+ SmallArray<TinyVector<1>> m_quadrature_positions;
+ PUGS_INLINE SmallArray<TinyVector<1>>
translateQuadraturePoints(const double& a, const double& b)
{
- Array<TinyVector<1>> points{number_points};
+ SmallArray<TinyVector<1>> points{number_points};
double length = b - a;
for (size_t i = 0; i < m_quadrature_positions.size(); i++) {
TinyVector<1> qpos{m_quadrature_positions[i]};
@@ -64,7 +409,7 @@ class IntegrationMethodLobatto : public IntegrationMethod<Order>
}
PUGS_INLINE
- constexpr void
+ void
fillArrayLobatto()
{
switch (Order) {
@@ -180,7 +525,7 @@ class IntegrationMethodLobatto : public IntegrationMethod<Order>
}
PUGS_INLINE
- Array<TinyVector<1>>
+ SmallArray<TinyVector<1>>
quadraturePoints(const double& a, const double& b)
{
return translateQuadraturePoints(a, b);
@@ -192,7 +537,7 @@ class IntegrationMethodLobatto : public IntegrationMethod<Order>
return number_points;
}
PUGS_INLINE
- Array<double>
+ SmallArray<double>
weights()
{
return m_weights;
@@ -209,16 +554,16 @@ class IntegrationMethodLobatto : public IntegrationMethod<Order>
};
template <size_t Order>
-class IntegrationMethodLegendre : public IntegrationMethod<Order>
+class IntegrationMethodLegendre : public IntegrationMethod
{
private:
static constexpr size_t number_points = (Order < 4 ? 2 : (Order < 6 ? 3 : (Order < 8 ? 4 : 5)));
- Array<double> m_weights;
- Array<TinyVector<1>> m_quadrature_positions;
- PUGS_INLINE Array<TinyVector<1>>
+ SmallArray<double> m_weights;
+ SmallArray<TinyVector<1>> m_quadrature_positions;
+ PUGS_INLINE SmallArray<TinyVector<1>>
translateQuadraturePoints(const double& a, const double& b)
{
- Array<TinyVector<1>> points{number_points};
+ SmallArray<TinyVector<1>> points{number_points};
double length = b - a;
for (size_t i = 0; i < m_quadrature_positions.size(); i++) {
TinyVector<1> qpos{m_quadrature_positions[i]};
@@ -228,78 +573,81 @@ class IntegrationMethodLegendre : public IntegrationMethod<Order>
}
PUGS_INLINE
- constexpr void
+ void
fillArrayLegendre()
{
- switch (Order) {
- case 0:
- case 1:
- case 2:
- case 3: {
- double oneovsqrt3 = 1. / std::sqrt(3);
- m_quadrature_positions[0] = -oneovsqrt3;
- m_quadrature_positions[1] = oneovsqrt3;
- m_weights[0] = 1;
- m_weights[1] = 1;
- // m_quadrature_positions.fill({-1, 0, 1});
- // m_weights.fill({oneov3, 4 * oneov3, oneov3});
- break;
- }
- case 4:
- case 5: {
- double oneover9 = 1. / 9;
- double coef = std::sqrt(3. / 5.);
- m_quadrature_positions[0] = -coef;
- m_quadrature_positions[1] = 0;
- m_quadrature_positions[2] = coef;
- m_weights[0] = 5 * oneover9;
- m_weights[1] = 8 * oneover9;
- m_weights[2] = 5 * oneover9;
- break;
- }
- case 6:
- case 7: {
- double oneov36 = 1. / 36.;
- double coef1 = std::sqrt(3. / 7. - 2. / 7. * std::sqrt(6. / 5.));
- double coef2 = std::sqrt(3. / 7. + 2. / 7. * std::sqrt(6. / 5.));
- double coef3 = 18. + std::sqrt(30.);
- double coef4 = 18. - std::sqrt(30.);
-
- m_quadrature_positions[0] = -coef2;
- m_quadrature_positions[1] = -coef1;
- m_quadrature_positions[2] = coef1;
- m_quadrature_positions[3] = coef2;
- m_weights[0] = coef4 * oneov36;
- m_weights[1] = coef3 * oneov36;
- m_weights[2] = coef3 * oneov36;
- m_weights[3] = coef4 * oneov36;
- break;
- }
- case 8:
- case 9: {
- double oneov900 = 1. / 900.;
- double oneov225 = 1. / 225.;
- double coef1 = 1. / 3. * std::sqrt(5. - 2. * std::sqrt(10. / 7.));
- double coef2 = 1. / 3. * std::sqrt(5. + 2. * std::sqrt(10. / 7.));
- double coef3 = 322. + 13. * std::sqrt(70.);
- double coef4 = 322. - 13. * std::sqrt(70.);
- m_quadrature_positions[0] = -coef2;
- m_quadrature_positions[1] = -coef1;
- m_quadrature_positions[2] = 0;
- m_quadrature_positions[3] = coef1;
- m_quadrature_positions[4] = coef2;
- m_weights[0] = coef4 * oneov900;
- m_weights[1] = coef3 * oneov900;
- m_weights[2] = 128. * oneov225;
- m_weights[3] = coef3 * oneov900;
- m_weights[4] = coef4 * oneov900;
- break;
- }
- default: {
- throw UnexpectedError("Gauss-Legendre quadratures handle orders up to 9.");
- break;
- }
- }
+ TensorialGaussLegendre<1> tensorial_gauss_legendre(Order);
+ copy_to(tensorial_gauss_legendre.pointList(), m_quadrature_positions);
+ copy_to(tensorial_gauss_legendre.weightList(), m_weights);
+ // switch (Order) {
+ // case 0:
+ // case 1:
+ // case 2:
+ // case 3: {
+ // double oneovsqrt3 = 1. / std::sqrt(3);
+ // m_quadrature_positions[0] = -oneovsqrt3;
+ // m_quadrature_positions[1] = oneovsqrt3;
+ // m_weights[0] = 1;
+ // m_weights[1] = 1;
+ // // m_quadrature_positions.fill({-1, 0, 1});
+ // // m_weights.fill({oneov3, 4 * oneov3, oneov3});
+ // break;
+ // }
+ // case 4:
+ // case 5: {
+ // double oneover9 = 1. / 9;
+ // double coef = std::sqrt(3. / 5.);
+ // m_quadrature_positions[0] = -coef;
+ // m_quadrature_positions[1] = 0;
+ // m_quadrature_positions[2] = coef;
+ // m_weights[0] = 5 * oneover9;
+ // m_weights[1] = 8 * oneover9;
+ // m_weights[2] = 5 * oneover9;
+ // break;
+ // }
+ // case 6:
+ // case 7: {
+ // double oneov36 = 1. / 36.;
+ // double coef1 = std::sqrt(3. / 7. - 2. / 7. * std::sqrt(6. / 5.));
+ // double coef2 = std::sqrt(3. / 7. + 2. / 7. * std::sqrt(6. / 5.));
+ // double coef3 = 18. + std::sqrt(30.);
+ // double coef4 = 18. - std::sqrt(30.);
+
+ // m_quadrature_positions[0] = -coef2;
+ // m_quadrature_positions[1] = -coef1;
+ // m_quadrature_positions[2] = coef1;
+ // m_quadrature_positions[3] = coef2;
+ // m_weights[0] = coef4 * oneov36;
+ // m_weights[1] = coef3 * oneov36;
+ // m_weights[2] = coef3 * oneov36;
+ // m_weights[3] = coef4 * oneov36;
+ // break;
+ // }
+ // case 8:
+ // case 9: {
+ // double oneov900 = 1. / 900.;
+ // double oneov225 = 1. / 225.;
+ // double coef1 = 1. / 3. * std::sqrt(5. - 2. * std::sqrt(10. / 7.));
+ // double coef2 = 1. / 3. * std::sqrt(5. + 2. * std::sqrt(10. / 7.));
+ // double coef3 = 322. + 13. * std::sqrt(70.);
+ // double coef4 = 322. - 13. * std::sqrt(70.);
+ // m_quadrature_positions[0] = -coef2;
+ // m_quadrature_positions[1] = -coef1;
+ // m_quadrature_positions[2] = 0;
+ // m_quadrature_positions[3] = coef1;
+ // m_quadrature_positions[4] = coef2;
+ // m_weights[0] = coef4 * oneov900;
+ // m_weights[1] = coef3 * oneov900;
+ // m_weights[2] = 128. * oneov225;
+ // m_weights[3] = coef3 * oneov900;
+ // m_weights[4] = coef4 * oneov900;
+ // break;
+ // }
+ // default: {
+ // throw UnexpectedError("Gauss-Legendre quadratures handle orders up to 9.");
+ // break;
+ // }
+ // }
}
public:
@@ -314,7 +662,7 @@ class IntegrationMethodLegendre : public IntegrationMethod<Order>
}
PUGS_INLINE
- Array<TinyVector<1>>
+ SmallArray<TinyVector<1>>
quadraturePoints(const double& a, const double& b)
{
return translateQuadraturePoints(a, b);
@@ -326,7 +674,7 @@ class IntegrationMethodLegendre : public IntegrationMethod<Order>
return number_points;
}
PUGS_INLINE
- Array<double>
+ SmallArray<double>
weights()
{
return m_weights;
@@ -350,17 +698,18 @@ class IntegrationTools
private:
// T m_values[N];
- std::shared_ptr<IntegrationMethod<Order>> m_integration_method;
+ std::shared_ptr<IntegrationMethod> m_integration_method;
size_t N;
public:
PUGS_INLINE T
integrate(const FunctionSymbolId& function, const double& a, const double& b) const
{
- T result = 0;
- Array<const TinyVector<1>> positions = quadraturePoints(a, b);
- Array<double> weights = m_integration_method->weights();
- Array<T> values = EvaluateAtPoints<double(const TinyVector<1>)>::evaluate(function, positions);
+ T result = 0;
+ SmallArray<const TinyVector<1>> positions = quadraturePoints(a, b);
+ SmallArray<double> weights = m_integration_method->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];
@@ -368,14 +717,15 @@ class IntegrationTools
return (b - a) / 2 * result;
}
- PUGS_INLINE Array<T>
- integrateFunction(const FunctionSymbolId& function, const Array<TinyVector<1>>& vertices) const
+ PUGS_INLINE SmallArray<T>
+ integrateFunction(const FunctionSymbolId& function, const SmallArray<TinyVector<1>>& vertices) const
{
- Array<const TinyVector<1>> positions = quadraturePoints(vertices);
- Array<T> result{vertices.size() - 1};
- Array<double> interval_size{vertices.size() - 1};
- Array<double> weights = m_integration_method->weights();
- Array<T> values = EvaluateAtPoints<double(const TinyVector<1>)>::evaluate(function, positions);
+ SmallArray<const TinyVector<1>> positions = quadraturePoints(vertices);
+ SmallArray<T> result{vertices.size() - 1};
+ SmallArray<double> interval_size{vertices.size() - 1};
+ SmallArray<double> weights = m_integration_method->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];
@@ -391,27 +741,31 @@ class IntegrationTools
PUGS_INLINE T
testIntegrate(const double& a, const double& b) const
{
- T result = 0;
- Array<TinyVector<1>> positions = quadraturePoints(a, b);
- Array<double> weights = m_integration_method->weights();
- Array<T> values(weights.size());
+ T result = 0;
+ SmallArray<TinyVector<1>> positions = quadraturePoints(a, b);
+
+ SmallArray<double> weights = m_integration_method->weights();
+
+ SmallArray<T> values(weights.size());
for (size_t j = 0; j < weights.size(); ++j) {
- values[j] = std::pow(positions[j][0], 3);
+ 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 (b - a) / 2 * result;
+
+ return 0.5 * (b - a) * result;
}
- PUGS_INLINE Array<T>
- testIntegrateFunction(const Array<TinyVector<1>>& vertices) const
+ PUGS_INLINE SmallArray<T>
+ testIntegrateFunction(const SmallArray<TinyVector<1>>& vertices) const
{
- Array<TinyVector<1>> positions = quadraturePoints(vertices);
- Array<double> interval_size{vertices.size() - 1};
- Array<double> weights = m_integration_method->weights();
- Array<T> values(positions.size());
+ SmallArray<TinyVector<1>> positions = quadraturePoints(vertices);
+ SmallArray<double> interval_size{vertices.size() - 1};
+ SmallArray<double> weights = m_integration_method->weights();
+ SmallArray<T> values(positions.size());
for (size_t j = 0; j < positions.size(); ++j) {
values[j] = std::pow(positions[j][0], 3);
}
@@ -419,7 +773,7 @@ class IntegrationTools
for (size_t i = 0; i < interval_size.size(); ++i) {
interval_size[i] = vertices[i + 1][0] - vertices[i][0];
}
- Array<T> result{vertices.size() - 1};
+ 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>) {
@@ -435,19 +789,19 @@ class IntegrationTools
return result;
}
- PUGS_INLINE Array<TinyVector<1>>
+ PUGS_INLINE SmallArray<TinyVector<1>>
quadraturePoints(const double& a, const double& b) const
{
return m_integration_method->quadraturePoints(a, b);
}
PUGS_INLINE
- Array<TinyVector<1>>
- quadraturePoints(const Array<TinyVector<1>>& positions) const
+ SmallArray<TinyVector<1>>
+ quadraturePoints(const SmallArray<TinyVector<1>>& positions) const
{
size_t number_of_intervals = positions.size() - 1;
size_t quadrature_size = m_integration_method->numberOfPoints();
- Array<TinyVector<1>> quadratures{quadrature_size * number_of_intervals};
+ 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];
@@ -459,7 +813,7 @@ class IntegrationTools
return quadratures;
}
- PUGS_INLINE Array<double>
+ PUGS_INLINE SmallArray<double>
weights() const
{
return m_integration_method->weights();
diff --git a/tests/test_FunctionEvaluation.cpp b/tests/test_FunctionEvaluation.cpp
index a7a79d156..88c9e552a 100644
--- a/tests/test_FunctionEvaluation.cpp
+++ b/tests/test_FunctionEvaluation.cpp
@@ -24,19 +24,76 @@ TEST_CASE("EvaluateAtPoints", "[integration]")
SECTION("Instantiation")
{
// positions[0] = 1;
- Array<double> result(1);
+ SmallArray<double> result(1);
// result = EvaluateAtPoints<double(TinyVector<1>)>::evaluate(f, positions);
// REQUIRE(result[0] == 1);
}
}
+TEST_CASE("QuadratureFormula", "[analysis]")
+{
+ SECTION("TensorialGaussLegendre")
+ {
+ auto integrate = [](auto f, auto quadrature_formula, const double a, const double b) {
+ auto point_list = quadrature_formula.pointList();
+ auto weight_list = quadrature_formula.weightList();
+
+ auto x = [&a, &b](auto x_hat) { return 0.5 * (b - a) * (x_hat + 1); };
+
+ auto value = weight_list[0] * f(x(point_list[0]));
+ for (size_t i = 1; i < weight_list.size(); ++i) {
+ value += weight_list[i] * f(x(point_list[i]));
+ }
+ return 0.5 * (b - a) * value;
+ };
+
+ auto p0 = [](const TinyVector<1>&) { return 3; };
+ auto p1 = [](const TinyVector<1>& X) {
+ const double x = X[0];
+ return 2 * x + 1;
+ };
+ auto p2 = [](const TinyVector<1>& X) {
+ const double x = X[0];
+ return 3 * x * x + 2 * x + 1;
+ };
+
+ SECTION("exact integration")
+ {
+ SECTION("order 0 and 1")
+ {
+ TensorialGaussLegendre<1> l0(0);
+ TensorialGaussLegendre<1> l1(1);
+
+ REQUIRE(l0.numberOfPoints() == 1);
+ REQUIRE(l1.numberOfPoints() == 1);
+
+ REQUIRE(integrate(p0, l1, 0.5, 1) == Catch::Approx(1.5));
+ REQUIRE(integrate(p1, l1, 0, 1) == Catch::Approx(2));
+ }
+
+ SECTION("order 2 and 3")
+ {
+ TensorialGaussLegendre<1> l2(2);
+ TensorialGaussLegendre<1> l3(3);
+
+ REQUIRE(l2.numberOfPoints() == 2);
+ REQUIRE(l3.numberOfPoints() == 2);
+
+ REQUIRE(integrate(p0, l3, 0.5, 1) == Catch::Approx(1.5));
+ REQUIRE(integrate(p1, l3, 0, 1) == Catch::Approx(2));
+ REQUIRE(integrate(p2, l3, 0, 1) == Catch::Approx(3));
+ }
+ }
+ }
+}
+
TEST_CASE("IntegrationTools", "[integration]")
{
IntegrationTools<2> integrationtools(QuadratureType::gausslobatto);
double a = 1;
double b = 2;
- Array<TinyVector<1>> vertices(4);
+ SmallArray<TinyVector<1>> vertices(4);
for (size_t i = 0; i < 4; i++) {
vertices[0] = {a};
vertices[1] = {b};
@@ -50,19 +107,21 @@ TEST_CASE("IntegrationTools", "[integration]")
}
SECTION("QuadraturePoints")
{
- Array<TinyVector<1>> quadratures = integrationtools.quadraturePoints(a, b);
- REQUIRE(((quadratures[0] == 1) and (quadratures[1] == 1.5) and (quadratures[2] == 2)));
+ SmallArray<TinyVector<1>> quadratures = integrationtools.quadraturePoints(a, b);
+ REQUIRE(quadratures[0] == 1);
+ REQUIRE(quadratures[1] == 1.5);
+ REQUIRE(quadratures[2] == 2);
}
SECTION("Weights")
{
- Array<double> weights = integrationtools.weights();
+ SmallArray<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));
- Array<double> results = integrationtools.testIntegrateFunction(vertices);
+ 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);
@@ -73,7 +132,7 @@ TEST_CASE("IntegrationTools", "[integration]")
{
double result = integrationtools5.testIntegrate(a, b);
REQUIRE(result == Catch::Approx((0.25 * (std::pow(2, 4) - std::pow(1, 4)))).epsilon(1E-14));
- Array<double> results = integrationtools.testIntegrateFunction(vertices);
+ 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));
@@ -83,7 +142,7 @@ TEST_CASE("IntegrationTools", "[integration]")
{
double result = integrationtools7.testIntegrate(a, b);
REQUIRE(result == Catch::Approx((0.25 * (std::pow(2, 4) - std::pow(1, 4)))).epsilon(1E-14));
- Array<double> results = integrationtools.testIntegrateFunction(vertices);
+ 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));
@@ -93,7 +152,7 @@ TEST_CASE("IntegrationTools", "[integration]")
{
double result = integrationtools9.testIntegrate(a, b);
REQUIRE(result == Catch::Approx((0.25 * (std::pow(2, 4) - std::pow(1, 4)))).epsilon(1E-14));
- Array<double> results = integrationtools.testIntegrateFunction(vertices);
+ 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));
@@ -103,7 +162,7 @@ TEST_CASE("IntegrationTools", "[integration]")
{
double result = integrationtools11.testIntegrate(a, b);
REQUIRE(result == Catch::Approx((0.25 * (std::pow(2, 4) - std::pow(1, 4)))).epsilon(1E-14));
- Array<double> results = integrationtools.testIntegrateFunction(vertices);
+ 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));
@@ -119,7 +178,7 @@ TEST_CASE("IntegrationTools", "[integration]")
{
double result = integrationtoolslegendre.testIntegrate(a, b);
REQUIRE(result == Catch::Approx((0.25 * (std::pow(2, 4) - std::pow(1, 4)))).epsilon(1E-14));
- Array<double> results = integrationtoolslegendre.testIntegrateFunction(vertices);
+ 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);
@@ -130,7 +189,7 @@ TEST_CASE("IntegrationTools", "[integration]")
{
double result = integrationtoolslegendre5.testIntegrate(a, b);
REQUIRE(result == Catch::Approx((0.25 * (std::pow(2, 4) - std::pow(1, 4)))).epsilon(1E-14));
- Array<double> results = integrationtoolslegendre5.testIntegrateFunction(vertices);
+ 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));
@@ -140,7 +199,7 @@ TEST_CASE("IntegrationTools", "[integration]")
{
double result = integrationtoolslegendre7.testIntegrate(a, b);
REQUIRE(result == Catch::Approx((0.25 * (std::pow(2, 4) - std::pow(1, 4)))).epsilon(1E-14));
- Array<double> results = integrationtoolslegendre7.testIntegrateFunction(vertices);
+ 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));
@@ -150,7 +209,7 @@ TEST_CASE("IntegrationTools", "[integration]")
{
double result = integrationtoolslegendre9.testIntegrate(a, b);
REQUIRE(result == Catch::Approx((0.25 * (std::pow(2, 4) - std::pow(1, 4)))).epsilon(1E-14));
- Array<double> results = integrationtoolslegendre9.testIntegrateFunction(vertices);
+ 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));
--
GitLab