test_DiamondDualMeshBuilder.cpp

#include <catch2/catch_test_macros.hpp>
#include <catch2/matchers/catch_matchers_all.hpp>

#include <MeshDataBaseForTests.hpp>
#include <mesh/DualMeshManager.hpp>
#include <mesh/MeshData.hpp>
#include <mesh/MeshDataManager.hpp>

#include <mesh/Connectivity.hpp>
#include <mesh/ItemValueUtils.hpp>
#include <mesh/Mesh.hpp>

// clazy:excludeall=non-pod-global-static

TEST_CASE("DiamondDualMeshBuilder", "[mesh]")
{
  SECTION("2D")
  {
    constexpr static size_t Dimension = 2;

    using ConnectivityType = Connectivity<Dimension>;
    using MeshType         = Mesh<ConnectivityType>;

    std::shared_ptr p_mesh      = MeshDataBaseForTests::get().hybrid2DMesh();
    const MeshType& primal_mesh = *p_mesh;

    REQUIRE(primal_mesh.numberOfNodes() == 53);
    REQUIRE(primal_mesh.numberOfFaces() == 110);
    REQUIRE(primal_mesh.numberOfCells() == 58);

    std::shared_ptr p_diamond_dual_mesh = DualMeshManager::instance().getDiamondDualMesh(primal_mesh);
    const MeshType& dual_mesh           = *p_diamond_dual_mesh;

    REQUIRE(dual_mesh.numberOfNodes() == 111);
    REQUIRE(dual_mesh.numberOfFaces() == 220);
    REQUIRE(dual_mesh.numberOfCells() == 110);

    auto dual_cell_volume = MeshDataManager::instance().getMeshData(dual_mesh).Vj();

    REQUIRE(sum(dual_cell_volume) == Catch::Approx(2));

    auto coord_comp = [](const TinyVector<Dimension>& a, const TinyVector<Dimension>& b) -> bool {
      return (a[0] < b[0]) or ((a[0] == b[0]) and (a[1] < b[1]));
    };

    auto dual_cell_to_node   = dual_mesh.connectivity().cellToNodeMatrix();
    auto primal_face_to_node = primal_mesh.connectivity().faceToNodeMatrix();
    auto primal_face_to_cell = primal_mesh.connectivity().faceToCellMatrix();

    auto dual_xr   = dual_mesh.xr();
    auto primal_xr = primal_mesh.xr();
    auto primal_xj = MeshDataManager::instance().getMeshData(primal_mesh).xj();

    using CoordSet = std::set<TinyVector<Dimension>,
                              std::function<bool(const TinyVector<Dimension>& a, const TinyVector<Dimension>& b)>>;

    CellId dual_cell_id   = 0;
    FaceId primal_face_id = 0;
    for (; dual_cell_id < dual_mesh.numberOfCells(); ++dual_cell_id, ++primal_face_id) {
      CoordSet dual_cell_node_coords{coord_comp};
      {
        auto dual_cell_node_list = dual_cell_to_node[dual_cell_id];
        for (size_t i_node = 0; i_node < dual_cell_node_list.size(); ++i_node) {
          dual_cell_node_coords.insert(dual_xr[dual_cell_node_list[i_node]]);
        }
      }

      CoordSet primal_coords{coord_comp};
      {
        auto primal_face_node_list = primal_face_to_node[primal_face_id];
        for (size_t i_node = 0; i_node < primal_face_node_list.size(); ++i_node) {
          primal_coords.insert(primal_xr[primal_face_node_list[i_node]]);
        }
        auto primal_face_cell_list = primal_face_to_cell[primal_face_id];
        for (size_t i_cell = 0; i_cell < primal_face_cell_list.size(); ++i_cell) {
          primal_coords.insert(primal_xj[primal_face_cell_list[i_cell]]);
        }
      }

      REQUIRE(primal_coords == dual_cell_node_coords);
    }
  }

  SECTION("3D")
  {
    constexpr static size_t Dimension = 3;

    using ConnectivityType = Connectivity<Dimension>;
    using MeshType         = Mesh<ConnectivityType>;

    std::shared_ptr p_mesh      = MeshDataBaseForTests::get().hybrid3DMesh();
    const MeshType& primal_mesh = *p_mesh;

    REQUIRE(primal_mesh.numberOfNodes() == 132);
    REQUIRE(primal_mesh.numberOfEdges() == 452);
    REQUIRE(primal_mesh.numberOfFaces() == 520);
    REQUIRE(primal_mesh.numberOfCells() == 199);

    std::shared_ptr p_diamond_dual_mesh = DualMeshManager::instance().getDiamondDualMesh(primal_mesh);
    const MeshType& dual_mesh           = *p_diamond_dual_mesh;

    REQUIRE(dual_mesh.numberOfNodes() == 331);
    REQUIRE(dual_mesh.numberOfEdges() == 1461);
    REQUIRE(dual_mesh.numberOfFaces() == 1651);
    REQUIRE(dual_mesh.numberOfCells() == 520);

    auto dual_cell_volume = MeshDataManager::instance().getMeshData(dual_mesh).Vj();

    REQUIRE(sum(dual_cell_volume) == Catch::Approx(2));

    auto coord_comp = [](const TinyVector<Dimension>& a, const TinyVector<Dimension>& b) -> bool {
      return (a[0] < b[0]) or ((a[0] == b[0]) and (a[1] < b[1])) or
             ((a[0] == b[0]) and (a[1] == b[1]) and (a[2] < b[2]));
    };

    auto dual_cell_to_node   = dual_mesh.connectivity().cellToNodeMatrix();
    auto primal_face_to_node = primal_mesh.connectivity().faceToNodeMatrix();
    auto primal_face_to_cell = primal_mesh.connectivity().faceToCellMatrix();

    auto dual_xr   = dual_mesh.xr();
    auto primal_xr = primal_mesh.xr();
    auto primal_xj = MeshDataManager::instance().getMeshData(primal_mesh).xj();

    using CoordSet = std::set<TinyVector<Dimension>,
                              std::function<bool(const TinyVector<Dimension>& a, const TinyVector<Dimension>& b)>>;

    CellId dual_cell_id   = 0;
    FaceId primal_face_id = 0;
    for (; dual_cell_id < dual_mesh.numberOfCells(); ++dual_cell_id, ++primal_face_id) {
      CoordSet dual_cell_node_coords{coord_comp};
      {
        auto dual_cell_node_list = dual_cell_to_node[dual_cell_id];
        for (size_t i_node = 0; i_node < dual_cell_node_list.size(); ++i_node) {
          dual_cell_node_coords.insert(dual_xr[dual_cell_node_list[i_node]]);
        }
      }

      CoordSet primal_coords{coord_comp};
      {
        auto primal_face_node_list = primal_face_to_node[primal_face_id];
        for (size_t i_node = 0; i_node < primal_face_node_list.size(); ++i_node) {
          primal_coords.insert(primal_xr[primal_face_node_list[i_node]]);
        }
        auto primal_face_cell_list = primal_face_to_cell[primal_face_id];
        for (size_t i_cell = 0; i_cell < primal_face_cell_list.size(); ++i_cell) {
          primal_coords.insert(primal_xj[primal_face_cell_list[i_cell]]);
        }
      }

      REQUIRE(primal_coords == dual_cell_node_coords);
    }
  }
}