Begin polynomial cell integration in 2D using divergence theorem

This allows to compute integral of centered canonical basis using
integrals over lines and thus valid for any polygonal cell types

src/language/modules/SchemeModule.cpp

@@ -682,6 +682,18 @@ SchemeModule::SchemeModule()
 
                               ));
 
+#warning REMOVE AFTER TESTS FINISHED
+  this->_addBuiltinFunction("test_reconstruction",
+                            std::function(
+
+                              [](const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>&
+                                   discrete_function_variant_list,
+                                 uint64_t degree) -> void {
+                                test_reconstruction(discrete_function_variant_list, degree);
+                              }
+
+                              ));
+
   MathFunctionRegisterForVh{*this};
 }
 

src/scheme/PolynomialReconstruction.cpp

@@ -270,8 +270,7 @@ _computeEjkMean(const QuadratureFormula<3>& quadrature,
 }
 
 template <typename NodeListT, size_t Dimension>
-void
-_computeEjkMean(const TinyVector<Dimension>& Xj,
+void _computeEjkMean(const TinyVector<Dimension>& Xj,
                      const NodeValue<const TinyVector<Dimension>>& xr,
                      const NodeListT& node_list,
                      const CellType& cell_type,
@@ -279,8 +278,7 @@ _computeEjkMean(const TinyVector<Dimension>& Xj,
                      const size_t degree,
                      const size_t basis_dimension,
                      SmallArray<double>& inv_Vj_wq_detJ_ek,
-                SmallArray<double>& mean_of_ejk)
-{}
+                     SmallArray<double>& mean_of_ejk);
 
 template <typename NodeListT>
 void
@@ -306,6 +304,86 @@ _computeEjkMean(const TinyVector<1>& Xj,
   }
 }
 
+template <typename ConformTransformationT>
+void
+_computeEjkBoundaryMean(const QuadratureFormula<1>& quadrature,
+                        const ConformTransformationT& T,
+                        const TinyVector<2>& Xj,
+                        const double Vi,
+                        const size_t degree,
+                        const size_t dimension,
+                        SmallArray<double>& inv_Vj_alpha_p_1_wq_X_prime_orth_ek,
+                        SmallArray<double>& mean_of_ejk)
+{
+  using Rd = TinyVector<2>;
+
+  const double velocity_perp_e1 = T.velocity()[1];
+
+  for (size_t i_q = 0; i_q < quadrature.numberOfPoints(); ++i_q) {
+    const double wq          = quadrature.weight(i_q);
+    const TinyVector<1> xi_q = quadrature.point(i_q);
+
+    const Rd X_Xj = T(xi_q) - Xj;
+
+    const double x_xj = X_Xj[0];
+    const double y_yj = X_Xj[1];
+
+    {
+      size_t k                                 = 0;
+      inv_Vj_alpha_p_1_wq_X_prime_orth_ek[k++] = x_xj * wq * velocity_perp_e1 / Vi;
+      for (; k <= degree; ++k) {
+        inv_Vj_alpha_p_1_wq_X_prime_orth_ek[k] =
+          x_xj * inv_Vj_alpha_p_1_wq_X_prime_orth_ek[k - 1] * ((1. * k) / (k + 1));
+      }
+
+      for (size_t i_y = 1; i_y <= degree; ++i_y) {
+        const size_t begin_i_y_1 = ((i_y - 1) * (2 * degree - i_y + 4)) / 2;
+        for (size_t l = 0; l <= degree - i_y; ++l) {
+          inv_Vj_alpha_p_1_wq_X_prime_orth_ek[k++] = y_yj * inv_Vj_alpha_p_1_wq_X_prime_orth_ek[begin_i_y_1 + l];
+        }
+      }
+    }
+
+    for (size_t k = 1; k < dimension; ++k) {
+      mean_of_ejk[k - 1] += inv_Vj_alpha_p_1_wq_X_prime_orth_ek[k];
+    }
+  }
+}
+
+void
+_computeEjkMeanByBoundary(const TinyVector<2>& Xj,
+                          const CellId& cell_id,
+                          const NodeValue<const TinyVector<2>>& xr,
+                          const auto& cell_to_face_matrix,
+                          const auto& face_to_node_matrix,
+                          const auto& cell_face_is_reversed,
+                          const CellValue<const double>& Vi,
+                          const size_t degree,
+                          const size_t basis_dimension,
+                          SmallArray<double>& inv_Vj_alpha_p_1_wq_X_prime_orth_ek,
+                          SmallArray<double>& mean_of_ejk)
+{
+  const auto& quadrature = QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor(degree + 1));
+
+  mean_of_ejk.fill(0);
+  auto cell_face_list = cell_to_face_matrix[cell_id];
+  for (size_t i_face = 0; i_face < cell_face_list.size(); ++i_face) {
+    bool is_reversed = cell_face_is_reversed[cell_id][i_face];
+
+    const FaceId face_id = cell_face_list[i_face];
+    auto face_node_list  = face_to_node_matrix[face_id];
+    if (is_reversed) {
+      const LineTransformation<2> T{xr[face_node_list[1]], xr[face_node_list[0]]};
+      _computeEjkBoundaryMean(quadrature, T, Xj, Vi[cell_id], degree, basis_dimension,
+                              inv_Vj_alpha_p_1_wq_X_prime_orth_ek, mean_of_ejk);
+    } else {
+      const LineTransformation<2> T{xr[face_node_list[0]], xr[face_node_list[1]]};
+      _computeEjkBoundaryMean(quadrature, T, Xj, Vi[cell_id], degree, basis_dimension,
+                              inv_Vj_alpha_p_1_wq_X_prime_orth_ek, mean_of_ejk);
+    }
+  }
+}
+
 template <typename NodeListT>
 void
 _computeEjkMean(const TinyVector<2>& Xj,
@@ -481,6 +559,7 @@ PolynomialReconstruction::_build(
   auto cell_is_owned       = mesh.connectivity().cellIsOwned();
   auto cell_type           = mesh.connectivity().cellType();
   auto cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+  auto cell_to_face_matrix = mesh.connectivity().cellToFaceMatrix();
 
   const size_t max_stencil_size = [&]() {
     size_t max_size = 0;
@@ -509,6 +588,9 @@ PolynomialReconstruction::_build(
   SmallArray<SmallArray<double>> mean_j_of_ejk_pool(Kokkos::DefaultExecutionSpace::concurrency());
   SmallArray<SmallArray<double>> mean_i_of_ejk_pool(Kokkos::DefaultExecutionSpace::concurrency());
 
+  SmallArray<SmallArray<double>> inv_Vj_alpha_p_1_wq_X_prime_orth_ek_pool(Kokkos::DefaultExecutionSpace::concurrency());
+  SmallArray<SmallArray<double>> mean_l_of_ejk_pool(Kokkos::DefaultExecutionSpace::concurrency());
+
   for (size_t i = 0; i < A_pool.size(); ++i) {
     A_pool[i] = SmallMatrix<double>(max_stencil_size, basis_dimension - 1);
     B_pool[i] = SmallMatrix<double>(max_stencil_size, number_of_columns);
@@ -518,15 +600,17 @@ PolynomialReconstruction::_build(
     inv_Vj_wq_detJ_ek_pool[i] = SmallArray<double>(basis_dimension);
     mean_j_of_ejk_pool[i]     = SmallArray<double>(basis_dimension - 1);
     mean_i_of_ejk_pool[i]     = SmallArray<double>(basis_dimension - 1);
+
+    inv_Vj_alpha_p_1_wq_X_prime_orth_ek_pool[i] = SmallArray<double>(basis_dimension);
   }
 
 #warning remove after rework
-  const bool FORCE_INTEGRATE = []() {
+  [[maybe_unused]] const int INTEGRATE = []() {
     auto value = std::getenv("INTEGRATE");
     if (value == nullptr) {
-      return false;
+      return 0;
     } else {
-      return true;
+      return std::atoi(value);
     }
   }();
 
@@ -596,7 +680,7 @@ PolynomialReconstruction::_build(
 
         ShrinkMatrixView A(A_pool[t], stencil_cell_list.size());
 
-        if ((m_descriptor.degree() == 1) and (not FORCE_INTEGRATE)) {
+        if ((m_descriptor.degree() == 1) and (INTEGRATE == 0)) {
           const Rd& Xj = xj[cell_j_id];
           for (size_t i = 0; i < stencil_cell_list.size(); ++i) {
             const CellId cell_i_id = stencil_cell_list[i];
@@ -605,6 +689,36 @@ PolynomialReconstruction::_build(
               A(i, l) = Xi_Xj[l];
             }
           }
+
+        } else if ((INTEGRATE == 2) or (INTEGRATE == 1)) {
+          if ((INTEGRATE == 2) and (MeshType::Dimension == 2)) {
+            if constexpr (MeshType::Dimension == 2) {
+              SmallArray<double>& inv_Vj_alpha_p_1_wq_X_prime_orth_ek = inv_Vj_alpha_p_1_wq_X_prime_orth_ek_pool[t];
+              SmallArray<double>& mean_j_of_ejk                       = mean_j_of_ejk_pool[t];
+              SmallArray<double>& mean_i_of_ejk                       = mean_i_of_ejk_pool[t];
+
+              auto face_to_node_matrix   = p_mesh->connectivity().faceToNodeMatrix();
+              auto cell_face_is_reversed = p_mesh->connectivity().cellFaceIsReversed();
+              const Rd& Xj               = xj[cell_j_id];
+
+              _computeEjkMeanByBoundary(Xj, cell_j_id, xr, cell_to_face_matrix, face_to_node_matrix,
+                                        cell_face_is_reversed, Vj, m_descriptor.degree(), basis_dimension,
+                                        inv_Vj_alpha_p_1_wq_X_prime_orth_ek, mean_j_of_ejk);
+
+              for (size_t i = 0; i < stencil_cell_list.size(); ++i) {
+                const CellId cell_i_id = stencil_cell_list[i];
+
+                _computeEjkMeanByBoundary(Xj, cell_i_id, xr, cell_to_face_matrix, face_to_node_matrix,
+                                          cell_face_is_reversed, Vj, m_descriptor.degree(), basis_dimension,
+                                          inv_Vj_alpha_p_1_wq_X_prime_orth_ek, mean_i_of_ejk);
+
+                for (size_t l = 0; l < basis_dimension - 1; ++l) {
+                  A(i, l) = mean_i_of_ejk[l] - mean_j_of_ejk[l];
+                }
+              }
+            } else {
+              throw UnexpectedError("invalid mesh dimension");
+            }
           } else {
             SmallArray<double>& inv_Vj_wq_detJ_ek = inv_Vj_wq_detJ_ek_pool[t];
             SmallArray<double>& mean_j_of_ejk     = mean_j_of_ejk_pool[t];
@@ -626,6 +740,9 @@ PolynomialReconstruction::_build(
               }
             }
           }
+        } else {
+          throw UnexpectedError("invalid integration strategy");
+        }
 
         if (m_descriptor.rowWeighting()) {
           // Add row weighting (give more importance to cells that are
@@ -679,12 +796,6 @@ PolynomialReconstruction::_build(
           Givens::solveCollectionInPlace(A, X, B);
         }
 
-        // Results are now written from the {ejk} basis to the {ek} basis
-        if (m_descriptor.degree() > 1) {
-          throw NotImplementedError(
-            "write the result from the {ejk} basis to the {ek} basis to get correct polynomial values");
-        }
-
         column_begin = 0;
         for (size_t i_dpk_variant = 0; i_dpk_variant < mutable_discrete_function_dpk_variant_list.size();
              ++i_dpk_variant) {

src/scheme/test_reconstruction.cpp

@@ -4,42 +4,17 @@
 #include <utils/Timer.hpp>
 
 void
-test_reconstruction(const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list)
+test_reconstruction(const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list,
+                    uint64_t degree)
 {
-  PolynomialReconstructionDescriptor descriptor{1};
-  // {
-  //   std::cout << "** variable list one by one (50 times)\n";
-  //   Timer t;
-  //   for (auto discrete_function_v : discrete_function_variant_list) {
-  //     std::visit(
-  //       [&](auto&& discrete_function) {
-  //         auto mesh_v = discrete_function.meshVariant();
-  //         std::visit(
-  //           [&](auto&& p_mesh) {
-  //             using DiscreteFunctionT = std::decay_t<decltype(discrete_function)>;
-  //             PolynomialReconstruction reconstruction{descriptor};
-  //             if constexpr (is_discrete_function_P0_v<DiscreteFunctionT>) {
-  //               for (size_t i = 0; i < 50; ++i) {
-  //                 auto res = reconstruction.build(*p_mesh, discrete_function);
-  //               }
-  //             } else if constexpr (is_discrete_function_P0_vector_v<DiscreteFunctionT>) {
-  //               std::cout << "Omitting discrete function p0 vector!\n";
-  //             }
-  //           },
-  //           mesh_v->variant());
-  //       },
-  //       discrete_function_v->discreteFunction());
-  //   }
-  //   t.pause();
-  //   std::cout << "t = " << t << '\n';
-  // }
-
+  PolynomialReconstructionDescriptor descriptor{degree};
   {
-    std::cout << "** variable list at once (50 times)\n";
+    const size_t nb_loops = 1;
+    std::cout << "** variable list at once (" << nb_loops << " times)\n";
 
     PolynomialReconstruction reconstruction{descriptor};
     Timer t;
-    for (size_t i = 0; i < 50; ++i) {
+    for (size_t i = 0; i < nb_loops; ++i) {
       auto res = reconstruction.build(discrete_function_variant_list);
     }
     t.pause();

src/scheme/test_reconstruction.hpp

@@ -8,6 +8,7 @@
 #include <vector>
 
 void test_reconstruction(
-  const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list);
+  const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list,
+  uint64_t degree = 1);
 
 #endif   // TEST_RECONSTRUCTION_HPP