From 68ef5ec9a822a8f7120e53e6fd51c9de76ec8bc7 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Del=20Pino?= <stephane.delpino44@gmail.com>
Date: Tue, 19 Nov 2024 19:23:03 +0100
Subject: [PATCH] [ci-skip] WIP: Add reconstrcution for
 DiscreteFunctionDPkVector of TinyVector

Almost works, remains to fix dispatching after reconstruction
---
 src/scheme/PolynomialReconstruction.cpp | 99 +++++++++++++++++++++----
 tests/test_PolynomialReconstruction.cpp | 80 ++++++++++++++++++++
 2 files changed, 164 insertions(+), 15 deletions(-)

diff --git a/src/scheme/PolynomialReconstruction.cpp b/src/scheme/PolynomialReconstruction.cpp
index d23644a72..783ec01ad 100644
--- a/src/scheme/PolynomialReconstruction.cpp
+++ b/src/scheme/PolynomialReconstruction.cpp
@@ -143,7 +143,13 @@ class PolynomialReconstruction::MutableDiscreteFunctionDPkVariant
   MutableDiscreteFunctionDPkVariant(const DiscreteFunctionDPkVector<Dimension, DataType>& discrete_function_dpk)
     : m_mutable_discrete_function_dpk{discrete_function_dpk}
   {
-    static_assert(std::is_same_v<DataType, double>,
+    static_assert(std::is_same_v<DataType, double> or                       //
+                    std::is_same_v<DataType, TinyVector<1, double>> or      //
+                    std::is_same_v<DataType, TinyVector<2, double>> or      //
+                    std::is_same_v<DataType, TinyVector<3, double>> or      //
+                    std::is_same_v<DataType, TinyMatrix<1, 1, double>> or   //
+                    std::is_same_v<DataType, TinyMatrix<2, 2, double>> or   //
+                    std::is_same_v<DataType, TinyMatrix<3, 3, double>>,
                   "DiscreteFunctionDPkVector with this DataType is not allowed in variant");
   }
 
@@ -718,7 +724,7 @@ PolynomialReconstruction::_getNumberOfColumns(
         using DiscreteFunctionT = std::decay_t<decltype(discrete_function)>;
         if constexpr (is_discrete_function_P0_v<DiscreteFunctionT>) {
           using data_type = std::decay_t<typename DiscreteFunctionT::data_type>;
-          if constexpr (std::is_same_v<data_type, double>) {
+          if constexpr (std::is_arithmetic_v<data_type>) {
             return 1;
           } else if constexpr (is_tiny_vector_v<data_type> or is_tiny_matrix_v<data_type>) {
             return data_type::Dimension;
@@ -728,7 +734,16 @@ PolynomialReconstruction::_getNumberOfColumns(
             // LCOV_EXCL_STOP
           }
         } else if constexpr (is_discrete_function_P0_vector_v<DiscreteFunctionT>) {
-          return discrete_function.size();
+          using data_type = std::decay_t<typename DiscreteFunctionT::data_type>;
+          if constexpr (std::is_arithmetic_v<data_type>) {
+            return discrete_function.size();
+          } else if constexpr (is_tiny_vector_v<data_type> or is_tiny_matrix_v<data_type>) {
+            return discrete_function.size() * data_type::Dimension;
+          } else {
+            // LCOV_EXCL_START
+            throw UnexpectedError("unexpected data type " + demangle<data_type>());
+            // LCOV_EXCL_STOP
+          }
         } else {
           // LCOV_EXCL_START
           throw UnexpectedError("unexpected discrete function type");
@@ -760,13 +775,9 @@ PolynomialReconstruction::_createMutableDiscreteFunctionDPKVariantList(
             DiscreteFunctionDPk<MeshType::Dimension, DataType>(p_mesh, m_descriptor.degree()));
         } else if constexpr (is_discrete_function_P0_vector_v<DiscreteFunctionT>) {
           using DataType = std::remove_const_t<std::decay_t<typename DiscreteFunctionT::data_type>>;
-          if constexpr (std::is_arithmetic_v<DataType>) {
-            mutable_discrete_function_dpk_variant_list.push_back(
-              DiscreteFunctionDPkVector<MeshType::Dimension, DataType>(p_mesh, m_descriptor.degree(),
-                                                                       discrete_function.size()));
-          } else {
-            throw NotImplementedError("reconstruction of DiscreteFunctionP0Vector of non arithmetic data type");
-          }
+          mutable_discrete_function_dpk_variant_list.push_back(
+            DiscreteFunctionDPkVector<MeshType::Dimension, DataType>(p_mesh, m_descriptor.degree(),
+                                                                     discrete_function.size()));
         } else {
           // LCOV_EXCL_START
           throw UnexpectedError("unexpected discrete function type");
@@ -951,7 +962,7 @@ PolynomialReconstruction::_build(
                     } else if constexpr (is_tiny_matrix_v<DataType>) {
                       if constexpr ((DataType::NumberOfColumns == DataType::NumberOfRows) and
                                     (DataType::NumberOfColumns == MeshType::Dimension)) {
-                        throw NotImplementedError("TinyMatrix symmetries for reconstruction");
+                        throw NotImplementedError("TinyMatrix symmetries for reconstruction of DiscreteFunctionP0");
                       }
                       const DataType& qi_qj = discrete_function[cell_i_id] - qj;
                       for (size_t k = 0; k < DataType::NumberOfRows; ++k) {
@@ -971,9 +982,10 @@ PolynomialReconstruction::_build(
               } else if constexpr (is_discrete_function_P0_vector_v<DiscreteFunctionT>) {
                 using DataType = std::decay_t<typename DiscreteFunctionT::data_type>;
 
+                const auto qj_vector = discrete_function[cell_j_id];
+
                 if constexpr (std::is_arithmetic_v<DataType>) {
-                  const auto qj_vector = discrete_function[cell_j_id];
-                  size_t index         = 0;
+                  size_t index = 0;
                   for (size_t i = 0; i < stencil_cell_list.size(); ++i, ++index) {
                     const CellId cell_i_id = stencil_cell_list[i];
                     for (size_t l = 0; l < qj_vector.size(); ++l) {
@@ -996,10 +1008,47 @@ PolynomialReconstruction::_build(
                       }
                     }
                   }
+                } else if constexpr (is_tiny_vector_v<DataType>) {
+                  size_t index = 0;
+                  for (size_t i = 0; i < stencil_cell_list.size(); ++i, ++index) {
+                    const CellId cell_i_id = stencil_cell_list[i];
+                    for (size_t l = 0; l < qj_vector.size(); ++l) {
+                      const DataType& qj    = qj_vector[l];
+                      const DataType& qi_qj = discrete_function[cell_i_id][l] - qj;
+                      for (size_t kB = column_begin + l * DataType::Dimension, k = 0; k < DataType::Dimension;
+                           ++k, ++kB) {
+                        B(index, kB) = qi_qj[k];
+                      }
+                    }
+                  }
+
+                  for (size_t i_symmetry = 0; i_symmetry < stencil_array.symmetryBoundaryStencilArrayList().size();
+                       ++i_symmetry) {
+                    auto& ghost_stencil  = stencil_array.symmetryBoundaryStencilArrayList()[i_symmetry].stencilArray();
+                    auto ghost_cell_list = ghost_stencil[cell_j_id];
+                    if (ghost_cell_list.size() > 0) {
+                      throw NotImplementedError("TinyVector symmetries for reconstruction of DiscreteFunctionP0Vector");
+                    }
+                  }
+                } else if constexpr (is_tiny_matrix_v<DataType>) {
+                  throw NotImplementedError("NIY TinyMatrix data");
+
+                  for (size_t i_symmetry = 0; i_symmetry < stencil_array.symmetryBoundaryStencilArrayList().size();
+                       ++i_symmetry) {
+                    auto& ghost_stencil  = stencil_array.symmetryBoundaryStencilArrayList()[i_symmetry].stencilArray();
+                    auto ghost_cell_list = ghost_stencil[cell_j_id];
+                    if (ghost_cell_list.size() > 0) {
+                      throw NotImplementedError("TinyMatrix symmetries for reconstruction of DiscreteFunctionP0Vector");
+                    }
+                  }
+                }
+
+                if constexpr (std::is_arithmetic_v<DataType>) {
                   column_begin += qj_vector.size();
-                } else {
-                  throw NotImplementedError("reconstruction of DiscreteFunctionP0Vector of non arithmetic data type");
+                } else if constexpr (is_tiny_vector_v<DataType> or is_tiny_matrix_v<DataType>) {
+                  column_begin += qj_vector.size() * DataType::Dimension;
                 }
+
               } else {
                 // LCOV_EXCL_START
                 throw UnexpectedError("invalid discrete function type");
@@ -1212,6 +1261,8 @@ PolynomialReconstruction::_build(
           Givens::solveCollectionInPlace(A, X, B);
         }
 
+        std::clog << "X = " << X << '\n';
+
         column_begin = 0;
         for (size_t i_dpk_variant = 0; i_dpk_variant < mutable_discrete_function_dpk_variant_list.size();
              ++i_dpk_variant) {
@@ -1318,6 +1369,24 @@ PolynomialReconstruction::_build(
                           dpk_j_ip1       = X(i, column_begin);
                         }
                         ++column_begin;
+                      } else if constexpr (is_tiny_vector_v<DataType>) {
+                        if (m_descriptor.degree() > 1) {
+                          auto& mean_j_of_ejk = mean_j_of_ejk_pool[t];
+                          for (size_t i = 0; i < basis_dimension - 1; ++i) {
+                            auto& dpk_j_0 = dpk_j[0];
+                            for (size_t k = 0; k < DataType::Dimension; ++k) {
+                              dpk_j_0[k] -= X(i, column_begin + k) * mean_j_of_ejk[i];
+                            }
+                          }
+                        }
+
+                        for (size_t i = 0; i < basis_dimension - 1; ++i) {
+                          auto& dpk_j_ip1 = dpk_j[i + 1];
+                          for (size_t k = 0; k < DataType::Dimension; ++k) {
+                            dpk_j_ip1[k] = X(i, column_begin + k);
+                          }
+                        }
+                        column_begin += DataType::Dimension;
                       } else {
                         // LCOV_EXCL_START
                         throw UnexpectedError("unexpected data type");
diff --git a/tests/test_PolynomialReconstruction.cpp b/tests/test_PolynomialReconstruction.cpp
index cdfb532f2..1c078129b 100644
--- a/tests/test_PolynomialReconstruction.cpp
+++ b/tests/test_PolynomialReconstruction.cpp
@@ -242,6 +242,86 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
             }
           }
 
+          SECTION("R3 vector data")
+          {
+            using R3 = TinyVector<3>;
+
+            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 vector_affine0 = [](const R1& x) -> R3 {
+                  return R3{+2.3 + 1.7 * x[0], -1.7 + 2.1 * x[0], +1.4 - 0.6 * x[0]};
+                };
+
+                auto vector_affine1 = [](const R1& x) -> R3 {
+                  return R3{+1.6 + 0.7 * x[0], -2.1 + 1.2 * x[0], +1.1 - 0.3 * x[0]};
+                };
+
+                auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
+
+                DiscreteFunctionP0Vector<R3> Vh{p_mesh, 2};
+
+                parallel_for(
+                  mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+                    Vh[cell_id][0] = vector_affine0(xj[cell_id]);
+                    Vh[cell_id][1] = vector_affine1(xj[cell_id]);
+                  });
+
+                auto reconstructions = PolynomialReconstruction{descriptor}.build(Vh);
+
+                auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<1, const R3>>();
+
+                {
+                  double max_mean_error = 0;
+                  for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+                    max_mean_error =
+                      std::max(max_mean_error, l2Norm(dpk_Vh(cell_id, 0)(xj[cell_id]) - vector_affine0(xj[cell_id])));
+                    max_mean_error =
+                      std::max(max_mean_error, l2Norm(dpk_Vh(cell_id, 1)(xj[cell_id]) - vector_affine1(xj[cell_id])));
+                  }
+                  REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-14));
+                }
+
+                {
+                  double max_slope_error = 0;
+                  {
+                    const TinyVector<3> slope0{+1.7, +2.1, -0.6};
+
+                    for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+                      for (size_t i = 0; i < R3::Dimension; ++i) {
+                        const double reconstructed_slope = (1 / 0.2) * (dpk_Vh(cell_id, 0)(R1{0.1} + xj[cell_id])[i] -
+                                                                        dpk_Vh(cell_id, 0)(xj[cell_id] - R1{0.1})[i]);
+
+                        std::clog << "reconstructed_slop0[" << i << "]=" << reconstructed_slope << " | expected slope0["
+                                  << i << "]=" << slope0[i] << '\n';
+                        max_slope_error = std::max(max_slope_error, std::abs(reconstructed_slope - slope0[i]));
+                      }
+                    }
+                  }
+
+                  {
+                    const TinyVector<3> slope1{+0.7, +1.2, -0.3};
+
+                    for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+                      for (size_t i = 0; i < R3::Dimension; ++i) {
+                        const double reconstructed_slope = (1 / 0.2) * (dpk_Vh(cell_id, 1)(R1{0.1} + xj[cell_id])[i] -
+                                                                        dpk_Vh(cell_id, 1)(xj[cell_id] - R1{0.1})[i]);
+
+                        std::clog << "reconstructed_slop1[" << i << "]=" << reconstructed_slope << " | expected slope1["
+                                  << i << "]=" << slope1[i] << '\n';
+                        max_slope_error = std::max(max_slope_error, std::abs(reconstructed_slope - slope1[i]));
+                      }
+                    }
+                  }
+                  REQUIRE(parallel::allReduceMax(max_slope_error) == Catch::Approx(0).margin(1E-14));
+                }
+              }
+            }
+          }
+
           SECTION("list of various types")
           {
             using R3x3 = TinyMatrix<3>;
-- 
GitLab