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

Begin quadrature formula rework 

Add bunch of Gauss-Legendre formula and begin testing them
parent 3cd6a765
Branches
Tags
1 merge request!124Add files for high order integration with quadratures
Show whitespace changes
Inline Side-by-side
src/algebra/IntegrationTools.hpp
+ 479
− 125
View file @ 546807e7
#ifndef INTEGRATION_TOOLS_HPP
#define INTEGRATION_TOOLS_HPP
#include <iostream>
#include <language/utils/EvaluateAtPoints.hpp>
#include <utils/Exceptions.hpp>
#include <utils/PugsAssert.hpp>
#include <utils/PugsMacros.hpp>
#include <utils/Types.hpp>
#include <cmath>
#include <language/utils/EvaluateAtPoints.hpp>
#include <utils/Types.hpp>
enum class QuadratureType : int8_t
{
......@@ -21,15 +19,362 @@ enum class QuadratureType : int8_t
QT__end
};
template <size_t Order>
class IntegrationMethod
template <size_t Dimension>
class QuadratureBase
{
protected:
size_t m_order;
SmallArray<const TinyVector<Dimension>> m_point_list;
SmallArray<const double> m_weight_list;
virtual void _buildPointAndWeightLists() = 0;
public:
// PUGS_INLINE virtual T integrate(const FunctionSymbolId& function, const double& a, const double& b) = 0;
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);
PUGS_INLINE virtual Array<TinyVector<1>> quadraturePoints(const double& a, const double& b) = 0;
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:
PUGS_INLINE virtual SmallArray<TinyVector<1>> quadraturePoints(const double& a, const double& b) = 0;
PUGS_INLINE virtual size_t numberOfPoints() = 0;
PUGS_INLINE virtual Array<double> weights() = 0;
PUGS_INLINE virtual SmallArray<double> weights() = 0;
// One does not use the '=default' constructor to avoid
// (zero-initialization) performances issues
PUGS_INLINE
......@@ -45,16 +390,16 @@ class IntegrationMethod
};
template <size_t Order>
class IntegrationMethodLobatto : public IntegrationMethod<Order>
class IntegrationMethodLobatto : public IntegrationMethod
{
private:
static constexpr size_t number_points = (Order < 4 ? 3 : (Order < 6 ? 4 : (Order < 8 ? 5 : Order < 10 ? 6 : 7)));
Array<double> m_weights;
Array<TinyVector<1>> m_quadrature_positions;
PUGS_INLINE Array<TinyVector<1>>
SmallArray<double> m_weights;
SmallArray<TinyVector<1>> m_quadrature_positions;
PUGS_INLINE SmallArray<TinyVector<1>>
translateQuadraturePoints(const double& a, const double& b)
{
Array<TinyVector<1>> points{number_points};
SmallArray<TinyVector<1>> points{number_points};
double length = b - a;
for (size_t i = 0; i < m_quadrature_positions.size(); i++) {
TinyVector<1> qpos{m_quadrature_positions[i]};
......@@ -64,7 +409,7 @@ class IntegrationMethodLobatto : public IntegrationMethod<Order>
}
PUGS_INLINE
constexpr void
void
fillArrayLobatto()
{
switch (Order) {
......@@ -180,7 +525,7 @@ class IntegrationMethodLobatto : public IntegrationMethod<Order>
}
PUGS_INLINE
Array<TinyVector<1>>
SmallArray<TinyVector<1>>
quadraturePoints(const double& a, const double& b)
{
return translateQuadraturePoints(a, b);
......@@ -192,7 +537,7 @@ class IntegrationMethodLobatto : public IntegrationMethod<Order>
return number_points;
}
PUGS_INLINE
Array<double>
SmallArray<double>
weights()
{
return m_weights;
......@@ -209,16 +554,16 @@ class IntegrationMethodLobatto : public IntegrationMethod<Order>
};
template <size_t Order>
class IntegrationMethodLegendre : public IntegrationMethod<Order>
class IntegrationMethodLegendre : public IntegrationMethod
{
private:
static constexpr size_t number_points = (Order < 4 ? 2 : (Order < 6 ? 3 : (Order < 8 ? 4 : 5)));
Array<double> m_weights;
Array<TinyVector<1>> m_quadrature_positions;
PUGS_INLINE Array<TinyVector<1>>
SmallArray<double> m_weights;
SmallArray<TinyVector<1>> m_quadrature_positions;
PUGS_INLINE SmallArray<TinyVector<1>>
translateQuadraturePoints(const double& a, const double& b)
{
Array<TinyVector<1>> points{number_points};
SmallArray<TinyVector<1>> points{number_points};
double length = b - a;
for (size_t i = 0; i < m_quadrature_positions.size(); i++) {
TinyVector<1> qpos{m_quadrature_positions[i]};
......@@ -228,78 +573,81 @@ class IntegrationMethodLegendre : public IntegrationMethod<Order>
}
PUGS_INLINE
constexpr void
void
fillArrayLegendre()
{
switch (Order) {
case 0:
case 1:
case 2:
case 3: {
double oneovsqrt3 = 1. / std::sqrt(3);
m_quadrature_positions[0] = -oneovsqrt3;
m_quadrature_positions[1] = oneovsqrt3;
m_weights[0] = 1;
m_weights[1] = 1;
// m_quadrature_positions.fill({-1, 0, 1});
// m_weights.fill({oneov3, 4 * oneov3, oneov3});
break;
}
case 4:
case 5: {
double oneover9 = 1. / 9;
double coef = std::sqrt(3. / 5.);
m_quadrature_positions[0] = -coef;
m_quadrature_positions[1] = 0;
m_quadrature_positions[2] = coef;
m_weights[0] = 5 * oneover9;
m_weights[1] = 8 * oneover9;
m_weights[2] = 5 * oneover9;
break;
}
case 6:
case 7: {
double oneov36 = 1. / 36.;
double coef1 = std::sqrt(3. / 7. - 2. / 7. * std::sqrt(6. / 5.));
double coef2 = std::sqrt(3. / 7. + 2. / 7. * std::sqrt(6. / 5.));
double coef3 = 18. + std::sqrt(30.);
double coef4 = 18. - std::sqrt(30.);
m_quadrature_positions[0] = -coef2;
m_quadrature_positions[1] = -coef1;
m_quadrature_positions[2] = coef1;
m_quadrature_positions[3] = coef2;
m_weights[0] = coef4 * oneov36;
m_weights[1] = coef3 * oneov36;
m_weights[2] = coef3 * oneov36;
m_weights[3] = coef4 * oneov36;
break;
}
case 8:
case 9: {
double oneov900 = 1. / 900.;
double oneov225 = 1. / 225.;
double coef1 = 1. / 3. * std::sqrt(5. - 2. * std::sqrt(10. / 7.));
double coef2 = 1. / 3. * std::sqrt(5. + 2. * std::sqrt(10. / 7.));
double coef3 = 322. + 13. * std::sqrt(70.);
double coef4 = 322. - 13. * std::sqrt(70.);
m_quadrature_positions[0] = -coef2;
m_quadrature_positions[1] = -coef1;
m_quadrature_positions[2] = 0;
m_quadrature_positions[3] = coef1;
m_quadrature_positions[4] = coef2;
m_weights[0] = coef4 * oneov900;
m_weights[1] = coef3 * oneov900;
m_weights[2] = 128. * oneov225;
m_weights[3] = coef3 * oneov900;
m_weights[4] = coef4 * oneov900;
break;
}
default: {
throw UnexpectedError("Gauss-Legendre quadratures handle orders up to 9.");
break;
}
}
TensorialGaussLegendre<1> tensorial_gauss_legendre(Order);
copy_to(tensorial_gauss_legendre.pointList(), m_quadrature_positions);
copy_to(tensorial_gauss_legendre.weightList(), m_weights);
// switch (Order) {
// case 0:
// case 1:
// case 2:
// case 3: {
// double oneovsqrt3 = 1. / std::sqrt(3);
// m_quadrature_positions[0] = -oneovsqrt3;
// m_quadrature_positions[1] = oneovsqrt3;
// m_weights[0] = 1;
// m_weights[1] = 1;
// // m_quadrature_positions.fill({-1, 0, 1});
// // m_weights.fill({oneov3, 4 * oneov3, oneov3});
// break;
// }
// case 4:
// case 5: {
// double oneover9 = 1. / 9;
// double coef = std::sqrt(3. / 5.);
// m_quadrature_positions[0] = -coef;
// m_quadrature_positions[1] = 0;
// m_quadrature_positions[2] = coef;
// m_weights[0] = 5 * oneover9;
// m_weights[1] = 8 * oneover9;
// m_weights[2] = 5 * oneover9;
// break;
// }
// case 6:
// case 7: {
// double oneov36 = 1. / 36.;
// double coef1 = std::sqrt(3. / 7. - 2. / 7. * std::sqrt(6. / 5.));
// double coef2 = std::sqrt(3. / 7. + 2. / 7. * std::sqrt(6. / 5.));
// double coef3 = 18. + std::sqrt(30.);
// double coef4 = 18. - std::sqrt(30.);
// m_quadrature_positions[0] = -coef2;
// m_quadrature_positions[1] = -coef1;
// m_quadrature_positions[2] = coef1;
// m_quadrature_positions[3] = coef2;
// m_weights[0] = coef4 * oneov36;
// m_weights[1] = coef3 * oneov36;
// m_weights[2] = coef3 * oneov36;
// m_weights[3] = coef4 * oneov36;
// break;
// }
// case 8:
// case 9: {
// double oneov900 = 1. / 900.;
// double oneov225 = 1. / 225.;
// double coef1 = 1. / 3. * std::sqrt(5. - 2. * std::sqrt(10. / 7.));
// double coef2 = 1. / 3. * std::sqrt(5. + 2. * std::sqrt(10. / 7.));
// double coef3 = 322. + 13. * std::sqrt(70.);
// double coef4 = 322. - 13. * std::sqrt(70.);
// m_quadrature_positions[0] = -coef2;
// m_quadrature_positions[1] = -coef1;
// m_quadrature_positions[2] = 0;
// m_quadrature_positions[3] = coef1;
// m_quadrature_positions[4] = coef2;
// m_weights[0] = coef4 * oneov900;
// m_weights[1] = coef3 * oneov900;
// m_weights[2] = 128. * oneov225;
// m_weights[3] = coef3 * oneov900;
// m_weights[4] = coef4 * oneov900;
// break;
// }
// default: {
// throw UnexpectedError("Gauss-Legendre quadratures handle orders up to 9.");
// break;
// }
// }
}
public:
......@@ -314,7 +662,7 @@ class IntegrationMethodLegendre : public IntegrationMethod<Order>
}
PUGS_INLINE
Array<TinyVector<1>>
SmallArray<TinyVector<1>>
quadraturePoints(const double& a, const double& b)
{
return translateQuadraturePoints(a, b);
......@@ -326,7 +674,7 @@ class IntegrationMethodLegendre : public IntegrationMethod<Order>
return number_points;
}
PUGS_INLINE
Array<double>
SmallArray<double>
weights()
{
return m_weights;
......@@ -350,7 +698,7 @@ class IntegrationTools
private:
// T m_values[N];
std::shared_ptr<IntegrationMethod<Order>> m_integration_method;
std::shared_ptr<IntegrationMethod> m_integration_method;
size_t N;
public:
......@@ -358,9 +706,10 @@ class IntegrationTools
integrate(const FunctionSymbolId& function, const double& a, const double& b) const
{
T result = 0;
Array<const TinyVector<1>> positions = quadraturePoints(a, b);
Array<double> weights = m_integration_method->weights();
Array<T> values = EvaluateAtPoints<double(const TinyVector<1>)>::evaluate(function, positions);
SmallArray<const TinyVector<1>> positions = quadraturePoints(a, b);
SmallArray<double> weights = m_integration_method->weights();
SmallArray<T> values(positions.size());
EvaluateAtPoints<double(const TinyVector<1>)>::evaluateTo(function, positions, values);
Assert(positions.size() == weights.size(), "Wrong number of quadrature points and/or weights in Gauss quadrature");
for (size_t i = 0; i < values.size(); ++i) {
result += weights[i] * values[i];
......@@ -368,14 +717,15 @@ class IntegrationTools
return (b - a) / 2 * result;
}
PUGS_INLINE Array<T>
integrateFunction(const FunctionSymbolId& function, const Array<TinyVector<1>>& vertices) const
PUGS_INLINE SmallArray<T>
integrateFunction(const FunctionSymbolId& function, const SmallArray<TinyVector<1>>& vertices) const
{
Array<const TinyVector<1>> positions = quadraturePoints(vertices);
Array<T> result{vertices.size() - 1};
Array<double> interval_size{vertices.size() - 1};
Array<double> weights = m_integration_method->weights();
Array<T> values = EvaluateAtPoints<double(const TinyVector<1>)>::evaluate(function, positions);
SmallArray<const TinyVector<1>> positions = quadraturePoints(vertices);
SmallArray<T> result{vertices.size() - 1};
SmallArray<double> interval_size{vertices.size() - 1};
SmallArray<double> weights = m_integration_method->weights();
SmallArray<T> values(positions.size());
EvaluateAtPoints<double(const TinyVector<1>)>::evaluateTo(function, positions, values);
for (size_t i = 0; i < interval_size.size(); ++i) {
interval_size[i] = vertices[i + 1][0] - vertices[i][0];
......@@ -392,26 +742,30 @@ class IntegrationTools
testIntegrate(const double& a, const double& b) const
{
T result = 0;
Array<TinyVector<1>> positions = quadraturePoints(a, b);
Array<double> weights = m_integration_method->weights();
Array<T> values(weights.size());
SmallArray<TinyVector<1>> positions = quadraturePoints(a, b);
SmallArray<double> weights = m_integration_method->weights();
SmallArray<T> values(weights.size());
for (size_t j = 0; j < weights.size(); ++j) {
values[j] = std::pow(positions[j][0], 3);
values[j] = std::pow(positions[j][0], 3.);
}
Assert(positions.size() == weights.size(), "Wrong number of quadrature points and/or weights in Gauss quadrature");
for (size_t i = 0; i < values.size(); ++i) {
result += weights[i] * values[i];
}
return (b - a) / 2 * result;
return 0.5 * (b - a) * result;
}
PUGS_INLINE Array<T>
testIntegrateFunction(const Array<TinyVector<1>>& vertices) const
PUGS_INLINE SmallArray<T>
testIntegrateFunction(const SmallArray<TinyVector<1>>& vertices) const
{
Array<TinyVector<1>> positions = quadraturePoints(vertices);
Array<double> interval_size{vertices.size() - 1};
Array<double> weights = m_integration_method->weights();
Array<T> values(positions.size());
SmallArray<TinyVector<1>> positions = quadraturePoints(vertices);
SmallArray<double> interval_size{vertices.size() - 1};
SmallArray<double> weights = m_integration_method->weights();
SmallArray<T> values(positions.size());
for (size_t j = 0; j < positions.size(); ++j) {
values[j] = std::pow(positions[j][0], 3);
}
......@@ -419,7 +773,7 @@ class IntegrationTools
for (size_t i = 0; i < interval_size.size(); ++i) {
interval_size[i] = vertices[i + 1][0] - vertices[i][0];
}
Array<T> result{vertices.size() - 1};
SmallArray<T> result{vertices.size() - 1};
if constexpr (std::is_arithmetic_v<T>) {
result.fill(0);
} else if constexpr (is_tiny_vector_v<T> or is_tiny_matrix_v<T>) {
......@@ -435,19 +789,19 @@ class IntegrationTools
return result;
}
PUGS_INLINE Array<TinyVector<1>>
PUGS_INLINE SmallArray<TinyVector<1>>
quadraturePoints(const double& a, const double& b) const
{
return m_integration_method->quadraturePoints(a, b);
}
PUGS_INLINE
Array<TinyVector<1>>
quadraturePoints(const Array<TinyVector<1>>& positions) const
SmallArray<TinyVector<1>>
quadraturePoints(const SmallArray<TinyVector<1>>& positions) const
{
size_t number_of_intervals = positions.size() - 1;
size_t quadrature_size = m_integration_method->numberOfPoints();
Array<TinyVector<1>> quadratures{quadrature_size * number_of_intervals};
SmallArray<TinyVector<1>> quadratures{quadrature_size * number_of_intervals};
for (size_t j = 0; j < number_of_intervals; ++j) {
double a = positions[j][0];
double b = positions[j + 1][0];
......@@ -459,7 +813,7 @@ class IntegrationTools
return quadratures;
}
PUGS_INLINE Array<double>
PUGS_INLINE SmallArray<double>
weights() const
{
return m_integration_method->weights();
......
tests/test_FunctionEvaluation.cpp
+ 73
− 14
View file @ 546807e7
......@@ -24,19 +24,76 @@ TEST_CASE("EvaluateAtPoints", "[integration]")
SECTION("Instantiation")
{
// positions[0] = 1;
Array<double> result(1);
SmallArray<double> result(1);
// result = EvaluateAtPoints<double(TinyVector<1>)>::evaluate(f, positions);
// REQUIRE(result[0] == 1);
}
}
TEST_CASE("QuadratureFormula", "[analysis]")
{
SECTION("TensorialGaussLegendre")
{
auto integrate = [](auto f, auto quadrature_formula, const double a, const double b) {
auto point_list = quadrature_formula.pointList();
auto weight_list = quadrature_formula.weightList();
auto x = [&a, &b](auto x_hat) { return 0.5 * (b - a) * (x_hat + 1); };
auto value = weight_list[0] * f(x(point_list[0]));
for (size_t i = 1; i < weight_list.size(); ++i) {
value += weight_list[i] * f(x(point_list[i]));
}
return 0.5 * (b - a) * value;
};
auto p0 = [](const TinyVector<1>&) { return 3; };
auto p1 = [](const TinyVector<1>& X) {
const double x = X[0];
return 2 * x + 1;
};
auto p2 = [](const TinyVector<1>& X) {
const double x = X[0];
return 3 * x * x + 2 * x + 1;
};
SECTION("exact integration")
{
SECTION("order 0 and 1")
{
TensorialGaussLegendre<1> l0(0);
TensorialGaussLegendre<1> l1(1);
REQUIRE(l0.numberOfPoints() == 1);
REQUIRE(l1.numberOfPoints() == 1);
REQUIRE(integrate(p0, l1, 0.5, 1) == Catch::Approx(1.5));
REQUIRE(integrate(p1, l1, 0, 1) == Catch::Approx(2));
}
SECTION("order 2 and 3")
{
TensorialGaussLegendre<1> l2(2);
TensorialGaussLegendre<1> l3(3);
REQUIRE(l2.numberOfPoints() == 2);
REQUIRE(l3.numberOfPoints() == 2);
REQUIRE(integrate(p0, l3, 0.5, 1) == Catch::Approx(1.5));
REQUIRE(integrate(p1, l3, 0, 1) == Catch::Approx(2));
REQUIRE(integrate(p2, l3, 0, 1) == Catch::Approx(3));
}
}
}
}
TEST_CASE("IntegrationTools", "[integration]")
{
IntegrationTools<2> integrationtools(QuadratureType::gausslobatto);
double a = 1;
double b = 2;
Array<TinyVector<1>> vertices(4);
SmallArray<TinyVector<1>> vertices(4);
for (size_t i = 0; i < 4; i++) {
vertices[0] = {a};
vertices[1] = {b};
......@@ -50,19 +107,21 @@ TEST_CASE("IntegrationTools", "[integration]")
}
SECTION("QuadraturePoints")
{
Array<TinyVector<1>> quadratures = integrationtools.quadraturePoints(a, b);
REQUIRE(((quadratures[0] == 1) and (quadratures[1] == 1.5) and (quadratures[2] == 2)));
SmallArray<TinyVector<1>> quadratures = integrationtools.quadraturePoints(a, b);
REQUIRE(quadratures[0] == 1);
REQUIRE(quadratures[1] == 1.5);
REQUIRE(quadratures[2] == 2);
}
SECTION("Weights")
{
Array<double> weights = integrationtools.weights();
SmallArray<double> weights = integrationtools.weights();
REQUIRE(((weights[0] == 1. / 3.) and (weights[1] == 4. / 3.) and (weights[2] == 1. / 3.)));
}
SECTION("TestOrder3")
{
double result = integrationtools.testIntegrate(a, b);
REQUIRE(result == Catch::Approx((0.25 * (std::pow(2, 4) - std::pow(1, 4)))).epsilon(1E-14));
Array<double> results = integrationtools.testIntegrateFunction(vertices);
SmallArray<double> results = integrationtools.testIntegrateFunction(vertices);
REQUIRE(results.size() == 3);
REQUIRE(results[0] == Catch::Approx((0.25 * (std::pow(2, 4) - std::pow(1, 4)))).epsilon(1E-14));
REQUIRE(results[1] == 0);
......@@ -73,7 +132,7 @@ TEST_CASE("IntegrationTools", "[integration]")
{
double result = integrationtools5.testIntegrate(a, b);
REQUIRE(result == Catch::Approx((0.25 * (std::pow(2, 4) - std::pow(1, 4)))).epsilon(1E-14));
Array<double> results = integrationtools.testIntegrateFunction(vertices);
SmallArray<double> results = integrationtools.testIntegrateFunction(vertices);
REQUIRE(results[0] == Catch::Approx((0.25 * (std::pow(2, 4) - std::pow(1, 4)))).epsilon(1E-14));
REQUIRE(results[1] == 0);
REQUIRE(results[2] == Catch::Approx((-0.25 * (std::pow(2, 4) - std::pow(1, 4)))).epsilon(1E-14));
......@@ -83,7 +142,7 @@ TEST_CASE("IntegrationTools", "[integration]")
{
double result = integrationtools7.testIntegrate(a, b);
REQUIRE(result == Catch::Approx((0.25 * (std::pow(2, 4) - std::pow(1, 4)))).epsilon(1E-14));
Array<double> results = integrationtools.testIntegrateFunction(vertices);
SmallArray<double> results = integrationtools.testIntegrateFunction(vertices);
REQUIRE(results[0] == Catch::Approx((0.25 * (std::pow(2, 4) - std::pow(1, 4)))).epsilon(1E-14));
REQUIRE(results[1] == 0);
REQUIRE(results[2] == Catch::Approx((-0.25 * (std::pow(2, 4) - std::pow(1, 4)))).epsilon(1E-14));
......@@ -93,7 +152,7 @@ TEST_CASE("IntegrationTools", "[integration]")
{
double result = integrationtools9.testIntegrate(a, b);
REQUIRE(result == Catch::Approx((0.25 * (std::pow(2, 4) - std::pow(1, 4)))).epsilon(1E-14));
Array<double> results = integrationtools.testIntegrateFunction(vertices);
SmallArray<double> results = integrationtools.testIntegrateFunction(vertices);
REQUIRE(results[0] == Catch::Approx((0.25 * (std::pow(2, 4) - std::pow(1, 4)))).epsilon(1E-14));
REQUIRE(results[1] == 0);
REQUIRE(results[2] == Catch::Approx((-0.25 * (std::pow(2, 4) - std::pow(1, 4)))).epsilon(1E-14));
......@@ -103,7 +162,7 @@ TEST_CASE("IntegrationTools", "[integration]")
{
double result = integrationtools11.testIntegrate(a, b);
REQUIRE(result == Catch::Approx((0.25 * (std::pow(2, 4) - std::pow(1, 4)))).epsilon(1E-14));
Array<double> results = integrationtools.testIntegrateFunction(vertices);
SmallArray<double> results = integrationtools.testIntegrateFunction(vertices);
REQUIRE(results[0] == Catch::Approx((0.25 * (std::pow(2, 4) - std::pow(1, 4)))).epsilon(1E-14));
REQUIRE(results[1] == 0);
REQUIRE(results[2] == Catch::Approx((-0.25 * (std::pow(2, 4) - std::pow(1, 4)))).epsilon(1E-14));
......@@ -119,7 +178,7 @@ TEST_CASE("IntegrationTools", "[integration]")
{
double result = integrationtoolslegendre.testIntegrate(a, b);
REQUIRE(result == Catch::Approx((0.25 * (std::pow(2, 4) - std::pow(1, 4)))).epsilon(1E-14));
Array<double> results = integrationtoolslegendre.testIntegrateFunction(vertices);
SmallArray<double> results = integrationtoolslegendre.testIntegrateFunction(vertices);
REQUIRE(results.size() == 3);
REQUIRE(results[0] == Catch::Approx((0.25 * (std::pow(2, 4) - std::pow(1, 4)))).epsilon(1E-14));
REQUIRE(results[1] == 0);
......@@ -130,7 +189,7 @@ TEST_CASE("IntegrationTools", "[integration]")
{
double result = integrationtoolslegendre5.testIntegrate(a, b);
REQUIRE(result == Catch::Approx((0.25 * (std::pow(2, 4) - std::pow(1, 4)))).epsilon(1E-14));
Array<double> results = integrationtoolslegendre5.testIntegrateFunction(vertices);
SmallArray<double> results = integrationtoolslegendre5.testIntegrateFunction(vertices);
REQUIRE(results[0] == Catch::Approx((0.25 * (std::pow(2, 4) - std::pow(1, 4)))).epsilon(1E-14));
REQUIRE(results[1] == 0);
REQUIRE(results[2] == Catch::Approx((-0.25 * (std::pow(2, 4) - std::pow(1, 4)))).epsilon(1E-14));
......@@ -140,7 +199,7 @@ TEST_CASE("IntegrationTools", "[integration]")
{
double result = integrationtoolslegendre7.testIntegrate(a, b);
REQUIRE(result == Catch::Approx((0.25 * (std::pow(2, 4) - std::pow(1, 4)))).epsilon(1E-14));
Array<double> results = integrationtoolslegendre7.testIntegrateFunction(vertices);
SmallArray<double> results = integrationtoolslegendre7.testIntegrateFunction(vertices);
REQUIRE(results[0] == Catch::Approx((0.25 * (std::pow(2, 4) - std::pow(1, 4)))).epsilon(1E-14));
REQUIRE(results[1] == 0);
REQUIRE(results[2] == Catch::Approx((-0.25 * (std::pow(2, 4) - std::pow(1, 4)))).epsilon(1E-14));
......@@ -150,7 +209,7 @@ TEST_CASE("IntegrationTools", "[integration]")
{
double result = integrationtoolslegendre9.testIntegrate(a, b);
REQUIRE(result == Catch::Approx((0.25 * (std::pow(2, 4) - std::pow(1, 4)))).epsilon(1E-14));
Array<double> results = integrationtoolslegendre9.testIntegrateFunction(vertices);
SmallArray<double> results = integrationtoolslegendre9.testIntegrateFunction(vertices);
REQUIRE(results[0] == Catch::Approx((0.25 * (std::pow(2, 4) - std::pow(1, 4)))).epsilon(1E-14));
REQUIRE(results[1] == 0);
REQUIRE(results[2] == Catch::Approx((-0.25 * (std::pow(2, 4) - std::pow(1, 4)))).epsilon(1E-14));
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please to comment