diff --git a/src/algebra/IntegrationTools.hpp b/src/algebra/IntegrationTools.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..1ffc6af1740426147e226404e35062a7ec749c7d
--- /dev/null
+++ b/src/algebra/IntegrationTools.hpp
@@ -0,0 +1,252 @@
+#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
diff --git a/src/language/utils/EvaluateAtPoints.hpp b/src/language/utils/EvaluateAtPoints.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..81d055dfbe113c50a16970b5b64df904e2810fbe
--- /dev/null
+++ b/src/language/utils/EvaluateAtPoints.hpp
@@ -0,0 +1,46 @@
+#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