diff --git a/src/mesh/CMakeLists.txt b/src/mesh/CMakeLists.txt
index 73d686fb20c158ad2693daa28124081902310f61..5baed9b1b0dd6886ea00efdae2e3775f4e0aa338 100644
--- a/src/mesh/CMakeLists.txt
+++ b/src/mesh/CMakeLists.txt
@@ -4,8 +4,9 @@ add_library(
PugsMesh
Connectivity.cpp
ConnectivityComputer.cpp
- GmshReader.cpp
ConnectivityDispatcher.cpp
+ GmshReader.cpp
+ MeshBuilderBase.cpp
SynchronizerManager.cpp)
# Additional dependencies
diff --git a/src/mesh/ConnectivityDispatcher.hpp b/src/mesh/ConnectivityDispatcher.hpp
index 3de9c1c592a0dcc447104ce032c9d6efa9daceca..73ebc548d98a01e4999611927272206df4af1c23 100644
--- a/src/mesh/ConnectivityDispatcher.hpp
+++ b/src/mesh/ConnectivityDispatcher.hpp
@@ -3,7 +3,6 @@
#include <mesh/ItemValue.hpp>
#include <mesh/ItemValueUtils.hpp>
-#include <mesh/Mesh.hpp>
#include <mesh/ConnectivityDescriptor.hpp>
diff --git a/src/mesh/GmshReader.cpp b/src/mesh/GmshReader.cpp
index 629f32801454edcf3e25a3c68aeeaa972a4ea6ff..e99c04722bff61b143d7fbf7c6b09ee881136320 100644
--- a/src/mesh/GmshReader.cpp
+++ b/src/mesh/GmshReader.cpp
@@ -23,218 +23,6 @@
#include <sstream>
#include <unordered_map>
-template <int Dimension>
-void
-GmshReader::_dispatch()
-{
- using ConnectivityType = Connectivity<Dimension>;
- using Rd = TinyVector<Dimension>;
- using MeshType = Mesh<ConnectivityType>;
-
- if (not m_mesh) {
- ConnectivityDescriptor descriptor;
- std::shared_ptr connectivity = ConnectivityType::build(descriptor);
- NodeValue<Rd> xr;
- m_mesh = std::make_shared<MeshType>(connectivity, xr);
- }
- const MeshType& mesh = static_cast<const MeshType&>(*m_mesh);
-
- ConnectivityDispatcher<Dimension> dispatcher(mesh.connectivity());
-
- std::shared_ptr dispatched_connectivity = dispatcher.dispatchedConnectivity();
- NodeValue<Rd> dispatched_xr = dispatcher.dispatch(mesh.xr());
-
- m_mesh = std::make_shared<MeshType>(dispatched_connectivity, dispatched_xr);
-}
-
-template <size_t Dimension>
-class ConnectivityFace;
-
-template <>
-class ConnectivityFace<2>
-{
- public:
- friend struct Hash;
- struct Hash
- {
- size_t
- operator()(const ConnectivityFace& f) const
- {
- size_t hash = 0;
- hash ^= std::hash<unsigned int>()(f.m_node0_id);
- hash ^= std::hash<unsigned int>()(f.m_node1_id) >> 1;
- return hash;
- }
- };
-
- private:
- const std::vector<int>& m_node_number_vector;
-
- unsigned int m_node0_id;
- unsigned int m_node1_id;
-
- bool m_reversed;
-
- public:
- std::vector<unsigned int>
- nodeIdList() const
- {
- return {m_node0_id, m_node1_id};
- }
-
- bool
- reversed() const
- {
- return m_reversed;
- }
-
- PUGS_INLINE
- bool
- operator==(const ConnectivityFace& f) const
- {
- return ((m_node0_id == f.m_node0_id) and (m_node1_id == f.m_node1_id));
- }
-
- PUGS_INLINE
- bool
- operator<(const ConnectivityFace& f) const
- {
- return ((m_node_number_vector[m_node0_id] < m_node_number_vector[f.m_node0_id]) or
- ((m_node_number_vector[m_node0_id] == m_node_number_vector[f.m_node0_id]) and
- (m_node_number_vector[m_node1_id] < m_node_number_vector[f.m_node1_id])));
- }
-
- PUGS_INLINE
- ConnectivityFace(const std::vector<unsigned int>& node_id_list, const std::vector<int>& node_number_vector)
- : m_node_number_vector(node_number_vector)
- {
- Assert(node_id_list.size() == 2);
-
- if (m_node_number_vector[node_id_list[0]] < m_node_number_vector[node_id_list[1]]) {
- m_node0_id = node_id_list[0];
- m_node1_id = node_id_list[1];
- m_reversed = false;
- } else {
- m_node0_id = node_id_list[1];
- m_node1_id = node_id_list[0];
- m_reversed = true;
- }
- }
-
- PUGS_INLINE
- ConnectivityFace(const ConnectivityFace&) = default;
-
- PUGS_INLINE
- ~ConnectivityFace() = default;
-};
-
-template <>
-class ConnectivityFace<3>
-{
- private:
- friend class GmshReader;
- friend struct Hash;
- struct Hash
- {
- size_t
- operator()(const ConnectivityFace& f) const
- {
- size_t hash = 0;
- for (size_t i = 0; i < f.m_node_id_list.size(); ++i) {
- hash ^= std::hash<unsigned int>()(f.m_node_id_list[i]) >> i;
- }
- return hash;
- }
- };
-
- private:
- bool m_reversed;
- std::vector<NodeId::base_type> m_node_id_list;
- const std::vector<int>& m_node_number_vector;
-
- PUGS_INLINE
- std::vector<unsigned int>
- _sort(const std::vector<unsigned int>& node_list)
- {
- const auto min_id = std::min_element(node_list.begin(), node_list.end());
- const int shift = std::distance(node_list.begin(), min_id);
-
- std::vector<unsigned int> rotated_node_list(node_list.size());
- if (node_list[(shift + 1) % node_list.size()] > node_list[(shift + node_list.size() - 1) % node_list.size()]) {
- for (size_t i = 0; i < node_list.size(); ++i) {
- rotated_node_list[i] = node_list[(shift + node_list.size() - i) % node_list.size()];
- m_reversed = true;
- }
- } else {
- for (size_t i = 0; i < node_list.size(); ++i) {
- rotated_node_list[i] = node_list[(shift + i) % node_list.size()];
- }
- }
-
- return rotated_node_list;
- }
-
- public:
- PUGS_INLINE
- const bool&
- reversed() const
- {
- return m_reversed;
- }
-
- PUGS_INLINE
- const std::vector<unsigned int>&
- nodeIdList() const
- {
- return m_node_id_list;
- }
-
- PUGS_INLINE
- ConnectivityFace(const std::vector<unsigned int>& given_node_id_list, const std::vector<int>& node_number_vector)
- : m_reversed(false), m_node_id_list(_sort(given_node_id_list)), m_node_number_vector(node_number_vector)
- {
- ;
- }
-
- public:
- bool
- operator==(const ConnectivityFace& f) const
- {
- if (m_node_id_list.size() == f.nodeIdList().size()) {
- for (size_t j = 0; j < m_node_id_list.size(); ++j) {
- if (m_node_id_list[j] != f.nodeIdList()[j]) {
- return false;
- }
- }
- return true;
- }
- return false;
- }
-
- PUGS_INLINE
- bool
- operator<(const ConnectivityFace& f) const
- {
- const size_t min_nb_nodes = std::min(f.m_node_id_list.size(), m_node_id_list.size());
- for (size_t i = 0; i < min_nb_nodes; ++i) {
- if (m_node_id_list[i] < f.m_node_id_list[i])
- return true;
- if (m_node_id_list[i] != f.m_node_id_list[i])
- return false;
- }
- return m_node_id_list.size() < f.m_node_id_list.size();
- }
-
- PUGS_INLINE
- ConnectivityFace(const ConnectivityFace&) = default;
-
- PUGS_INLINE
- ConnectivityFace() = delete;
-
- PUGS_INLINE
- ~ConnectivityFace() = default;
-};
-
GmshReader::GmshReader(const std::string& filename) : m_filename(filename)
{
if (parallel::rank() == 0) {
@@ -592,293 +380,6 @@ GmshReader::__readPhysicalNames2_2()
}
}
-template <size_t Dimension>
-void
-GmshReader::__computeCellFaceAndFaceNodeConnectivities(ConnectivityDescriptor& descriptor)
-{
- static_assert((Dimension == 2) or (Dimension == 3), "Invalid dimension to compute cell-face connectivities");
- using CellFaceInfo = std::tuple<CellId, unsigned short, bool>;
- using Face = ConnectivityFace<Dimension>;
-
- const auto& node_number_vector = descriptor.node_number_vector;
- Array<unsigned short> cell_nb_faces(descriptor.cell_to_node_vector.size());
- std::map<Face, std::vector<CellFaceInfo>> face_cells_map;
- for (CellId j = 0; j < descriptor.cell_to_node_vector.size(); ++j) {
- const auto& cell_nodes = descriptor.cell_to_node_vector[j];
-
- if constexpr (Dimension == 2) {
- switch (descriptor.cell_type_vector[j]) {
- case CellType::Triangle: {
- cell_nb_faces[j] = 3;
- // face 0
- Face f0({cell_nodes[1], cell_nodes[2]}, node_number_vector);
- face_cells_map[f0].emplace_back(std::make_tuple(j, 0, f0.reversed()));
-
- // face 1
- Face f1({cell_nodes[2], cell_nodes[0]}, node_number_vector);
- face_cells_map[f1].emplace_back(std::make_tuple(j, 1, f1.reversed()));
-
- // face 2
- Face f2({cell_nodes[0], cell_nodes[1]}, node_number_vector);
- face_cells_map[f2].emplace_back(std::make_tuple(j, 2, f2.reversed()));
- break;
- }
- case CellType::Quadrangle: {
- cell_nb_faces[j] = 4;
- // face 0
- Face f0({cell_nodes[0], cell_nodes[1]}, node_number_vector);
- face_cells_map[f0].emplace_back(std::make_tuple(j, 0, f0.reversed()));
-
- // face 1
- Face f1({cell_nodes[1], cell_nodes[2]}, node_number_vector);
- face_cells_map[f1].emplace_back(std::make_tuple(j, 1, f1.reversed()));
-
- // face 2
- Face f2({cell_nodes[2], cell_nodes[3]}, node_number_vector);
- face_cells_map[f2].emplace_back(std::make_tuple(j, 2, f2.reversed()));
-
- // face 3
- Face f3({cell_nodes[3], cell_nodes[0]}, node_number_vector);
- face_cells_map[f3].emplace_back(std::make_tuple(j, 3, f3.reversed()));
- break;
- }
- default: {
- std::ostringstream error_msg;
- error_msg << name(descriptor.cell_type_vector[j]) << ": unexpected cell type in dimension 2";
- throw UnexpectedError(error_msg.str());
- }
- }
- } else if constexpr (Dimension == 3) {
- switch (descriptor.cell_type_vector[j]) {
- case CellType::Tetrahedron: {
- cell_nb_faces[j] = 4;
- // face 0
- Face f0({cell_nodes[1], cell_nodes[2], cell_nodes[3]}, node_number_vector);
- face_cells_map[f0].emplace_back(std::make_tuple(j, 0, f0.reversed()));
-
- // face 1
- Face f1({cell_nodes[0], cell_nodes[3], cell_nodes[2]}, node_number_vector);
- face_cells_map[f1].emplace_back(std::make_tuple(j, 1, f1.reversed()));
-
- // face 2
- Face f2({cell_nodes[0], cell_nodes[1], cell_nodes[3]}, node_number_vector);
- face_cells_map[f2].emplace_back(std::make_tuple(j, 2, f2.reversed()));
-
- // face 3
- Face f3({cell_nodes[0], cell_nodes[2], cell_nodes[1]}, node_number_vector);
- face_cells_map[f3].emplace_back(std::make_tuple(j, 3, f3.reversed()));
- break;
- }
- case CellType::Hexahedron: {
- // face 0
- Face f0({cell_nodes[3], cell_nodes[2], cell_nodes[1], cell_nodes[0]}, node_number_vector);
- face_cells_map[f0].emplace_back(std::make_tuple(j, 0, f0.reversed()));
-
- // face 1
- Face f1({cell_nodes[4], cell_nodes[5], cell_nodes[6], cell_nodes[7]}, node_number_vector);
- face_cells_map[f1].emplace_back(std::make_tuple(j, 1, f1.reversed()));
-
- // face 2
- Face f2({cell_nodes[0], cell_nodes[4], cell_nodes[7], cell_nodes[3]}, node_number_vector);
- face_cells_map[f2].emplace_back(std::make_tuple(j, 2, f2.reversed()));
-
- // face 3
- Face f3({cell_nodes[1], cell_nodes[2], cell_nodes[6], cell_nodes[5]}, node_number_vector);
- face_cells_map[f3].emplace_back(std::make_tuple(j, 3, f3.reversed()));
-
- // face 4
- Face f4({cell_nodes[0], cell_nodes[1], cell_nodes[5], cell_nodes[4]}, node_number_vector);
- face_cells_map[f4].emplace_back(std::make_tuple(j, 4, f4.reversed()));
-
- // face 5
- Face f5({cell_nodes[3], cell_nodes[7], cell_nodes[6], cell_nodes[2]}, node_number_vector);
- face_cells_map[f5].emplace_back(std::make_tuple(j, 5, f5.reversed()));
-
- cell_nb_faces[j] = 6;
- break;
- }
- default: {
- std::ostringstream error_msg;
- error_msg << name(descriptor.cell_type_vector[j]) << ": unexpected cell type in dimension 3";
- throw UnexpectedError(error_msg.str());
- }
- }
- }
- }
-
- {
- descriptor.cell_to_face_vector.resize(descriptor.cell_to_node_vector.size());
- for (CellId j = 0; j < descriptor.cell_to_face_vector.size(); ++j) {
- descriptor.cell_to_face_vector[j].resize(cell_nb_faces[j]);
- }
- FaceId l = 0;
- for (const auto& face_cells_vector : face_cells_map) {
- const auto& cells_vector = face_cells_vector.second;
- for (unsigned short lj = 0; lj < cells_vector.size(); ++lj) {
- const auto& [cell_number, cell_local_face, reversed] = cells_vector[lj];
- descriptor.cell_to_face_vector[cell_number][cell_local_face] = l;
- }
- ++l;
- }
- }
-
- {
- descriptor.cell_face_is_reversed_vector.resize(descriptor.cell_to_node_vector.size());
- for (CellId j = 0; j < descriptor.cell_face_is_reversed_vector.size(); ++j) {
- descriptor.cell_face_is_reversed_vector[j] = Array<bool>(cell_nb_faces[j]);
- }
- for (const auto& face_cells_vector : face_cells_map) {
- const auto& cells_vector = face_cells_vector.second;
- for (unsigned short lj = 0; lj < cells_vector.size(); ++lj) {
- const auto& [cell_number, cell_local_face, reversed] = cells_vector[lj];
- descriptor.cell_face_is_reversed_vector[cell_number][cell_local_face] = reversed;
- }
- }
- }
-
- {
- descriptor.face_to_node_vector.resize(face_cells_map.size());
- int l = 0;
- for (const auto& face_info : face_cells_map) {
- const Face& face = face_info.first;
- descriptor.face_to_node_vector[l] = face.nodeIdList();
- ++l;
- }
- }
-
- {
- // Face numbers may change if numbers are provided in the file
- descriptor.face_number_vector.resize(face_cells_map.size());
- for (size_t l = 0; l < face_cells_map.size(); ++l) {
- descriptor.face_number_vector[l] = l;
- }
- }
-}
-
-template <size_t Dimension>
-void
-GmshReader::__computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities(ConnectivityDescriptor& descriptor)
-{
- static_assert(Dimension == 3, "Invalid dimension to compute face-edge connectivities");
- using FaceEdgeInfo = std::tuple<FaceId, unsigned short, bool>;
- using Edge = ConnectivityFace<2>;
-
- const auto& node_number_vector = descriptor.node_number_vector;
- Array<unsigned short> face_nb_edges(descriptor.face_to_node_vector.size());
- std::map<Edge, std::vector<FaceEdgeInfo>> edge_faces_map;
- for (FaceId l = 0; l < descriptor.face_to_node_vector.size(); ++l) {
- const auto& face_nodes = descriptor.face_to_node_vector[l];
-
- face_nb_edges[l] = face_nodes.size();
- for (size_t r = 0; r < face_nodes.size() - 1; ++r) {
- Edge e({face_nodes[r], face_nodes[r + 1]}, node_number_vector);
- edge_faces_map[e].emplace_back(std::make_tuple(l, r, e.reversed()));
- }
- {
- Edge e({face_nodes[face_nodes.size() - 1], face_nodes[0]}, node_number_vector);
- edge_faces_map[e].emplace_back(std::make_tuple(l, face_nodes.size() - 1, e.reversed()));
- }
- }
-
- std::unordered_map<Edge, EdgeId, typename Edge::Hash> edge_id_map;
- {
- descriptor.face_to_edge_vector.resize(descriptor.face_to_node_vector.size());
- for (FaceId l = 0; l < descriptor.face_to_node_vector.size(); ++l) {
- descriptor.face_to_edge_vector[l].resize(face_nb_edges[l]);
- }
- EdgeId e = 0;
- for (const auto& edge_faces_vector : edge_faces_map) {
- const auto& faces_vector = edge_faces_vector.second;
- for (unsigned short l = 0; l < faces_vector.size(); ++l) {
- const auto& [face_number, face_local_edge, reversed] = faces_vector[l];
- descriptor.face_to_edge_vector[face_number][face_local_edge] = e;
- }
- edge_id_map[edge_faces_vector.first] = e;
- ++e;
- }
- }
-
- {
- descriptor.face_edge_is_reversed_vector.resize(descriptor.face_to_node_vector.size());
- for (FaceId j = 0; j < descriptor.face_edge_is_reversed_vector.size(); ++j) {
- descriptor.face_edge_is_reversed_vector[j] = Array<bool>(face_nb_edges[j]);
- }
- for (const auto& edge_faces_vector : edge_faces_map) {
- const auto& faces_vector = edge_faces_vector.second;
- for (unsigned short lj = 0; lj < faces_vector.size(); ++lj) {
- const auto& [face_number, face_local_edge, reversed] = faces_vector[lj];
- descriptor.face_edge_is_reversed_vector[face_number][face_local_edge] = reversed;
- }
- }
- }
-
- {
- descriptor.edge_to_node_vector.resize(edge_faces_map.size());
- int e = 0;
- for (const auto& edge_info : edge_faces_map) {
- const Edge& edge = edge_info.first;
- descriptor.edge_to_node_vector[e] = edge.nodeIdList();
- ++e;
- }
- }
-
- {
- // Edge numbers may change if numbers are provided in the file
- descriptor.edge_number_vector.resize(edge_faces_map.size());
- for (size_t e = 0; e < edge_faces_map.size(); ++e) {
- descriptor.edge_number_vector[e] = e;
- }
- }
-
- {
- descriptor.cell_to_node_vector.reserve(descriptor.cell_to_node_vector.size());
- for (CellId j = 0; j < descriptor.cell_to_node_vector.size(); ++j) {
- const auto& cell_nodes = descriptor.cell_to_node_vector[j];
-
- switch (descriptor.cell_type_vector[j]) {
- case CellType::Tetrahedron: {
- constexpr int local_edge[6][2] = {{0, 1}, {0, 2}, {0, 3}, {1, 2}, {2, 3}, {3, 1}};
- std::vector<unsigned int> cell_edge_vector;
- cell_edge_vector.reserve(6);
- for (int i_edge = 0; i_edge < 6; ++i_edge) {
- const auto e = local_edge[i_edge];
- Edge edge{{cell_nodes[e[0]], cell_nodes[e[1]]}, node_number_vector};
- auto i = edge_id_map.find(edge);
- if (i == edge_id_map.end()) {
- throw NormalError("could not find this edge");
- }
- cell_edge_vector.push_back(i->second);
- }
- descriptor.cell_to_edge_vector.emplace_back(cell_edge_vector);
- break;
- }
- case CellType::Hexahedron: {
- constexpr int local_edge[12][2] = {{0, 1}, {1, 2}, {2, 3}, {3, 0}, {4, 5}, {5, 6},
- {6, 7}, {7, 4}, {0, 4}, {1, 5}, {2, 6}, {3, 7}};
- std::vector<unsigned int> cell_edge_vector;
- cell_edge_vector.reserve(12);
- for (int i_edge = 0; i_edge < 12; ++i_edge) {
- const auto e = local_edge[i_edge];
- Edge edge{{cell_nodes[e[0]], cell_nodes[e[1]]}, node_number_vector};
- auto i = edge_id_map.find(edge);
- if (i == edge_id_map.end()) {
- throw NormalError("could not find this edge");
- }
- cell_edge_vector.push_back(i->second);
- }
- descriptor.cell_to_edge_vector.emplace_back(cell_edge_vector);
- break;
- }
- default: {
- std::stringstream error_msg;
- error_msg << name(descriptor.cell_type_vector[j]) << ": unexpected cell type in dimension 3";
- throw UnexpectedError(error_msg.str());
- }
- }
- }
- }
-}
-
void
GmshReader::__proceedData()
{
@@ -1177,7 +678,7 @@ GmshReader::__proceedData()
descriptor.addRefItemList(RefCellList(physical_ref_id.refId(), cell_list));
}
- this->__computeCellFaceAndFaceNodeConnectivities<3>(descriptor);
+ MeshBuilderBase::_computeCellFaceAndFaceNodeConnectivities<3>(descriptor);
const auto& node_number_vector = descriptor.node_number_vector;
@@ -1274,7 +775,7 @@ GmshReader::__proceedData()
descriptor.addRefItemList(RefFaceList{physical_ref_id.refId(), face_list});
}
}
- this->__computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities<3>(descriptor);
+ MeshBuilderBase::_computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities<3>(descriptor);
{
using Edge = ConnectivityFace<2>;
@@ -1434,7 +935,7 @@ GmshReader::__proceedData()
descriptor.addRefItemList(RefCellList(physical_ref_id.refId(), cell_list));
}
- this->__computeCellFaceAndFaceNodeConnectivities<2>(descriptor);
+ MeshBuilderBase::_computeCellFaceAndFaceNodeConnectivities<2>(descriptor);
using Face = ConnectivityFace<2>;
const auto& node_number_vector = descriptor.node_number_vector;
diff --git a/src/mesh/GmshReader.hpp b/src/mesh/GmshReader.hpp
index 1d2d70e4f7fa4d9b4812ab50ae31605bf44f4cf8..4edd1fc7b081d67bf74fa132ed4407f43896346d 100644
--- a/src/mesh/GmshReader.hpp
+++ b/src/mesh/GmshReader.hpp
@@ -2,21 +2,19 @@
#define GMSH_READER_HPP
#include <algebra/TinyVector.hpp>
-
-#include <utils/Array.hpp>
-
#include <mesh/IMesh.hpp>
+#include <mesh/MeshBuilderBase.hpp>
#include <mesh/RefId.hpp>
+#include <utils/Array.hpp>
#include <array>
#include <fstream>
#include <map>
+#include <memory>
#include <string>
#include <vector>
-class ConnectivityDescriptor;
-
-class GmshReader
+class GmshReader : public MeshBuilderBase
{
public:
using R3 = TinyVector<3, double>;
@@ -263,24 +261,7 @@ class GmshReader
*/
void __readPhysicalNames2_2();
- template <size_t Dimension>
- void __computeCellFaceAndFaceNodeConnectivities(ConnectivityDescriptor& descriptor);
-
- template <size_t Dimension>
- void __computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities(ConnectivityDescriptor& descriptor);
-
- template <int Dimension>
- void _dispatch();
-
- std::shared_ptr<IMesh> m_mesh;
-
public:
- std::shared_ptr<IMesh>
- mesh() const
- {
- return m_mesh;
- }
-
GmshReader(const std::string& filename);
~GmshReader() = default;
};
diff --git a/src/mesh/MeshBuilderBase.cpp b/src/mesh/MeshBuilderBase.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..74b2e0a2d71e9a56f8447b8b93e7578f72c93702
--- /dev/null
+++ b/src/mesh/MeshBuilderBase.cpp
@@ -0,0 +1,336 @@
+#include <mesh/MeshBuilderBase.hpp>
+
+#include <mesh/Connectivity.hpp>
+#include <mesh/ConnectivityDescriptor.hpp>
+#include <mesh/ConnectivityDispatcher.hpp>
+#include <mesh/ItemId.hpp>
+#include <mesh/Mesh.hpp>
+#include <utils/PugsAssert.hpp>
+#include <utils/PugsMacros.hpp>
+
+#include <vector>
+
+template <int Dimension>
+void
+MeshBuilderBase::_dispatch()
+{
+ if (parallel::size() == 1) {
+ return;
+ }
+
+ using ConnectivityType = Connectivity<Dimension>;
+ using Rd = TinyVector<Dimension>;
+ using MeshType = Mesh<ConnectivityType>;
+
+ if (not m_mesh) {
+ ConnectivityDescriptor descriptor;
+ std::shared_ptr connectivity = ConnectivityType::build(descriptor);
+ NodeValue<Rd> xr;
+ m_mesh = std::make_shared<MeshType>(connectivity, xr);
+ }
+ const MeshType& mesh = static_cast<const MeshType&>(*m_mesh);
+
+ ConnectivityDispatcher<Dimension> dispatcher(mesh.connectivity());
+
+ std::shared_ptr dispatched_connectivity = dispatcher.dispatchedConnectivity();
+ NodeValue<Rd> dispatched_xr = dispatcher.dispatch(mesh.xr());
+
+ m_mesh = std::make_shared<MeshType>(dispatched_connectivity, dispatched_xr);
+}
+
+template <size_t Dimension>
+void
+MeshBuilderBase::_computeCellFaceAndFaceNodeConnectivities(ConnectivityDescriptor& descriptor)
+{
+ static_assert((Dimension == 2) or (Dimension == 3), "Invalid dimension to compute cell-face connectivities");
+ using CellFaceInfo = std::tuple<CellId, unsigned short, bool>;
+ using Face = ConnectivityFace<Dimension>;
+
+ const auto& node_number_vector = descriptor.node_number_vector;
+ Array<unsigned short> cell_nb_faces(descriptor.cell_to_node_vector.size());
+ std::map<Face, std::vector<CellFaceInfo>> face_cells_map;
+ for (CellId j = 0; j < descriptor.cell_to_node_vector.size(); ++j) {
+ const auto& cell_nodes = descriptor.cell_to_node_vector[j];
+
+ if constexpr (Dimension == 2) {
+ switch (descriptor.cell_type_vector[j]) {
+ case CellType::Triangle: {
+ cell_nb_faces[j] = 3;
+ // face 0
+ Face f0({cell_nodes[1], cell_nodes[2]}, node_number_vector);
+ face_cells_map[f0].emplace_back(std::make_tuple(j, 0, f0.reversed()));
+
+ // face 1
+ Face f1({cell_nodes[2], cell_nodes[0]}, node_number_vector);
+ face_cells_map[f1].emplace_back(std::make_tuple(j, 1, f1.reversed()));
+
+ // face 2
+ Face f2({cell_nodes[0], cell_nodes[1]}, node_number_vector);
+ face_cells_map[f2].emplace_back(std::make_tuple(j, 2, f2.reversed()));
+ break;
+ }
+ case CellType::Quadrangle: {
+ cell_nb_faces[j] = 4;
+ // face 0
+ Face f0({cell_nodes[0], cell_nodes[1]}, node_number_vector);
+ face_cells_map[f0].emplace_back(std::make_tuple(j, 0, f0.reversed()));
+
+ // face 1
+ Face f1({cell_nodes[1], cell_nodes[2]}, node_number_vector);
+ face_cells_map[f1].emplace_back(std::make_tuple(j, 1, f1.reversed()));
+
+ // face 2
+ Face f2({cell_nodes[2], cell_nodes[3]}, node_number_vector);
+ face_cells_map[f2].emplace_back(std::make_tuple(j, 2, f2.reversed()));
+
+ // face 3
+ Face f3({cell_nodes[3], cell_nodes[0]}, node_number_vector);
+ face_cells_map[f3].emplace_back(std::make_tuple(j, 3, f3.reversed()));
+ break;
+ }
+ default: {
+ std::ostringstream error_msg;
+ error_msg << name(descriptor.cell_type_vector[j]) << ": unexpected cell type in dimension 2";
+ throw UnexpectedError(error_msg.str());
+ }
+ }
+ } else if constexpr (Dimension == 3) {
+ switch (descriptor.cell_type_vector[j]) {
+ case CellType::Tetrahedron: {
+ cell_nb_faces[j] = 4;
+ // face 0
+ Face f0({cell_nodes[1], cell_nodes[2], cell_nodes[3]}, node_number_vector);
+ face_cells_map[f0].emplace_back(std::make_tuple(j, 0, f0.reversed()));
+
+ // face 1
+ Face f1({cell_nodes[0], cell_nodes[3], cell_nodes[2]}, node_number_vector);
+ face_cells_map[f1].emplace_back(std::make_tuple(j, 1, f1.reversed()));
+
+ // face 2
+ Face f2({cell_nodes[0], cell_nodes[1], cell_nodes[3]}, node_number_vector);
+ face_cells_map[f2].emplace_back(std::make_tuple(j, 2, f2.reversed()));
+
+ // face 3
+ Face f3({cell_nodes[0], cell_nodes[2], cell_nodes[1]}, node_number_vector);
+ face_cells_map[f3].emplace_back(std::make_tuple(j, 3, f3.reversed()));
+ break;
+ }
+ case CellType::Hexahedron: {
+ // face 0
+ Face f0({cell_nodes[3], cell_nodes[2], cell_nodes[1], cell_nodes[0]}, node_number_vector);
+ face_cells_map[f0].emplace_back(std::make_tuple(j, 0, f0.reversed()));
+
+ // face 1
+ Face f1({cell_nodes[4], cell_nodes[5], cell_nodes[6], cell_nodes[7]}, node_number_vector);
+ face_cells_map[f1].emplace_back(std::make_tuple(j, 1, f1.reversed()));
+
+ // face 2
+ Face f2({cell_nodes[0], cell_nodes[4], cell_nodes[7], cell_nodes[3]}, node_number_vector);
+ face_cells_map[f2].emplace_back(std::make_tuple(j, 2, f2.reversed()));
+
+ // face 3
+ Face f3({cell_nodes[1], cell_nodes[2], cell_nodes[6], cell_nodes[5]}, node_number_vector);
+ face_cells_map[f3].emplace_back(std::make_tuple(j, 3, f3.reversed()));
+
+ // face 4
+ Face f4({cell_nodes[0], cell_nodes[1], cell_nodes[5], cell_nodes[4]}, node_number_vector);
+ face_cells_map[f4].emplace_back(std::make_tuple(j, 4, f4.reversed()));
+
+ // face 5
+ Face f5({cell_nodes[3], cell_nodes[7], cell_nodes[6], cell_nodes[2]}, node_number_vector);
+ face_cells_map[f5].emplace_back(std::make_tuple(j, 5, f5.reversed()));
+
+ cell_nb_faces[j] = 6;
+ break;
+ }
+ default: {
+ std::ostringstream error_msg;
+ error_msg << name(descriptor.cell_type_vector[j]) << ": unexpected cell type in dimension 3";
+ throw UnexpectedError(error_msg.str());
+ }
+ }
+ }
+ }
+
+ {
+ descriptor.cell_to_face_vector.resize(descriptor.cell_to_node_vector.size());
+ for (CellId j = 0; j < descriptor.cell_to_face_vector.size(); ++j) {
+ descriptor.cell_to_face_vector[j].resize(cell_nb_faces[j]);
+ }
+ FaceId l = 0;
+ for (const auto& face_cells_vector : face_cells_map) {
+ const auto& cells_vector = face_cells_vector.second;
+ for (unsigned short lj = 0; lj < cells_vector.size(); ++lj) {
+ const auto& [cell_number, cell_local_face, reversed] = cells_vector[lj];
+ descriptor.cell_to_face_vector[cell_number][cell_local_face] = l;
+ }
+ ++l;
+ }
+ }
+
+ {
+ descriptor.cell_face_is_reversed_vector.resize(descriptor.cell_to_node_vector.size());
+ for (CellId j = 0; j < descriptor.cell_face_is_reversed_vector.size(); ++j) {
+ descriptor.cell_face_is_reversed_vector[j] = Array<bool>(cell_nb_faces[j]);
+ }
+ for (const auto& face_cells_vector : face_cells_map) {
+ const auto& cells_vector = face_cells_vector.second;
+ for (unsigned short lj = 0; lj < cells_vector.size(); ++lj) {
+ const auto& [cell_number, cell_local_face, reversed] = cells_vector[lj];
+ descriptor.cell_face_is_reversed_vector[cell_number][cell_local_face] = reversed;
+ }
+ }
+ }
+
+ {
+ descriptor.face_to_node_vector.resize(face_cells_map.size());
+ int l = 0;
+ for (const auto& face_info : face_cells_map) {
+ const Face& face = face_info.first;
+ descriptor.face_to_node_vector[l] = face.nodeIdList();
+ ++l;
+ }
+ }
+
+ {
+ // Face numbers may change if numbers are provided in the file
+ descriptor.face_number_vector.resize(face_cells_map.size());
+ for (size_t l = 0; l < face_cells_map.size(); ++l) {
+ descriptor.face_number_vector[l] = l;
+ }
+ }
+}
+
+template <size_t Dimension>
+void
+MeshBuilderBase::_computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities(ConnectivityDescriptor& descriptor)
+{
+ static_assert(Dimension == 3, "Invalid dimension to compute face-edge connectivities");
+ using FaceEdgeInfo = std::tuple<FaceId, unsigned short, bool>;
+ using Edge = ConnectivityFace<2>;
+
+ const auto& node_number_vector = descriptor.node_number_vector;
+ Array<unsigned short> face_nb_edges(descriptor.face_to_node_vector.size());
+ std::map<Edge, std::vector<FaceEdgeInfo>> edge_faces_map;
+ for (FaceId l = 0; l < descriptor.face_to_node_vector.size(); ++l) {
+ const auto& face_nodes = descriptor.face_to_node_vector[l];
+
+ face_nb_edges[l] = face_nodes.size();
+ for (size_t r = 0; r < face_nodes.size() - 1; ++r) {
+ Edge e({face_nodes[r], face_nodes[r + 1]}, node_number_vector);
+ edge_faces_map[e].emplace_back(std::make_tuple(l, r, e.reversed()));
+ }
+ {
+ Edge e({face_nodes[face_nodes.size() - 1], face_nodes[0]}, node_number_vector);
+ edge_faces_map[e].emplace_back(std::make_tuple(l, face_nodes.size() - 1, e.reversed()));
+ }
+ }
+
+ std::unordered_map<Edge, EdgeId, typename Edge::Hash> edge_id_map;
+ {
+ descriptor.face_to_edge_vector.resize(descriptor.face_to_node_vector.size());
+ for (FaceId l = 0; l < descriptor.face_to_node_vector.size(); ++l) {
+ descriptor.face_to_edge_vector[l].resize(face_nb_edges[l]);
+ }
+ EdgeId e = 0;
+ for (const auto& edge_faces_vector : edge_faces_map) {
+ const auto& faces_vector = edge_faces_vector.second;
+ for (unsigned short l = 0; l < faces_vector.size(); ++l) {
+ const auto& [face_number, face_local_edge, reversed] = faces_vector[l];
+ descriptor.face_to_edge_vector[face_number][face_local_edge] = e;
+ }
+ edge_id_map[edge_faces_vector.first] = e;
+ ++e;
+ }
+ }
+
+ {
+ descriptor.face_edge_is_reversed_vector.resize(descriptor.face_to_node_vector.size());
+ for (FaceId j = 0; j < descriptor.face_edge_is_reversed_vector.size(); ++j) {
+ descriptor.face_edge_is_reversed_vector[j] = Array<bool>(face_nb_edges[j]);
+ }
+ for (const auto& edge_faces_vector : edge_faces_map) {
+ const auto& faces_vector = edge_faces_vector.second;
+ for (unsigned short lj = 0; lj < faces_vector.size(); ++lj) {
+ const auto& [face_number, face_local_edge, reversed] = faces_vector[lj];
+ descriptor.face_edge_is_reversed_vector[face_number][face_local_edge] = reversed;
+ }
+ }
+ }
+
+ {
+ descriptor.edge_to_node_vector.resize(edge_faces_map.size());
+ int e = 0;
+ for (const auto& edge_info : edge_faces_map) {
+ const Edge& edge = edge_info.first;
+ descriptor.edge_to_node_vector[e] = edge.nodeIdList();
+ ++e;
+ }
+ }
+
+ {
+ // Edge numbers may change if numbers are provided in the file
+ descriptor.edge_number_vector.resize(edge_faces_map.size());
+ for (size_t e = 0; e < edge_faces_map.size(); ++e) {
+ descriptor.edge_number_vector[e] = e;
+ }
+ }
+
+ {
+ descriptor.cell_to_node_vector.reserve(descriptor.cell_to_node_vector.size());
+ for (CellId j = 0; j < descriptor.cell_to_node_vector.size(); ++j) {
+ const auto& cell_nodes = descriptor.cell_to_node_vector[j];
+
+ switch (descriptor.cell_type_vector[j]) {
+ case CellType::Tetrahedron: {
+ constexpr int local_edge[6][2] = {{0, 1}, {0, 2}, {0, 3}, {1, 2}, {2, 3}, {3, 1}};
+ std::vector<unsigned int> cell_edge_vector;
+ cell_edge_vector.reserve(6);
+ for (int i_edge = 0; i_edge < 6; ++i_edge) {
+ const auto e = local_edge[i_edge];
+ Edge edge{{cell_nodes[e[0]], cell_nodes[e[1]]}, node_number_vector};
+ auto i = edge_id_map.find(edge);
+ if (i == edge_id_map.end()) {
+ throw NormalError("could not find this edge");
+ }
+ cell_edge_vector.push_back(i->second);
+ }
+ descriptor.cell_to_edge_vector.emplace_back(cell_edge_vector);
+ break;
+ }
+ case CellType::Hexahedron: {
+ constexpr int local_edge[12][2] = {{0, 1}, {1, 2}, {2, 3}, {3, 0}, {4, 5}, {5, 6},
+ {6, 7}, {7, 4}, {0, 4}, {1, 5}, {2, 6}, {3, 7}};
+ std::vector<unsigned int> cell_edge_vector;
+ cell_edge_vector.reserve(12);
+ for (int i_edge = 0; i_edge < 12; ++i_edge) {
+ const auto e = local_edge[i_edge];
+ Edge edge{{cell_nodes[e[0]], cell_nodes[e[1]]}, node_number_vector};
+ auto i = edge_id_map.find(edge);
+ if (i == edge_id_map.end()) {
+ throw NormalError("could not find this edge");
+ }
+ cell_edge_vector.push_back(i->second);
+ }
+ descriptor.cell_to_edge_vector.emplace_back(cell_edge_vector);
+ break;
+ }
+ default: {
+ std::stringstream error_msg;
+ error_msg << name(descriptor.cell_type_vector[j]) << ": unexpected cell type in dimension 3";
+ throw UnexpectedError(error_msg.str());
+ }
+ }
+ }
+ }
+}
+
+template void MeshBuilderBase::_dispatch<1>();
+template void MeshBuilderBase::_dispatch<2>();
+template void MeshBuilderBase::_dispatch<3>();
+
+template void MeshBuilderBase::_computeCellFaceAndFaceNodeConnectivities<2>(ConnectivityDescriptor& descriptor);
+template void MeshBuilderBase::_computeCellFaceAndFaceNodeConnectivities<3>(ConnectivityDescriptor& descriptor);
+
+template void MeshBuilderBase::_computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities<3>(
+ ConnectivityDescriptor& descriptor);
diff --git a/src/mesh/MeshBuilderBase.hpp b/src/mesh/MeshBuilderBase.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..dcca07f5ab32ec50e19d4909ce0f7b5a68d6cf8f
--- /dev/null
+++ b/src/mesh/MeshBuilderBase.hpp
@@ -0,0 +1,230 @@
+#ifndef MESH_BUILDER_BASE_HPP
+#define MESH_BUILDER_BASE_HPP
+
+#include <mesh/IMesh.hpp>
+#include <mesh/ItemId.hpp>
+#include <utils/PugsAssert.hpp>
+#include <utils/PugsMacros.hpp>
+
+#include <memory>
+#include <vector>
+
+class ConnectivityDescriptor;
+
+class MeshBuilderBase
+{
+ protected:
+ template <size_t Dimension>
+ class ConnectivityFace;
+
+ std::shared_ptr<IMesh> m_mesh;
+
+ template <size_t Dimension>
+ static void _computeCellFaceAndFaceNodeConnectivities(ConnectivityDescriptor& descriptor);
+
+ template <size_t Dimension>
+ static void _computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities(ConnectivityDescriptor& descriptor);
+
+ template <int Dimension>
+ void _dispatch();
+
+ public:
+ std::shared_ptr<IMesh>
+ mesh() const
+ {
+ return m_mesh;
+ }
+
+ MeshBuilderBase() = default;
+ ~MeshBuilderBase() = default;
+};
+
+template <>
+class MeshBuilderBase::ConnectivityFace<2>
+{
+ public:
+ friend struct Hash;
+ friend class CartesianMeshBuilder;
+
+ struct Hash
+ {
+ size_t
+ operator()(const ConnectivityFace& f) const
+ {
+ size_t hash = 0;
+ hash ^= std::hash<unsigned int>()(f.m_node0_id);
+ hash ^= std::hash<unsigned int>()(f.m_node1_id) >> 1;
+ return hash;
+ }
+ };
+
+ private:
+ const std::vector<int>& m_node_number_vector;
+
+ unsigned int m_node0_id;
+ unsigned int m_node1_id;
+
+ bool m_reversed;
+
+ public:
+ std::vector<unsigned int>
+ nodeIdList() const
+ {
+ return {m_node0_id, m_node1_id};
+ }
+
+ bool
+ reversed() const
+ {
+ return m_reversed;
+ }
+
+ PUGS_INLINE
+ bool
+ operator==(const ConnectivityFace& f) const
+ {
+ return ((m_node0_id == f.m_node0_id) and (m_node1_id == f.m_node1_id));
+ }
+
+ PUGS_INLINE
+ bool
+ operator<(const ConnectivityFace& f) const
+ {
+ return ((m_node_number_vector[m_node0_id] < m_node_number_vector[f.m_node0_id]) or
+ ((m_node_number_vector[m_node0_id] == m_node_number_vector[f.m_node0_id]) and
+ (m_node_number_vector[m_node1_id] < m_node_number_vector[f.m_node1_id])));
+ }
+
+ PUGS_INLINE
+ ConnectivityFace(const std::vector<unsigned int>& node_id_list, const std::vector<int>& node_number_vector)
+ : m_node_number_vector(node_number_vector)
+ {
+ Assert(node_id_list.size() == 2);
+
+ if (m_node_number_vector[node_id_list[0]] < m_node_number_vector[node_id_list[1]]) {
+ m_node0_id = node_id_list[0];
+ m_node1_id = node_id_list[1];
+ m_reversed = false;
+ } else {
+ m_node0_id = node_id_list[1];
+ m_node1_id = node_id_list[0];
+ m_reversed = true;
+ }
+ }
+
+ PUGS_INLINE
+ ConnectivityFace(const ConnectivityFace&) = default;
+
+ PUGS_INLINE
+ ~ConnectivityFace() = default;
+};
+
+template <>
+class MeshBuilderBase::ConnectivityFace<3>
+{
+ private:
+ friend class GmshReader;
+ friend class CartesianMeshBuilder;
+ friend struct Hash;
+ struct Hash
+ {
+ size_t
+ operator()(const ConnectivityFace& f) const
+ {
+ size_t hash = 0;
+ for (size_t i = 0; i < f.m_node_id_list.size(); ++i) {
+ hash ^= std::hash<unsigned int>()(f.m_node_id_list[i]) >> i;
+ }
+ return hash;
+ }
+ };
+
+ private:
+ bool m_reversed;
+ std::vector<NodeId::base_type> m_node_id_list;
+ const std::vector<int>& m_node_number_vector;
+
+ PUGS_INLINE
+ std::vector<unsigned int>
+ _sort(const std::vector<unsigned int>& node_list)
+ {
+ const auto min_id = std::min_element(node_list.begin(), node_list.end());
+ const int shift = std::distance(node_list.begin(), min_id);
+
+ std::vector<unsigned int> rotated_node_list(node_list.size());
+ if (node_list[(shift + 1) % node_list.size()] > node_list[(shift + node_list.size() - 1) % node_list.size()]) {
+ for (size_t i = 0; i < node_list.size(); ++i) {
+ rotated_node_list[i] = node_list[(shift + node_list.size() - i) % node_list.size()];
+ m_reversed = true;
+ }
+ } else {
+ for (size_t i = 0; i < node_list.size(); ++i) {
+ rotated_node_list[i] = node_list[(shift + i) % node_list.size()];
+ }
+ }
+
+ return rotated_node_list;
+ }
+
+ public:
+ PUGS_INLINE
+ const bool&
+ reversed() const
+ {
+ return m_reversed;
+ }
+
+ PUGS_INLINE
+ const std::vector<unsigned int>&
+ nodeIdList() const
+ {
+ return m_node_id_list;
+ }
+
+ PUGS_INLINE
+ ConnectivityFace(const std::vector<unsigned int>& given_node_id_list, const std::vector<int>& node_number_vector)
+ : m_reversed(false), m_node_id_list(_sort(given_node_id_list)), m_node_number_vector(node_number_vector)
+ {
+ ;
+ }
+
+ public:
+ bool
+ operator==(const ConnectivityFace& f) const
+ {
+ if (m_node_id_list.size() == f.nodeIdList().size()) {
+ for (size_t j = 0; j < m_node_id_list.size(); ++j) {
+ if (m_node_id_list[j] != f.nodeIdList()[j]) {
+ return false;
+ }
+ }
+ return true;
+ }
+ return false;
+ }
+
+ PUGS_INLINE
+ bool
+ operator<(const ConnectivityFace& f) const
+ {
+ const size_t min_nb_nodes = std::min(f.m_node_id_list.size(), m_node_id_list.size());
+ for (size_t i = 0; i < min_nb_nodes; ++i) {
+ if (m_node_id_list[i] < f.m_node_id_list[i])
+ return true;
+ if (m_node_id_list[i] != f.m_node_id_list[i])
+ return false;
+ }
+ return m_node_id_list.size() < f.m_node_id_list.size();
+ }
+
+ PUGS_INLINE
+ ConnectivityFace(const ConnectivityFace&) = default;
+
+ PUGS_INLINE
+ ConnectivityFace() = delete;
+
+ PUGS_INLINE
+ ~ConnectivityFace() = default;
+};
+
+#endif // MESH_BUILDER_BASE_HPP