From 95f2c6f0dc79133359da0d662735d6aa5ec13bc4 Mon Sep 17 00:00:00 2001
From: Stephane Del Pino <stephane.delpino44@gmail.com>
Date: Tue, 2 Jun 2020 16:22:45 +0200
Subject: [PATCH] Add regular Cartesian mesh generation
The mesh is defined by two opposite corners and the number of cells in
each direction.
---
CMakeLists.txt | 1 +
src/language/modules/MeshModule.cpp | 87 ++++-
src/mesh/CMakeLists.txt | 1 +
src/mesh/CartesianMeshBuilder.cpp | 549 ++++++++++++++++++++++++++++
src/mesh/CartesianMeshBuilder.hpp | 32 ++
5 files changed, 663 insertions(+), 7 deletions(-)
create mode 100644 src/mesh/CartesianMeshBuilder.cpp
create mode 100644 src/mesh/CartesianMeshBuilder.hpp
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 168fe3a82..33f0df0c8 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -397,6 +397,7 @@ target_link_libraries(
PugsLanguage
PugsLanguageAST
PugsLanguageModules
+ PugsMesh
PugsLanguageUtils
kokkos
${PARMETIS_LIBRARIES}
diff --git a/src/language/modules/MeshModule.cpp b/src/language/modules/MeshModule.cpp
index f4472b749..417737777 100644
--- a/src/language/modules/MeshModule.cpp
+++ b/src/language/modules/MeshModule.cpp
@@ -1,11 +1,13 @@
#include <language/modules/MeshModule.hpp>
+#include <algebra/TinyVector.hpp>
#include <language/node_processor/ExecutionPolicy.hpp>
#include <language/utils/BuiltinFunctionEmbedder.hpp>
#include <language/utils/FunctionTable.hpp>
#include <language/utils/PugsFunctionAdapter.hpp>
#include <language/utils/SymbolTable.hpp>
#include <language/utils/TypeDescriptor.hpp>
+#include <mesh/CartesianMeshBuilder.hpp>
#include <mesh/Connectivity.hpp>
#include <mesh/GmshReader.hpp>
#include <mesh/Mesh.hpp>
@@ -104,6 +106,66 @@ MeshModule::MeshModule()
));
+ this->_addBuiltinFunction("cartesian1dMesh",
+ std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<const IMesh>, TinyVector<1>,
+ TinyVector<1>, std::vector<uint64_t>>>(
+
+ std::function<std::shared_ptr<const IMesh>(TinyVector<1>, TinyVector<1>,
+ std::vector<uint64_t>)>{
+
+ [&](const TinyVector<1> a, const TinyVector<1> b,
+ const std::vector<uint64_t>& box_sizes) -> std::shared_ptr<const IMesh> {
+ constexpr uint64_t dimension = 1;
+
+ if (box_sizes.size() != dimension) {
+ throw NormalError("expecting " + std::to_string(dimension) +
+ " dimensions, provided " + std::to_string(box_sizes.size()));
+ }
+
+ const TinyVector<dimension, uint64_t> sizes = [&]() {
+ TinyVector<dimension, uint64_t> s;
+ for (size_t i = 0; i < dimension; ++i) {
+ s[i] = box_sizes[i];
+ }
+ return s;
+ }();
+
+ CartesianMeshBuilder builder{a, b, sizes};
+ return builder.mesh();
+ }}
+
+ ));
+
+ this->_addBuiltinFunction("cartesian2dMesh",
+ std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<const IMesh>, TinyVector<2>,
+ TinyVector<2>, std::vector<uint64_t>>>(
+
+ std::function<std::shared_ptr<const IMesh>(TinyVector<2>, TinyVector<2>,
+ std::vector<uint64_t>)>{
+
+ [&](const TinyVector<2> a, const TinyVector<2> b,
+ const std::vector<uint64_t>& box_sizes) -> std::shared_ptr<const IMesh> {
+ constexpr uint64_t dimension = 2;
+
+ if (box_sizes.size() != dimension) {
+ throw NormalError("expecting " + std::to_string(dimension) +
+ " dimensions, provided " + std::to_string(box_sizes.size()));
+ }
+
+ const TinyVector<dimension, uint64_t> sizes = [&]() {
+ TinyVector<dimension, uint64_t> s;
+ for (size_t i = 0; i < dimension; ++i) {
+ s[i] = box_sizes[i];
+ }
+ return s;
+ }();
+
+ CartesianMeshBuilder builder{a, b, sizes};
+ return builder.mesh();
+ }}
+
+ ));
+
this->_addBuiltinFunction("cartesian3dMesh",
std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<const IMesh>, TinyVector<3>,
TinyVector<3>, std::vector<uint64_t>>>(
@@ -111,14 +173,25 @@ MeshModule::MeshModule()
std::function<std::shared_ptr<const IMesh>(TinyVector<3>, TinyVector<3>,
std::vector<uint64_t>)>{
- [](const TinyVector<3>& a, const TinyVector<3>& b,
- const std::vector<uint64_t>& box_sizes) -> std::shared_ptr<const IMesh> {
- std::cout << a << " -> " << b << " ";
- for (auto i : box_sizes) {
- std::cout << i << ' ';
+ [&](const TinyVector<3>& a, const TinyVector<3>& b,
+ const std::vector<uint64_t>& box_sizes) -> std::shared_ptr<const IMesh> {
+ constexpr uint64_t dimension = 3;
+
+ if (box_sizes.size() != dimension) {
+ throw NormalError("expecting " + std::to_string(dimension) +
+ " dimensions, provided " + std::to_string(box_sizes.size()));
}
- std::cout << '\n';
- return nullptr;
+
+ const TinyVector<dimension, uint64_t> sizes = [&]() {
+ TinyVector<dimension, uint64_t> s;
+ for (size_t i = 0; i < dimension; ++i) {
+ s[i] = box_sizes[i];
+ }
+ return s;
+ }();
+
+ CartesianMeshBuilder builder{a, b, sizes};
+ return builder.mesh();
}}
));
diff --git a/src/mesh/CMakeLists.txt b/src/mesh/CMakeLists.txt
index 5baed9b1b..3b92ca14e 100644
--- a/src/mesh/CMakeLists.txt
+++ b/src/mesh/CMakeLists.txt
@@ -2,6 +2,7 @@
add_library(
PugsMesh
+ CartesianMeshBuilder.cpp
Connectivity.cpp
ConnectivityComputer.cpp
ConnectivityDispatcher.cpp
diff --git a/src/mesh/CartesianMeshBuilder.cpp b/src/mesh/CartesianMeshBuilder.cpp
new file mode 100644
index 000000000..f422eb93f
--- /dev/null
+++ b/src/mesh/CartesianMeshBuilder.cpp
@@ -0,0 +1,549 @@
+#include <mesh/CartesianMeshBuilder.hpp>
+
+#include <mesh/Connectivity.hpp>
+#include <mesh/ConnectivityDescriptor.hpp>
+#include <mesh/ConnectivityDispatcher.hpp>
+#include <mesh/RefId.hpp>
+#include <utils/Array.hpp>
+#include <utils/Messenger.hpp>
+
+template <size_t Dimension>
+void
+CartesianMeshBuilder::_buildBoundaryNodeList(const TinyVector<Dimension, uint64_t>&, ConnectivityDescriptor&)
+{
+ static_assert(Dimension <= 1, "unexpected dimension");
+}
+
+template <>
+void
+CartesianMeshBuilder::_buildBoundaryNodeList(const TinyVector<1, uint64_t>& cell_size,
+ ConnectivityDescriptor& descriptor)
+{
+ { // xmin
+ Array<NodeId> boundary_faces(1);
+ boundary_faces[0] = 0;
+ descriptor.addRefItemList(RefNodeList{RefId{0, "XMIN"}, boundary_faces});
+ }
+ { // xmax
+ Array<NodeId> boundary_faces(1);
+ boundary_faces[0] = cell_size[0];
+ descriptor.addRefItemList(RefNodeList{RefId{1, "XMAX"}, boundary_faces});
+ }
+}
+
+template <size_t Dimension>
+void
+CartesianMeshBuilder::_buildBoundaryFaceList(const TinyVector<Dimension, uint64_t>&, ConnectivityDescriptor&)
+{
+ static_assert(Dimension >= 2 and Dimension <= 3, "unexpected dimension");
+}
+
+template <>
+void
+CartesianMeshBuilder::_buildBoundaryFaceList(const TinyVector<2, uint64_t>& cell_size,
+ ConnectivityDescriptor& descriptor)
+{
+ auto cell_number = [&](const TinyVector<2, uint64_t>& cell_logic_id) {
+ return cell_logic_id[0] * cell_size[1] + cell_logic_id[1];
+ };
+
+ { // xmin
+ const size_t i = 0;
+ Array<FaceId> boundary_faces(cell_size[1]);
+ for (size_t j = 0; j < cell_size[1]; ++j) {
+ constexpr size_t left_face = 3;
+
+ const size_t cell_id = cell_number(TinyVector<2, uint64_t>{i, j});
+ const size_t face_id = descriptor.cell_to_face_vector[cell_id][left_face];
+
+ boundary_faces[j] = face_id;
+ }
+ descriptor.addRefItemList(RefFaceList{RefId{0, "XMIN"}, boundary_faces});
+ }
+
+ { // xmax
+ const size_t i = cell_size[0] - 1;
+ Array<FaceId> boundary_faces(cell_size[1]);
+ for (size_t j = 0; j < cell_size[1]; ++j) {
+ constexpr size_t right_face = 1;
+
+ const size_t cell_id = cell_number(TinyVector<2, uint64_t>{i, j});
+ const size_t face_id = descriptor.cell_to_face_vector[cell_id][right_face];
+
+ boundary_faces[j] = face_id;
+ }
+ descriptor.addRefItemList(RefFaceList{RefId{1, "XMAX"}, boundary_faces});
+ }
+
+ { // ymin
+ const size_t j = 0;
+ Array<FaceId> boundary_faces(cell_size[0]);
+ for (size_t i = 0; i < cell_size[0]; ++i) {
+ constexpr size_t bottom_face = 0;
+
+ const size_t cell_id = cell_number(TinyVector<2, uint64_t>{i, j});
+ const size_t face_id = descriptor.cell_to_face_vector[cell_id][bottom_face];
+
+ boundary_faces[i] = face_id;
+ }
+ descriptor.addRefItemList(RefFaceList{RefId{2, "YMIN"}, boundary_faces});
+ }
+
+ { // ymax
+ const size_t j = cell_size[1] - 1;
+ Array<FaceId> boundary_faces(cell_size[0]);
+ for (size_t i = 0; i < cell_size[0]; ++i) {
+ constexpr size_t top_face = 2;
+
+ const size_t cell_id = cell_number(TinyVector<2, uint64_t>{i, j});
+ const size_t face_id = descriptor.cell_to_face_vector[cell_id][top_face];
+
+ boundary_faces[i] = face_id;
+ }
+ descriptor.addRefItemList(RefFaceList{RefId{3, "YMAX"}, boundary_faces});
+ }
+}
+
+template <>
+void
+CartesianMeshBuilder::_buildBoundaryFaceList(const TinyVector<3, uint64_t>& cell_size,
+ ConnectivityDescriptor& descriptor)
+{
+ using Face = ConnectivityFace<3>;
+
+ const std::unordered_map<Face, FaceId, typename Face::Hash> face_to_id_map = [&] {
+ std::unordered_map<Face, FaceId, typename Face::Hash> face_to_id_map;
+ for (FaceId l = 0; l < descriptor.face_to_node_vector.size(); ++l) {
+ const auto& node_vector = descriptor.face_to_node_vector[l];
+ face_to_id_map[Face(node_vector, descriptor.node_number_vector)] = l;
+ }
+ return face_to_id_map;
+ }();
+
+ const std::unordered_map<int, FaceId> face_number_id_map = [&] {
+ std::unordered_map<int, FaceId> face_number_id_map;
+ for (size_t l = 0; l < descriptor.face_number_vector.size(); ++l) {
+ face_number_id_map[descriptor.face_number_vector[l]] = l;
+ }
+ Assert(face_number_id_map.size() == descriptor.face_number_vector.size());
+ return face_number_id_map;
+ }();
+
+ const TinyVector<3, uint64_t> node_size{cell_size[0] + 1, cell_size[1] + 1, cell_size[2] + 1};
+ auto node_number = [&](const TinyVector<3, uint64_t>& node_logic_id) {
+ return (node_logic_id[0] * node_size[1] + node_logic_id[1]) * node_size[2] + node_logic_id[2];
+ };
+
+ { // xmin
+ const size_t i = 0;
+ Array<FaceId> boundary_faces(cell_size[1] * cell_size[2]);
+ size_t l = 0;
+ for (size_t j = 0; j < cell_size[1]; ++j) {
+ for (size_t k = 0; k < cell_size[2]; ++k) {
+ const uint32_t node_0_id = node_number(TinyVector<3, uint64_t>{i, j, k});
+ const uint32_t node_1_id = node_number(TinyVector<3, uint64_t>{i, j + 1, k});
+ const uint32_t node_2_id = node_number(TinyVector<3, uint64_t>{i, j + 1, k + 1});
+ const uint32_t node_3_id = node_number(TinyVector<3, uint64_t>{i, j, k + 1});
+
+ auto i_face =
+ face_to_id_map.find(Face{{node_0_id, node_1_id, node_2_id, node_3_id}, descriptor.node_number_vector});
+ Assert(i_face != face_to_id_map.end());
+
+ boundary_faces[l++] = i_face->second;
+ }
+ }
+ Assert(l == cell_size[1] * cell_size[2]);
+ descriptor.addRefItemList(RefFaceList{RefId{0, "XMIN"}, boundary_faces});
+ }
+
+ { // xmax
+ const size_t i = cell_size[0] - 1;
+ Array<FaceId> boundary_faces(cell_size[1] * cell_size[2]);
+ size_t l = 0;
+ for (size_t j = 0; j < cell_size[1]; ++j) {
+ for (size_t k = 0; k < cell_size[2]; ++k) {
+ const uint32_t node_0_id = node_number(TinyVector<3, uint64_t>{i + 1, j, k});
+ const uint32_t node_1_id = node_number(TinyVector<3, uint64_t>{i + 1, j + 1, k});
+ const uint32_t node_2_id = node_number(TinyVector<3, uint64_t>{i + 1, j + 1, k + 1});
+ const uint32_t node_3_id = node_number(TinyVector<3, uint64_t>{i + 1, j, k + 1});
+
+ auto i_face =
+ face_to_id_map.find(Face{{node_0_id, node_1_id, node_2_id, node_3_id}, descriptor.node_number_vector});
+ Assert(i_face != face_to_id_map.end());
+
+ boundary_faces[l++] = i_face->second;
+ }
+ }
+ Assert(l == cell_size[1] * cell_size[2]);
+ descriptor.addRefItemList(RefFaceList{RefId{1, "XMAX"}, boundary_faces});
+ }
+
+ { // ymin
+ const size_t j = 0;
+ Array<FaceId> boundary_faces(cell_size[0] * cell_size[2]);
+ size_t l = 0;
+ for (size_t i = 0; i < cell_size[0]; ++i) {
+ for (size_t k = 0; k < cell_size[2]; ++k) {
+ const uint32_t node_0_id = node_number(TinyVector<3, uint64_t>{i, j, k});
+ const uint32_t node_1_id = node_number(TinyVector<3, uint64_t>{i + 1, j, k});
+ const uint32_t node_2_id = node_number(TinyVector<3, uint64_t>{i + 1, j, k + 1});
+ const uint32_t node_3_id = node_number(TinyVector<3, uint64_t>{i, j, k + 1});
+
+ auto i_face =
+ face_to_id_map.find(Face{{node_0_id, node_1_id, node_2_id, node_3_id}, descriptor.node_number_vector});
+ Assert(i_face != face_to_id_map.end());
+
+ boundary_faces[l++] = i_face->second;
+ }
+ }
+ Assert(l == cell_size[0] * cell_size[2]);
+ descriptor.addRefItemList(RefFaceList{RefId{2, "YMIN"}, boundary_faces});
+ }
+
+ {
+ // ymax
+ const size_t j = cell_size[1] - 1;
+ Array<FaceId> boundary_faces(cell_size[0] * cell_size[2]);
+ size_t l = 0;
+ for (size_t i = 0; i < cell_size[0]; ++i) {
+ for (size_t k = 0; k < cell_size[2]; ++k) {
+ const uint32_t node_0_id = node_number(TinyVector<3, uint64_t>{i, j + 1, k});
+ const uint32_t node_1_id = node_number(TinyVector<3, uint64_t>{i + 1, j + 1, k});
+ const uint32_t node_2_id = node_number(TinyVector<3, uint64_t>{i + 1, j + 1, k + 1});
+ const uint32_t node_3_id = node_number(TinyVector<3, uint64_t>{i, j + 1, k + 1});
+
+ auto i_face =
+ face_to_id_map.find(Face{{node_0_id, node_1_id, node_2_id, node_3_id}, descriptor.node_number_vector});
+ Assert(i_face != face_to_id_map.end());
+
+ boundary_faces[l++] = i_face->second;
+ }
+ }
+ Assert(l == cell_size[0] * cell_size[2]);
+ descriptor.addRefItemList(RefFaceList{RefId{3, "YMAX"}, boundary_faces});
+ }
+
+ { // zmin
+ const size_t k = 0;
+ Array<FaceId> boundary_faces(cell_size[0] * cell_size[1]);
+ size_t l = 0;
+ for (size_t i = 0; i < cell_size[0]; ++i) {
+ for (size_t j = 0; j < cell_size[1]; ++j) {
+ const uint32_t node_0_id = node_number(TinyVector<3, uint64_t>{i, j, k});
+ const uint32_t node_1_id = node_number(TinyVector<3, uint64_t>{i + 1, j, k});
+ const uint32_t node_2_id = node_number(TinyVector<3, uint64_t>{i + 1, j + 1, k});
+ const uint32_t node_3_id = node_number(TinyVector<3, uint64_t>{i, j + 1, k});
+
+ auto i_face =
+ face_to_id_map.find(Face{{node_0_id, node_1_id, node_2_id, node_3_id}, descriptor.node_number_vector});
+ Assert(i_face != face_to_id_map.end());
+
+ boundary_faces[l++] = i_face->second;
+ }
+ }
+ Assert(l == cell_size[0] * cell_size[1]);
+ descriptor.addRefItemList(RefFaceList{RefId{4, "ZMIN"}, boundary_faces});
+ }
+
+ { // zmax
+ const size_t k = cell_size[2] - 1;
+ Array<FaceId> boundary_faces(cell_size[0] * cell_size[1]);
+ size_t l = 0;
+ for (size_t i = 0; i < cell_size[0]; ++i) {
+ for (size_t j = 0; j < cell_size[1]; ++j) {
+ const uint32_t node_0_id = node_number(TinyVector<3, uint64_t>{i, j, k + 1});
+ const uint32_t node_1_id = node_number(TinyVector<3, uint64_t>{i + 1, j, k + 1});
+ const uint32_t node_2_id = node_number(TinyVector<3, uint64_t>{i + 1, j + 1, k + 1});
+ const uint32_t node_3_id = node_number(TinyVector<3, uint64_t>{i, j + 1, k + 1});
+
+ auto i_face =
+ face_to_id_map.find(Face{{node_0_id, node_1_id, node_2_id, node_3_id}, descriptor.node_number_vector});
+ Assert(i_face != face_to_id_map.end());
+
+ boundary_faces[l++] = i_face->second;
+ }
+ }
+ Assert(l == cell_size[0] * cell_size[1]);
+ descriptor.addRefItemList(RefFaceList{RefId{5, "ZMAX"}, boundary_faces});
+ }
+}
+
+template <>
+void
+CartesianMeshBuilder::_buildCartesianMesh(const TinyVector<1>& a,
+ const TinyVector<1>& b,
+ const TinyVector<1, uint64_t>& cell_size)
+{
+ const size_t number_of_cells = cell_size[0];
+ const size_t number_of_nodes = cell_size[0] + 1;
+
+ ConnectivityDescriptor descriptor;
+ descriptor.node_number_vector.resize(number_of_nodes);
+ for (size_t i = 0; i < number_of_nodes; ++i) {
+ descriptor.node_number_vector[i] = i;
+ }
+
+ descriptor.cell_number_vector.resize(number_of_cells);
+ for (size_t i = 0; i < number_of_cells; ++i) {
+ descriptor.cell_number_vector[i] = i;
+ }
+
+ descriptor.cell_type_vector.resize(number_of_cells);
+ std::fill(descriptor.cell_type_vector.begin(), descriptor.cell_type_vector.end(), CellType::Line);
+
+ descriptor.cell_to_node_vector.resize(number_of_cells);
+ constexpr size_t nb_node_per_cell = 2;
+ for (size_t j = 0; j < number_of_cells; ++j) {
+ descriptor.cell_to_node_vector[j].resize(nb_node_per_cell);
+ for (size_t r = 0; r < nb_node_per_cell; ++r) {
+ descriptor.cell_to_node_vector[j][0] = j;
+ descriptor.cell_to_node_vector[j][1] = j + 1;
+ }
+ }
+
+ this->_buildBoundaryNodeList(cell_size, descriptor);
+
+ descriptor.cell_owner_vector.resize(number_of_cells);
+ std::fill(descriptor.cell_owner_vector.begin(), descriptor.cell_owner_vector.end(), parallel::rank());
+
+ descriptor.node_owner_vector.resize(descriptor.node_number_vector.size());
+ std::fill(descriptor.node_owner_vector.begin(), descriptor.node_owner_vector.end(), parallel::rank());
+
+ std::shared_ptr p_connectivity = Connectivity1D::build(descriptor);
+ Connectivity1D& connectivity = *p_connectivity;
+
+ const double h = (b[0] - a[0]) / cell_size[0];
+
+ NodeValue<TinyVector<1>> xr(connectivity);
+ for (NodeId r = 0; r < number_of_nodes; ++r) {
+ xr[r] = a + r * h;
+ }
+
+ m_mesh = std::make_shared<Mesh<Connectivity1D>>(p_connectivity, xr);
+}
+
+template <>
+void
+CartesianMeshBuilder::_buildCartesianMesh(const TinyVector<2>& a,
+ const TinyVector<2>& b,
+ const TinyVector<2, uint64_t>& cell_size)
+{
+ constexpr size_t Dimension = 2;
+
+ const TinyVector<Dimension, uint64_t> node_size{cell_size[0] + 1, cell_size[1] + 1};
+
+ const size_t number_of_cells = cell_size[0] * cell_size[1];
+ const size_t number_of_nodes = node_size[0] * node_size[1];
+
+ ConnectivityDescriptor descriptor;
+ descriptor.node_number_vector.resize(number_of_nodes);
+ for (size_t i = 0; i < number_of_nodes; ++i) {
+ descriptor.node_number_vector[i] = i;
+ }
+
+ descriptor.cell_number_vector.resize(number_of_cells);
+ for (size_t i = 0; i < number_of_cells; ++i) {
+ descriptor.cell_number_vector[i] = i;
+ }
+
+ descriptor.cell_type_vector.resize(number_of_cells);
+ std::fill(descriptor.cell_type_vector.begin(), descriptor.cell_type_vector.end(), CellType::Quadrangle);
+
+ auto node_logic_id = [&](size_t r) {
+ const uint64_t r0 = r / node_size[1];
+ const uint64_t r1 = r % node_size[1];
+ return TinyVector<Dimension, uint64_t>{r0, r1};
+ };
+
+ auto node_number = [&](const TinyVector<Dimension, uint64_t>& node_logic_id) {
+ return node_logic_id[0] * node_size[1] + node_logic_id[1];
+ };
+
+ auto cell_logic_id = [&](size_t j) {
+ const uint64_t j0 = j / cell_size[1];
+ const uint64_t j1 = j % cell_size[1];
+ return TinyVector<Dimension, uint64_t>{j0, j1};
+ };
+
+ descriptor.cell_to_node_vector.resize(number_of_cells);
+ constexpr size_t nb_node_per_cell = 1 << Dimension;
+ for (size_t j = 0; j < number_of_cells; ++j) {
+ TinyVector<Dimension, size_t> cell_index = cell_logic_id(j);
+ descriptor.cell_to_node_vector[j].resize(nb_node_per_cell);
+ for (size_t r = 0; r < nb_node_per_cell; ++r) {
+ descriptor.cell_to_node_vector[j][0] = node_number(cell_index + TinyVector<Dimension, uint64_t>{0, 0});
+ descriptor.cell_to_node_vector[j][1] = node_number(cell_index + TinyVector<Dimension, uint64_t>{1, 0});
+ descriptor.cell_to_node_vector[j][2] = node_number(cell_index + TinyVector<Dimension, uint64_t>{1, 1});
+ descriptor.cell_to_node_vector[j][3] = node_number(cell_index + TinyVector<Dimension, uint64_t>{0, 1});
+ }
+ }
+
+ MeshBuilderBase::_computeCellFaceAndFaceNodeConnectivities<Dimension>(descriptor);
+
+ this->_buildBoundaryFaceList(cell_size, descriptor);
+
+ descriptor.cell_owner_vector.resize(number_of_cells);
+ std::fill(descriptor.cell_owner_vector.begin(), descriptor.cell_owner_vector.end(), parallel::rank());
+
+ descriptor.face_owner_vector.resize(descriptor.face_number_vector.size());
+ std::fill(descriptor.face_owner_vector.begin(), descriptor.face_owner_vector.end(), parallel::rank());
+
+ descriptor.node_owner_vector.resize(descriptor.node_number_vector.size());
+ std::fill(descriptor.node_owner_vector.begin(), descriptor.node_owner_vector.end(), parallel::rank());
+
+ std::shared_ptr p_connectivity = Connectivity<Dimension>::build(descriptor);
+ Connectivity<Dimension>& connectivity = *p_connectivity;
+
+ const TinyVector<Dimension> h{(b[0] - a[0]) / cell_size[0], (b[1] - a[1]) / cell_size[1]};
+
+ NodeValue<TinyVector<Dimension>> xr(connectivity);
+ for (NodeId r = 0; r < number_of_nodes; ++r) {
+ const TinyVector<Dimension, uint64_t> node_index = node_logic_id(r);
+ for (size_t i = 0; i < Dimension; ++i) {
+ xr[r][i] = a[i] + node_index[i] * h[i];
+ }
+ }
+
+ m_mesh = std::make_shared<Mesh<Connectivity<Dimension>>>(p_connectivity, xr);
+}
+
+template <>
+void
+CartesianMeshBuilder::_buildCartesianMesh(const TinyVector<3>& a,
+ const TinyVector<3>& b,
+ const TinyVector<3, uint64_t>& cell_size)
+{
+ constexpr size_t Dimension = 3;
+
+ const TinyVector<Dimension, uint64_t> node_size = [&] {
+ TinyVector node_size{cell_size};
+ for (size_t i = 0; i < Dimension; ++i) {
+ node_size[i] += 1;
+ }
+ return node_size;
+ }();
+
+ auto count_items = [](auto&& v) {
+ using V_Type = std::decay_t<decltype(v)>;
+ static_assert(is_tiny_vector_v<V_Type>);
+
+ size_t sum = v[0];
+ for (size_t i = 1; i < Dimension; ++i) {
+ sum *= v[i];
+ }
+ return sum;
+ };
+
+ const size_t number_of_cells = count_items(cell_size);
+ const size_t number_of_nodes = count_items(node_size);
+
+ ConnectivityDescriptor descriptor;
+ descriptor.node_number_vector.resize(number_of_nodes);
+ for (size_t i = 0; i < number_of_nodes; ++i) {
+ descriptor.node_number_vector[i] = i;
+ }
+
+ descriptor.cell_number_vector.resize(number_of_cells);
+ for (size_t i = 0; i < number_of_cells; ++i) {
+ descriptor.cell_number_vector[i] = i;
+ }
+
+ descriptor.cell_type_vector.resize(number_of_cells);
+ std::fill(descriptor.cell_type_vector.begin(), descriptor.cell_type_vector.end(), CellType::Hexahedron);
+
+ auto cell_logic_id = [&](size_t j) {
+ const size_t slice1 = cell_size[1] * cell_size[2];
+ const size_t& slice2 = cell_size[2];
+ const uint64_t j0 = j / slice1;
+ const uint64_t j1 = (j - j0 * slice1) / slice2;
+ const uint64_t j2 = j - (j0 * slice1 + j1 * slice2);
+ return TinyVector<Dimension, uint64_t>{j0, j1, j2};
+ };
+
+ auto node_logic_id = [&](size_t r) {
+ const size_t slice1 = node_size[1] * node_size[2];
+ const size_t& slice2 = node_size[2];
+ const uint64_t r0 = r / slice1;
+ const uint64_t r1 = (r - r0 * slice1) / slice2;
+ const uint64_t r2 = r - (r0 * slice1 + r1 * slice2);
+ return TinyVector<Dimension, uint64_t>{r0, r1, r2};
+ };
+
+ auto node_number = [&](const TinyVector<Dimension, uint64_t>& node_logic_id) {
+ return (node_logic_id[0] * node_size[1] + node_logic_id[1]) * node_size[2] + node_logic_id[2];
+ };
+
+ descriptor.cell_to_node_vector.resize(number_of_cells);
+ constexpr size_t nb_node_per_cell = 1 << Dimension;
+ for (size_t j = 0; j < number_of_cells; ++j) {
+ TinyVector<Dimension, size_t> cell_index = cell_logic_id(j);
+ descriptor.cell_to_node_vector[j].resize(nb_node_per_cell);
+ for (size_t r = 0; r < nb_node_per_cell; ++r) {
+ static_assert(Dimension == 3, "unexpected dimension");
+ descriptor.cell_to_node_vector[j][0] = node_number(cell_index + TinyVector<Dimension, uint64_t>{0, 0, 0});
+ descriptor.cell_to_node_vector[j][1] = node_number(cell_index + TinyVector<Dimension, uint64_t>{1, 0, 0});
+ descriptor.cell_to_node_vector[j][2] = node_number(cell_index + TinyVector<Dimension, uint64_t>{1, 1, 0});
+ descriptor.cell_to_node_vector[j][3] = node_number(cell_index + TinyVector<Dimension, uint64_t>{0, 1, 0});
+ descriptor.cell_to_node_vector[j][4] = node_number(cell_index + TinyVector<Dimension, uint64_t>{0, 0, 1});
+ descriptor.cell_to_node_vector[j][5] = node_number(cell_index + TinyVector<Dimension, uint64_t>{1, 0, 1});
+ descriptor.cell_to_node_vector[j][6] = node_number(cell_index + TinyVector<Dimension, uint64_t>{1, 1, 1});
+ descriptor.cell_to_node_vector[j][7] = node_number(cell_index + TinyVector<Dimension, uint64_t>{0, 1, 1});
+ }
+ }
+
+ MeshBuilderBase::_computeCellFaceAndFaceNodeConnectivities<Dimension>(descriptor);
+ MeshBuilderBase::_computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities<Dimension>(descriptor);
+
+ this->_buildBoundaryFaceList(cell_size, descriptor);
+
+ descriptor.cell_owner_vector.resize(number_of_cells);
+ std::fill(descriptor.cell_owner_vector.begin(), descriptor.cell_owner_vector.end(), parallel::rank());
+
+ descriptor.face_owner_vector.resize(descriptor.face_number_vector.size());
+ std::fill(descriptor.face_owner_vector.begin(), descriptor.face_owner_vector.end(), parallel::rank());
+
+ descriptor.edge_owner_vector.resize(descriptor.edge_number_vector.size());
+ std::fill(descriptor.edge_owner_vector.begin(), descriptor.edge_owner_vector.end(), parallel::rank());
+
+ descriptor.node_owner_vector.resize(descriptor.node_number_vector.size());
+ std::fill(descriptor.node_owner_vector.begin(), descriptor.node_owner_vector.end(), parallel::rank());
+
+ std::shared_ptr p_connectivity = Connectivity<Dimension>::build(descriptor);
+ Connectivity<Dimension>& connectivity = *p_connectivity;
+
+ const TinyVector<Dimension> h{(b[0] - a[0]) / cell_size[0], (b[1] - a[1]) / cell_size[1],
+ (b[2] - a[2]) / cell_size[2]};
+
+ NodeValue<TinyVector<Dimension>> xr(connectivity);
+ for (NodeId r = 0; r < number_of_nodes; ++r) {
+ const TinyVector<Dimension, uint64_t> node_index = node_logic_id(r);
+ for (size_t i = 0; i < Dimension; ++i) {
+ xr[r][i] = a[i] + node_index[i] * h[i];
+ }
+ }
+
+ m_mesh = std::make_shared<Mesh<Connectivity<Dimension>>>(p_connectivity, xr);
+}
+
+template <size_t Dimension>
+CartesianMeshBuilder::CartesianMeshBuilder(const TinyVector<Dimension>& a,
+ const TinyVector<Dimension>& b,
+ const TinyVector<Dimension, uint64_t>& size)
+{
+ if (parallel::rank() == 0) {
+ this->_buildCartesianMesh(a, b, size);
+ }
+ this->_dispatch<Dimension>();
+}
+
+template CartesianMeshBuilder::CartesianMeshBuilder(const TinyVector<1>&,
+ const TinyVector<1>&,
+ const TinyVector<1, uint64_t>&);
+
+template CartesianMeshBuilder::CartesianMeshBuilder(const TinyVector<2>&,
+ const TinyVector<2>&,
+ const TinyVector<2, uint64_t>&);
+
+template CartesianMeshBuilder::CartesianMeshBuilder(const TinyVector<3>&,
+ const TinyVector<3>&,
+ const TinyVector<3, uint64_t>&);
diff --git a/src/mesh/CartesianMeshBuilder.hpp b/src/mesh/CartesianMeshBuilder.hpp
new file mode 100644
index 000000000..a4a9ff9e7
--- /dev/null
+++ b/src/mesh/CartesianMeshBuilder.hpp
@@ -0,0 +1,32 @@
+#ifndef CARTESIAN_MESH_BUILDER_HPP
+#define CARTESIAN_MESH_BUILDER_HPP
+
+#include <algebra/TinyVector.hpp>
+#include <mesh/Mesh.hpp>
+#include <mesh/MeshBuilderBase.hpp>
+
+#include <memory>
+
+class CartesianMeshBuilder : public MeshBuilderBase
+{
+ private:
+ template <size_t Dimension>
+ void _buildBoundaryNodeList(const TinyVector<Dimension, uint64_t>& cell_size, ConnectivityDescriptor& descriptor);
+
+ template <size_t Dimension>
+ void _buildBoundaryFaceList(const TinyVector<Dimension, uint64_t>& cell_size, ConnectivityDescriptor& descriptor);
+
+ template <size_t Dimension>
+ void _buildCartesianMesh(const TinyVector<Dimension>& a,
+ const TinyVector<Dimension>& b,
+ const TinyVector<Dimension, uint64_t>& size);
+
+ public:
+ template <size_t Dimension>
+ CartesianMeshBuilder(const TinyVector<Dimension>& a,
+ const TinyVector<Dimension>& b,
+ const TinyVector<Dimension, uint64_t>& size);
+ ~CartesianMeshBuilder() = default;
+};
+
+#endif // CARTESIAN_MESH_BUILDER_HPP
--
GitLab