Skip to content
Snippets Groups Projects
Commit 6bc49ab3 authored by Emmanuel Labourasse's avatar Emmanuel Labourasse
Browse files

Add files for high order integration with quadratures

warning: Untested yet
parent b5d624f9
Branches
Tags
1 merge request!124Add files for high order integration with quadratures
#ifndef INTEGRATION_TOOLS_HPP
#define INTEGRATION_TOOLS_HPP
#include <iostream>
#include <utils/Exceptions.hpp>
#include <utils/PugsAssert.hpp>
#include <utils/PugsMacros.hpp>
#include <cmath>
#include <language/utils/EvaluateAtPoints.hpp>
#include <language/utils/PugsFunctionAdapter.hpp>
#include <utils/Types.hpp>
enum class QuadratureType : int8_t
{
QT__begin = 0,
//
gausslobatto = QT__begin,
//
QT__end
};
template <size_t Order, typename T = double>
class IntegrationMethod
{
public:
PUGS_INLINE virtual T integrate(const FunctionSymbolId& function, const double& a, const double& b) = 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, typename T = double>
class IntegrationMethodGL : public IntegrationMethod<Order, T>
{
private:
Array<double> weights;
Array<double> quadrature_positions;
PUGS_INLINE Array<TinyVector<1>>
translateQuadraturePoints(const double& a, const double& b)
{
Array<TinyVector<1>> points;
double length = b - a;
for (size_t i = 0; i < quadrature_positions.size(); i++) {
TinyVector<1> qpos{quadrature_positions[i]};
points[i] = 0.5 * (length * qpos + a + b);
}
return points;
}
PUGS_INLINE
constexpr void
fillArrayGL() noexcept
{
switch (Order) {
case 0:
case 1:
case 2:
case 3: {
double oneov3 = 1. / 3.;
quadrature_positions[0] = -1.;
quadrature_positions[1] = 0.;
quadrature_positions[2] = 1.;
weights[0] = oneov3;
weights[1] = 4 * oneov3;
weights[2] = oneov3;
// quadrature_positions.fill({-1, 0, 1});
// weights.fill({oneov3, 4 * oneov3, oneov3});
break;
}
case 4:
case 5: {
double oneov6 = 1. / 6.;
double coef = std::sqrt(1. / 5.);
quadrature_positions[0] = -1;
quadrature_positions[1] = -coef;
quadrature_positions[2] = coef;
quadrature_positions[3] = 1;
weights[0] = oneov6;
weights[1] = 5 * oneov6;
weights[2] = 5 * oneov6;
weights[3] = oneov6;
break;
}
case 6:
case 7: {
double oneov45 = 1. / 45.;
double oneov10 = 1. / 10.;
double coef = std::sqrt(3. / 7.);
quadrature_positions[0] = -1;
quadrature_positions[1] = -coef;
quadrature_positions[2] = 0;
quadrature_positions[3] = coef;
quadrature_positions[4] = 1;
weights[0] = oneov10;
weights[1] = 98 * oneov45;
weights[2] = 32 * oneov45;
weights[3] = 98 * oneov45;
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 = std::sqrt(14. + std::sqrt(7)) * oneov30;
double coef4 = std::sqrt(14. - std::sqrt(7)) * oneov30;
quadrature_positions[0] = -1;
quadrature_positions[1] = -coef2;
quadrature_positions[2] = -coef1;
quadrature_positions[3] = coef1;
quadrature_positions[4] = coef2;
quadrature_positions[5] = 1;
weights[0] = 2. * oneov30;
weights[1] = coef4 * oneov30;
weights[2] = coef3 * oneov30;
weights[3] = coef3 * oneov30;
weights[4] = coef4 * oneov30;
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. / 11. - 2 * std::sqrt(5. / 3.) / 11.);
double coef2 = std::sqrt(5. / 11. + 2 * std::sqrt(5. / 3.) / 11.);
double coef3 = 124. + 7 * std::sqrt(15.) * oneov350;
double coef4 = 124. - 7 * std::sqrt(15.) * oneov350;
quadrature_positions[0] = -1;
quadrature_positions[1] = -coef2;
quadrature_positions[2] = -coef1;
quadrature_positions[3] = 0;
quadrature_positions[4] = coef1;
quadrature_positions[5] = coef2;
quadrature_positions[6] = 1;
weights[0] = oneov21;
weights[1] = coef4 * oneov350;
weights[2] = coef3 * oneov350;
weights[3] = coef3 * oneov350;
weights[4] = coef0;
weights[5] = coef4 * oneov350;
weights[6] = oneov21;
break;
}
default: {
throw UnexpectedError("Gauss-Lobatto quadratures handle orders up to 11.");
break;
}
}
}
public:
PUGS_INLINE T
integrate(const FunctionSymbolId& function, const double& a, const double& b)
{
T result = 0;
Array<TinyVector<1>> positions = translateQuadraturePoints(a, b);
Array<T> values = EvaluateAtPoints<double(TinyVector<1>)>::evaluate(function, positions);
Assert(positions.size() == weights.size(),
"Wrong number of quadrature points and/or weights in Gauss-Lobatto quadrature");
for (size_t i = 0; i < values.size(); ++i) {
result += weights[i] * values[i];
}
return result;
};
// One does not use the '=default' constructor to avoid
// (zero-initialization) performances issues
PUGS_INLINE
constexpr IntegrationMethodGL() noexcept
{
fillArrayGL();
}
PUGS_INLINE
constexpr IntegrationMethodGL(const IntegrationMethodGL&) noexcept = default;
PUGS_INLINE
constexpr IntegrationMethodGL(IntegrationMethodGL&& v) noexcept = default;
PUGS_INLINE
~IntegrationMethodGL() 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<Order, T>> m_integration_method;
size_t N;
public:
PUGS_INLINE T
integrate(const FunctionSymbolId& function, const double& a, const double& b) const
{
return m_integration_method->integrate(function, a, b);
}
PUGS_INLINE
constexpr size_t
order() const
{
return Order;
}
PUGS_INLINE constexpr IntegrationTools(QuadratureType quadrature)
{
switch (quadrature) {
case QuadratureType::gausslobatto: {
m_integration_method = std::make_shared<IntegrationMethodGL<Order, T>>();
break;
}
case QuadratureType::QT__end: {
throw UnexpectedError("QuadratureType is not defined!");
}
}
}
// One does not use the '=default' constructor to avoid
// (zero-initialization) performances issues
PUGS_INLINE
constexpr IntegrationTools() noexcept {}
PUGS_INLINE
constexpr IntegrationTools(const IntegrationTools&) noexcept = default;
PUGS_INLINE
constexpr IntegrationTools(IntegrationTools&& v) noexcept = default;
PUGS_INLINE
~IntegrationTools() noexcept = default;
};
#endif // INTEGRATION_TOOLS_HPP
#ifndef EVALUATE_AT_POINTS_HPP
#define EVALUATE_AT_POINTS_HPP
#include <language/utils/PugsFunctionAdapter.hpp>
template <typename T>
class EvaluateAtPoints;
template <typename OutputType, typename InputType>
class EvaluateAtPoints<OutputType(InputType)> : public PugsFunctionAdapter<OutputType(InputType)>
{
// static constexpr size_t Dimension = OutputType::Dimension;
using Adapter = PugsFunctionAdapter<OutputType(InputType)>;
public:
// template <size_t vectorlength>
static inline Array<OutputType>
evaluate(const FunctionSymbolId& function_symbol_id, const Array<const InputType>& position)
{
// static_assert(position.size()==vectorlength, "The length of the output must be equal to the template
// argument");
auto& expression = Adapter::getFunctionExpression(function_symbol_id);
auto convert_result = Adapter::getResultConverter(expression.m_data_type);
Array<ExecutionPolicy> context_list = Adapter::getContextList(expression);
using execution_space = typename Kokkos::DefaultExecutionSpace::execution_space;
Kokkos::Experimental::UniqueToken<execution_space, Kokkos::Experimental::UniqueTokenScope::Global> tokens;
Array<OutputType> value(position.size());
parallel_for(position.size(), [=, &expression, &tokens](size_t i) {
const int32_t t = tokens.acquire();
auto& execution_policy = context_list[t];
Adapter::convertArgs(execution_policy.currentContext(), position[i]);
auto result = expression.execute(execution_policy);
value[i] = convert_result(std::move(result));
tokens.release(t);
});
return value;
}
};
#endif // INTERPOLATE_ITEM_VALUE_HPP
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment