Fix degree 1 reconstruction for discrete function p0 vectors

Also add tests

src/scheme/PolynomialReconstruction.cpp

@@ -122,12 +122,16 @@ PolynomialReconstruction::_build(
             } else if constexpr (is_tiny_vector_v<data_type> or is_tiny_matrix_v<data_type>) {
               return data_type::Dimension;
             } else {
+              // LCOV_EXCL_START
               throw UnexpectedError("unexpected data type " + demangle<data_type>());
+              // LCOV_EXCL_STOP
             }
           } else if constexpr (is_discrete_function_P0_vector_v<DiscreteFunctionT>) {
             return discrete_function.size();
           } else {
+            // LCOV_EXCL_START
             throw UnexpectedError("unexpected discrete function type");
+            // LCOV_EXCL_STOP
           }
         },
         discrete_function_variant_p->discreteFunction());
@@ -171,9 +175,11 @@ PolynomialReconstruction::_build(
         } else if constexpr (is_discrete_function_P0_vector_v<DiscreteFunctionT>) {
           using DataType = std::remove_const_t<std::decay_t<typename DiscreteFunctionT::data_type>>;
           mutable_discrete_function_dpk_variant_list.push_back(
-            DiscreteFunctionDPkVector<MeshType::Dimension, DataType>(p_mesh, discrete_function.size(), degree));
+            DiscreteFunctionDPkVector<MeshType::Dimension, DataType>(p_mesh, degree, discrete_function.size()));
         } else {
+          // LCOV_EXCL_START
           throw UnexpectedError("unexpected discrete function type");
+          // LCOV_EXCL_STOP
         }
       },
       discrete_function_variant->discreteFunction());
@@ -245,7 +251,9 @@ PolynomialReconstruction::_build(
                 }
                 column_begin += pj_vector.size();
               } else {
+                // LCOV_EXCL_START
                 throw UnexpectedError("invalid discrete function type");
+                // LCOV_EXCL_STOP
               }
             },
             discrete_function_variant->discreteFunction());
@@ -282,6 +290,7 @@ PolynomialReconstruction::_build(
 
               if constexpr (std::is_same_v<DataType, P0DataType>) {
                 if constexpr (is_discrete_function_P0_v<P0FunctionT>) {
+                  if constexpr (is_discrete_function_dpk_scalar_v<DPkFunctionT>) {
                     auto dpk_j = dpk_function.coefficients(cell_j_id);
                     dpk_j[0]   = p0_function[cell_j_id];
 
@@ -310,9 +319,17 @@ PolynomialReconstruction::_build(
                       }
                       column_begin += DataType::Dimension;
                     } else {
+                      // LCOV_EXCL_START
                       throw UnexpectedError("unexpected data type");
+                      // LCOV_EXCL_STOP
+                    }
+                  } else {
+                    // LCOV_EXCL_START
+                    throw UnexpectedError("unexpected discrete dpk function type");
+                    // LCOV_EXCL_STOP
                   }
                 } else if constexpr (is_discrete_function_P0_vector_v<P0FunctionT>) {
+                  if constexpr (is_discrete_function_dpk_vector_v<DPkFunctionT>) {
                     auto dpk_j        = dpk_function.coefficients(cell_j_id);
                     auto cell_vector  = p0_function[cell_j_id];
                     const size_t size = MeshType::Dimension + 1;
@@ -327,14 +344,25 @@ PolynomialReconstruction::_build(
                         }
                         ++column_begin;
                       } else {
+                        // LCOV_EXCL_START
                         throw UnexpectedError("unexpected data type");
+                        // LCOV_EXCL_STOP
                       }
                     }
                   } else {
+                    // LCOV_EXCL_START
+                    throw UnexpectedError("unexpected discrete dpk function type");
+                    // LCOV_EXCL_STOP
+                  }
+                } else {
+                  // LCOV_EXCL_START
                   throw UnexpectedError("unexpected discrete function type");
+                  // LCOV_EXCL_STOP
                 }
               } else {
+                // LCOV_EXCL_START
                 throw UnexpectedError("incompatible data types");
+                // LCOV_EXCL_STOP
               }
             },
             dpk_variant.mutableDiscreteFunctionDPk(), discrete_function_variant->discreteFunction());
@@ -364,7 +392,7 @@ PolynomialReconstruction::build(
   const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const
 {
   if (not hasSameMesh(discrete_function_variant_list)) {
-    throw NormalError("cannot reconstruct functions living of different mesh simultaneously");
+    throw NormalError("cannot reconstruct functions living of different meshes simultaneously");
   }
 
   auto mesh_v = getCommonMesh(discrete_function_variant_list);

tests/test_PolynomialReconstruction.cpp

@@ -171,6 +171,64 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
       }
     }
 
+    SECTION("R 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_affine = [](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 xj = MeshDataManager::instance().getMeshData(mesh).xj();
+
+          DiscreteFunctionP0Vector<double> Vh{p_mesh, 3};
+
+          parallel_for(
+            mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+              auto vector = vector_affine(xj[cell_id]);
+              for (size_t i = 0; i < vector.dimension(); ++i) {
+                Vh[cell_id][i] = vector[i];
+              }
+            });
+
+          auto reconstructions = PolynomialReconstruction{}.build(Vh);
+
+          auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<1, const double>>();
+
+          {
+            double max_mean_error = 0;
+            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+              auto vector = vector_affine(xj[cell_id]);
+              for (size_t i = 0; i < vector.dimension(); ++i) {
+                max_mean_error = std::max(max_mean_error, std::abs(dpk_Vh(cell_id, i)(xj[cell_id]) - vector[i]));
+              }
+            }
+            REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-13));
+          }
+
+          {
+            double max_slope_error = 0;
+            const TinyVector<3> slope{+1.7, +2.1, -0.6};
+
+            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+              for (size_t i = 0; i < slope.dimension(); ++i) {
+                const double reconstructed_slope =
+                  (1 / 0.2) * (dpk_Vh(cell_id, i)(R1{0.1} + xj[cell_id]) - dpk_Vh(cell_id, i)(xj[cell_id] - R1{0.1}));
+
+                max_slope_error = std::max(max_slope_error, std::abs(reconstructed_slope - slope[i]));
+              }
+              REQUIRE(max_slope_error == Catch::Approx(0).margin(1E-13));
+            }
+          }
+        }
+      }
+    }
+
     SECTION("list of various types")
     {
       using R3x3 = TinyMatrix<3>;
@@ -198,6 +256,10 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
             };
           };
 
+          auto vector_affine = [](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 xj = MeshDataManager::instance().getMeshData(mesh).xj();
 
           DiscreteFunctionP0<double> fh{p_mesh};
@@ -212,8 +274,18 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
           parallel_for(
             mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { Ah[cell_id] = R3x3_affine(xj[cell_id]); });
 
+          DiscreteFunctionP0Vector<double> Vh{p_mesh, 3};
+          parallel_for(
+            mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+              auto vector = vector_affine(xj[cell_id]);
+              for (size_t i = 0; i < vector.dimension(); ++i) {
+                Vh[cell_id][i] = vector[i];
+              }
+            });
+
           auto reconstructions = PolynomialReconstruction{}.build(std::make_shared<DiscreteFunctionVariant>(fh), uh,
-                                                                  std::make_shared<DiscreteFunctionP0<R3x3>>(Ah));
+                                                                  std::make_shared<DiscreteFunctionP0<R3x3>>(Ah),
+                                                                  DiscreteFunctionVariant(Vh));
 
           auto dpk_fh = reconstructions[0]->get<DiscreteFunctionDPk<1, const double>>();
 
@@ -283,6 +355,34 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
             }
             REQUIRE(max_slope_error == Catch::Approx(0).margin(1E-13));
           }
+
+          auto dpk_Vh = reconstructions[3]->get<DiscreteFunctionDPkVector<1, const double>>();
+
+          {
+            double max_mean_error = 0;
+            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+              auto vector = vector_affine(xj[cell_id]);
+              for (size_t i = 0; i < vector.dimension(); ++i) {
+                max_mean_error = std::max(max_mean_error, std::abs(dpk_Vh(cell_id, i)(xj[cell_id]) - vector[i]));
+              }
+            }
+            REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-13));
+          }
+
+          {
+            double max_slope_error = 0;
+            const TinyVector<3> slope{+1.7, +2.1, -0.6};
+
+            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+              for (size_t i = 0; i < slope.dimension(); ++i) {
+                const double reconstructed_slope =
+                  (1 / 0.2) * (dpk_Vh(cell_id, i)(R1{0.1} + xj[cell_id]) - dpk_Vh(cell_id, i)(xj[cell_id] - R1{0.1}));
+
+                max_slope_error = std::max(max_slope_error, std::abs(reconstructed_slope - slope[i]));
+              }
+              REQUIRE(max_slope_error == Catch::Approx(0).margin(1E-13));
+            }
+          }
         }
       }
     }
@@ -464,6 +564,78 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
         }
       }
     }
+
+    SECTION("vector data")
+    {
+      using R4 = TinyVector<4>;
+
+      for (auto named_mesh : MeshDataBaseForTests::get().all2DMeshes()) {
+        SECTION(named_mesh.name())
+        {
+          auto p_mesh = named_mesh.mesh()->get<Mesh<2>>();
+          auto& mesh  = *p_mesh;
+
+          auto vector_affine = [](const R2& x) -> R4 {
+            return R4{+2.3 + 1.7 * x[0] + 1.2 * x[1], -1.7 + 2.1 * x[0] - 2.2 * x[1],   //
+                      +1.4 - 0.6 * x[0] - 2.1 * x[1], +2.4 - 2.3 * x[0] + 1.3 * x[1]};
+          };
+          auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
+
+          DiscreteFunctionP0Vector<double> Vh{p_mesh, 4};
+
+          parallel_for(
+            mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+              auto vector = vector_affine(xj[cell_id]);
+              for (size_t i = 0; i < vector.dimension(); ++i) {
+                Vh[cell_id][i] = vector[i];
+              }
+            });
+
+          auto reconstructions = PolynomialReconstruction{}.build(Vh);
+
+          auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<2, const double>>();
+
+          {
+            double max_mean_error = 0;
+            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+              auto vector = vector_affine(xj[cell_id]);
+              for (size_t i = 0; i < vector.dimension(); ++i) {
+                max_mean_error = std::max(max_mean_error, std::abs(dpk_Vh(cell_id, i)(xj[cell_id]) - vector[i]));
+              }
+            }
+            REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-13));
+          }
+
+          {
+            double max_x_slope_error = 0;
+            const R4 slope{+1.7, +2.1, -0.6, -2.3};
+            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+              for (size_t i = 0; i < slope.dimension(); ++i) {
+                const double reconstructed_slope = (1 / 0.2) * (dpk_Vh(cell_id, i)(R2{0.1, 0} + xj[cell_id]) -
+                                                                dpk_Vh(cell_id, i)(xj[cell_id] - R2{0.1, 0}));
+
+                max_x_slope_error = std::max(max_x_slope_error, std::abs(reconstructed_slope - slope[i]));
+              }
+            }
+            REQUIRE(max_x_slope_error == Catch::Approx(0).margin(1E-12));
+          }
+
+          {
+            double max_y_slope_error = 0;
+            const R4 slope{+1.2, -2.2, -2.1, +1.3};
+            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+              for (size_t i = 0; i < slope.dimension(); ++i) {
+                const double reconstructed_slope = (1 / 0.2) * (dpk_Vh(cell_id, i)(R2{0, 0.1} + xj[cell_id]) -
+                                                                dpk_Vh(cell_id, i)(xj[cell_id] - R2{0, 0.1}));
+
+                max_y_slope_error = std::max(max_y_slope_error, std::abs(reconstructed_slope - slope[i]));
+              }
+            }
+            REQUIRE(max_y_slope_error == Catch::Approx(0).margin(1E-12));
+          }
+        }
+      }
+    }
   }
 
   SECTION("3D")
@@ -675,5 +847,104 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
         }
       }
     }
+
+    SECTION("vector data")
+    {
+      using R4 = TinyVector<4>;
+
+      for (auto named_mesh : MeshDataBaseForTests::get().all3DMeshes()) {
+        SECTION(named_mesh.name())
+        {
+          auto p_mesh = named_mesh.mesh()->get<Mesh<3>>();
+          auto& mesh  = *p_mesh;
+
+          auto vector_affine = [](const R3& x) -> R4 {
+            return R4{+2.3 + 1.7 * x[0] + 1.2 * x[1] - 1.3 * x[2], -1.7 + 2.1 * x[0] - 2.2 * x[1] - 2.4 * x[2],
+                      //
+                      +2.4 - 2.3 * x[0] + 1.3 * x[1] + 1.4 * x[2], -0.2 + 3.1 * x[0] + 0.8 * x[1] - 1.8 * x[2]};
+          };
+          auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
+
+          DiscreteFunctionP0Vector<double> Vh{p_mesh, 4};
+
+          parallel_for(
+            mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+              auto vector = vector_affine(xj[cell_id]);
+              for (size_t i = 0; i < vector.dimension(); ++i) {
+                Vh[cell_id][i] = vector[i];
+              }
+            });
+
+          auto reconstructions = PolynomialReconstruction{}.build(Vh);
+
+          auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<3, const double>>();
+
+          {
+            double max_mean_error = 0;
+            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+              auto vector = vector_affine(xj[cell_id]);
+              for (size_t i = 0; i < vector.dimension(); ++i) {
+                max_mean_error = std::max(max_mean_error, std::abs(dpk_Vh(cell_id, i)(xj[cell_id]) - vector[i]));
+              }
+            }
+            REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-13));
+          }
+
+          {
+            double max_x_slope_error = 0;
+            const R4 slope{+1.7, 2.1, -2.3, +3.1};
+            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+              for (size_t i = 0; i < slope.dimension(); ++i) {
+                const double reconstructed_slope = (1 / 0.2) * (dpk_Vh(cell_id, i)(R3{0.1, 0, 0} + xj[cell_id]) -
+                                                                dpk_Vh(cell_id, i)(xj[cell_id] - R3{0.1, 0, 0}));
+
+                max_x_slope_error = std::max(max_x_slope_error, std::abs(reconstructed_slope - slope[i]));
+              }
+            }
+            REQUIRE(max_x_slope_error == Catch::Approx(0).margin(1E-12));
+          }
+
+          {
+            double max_y_slope_error = 0;
+            const R4 slope{+1.2, -2.2, 1.3, +0.8};
+            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+              for (size_t i = 0; i < slope.dimension(); ++i) {
+                const double reconstructed_slope = (1 / 0.2) * (dpk_Vh(cell_id, i)(R3{0, 0.1, 0} + xj[cell_id]) -
+                                                                dpk_Vh(cell_id, i)(xj[cell_id] - R3{0, 0.1, 0}));
+
+                max_y_slope_error = std::max(max_y_slope_error, std::abs(reconstructed_slope - slope[i]));
+              }
+            }
+            REQUIRE(max_y_slope_error == Catch::Approx(0).margin(1E-12));
+          }
+
+          {
+            double max_z_slope_error = 0;
+            const R4 slope{-1.3, -2.4, +1.4, -1.8};
+            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+              for (size_t i = 0; i < slope.dimension(); ++i) {
+                const double reconstructed_slope = (1 / 0.2) * (dpk_Vh(cell_id, i)(R3{0, 0, 0.1} + xj[cell_id]) -
+                                                                dpk_Vh(cell_id, i)(xj[cell_id] - R3{0, 0, 0.1}));
+
+                max_z_slope_error = std::max(max_z_slope_error, std::abs(reconstructed_slope - slope[i]));
+              }
+            }
+            REQUIRE(max_z_slope_error == Catch::Approx(0).margin(1E-12));
+          }
+        }
+      }
+    }
+  }
+
+  SECTION("errors")
+  {
+    auto p_mesh1 = MeshDataBaseForTests::get().unordered1DMesh()->get<Mesh<1>>();
+    DiscreteFunctionP0<double> f1{p_mesh1};
+
+    auto p_mesh2 = MeshDataBaseForTests::get().cartesian1DMesh()->get<Mesh<1>>();
+    DiscreteFunctionP0<double> f2{p_mesh2};
+
+    REQUIRE_THROWS_WITH(PolynomialReconstruction{}.build(f1, f2),
+                        "error: cannot reconstruct functions living of different meshes simultaneously");
   }
 }