From a5b213029ba3bf19656c9c44adf334f9c183e87d Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Del=20Pino?= <stephane.delpino44@gmail.com>
Date: Mon, 4 Oct 2021 15:27:57 +0200
Subject: [PATCH] Displace Gauss Legendre quadrature formula in its own file
All orders are tested in 1d.
It is quite tricky to test correctly these formulae since starting
from order 10 formula give a very (very actually) precise evaluation
of common functions.
---
CMakeLists.txt | 3 +-
src/CMakeLists.txt | 3 +
src/algebra/IntegrationTools.hpp | 354 +---------------
src/analysis/CMakeLists.txt | 5 +
src/analysis/GaussLegendreQuadrature.cpp | 291 +++++++++++++
src/analysis/GaussLegendreQuadrature.hpp | 25 ++
src/analysis/QuadratureFormula.hpp | 52 +++
tests/CMakeLists.txt | 3 +
tests/test_FunctionEvaluation.cpp | 105 -----
tests/test_GaussLegendreQuadrature.cpp | 503 +++++++++++++++++++++++
10 files changed, 887 insertions(+), 457 deletions(-)
create mode 100644 src/analysis/CMakeLists.txt
create mode 100644 src/analysis/GaussLegendreQuadrature.cpp
create mode 100644 src/analysis/GaussLegendreQuadrature.hpp
create mode 100644 src/analysis/QuadratureFormula.hpp
create mode 100644 tests/test_GaussLegendreQuadrature.cpp
diff --git a/CMakeLists.txt b/CMakeLists.txt
index e9255bff6..4b302fe37 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -573,6 +573,7 @@ target_link_libraries(
pugs
PugsMesh
PugsAlgebra
+ PugsAnalysis
PugsUtils
PugsLanguage
PugsLanguageAST
@@ -606,7 +607,7 @@ if(${CMAKE_VERSION} VERSION_LESS "3.13.0")
LIBRARY DESTINATION lib
ARCHIVE DESTINATION lib)
else()
- install(TARGETS pugs PugsAlgebra PugsUtils PugsLanguage PugsLanguageAST PugsLanguageModules PugsLanguageAlgorithms PugsLanguageUtils PugsMesh PugsScheme PugsOutput
+ install(TARGETS pugs PugsAlgebra PugsAnalysis PugsUtils PugsLanguage PugsLanguageAST PugsLanguageModules PugsLanguageAlgorithms PugsLanguageUtils PugsMesh PugsScheme PugsOutput
RUNTIME DESTINATION bin
LIBRARY DESTINATION lib
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 16904a7be..54cae625b 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -12,6 +12,9 @@ add_subdirectory(language)
# Pugs algebra
add_subdirectory(algebra)
+# Pugs analysis
+add_subdirectory(analysis)
+
# Pugs mesh
add_subdirectory(mesh)
diff --git a/src/algebra/IntegrationTools.hpp b/src/algebra/IntegrationTools.hpp
index 756fe09b9..08642b794 100644
--- a/src/algebra/IntegrationTools.hpp
+++ b/src/algebra/IntegrationTools.hpp
@@ -7,6 +7,8 @@
#include <utils/PugsMacros.hpp>
#include <utils/Types.hpp>
+#include <analysis/GaussLegendreQuadrature.hpp>
+
#include <cmath>
enum class QuadratureType : int8_t
@@ -19,356 +21,6 @@ enum class QuadratureType : int8_t
QT__end
};
-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:
- 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;
-}
-
-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:
@@ -576,7 +228,7 @@ class IntegrationMethodLegendre : public IntegrationMethod
void
fillArrayLegendre()
{
- TensorialGaussLegendre<1> tensorial_gauss_legendre(Order);
+ GaussLegendreQuadrature<1> tensorial_gauss_legendre(Order);
copy_to(tensorial_gauss_legendre.pointList(), m_quadrature_positions);
copy_to(tensorial_gauss_legendre.weightList(), m_weights);
// switch (Order) {
diff --git a/src/analysis/CMakeLists.txt b/src/analysis/CMakeLists.txt
new file mode 100644
index 000000000..b2e2f4c2f
--- /dev/null
+++ b/src/analysis/CMakeLists.txt
@@ -0,0 +1,5 @@
+# ------------------- Source files --------------------
+
+add_library(
+ PugsAnalysis
+ GaussLegendreQuadrature.cpp)
diff --git a/src/analysis/GaussLegendreQuadrature.cpp b/src/analysis/GaussLegendreQuadrature.cpp
new file mode 100644
index 000000000..e9cc2b963
--- /dev/null
+++ b/src/analysis/GaussLegendreQuadrature.cpp
@@ -0,0 +1,291 @@
+#include <analysis/GaussLegendreQuadrature.hpp>
+#include <utils/Exceptions.hpp>
+
+template <>
+void
+GaussLegendreQuadrature<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 NormalError("Gauss-Legendre quadrature formulae handle orders up to 23");
+ }
+ }
+
+ m_point_list = point_list;
+ m_weight_list = weight_list;
+}
+
+template <>
+void
+GaussLegendreQuadrature<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);
+
+ GaussLegendreQuadrature<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
+GaussLegendreQuadrature<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);
+
+ GaussLegendreQuadrature<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;
+}
diff --git a/src/analysis/GaussLegendreQuadrature.hpp b/src/analysis/GaussLegendreQuadrature.hpp
new file mode 100644
index 000000000..3e4a8b0cc
--- /dev/null
+++ b/src/analysis/GaussLegendreQuadrature.hpp
@@ -0,0 +1,25 @@
+#ifndef GAUSS_LEGENDRE_QUADRATURE_HPP
+#define GAUSS_LEGENDRE_QUADRATURE_HPP
+
+#include <analysis/QuadratureFormula.hpp>
+
+template <size_t Dimension>
+class GaussLegendreQuadrature final : public QuadratureForumla<Dimension>
+{
+ private:
+ void _buildPointAndWeightLists();
+
+ public:
+ GaussLegendreQuadrature(GaussLegendreQuadrature&&) = default;
+ GaussLegendreQuadrature(const GaussLegendreQuadrature&) = default;
+
+ explicit GaussLegendreQuadrature(const size_t order) : QuadratureForumla<Dimension>(order)
+ {
+ this->_buildPointAndWeightLists();
+ }
+
+ GaussLegendreQuadrature() = delete;
+ ~GaussLegendreQuadrature() = default;
+};
+
+#endif // GAUSS_LEGENDRE_QUADRATURE_HPP
diff --git a/src/analysis/QuadratureFormula.hpp b/src/analysis/QuadratureFormula.hpp
new file mode 100644
index 000000000..c92d31872
--- /dev/null
+++ b/src/analysis/QuadratureFormula.hpp
@@ -0,0 +1,52 @@
+#ifndef QUADRATURE_FORMULA_HPP
+#define QUADRATURE_FORMULA_HPP
+
+#include <algebra/TinyVector.hpp>
+#include <utils/PugsAssert.hpp>
+#include <utils/PugsMacros.hpp>
+#include <utils/SmallArray.hpp>
+
+template <size_t Dimension>
+class QuadratureForumla
+{
+ protected:
+ size_t m_order;
+ SmallArray<const TinyVector<Dimension>> m_point_list;
+ SmallArray<const double> m_weight_list;
+
+ virtual void _buildPointAndWeightLists() = 0;
+
+ public:
+ PUGS_INLINE size_t
+ numberOfPoints() const
+ {
+ Assert(m_point_list.size() == m_weight_list.size());
+ return m_point_list.size();
+ }
+
+ PUGS_INLINE
+ const SmallArray<const TinyVector<Dimension>>&
+ pointList() const
+ {
+ return m_point_list;
+ }
+
+ PUGS_INLINE
+ const SmallArray<const double>&
+ weightList() const
+ {
+ return m_weight_list;
+ }
+
+ QuadratureForumla(QuadratureForumla&&) = default;
+ QuadratureForumla(const QuadratureForumla&) = default;
+
+ protected:
+ explicit QuadratureForumla(const size_t order) : m_order{order} {}
+
+ public:
+ QuadratureForumla() = delete;
+ virtual ~QuadratureForumla() = default;
+};
+
+#endif // QUADRATURE_FORMULA_HPP
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 5cc4fee50..b0b2b7328 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -80,6 +80,7 @@ add_executable (unit_tests
test_FunctionProcessor.cpp
test_FunctionSymbolId.cpp
test_FunctionTable.cpp
+ test_GaussLegendreQuadrature.cpp
test_IfProcessor.cpp
test_IncDecExpressionProcessor.cpp
test_INodeProcessor.cpp
@@ -145,6 +146,7 @@ target_link_libraries (unit_tests
PugsLanguage
PugsMesh
PugsAlgebra
+ PugsAnalysis
PugsScheme
PugsOutput
PugsUtils
@@ -160,6 +162,7 @@ target_link_libraries (unit_tests
target_link_libraries (mpi_unit_tests
test_Pugs_MeshDataBase
PugsAlgebra
+ PugsAnalysis
PugsUtils
PugsLanguage
PugsLanguageAST
diff --git a/tests/test_FunctionEvaluation.cpp b/tests/test_FunctionEvaluation.cpp
index 7893ba07e..2e290556e 100644
--- a/tests/test_FunctionEvaluation.cpp
+++ b/tests/test_FunctionEvaluation.cpp
@@ -7,114 +7,9 @@
#include <language/utils/FunctionSymbolId.hpp>
#include <language/utils/SymbolTable.hpp>
#include <utils/PugsAssert.hpp>
-// Instantiate to ensure full coverage is performed
-// template class EvaluateAtPoints<double(TinyVector<1>)>;
-template class TinyVector<1>;
-template class EvaluateAtPoints<double(TinyVector<1>)>;
-template class IntegrationTools<2>;
// clazy:excludeall=non-pod-global-static
-TEST_CASE("EvaluateAtPoints", "[integration]")
-{
- std::shared_ptr<SymbolTable> s = std::make_shared<SymbolTable>();
- const FunctionSymbolId f{2, s};
- const Array<TinyVector<1>> positions(1);
- positions[0] = 1;
- SECTION("Instantiation")
- {
- // positions[0] = 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();
-
- double alpha = 0.5 * (b - a);
- double beta = 0.5 * (a + b);
-
- auto x = [&alpha, &beta](auto x_hat) { return alpha * x_hat + beta; };
-
- 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 alpha * value;
- };
-
- auto get_order = [&integrate](auto f, auto quadrature_formula, const double a, const double b,
- const double exact_value) {
- double int_ab = integrate(f, quadrature_formula, a, b);
- double int_first_half = integrate(f, quadrature_formula, a, 0.5 * (a + b));
- double int_second_half = integrate(f, quadrature_formula, 0.5 * (a + b), b);
-
- return -std::log((int_first_half + int_second_half - exact_value) / (int_ab - exact_value)) / std::log(2);
- };
-
- auto p0 = [](const TinyVector<1>&) { return 1; };
- 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;
- };
- auto p3 = [](const TinyVector<1>& X) {
- const double x = X[0];
- return 4 * std::pow(x, 3) + 3 * x * x + 2 * x + 1;
- };
- auto p4 = [](const TinyVector<1>& X) {
- const double x = X[0];
- return 5 * std::pow(x, 4) + 4 * std::pow(x, 3) + 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(0.5));
- REQUIRE(integrate(p1, l1, 0, 1) == Catch::Approx(2));
- REQUIRE(integrate(p2, l1, 0, 1) != Catch::Approx(3));
-
- REQUIRE(get_order(p2, l1, 0, 1, 3) == 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(0.5));
- REQUIRE(integrate(p1, l3, 0, 1) == Catch::Approx(2));
- REQUIRE(integrate(p2, l3, 0, 1) == Catch::Approx(3));
- REQUIRE(integrate(p3, l3, 0, 1) == Catch::Approx(4));
- REQUIRE(integrate(p4, l3, 0, 1) != Catch::Approx(5));
-
- REQUIRE(get_order(p4, l3, 0, 1, 5) == Catch::Approx(4));
- }
- }
- }
-}
-
TEST_CASE("IntegrationTools", "[integration]")
{
IntegrationTools<2> integrationtools(QuadratureType::gausslobatto);
diff --git a/tests/test_GaussLegendreQuadrature.cpp b/tests/test_GaussLegendreQuadrature.cpp
new file mode 100644
index 000000000..f20013560
--- /dev/null
+++ b/tests/test_GaussLegendreQuadrature.cpp
@@ -0,0 +1,503 @@
+#include <catch2/catch_approx.hpp>
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_all.hpp>
+
+#include <analysis/GaussLegendreQuadrature.hpp>
+#include <utils/Exceptions.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("GaussLegendreQuadrature", "[analysis]")
+{
+ auto same_formulas = [](auto quadrature_formula0, auto quadrature_formula1) {
+ auto point_list0 = quadrature_formula0.pointList();
+ auto weight_list0 = quadrature_formula0.weightList();
+
+ auto point_list1 = quadrature_formula1.pointList();
+ auto weight_list1 = quadrature_formula1.weightList();
+
+ bool is_same = true;
+
+ for (size_t i = 0; i < point_list0.size(); ++i) {
+ if (point_list0[i] != point_list1[i]) {
+ return false;
+ }
+ }
+ for (size_t i = 0; i < weight_list0.size(); ++i) {
+ if (weight_list0[i] != weight_list1[i]) {
+ return false;
+ }
+ }
+
+ return is_same;
+ };
+
+ SECTION("1D")
+ {
+ 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();
+
+ double alpha = 0.5 * (b - a);
+ double beta = 0.5 * (a + b);
+
+ auto x = [&alpha, &beta](auto x_hat) { return alpha * x_hat + beta; };
+
+ 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 alpha * value;
+ };
+
+ auto get_order = [&integrate](auto f, auto quadrature_formula, const double a, const double b,
+ const double exact_value) {
+ double int_ab = integrate(f, quadrature_formula, a, b);
+ double int_first_half = integrate(f, quadrature_formula, a, 0.5 * (a + b));
+ double int_second_half = integrate(f, quadrature_formula, 0.5 * (a + b), b);
+
+ return -std::log((int_first_half + int_second_half - exact_value) / (int_ab - exact_value)) / std::log(2);
+ };
+
+ auto p0 = [](const TinyVector<1>&) { return 1; };
+ auto p1 = [&p0](const TinyVector<1>& X) {
+ const double x = X[0];
+ return 2 * x + p0(X);
+ };
+ auto p2 = [&p1](const TinyVector<1>& X) {
+ const double x = X[0];
+ return 3 * x * x + p1(X);
+ };
+ auto p3 = [&p2](const TinyVector<1>& X) {
+ const double x = X[0];
+ return 4 * std::pow(x, 3) + p2(X);
+ };
+ auto p4 = [&p3](const TinyVector<1>& X) {
+ const double x = X[0];
+ return 5 * std::pow(x, 4) + p3(X);
+ };
+ auto p5 = [&p4](const TinyVector<1>& X) {
+ const double x = X[0];
+ return 6 * std::pow(x, 5) + p4(X);
+ };
+ auto p6 = [&p5](const TinyVector<1>& X) {
+ const double x = X[0];
+ return 7 * std::pow(x, 6) + p5(X);
+ };
+ auto p7 = [&p6](const TinyVector<1>& X) {
+ const double x = X[0];
+ return 8 * std::pow(x, 7) + p6(X);
+ };
+ auto p8 = [&p7](const TinyVector<1>& X) {
+ const double x = X[0];
+ return 9 * std::pow(x, 8) + p7(X);
+ };
+ auto p9 = [&p8](const TinyVector<1>& X) {
+ const double x = X[0];
+ return 10 * std::pow(x, 9) + p8(X);
+ };
+ auto p10 = [&p9](const TinyVector<1>& X) {
+ const double x = X[0];
+ return 11 * std::pow(x, 10) + p9(X);
+ };
+ auto p11 = [&p10](const TinyVector<1>& X) {
+ const double x = X[0];
+ return 12 * std::pow(x, 11) + p10(X);
+ };
+ auto p12 = [&p11](const TinyVector<1>& X) {
+ const double x = X[0];
+ return 13 * std::pow(x, 12) + p11(X);
+ };
+ auto p13 = [&p12](const TinyVector<1>& X) {
+ const double x = X[0];
+ return 14 * std::pow(x, 13) + p12(X);
+ };
+ auto p14 = [&p13](const TinyVector<1>& X) {
+ const double x = X[0];
+ return 15 * std::pow(x, 14) + p13(X);
+ };
+ auto p15 = [&p14](const TinyVector<1>& X) {
+ const double x = X[0];
+ return 16 * std::pow(x, 15) + p14(X);
+ };
+ auto p16 = [&p15](const TinyVector<1>& X) {
+ const double x = X[0];
+ return 17 * std::pow(x, 16) + p15(X);
+ };
+ auto p17 = [&p16](const TinyVector<1>& X) {
+ const double x = X[0];
+ return 18 * std::pow(x, 17) + p16(X);
+ };
+ auto p18 = [&p17](const TinyVector<1>& X) {
+ const double x = X[0];
+ return 19 * std::pow(x, 18) + p17(X);
+ };
+ auto p19 = [&p18](const TinyVector<1>& X) {
+ const double x = X[0];
+ return 20 * std::pow(x, 19) + p18(X);
+ };
+ auto p20 = [&p19](const TinyVector<1>& X) {
+ const double x = X[0];
+ return 21 * std::pow(x, 20) + p19(X);
+ };
+ auto p21 = [&p20](const TinyVector<1>& X) {
+ const double x = X[0];
+ return 22 * std::pow(x, 21) + p20(X);
+ };
+ auto p22 = [&p21](const TinyVector<1>& X) {
+ const double x = X[0];
+ return 23 * std::pow(x, 22) + p21(X);
+ };
+ auto p23 = [&p22](const TinyVector<1>& X) {
+ const double x = X[0];
+ return 24 * std::pow(x, 23) + p22(X);
+ };
+ auto p24 = [&p23](const TinyVector<1>& X) {
+ const double x = X[0];
+ return 25 * std::pow(x, 24) + p23(X);
+ };
+
+ SECTION("order 0 and 1")
+ {
+ GaussLegendreQuadrature<1> l0(0);
+ GaussLegendreQuadrature<1> l1(1);
+
+ REQUIRE(same_formulas(l0, l1));
+
+ REQUIRE(l0.numberOfPoints() == 1);
+ REQUIRE(l1.numberOfPoints() == 1);
+
+ REQUIRE(integrate(p0, l1, 0.5, 1) == Catch::Approx(0.5));
+ REQUIRE(integrate(p1, l1, 0, 1) == Catch::Approx(2));
+ REQUIRE(integrate(p2, l1, 0, 1) != Catch::Approx(3));
+
+ REQUIRE(get_order(p2, l1, -1, 1, 4) == Catch::Approx(2));
+ }
+
+ SECTION("order 2 and 3")
+ {
+ GaussLegendreQuadrature<1> l2(2);
+ GaussLegendreQuadrature<1> l3(3);
+
+ REQUIRE(same_formulas(l2, l3));
+
+ REQUIRE(l2.numberOfPoints() == 2);
+ REQUIRE(l3.numberOfPoints() == 2);
+
+ REQUIRE(integrate(p0, l3, 0.5, 1) == Catch::Approx(0.5));
+ REQUIRE(integrate(p1, l3, 0, 1) == Catch::Approx(2));
+ REQUIRE(integrate(p2, l3, 0, 1) == Catch::Approx(3));
+ REQUIRE(integrate(p3, l3, 0, 1) == Catch::Approx(4));
+ REQUIRE(integrate(p4, l3, 0, 1) != Catch::Approx(5));
+
+ REQUIRE(get_order(p4, l3, -1, 1, 6) == Catch::Approx(4));
+ }
+
+ SECTION("order 4 and 5")
+ {
+ GaussLegendreQuadrature<1> l4(4);
+ GaussLegendreQuadrature<1> l5(5);
+
+ REQUIRE(same_formulas(l4, l5));
+
+ REQUIRE(l4.numberOfPoints() == 3);
+ REQUIRE(l5.numberOfPoints() == 3);
+
+ REQUIRE(integrate(p0, l5, 0.5, 1) == Catch::Approx(0.5));
+ REQUIRE(integrate(p1, l5, 0, 1) == Catch::Approx(2));
+ REQUIRE(integrate(p2, l5, 0, 1) == Catch::Approx(3));
+ REQUIRE(integrate(p3, l5, 0, 1) == Catch::Approx(4));
+ REQUIRE(integrate(p4, l5, 0, 1) == Catch::Approx(5));
+ REQUIRE(integrate(p5, l5, 0, 1) == Catch::Approx(6));
+ REQUIRE(integrate(p6, l5, 0, 1) != Catch::Approx(7));
+
+ REQUIRE(get_order(p6, l5, -1, 1, 8) == Catch::Approx(6));
+ }
+
+ SECTION("order 6 and 7")
+ {
+ GaussLegendreQuadrature<1> l6(6);
+ GaussLegendreQuadrature<1> l7(7);
+
+ REQUIRE(same_formulas(l6, l7));
+
+ REQUIRE(l6.numberOfPoints() == 4);
+ REQUIRE(l7.numberOfPoints() == 4);
+
+ REQUIRE(integrate(p0, l7, 0.5, 1) == Catch::Approx(0.5));
+ REQUIRE(integrate(p1, l7, 0, 1) == Catch::Approx(2));
+ REQUIRE(integrate(p2, l7, 0, 1) == Catch::Approx(3));
+ REQUIRE(integrate(p3, l7, 0, 1) == Catch::Approx(4));
+ REQUIRE(integrate(p4, l7, 0, 1) == Catch::Approx(5));
+ REQUIRE(integrate(p5, l7, 0, 1) == Catch::Approx(6));
+ REQUIRE(integrate(p6, l7, 0, 1) == Catch::Approx(7));
+ REQUIRE(integrate(p7, l7, 0, 1) == Catch::Approx(8));
+ REQUIRE(integrate(p8, l7, 0, 1) != Catch::Approx(9));
+
+ REQUIRE(get_order(p8, l7, -1, 1, 10) == Catch::Approx(8));
+ }
+
+ SECTION("order 8 and 9")
+ {
+ GaussLegendreQuadrature<1> l8(8);
+ GaussLegendreQuadrature<1> l9(9);
+
+ REQUIRE(same_formulas(l8, l9));
+
+ REQUIRE(l8.numberOfPoints() == 5);
+ REQUIRE(l9.numberOfPoints() == 5);
+
+ REQUIRE(integrate(p0, l9, 0.5, 1) == Catch::Approx(0.5));
+ REQUIRE(integrate(p1, l9, 0, 1) == Catch::Approx(2));
+ REQUIRE(integrate(p2, l9, 0, 1) == Catch::Approx(3));
+ REQUIRE(integrate(p3, l9, 0, 1) == Catch::Approx(4));
+ REQUIRE(integrate(p4, l9, 0, 1) == Catch::Approx(5));
+ REQUIRE(integrate(p5, l9, 0, 1) == Catch::Approx(6));
+ REQUIRE(integrate(p6, l9, 0, 1) == Catch::Approx(7));
+ REQUIRE(integrate(p7, l9, 0, 1) == Catch::Approx(8));
+ REQUIRE(integrate(p8, l9, 0, 1) == Catch::Approx(9));
+ REQUIRE(integrate(p9, l9, 0, 1) == Catch::Approx(10));
+ REQUIRE(integrate(p10, l9, 0, 1) != Catch::Approx(11).epsilon(1E-13));
+
+ REQUIRE(get_order(p10, l9, -1, 1, 12) == Catch::Approx(10));
+ }
+
+ SECTION("order 10 and 11")
+ {
+ GaussLegendreQuadrature<1> l10(10);
+ GaussLegendreQuadrature<1> l11(11);
+
+ REQUIRE(same_formulas(l10, l11));
+
+ REQUIRE(l10.numberOfPoints() == 6);
+ REQUIRE(l11.numberOfPoints() == 6);
+
+ REQUIRE(integrate(p0, l11, 0.5, 1) == Catch::Approx(0.5));
+ REQUIRE(integrate(p1, l11, 0, 1) == Catch::Approx(2));
+ REQUIRE(integrate(p2, l11, 0, 1) == Catch::Approx(3));
+ REQUIRE(integrate(p3, l11, 0, 1) == Catch::Approx(4));
+ REQUIRE(integrate(p4, l11, 0, 1) == Catch::Approx(5));
+ REQUIRE(integrate(p5, l11, 0, 1) == Catch::Approx(6));
+ REQUIRE(integrate(p6, l11, 0, 1) == Catch::Approx(7));
+ REQUIRE(integrate(p7, l11, 0, 1) == Catch::Approx(8));
+ REQUIRE(integrate(p8, l11, 0, 1) == Catch::Approx(9));
+ REQUIRE(integrate(p9, l11, 0, 1) == Catch::Approx(10));
+ REQUIRE(integrate(p10, l11, 0, 1) == Catch::Approx(11));
+ REQUIRE(integrate(p11, l11, 0, 1) == Catch::Approx(12));
+ REQUIRE(integrate(p12, l11, 0, 1) != Catch::Approx(13).epsilon(1E-13));
+
+ REQUIRE(get_order(p12, l11, -1, 1, 14) == Catch::Approx(12));
+ }
+
+ SECTION("order 12 and 13")
+ {
+ GaussLegendreQuadrature<1> l12(12);
+ GaussLegendreQuadrature<1> l13(13);
+
+ REQUIRE(same_formulas(l12, l13));
+
+ REQUIRE(l12.numberOfPoints() == 7);
+ REQUIRE(l13.numberOfPoints() == 7);
+
+ REQUIRE(integrate(p0, l13, 0.5, 1) == Catch::Approx(0.5));
+ REQUIRE(integrate(p1, l13, 0, 1) == Catch::Approx(2));
+ REQUIRE(integrate(p2, l13, 0, 1) == Catch::Approx(3));
+ REQUIRE(integrate(p3, l13, 0, 1) == Catch::Approx(4));
+ REQUIRE(integrate(p4, l13, 0, 1) == Catch::Approx(5));
+ REQUIRE(integrate(p5, l13, 0, 1) == Catch::Approx(6));
+ REQUIRE(integrate(p6, l13, 0, 1) == Catch::Approx(7));
+ REQUIRE(integrate(p7, l13, 0, 1) == Catch::Approx(8));
+ REQUIRE(integrate(p8, l13, 0, 1) == Catch::Approx(9));
+ REQUIRE(integrate(p9, l13, 0, 1) == Catch::Approx(10));
+ REQUIRE(integrate(p10, l13, 0, 1) == Catch::Approx(11));
+ REQUIRE(integrate(p11, l13, 0, 1) == Catch::Approx(12));
+ REQUIRE(integrate(p12, l13, 0, 1) == Catch::Approx(13));
+ REQUIRE(integrate(p13, l13, 0, 1) == Catch::Approx(14));
+ REQUIRE(integrate(p14, l13, 0, 1) != Catch::Approx(15).epsilon(1E-13));
+
+ REQUIRE(get_order(p14, l13, -1, 1, 16) == Catch::Approx(14));
+ }
+
+ SECTION("order 14 and 15")
+ {
+ GaussLegendreQuadrature<1> l14(14);
+ GaussLegendreQuadrature<1> l15(15);
+
+ REQUIRE(same_formulas(l14, l15));
+
+ REQUIRE(l14.numberOfPoints() == 8);
+ REQUIRE(l15.numberOfPoints() == 8);
+
+ REQUIRE(integrate(p0, l15, 0.5, 1) == Catch::Approx(0.5));
+ REQUIRE(integrate(p1, l15, 0, 1) == Catch::Approx(2));
+ REQUIRE(integrate(p2, l15, 0, 1) == Catch::Approx(3));
+ REQUIRE(integrate(p3, l15, 0, 1) == Catch::Approx(4));
+ REQUIRE(integrate(p4, l15, 0, 1) == Catch::Approx(5));
+ REQUIRE(integrate(p5, l15, 0, 1) == Catch::Approx(6));
+ REQUIRE(integrate(p6, l15, 0, 1) == Catch::Approx(7));
+ REQUIRE(integrate(p7, l15, 0, 1) == Catch::Approx(8));
+ REQUIRE(integrate(p8, l15, 0, 1) == Catch::Approx(9));
+ REQUIRE(integrate(p9, l15, 0, 1) == Catch::Approx(10));
+ REQUIRE(integrate(p10, l15, 0, 1) == Catch::Approx(11));
+ REQUIRE(integrate(p11, l15, 0, 1) == Catch::Approx(12));
+ REQUIRE(integrate(p12, l15, 0, 1) == Catch::Approx(13));
+ REQUIRE(integrate(p13, l15, 0, 1) == Catch::Approx(14));
+ REQUIRE(integrate(p14, l15, 0, 1) == Catch::Approx(15));
+ REQUIRE(integrate(p15, l15, 0, 1) == Catch::Approx(16));
+ REQUIRE(integrate(p16, l15, 0, 1) != Catch::Approx(17).epsilon(1E-13));
+
+ REQUIRE(get_order(p16, l15, -1, 1, 18) == Catch::Approx(16));
+ }
+
+ SECTION("order 16 and 17")
+ {
+ GaussLegendreQuadrature<1> l16(16);
+ GaussLegendreQuadrature<1> l17(17);
+
+ REQUIRE(same_formulas(l16, l17));
+
+ REQUIRE(l16.numberOfPoints() == 9);
+ REQUIRE(l17.numberOfPoints() == 9);
+
+ REQUIRE(integrate(p0, l17, 0.5, 1) == Catch::Approx(0.5));
+ REQUIRE(integrate(p1, l17, 0, 1) == Catch::Approx(2));
+ REQUIRE(integrate(p2, l17, 0, 1) == Catch::Approx(3));
+ REQUIRE(integrate(p3, l17, 0, 1) == Catch::Approx(4));
+ REQUIRE(integrate(p4, l17, 0, 1) == Catch::Approx(5));
+ REQUIRE(integrate(p5, l17, 0, 1) == Catch::Approx(6));
+ REQUIRE(integrate(p6, l17, 0, 1) == Catch::Approx(7));
+ REQUIRE(integrate(p7, l17, 0, 1) == Catch::Approx(8));
+ REQUIRE(integrate(p8, l17, 0, 1) == Catch::Approx(9));
+ REQUIRE(integrate(p9, l17, 0, 1) == Catch::Approx(10));
+ REQUIRE(integrate(p10, l17, 0, 1) == Catch::Approx(11));
+ REQUIRE(integrate(p11, l17, 0, 1) == Catch::Approx(12));
+ REQUIRE(integrate(p12, l17, 0, 1) == Catch::Approx(13));
+ REQUIRE(integrate(p13, l17, 0, 1) == Catch::Approx(14));
+ REQUIRE(integrate(p14, l17, 0, 1) == Catch::Approx(15));
+ REQUIRE(integrate(p15, l17, 0, 1) == Catch::Approx(16));
+ REQUIRE(integrate(p16, l17, 0, 1) == Catch::Approx(17));
+ REQUIRE(integrate(p17, l17, 0, 1) == Catch::Approx(18));
+ REQUIRE(integrate(p18, l17, 0, 1) != Catch::Approx(19).epsilon(1E-13));
+
+ REQUIRE(get_order(p18, l17, -1, 1, 20) == Catch::Approx(18));
+ }
+
+ SECTION("order 18 and 19")
+ {
+ GaussLegendreQuadrature<1> l18(18);
+ GaussLegendreQuadrature<1> l19(19);
+
+ REQUIRE(same_formulas(l18, l19));
+
+ REQUIRE(l18.numberOfPoints() == 10);
+ REQUIRE(l19.numberOfPoints() == 10);
+
+ REQUIRE(integrate(p0, l19, 0.5, 1) == Catch::Approx(0.5));
+ REQUIRE(integrate(p1, l19, 0, 1) == Catch::Approx(2));
+ REQUIRE(integrate(p2, l19, 0, 1) == Catch::Approx(3));
+ REQUIRE(integrate(p3, l19, 0, 1) == Catch::Approx(4));
+ REQUIRE(integrate(p4, l19, 0, 1) == Catch::Approx(5));
+ REQUIRE(integrate(p5, l19, 0, 1) == Catch::Approx(6));
+ REQUIRE(integrate(p6, l19, 0, 1) == Catch::Approx(7));
+ REQUIRE(integrate(p7, l19, 0, 1) == Catch::Approx(8));
+ REQUIRE(integrate(p8, l19, 0, 1) == Catch::Approx(9));
+ REQUIRE(integrate(p9, l19, 0, 1) == Catch::Approx(10));
+ REQUIRE(integrate(p10, l19, 0, 1) == Catch::Approx(11));
+ REQUIRE(integrate(p11, l19, 0, 1) == Catch::Approx(12));
+ REQUIRE(integrate(p12, l19, 0, 1) == Catch::Approx(13));
+ REQUIRE(integrate(p13, l19, 0, 1) == Catch::Approx(14));
+ REQUIRE(integrate(p14, l19, 0, 1) == Catch::Approx(15));
+ REQUIRE(integrate(p15, l19, 0, 1) == Catch::Approx(16));
+ REQUIRE(integrate(p16, l19, 0, 1) == Catch::Approx(17));
+ REQUIRE(integrate(p17, l19, 0, 1) == Catch::Approx(18));
+ REQUIRE(integrate(p18, l19, 0, 1) == Catch::Approx(19));
+ REQUIRE(integrate(p19, l19, 0, 1) == Catch::Approx(20));
+ REQUIRE(integrate(p20, l19, 0, 1) != Catch::Approx(21).epsilon(1E-13));
+
+ REQUIRE(get_order(p20, l19, -1, 1, 22) == Catch::Approx(20));
+ }
+
+ SECTION("order 20 and 21")
+ {
+ GaussLegendreQuadrature<1> l20(20);
+ GaussLegendreQuadrature<1> l21(21);
+
+ REQUIRE(same_formulas(l20, l21));
+
+ REQUIRE(l20.numberOfPoints() == 11);
+ REQUIRE(l21.numberOfPoints() == 11);
+
+ REQUIRE(integrate(p0, l21, 0.5, 1) == Catch::Approx(0.5));
+ REQUIRE(integrate(p1, l21, 0, 1) == Catch::Approx(2));
+ REQUIRE(integrate(p2, l21, 0, 1) == Catch::Approx(3));
+ REQUIRE(integrate(p3, l21, 0, 1) == Catch::Approx(4));
+ REQUIRE(integrate(p4, l21, 0, 1) == Catch::Approx(5));
+ REQUIRE(integrate(p5, l21, 0, 1) == Catch::Approx(6));
+ REQUIRE(integrate(p6, l21, 0, 1) == Catch::Approx(7));
+ REQUIRE(integrate(p7, l21, 0, 1) == Catch::Approx(8));
+ REQUIRE(integrate(p8, l21, 0, 1) == Catch::Approx(9));
+ REQUIRE(integrate(p9, l21, 0, 1) == Catch::Approx(10));
+ REQUIRE(integrate(p10, l21, 0, 1) == Catch::Approx(11));
+ REQUIRE(integrate(p11, l21, 0, 1) == Catch::Approx(12));
+ REQUIRE(integrate(p12, l21, 0, 1) == Catch::Approx(13));
+ REQUIRE(integrate(p13, l21, 0, 1) == Catch::Approx(14));
+ REQUIRE(integrate(p14, l21, 0, 1) == Catch::Approx(15));
+ REQUIRE(integrate(p15, l21, 0, 1) == Catch::Approx(16));
+ REQUIRE(integrate(p16, l21, 0, 1) == Catch::Approx(17));
+ REQUIRE(integrate(p17, l21, 0, 1) == Catch::Approx(18));
+ REQUIRE(integrate(p18, l21, 0, 1) == Catch::Approx(19));
+ REQUIRE(integrate(p19, l21, 0, 1) == Catch::Approx(20));
+ REQUIRE(integrate(p20, l21, 0, 1) == Catch::Approx(21));
+ REQUIRE(integrate(p21, l21, 0, 1) == Catch::Approx(22));
+ REQUIRE(integrate(p22, l21, 0, 1) != Catch::Approx(23).epsilon(1E-14));
+
+ REQUIRE(get_order(p22, l21, -1, 1, 24) == Catch::Approx(22).margin(0.02));
+ }
+
+ SECTION("order 22 and 23")
+ {
+ GaussLegendreQuadrature<1> l22(22);
+ GaussLegendreQuadrature<1> l23(23);
+
+ REQUIRE(same_formulas(l22, l23));
+
+ REQUIRE(l22.numberOfPoints() == 12);
+ REQUIRE(l23.numberOfPoints() == 12);
+
+ REQUIRE(integrate(p0, l23, 0.5, 1) == Catch::Approx(0.5));
+ REQUIRE(integrate(p1, l23, 0, 1) == Catch::Approx(2));
+ REQUIRE(integrate(p2, l23, 0, 1) == Catch::Approx(3));
+ REQUIRE(integrate(p3, l23, 0, 1) == Catch::Approx(4));
+ REQUIRE(integrate(p4, l23, 0, 1) == Catch::Approx(5));
+ REQUIRE(integrate(p5, l23, 0, 1) == Catch::Approx(6));
+ REQUIRE(integrate(p6, l23, 0, 1) == Catch::Approx(7));
+ REQUIRE(integrate(p7, l23, 0, 1) == Catch::Approx(8));
+ REQUIRE(integrate(p8, l23, 0, 1) == Catch::Approx(9));
+ REQUIRE(integrate(p9, l23, 0, 1) == Catch::Approx(10));
+ REQUIRE(integrate(p10, l23, 0, 1) == Catch::Approx(11));
+ REQUIRE(integrate(p11, l23, 0, 1) == Catch::Approx(12));
+ REQUIRE(integrate(p12, l23, 0, 1) == Catch::Approx(13));
+ REQUIRE(integrate(p13, l23, 0, 1) == Catch::Approx(14));
+ REQUIRE(integrate(p14, l23, 0, 1) == Catch::Approx(15));
+ REQUIRE(integrate(p15, l23, 0, 1) == Catch::Approx(16));
+ REQUIRE(integrate(p16, l23, 0, 1) == Catch::Approx(17));
+ REQUIRE(integrate(p17, l23, 0, 1) == Catch::Approx(18));
+ REQUIRE(integrate(p18, l23, 0, 1) == Catch::Approx(19));
+ REQUIRE(integrate(p19, l23, 0, 1) == Catch::Approx(20));
+ REQUIRE(integrate(p20, l23, 0, 1) == Catch::Approx(21));
+ REQUIRE(integrate(p21, l23, 0, 1) == Catch::Approx(22));
+ REQUIRE(integrate(p22, l23, 0, 1) == Catch::Approx(23));
+ REQUIRE(integrate(p23, l23, 0, 1) == Catch::Approx(24));
+ REQUIRE(integrate(p24, l23, 0, 1) != Catch::Approx(25).epsilon(1E-15));
+
+ REQUIRE(get_order(p24, l23, -1, 1, 26) == Catch::Approx(24).margin(0.05));
+ }
+
+ SECTION("not implemented formulae")
+ {
+ REQUIRE_THROWS_WITH(GaussLegendreQuadrature<1>(24),
+ "error: Gauss-Legendre quadrature formulae handle orders up to 23");
+ }
+ }
+}
--
GitLab