diff --git a/src/algebra/Givens.hpp b/src/algebra/Givens.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..02a71d0aa97d646689d2ef0b6f059131478f38f9
--- /dev/null
+++ b/src/algebra/Givens.hpp
@@ -0,0 +1,99 @@
+#ifndef GIVENS_HPP
+#define GIVENS_HPP
+#include <cmath>
+#include <iomanip>
+#include <iostream>
+#include <rang.hpp>
+#include <utils/Exceptions.hpp>
+#include <utils/PugsAssert.hpp>
+
+class Givens
+{
+ private:
+ static void
+ _givens(const double a, const double b, double& c, double& s)
+ {
+ if (b == 0) {
+ c = 1;
+ s = 0;
+ } else {
+ if (std::abs(b) > std::abs(a)) {
+ const double tau = -a / b;
+ s = 1. / sqrt(1 + tau * tau);
+ c = s * tau;
+ } else {
+ Assert(a != 0);
+ const double tau = -b / a;
+ c = 1 / sqrt(1 + tau * tau);
+ s = c * tau;
+ }
+ }
+ }
+
+ static void
+ _rotate(const double Ai_1j, const double Aij, const double c, const double s, double& ui_1, double& ui)
+ {
+ ui_1 = c * Ai_1j - s * Aij;
+ ui = s * Ai_1j + c * Aij;
+ }
+
+ template <typename MatrixType, typename VectorType, typename RHSVectorType>
+ static void
+ _lift(const MatrixType& A, const RHSVectorType& b, VectorType& x)
+ {
+ x.fill(0);
+ for (size_t i = x.size(); i >= 1; --i) {
+ double sum = 0;
+ for (size_t j = i; j < x.size(); ++j) {
+ sum += A(i - 1, j) * x[j];
+ }
+ x[i - 1] = (b[i - 1] - sum) / A(i - 1, i - 1);
+ }
+ }
+
+ public:
+ template <typename MatrixType, typename VectorType, typename RHSVectorType>
+ static void
+ solveInPlace(MatrixType& A, VectorType& x, RHSVectorType& b)
+ {
+ Assert(x.size() == A.numberOfColumns(), "The number of columns of A must be the size of x");
+ Assert(b.size() == A.numberOfRows(), "The number of rows of A must be the size of b");
+
+ for (size_t j = 0; j < A.numberOfColumns(); ++j) {
+ for (size_t i = A.numberOfRows() - 1; i > j; --i) {
+ double c(0), s(0);
+ _givens(A(i - 1, j), A(i, j), c, s);
+ for (size_t k = j; k < A.numberOfColumns(); ++k) {
+ double xi_1(0), xi(0);
+ _rotate(A(i - 1, k), A(i, k), c, s, xi_1, xi);
+ A(i - 1, k) = xi_1;
+ A(i, k) = xi;
+ }
+ double xi_1(0), xi(0);
+ _rotate(b[i - 1], b[i], c, s, xi_1, xi);
+ b[i - 1] = xi_1;
+ b[i] = xi;
+ }
+ }
+
+ _lift(A, b, x);
+ }
+
+ template <typename MatrixType, typename VectorType, typename RHSVectorType>
+ static void
+ solve(const MatrixType& A, VectorType& x, const RHSVectorType& b)
+ {
+ Assert(x.size() == A.numberOfColumns(), "The number of columns of A must be the size of x");
+ Assert(b.size() == A.numberOfRows(), "The number of rows of A must be the size of b");
+
+ MatrixType rotateA = copy(A);
+ RHSVectorType rotateb = copy(b);
+
+ solveInPlace(rotateA, x, rotateb);
+ }
+
+ Givens() = delete;
+ ~Givens() = delete;
+};
+
+#endif // GIVENS_HPP
diff --git a/src/scheme/PolynomialReconstruction.hpp b/src/scheme/PolynomialReconstruction.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..7010d3d196b012b2f39a434c0c81111844f4b608
--- /dev/null
+++ b/src/scheme/PolynomialReconstruction.hpp
@@ -0,0 +1,81 @@
+#ifndef POLYNOMIAL_RECONSTRUCTION_HPP
+#define POLYNOMIAL_RECONSTRUCTION_HPP
+
+#include <mesh/MeshTraits.hpp>
+#include <mesh/StencilBuilder.hpp>
+#include <scheme/DiscreteFunctionDPk.hpp>
+
+#include <algebra/SmallMatrix.hpp>
+#include <algebra/SmallVector.hpp>
+
+#include <algebra/Givens.hpp>
+#include <analysis/GaussLegendreQuadratureDescriptor.hpp>
+#include <analysis/QuadratureFormula.hpp>
+#include <analysis/QuadratureManager.hpp>
+#include <geometry/LineTransformation.hpp>
+#include <mesh/MeshData.hpp>
+#include <mesh/MeshDataManager.hpp>
+
+class PolynomialReconstruction
+{
+ public:
+ template <MeshConcept MeshType, typename DataType>
+ PUGS_INLINE DiscreteFunctionDPk<MeshType::Dimension, std::remove_const_t<DataType>>
+ build(const MeshType& mesh, const DiscreteFunctionP0<DataType> p0_function)
+ {
+ const size_t degree = 1;
+ const Stencil stencil = StencilBuilder{}.build(mesh.connectivity());
+
+ auto xr = mesh.xr();
+ auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
+
+ auto cell_is_owned = mesh.connectivity().cellIsOwned();
+
+ DiscreteFunctionDPk<MeshType::Dimension, std::remove_const_t<DataType>> dpk{mesh.meshVariant(), degree};
+
+ if constexpr ((is_polygonal_mesh_v<MeshType>)and(MeshType::Dimension == 1)) {
+ auto qf = QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor(1));
+
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_j_id) {
+ if (cell_is_owned[cell_j_id]) {
+ auto stencil_cell_list = stencil[cell_j_id];
+ SmallVector<double> b(stencil_cell_list.size());
+
+ const double p_j = p0_function[cell_j_id];
+ for (size_t i = 0; i < stencil_cell_list.size(); ++i) {
+ const CellId cell_i_id = stencil_cell_list[i];
+
+ b[i] = p0_function[cell_i_id] - p_j;
+ }
+
+ SmallMatrix<double> A(stencil_cell_list.size(), degree);
+
+ using R1 = TinyVector<1>;
+
+ const R1& Xj = xj[cell_j_id];
+ auto X_Xj = [&](const R1& X) { return X - Xj; };
+
+ for (size_t i = 0; i < stencil_cell_list.size(); ++i) {
+ const CellId cell_i_id = stencil_cell_list[i];
+
+ A(i, 0) = X_Xj(xj[cell_i_id])[0];
+ }
+
+ SmallVector<double> x(1);
+ Givens::solveInPlace(A, x, b);
+
+ dpk.coefficients(cell_j_id)[0] = p_j;
+ dpk.coefficients(cell_j_id)[1] = x[0];
+ }
+ });
+ }
+
+ synchronize(dpk.cellArrays());
+ return dpk;
+ }
+
+ PolynomialReconstruction() {}
+};
+
+#endif // POLYNOMIAL_RECONSTRUCTION_HPP
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 68da0dc88f3ab393e8e7623bc4580fff21161b6e..fccafd015e31ced8c3cf5981b216a4015e15b286 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -96,6 +96,7 @@ add_executable (unit_tests
test_GaussLegendreQuadratureDescriptor.cpp
test_GaussLobattoQuadratureDescriptor.cpp
test_GaussQuadratureDescriptor.cpp
+ test_Givens.cpp
test_IfProcessor.cpp
test_IncDecExpressionProcessor.cpp
test_IntegrateCellArray.cpp
@@ -203,6 +204,7 @@ add_executable (mpi_unit_tests
test_OFStream.cpp
test_ParallelChecker_read.cpp
test_Partitioner.cpp
+ test_PolynomialReconstruction.cpp
test_RandomEngine.cpp
test_StencilBuilder.cpp
test_SubItemArrayPerItemVariant.cpp
diff --git a/tests/test_Givens.cpp b/tests/test_Givens.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..818d8222278258fa9323e31925c7d584f471777a
--- /dev/null
+++ b/tests/test_Givens.cpp
@@ -0,0 +1,96 @@
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_all.hpp>
+
+#include <algebra/Givens.hpp>
+#include <algebra/SmallMatrix.hpp>
+#include <algebra/SmallVector.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("Givens", "[algebra]")
+{
+ SECTION("square matrix")
+ {
+ SmallMatrix<double> A{5, 5};
+ A.fill(0);
+ A(0, 0) = 2;
+ A(0, 1) = -1;
+
+ A(1, 0) = -0.2;
+ A(1, 1) = 2;
+ A(1, 2) = -1;
+
+ A(2, 1) = -1;
+ A(2, 2) = 4;
+ A(2, 3) = -2;
+
+ A(3, 2) = -1;
+ A(3, 3) = 2;
+ A(3, 4) = -0.1;
+
+ A(4, 3) = 1;
+ A(4, 4) = 3;
+
+ // SmallMatrix<double> A{S.getMatrix()};
+
+ SmallVector<const double> x_exact = [] {
+ SmallVector<double> y{5};
+ y[0] = 1;
+ y[1] = 3;
+ y[2] = 2;
+ y[3] = 4;
+ y[4] = 5;
+ return y;
+ }();
+
+ SmallVector<double> b = A * x_exact;
+
+ SmallVector<double> x{5};
+ x = zero;
+
+ Givens::solve(A, x, b);
+ SmallVector<double> error = x - x_exact;
+
+ REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x, x)));
+ }
+
+ SECTION("rectangular matrix")
+ {
+ SmallMatrix<double> A{5, 3};
+ A.fill(0);
+ A(0, 0) = 2;
+ A(0, 1) = -1;
+
+ A(1, 0) = -0.2;
+ A(1, 1) = 2;
+ A(1, 2) = -1;
+
+ A(2, 1) = -1;
+ A(2, 2) = 4;
+
+ A(3, 2) = -1;
+ A(4, 0) = 0.5;
+ A(4, 1) = 1;
+ A(4, 2) = 1;
+
+ // SmallMatrix<double> A{S.getMatrix()};
+
+ SmallVector<const double> x_exact = [] {
+ SmallVector<double> y{3};
+ y[0] = 1;
+ y[1] = 3;
+ y[2] = 2;
+ return y;
+ }();
+
+ SmallVector<double> b = A * x_exact;
+
+ SmallVector<double> x{3};
+ x = zero;
+
+ Givens::solve(A, x, b);
+ SmallVector<double> error = x - x_exact;
+
+ REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x, x)));
+ }
+}
diff --git a/tests/test_PolynomialReconstruction.cpp b/tests/test_PolynomialReconstruction.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..9a1df325788e1277b4e0996c0316bb0fe9081ce1
--- /dev/null
+++ b/tests/test_PolynomialReconstruction.cpp
@@ -0,0 +1,65 @@
+#include <catch2/catch_approx.hpp>
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_all.hpp>
+
+#include <Kokkos_Core.hpp>
+
+#include <utils/PugsAssert.hpp>
+#include <utils/Types.hpp>
+
+#include <algebra/SmallMatrix.hpp>
+#include <algebra/SmallVector.hpp>
+#include <mesh/Mesh.hpp>
+#include <mesh/MeshDataManager.hpp>
+#include <scheme/DiscreteFunctionP0.hpp>
+#include <scheme/PolynomialReconstruction.hpp>
+
+#include <MeshDataBaseForTests.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("PolynomialReconstruction", "[scheme]")
+{
+ SECTION("1D")
+ {
+ using R1 = TinyVector<1>;
+
+ for (auto named_mesh : MeshDataBaseForTests::get().all1DMeshes()) {
+ SECTION(named_mesh.name())
+ {
+ auto p_mesh = named_mesh.mesh()->get<Mesh<1>>();
+ auto& mesh = *p_mesh;
+
+ auto affine = [](const R1& x) { return 2.3 + 1.7 * x[0]; };
+ auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
+
+ DiscreteFunctionP0<double> fh{p_mesh};
+
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = affine(xj[cell_id]); });
+
+ PolynomialReconstruction reconstruction;
+ auto dpk = reconstruction.build(mesh, fh);
+
+ {
+ double max_mean_error = 0;
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ max_mean_error = std::max(max_mean_error, std::abs(dpk[cell_id](xj[cell_id]) - affine(xj[cell_id])));
+ }
+ REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-14));
+ }
+
+ {
+ double max_slope_error = 0;
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ const double reconstructed_slope =
+ (dpk[cell_id](R1{0.1} + xj[cell_id]) - dpk[cell_id](xj[cell_id] - R1{0.1})) / 0.2;
+
+ max_slope_error = std::max(max_slope_error, std::abs(reconstructed_slope - 1.7));
+ }
+ REQUIRE(max_slope_error == Catch::Approx(0).margin(1E-14));
+ }
+ }
+ }
+ }
+}