diff --git a/CMakeLists.txt b/CMakeLists.txt
index e9255bff6ac523e8bfdf61fe42044e4366d3993f..4b302fe373681cdc862851e871c53d3842d0d842 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 16904a7be0cddc9aa54206638ab4b5c759cd0de8..54cae625bc7c899c433bfdf4021ef6242dd5aaf4 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 756fe09b9b4a2428a5649f37f6007ae9404298ca..08642b794df8470eb1e4598c75d643ccdf1de364 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 0000000000000000000000000000000000000000..b2e2f4c2f472639c54e013de1495a33a7433d021
--- /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 0000000000000000000000000000000000000000..e9cc2b9638ee95045444efad116dccad0f117fc7
--- /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 0000000000000000000000000000000000000000..3e4a8b0cc3538f1539464dc2d81d478c1e6b1463
--- /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 0000000000000000000000000000000000000000..c92d3187246268ec61cf15356d57f7ea39111df8
--- /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 5cc4fee50ca3e8d71899a58083f1c251fa30068c..b0b2b73281d3c7473cdd0636bd2c1490fd21329c 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 7893ba07eec45c4b33c842454172dde6446e46ac..2e290556eac3a8d46d76cd9fd114cfa025e6ca2f 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 0000000000000000000000000000000000000000..f2001356081242c33c0412c7d56c174f296f71aa
--- /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");
+ }
+ }
+}