Add tests for load balancer

src/scheme/LoadBalancer.cpp

@@ -44,6 +44,8 @@ LoadBalancer::balance(const std::vector<std::shared_ptr<const DiscreteFunctionVa
       },
       p_balanced_mesh_v->variant());
   } else {
+    // This macro test is used to avoid none reachable code for serial builds of pugs
+#ifdef PUGS_HAS_MPI
     std::visit(
       [&mesh_balancer, &balanced_discrete_function_list, &discrete_function_list](auto&& p_balanced_mesh) {
         using MeshType             = mesh_type_t<decltype(p_balanced_mesh)>;
@@ -69,6 +71,7 @@ LoadBalancer::balance(const std::vector<std::shared_ptr<const DiscreteFunctionVa
         }
       },
       p_balanced_mesh_v->variant());
+#endif   // PUGS_HAS_MPI
   }
 
   return balanced_discrete_function_list;

tests/CMakeLists.txt

@@ -134,6 +134,7 @@ add_executable (unit_tests
   test_LinearSolverOptions.cpp
   test_LineTransformation.cpp
   test_ListAffectationProcessor.cpp
+  test_LoadBalancer_serial.cpp
   test_MathModule.cpp
   test_MedianDualConnectivityBuilder.cpp
   test_MedianDualMeshBuilder.cpp
@@ -255,6 +256,7 @@ add_executable (mpi_unit_tests
   test_ItemValueUtils.cpp
   test_ItemValueVariant.cpp
   test_ItemValueVariantFunctionInterpoler.cpp
+  test_LoadBalancer_parallel.cpp
   test_MeshEdgeBoundary.cpp
   test_MeshEdgeInterface.cpp
   test_MeshFaceBoundary.cpp

tests/test_LoadBalancer_parallel.cpp

+#include <catch2/catch_all.hpp>
+#include <catch2/catch_test_macros.hpp>
+
+#include <scheme/LoadBalancer.hpp>
+
+#include <MeshDataBaseForTests.hpp>
+#include <mesh/Mesh.hpp>
+#include <mesh/MeshVariant.hpp>
+#include <scheme/DiscreteFunctionP0.hpp>
+#include <scheme/DiscreteFunctionP0Vector.hpp>
+#include <scheme/DiscreteFunctionVariant.hpp>
+#include <utils/PartitionerOptions.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("LoadBalancer parallel", "[scheme]")
+{
+  auto pugs_has_partitioner = [&](const PartitionerLibrary& partitioner_library) {
+    switch (partitioner_library) {
+    case PartitionerLibrary::parmetis: {
+#ifdef PUGS_HAS_PARMETIS
+      return true;
+#else    // PUGS_HAS_PARMETIS
+      return false;
+#endif   // PUGS_HAS_PARMETIS
+    }
+    case PartitionerLibrary::ptscotch: {
+#ifdef PUGS_HAS_PARMETIS
+      return true;
+#else    // PUGS_HAS_PARMETIS
+      return false;
+#endif   // PUGS_HAS_PARMETIS
+    }
+    default: {
+      return false;
+    }
+    }
+  };
+
+  constexpr std::array partitioner_libs = {PartitionerLibrary::parmetis, PartitionerLibrary::ptscotch};
+
+  const auto initial_library = PartitionerOptions::default_options.library();
+
+  for (auto lib : partitioner_libs) {
+    if (pugs_has_partitioner(lib)) {
+      SECTION(name(lib))
+      {
+        SECTION("1D")
+        {
+          for (auto&& named_mesh : MeshDataBaseForTests::get().all1DMeshes()) {
+            std::shared_ptr<const Mesh<1>> initial_mesh_1d = named_mesh.mesh()->get<Mesh<1>>();
+            SECTION(named_mesh.name())
+            {
+              DiscreteFunctionP0<double> initial_fh{initial_mesh_1d};
+              DiscreteFunctionP0Vector<double> initial_vh{initial_mesh_1d, 3};
+
+              {
+                auto initial_xj = MeshDataManager::instance().getMeshData(*initial_mesh_1d).xj();
+
+                for (CellId cell_id = 0; cell_id < initial_mesh_1d->numberOfCells(); ++cell_id) {
+                  initial_fh[cell_id] = 2 * dot(initial_xj[cell_id], initial_xj[cell_id]);
+
+                  initial_vh[cell_id][0] = 3 * initial_xj[cell_id][0] - 1;
+                  initial_vh[cell_id][1] = std::sin(initial_xj[cell_id][0]);
+                  initial_vh[cell_id][2] = -2 * initial_xj[cell_id][0] + 4;
+                }
+              }
+
+              LoadBalancer load_balancer;
+              auto balanced_discrete_function_list =
+                load_balancer.balance({std::make_shared<DiscreteFunctionVariant>(initial_fh),
+                                       std::make_shared<DiscreteFunctionVariant>(initial_vh)});
+
+              REQUIRE(balanced_discrete_function_list.size() == 2);
+
+              DiscreteFunctionP0<const double> final_fh =
+                balanced_discrete_function_list[0]->get<DiscreteFunctionP0<const double>>();
+              DiscreteFunctionP0Vector<const double> final_vh =
+                balanced_discrete_function_list[1]->get<DiscreteFunctionP0Vector<const double>>();
+
+              REQUIRE(final_vh.size() == 3);
+              {
+                auto final_mesh_1d_v = final_fh.meshVariant();
+                auto final_mesh_1d   = final_mesh_1d_v->get<Mesh<1>>();
+
+                REQUIRE(final_mesh_1d->id() > initial_mesh_1d->id());
+
+                REQUIRE_THROWS_WITH(load_balancer.balance({std::make_shared<DiscreteFunctionVariant>(initial_fh),
+                                                           std::make_shared<DiscreteFunctionVariant>(final_fh)}),
+                                    "error: discrete functions are not defined on the same mesh");
+
+                {
+                  auto final_xj = MeshDataManager::instance().getMeshData(*final_mesh_1d).xj();
+
+                  bool final_fh_is_same = true;
+                  bool final_vh_is_same = true;
+                  for (CellId cell_id = 0; cell_id < final_mesh_1d->numberOfCells(); ++cell_id) {
+                    if (final_fh[cell_id] != 2 * dot(final_xj[cell_id], final_xj[cell_id])) {
+                      final_fh_is_same = false;
+                    }
+
+                    if (final_vh[cell_id][0] != 3 * final_xj[cell_id][0] - 1) {
+                      final_vh_is_same = false;
+                    }
+                    if (final_vh[cell_id][1] != std::sin(final_xj[cell_id][0])) {
+                      final_vh_is_same = false;
+                    }
+                    if (final_vh[cell_id][2] != -2 * final_xj[cell_id][0] + 4) {
+                      final_vh_is_same = false;
+                    }
+                  }
+
+                  REQUIRE(parallel::allReduceAnd(final_fh_is_same));
+                  REQUIRE(parallel::allReduceAnd(final_vh_is_same));
+                }
+              }
+            }
+          }
+        }
+
+        SECTION("2D")
+        {
+          for (auto&& named_mesh : MeshDataBaseForTests::get().all2DMeshes()) {
+            std::shared_ptr<const Mesh<2>> initial_mesh_2d = named_mesh.mesh()->get<Mesh<2>>();
+            SECTION(named_mesh.name())
+            {
+              DiscreteFunctionP0<double> initial_fh{initial_mesh_2d};
+              DiscreteFunctionP0Vector<double> initial_vh{initial_mesh_2d, 3};
+
+              {
+                auto initial_xj = MeshDataManager::instance().getMeshData(*initial_mesh_2d).xj();
+
+                for (CellId cell_id = 0; cell_id < initial_mesh_2d->numberOfCells(); ++cell_id) {
+                  initial_fh[cell_id] = 2 * dot(initial_xj[cell_id], initial_xj[cell_id]);
+
+                  initial_vh[cell_id][0] = 3 * initial_xj[cell_id][0] - 1;
+                  initial_vh[cell_id][1] = std::sin(initial_xj[cell_id][0] + initial_xj[cell_id][1]);
+                  initial_vh[cell_id][2] = -2 * initial_xj[cell_id][0] * initial_xj[cell_id][1] + 4;
+                }
+              }
+
+              LoadBalancer load_balancer;
+              auto balanced_discrete_function_list =
+                load_balancer.balance({std::make_shared<DiscreteFunctionVariant>(initial_fh),
+                                       std::make_shared<DiscreteFunctionVariant>(initial_vh)});
+
+              REQUIRE(balanced_discrete_function_list.size() == 2);
+
+              DiscreteFunctionP0<const double> final_fh =
+                balanced_discrete_function_list[0]->get<DiscreteFunctionP0<const double>>();
+              DiscreteFunctionP0Vector<const double> final_vh =
+                balanced_discrete_function_list[1]->get<DiscreteFunctionP0Vector<const double>>();
+
+              REQUIRE(final_vh.size() == 3);
+              {
+                auto final_mesh_2d_v = final_fh.meshVariant();
+                auto final_mesh_2d   = final_mesh_2d_v->get<Mesh<2>>();
+
+                REQUIRE(final_mesh_2d->id() > initial_mesh_2d->id());
+
+                REQUIRE_THROWS_WITH(load_balancer.balance({std::make_shared<DiscreteFunctionVariant>(initial_fh),
+                                                           std::make_shared<DiscreteFunctionVariant>(final_fh)}),
+                                    "error: discrete functions are not defined on the same mesh");
+
+                {
+                  auto final_xj = MeshDataManager::instance().getMeshData(*final_mesh_2d).xj();
+
+                  bool final_fh_is_same = true;
+                  bool final_vh_is_same = true;
+                  for (CellId cell_id = 0; cell_id < final_mesh_2d->numberOfCells(); ++cell_id) {
+                    if (final_fh[cell_id] != 2 * dot(final_xj[cell_id], final_xj[cell_id])) {
+                      final_fh_is_same = false;
+                    }
+
+                    if (final_vh[cell_id][0] != 3 * final_xj[cell_id][0] - 1) {
+                      final_vh_is_same = false;
+                    }
+                    if (final_vh[cell_id][1] != std::sin(final_xj[cell_id][0] + final_xj[cell_id][1])) {
+                      final_vh_is_same = false;
+                    }
+                    if (final_vh[cell_id][2] != -2 * final_xj[cell_id][0] * final_xj[cell_id][1] + 4) {
+                      final_vh_is_same = false;
+                    }
+                  }
+
+                  REQUIRE(parallel::allReduceAnd(final_fh_is_same));
+                  REQUIRE(parallel::allReduceAnd(final_vh_is_same));
+                }
+              }
+            }
+          }
+        }
+
+        SECTION("3D")
+        {
+          for (auto&& named_mesh : MeshDataBaseForTests::get().all3DMeshes()) {
+            std::shared_ptr<const Mesh<3>> initial_mesh_3d = named_mesh.mesh()->get<Mesh<3>>();
+            SECTION(named_mesh.name())
+            {
+              DiscreteFunctionP0<double> initial_fh{initial_mesh_3d};
+              DiscreteFunctionP0Vector<double> initial_vh{initial_mesh_3d, 3};
+
+              {
+                auto initial_xj = MeshDataManager::instance().getMeshData(*initial_mesh_3d).xj();
+
+                for (CellId cell_id = 0; cell_id < initial_mesh_3d->numberOfCells(); ++cell_id) {
+                  initial_fh[cell_id] = 2 * dot(initial_xj[cell_id], initial_xj[cell_id]);
+
+                  initial_vh[cell_id][0] = 3 * initial_xj[cell_id][0] - 1;
+                  initial_vh[cell_id][1] = std::sin(initial_xj[cell_id][1] + initial_xj[cell_id][2]);
+                  initial_vh[cell_id][2] = -2 * initial_xj[cell_id][0] * initial_xj[cell_id][1] + 4;
+                }
+              }
+
+              LoadBalancer load_balancer;
+              auto balanced_discrete_function_list =
+                load_balancer.balance({std::make_shared<DiscreteFunctionVariant>(initial_fh),
+                                       std::make_shared<DiscreteFunctionVariant>(initial_vh)});
+
+              REQUIRE(balanced_discrete_function_list.size() == 2);
+
+              DiscreteFunctionP0<const double> final_fh =
+                balanced_discrete_function_list[0]->get<DiscreteFunctionP0<const double>>();
+              DiscreteFunctionP0Vector<const double> final_vh =
+                balanced_discrete_function_list[1]->get<DiscreteFunctionP0Vector<const double>>();
+
+              REQUIRE(final_vh.size() == 3);
+              {
+                auto final_mesh_3d_v = final_fh.meshVariant();
+                auto final_mesh_3d   = final_mesh_3d_v->get<Mesh<3>>();
+
+                REQUIRE(final_mesh_3d->id() > initial_mesh_3d->id());
+
+                REQUIRE_THROWS_WITH(load_balancer.balance({std::make_shared<DiscreteFunctionVariant>(initial_fh),
+                                                           std::make_shared<DiscreteFunctionVariant>(final_fh)}),
+                                    "error: discrete functions are not defined on the same mesh");
+
+                {
+                  auto final_xj = MeshDataManager::instance().getMeshData(*final_mesh_3d).xj();
+
+                  bool final_fh_is_same = true;
+                  bool final_vh_is_same = true;
+                  for (CellId cell_id = 0; cell_id < final_mesh_3d->numberOfCells(); ++cell_id) {
+                    if (final_fh[cell_id] != 2 * dot(final_xj[cell_id], final_xj[cell_id])) {
+                      final_fh_is_same = false;
+                    }
+
+                    if (final_vh[cell_id][0] != 3 * final_xj[cell_id][0] - 1) {
+                      final_vh_is_same = false;
+                    }
+                    if (final_vh[cell_id][1] != std::sin(final_xj[cell_id][1] + final_xj[cell_id][2])) {
+                      final_vh_is_same = false;
+                    }
+                    if (final_vh[cell_id][2] != -2 * final_xj[cell_id][0] * final_xj[cell_id][1] + 4) {
+                      final_vh_is_same = false;
+                    }
+                  }
+
+                  REQUIRE(parallel::allReduceAnd(final_fh_is_same));
+                  REQUIRE(parallel::allReduceAnd(final_vh_is_same));
+                }
+              }
+            }
+          }
+        }
+      }
+    }
+  }
+
+  PartitionerOptions::default_options.library() = initial_library;
+}

tests/test_LoadBalancer_serial.cpp

+#include <catch2/catch_all.hpp>
+#include <catch2/catch_test_macros.hpp>
+
+#include <scheme/LoadBalancer.hpp>
+
+#include <MeshDataBaseForTests.hpp>
+#include <mesh/Mesh.hpp>
+#include <mesh/MeshVariant.hpp>
+#include <scheme/DiscreteFunctionP0.hpp>
+#include <scheme/DiscreteFunctionP0Vector.hpp>
+#include <scheme/DiscreteFunctionVariant.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("LoadBalancer serial", "[scheme]")
+{
+  SECTION("1D")
+  {
+    for (auto&& named_mesh : MeshDataBaseForTests::get().all1DMeshes()) {
+      std::shared_ptr<const Mesh<1>> initial_mesh_1d = named_mesh.mesh()->get<Mesh<1>>();
+      SECTION(named_mesh.name())
+      {
+        DiscreteFunctionP0<double> initial_fh{initial_mesh_1d};
+        DiscreteFunctionP0Vector<double> initial_vh{initial_mesh_1d, 3};
+
+        {
+          auto initial_xj = MeshDataManager::instance().getMeshData(*initial_mesh_1d).xj();
+
+          for (CellId cell_id = 0; cell_id < initial_mesh_1d->numberOfCells(); ++cell_id) {
+            initial_fh[cell_id] = 2 * dot(initial_xj[cell_id], initial_xj[cell_id]);
+
+            initial_vh[cell_id][0] = 3 * initial_xj[cell_id][0] - 1;
+            initial_vh[cell_id][1] = std::sin(initial_xj[cell_id][0]);
+            initial_vh[cell_id][2] = -2 * initial_xj[cell_id][0] + 4;
+          }
+        }
+
+        LoadBalancer load_balancer;
+        auto balanced_discrete_function_list =
+          load_balancer.balance({std::make_shared<DiscreteFunctionVariant>(initial_fh),
+                                 std::make_shared<DiscreteFunctionVariant>(initial_vh)});
+
+        REQUIRE(balanced_discrete_function_list.size() == 2);
+
+        DiscreteFunctionP0<const double> final_fh =
+          balanced_discrete_function_list[0]->get<DiscreteFunctionP0<const double>>();
+        DiscreteFunctionP0Vector<const double> final_vh =
+          balanced_discrete_function_list[1]->get<DiscreteFunctionP0Vector<const double>>();
+
+        REQUIRE(final_vh.size() == 3);
+        {
+          auto final_mesh_1d_v = final_fh.meshVariant();
+          auto final_mesh_1d   = final_mesh_1d_v->get<Mesh<1>>();
+
+          REQUIRE(final_mesh_1d->id() > initial_mesh_1d->id());
+
+          REQUIRE_THROWS_WITH(load_balancer.balance({std::make_shared<DiscreteFunctionVariant>(initial_fh),
+                                                     std::make_shared<DiscreteFunctionVariant>(final_fh)}),
+                              "error: discrete functions are not defined on the same mesh");
+
+          {
+            auto final_xj = MeshDataManager::instance().getMeshData(*final_mesh_1d).xj();
+
+            bool final_fh_is_same = true;
+            bool final_vh_is_same = true;
+            for (CellId cell_id = 0; cell_id < final_mesh_1d->numberOfCells(); ++cell_id) {
+              if (final_fh[cell_id] != 2 * dot(final_xj[cell_id], final_xj[cell_id])) {
+                final_fh_is_same = false;
+              }
+
+              if (final_vh[cell_id][0] != 3 * final_xj[cell_id][0] - 1) {
+                final_vh_is_same = false;
+              }
+              if (final_vh[cell_id][1] != std::sin(final_xj[cell_id][0])) {
+                final_vh_is_same = false;
+              }
+              if (final_vh[cell_id][2] != -2 * final_xj[cell_id][0] + 4) {
+                final_vh_is_same = false;
+              }
+            }
+
+            REQUIRE(final_fh_is_same);
+            REQUIRE(final_vh_is_same);
+          }
+        }
+      }
+    }
+  }
+
+  SECTION("2D")
+  {
+    for (auto&& named_mesh : MeshDataBaseForTests::get().all2DMeshes()) {
+      std::shared_ptr<const Mesh<2>> initial_mesh_2d = named_mesh.mesh()->get<Mesh<2>>();
+      SECTION(named_mesh.name())
+      {
+        DiscreteFunctionP0<double> initial_fh{initial_mesh_2d};
+        DiscreteFunctionP0Vector<double> initial_vh{initial_mesh_2d, 3};
+
+        {
+          auto initial_xj = MeshDataManager::instance().getMeshData(*initial_mesh_2d).xj();
+
+          for (CellId cell_id = 0; cell_id < initial_mesh_2d->numberOfCells(); ++cell_id) {
+            initial_fh[cell_id] = 2 * dot(initial_xj[cell_id], initial_xj[cell_id]);
+
+            initial_vh[cell_id][0] = 3 * initial_xj[cell_id][0] - 1;
+            initial_vh[cell_id][1] = std::sin(initial_xj[cell_id][0] + initial_xj[cell_id][1]);
+            initial_vh[cell_id][2] = -2 * initial_xj[cell_id][0] * initial_xj[cell_id][1] + 4;
+          }
+        }
+
+        LoadBalancer load_balancer;
+        auto balanced_discrete_function_list =
+          load_balancer.balance({std::make_shared<DiscreteFunctionVariant>(initial_fh),
+                                 std::make_shared<DiscreteFunctionVariant>(initial_vh)});
+
+        REQUIRE(balanced_discrete_function_list.size() == 2);
+
+        DiscreteFunctionP0<const double> final_fh =
+          balanced_discrete_function_list[0]->get<DiscreteFunctionP0<const double>>();
+        DiscreteFunctionP0Vector<const double> final_vh =
+          balanced_discrete_function_list[1]->get<DiscreteFunctionP0Vector<const double>>();
+
+        REQUIRE(final_vh.size() == 3);
+        {
+          auto final_mesh_2d_v = final_fh.meshVariant();
+          auto final_mesh_2d   = final_mesh_2d_v->get<Mesh<2>>();
+
+          REQUIRE(final_mesh_2d->id() > initial_mesh_2d->id());
+
+          REQUIRE_THROWS_WITH(load_balancer.balance({std::make_shared<DiscreteFunctionVariant>(initial_fh),
+                                                     std::make_shared<DiscreteFunctionVariant>(final_fh)}),
+                              "error: discrete functions are not defined on the same mesh");
+
+          {
+            auto final_xj = MeshDataManager::instance().getMeshData(*final_mesh_2d).xj();
+
+            bool final_fh_is_same = true;
+            bool final_vh_is_same = true;
+            for (CellId cell_id = 0; cell_id < final_mesh_2d->numberOfCells(); ++cell_id) {
+              if (final_fh[cell_id] != 2 * dot(final_xj[cell_id], final_xj[cell_id])) {
+                final_fh_is_same = false;
+              }
+
+              if (final_vh[cell_id][0] != 3 * final_xj[cell_id][0] - 1) {
+                final_vh_is_same = false;
+              }
+              if (final_vh[cell_id][1] != std::sin(final_xj[cell_id][0] + final_xj[cell_id][1])) {
+                final_vh_is_same = false;
+              }
+              if (final_vh[cell_id][2] != -2 * final_xj[cell_id][0] * final_xj[cell_id][1] + 4) {
+                final_vh_is_same = false;
+              }
+            }
+
+            REQUIRE(final_fh_is_same);
+            REQUIRE(final_vh_is_same);
+          }
+        }
+      }
+    }
+  }
+
+  SECTION("3D")
+  {
+    for (auto&& named_mesh : MeshDataBaseForTests::get().all3DMeshes()) {
+      std::shared_ptr<const Mesh<3>> initial_mesh_3d = named_mesh.mesh()->get<Mesh<3>>();
+      SECTION(named_mesh.name())
+      {
+        DiscreteFunctionP0<double> initial_fh{initial_mesh_3d};
+        DiscreteFunctionP0Vector<double> initial_vh{initial_mesh_3d, 3};
+
+        {
+          auto initial_xj = MeshDataManager::instance().getMeshData(*initial_mesh_3d).xj();
+
+          for (CellId cell_id = 0; cell_id < initial_mesh_3d->numberOfCells(); ++cell_id) {
+            initial_fh[cell_id] = 2 * dot(initial_xj[cell_id], initial_xj[cell_id]);
+
+            initial_vh[cell_id][0] = 3 * initial_xj[cell_id][0] - 1;
+            initial_vh[cell_id][1] = std::sin(initial_xj[cell_id][1] + initial_xj[cell_id][2]);
+            initial_vh[cell_id][2] = -2 * initial_xj[cell_id][0] * initial_xj[cell_id][1] + 4;
+          }
+        }
+
+        LoadBalancer load_balancer;
+        auto balanced_discrete_function_list =
+          load_balancer.balance({std::make_shared<DiscreteFunctionVariant>(initial_fh),
+                                 std::make_shared<DiscreteFunctionVariant>(initial_vh)});
+
+        REQUIRE(balanced_discrete_function_list.size() == 2);
+
+        DiscreteFunctionP0<const double> final_fh =
+          balanced_discrete_function_list[0]->get<DiscreteFunctionP0<const double>>();
+        DiscreteFunctionP0Vector<const double> final_vh =
+          balanced_discrete_function_list[1]->get<DiscreteFunctionP0Vector<const double>>();
+
+        REQUIRE(final_vh.size() == 3);
+        {
+          auto final_mesh_3d_v = final_fh.meshVariant();
+          auto final_mesh_3d   = final_mesh_3d_v->get<Mesh<3>>();
+
+          REQUIRE(final_mesh_3d->id() > initial_mesh_3d->id());
+
+          REQUIRE_THROWS_WITH(load_balancer.balance({std::make_shared<DiscreteFunctionVariant>(initial_fh),
+                                                     std::make_shared<DiscreteFunctionVariant>(final_fh)}),
+                              "error: discrete functions are not defined on the same mesh");
+
+          {
+            auto final_xj = MeshDataManager::instance().getMeshData(*final_mesh_3d).xj();
+
+            bool final_fh_is_same = true;
+            bool final_vh_is_same = true;
+            for (CellId cell_id = 0; cell_id < final_mesh_3d->numberOfCells(); ++cell_id) {
+              if (final_fh[cell_id] != 2 * dot(final_xj[cell_id], final_xj[cell_id])) {
+                final_fh_is_same = false;
+              }
+
+              if (final_vh[cell_id][0] != 3 * final_xj[cell_id][0] - 1) {
+                final_vh_is_same = false;
+              }
+              if (final_vh[cell_id][1] != std::sin(final_xj[cell_id][1] + final_xj[cell_id][2])) {
+                final_vh_is_same = false;
+              }
+              if (final_vh[cell_id][2] != -2 * final_xj[cell_id][0] * final_xj[cell_id][1] + 4) {
+                final_vh_is_same = false;
+              }
+            }
+
+            REQUIRE(final_fh_is_same);
+            REQUIRE(final_vh_is_same);
+          }
+        }
+      }
+    }
+  }
+}