Add row weighting for least-square resolution

The idea is to add weight for the values that are provided by cells
that are close to the reconstruction

src/scheme/PolynomialReconstruction.cpp

@@ -278,6 +278,17 @@ PolynomialReconstruction::_build(
           }
         }
 
+        for (size_t i = 0; i < stencil_cell_list.size(); ++i) {
+          const CellId cell_i_id = stencil_cell_list[i];
+          const double weight    = 1. / l2Norm(xj[cell_i_id] - Xj);
+          for (size_t l = 0; l < basis_dimension - 1; ++l) {
+            A(i, l) *= weight;
+          }
+          for (size_t l = 0; l < number_of_columns; ++l) {
+            B(i, l) *= weight;
+          }
+        }
+
         const SmallMatrix<double>& X = X_pool[t];
         Givens::solveCollectionInPlace(A, X, B);
 

src/scheme/PolynomialReconstruction.hpp

@@ -100,6 +100,15 @@ class PolynomialReconstruction
                 }
               }
 
+              for (size_t i = 0; i < stencil_cell_list.size(); ++i) {
+                const CellId cell_i_id = stencil_cell_list[i];
+                const double weight    = 1. / l2Norm(xj[cell_i_id] - Xj);
+                for (size_t l = 0; l < MeshType::Dimension; ++l) {
+                  A(i, l) *= weight;
+                }
+                b[i] *= weight;
+              }
+
               SmallVector<double> x(MeshType::Dimension);
               Givens::solveInPlace(A, x, b);
 
@@ -150,6 +159,17 @@ class PolynomialReconstruction
                 }
               }
 
+              for (size_t i = 0; i < stencil_cell_list.size(); ++i) {
+                const CellId cell_i_id = stencil_cell_list[i];
+                const double weight    = 1. / l2Norm(xj[cell_i_id] - Xj);
+                for (size_t l = 0; l < MeshType::Dimension - 1; ++l) {
+                  A(i, l) *= weight;
+                }
+                for (size_t l = 0; l < DataType::Dimension; ++l) {
+                  B(i, l) *= weight;
+                }
+              }
+
               TinyMatrix<MeshType::Dimension, DataType::Dimension> X;
               Givens::solveCollectionInPlace(A, X, B);
 
@@ -205,6 +225,17 @@ class PolynomialReconstruction
                 }
               }
 
+              for (size_t i = 0; i < stencil_cell_list.size(); ++i) {
+                const CellId cell_i_id = stencil_cell_list[i];
+                const double weight    = 1. / l2Norm(xj[cell_i_id] - Xj);
+                for (size_t l = 0; l < MeshType::Dimension - 1; ++l) {
+                  A(i, l) *= weight;
+                }
+                for (size_t l = 0; l < DataType::Dimension; ++l) {
+                  B(i, l) *= weight;
+                }
+              }
+
               TinyMatrix<MeshType::Dimension, DataType::Dimension> X;
               Givens::solveCollectionInPlace(A, X, B);
 

tests/test_PolynomialReconstruction.cpp

@@ -506,7 +506,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
 
                 max_y_slope_error = std::max(max_y_slope_error, l2Norm(reconstructed_slope - R3{-2.2, 1.3, -1.1}));
               }
-              REQUIRE(max_y_slope_error == Catch::Approx(0).margin(1E-14));
+              REQUIRE(max_y_slope_error == Catch::Approx(0).margin(1E-13));
             }
           }
         }