diff --git a/src/language/modules/MeshModule.cpp b/src/language/modules/MeshModule.cpp
index da8eb6a69b00159c816b8fa2624071df6104bbab..1bb96c9ac41d1b1fd2bcb85c248fb66aac39d1bf 100644
--- a/src/language/modules/MeshModule.cpp
+++ b/src/language/modules/MeshModule.cpp
@@ -167,6 +167,37 @@ MeshModule::MeshModule()
}
));
+
+ this->_addBuiltinFunction("dual", std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<const IMesh>(
+ const std::shared_ptr<const IMesh>&)>>(
+
+ [](const std::shared_ptr<const IMesh>& i_mesh) -> std::shared_ptr<const IMesh> {
+ switch (i_mesh->dimension()) {
+ case 1: {
+ using MeshType = Mesh<Connectivity<1>>;
+
+ std::shared_ptr p_mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
+ return DualMeshManager::instance().getDualMesh(p_mesh);
+ }
+ case 2: {
+ using MeshType = Mesh<Connectivity<2>>;
+
+ std::shared_ptr p_mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
+ return DualMeshManager::instance().getDualMesh(p_mesh);
+ }
+ case 3: {
+ using MeshType = Mesh<Connectivity<3>>;
+
+ std::shared_ptr p_mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
+ return DualMeshManager::instance().getDualMesh(p_mesh);
+ }
+ default: {
+ throw UnexpectedError("invalid dimension");
+ }
+ }
+ }
+
+ ));
}
void
diff --git a/src/mesh/CMakeLists.txt b/src/mesh/CMakeLists.txt
index 8f75c57d649b22f7ee7e63a6fb1af6f24aaefd61..eeef93dd8915ae3f413c27e36c71ae85043d168d 100644
--- a/src/mesh/CMakeLists.txt
+++ b/src/mesh/CMakeLists.txt
@@ -9,7 +9,9 @@ add_library(
ConnectivityDispatcher.cpp
DiamondDualConnectivityBuilder.cpp
DiamondDualMeshBuilder.cpp
+ DualConnectivityBuilder.cpp
DualConnectivityManager.cpp
+ DualMeshBuilder.cpp
DualMeshManager.cpp
GmshReader.cpp
IConnectivity.cpp
diff --git a/src/mesh/ConnectivityToDualConnectivityDataMapper.hpp b/src/mesh/ConnectivityToDualConnectivityDataMapper.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..49f180abb19e1cf1a33943bea1d653a2b30dffc1
--- /dev/null
+++ b/src/mesh/ConnectivityToDualConnectivityDataMapper.hpp
@@ -0,0 +1,144 @@
+#ifndef CONNECTIVITY_TO_DUAL_CONNECTIVITY_DATA_MAPPER_HPP
+#define CONNECTIVITY_TO_DUAL_CONNECTIVITY_DATA_MAPPER_HPP
+
+#include <mesh/Connectivity.hpp>
+#include <mesh/IConnectivityToDualConnectivityDataMapper.hpp>
+#include <mesh/ItemIdToItemIdMap.hpp>
+#include <mesh/ItemValue.hpp>
+#include <utils/Array.hpp>
+#include <utils/PugsAssert.hpp>
+
+template <size_t Dimension>
+class ConnectivityToDualConnectivityDataMapper : public IConnectivityToDualConnectivityDataMapper
+{
+ private:
+ const IConnectivity* m_primal_connectivity;
+ const IConnectivity* m_dual_connectivity;
+
+ NodeIdToNodeIdMap m_primal_node_to_dual_node_map;
+ CellIdToNodeIdMap m_primal_cell_to_dual_node_map;
+ FaceIdToCellIdMap m_primal_face_to_dual_cell_map;
+
+ public:
+ template <typename OriginDataType1, typename OriginDataType2, typename DestinationDataType>
+ void
+ toDualNode(const NodeValue<OriginDataType1>& primal_node_value,
+ const CellValue<OriginDataType2>& primal_cell_value,
+ const NodeValue<DestinationDataType>& dual_node_value) const
+ {
+ static_assert(not std::is_const_v<DestinationDataType>, "destination data type must not be constant");
+ static_assert(std::is_same_v<std::remove_const_t<OriginDataType1>, DestinationDataType>, "incompatible types");
+ static_assert(std::is_same_v<std::remove_const_t<OriginDataType2>, DestinationDataType>, "incompatible types");
+
+ Assert(m_primal_connectivity == primal_cell_value.connectivity_ptr().get());
+ Assert(m_primal_connectivity == primal_node_value.connectivity_ptr().get());
+ Assert(m_dual_connectivity == dual_node_value.connectivity_ptr().get());
+
+ parallel_for(
+ m_primal_node_to_dual_node_map.size(), PUGS_LAMBDA(size_t i) {
+ const auto [primal_node_id, dual_node_id] = m_primal_node_to_dual_node_map[i];
+
+ dual_node_value[dual_node_id] = primal_node_value[primal_node_id];
+ });
+
+ parallel_for(
+ m_primal_cell_to_dual_node_map.size(), PUGS_LAMBDA(size_t i) {
+ const auto [primal_cell_id, dual_node_id] = m_primal_cell_to_dual_node_map[i];
+ dual_node_value[dual_node_id] = primal_cell_value[primal_cell_id];
+ });
+ }
+
+ template <typename OriginDataType, typename DestinationDataType1, typename DestinationDataType2>
+ void
+ fromDualNode(const NodeValue<OriginDataType>& dual_node_value,
+ const NodeValue<DestinationDataType1>& primal_node_value,
+ const CellValue<DestinationDataType2>& primal_cell_value) const
+ {
+ static_assert(not std::is_const_v<DestinationDataType1>, "destination data type must not be constant");
+ static_assert(not std::is_const_v<DestinationDataType2>, "destination data type must not be constant");
+ static_assert(std::is_same_v<std::remove_const_t<OriginDataType>, DestinationDataType1>, "incompatible types");
+ static_assert(std::is_same_v<std::remove_const_t<OriginDataType>, DestinationDataType2>, "incompatible types");
+
+ Assert(m_primal_connectivity == primal_cell_value.connectivity_ptr().get());
+ Assert(m_primal_connectivity == primal_node_value.connectivity_ptr().get());
+ Assert(m_dual_connectivity == dual_node_value.connectivity_ptr().get());
+
+ parallel_for(
+ m_primal_node_to_dual_node_map.size(), PUGS_LAMBDA(size_t i) {
+ const auto [primal_node_id, dual_node_id] = m_primal_node_to_dual_node_map[i];
+
+ primal_node_value[primal_node_id] = dual_node_value[dual_node_id];
+ });
+
+ parallel_for(
+ m_primal_cell_to_dual_node_map.size(), PUGS_LAMBDA(size_t i) {
+ const auto [primal_cell_id, dual_node_id] = m_primal_cell_to_dual_node_map[i];
+ primal_cell_value[primal_cell_id] = dual_node_value[dual_node_id];
+ });
+ }
+
+ template <typename OriginDataType, typename DestinationDataType>
+ void
+ toDualCell(const FaceValue<OriginDataType>& primal_face_value,
+ const CellValue<DestinationDataType>& dual_cell_value) const
+ {
+ static_assert(not std::is_const_v<DestinationDataType>, "destination data type must not be constant");
+ static_assert(std::is_same_v<std::remove_const_t<OriginDataType>, DestinationDataType>, "incompatible types");
+
+ Assert(m_primal_connectivity == primal_face_value.connectivity_ptr().get());
+ Assert(m_dual_connectivity == dual_cell_value.connectivity_ptr().get());
+
+ parallel_for(
+ m_primal_face_to_dual_cell_map.size(), PUGS_LAMBDA(size_t i) {
+ const auto [primal_face_id, dual_cell_id] = m_primal_face_to_dual_cell_map[i];
+
+ dual_cell_value[dual_cell_id] = primal_face_value[primal_face_id];
+ });
+
+ parallel_for(
+ m_primal_cell_to_dual_node_map.size(), PUGS_LAMBDA(size_t i) {
+ const auto [primal_face_id, dual_cell_id] = m_primal_face_to_dual_cell_map[i];
+
+ dual_cell_value[dual_cell_id] = primal_face_value[primal_face_id];
+ });
+ }
+
+ template <typename OriginDataType, typename DestinationDataType>
+ void
+ fromDualCell(const CellValue<DestinationDataType>& dual_cell_value,
+ const FaceValue<OriginDataType>& primal_face_value) const
+ {
+ static_assert(not std::is_const_v<DestinationDataType>, "destination data type must not be constant");
+ static_assert(std::is_same_v<std::remove_const_t<OriginDataType>, DestinationDataType>, "incompatible types");
+
+ Assert(m_primal_connectivity == primal_face_value.connectivity_ptr().get());
+ Assert(m_dual_connectivity == dual_cell_value.connectivity_ptr().get());
+
+ parallel_for(
+ m_primal_face_to_dual_cell_map.size(), PUGS_LAMBDA(size_t i) {
+ const auto [primal_face_id, dual_cell_id] = m_primal_face_to_dual_cell_map[i];
+
+ primal_face_value[primal_face_id] = dual_cell_value[dual_cell_id];
+ });
+
+ parallel_for(
+ m_primal_cell_to_dual_node_map.size(), PUGS_LAMBDA(size_t i) {
+ const auto [primal_face_id, dual_cell_id] = m_primal_face_to_dual_cell_map[i];
+ primal_face_value[primal_face_id] = dual_cell_value[dual_cell_id];
+ });
+ }
+
+ ConnectivityToDualConnectivityDataMapper(const Connectivity<Dimension>& primal_connectivity,
+ const Connectivity<Dimension>& dual_connectivity,
+ const NodeIdToNodeIdMap& primal_node_to_dual_node_map,
+ const CellIdToNodeIdMap& primal_cell_to_dual_node_map,
+ const FaceIdToCellIdMap& primal_face_to_dual_cell_map)
+ : m_primal_connectivity{&primal_connectivity},
+ m_dual_connectivity{&dual_connectivity},
+ m_primal_node_to_dual_node_map{primal_node_to_dual_node_map},
+ m_primal_cell_to_dual_node_map{primal_cell_to_dual_node_map},
+ m_primal_face_to_dual_cell_map{primal_face_to_dual_cell_map}
+ {}
+};
+
+#endif // CONNECTIVITY_TO_DUAL_CONNECTIVITY_DATA_MAPPER_HPP
diff --git a/src/mesh/DiamondDualMeshBuilder.cpp b/src/mesh/DiamondDualMeshBuilder.cpp
index 175aadfa4562f8bb76967a4f26e6ab090c63acaa..6276012f66faf8101a073fcf22e698a6d70fb3ba 100644
--- a/src/mesh/DiamondDualMeshBuilder.cpp
+++ b/src/mesh/DiamondDualMeshBuilder.cpp
@@ -2,7 +2,6 @@
#include <mesh/Connectivity.hpp>
#include <mesh/ConnectivityToDiamondDualConnectivityDataMapper.hpp>
-#include <mesh/DiamondDualConnectivityBuilder.hpp>
#include <mesh/DualConnectivityManager.hpp>
#include <mesh/ItemValueUtils.hpp>
#include <mesh/Mesh.hpp>
diff --git a/src/mesh/DualConnectivityBuilder.cpp b/src/mesh/DualConnectivityBuilder.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..77b00b0f1f62bc03f741ac97cb2d77a9565b2226
--- /dev/null
+++ b/src/mesh/DualConnectivityBuilder.cpp
@@ -0,0 +1,585 @@
+#include <mesh/DualConnectivityBuilder.hpp>
+
+#include <mesh/Connectivity.hpp>
+#include <mesh/ConnectivityDescriptor.hpp>
+#include <mesh/ConnectivityDispatcher.hpp>
+#include <mesh/ConnectivityToDualConnectivityDataMapper.hpp>
+#include <mesh/ItemValueUtils.hpp>
+#include <mesh/Mesh.hpp>
+#include <mesh/RefId.hpp>
+#include <utils/Array.hpp>
+#include <utils/Messenger.hpp>
+
+#include <vector>
+
+template <size_t Dimension>
+void
+DualConnectivityBuilder::_buildConnectivityDescriptor(const Connectivity<Dimension>& primal_connectivity,
+ ConnectivityDescriptor& dual_descriptor)
+{
+ const size_t primal_number_of_nodes = primal_connectivity.numberOfNodes();
+ const size_t primal_number_of_faces = primal_connectivity.numberOfFaces();
+ const size_t primal_number_of_cells = primal_connectivity.numberOfCells();
+
+ const auto& primal_node_number = primal_connectivity.nodeNumber();
+ const auto& primal_cell_number = primal_connectivity.cellNumber();
+
+ const size_t dual_number_of_nodes = primal_number_of_cells + primal_number_of_nodes;
+ const size_t dual_number_of_cells = primal_number_of_faces;
+
+ {
+ m_primal_node_to_dual_node_map = NodeIdToNodeIdMap{primal_number_of_nodes};
+
+ NodeId dual_node_id = 0;
+ for (NodeId primal_node_id = 0; primal_node_id < primal_number_of_nodes; ++primal_node_id) {
+ m_primal_node_to_dual_node_map[primal_node_id] = std::make_pair(primal_node_id, dual_node_id++);
+ }
+
+ m_primal_cell_to_dual_node_map = CellIdToNodeIdMap{primal_number_of_cells};
+ for (CellId primal_cell_id = 0; primal_cell_id < primal_number_of_cells; ++primal_cell_id) {
+ m_primal_cell_to_dual_node_map[primal_cell_id] = std::make_pair(primal_cell_id, dual_node_id++);
+ }
+ }
+
+ dual_descriptor.node_number_vector.resize(dual_number_of_nodes);
+
+ parallel_for(m_primal_node_to_dual_node_map.size(), [&](size_t i) {
+ const auto [primal_node_id, dual_dual_node_id] = m_primal_node_to_dual_node_map[i];
+
+ dual_descriptor.node_number_vector[dual_dual_node_id] = primal_node_number[primal_node_id];
+ });
+
+ const size_t cell_number_shift = max(primal_node_number) + 1;
+ parallel_for(primal_number_of_cells, [&](size_t i) {
+ const auto [primal_cell_id, dual_dual_node_id] = m_primal_cell_to_dual_node_map[i];
+
+ dual_descriptor.node_number_vector[dual_dual_node_id] = primal_cell_number[primal_cell_id] + cell_number_shift;
+ });
+
+ {
+ m_primal_face_to_dual_cell_map = FaceIdToCellIdMap{primal_number_of_faces};
+ CellId dual_cell_id = 0;
+ for (FaceId primal_face_id = 0; primal_face_id < primal_number_of_faces; ++primal_face_id) {
+ m_primal_face_to_dual_cell_map[primal_face_id] = std::make_pair(primal_face_id, dual_cell_id++);
+ }
+ }
+
+ dual_descriptor.cell_number_vector.resize(dual_number_of_cells);
+ const auto& primal_face_number = primal_connectivity.faceNumber();
+ parallel_for(dual_number_of_cells, [&](size_t i) {
+ const auto [primal_face_id, dual_cell_id] = m_primal_face_to_dual_cell_map[i];
+ dual_descriptor.cell_number_vector[dual_cell_id] = primal_face_number[primal_face_id];
+ });
+
+ if constexpr (Dimension == 3) {
+ const size_t number_of_edges = dual_descriptor.edge_to_node_vector.size();
+ dual_descriptor.edge_number_vector.resize(number_of_edges);
+ for (size_t i_edge = 0; i_edge < number_of_edges; ++i_edge) {
+ dual_descriptor.edge_number_vector[i_edge] = i_edge;
+ }
+ if (parallel::size() > 1) {
+ throw NotImplementedError("parallel edge numbering is undefined");
+ }
+ }
+
+ dual_descriptor.cell_type_vector.resize(dual_number_of_cells);
+
+ const auto& primal_face_to_cell_matrix = primal_connectivity.faceToCellMatrix();
+ const auto& primal_face_to_node_matrix = primal_connectivity.faceToNodeMatrix();
+
+ static_assert(Dimension > 1);
+ parallel_for(primal_number_of_faces, [&](FaceId face_id) {
+ const size_t i_cell = face_id;
+ const auto& primal_face_cell_list = primal_face_to_cell_matrix[face_id];
+
+ if constexpr (Dimension == 2) {
+ if (primal_face_cell_list.size() == 1) {
+ dual_descriptor.cell_type_vector[i_cell] = CellType::Triangle;
+ } else {
+ Assert(primal_face_cell_list.size() == 2);
+ dual_descriptor.cell_type_vector[i_cell] = CellType::Quadrangle;
+ }
+ } else if constexpr (Dimension == 3) {
+ if (primal_face_cell_list.size() == 1) {
+ if (primal_face_to_node_matrix[face_id].size() == 3) {
+ dual_descriptor.cell_type_vector[i_cell] = CellType::Tetrahedron;
+ } else {
+ dual_descriptor.cell_type_vector[i_cell] = CellType::Pyramid;
+ }
+ } else {
+ Assert(primal_face_cell_list.size() == 2);
+ dual_descriptor.cell_type_vector[i_cell] = CellType::Diamond;
+ }
+ }
+ });
+
+ dual_descriptor.cell_to_node_vector.resize(dual_number_of_cells);
+
+ const auto& primal_face_local_number_in_their_cells = primal_connectivity.faceLocalNumbersInTheirCells();
+ const auto& cell_face_is_reversed = primal_connectivity.cellFaceIsReversed();
+ parallel_for(primal_number_of_faces, [&](FaceId face_id) {
+ const size_t& i_dual_cell = face_id;
+ const auto& primal_face_cell_list = primal_face_to_cell_matrix[face_id];
+ const auto& primal_face_node_list = primal_face_to_node_matrix[face_id];
+ if (primal_face_cell_list.size() == 1) {
+ dual_descriptor.cell_to_node_vector[i_dual_cell].resize(primal_face_node_list.size() + 1);
+
+ const CellId cell_id = primal_face_cell_list[0];
+ const auto i_face_in_cell = primal_face_local_number_in_their_cells(face_id, 0);
+
+ for (size_t i_node = 0; i_node < primal_face_node_list.size(); ++i_node) {
+ dual_descriptor.cell_to_node_vector[i_dual_cell][i_node] = primal_face_node_list[i_node];
+ }
+ dual_descriptor.cell_to_node_vector[i_dual_cell][primal_face_node_list.size()] = primal_number_of_nodes + cell_id;
+
+ if constexpr (Dimension == 2) {
+ if (cell_face_is_reversed(cell_id, i_face_in_cell)) {
+ std::swap(dual_descriptor.cell_to_node_vector[i_dual_cell][0],
+ dual_descriptor.cell_to_node_vector[i_dual_cell][1]);
+ }
+ } else {
+ if (not cell_face_is_reversed(cell_id, i_face_in_cell)) {
+ // In 3D the basis of the pyramid is described in the
+ // indirect way IF the face is not reversed. In other words
+ // the "topological normal" must point to the "top" of the
+ // pyramid.
+ for (size_t i_node = 0; i_node < primal_face_node_list.size() / 2; ++i_node) {
+ std::swap(dual_descriptor.cell_to_node_vector[i_dual_cell][i_node],
+ dual_descriptor.cell_to_node_vector[i_dual_cell][primal_face_node_list.size() - 1 - i_node]);
+ }
+ }
+ }
+ } else {
+ Assert(primal_face_cell_list.size() == 2);
+ dual_descriptor.cell_to_node_vector[i_dual_cell].resize(primal_face_node_list.size() + 2);
+
+ const CellId cell0_id = primal_face_cell_list[0];
+ const CellId cell1_id = primal_face_cell_list[1];
+ const auto i_face_in_cell0 = primal_face_local_number_in_their_cells(face_id, 0);
+
+ if constexpr (Dimension == 2) {
+ Assert(primal_face_node_list.size() == 2);
+ dual_descriptor.cell_to_node_vector[i_dual_cell][0] = primal_number_of_nodes + cell0_id;
+ dual_descriptor.cell_to_node_vector[i_dual_cell][1] = primal_face_node_list[0];
+ dual_descriptor.cell_to_node_vector[i_dual_cell][2] = primal_number_of_nodes + cell1_id;
+ dual_descriptor.cell_to_node_vector[i_dual_cell][3] = primal_face_node_list[1];
+
+ if (cell_face_is_reversed(cell0_id, i_face_in_cell0)) {
+ std::swap(dual_descriptor.cell_to_node_vector[i_dual_cell][1],
+ dual_descriptor.cell_to_node_vector[i_dual_cell][3]);
+ }
+ } else {
+ dual_descriptor.cell_to_node_vector[i_dual_cell][0] = primal_number_of_nodes + cell0_id;
+ for (size_t i_node = 0; i_node < primal_face_node_list.size(); ++i_node) {
+ dual_descriptor.cell_to_node_vector[i_dual_cell][i_node + 1] = primal_face_node_list[i_node];
+ }
+ dual_descriptor.cell_to_node_vector[i_dual_cell][primal_face_node_list.size() + 1] =
+ primal_number_of_nodes + cell1_id;
+
+ if (cell_face_is_reversed(cell0_id, i_face_in_cell0)) {
+ std::swap(dual_descriptor.cell_to_node_vector[i_dual_cell][0],
+ dual_descriptor.cell_to_node_vector[i_dual_cell][primal_face_node_list.size() + 1]);
+ }
+ }
+ }
+ });
+}
+
+template <>
+void
+DualConnectivityBuilder::_buildConnectivityDescriptor<1>(const Connectivity<1>& primal_connectivity,
+ ConnectivityDescriptor& dual_descriptor)
+{
+ const size_t primal_number_of_nodes = primal_connectivity.numberOfNodes();
+ const size_t primal_number_of_faces = primal_connectivity.numberOfFaces();
+ const size_t primal_number_of_cells = primal_connectivity.numberOfCells();
+
+ const size_t dual_number_of_nodes = 2 * primal_number_of_nodes - primal_number_of_cells;
+
+ const size_t dual_number_of_cells = primal_number_of_faces;
+ const size_t number_of_kept_nodes = 2 * (dual_number_of_nodes - dual_number_of_cells);
+
+ const auto& primal_node_to_cell_matrix = primal_connectivity.nodeToCellMatrix();
+ size_t next_kept_node_id = 0;
+
+ dual_descriptor.node_number_vector.resize(dual_number_of_nodes);
+
+ const auto& primal_node_number = primal_connectivity.nodeNumber();
+
+ for (NodeId primal_node_id = 0; primal_node_id < primal_connectivity.numberOfNodes(); ++primal_node_id) {
+ const auto& primal_node_cell_list = primal_node_to_cell_matrix[primal_node_id];
+ if (primal_node_cell_list.size() == 1) {
+ dual_descriptor.node_number_vector[next_kept_node_id++] = primal_node_number[primal_node_id];
+ }
+ }
+
+ const size_t primal_number_of_kept_nodes = next_kept_node_id;
+
+ const auto& primal_cell_number = primal_connectivity.cellNumber();
+
+ const size_t cell_number_shift = max(primal_node_number) + 1;
+
+ for (CellId primal_cell_id = 0; primal_cell_id < primal_number_of_cells; ++primal_cell_id) {
+ dual_descriptor.node_number_vector[primal_number_of_kept_nodes + primal_cell_id] =
+ primal_cell_number[primal_cell_id] + cell_number_shift;
+ }
+
+ if (number_of_kept_nodes != next_kept_node_id) {
+ throw UnexpectedError("unexpected number of kept node" + std::to_string(next_kept_node_id) +
+ " != " + std::to_string(number_of_kept_nodes));
+ }
+
+ dual_descriptor.cell_number_vector.resize(dual_number_of_cells);
+
+ for (NodeId primal_node_id = 0; primal_node_id < primal_number_of_nodes; ++primal_node_id) {
+ dual_descriptor.cell_number_vector[primal_node_id] = primal_node_number[primal_node_id];
+ }
+
+ dual_descriptor.cell_type_vector.resize(dual_number_of_cells);
+
+ for (NodeId node_id = 0; node_id < primal_connectivity.numberOfNodes(); ++node_id) {
+ const size_t i_cell = node_id;
+
+ dual_descriptor.cell_type_vector[i_cell] = CellType::Line;
+ }
+
+ dual_descriptor.cell_to_node_vector.resize(dual_number_of_cells);
+
+ const auto& primal_node_local_number_in_their_cells = primal_connectivity.nodeLocalNumbersInTheirCells();
+
+ {
+ size_t next_kept_node_id = 0;
+ for (NodeId i_node = 0; i_node < primal_connectivity.numberOfNodes(); ++i_node) {
+ const size_t& i_dual_cell = i_node;
+ const auto& primal_node_cell_list = primal_node_to_cell_matrix[i_node];
+ dual_descriptor.cell_to_node_vector[i_dual_cell].resize(2);
+ if (primal_node_cell_list.size() == 1) {
+ const auto i_node_in_cell = primal_node_local_number_in_their_cells(i_node, 0);
+
+ dual_descriptor.cell_to_node_vector[i_dual_cell][i_node_in_cell] = next_kept_node_id++;
+ dual_descriptor.cell_to_node_vector[i_dual_cell][1 - i_node_in_cell] =
+ number_of_kept_nodes + primal_node_cell_list[0];
+ } else {
+ dual_descriptor.cell_to_node_vector[i_dual_cell][0] = number_of_kept_nodes + primal_node_cell_list[0];
+ dual_descriptor.cell_to_node_vector[i_dual_cell][1] = number_of_kept_nodes + primal_node_cell_list[1];
+ }
+ }
+ }
+}
+
+template <size_t Dimension>
+void
+DualConnectivityBuilder::_buildConnectivityFrom(const IConnectivity& i_primal_connectivity)
+{
+ static_assert(Dimension <= 3, "invalid connectivity dimension");
+
+ if constexpr (Dimension != 2) {
+ throw NotImplementedError("dual mesh in dimension " + std::to_string(Dimension));
+ }
+
+ using ConnectivityType = Connectivity<Dimension>;
+
+ const ConnectivityType& primal_connectivity = dynamic_cast<const ConnectivityType&>(i_primal_connectivity);
+
+ ConnectivityDescriptor dual_descriptor;
+
+ this->_buildConnectivityDescriptor(primal_connectivity, dual_descriptor);
+
+ if constexpr (Dimension > 1) {
+ ConnectivityBuilderBase::_computeCellFaceAndFaceNodeConnectivities<Dimension>(dual_descriptor);
+ if constexpr (Dimension > 2) {
+ ConnectivityBuilderBase::_computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities<Dimension>(dual_descriptor);
+ }
+ }
+
+ {
+ const std::unordered_map<unsigned int, NodeId> node_to_id_map = [&] {
+ std::unordered_map<unsigned int, NodeId> node_to_id_map;
+ for (size_t i_node = 0; i_node < dual_descriptor.node_number_vector.size(); ++i_node) {
+ node_to_id_map[dual_descriptor.node_number_vector[i_node]] = i_node;
+ }
+ return node_to_id_map;
+ }();
+
+ for (size_t i_node_list = 0; i_node_list < primal_connectivity.template numberOfRefItemList<ItemType::node>();
+ ++i_node_list) {
+ const auto& primal_ref_node_list = primal_connectivity.template refItemList<ItemType::node>(i_node_list);
+ const auto& primal_node_list = primal_ref_node_list.list();
+
+ const std::vector<NodeId> dual_node_list = [&]() {
+ std::vector<NodeId> dual_node_list;
+
+ for (size_t i_primal_node = 0; i_primal_node < primal_node_list.size(); ++i_primal_node) {
+ auto primal_node_id = primal_node_list[i_primal_node];
+ const auto i_dual_node = node_to_id_map.find(primal_connectivity.nodeNumber()[primal_node_id]);
+ if (i_dual_node != node_to_id_map.end()) {
+ dual_node_list.push_back(i_dual_node->second);
+ }
+ }
+
+ return dual_node_list;
+ }();
+
+ if (parallel::allReduceOr(dual_node_list.size() > 0)) {
+ Array<NodeId> node_array(dual_node_list.size());
+ for (size_t i = 0; i < dual_node_list.size(); ++i) {
+ node_array[i] = dual_node_list[i];
+ }
+ dual_descriptor.addRefItemList(RefNodeList{primal_ref_node_list.refId(), node_array});
+ }
+ }
+ }
+
+ if constexpr (Dimension > 1) {
+ const auto& primal_face_to_node_matrix = primal_connectivity.faceToNodeMatrix();
+
+ using Face = ConnectivityFace<Dimension>;
+
+ 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 < dual_descriptor.face_to_node_vector.size(); ++l) {
+ const auto& node_vector = dual_descriptor.face_to_node_vector[l];
+
+ face_to_id_map[Face(node_vector, dual_descriptor.node_number_vector)] = l;
+ }
+ return face_to_id_map;
+ }();
+
+ for (size_t i_face_list = 0; i_face_list < primal_connectivity.template numberOfRefItemList<ItemType::face>();
+ ++i_face_list) {
+ const auto& primal_ref_face_list = primal_connectivity.template refItemList<ItemType::face>(i_face_list);
+ const auto& primal_face_list = primal_ref_face_list.list();
+
+ const std::vector<FaceId> dual_face_list = [&]() {
+ std::vector<FaceId> dual_face_list;
+ dual_face_list.reserve(primal_face_list.size());
+ for (size_t i = 0; i < primal_face_list.size(); ++i) {
+ FaceId primal_face_id = primal_face_list[i];
+
+ const auto& primal_face_node_list = primal_face_to_node_matrix[primal_face_id];
+
+ const auto i_dual_face = [&]() {
+ std::vector<unsigned int> node_list(primal_face_node_list.size());
+ for (size_t i = 0; i < primal_face_node_list.size(); ++i) {
+ node_list[i] = primal_face_node_list[i];
+ }
+ return face_to_id_map.find(Face(node_list, dual_descriptor.node_number_vector));
+ }();
+
+ if (i_dual_face != face_to_id_map.end()) {
+ dual_face_list.push_back(i_dual_face->second);
+ }
+ }
+ return dual_face_list;
+ }();
+
+ if (parallel::allReduceOr(dual_face_list.size() > 0)) {
+ Array<FaceId> face_array(dual_face_list.size());
+ for (size_t i = 0; i < dual_face_list.size(); ++i) {
+ face_array[i] = dual_face_list[i];
+ }
+ dual_descriptor.addRefItemList(RefFaceList{primal_ref_face_list.refId(), face_array});
+ }
+ }
+ }
+
+ if constexpr (Dimension > 2) {
+ const auto& primal_edge_to_node_matrix = primal_connectivity.edgeToNodeMatrix();
+ using Edge = ConnectivityFace<2>;
+
+ const std::unordered_map<Edge, EdgeId, typename Edge::Hash> edge_to_id_map = [&] {
+ std::unordered_map<Edge, EdgeId, typename Edge::Hash> edge_to_id_map;
+ for (EdgeId l = 0; l < dual_descriptor.edge_to_node_vector.size(); ++l) {
+ const auto& node_vector = dual_descriptor.edge_to_node_vector[l];
+ edge_to_id_map[Edge(node_vector, dual_descriptor.node_number_vector)] = l;
+ }
+ return edge_to_id_map;
+ }();
+
+ for (size_t i_edge_list = 0; i_edge_list < primal_connectivity.template numberOfRefItemList<ItemType::edge>();
+ ++i_edge_list) {
+ const auto& primal_ref_edge_list = primal_connectivity.template refItemList<ItemType::edge>(i_edge_list);
+ const auto& primal_edge_list = primal_ref_edge_list.list();
+
+ const std::vector<EdgeId> dual_edge_list = [&]() {
+ std::vector<EdgeId> dual_edge_list;
+ dual_edge_list.reserve(primal_edge_list.size());
+ for (size_t i = 0; i < primal_edge_list.size(); ++i) {
+ EdgeId primal_edge_id = primal_edge_list[i];
+
+ const auto& primal_edge_node_list = primal_edge_to_node_matrix[primal_edge_id];
+
+ const auto i_dual_edge = [&]() {
+ std::vector<unsigned int> node_list(primal_edge_node_list.size());
+ for (size_t i = 0; i < primal_edge_node_list.size(); ++i) {
+ node_list[i] = primal_edge_node_list[i];
+ }
+ return edge_to_id_map.find(Edge(node_list, dual_descriptor.node_number_vector));
+ }();
+
+ if (i_dual_edge != edge_to_id_map.end()) {
+ dual_edge_list.push_back(i_dual_edge->second);
+ }
+ }
+ return dual_edge_list;
+ }();
+
+ if (parallel::allReduceOr(dual_edge_list.size() > 0)) {
+ Array<EdgeId> edge_array(dual_edge_list.size());
+ for (size_t i = 0; i < dual_edge_list.size(); ++i) {
+ edge_array[i] = dual_edge_list[i];
+ }
+ dual_descriptor.addRefItemList(RefEdgeList{primal_ref_edge_list.refId(), edge_array});
+ }
+ }
+ }
+
+ const size_t primal_number_of_nodes = primal_connectivity.numberOfNodes();
+ const size_t primal_number_of_cells = primal_connectivity.numberOfCells();
+
+ dual_descriptor.node_owner_vector.resize(dual_descriptor.node_number_vector.size());
+
+ if constexpr (Dimension == 1) {
+ const auto& node_to_cell_matrix = primal_connectivity.nodeToCellMatrix();
+ const auto& primal_node_owner = primal_connectivity.nodeOwner();
+ size_t next_kept_node_id = 0;
+ for (NodeId primal_node_id = 0; primal_node_id < primal_connectivity.numberOfNodes(); ++primal_node_id) {
+ if (node_to_cell_matrix[primal_node_id].size() == 1) {
+ dual_descriptor.node_owner_vector[next_kept_node_id++] = primal_node_owner[primal_node_id];
+ }
+ }
+ const size_t number_of_kept_nodes = next_kept_node_id;
+ const auto& primal_cell_owner = primal_connectivity.cellOwner();
+ for (CellId primal_cell_id = 0; primal_cell_id < primal_number_of_cells; ++primal_cell_id) {
+ dual_descriptor.node_owner_vector[number_of_kept_nodes + primal_cell_id] = primal_cell_owner[primal_cell_id];
+ }
+ } else {
+ const auto& primal_node_owner = primal_connectivity.nodeOwner();
+ for (NodeId primal_node_id = 0; primal_node_id < primal_connectivity.numberOfNodes(); ++primal_node_id) {
+ dual_descriptor.node_owner_vector[primal_node_id] = primal_node_owner[primal_node_id];
+ }
+ const auto& primal_cell_owner = primal_connectivity.cellOwner();
+ for (CellId primal_cell_id = 0; primal_cell_id < primal_number_of_cells; ++primal_cell_id) {
+ dual_descriptor.node_owner_vector[primal_number_of_nodes + primal_cell_id] = primal_cell_owner[primal_cell_id];
+ }
+ }
+
+ if constexpr (Dimension == 1) {
+ dual_descriptor.cell_owner_vector.resize(dual_descriptor.cell_number_vector.size());
+ const auto& primal_node_owner = primal_connectivity.nodeOwner();
+ for (NodeId primal_node_id = 0; primal_node_id < primal_number_of_nodes; ++primal_node_id) {
+ dual_descriptor.cell_owner_vector[primal_node_id] = primal_node_owner[primal_node_id];
+ }
+ } else {
+ dual_descriptor.cell_owner_vector.resize(dual_descriptor.cell_number_vector.size());
+ const size_t primal_number_of_faces = primal_connectivity.numberOfFaces();
+ const auto& primal_face_owner = primal_connectivity.faceOwner();
+ for (FaceId primal_face_id = 0; primal_face_id < primal_number_of_faces; ++primal_face_id) {
+ dual_descriptor.cell_owner_vector[primal_face_id] = primal_face_owner[primal_face_id];
+ }
+ }
+
+ {
+ std::vector<int> face_cell_owner(dual_descriptor.face_number_vector.size());
+ std::fill(std::begin(face_cell_owner), std::end(face_cell_owner), parallel::size());
+
+ for (size_t i_cell = 0; i_cell < dual_descriptor.cell_to_face_vector.size(); ++i_cell) {
+ const auto& cell_face_list = dual_descriptor.cell_to_face_vector[i_cell];
+ for (size_t i_face = 0; i_face < cell_face_list.size(); ++i_face) {
+ const size_t face_id = cell_face_list[i_face];
+ face_cell_owner[face_id] = std::min(face_cell_owner[face_id], dual_descriptor.cell_number_vector[i_cell]);
+ }
+ }
+
+ dual_descriptor.face_owner_vector.resize(face_cell_owner.size());
+ for (size_t i_face = 0; i_face < face_cell_owner.size(); ++i_face) {
+ dual_descriptor.face_owner_vector[i_face] = dual_descriptor.cell_owner_vector[face_cell_owner[i_face]];
+ }
+ }
+
+ if constexpr (Dimension == 3) {
+ std::vector<int> edge_cell_owner(dual_descriptor.edge_number_vector.size());
+ std::fill(std::begin(edge_cell_owner), std::end(edge_cell_owner), parallel::size());
+
+ for (size_t i_cell = 0; i_cell < dual_descriptor.cell_to_face_vector.size(); ++i_cell) {
+ const auto& cell_edge_list = dual_descriptor.cell_to_edge_vector[i_cell];
+ for (size_t i_edge = 0; i_edge < cell_edge_list.size(); ++i_edge) {
+ const size_t edge_id = cell_edge_list[i_edge];
+ edge_cell_owner[edge_id] = std::min(edge_cell_owner[edge_id], dual_descriptor.cell_number_vector[i_cell]);
+ }
+ }
+
+ dual_descriptor.edge_owner_vector.resize(edge_cell_owner.size());
+ for (size_t i_edge = 0; i_edge < edge_cell_owner.size(); ++i_edge) {
+ dual_descriptor.edge_owner_vector[i_edge] = dual_descriptor.cell_owner_vector[edge_cell_owner[i_edge]];
+ }
+ }
+
+ m_connectivity = ConnectivityType::build(dual_descriptor);
+
+ if constexpr (Dimension == 1) {
+ const auto& node_to_cell_matrix = primal_connectivity.nodeToCellMatrix();
+
+ NodeId dual_node_id = 0;
+ m_primal_node_to_dual_node_map = [&]() {
+ std::vector<std::pair<NodeId, NodeId>> primal_node_to_dual_node_vector;
+ for (NodeId primal_node_id = 0; primal_node_id < primal_connectivity.numberOfNodes(); ++primal_node_id) {
+ if (node_to_cell_matrix[primal_node_id].size() == 1) {
+ primal_node_to_dual_node_vector.push_back(std::make_pair(primal_node_id, dual_node_id++));
+ }
+ }
+ return convert_to_array(primal_node_to_dual_node_vector);
+ }();
+
+ m_primal_cell_to_dual_node_map = [&]() {
+ CellIdToNodeIdMap primal_cell_to_dual_node_map{primal_number_of_cells};
+ for (CellId primal_cell_id = 0; primal_cell_id < primal_cell_to_dual_node_map.size(); ++primal_cell_id) {
+ primal_cell_to_dual_node_map[primal_cell_id] = std::make_pair(primal_cell_id, dual_node_id++);
+ }
+ return primal_cell_to_dual_node_map;
+ }();
+
+ m_primal_face_to_dual_cell_map = [&]() {
+ FaceIdToCellIdMap primal_face_to_dual_cell_map{primal_connectivity.numberOfFaces()};
+ for (size_t id = 0; id < primal_face_to_dual_cell_map.size(); ++id) {
+ const CellId dual_cell_id = id;
+ const FaceId primal_face_id = id;
+
+ primal_face_to_dual_cell_map[id] = std::make_pair(primal_face_id, dual_cell_id);
+ }
+ return primal_face_to_dual_cell_map;
+ }();
+ }
+
+ m_mapper =
+ std::make_shared<ConnectivityToDualConnectivityDataMapper<Dimension>>(primal_connectivity,
+ dynamic_cast<const ConnectivityType&>(
+ *m_connectivity),
+ m_primal_node_to_dual_node_map,
+ m_primal_cell_to_dual_node_map,
+ m_primal_face_to_dual_cell_map);
+}
+
+DualConnectivityBuilder::DualConnectivityBuilder(const IConnectivity& connectivity)
+{
+ if (parallel::size() > 1) {
+ throw NotImplementedError("Construction of diamond dual mesh is not implemented in parallel");
+ }
+ switch (connectivity.dimension()) {
+ case 1: {
+ this->_buildConnectivityFrom<1>(connectivity);
+ break;
+ }
+ case 2: {
+ this->_buildConnectivityFrom<2>(connectivity);
+ break;
+ }
+ case 3: {
+ this->_buildConnectivityFrom<3>(connectivity);
+ break;
+ }
+ default: {
+ throw UnexpectedError("invalid connectivity dimension: " + std::to_string(connectivity.dimension()));
+ }
+ }
+}
diff --git a/src/mesh/DualConnectivityBuilder.hpp b/src/mesh/DualConnectivityBuilder.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..2edb39c6c382e79199e9cfd39f57aaab30aef5d8
--- /dev/null
+++ b/src/mesh/DualConnectivityBuilder.hpp
@@ -0,0 +1,42 @@
+#ifndef DUAL_CONNECTIVITY_BUILDER_HPP
+#define DUAL_CONNECTIVITY_BUILDER_HPP
+
+#include <mesh/ConnectivityBuilderBase.hpp>
+#include <mesh/IConnectivityToDualConnectivityDataMapper.hpp>
+#include <mesh/ItemIdToItemIdMap.hpp>
+
+#include <memory>
+
+template <size_t>
+class Connectivity;
+class ConnectivityDescriptor;
+
+class DualConnectivityBuilder : public ConnectivityBuilderBase
+{
+ private:
+ NodeIdToNodeIdMap m_primal_node_to_dual_node_map;
+ CellIdToNodeIdMap m_primal_cell_to_dual_node_map;
+ FaceIdToCellIdMap m_primal_face_to_dual_cell_map;
+
+ std::shared_ptr<IConnectivityToDualConnectivityDataMapper> m_mapper;
+
+ template <size_t Dimension>
+ void _buildConnectivityDescriptor(const Connectivity<Dimension>&, ConnectivityDescriptor&);
+
+ template <size_t Dimension>
+ void _buildConnectivityFrom(const IConnectivity&);
+
+ friend class DualConnectivityManager;
+ DualConnectivityBuilder(const IConnectivity&);
+
+ public:
+ std::shared_ptr<IConnectivityToDualConnectivityDataMapper>
+ mapper() const
+ {
+ return m_mapper;
+ }
+
+ ~DualConnectivityBuilder() = default;
+};
+
+#endif // DUAL_CONNECTIVITY_BUILDER_HPP
diff --git a/src/mesh/DualConnectivityManager.cpp b/src/mesh/DualConnectivityManager.cpp
index 60a4279810ccdb19c3ff6f160a87839e09fb872b..b923cc0205888b25da08eb920100396a60622ef4 100644
--- a/src/mesh/DualConnectivityManager.cpp
+++ b/src/mesh/DualConnectivityManager.cpp
@@ -2,7 +2,9 @@
#include <mesh/Connectivity.hpp>
#include <mesh/ConnectivityToDiamondDualConnectivityDataMapper.hpp>
+#include <mesh/ConnectivityToDualConnectivityDataMapper.hpp>
#include <mesh/DiamondDualConnectivityBuilder.hpp>
+#include <mesh/DualConnectivityBuilder.hpp>
#include <mesh/DualConnectivityManager.hpp>
#include <utils/Exceptions.hpp>
@@ -53,6 +55,16 @@ DualConnectivityManager::deleteConnectivity(const IConnectivity* p_connectivity)
} while (has_removed);
}
+template <typename DualConnectivityBuilderType>
+DualConnectivityManager::DualConnectivityInfo
+DualConnectivityManager::_buildDualConnectivity(const Key& key, const IConnectivity& connectivity)
+{
+ DualConnectivityBuilderType builder{connectivity};
+ DualConnectivityInfo connectivity_info{builder.connectivity(), builder.mapper()};
+ m_connectivity_to_dual_connectivity_info_map[key] = connectivity_info;
+ return connectivity_info;
+}
+
DualConnectivityManager::DualConnectivityInfo
DualConnectivityManager::_getDualConnectivityInfo(const DualMeshType& type, const IConnectivity& connectivity)
{
@@ -63,13 +75,17 @@ DualConnectivityManager::_getDualConnectivityInfo(const DualMeshType& type, cons
i_connectivity != m_connectivity_to_dual_connectivity_info_map.end()) {
return i_connectivity->second;
} else {
- DiamondDualConnectivityBuilder builder{connectivity};
-
- DualConnectivityInfo connectivity_info{builder.connectivity(), builder.mapper()};
-
- m_connectivity_to_dual_connectivity_info_map[key] = connectivity_info;
-
- return connectivity_info;
+ switch (type) {
+ case DualMeshType::Classic: {
+ return this->_buildDualConnectivity<DualConnectivityBuilder>(key, connectivity);
+ }
+ case DualMeshType::Diamond: {
+ return this->_buildDualConnectivity<DiamondDualConnectivityBuilder>(key, connectivity);
+ }
+ default: {
+ throw UnexpectedError("invalid dual mesh type");
+ }
+ }
}
}
@@ -97,6 +113,30 @@ DualConnectivityManager::getConnectivityToDiamondDualConnectivityDataMapper(cons
}
}
+template <size_t Dimension>
+std::shared_ptr<const Connectivity<Dimension>>
+DualConnectivityManager::getDualConnectivity(const Connectivity<Dimension>& connectivity)
+{
+ return std::dynamic_pointer_cast<const Connectivity<Dimension>>(
+ this->_getDualConnectivityInfo(DualMeshType::Classic, connectivity).dualConnectivity());
+}
+
+template <size_t Dimension>
+std::shared_ptr<const ConnectivityToDualConnectivityDataMapper<Dimension>>
+DualConnectivityManager::getConnectivityToDualConnectivityDataMapper(const Connectivity<Dimension>& connectivity)
+{
+ auto i_data_mapper =
+ this->_getDualConnectivityInfo(DualMeshType::Classic, connectivity).connectivityToDualConnectivityDataMapper();
+ auto data_mapper =
+ std::dynamic_pointer_cast<const ConnectivityToDualConnectivityDataMapper<Dimension>>(i_data_mapper);
+
+ if (data_mapper.use_count() > 0) {
+ return data_mapper;
+ } else {
+ throw UnexpectedError("invalid connectivity data mapper type");
+ }
+}
+
template std::shared_ptr<const Connectivity<1>> DualConnectivityManager::getDiamondDualConnectivity(
const Connectivity<1>& connectivity);
@@ -114,3 +154,21 @@ DualConnectivityManager::getConnectivityToDiamondDualConnectivityDataMapper(cons
template std::shared_ptr<const ConnectivityToDiamondDualConnectivityDataMapper<3>>
DualConnectivityManager::getConnectivityToDiamondDualConnectivityDataMapper(const Connectivity<3>&);
+
+template std::shared_ptr<const Connectivity<1>> DualConnectivityManager::getDualConnectivity(
+ const Connectivity<1>& connectivity);
+
+template std::shared_ptr<const Connectivity<2>> DualConnectivityManager::getDualConnectivity(
+ const Connectivity<2>& connectivity);
+
+template std::shared_ptr<const Connectivity<3>> DualConnectivityManager::getDualConnectivity(
+ const Connectivity<3>& connectivity);
+
+template std::shared_ptr<const ConnectivityToDualConnectivityDataMapper<1>>
+DualConnectivityManager::getConnectivityToDualConnectivityDataMapper(const Connectivity<1>&);
+
+template std::shared_ptr<const ConnectivityToDualConnectivityDataMapper<2>>
+DualConnectivityManager::getConnectivityToDualConnectivityDataMapper(const Connectivity<2>&);
+
+template std::shared_ptr<const ConnectivityToDualConnectivityDataMapper<3>>
+DualConnectivityManager::getConnectivityToDualConnectivityDataMapper(const Connectivity<3>&);
diff --git a/src/mesh/DualConnectivityManager.hpp b/src/mesh/DualConnectivityManager.hpp
index eb6e0cad10862df32bf330deefab00b5c7aaf658..1f6819a89842be7c1973587def875f835befd9f9 100644
--- a/src/mesh/DualConnectivityManager.hpp
+++ b/src/mesh/DualConnectivityManager.hpp
@@ -14,6 +14,9 @@ class Connectivity;
template <size_t Dimension>
class ConnectivityToDiamondDualConnectivityDataMapper;
+template <size_t Dimension>
+class ConnectivityToDualConnectivityDataMapper;
+
class DualConnectivityManager
{
private:
@@ -55,8 +58,6 @@ class DualConnectivityManager
~DualConnectivityInfo() = default;
};
- DualConnectivityInfo _getDualConnectivityInfo(const DualMeshType& type, const IConnectivity& connectivity);
-
using Key = std::pair<DualMeshType, const IConnectivity*>;
struct HashKey
{
@@ -67,6 +68,11 @@ class DualConnectivityManager
}
};
+ template <typename DualConnectivityBuilderType>
+ DualConnectivityInfo _buildDualConnectivity(const Key& key, const IConnectivity& connectivity);
+
+ DualConnectivityInfo _getDualConnectivityInfo(const DualMeshType& type, const IConnectivity& connectivity);
+
std::unordered_map<Key, DualConnectivityInfo, HashKey> m_connectivity_to_dual_connectivity_info_map;
static DualConnectivityManager* m_instance;
@@ -91,6 +97,13 @@ class DualConnectivityManager
void deleteConnectivity(const IConnectivity*);
+ template <size_t Dimension>
+ std::shared_ptr<const Connectivity<Dimension>> getDualConnectivity(const Connectivity<Dimension>&);
+
+ template <size_t Dimension>
+ std::shared_ptr<const ConnectivityToDualConnectivityDataMapper<Dimension>>
+ getConnectivityToDualConnectivityDataMapper(const Connectivity<Dimension>&);
+
template <size_t Dimension>
std::shared_ptr<const Connectivity<Dimension>> getDiamondDualConnectivity(const Connectivity<Dimension>&);
diff --git a/src/mesh/DualMeshBuilder.cpp b/src/mesh/DualMeshBuilder.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..47cac173e0a04500a3f7e4a83b8252375a306f0a
--- /dev/null
+++ b/src/mesh/DualMeshBuilder.cpp
@@ -0,0 +1,62 @@
+#include <mesh/DualMeshBuilder.hpp>
+
+#include <mesh/Connectivity.hpp>
+#include <mesh/ConnectivityToDualConnectivityDataMapper.hpp>
+#include <mesh/DualConnectivityManager.hpp>
+#include <mesh/ItemValueUtils.hpp>
+#include <mesh/Mesh.hpp>
+#include <mesh/MeshData.hpp>
+#include <mesh/MeshDataManager.hpp>
+
+template <size_t Dimension>
+void
+DualMeshBuilder::_buildDualMeshFrom(const std::shared_ptr<const IMesh>& p_i_mesh)
+{
+ using ConnectivityType = Connectivity<Dimension>;
+ using MeshType = Mesh<Connectivity<Dimension>>;
+
+ std::shared_ptr p_primal_mesh = std::dynamic_pointer_cast<const MeshType>(p_i_mesh);
+ const MeshType& primal_mesh = *p_primal_mesh;
+
+ DualConnectivityManager& manager = DualConnectivityManager::instance();
+
+ std::shared_ptr<const ConnectivityType> p_dual_connectivity = manager.getDualConnectivity(primal_mesh.connectivity());
+
+ const ConnectivityType& dual_connectivity = *p_dual_connectivity;
+
+ const NodeValue<const TinyVector<Dimension>> primal_xr = primal_mesh.xr();
+
+ MeshData<Dimension>& primal_mesh_data = MeshDataManager::instance().getMeshData(primal_mesh);
+ const CellValue<const TinyVector<Dimension>> primal_xj = primal_mesh_data.xj();
+
+ std::shared_ptr connectivity_to_dual_connectivity_data_mapper =
+ manager.getConnectivityToDualConnectivityDataMapper(primal_mesh.connectivity());
+
+ NodeValue<TinyVector<Dimension>> dual_xr{dual_connectivity};
+ connectivity_to_dual_connectivity_data_mapper->toDualNode(primal_xr, primal_xj, dual_xr);
+
+ m_mesh = std::make_shared<MeshType>(p_dual_connectivity, dual_xr);
+}
+
+DualMeshBuilder::DualMeshBuilder(const std::shared_ptr<const IMesh>& p_mesh)
+{
+ std::cout << "building DualMesh\n";
+
+ switch (p_mesh->dimension()) {
+ case 1: {
+ this->_buildDualMeshFrom<1>(p_mesh);
+ break;
+ }
+ case 2: {
+ this->_buildDualMeshFrom<2>(p_mesh);
+ break;
+ }
+ case 3: {
+ this->_buildDualMeshFrom<3>(p_mesh);
+ break;
+ }
+ default: {
+ throw UnexpectedError("invalid mesh dimension: " + std::to_string(p_mesh->dimension()));
+ }
+ }
+}
diff --git a/src/mesh/DualMeshBuilder.hpp b/src/mesh/DualMeshBuilder.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..1756d7e92ad7e8bfc4fffb2ca2e41ac1bbcb78f1
--- /dev/null
+++ b/src/mesh/DualMeshBuilder.hpp
@@ -0,0 +1,21 @@
+#ifndef DUAL_MESH_BUILDER_HPP
+#define DUAL_MESH_BUILDER_HPP
+
+#include <mesh/MeshBuilderBase.hpp>
+
+#include <memory>
+
+class DualMeshBuilder : public MeshBuilderBase
+{
+ private:
+ template <size_t Dimension>
+ void _buildDualMeshFrom(const std::shared_ptr<const IMesh>&);
+
+ friend class DualMeshManager;
+ DualMeshBuilder(const std::shared_ptr<const IMesh>&);
+
+ public:
+ ~DualMeshBuilder() = default;
+};
+
+#endif // DUAL_MESH_BUILDER_HPP
diff --git a/src/mesh/DualMeshManager.cpp b/src/mesh/DualMeshManager.cpp
index cf903145b65709d53aec0210b731c57e51fe878d..4fa108cf461e1a070725cf14866b41f346e8b20f 100644
--- a/src/mesh/DualMeshManager.cpp
+++ b/src/mesh/DualMeshManager.cpp
@@ -2,6 +2,7 @@
#include <mesh/Connectivity.hpp>
#include <mesh/DiamondDualMeshBuilder.hpp>
+#include <mesh/DualMeshBuilder.hpp>
#include <mesh/Mesh.hpp>
#include <utils/Exceptions.hpp>
#include <utils/PugsAssert.hpp>
@@ -52,6 +53,23 @@ DualMeshManager::deleteMesh(const IMesh* p_mesh)
} while (has_removed);
}
+template <size_t Dimension>
+std::shared_ptr<const Mesh<Connectivity<Dimension>>>
+DualMeshManager::getDualMesh(std::shared_ptr<const Mesh<Connectivity<Dimension>>> mesh)
+{
+ const IMesh* p_mesh = mesh.get();
+
+ auto key = std::make_pair(DualMeshType::Classic, p_mesh);
+ if (auto i_mesh_data = m_mesh_to_dual_mesh_map.find(key); i_mesh_data != m_mesh_to_dual_mesh_map.end()) {
+ return std::dynamic_pointer_cast<const Mesh<Connectivity<Dimension>>>(i_mesh_data->second);
+ } else {
+ DualMeshBuilder builder{mesh};
+
+ m_mesh_to_dual_mesh_map[key] = builder.mesh();
+ return std::dynamic_pointer_cast<const Mesh<Connectivity<Dimension>>>(builder.mesh());
+ }
+}
+
template <size_t Dimension>
std::shared_ptr<const Mesh<Connectivity<Dimension>>>
DualMeshManager::getDiamondDualMesh(std::shared_ptr<const Mesh<Connectivity<Dimension>>> mesh)
@@ -69,6 +87,13 @@ DualMeshManager::getDiamondDualMesh(std::shared_ptr<const Mesh<Connectivity<Dime
}
}
+template std::shared_ptr<const Mesh<Connectivity<1>>> DualMeshManager::getDualMesh(
+ std::shared_ptr<const Mesh<Connectivity<1>>>);
+template std::shared_ptr<const Mesh<Connectivity<2>>> DualMeshManager::getDualMesh(
+ std::shared_ptr<const Mesh<Connectivity<2>>>);
+template std::shared_ptr<const Mesh<Connectivity<3>>> DualMeshManager::getDualMesh(
+ std::shared_ptr<const Mesh<Connectivity<3>>>);
+
template std::shared_ptr<const Mesh<Connectivity<1>>> DualMeshManager::getDiamondDualMesh(
std::shared_ptr<const Mesh<Connectivity<1>>>);
template std::shared_ptr<const Mesh<Connectivity<2>>> DualMeshManager::getDiamondDualMesh(
diff --git a/src/mesh/DualMeshManager.hpp b/src/mesh/DualMeshManager.hpp
index e2b92b599254bbac404478fd590df6747281ae29..dd8af5890564f459034905158d886ee600811ce3 100644
--- a/src/mesh/DualMeshManager.hpp
+++ b/src/mesh/DualMeshManager.hpp
@@ -55,6 +55,10 @@ class DualMeshManager
template <size_t Dimension>
std::shared_ptr<const Mesh<Connectivity<Dimension>>> getDiamondDualMesh(
std::shared_ptr<const Mesh<Connectivity<Dimension>>>);
+
+ template <size_t Dimension>
+ std::shared_ptr<const Mesh<Connectivity<Dimension>>> getDualMesh(
+ std::shared_ptr<const Mesh<Connectivity<Dimension>>>);
};
#endif // DUAL_MESH_MANAGER_HPP