diff --git a/src/algebra/IntegrationTools.hpp b/src/algebra/IntegrationTools.hpp
index 08642b794df8470eb1e4598c75d643ccdf1de364..954924d6ec9f6db78127d3e6afe4fae28db6cb3b 100644
--- a/src/algebra/IntegrationTools.hpp
+++ b/src/algebra/IntegrationTools.hpp
@@ -8,6 +8,7 @@
#include <utils/Types.hpp>
#include <analysis/GaussLegendreQuadrature.hpp>
+#include <analysis/GaussLobattoQuadrature.hpp>
#include <cmath>
@@ -21,337 +22,13 @@ enum class QuadratureType : int8_t
QT__end
};
-class IntegrationMethod
-{
- public:
- PUGS_INLINE virtual SmallArray<TinyVector<1>> quadraturePoints(const double& a, const double& b) = 0;
- PUGS_INLINE virtual size_t numberOfPoints() = 0;
- PUGS_INLINE virtual SmallArray<double> weights() = 0;
- // One does not use the '=default' constructor to avoid
- // (zero-initialization) performances issues
- PUGS_INLINE
- constexpr IntegrationMethod() noexcept {}
-
- PUGS_INLINE
- constexpr IntegrationMethod(const IntegrationMethod&) noexcept = default;
-
- PUGS_INLINE
- constexpr IntegrationMethod(IntegrationMethod&& v) noexcept = default;
-
- virtual PUGS_INLINE ~IntegrationMethod() noexcept = default;
-};
-
-template <size_t Order>
-class IntegrationMethodLobatto : public IntegrationMethod
-{
- private:
- static constexpr size_t number_points = (Order < 4 ? 3 : (Order < 6 ? 4 : (Order < 8 ? 5 : Order < 10 ? 6 : 7)));
- SmallArray<double> m_weights;
- SmallArray<TinyVector<1>> m_quadrature_positions;
- PUGS_INLINE SmallArray<TinyVector<1>>
- translateQuadraturePoints(const double& a, const double& b)
- {
- SmallArray<TinyVector<1>> points{number_points};
- double length = b - a;
- for (size_t i = 0; i < m_quadrature_positions.size(); i++) {
- TinyVector<1> qpos{m_quadrature_positions[i]};
- points[i] = 0.5 * (length * qpos + a + b);
- }
- return points;
- }
-
- PUGS_INLINE
- void
- fillArrayLobatto()
- {
- switch (Order) {
- case 0:
- case 1:
- case 2:
- case 3: {
- double oneov3 = 1. / 3.;
- m_quadrature_positions[0] = -1.;
- m_quadrature_positions[1] = 0.;
- m_quadrature_positions[2] = 1.;
- m_weights[0] = oneov3;
- m_weights[1] = 4 * oneov3;
- m_weights[2] = oneov3;
- // m_quadrature_positions.fill({-1, 0, 1});
- // m_weights.fill({oneov3, 4 * oneov3, oneov3});
- break;
- }
- case 4:
- case 5: {
- double oneov6 = 1. / 6.;
- double coef = std::sqrt(1. / 5.);
- m_quadrature_positions[0] = -1;
- m_quadrature_positions[1] = -coef;
- m_quadrature_positions[2] = coef;
- m_quadrature_positions[3] = 1;
- m_weights[0] = oneov6;
- m_weights[1] = 5 * oneov6;
- m_weights[2] = 5 * oneov6;
- m_weights[3] = oneov6;
- break;
- }
- case 6:
- case 7: {
- double oneov90 = 1. / 90.;
- double oneov10 = 1. / 10.;
- double coef = std::sqrt(3. / 7.);
- m_quadrature_positions[0] = -1;
- m_quadrature_positions[1] = -coef;
- m_quadrature_positions[2] = 0;
- m_quadrature_positions[3] = coef;
- m_quadrature_positions[4] = 1;
- m_weights[0] = oneov10;
- m_weights[1] = 49 * oneov90;
- m_weights[2] = 64 * oneov90;
- m_weights[3] = 49 * oneov90;
- m_weights[4] = oneov10;
- break;
- }
- case 8:
- case 9: {
- double oneov30 = 1. / 30.;
- double coef1 = std::sqrt(1. / 3. - 2 * std::sqrt(7) / 21.);
- double coef2 = std::sqrt(1. / 3. + 2 * std::sqrt(7) / 21.);
- double coef3 = (14. + std::sqrt(7)) * oneov30;
- double coef4 = (14. - std::sqrt(7)) * oneov30;
- m_quadrature_positions[0] = -1;
- m_quadrature_positions[1] = -coef2;
- m_quadrature_positions[2] = -coef1;
- m_quadrature_positions[3] = coef1;
- m_quadrature_positions[4] = coef2;
- m_quadrature_positions[5] = 1;
- m_weights[0] = 2. * oneov30;
- m_weights[1] = coef4;
- m_weights[2] = coef3;
- m_weights[3] = coef3;
- m_weights[4] = coef4;
- m_weights[5] = 2. * oneov30;
- break;
- }
- case 10:
- case 11: {
- double oneov350 = 1. / 350.;
- double oneov21 = 1. / 21.;
- double coef0 = 256. / 525.;
- double coef1 = std::sqrt((5. - 2 * std::sqrt(5. / 3.)) / 11.);
- double coef2 = std::sqrt((5. + 2 * std::sqrt(5. / 3.)) / 11.);
- double coef3 = (124. + 7 * std::sqrt(15.)) * oneov350;
- double coef4 = (124. - 7 * std::sqrt(15.)) * oneov350;
- m_quadrature_positions[0] = -1;
- m_quadrature_positions[1] = -coef2;
- m_quadrature_positions[2] = -coef1;
- m_quadrature_positions[3] = 0;
- m_quadrature_positions[4] = coef1;
- m_quadrature_positions[5] = coef2;
- m_quadrature_positions[6] = 1;
- m_weights[0] = oneov21;
- m_weights[1] = coef4;
- m_weights[2] = coef3;
- m_weights[3] = coef0;
- m_weights[4] = coef3;
- m_weights[5] = coef4;
- m_weights[6] = oneov21;
- break;
- }
-
- default: {
- throw UnexpectedError("Gauss-Lobatto quadratures handle orders up to 11.");
- break;
- }
- }
- }
-
- public:
- // One does not use the '=default' constructor to avoid
- // (zero-initialization) performances issues
- PUGS_INLINE
- constexpr IntegrationMethodLobatto() noexcept : m_weights{number_points}, m_quadrature_positions{number_points}
- {
- m_weights.fill(0.);
- m_quadrature_positions.fill(TinyVector<1, double>{0});
- fillArrayLobatto();
- }
-
- PUGS_INLINE
- SmallArray<TinyVector<1>>
- quadraturePoints(const double& a, const double& b)
- {
- return translateQuadraturePoints(a, b);
- }
-
- PUGS_INLINE size_t
- numberOfPoints()
- {
- return number_points;
- }
- PUGS_INLINE
- SmallArray<double>
- weights()
- {
- return m_weights;
- }
-
- PUGS_INLINE
- constexpr IntegrationMethodLobatto(const IntegrationMethodLobatto&) noexcept = default;
-
- PUGS_INLINE
- constexpr IntegrationMethodLobatto(IntegrationMethodLobatto&& v) noexcept = default;
-
- PUGS_INLINE
- ~IntegrationMethodLobatto() noexcept = default;
-};
-
-template <size_t Order>
-class IntegrationMethodLegendre : public IntegrationMethod
-{
- private:
- static constexpr size_t number_points = (Order < 4 ? 2 : (Order < 6 ? 3 : (Order < 8 ? 4 : 5)));
- SmallArray<double> m_weights;
- SmallArray<TinyVector<1>> m_quadrature_positions;
- PUGS_INLINE SmallArray<TinyVector<1>>
- translateQuadraturePoints(const double& a, const double& b)
- {
- SmallArray<TinyVector<1>> points{number_points};
- double length = b - a;
- for (size_t i = 0; i < m_quadrature_positions.size(); i++) {
- TinyVector<1> qpos{m_quadrature_positions[i]};
- points[i] = 0.5 * (length * qpos + a + b);
- }
- return points;
- }
-
- PUGS_INLINE
- void
- fillArrayLegendre()
- {
- 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) {
- // case 0:
- // case 1:
- // case 2:
- // case 3: {
- // double oneovsqrt3 = 1. / std::sqrt(3);
- // m_quadrature_positions[0] = -oneovsqrt3;
- // m_quadrature_positions[1] = oneovsqrt3;
- // m_weights[0] = 1;
- // m_weights[1] = 1;
- // // m_quadrature_positions.fill({-1, 0, 1});
- // // m_weights.fill({oneov3, 4 * oneov3, oneov3});
- // break;
- // }
- // case 4:
- // case 5: {
- // double oneover9 = 1. / 9;
- // double coef = std::sqrt(3. / 5.);
- // m_quadrature_positions[0] = -coef;
- // m_quadrature_positions[1] = 0;
- // m_quadrature_positions[2] = coef;
- // m_weights[0] = 5 * oneover9;
- // m_weights[1] = 8 * oneover9;
- // m_weights[2] = 5 * oneover9;
- // break;
- // }
- // case 6:
- // case 7: {
- // double oneov36 = 1. / 36.;
- // double coef1 = std::sqrt(3. / 7. - 2. / 7. * std::sqrt(6. / 5.));
- // double coef2 = std::sqrt(3. / 7. + 2. / 7. * std::sqrt(6. / 5.));
- // double coef3 = 18. + std::sqrt(30.);
- // double coef4 = 18. - std::sqrt(30.);
-
- // m_quadrature_positions[0] = -coef2;
- // m_quadrature_positions[1] = -coef1;
- // m_quadrature_positions[2] = coef1;
- // m_quadrature_positions[3] = coef2;
- // m_weights[0] = coef4 * oneov36;
- // m_weights[1] = coef3 * oneov36;
- // m_weights[2] = coef3 * oneov36;
- // m_weights[3] = coef4 * oneov36;
- // break;
- // }
- // case 8:
- // case 9: {
- // double oneov900 = 1. / 900.;
- // double oneov225 = 1. / 225.;
- // double coef1 = 1. / 3. * std::sqrt(5. - 2. * std::sqrt(10. / 7.));
- // double coef2 = 1. / 3. * std::sqrt(5. + 2. * std::sqrt(10. / 7.));
- // double coef3 = 322. + 13. * std::sqrt(70.);
- // double coef4 = 322. - 13. * std::sqrt(70.);
- // m_quadrature_positions[0] = -coef2;
- // m_quadrature_positions[1] = -coef1;
- // m_quadrature_positions[2] = 0;
- // m_quadrature_positions[3] = coef1;
- // m_quadrature_positions[4] = coef2;
- // m_weights[0] = coef4 * oneov900;
- // m_weights[1] = coef3 * oneov900;
- // m_weights[2] = 128. * oneov225;
- // m_weights[3] = coef3 * oneov900;
- // m_weights[4] = coef4 * oneov900;
- // break;
- // }
- // default: {
- // throw UnexpectedError("Gauss-Legendre quadratures handle orders up to 9.");
- // break;
- // }
- // }
- }
-
- public:
- // One does not use the '=default' constructor to avoid
- // (zero-initialization) performances issues
- PUGS_INLINE
- constexpr IntegrationMethodLegendre() noexcept : m_weights{number_points}, m_quadrature_positions{number_points}
- {
- m_weights.fill(0.);
- m_quadrature_positions.fill(TinyVector<1, double>{0});
- fillArrayLegendre();
- }
-
- PUGS_INLINE
- SmallArray<TinyVector<1>>
- quadraturePoints(const double& a, const double& b)
- {
- return translateQuadraturePoints(a, b);
- }
-
- PUGS_INLINE size_t
- numberOfPoints()
- {
- return number_points;
- }
- PUGS_INLINE
- SmallArray<double>
- weights()
- {
- return m_weights;
- }
-
- PUGS_INLINE
- constexpr IntegrationMethodLegendre(const IntegrationMethodLegendre&) noexcept = default;
-
- PUGS_INLINE
- constexpr IntegrationMethodLegendre(IntegrationMethodLegendre&& v) noexcept = default;
-
- PUGS_INLINE
- ~IntegrationMethodLegendre() noexcept = default;
-};
-
template <size_t Order, typename T = double>
class IntegrationTools
{
public:
using data_type = T;
-
- private:
- // T m_values[N];
- std::shared_ptr<IntegrationMethod> m_integration_method;
- size_t N;
+ SmallArray<const double> m_weights;
+ SmallArray<const TinyVector<1>> m_quadrature_positions;
public:
PUGS_INLINE T
@@ -359,7 +36,7 @@ class IntegrationTools
{
T result = 0;
SmallArray<const TinyVector<1>> positions = quadraturePoints(a, b);
- SmallArray<double> weights = m_integration_method->weights();
+ SmallArray<const double> weights = m_weights;
SmallArray<T> values(positions.size());
EvaluateAtPoints<double(const TinyVector<1>)>::evaluateTo(function, positions, values);
Assert(positions.size() == weights.size(), "Wrong number of quadrature points and/or weights in Gauss quadrature");
@@ -375,7 +52,7 @@ class IntegrationTools
SmallArray<const TinyVector<1>> positions = quadraturePoints(vertices);
SmallArray<T> result{vertices.size() - 1};
SmallArray<double> interval_size{vertices.size() - 1};
- SmallArray<double> weights = m_integration_method->weights();
+ SmallArray<const double> weights = m_weights;
SmallArray<T> values(positions.size());
EvaluateAtPoints<double(const TinyVector<1>)>::evaluateTo(function, positions, values);
@@ -393,10 +70,10 @@ class IntegrationTools
PUGS_INLINE T
testIntegrate(const double& a, const double& b) const
{
- T result = 0;
- SmallArray<TinyVector<1>> positions = quadraturePoints(a, b);
+ T result = 0;
+ SmallArray<const TinyVector<1>> positions = quadraturePoints(a, b);
- SmallArray<double> weights = m_integration_method->weights();
+ SmallArray<const double> weights = m_weights;
SmallArray<T> values(weights.size());
for (size_t j = 0; j < weights.size(); ++j) {
@@ -414,9 +91,9 @@ class IntegrationTools
PUGS_INLINE SmallArray<T>
testIntegrateFunction(const SmallArray<TinyVector<1>>& vertices) const
{
- SmallArray<TinyVector<1>> positions = quadraturePoints(vertices);
+ SmallArray<const TinyVector<1>> positions = quadraturePoints(vertices);
SmallArray<double> interval_size{vertices.size() - 1};
- SmallArray<double> weights = m_integration_method->weights();
+ SmallArray<const double> weights = this->weights();
SmallArray<T> values(positions.size());
for (size_t j = 0; j < positions.size(); ++j) {
values[j] = std::pow(positions[j][0], 3);
@@ -441,10 +118,16 @@ class IntegrationTools
return result;
}
- PUGS_INLINE SmallArray<TinyVector<1>>
+ PUGS_INLINE SmallArray<const TinyVector<1>>
quadraturePoints(const double& a, const double& b) const
{
- return m_integration_method->quadraturePoints(a, b);
+ SmallArray<TinyVector<1>> points{m_quadrature_positions.size()};
+ double length = b - a;
+ for (size_t i = 0; i < m_quadrature_positions.size(); i++) {
+ TinyVector<1> qpos{m_quadrature_positions[i]};
+ points[i] = 0.5 * (length * qpos + a + b);
+ }
+ return points;
}
PUGS_INLINE
@@ -452,12 +135,12 @@ class IntegrationTools
quadraturePoints(const SmallArray<TinyVector<1>>& positions) const
{
size_t number_of_intervals = positions.size() - 1;
- size_t quadrature_size = m_integration_method->numberOfPoints();
+ size_t quadrature_size = m_quadrature_positions.size();
SmallArray<TinyVector<1>> quadratures{quadrature_size * number_of_intervals};
for (size_t j = 0; j < number_of_intervals; ++j) {
double a = positions[j][0];
double b = positions[j + 1][0];
- auto intermediate_positions = m_integration_method->quadraturePoints(a, b);
+ auto intermediate_positions = this->quadraturePoints(a, b);
for (size_t k = 0; k < quadrature_size; k++) {
quadratures[j * quadrature_size + k] = intermediate_positions[k];
}
@@ -465,10 +148,10 @@ class IntegrationTools
return quadratures;
}
- PUGS_INLINE SmallArray<double>
+ PUGS_INLINE SmallArray<const double>
weights() const
{
- return m_integration_method->weights();
+ return m_weights;
}
PUGS_INLINE
@@ -478,15 +161,20 @@ class IntegrationTools
return Order;
}
- PUGS_INLINE constexpr IntegrationTools(QuadratureType quadrature)
+ PUGS_INLINE
+ IntegrationTools(QuadratureType quadrature)
{
switch (quadrature) {
case QuadratureType::gausslobatto: {
- m_integration_method = std::make_shared<IntegrationMethodLobatto<Order>>();
+ GaussLobattoQuadrature<1> tensorial_gauss_lobatto(Order);
+ m_quadrature_positions = tensorial_gauss_lobatto.pointList();
+ m_weights = tensorial_gauss_lobatto.weightList();
break;
}
case QuadratureType::gausslegendre: {
- m_integration_method = std::make_shared<IntegrationMethodLegendre<Order>>();
+ GaussLegendreQuadrature<1> tensorial_gauss_legendre(Order);
+ m_quadrature_positions = tensorial_gauss_legendre.pointList();
+ m_weights = tensorial_gauss_legendre.weightList();
break;
}
case QuadratureType::QT__end: {
diff --git a/src/analysis/CMakeLists.txt b/src/analysis/CMakeLists.txt
index b2e2f4c2f472639c54e013de1495a33a7433d021..03c39ec899321a01c4ec569a44bd3c3791e39589 100644
--- a/src/analysis/CMakeLists.txt
+++ b/src/analysis/CMakeLists.txt
@@ -2,4 +2,5 @@
add_library(
PugsAnalysis
- GaussLegendreQuadrature.cpp)
+ GaussLegendreQuadrature.cpp
+ GaussLobattoQuadrature.cpp)
diff --git a/src/analysis/GaussLobattoQuadrature.cpp b/src/analysis/GaussLobattoQuadrature.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..267b91c36c2df63bdb812d307777eef5f47ff9d2
--- /dev/null
+++ b/src/analysis/GaussLobattoQuadrature.cpp
@@ -0,0 +1,180 @@
+#include <analysis/GaussLobattoQuadrature.hpp>
+#include <utils/Exceptions.hpp>
+
+template <>
+void
+GaussLobattoQuadrature<1>::_buildPointAndWeightLists()
+{
+ const size_t nb_points = m_order / 2 + 2;
+
+ SmallArray<TinyVector<1>> point_list(nb_points);
+ SmallArray<double> weight_list(nb_points);
+
+ switch (m_order) {
+ case 0:
+ case 1: {
+ point_list[0] = -1;
+ point_list[1] = +1;
+
+ weight_list[0] = 1;
+ weight_list[1] = 1;
+ break;
+ }
+ case 2:
+ case 3: {
+ point_list[0] = -1;
+ point_list[1] = 0;
+ point_list[2] = +1;
+
+ weight_list[0] = 0.3333333333333333333333333;
+ weight_list[1] = 1.3333333333333333333333333;
+ weight_list[2] = 0.3333333333333333333333333;
+ break;
+ }
+ case 4:
+ case 5: {
+ point_list[0] = -1;
+ point_list[1] = -0.4472135954999579392818347;
+ point_list[2] = +0.4472135954999579392818347;
+ point_list[3] = +1;
+
+ weight_list[0] = 0.1666666666666666666666667;
+ weight_list[1] = 0.8333333333333333333333333;
+ weight_list[2] = 0.8333333333333333333333333;
+ weight_list[3] = 0.1666666666666666666666667;
+ break;
+ }
+ case 6:
+ case 7: {
+ point_list[0] = -1;
+ point_list[1] = -0.6546536707079771437982925;
+ point_list[2] = 0;
+ point_list[3] = +0.6546536707079771437982925;
+ point_list[4] = +1;
+
+ weight_list[0] = 0.1;
+ weight_list[1] = 0.5444444444444444444444444;
+ weight_list[2] = 0.7111111111111111111111111;
+ weight_list[3] = 0.5444444444444444444444444;
+ weight_list[4] = 0.1;
+ break;
+ }
+ case 8:
+ case 9: {
+ point_list[0] = -1;
+ point_list[1] = -0.7650553239294646928510030;
+ point_list[2] = -0.2852315164806450963141510;
+ point_list[3] = +0.2852315164806450963141510;
+ point_list[4] = +0.7650553239294646928510030;
+ point_list[5] = +1;
+
+ weight_list[0] = 0.0666666666666666666666667;
+ weight_list[1] = 0.3784749562978469803166128;
+ weight_list[2] = 0.5548583770354863530167205;
+ weight_list[3] = 0.5548583770354863530167205;
+ weight_list[4] = 0.3784749562978469803166128;
+ weight_list[5] = 0.0666666666666666666666667;
+ break;
+ }
+ case 10:
+ case 11: {
+ point_list[0] = -1;
+ point_list[1] = -0.8302238962785669298720322;
+ point_list[2] = -0.4688487934707142138037719;
+ point_list[3] = 0;
+ point_list[4] = +0.4688487934707142138037719;
+ point_list[5] = +0.8302238962785669298720322;
+ point_list[6] = +1;
+
+ weight_list[0] = 0.0476190476190476190476190;
+ weight_list[1] = 0.2768260473615659480107004;
+ weight_list[2] = 0.4317453812098626234178710;
+ weight_list[3] = 0.4876190476190476190476190;
+ weight_list[4] = 0.4317453812098626234178710;
+ weight_list[5] = 0.2768260473615659480107004;
+ weight_list[6] = 0.0476190476190476190476190;
+ break;
+ }
+ case 12:
+ case 13: {
+ point_list[0] = -1;
+ point_list[1] = -0.8717401485096066153374457;
+ point_list[2] = -0.5917001814331423021445107;
+ point_list[3] = -0.2092992179024788687686573;
+ point_list[4] = +0.2092992179024788687686573;
+ point_list[5] = +0.5917001814331423021445107;
+ point_list[6] = +0.8717401485096066153374457;
+ point_list[7] = +1;
+
+ weight_list[0] = 0.0357142857142857142857143;
+ weight_list[1] = 0.2107042271435060393829921;
+ weight_list[2] = 0.3411226924835043647642407;
+ weight_list[3] = 0.4124587946587038815670530;
+ weight_list[4] = 0.4124587946587038815670530;
+ weight_list[5] = 0.3411226924835043647642407;
+ weight_list[6] = 0.2107042271435060393829921;
+ weight_list[7] = 0.0357142857142857142857143;
+ break;
+ }
+ default: {
+ throw NormalError("Gauss-Lobatto quadrature formulae handle orders up to 13");
+ }
+ }
+
+ m_point_list = point_list;
+ m_weight_list = weight_list;
+}
+
+template <>
+void
+GaussLobattoQuadrature<2>::_buildPointAndWeightLists()
+{
+ const size_t nb_points_1d = this->m_order / 2 + 2;
+ const size_t nb_points = nb_points_1d * nb_points_1d;
+
+ SmallArray<TinyVector<2>> point_list(nb_points);
+ SmallArray<double> weight_list(nb_points);
+
+ GaussLobattoQuadrature<1> gauss_lobatto_1d(m_order);
+ const auto& point_list_1d = gauss_lobatto_1d.pointList();
+ const auto& weight_list_1d = gauss_lobatto_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
+GaussLobattoQuadrature<3>::_buildPointAndWeightLists()
+{
+ const size_t nb_points_1d = this->m_order / 2 + 2;
+ 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);
+
+ GaussLobattoQuadrature<1> gauss_lobatto_1d(m_order);
+ const auto& point_list_1d = gauss_lobatto_1d.pointList();
+ const auto& weight_list_1d = gauss_lobatto_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/GaussLobattoQuadrature.hpp b/src/analysis/GaussLobattoQuadrature.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..89dbb2e8b11da89bc3c326fd8a39e7a3412c2e4f
--- /dev/null
+++ b/src/analysis/GaussLobattoQuadrature.hpp
@@ -0,0 +1,25 @@
+#ifndef GAUSS_LOBATTO_QUADRATURE_HPP
+#define GAUSS_LOBATTO_QUADRATURE_HPP
+
+#include <analysis/QuadratureFormula.hpp>
+
+template <size_t Dimension>
+class GaussLobattoQuadrature final : public QuadratureForumla<Dimension>
+{
+ private:
+ void _buildPointAndWeightLists();
+
+ public:
+ GaussLobattoQuadrature(GaussLobattoQuadrature&&) = default;
+ GaussLobattoQuadrature(const GaussLobattoQuadrature&) = default;
+
+ explicit GaussLobattoQuadrature(const size_t order) : QuadratureForumla<Dimension>(order)
+ {
+ this->_buildPointAndWeightLists();
+ }
+
+ GaussLobattoQuadrature() = delete;
+ ~GaussLobattoQuadrature() = default;
+};
+
+#endif // GAUSS_LOBATTO_QUADRATURE_HPP
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index b0b2b73281d3c7473cdd0636bd2c1490fd21329c..b5d57ecad09204171ec829e250edecf565740e48 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -81,6 +81,7 @@ add_executable (unit_tests
test_FunctionSymbolId.cpp
test_FunctionTable.cpp
test_GaussLegendreQuadrature.cpp
+ test_GaussLobattoQuadrature.cpp
test_IfProcessor.cpp
test_IncDecExpressionProcessor.cpp
test_INodeProcessor.cpp
diff --git a/tests/test_FunctionEvaluation.cpp b/tests/test_FunctionEvaluation.cpp
index 2e290556eac3a8d46d76cd9fd114cfa025e6ca2f..9a8839f5aada5171b4071d7cd787c1ddece8118c 100644
--- a/tests/test_FunctionEvaluation.cpp
+++ b/tests/test_FunctionEvaluation.cpp
@@ -30,14 +30,14 @@ TEST_CASE("IntegrationTools", "[integration]")
}
SECTION("QuadraturePoints")
{
- SmallArray<TinyVector<1>> quadratures = integrationtools.quadraturePoints(a, b);
+ SmallArray<const TinyVector<1>> quadratures = integrationtools.quadraturePoints(a, b);
REQUIRE(quadratures[0] == 1);
REQUIRE(quadratures[1] == 1.5);
REQUIRE(quadratures[2] == 2);
}
SECTION("Weights")
{
- SmallArray<double> weights = integrationtools.weights();
+ SmallArray<const double> weights = integrationtools.weights();
REQUIRE(((weights[0] == 1. / 3.) and (weights[1] == 4. / 3.) and (weights[2] == 1. / 3.)));
}
SECTION("TestOrder3")
diff --git a/tests/test_GaussLobattoQuadrature.cpp b/tests/test_GaussLobattoQuadrature.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..d60a0653b240b1c45077580d74a2f56e85b8da4a
--- /dev/null
+++ b/tests/test_GaussLobattoQuadrature.cpp
@@ -0,0 +1,409 @@
+#include <catch2/catch_approx.hpp>
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_all.hpp>
+
+#include <analysis/GaussLobattoQuadrature.hpp>
+#include <utils/Exceptions.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("GaussLobattoQuadrature", "[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 is_symmetric_formula = [](auto quadrature_formula) {
+ auto point_list = quadrature_formula.pointList();
+ auto weight_list = quadrature_formula.weightList();
+
+ bool is_symmetric = true;
+
+ const size_t middle_index = point_list.size() / 2;
+ for (size_t i = 0; i <= middle_index; ++i) {
+ if (point_list[i] != -point_list[point_list.size() - 1 - i]) {
+ return false;
+ }
+ }
+ for (size_t i = 0; i <= middle_index; ++i) {
+ if (weight_list[i] != weight_list[point_list.size() - 1 - i]) {
+ return false;
+ }
+ }
+
+ return is_symmetric;
+ };
+
+ 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);
+ };
+
+ SECTION("order 0 and 1")
+ {
+ GaussLobattoQuadrature<1> l0(0);
+ GaussLobattoQuadrature<1> l1(1);
+
+ REQUIRE(same_formulas(l0, l1));
+ REQUIRE(is_symmetric_formula(l1));
+
+ REQUIRE(l1.numberOfPoints() == 2);
+
+ 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")
+ {
+ GaussLobattoQuadrature<1> l2(2);
+ GaussLobattoQuadrature<1> l3(3);
+
+ REQUIRE(same_formulas(l2, l3));
+ REQUIRE(is_symmetric_formula(l3));
+
+ REQUIRE(l3.numberOfPoints() == 3);
+
+ 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")
+ {
+ GaussLobattoQuadrature<1> l4(4);
+ GaussLobattoQuadrature<1> l5(5);
+
+ REQUIRE(same_formulas(l4, l5));
+ REQUIRE(is_symmetric_formula(l5));
+
+ REQUIRE(l5.numberOfPoints() == 4);
+
+ 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")
+ {
+ GaussLobattoQuadrature<1> l6(6);
+ GaussLobattoQuadrature<1> l7(7);
+
+ REQUIRE(same_formulas(l6, l7));
+ REQUIRE(is_symmetric_formula(l7));
+
+ REQUIRE(l7.numberOfPoints() == 5);
+
+ 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")
+ {
+ GaussLobattoQuadrature<1> l8(8);
+ GaussLobattoQuadrature<1> l9(9);
+
+ REQUIRE(same_formulas(l8, l9));
+ REQUIRE(is_symmetric_formula(l9));
+
+ REQUIRE(l9.numberOfPoints() == 6);
+
+ 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")
+ {
+ GaussLobattoQuadrature<1> l10(10);
+ GaussLobattoQuadrature<1> l11(11);
+
+ REQUIRE(same_formulas(l10, l11));
+ REQUIRE(is_symmetric_formula(l11));
+
+ REQUIRE(l11.numberOfPoints() == 7);
+
+ 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")
+ {
+ GaussLobattoQuadrature<1> l12(12);
+ GaussLobattoQuadrature<1> l13(13);
+
+ REQUIRE(same_formulas(l12, l13));
+ REQUIRE(is_symmetric_formula(l13));
+
+ REQUIRE(l13.numberOfPoints() == 8);
+
+ 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("not implemented formulae")
+ {
+ REQUIRE_THROWS_WITH(GaussLobattoQuadrature<1>(14),
+ "error: Gauss-Lobatto quadrature formulae handle orders up to 13");
+ }
+ }
+
+ SECTION("2D")
+ {
+ auto integrate = [](auto f, auto quadrature_formula, const double xa, const double xb, const double ya,
+ const double yb) {
+ auto point_list = quadrature_formula.pointList();
+ auto weight_list = quadrature_formula.weightList();
+
+ const double alphax = 0.5 * (xb - xa);
+ const double betax = 0.5 * (xa + xb);
+
+ const double alphay = 0.5 * (yb - ya);
+ const double betay = 0.5 * (ya + yb);
+
+ auto x = [&alphax, &betax, &alphay, &betay](auto x_hat) {
+ return TinyVector<2>{alphax * x_hat[0] + betax, alphay * x_hat[1] + betay};
+ };
+
+ 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 alphax * alphay * value;
+ };
+
+ auto p6 = [](const double x) { return 7 * std::pow(x, 6); };
+ auto p7 = [](const double x) { return 8 * std::pow(x, 7); };
+
+ auto px7y6 = [&p6, &p7](const TinyVector<2>& X) {
+ const double x = X[0];
+ const double y = X[1];
+ return p7(x) * p6(y);
+ };
+
+ SECTION("order 6 and 7")
+ {
+ GaussLobattoQuadrature<2> l6(6);
+ GaussLobattoQuadrature<2> l7(7);
+
+ REQUIRE(same_formulas(l6, l7));
+
+ REQUIRE(l7.numberOfPoints() == 5 * 5);
+
+ REQUIRE(integrate(px7y6, l7, 0, 1, 0.2, 0.8) == Catch::Approx(std::pow(0.8, 7) - std::pow(0.2, 7)));
+ }
+ }
+
+ SECTION("3D")
+ {
+ auto integrate = [](auto f, auto quadrature_formula, const double xa, const double xb, const double ya,
+ const double yb, const double za, const double zb) {
+ auto point_list = quadrature_formula.pointList();
+ auto weight_list = quadrature_formula.weightList();
+
+ const double alphax = 0.5 * (xb - xa);
+ const double betax = 0.5 * (xa + xb);
+
+ const double alphay = 0.5 * (yb - ya);
+ const double betay = 0.5 * (ya + yb);
+
+ const double alphaz = 0.5 * (zb - za);
+ const double betaz = 0.5 * (za + zb);
+
+ auto x = [&alphax, &betax, &alphay, &betay, &alphaz, &betaz](auto x_hat) {
+ return TinyVector<3>{alphax * x_hat[0] + betax, alphay * x_hat[1] + betay, alphaz * x_hat[2] + betaz};
+ };
+
+ 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 alphax * alphay * alphaz * value;
+ };
+
+ auto p6 = [](const double x) { return 7 * std::pow(x, 6); };
+ auto p7 = [](const double x) { return 8 * std::pow(x, 7); };
+
+ auto px7y6 = [&p6, &p7](const TinyVector<3>& X) {
+ const double x = X[0];
+ const double y = X[1];
+ const double z = X[2];
+ return p7(x) * p6(y) + p7(z);
+ };
+
+ SECTION("order 6 and 7")
+ {
+ GaussLobattoQuadrature<3> l6(6);
+ GaussLobattoQuadrature<3> l7(7);
+
+ REQUIRE(same_formulas(l6, l7));
+
+ REQUIRE(l7.numberOfPoints() == 5 * 5 * 5);
+
+ REQUIRE(integrate(px7y6, l7, 0, 1, 0.2, 0.8, -0.1, 0.7) ==
+ Catch::Approx((std::pow(0.8, 7) - std::pow(0.2, 7)) * (0.7 - -0.1) +
+ (0.8 - 0.2) * (std::pow(0.7, 8) - std::pow(-0.1, 8))));
+ }
+ }
+}