diff --git a/src/language/modules/MeshModule.cpp b/src/language/modules/MeshModule.cpp
index b9a65729637128081a7206ad6fa5d98ad95e9978..89e678c3fb4095350842b19a9145502cd7f558a8 100644
--- a/src/language/modules/MeshModule.cpp
+++ b/src/language/modules/MeshModule.cpp
@@ -9,7 +9,7 @@
#include <language/utils/TypeDescriptor.hpp>
#include <mesh/CartesianMeshBuilder.hpp>
#include <mesh/Connectivity.hpp>
-#include <mesh/DiamondDualMeshBuilder.hpp>
+#include <mesh/DiamondDualMeshManager.hpp>
#include <mesh/GmshReader.hpp>
#include <mesh/Mesh.hpp>
#include <utils/Exceptions.hpp>
@@ -184,9 +184,30 @@ MeshModule::MeshModule()
std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<const IMesh>(
const std::shared_ptr<const IMesh>&)>>(
- [](const std::shared_ptr<const IMesh>& p_mesh) -> std::shared_ptr<const IMesh> {
- DiamondDualMeshBuilder builder{p_mesh};
- return builder.mesh();
+ [](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 DiamondDualMeshManager::instance().getDiamondDualMesh(p_mesh);
+ }
+ case 2: {
+ using MeshType = Mesh<Connectivity<2>>;
+
+ std::shared_ptr p_mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
+ return DiamondDualMeshManager::instance().getDiamondDualMesh(p_mesh);
+ }
+ case 3: {
+ using MeshType = Mesh<Connectivity<3>>;
+
+ std::shared_ptr p_mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
+ return DiamondDualMeshManager::instance().getDiamondDualMesh(p_mesh);
+ }
+ default: {
+ throw UnexpectedError("invalid dimension");
+ }
+ }
}
));
diff --git a/src/language/modules/VTKModule.cpp b/src/language/modules/VTKModule.cpp
index f815a78e763b387ee5707d9ac65ec1f038a36b14..ce84141ee866f7ac6fdbb8448db78d0d1b64de88 100644
--- a/src/language/modules/VTKModule.cpp
+++ b/src/language/modules/VTKModule.cpp
@@ -24,7 +24,10 @@ VTKModule::VTKModule()
const std::shared_ptr<const MeshType> mesh =
std::dynamic_pointer_cast<const MeshType>(p_mesh);
- writer.write(mesh, OutputNamedItemValueSet{}, time, true);
+ writer.write(mesh,
+ OutputNamedItemValueSet{
+ NamedItemValue{"cell_number", mesh->connectivity().cellNumber()}},
+ time, true);
break;
}
case 2: {
@@ -32,7 +35,10 @@ VTKModule::VTKModule()
const std::shared_ptr<const MeshType> mesh =
std::dynamic_pointer_cast<const MeshType>(p_mesh);
- writer.write(mesh, OutputNamedItemValueSet{}, time, true);
+ writer.write(mesh,
+ OutputNamedItemValueSet{
+ NamedItemValue{"cell_number", mesh->connectivity().cellNumber()}},
+ time, true);
break;
}
case 3: {
@@ -40,7 +46,10 @@ VTKModule::VTKModule()
const std::shared_ptr<const MeshType> mesh =
std::dynamic_pointer_cast<const MeshType>(p_mesh);
- writer.write(mesh, OutputNamedItemValueSet{}, time, true);
+ writer.write(mesh,
+ OutputNamedItemValueSet{
+ NamedItemValue{"cell_number", mesh->connectivity().cellNumber()}},
+ time, true);
break;
}
}
diff --git a/src/main.cpp b/src/main.cpp
index 9b853028133fd9035f67f743c62ee2a0a298202a..aad67d58a5d228cdc502e9b024cfc497fef23291 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -1,6 +1,8 @@
#include <utils/PugsUtils.hpp>
#include <language/PugsParser.hpp>
+#include <mesh/DiamondDualConnectivityManager.hpp>
+#include <mesh/DiamondDualMeshManager.hpp>
#include <mesh/MeshDataManager.hpp>
#include <mesh/SynchronizerManager.hpp>
@@ -11,9 +13,13 @@ main(int argc, char* argv[])
SynchronizerManager::create();
MeshDataManager::create();
+ DiamondDualConnectivityManager::create();
+ DiamondDualMeshManager::create();
parser(filename);
+ DiamondDualMeshManager::destroy();
+ DiamondDualConnectivityManager::destroy();
MeshDataManager::destroy();
SynchronizerManager::destroy();
finalize();
diff --git a/src/mesh/CMakeLists.txt b/src/mesh/CMakeLists.txt
index 92c865c599767db0d5e5fed84e7d16d0cbec8579..6c0088831c86023e31914b953c9ea137e5411d21 100644
--- a/src/mesh/CMakeLists.txt
+++ b/src/mesh/CMakeLists.txt
@@ -8,8 +8,11 @@ add_library(
ConnectivityComputer.cpp
ConnectivityDispatcher.cpp
DiamondDualConnectivityBuilder.cpp
+ DiamondDualConnectivityManager.cpp
DiamondDualMeshBuilder.cpp
+ DiamondDualMeshManager.cpp
GmshReader.cpp
+ IConnectivity.cpp
IMesh.cpp
LogicalConnectivityBuilder.cpp
MeshBuilderBase.cpp
diff --git a/src/mesh/Connectivity.hpp b/src/mesh/Connectivity.hpp
index 5a888190841e7886d4548a5405dc3126223318ac..3bf5cd4201c574749e5b6a67a2b6c4d667bd88f1 100644
--- a/src/mesh/Connectivity.hpp
+++ b/src/mesh/Connectivity.hpp
@@ -12,7 +12,6 @@
#include <mesh/RefId.hpp>
#include <mesh/RefItemList.hpp>
#include <mesh/SubItemValuePerItem.hpp>
-#include <mesh/SynchronizerManager.hpp>
#include <utils/CSRGraph.hpp>
#include <utils/Exceptions.hpp>
#include <utils/PugsAssert.hpp>
@@ -651,11 +650,7 @@ class Connectivity final : public IConnectivity
void _buildFrom(const ConnectivityDescriptor& descriptor);
public:
- ~Connectivity()
- {
- auto& manager = SynchronizerManager::instance();
- manager.deleteConnectivitySynchronizer(this);
- }
+ ~Connectivity() = default;
};
template <size_t Dim>
diff --git a/src/mesh/ConnectivityBuilderBase.cpp b/src/mesh/ConnectivityBuilderBase.cpp
index f5d47e25003a710e1e2c18ba6b772ffd04bb532f..60c38bfb84fb975c0e2cdc3b07d1d97e39eb0021 100644
--- a/src/mesh/ConnectivityBuilderBase.cpp
+++ b/src/mesh/ConnectivityBuilderBase.cpp
@@ -2,12 +2,13 @@
#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 <map>
+#include <unordered_map>
#include <vector>
template <size_t Dimension>
diff --git a/src/mesh/ConnectivityToDiamondDualConnectivityDataMapper.hpp b/src/mesh/ConnectivityToDiamondDualConnectivityDataMapper.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..8d35d2ff9541fbfab7e8e6e30976bd5f88970f56
--- /dev/null
+++ b/src/mesh/ConnectivityToDiamondDualConnectivityDataMapper.hpp
@@ -0,0 +1,153 @@
+#ifndef CONNECTIVITY_TO_DIAMOND_DUAL_CONNECTIVITY_DATA_MAPPER_HPP
+#define CONNECTIVITY_TO_DIAMOND_DUAL_CONNECTIVITY_DATA_MAPPER_HPP
+
+#include <mesh/Connectivity.hpp>
+#include <mesh/ItemIdToItemIdMap.hpp>
+#include <mesh/ItemValue.hpp>
+#include <utils/Array.hpp>
+#include <utils/PugsAssert.hpp>
+
+class IConnectivityToDiamondDualConnectivityDataMapper
+{
+ public:
+ IConnectivityToDiamondDualConnectivityDataMapper(const IConnectivityToDiamondDualConnectivityDataMapper&) = delete;
+ IConnectivityToDiamondDualConnectivityDataMapper(IConnectivityToDiamondDualConnectivityDataMapper&&) = delete;
+
+ IConnectivityToDiamondDualConnectivityDataMapper() = default;
+ virtual ~IConnectivityToDiamondDualConnectivityDataMapper() = default;
+};
+
+template <size_t Dimension>
+class ConnectivityToDiamondDualConnectivityDataMapper : public IConnectivityToDiamondDualConnectivityDataMapper
+{
+ 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];
+ });
+ }
+
+ ConnectivityToDiamondDualConnectivityDataMapper(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_DIAMOND_DUAL_CONNECTIVITY_DATA_MAPPER_HPP
diff --git a/src/mesh/DiamondDualConnectivityBuilder.cpp b/src/mesh/DiamondDualConnectivityBuilder.cpp
index 8bcbe9bad8e15cbef1dfa449c0d6defc115f126f..76adb60de4b047cd5d456662db96d10c6e1677e7 100644
--- a/src/mesh/DiamondDualConnectivityBuilder.cpp
+++ b/src/mesh/DiamondDualConnectivityBuilder.cpp
@@ -3,12 +3,16 @@
#include <mesh/Connectivity.hpp>
#include <mesh/ConnectivityDescriptor.hpp>
#include <mesh/ConnectivityDispatcher.hpp>
+#include <mesh/ConnectivityToDiamondDualConnectivityDataMapper.hpp>
#include <mesh/ItemValueUtils.hpp>
#include <mesh/Mesh.hpp>
#include <mesh/RefId.hpp>
#include <utils/Array.hpp>
+#include <utils/ArrayUtils.hpp>
#include <utils/Messenger.hpp>
+#include <vector>
+
template <size_t Dimension>
void
DiamondDualConnectivityBuilder::_buildDiamondConnectivityDescriptor(const Connectivity<Dimension>& primal_connectivity,
@@ -18,31 +22,55 @@ DiamondDualConnectivityBuilder::_buildDiamondConnectivityDescriptor(const Connec
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 diamond_number_of_nodes = primal_number_of_cells + primal_number_of_nodes;
+ const size_t diamond_number_of_cells = primal_number_of_faces;
+
+ {
+ m_primal_node_to_dual_node_map = NodeIdToNodeIdMap{primal_number_of_nodes};
+
+ NodeId diamond_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, diamond_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, diamond_node_id++);
+ }
+ }
diamond_descriptor.node_number_vector.resize(diamond_number_of_nodes);
- const auto& primal_node_number = primal_connectivity.nodeNumber();
+ for (size_t i = 0; i < m_primal_node_to_dual_node_map.size(); ++i) {
+ const auto [primal_node_id, diamond_dual_node_id] = m_primal_node_to_dual_node_map[i];
- for (NodeId primal_node_id = 0; primal_node_id < primal_connectivity.numberOfNodes(); ++primal_node_id) {
- diamond_descriptor.node_number_vector[primal_node_id] = primal_node_number[primal_node_id];
+ diamond_descriptor.node_number_vector[diamond_dual_node_id] = primal_node_number[primal_node_id];
}
- const auto& primal_cell_number = primal_connectivity.cellNumber();
+ const size_t cell_number_shift = max(primal_node_number) + 1;
+ for (size_t i = 0; i < primal_number_of_cells; ++i) {
+ const auto [primal_cell_id, diamond_dual_node_id] = m_primal_cell_to_dual_node_map[i];
- const size_t max_node_number = max(primal_node_number);
+ diamond_descriptor.node_number_vector[diamond_dual_node_id] =
+ primal_cell_number[primal_cell_id] + cell_number_shift;
+ }
- for (CellId primal_cell_id = 0; primal_cell_id < primal_number_of_cells; ++primal_cell_id) {
- diamond_descriptor.node_number_vector[primal_number_of_nodes + primal_cell_id] =
- primal_cell_number[primal_cell_id] + max_node_number;
+ {
+ m_primal_face_to_dual_cell_map = FaceIdToCellIdMap{primal_number_of_faces};
+ CellId diamond_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, diamond_cell_id++);
+ }
}
- const size_t diamond_number_of_cells = primal_number_of_faces;
diamond_descriptor.cell_number_vector.resize(diamond_number_of_cells);
-
const auto& primal_face_number = primal_connectivity.faceNumber();
- for (FaceId i_primal_face = 0; i_primal_face < primal_number_of_faces; ++i_primal_face) {
- diamond_descriptor.cell_number_vector[i_primal_face] = primal_face_number[i_primal_face];
+ for (size_t i = 0; i < diamond_number_of_cells; ++i) {
+ const auto [primal_face_id, dual_cell_id] = m_primal_face_to_dual_cell_map[i];
+ diamond_descriptor.cell_number_vector[dual_cell_id] = primal_face_number[primal_face_id];
}
if constexpr (Dimension == 3) {
@@ -60,7 +88,7 @@ DiamondDualConnectivityBuilder::_buildDiamondConnectivityDescriptor(const Connec
const auto& primal_face_to_cell_matrix = primal_connectivity.faceToCellMatrix();
- for (FaceId i_face = 0; i_face < primal_connectivity.numberOfFaces(); ++i_face) {
+ for (FaceId i_face = 0; i_face < primal_number_of_faces; ++i_face) {
const size_t i_cell = i_face;
const auto& primal_face_cell_list = primal_face_to_cell_matrix[i_face];
@@ -86,7 +114,7 @@ DiamondDualConnectivityBuilder::_buildDiamondConnectivityDescriptor(const Connec
const auto& primal_face_to_node_matrix = primal_connectivity.faceToNodeMatrix();
const auto& primal_face_local_number_in_their_cells = primal_connectivity.faceLocalNumbersInTheirCells();
const auto& cell_face_is_reversed = primal_connectivity.cellFaceIsReversed();
- for (FaceId i_face = 0; i_face < primal_connectivity.numberOfFaces(); ++i_face) {
+ for (FaceId i_face = 0; i_face < primal_number_of_faces; ++i_face) {
const size_t& i_diamond_cell = i_face;
const auto& primal_face_cell_list = primal_face_to_cell_matrix[i_face];
const auto& primal_face_node_list = primal_face_to_node_matrix[i_face];
@@ -162,23 +190,34 @@ DiamondDualConnectivityBuilder::_buildDiamondConnectivityDescriptor<1>(const Con
const size_t diamond_number_of_nodes = 2 * primal_number_of_nodes - primal_number_of_cells;
- diamond_descriptor.node_number_vector.resize(diamond_number_of_nodes);
-
- const auto& primal_node_number = primal_connectivity.nodeNumber();
-
const size_t diamond_number_of_cells = primal_number_of_faces;
const size_t number_of_kept_nodes = 2 * (diamond_number_of_nodes - diamond_number_of_cells);
const auto& primal_node_to_cell_matrix = primal_connectivity.nodeToCellMatrix();
size_t next_kept_node_id = 0;
- for (NodeId node_id = 0; node_id < primal_connectivity.numberOfNodes(); ++node_id) {
- const auto& primal_node_cell_list = primal_node_to_cell_matrix[node_id];
+ diamond_descriptor.node_number_vector.resize(diamond_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) {
- diamond_descriptor.node_number_vector[next_kept_node_id++] = primal_node_number[node_id];
+ diamond_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) {
+ diamond_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));
@@ -200,8 +239,6 @@ DiamondDualConnectivityBuilder::_buildDiamondConnectivityDescriptor<1>(const Con
diamond_descriptor.cell_to_node_vector.resize(diamond_number_of_cells);
- ;
-
const auto& primal_node_local_number_in_their_cells = primal_connectivity.nodeLocalNumbersInTheirCells();
{
@@ -278,7 +315,6 @@ DiamondDualConnectivityBuilder::_buildDiamondConnectivityFrom(const IConnectivit
node_array[i] = diamond_node_list[i];
}
diamond_descriptor.addRefItemList(RefNodeList{primal_ref_node_list.refId(), node_array});
- std::cout << "stored " << primal_ref_node_list.refId() << '\n';
}
}
}
@@ -466,15 +502,62 @@ DiamondDualConnectivityBuilder::_buildDiamondConnectivityFrom(const IConnectivit
diamond_descriptor.edge_owner_vector.resize(edge_cell_owner.size());
for (size_t i_edge = 0; i_edge < edge_cell_owner.size(); ++i_edge) {
- diamond_descriptor.face_owner_vector[i_edge] = diamond_descriptor.cell_owner_vector[edge_cell_owner[i_edge]];
+ diamond_descriptor.edge_owner_vector[i_edge] = diamond_descriptor.cell_owner_vector[edge_cell_owner[i_edge]];
}
}
m_connectivity = ConnectivityType::build(diamond_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<ConnectivityToDiamondDualConnectivityDataMapper<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);
}
DiamondDualConnectivityBuilder::DiamondDualConnectivityBuilder(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->_buildDiamondConnectivityFrom<1>(connectivity);
diff --git a/src/mesh/DiamondDualConnectivityBuilder.hpp b/src/mesh/DiamondDualConnectivityBuilder.hpp
index 670372df9e53d3bfa8dd3f1b03545d9264b41dc4..36b120811ad6478af90b89112740d272014b2e7a 100644
--- a/src/mesh/DiamondDualConnectivityBuilder.hpp
+++ b/src/mesh/DiamondDualConnectivityBuilder.hpp
@@ -2,6 +2,7 @@
#define DIAMOND_DUAL_CONNECTIVITY_BUILDER_HPP
#include <mesh/ConnectivityBuilderBase.hpp>
+#include <mesh/ItemIdToItemIdMap.hpp>
#include <memory>
@@ -9,17 +10,33 @@ template <size_t>
class Connectivity;
class ConnectivityDescriptor;
+class IConnectivityToDiamondDualConnectivityDataMapper;
+
class DiamondDualConnectivityBuilder : 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<IConnectivityToDiamondDualConnectivityDataMapper> m_mapper;
+
template <size_t Dimension>
void _buildDiamondConnectivityDescriptor(const Connectivity<Dimension>&, ConnectivityDescriptor&);
template <size_t Dimension>
void _buildDiamondConnectivityFrom(const IConnectivity&);
- public:
+ friend class DiamondDualConnectivityManager;
DiamondDualConnectivityBuilder(const IConnectivity&);
+
+ public:
+ std::shared_ptr<IConnectivityToDiamondDualConnectivityDataMapper>
+ mapper() const
+ {
+ return m_mapper;
+ }
+
~DiamondDualConnectivityBuilder() = default;
};
diff --git a/src/mesh/DiamondDualConnectivityManager.cpp b/src/mesh/DiamondDualConnectivityManager.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..9459010b69e89acaf3f247866d2f6e81ef00cd81
--- /dev/null
+++ b/src/mesh/DiamondDualConnectivityManager.cpp
@@ -0,0 +1,96 @@
+#include <utils/PugsAssert.hpp>
+
+#include <mesh/Connectivity.hpp>
+#include <mesh/ConnectivityToDiamondDualConnectivityDataMapper.hpp>
+#include <mesh/DiamondDualConnectivityBuilder.hpp>
+#include <mesh/DiamondDualConnectivityManager.hpp>
+#include <utils/Exceptions.hpp>
+
+#include <sstream>
+
+DiamondDualConnectivityManager* DiamondDualConnectivityManager::m_instance{nullptr};
+
+void
+DiamondDualConnectivityManager::create()
+{
+ Assert(m_instance == nullptr, "DiamondDualConnectivityManager is already created");
+ m_instance = new DiamondDualConnectivityManager;
+}
+
+void
+DiamondDualConnectivityManager::destroy()
+{
+ Assert(m_instance != nullptr, "DiamondDualConnectivityManager was not created!");
+
+ if (m_instance->m_connectivity_to_diamond_dual_connectivity_info_map.size() > 0) {
+ std::stringstream error;
+ error << ": some connectivities are still registered\n";
+ for (const auto& [connectivity, diamond_dual_connectivity_info] :
+ m_instance->m_connectivity_to_diamond_dual_connectivity_info_map) {
+ error << " - connectivity " << rang::fgB::magenta << connectivity << rang::style::reset << '\n';
+ }
+ throw UnexpectedError(error.str());
+ }
+ delete m_instance;
+ m_instance = nullptr;
+}
+
+void
+DiamondDualConnectivityManager::deleteConnectivity(const IConnectivity* p_connectivity)
+{
+ m_connectivity_to_diamond_dual_connectivity_info_map.erase(p_connectivity);
+}
+
+DiamondDualConnectivityManager::DiamondDualConnectivityInfo
+DiamondDualConnectivityManager::_getDiamondDualConnectivityInfo(const IConnectivity& connectivity)
+{
+ const IConnectivity* p_connectivity = &connectivity;
+
+ if (auto i_connectivity = m_connectivity_to_diamond_dual_connectivity_info_map.find(p_connectivity);
+ i_connectivity != m_connectivity_to_diamond_dual_connectivity_info_map.end()) {
+ return i_connectivity->second;
+ } else {
+ DiamondDualConnectivityBuilder builder{connectivity};
+
+ DiamondDualConnectivityInfo connectivity_info{builder.connectivity(), builder.mapper()};
+
+ m_connectivity_to_diamond_dual_connectivity_info_map[p_connectivity] = connectivity_info;
+
+ return connectivity_info;
+ }
+}
+
+template <size_t Dimension>
+std::shared_ptr<const Connectivity<Dimension>>
+DiamondDualConnectivityManager::getDiamondDualConnectivity(const Connectivity<Dimension>& connectivity)
+{
+ return std::dynamic_pointer_cast<const Connectivity<Dimension>>(
+ this->_getDiamondDualConnectivityInfo(connectivity).diamondDualConnectivity());
+}
+
+template <size_t Dimension>
+std::shared_ptr<const ConnectivityToDiamondDualConnectivityDataMapper<Dimension>>
+DiamondDualConnectivityManager::getConnectivityToDiamondDualConnectivityDataMapper(
+ const Connectivity<Dimension>& connectivity)
+{
+ return std::dynamic_pointer_cast<const ConnectivityToDiamondDualConnectivityDataMapper<Dimension>>(
+ this->_getDiamondDualConnectivityInfo(connectivity).connectivityToDiamondDualConnectivityDataMapper());
+}
+
+template std::shared_ptr<const Connectivity<1>> DiamondDualConnectivityManager::getDiamondDualConnectivity(
+ const Connectivity<1>& connectivity);
+
+template std::shared_ptr<const Connectivity<2>> DiamondDualConnectivityManager::getDiamondDualConnectivity(
+ const Connectivity<2>& connectivity);
+
+template std::shared_ptr<const Connectivity<3>> DiamondDualConnectivityManager::getDiamondDualConnectivity(
+ const Connectivity<3>& connectivity);
+
+template std::shared_ptr<const ConnectivityToDiamondDualConnectivityDataMapper<1>>
+DiamondDualConnectivityManager::getConnectivityToDiamondDualConnectivityDataMapper(const Connectivity<1>&);
+
+template std::shared_ptr<const ConnectivityToDiamondDualConnectivityDataMapper<2>>
+DiamondDualConnectivityManager::getConnectivityToDiamondDualConnectivityDataMapper(const Connectivity<2>&);
+
+template std::shared_ptr<const ConnectivityToDiamondDualConnectivityDataMapper<3>>
+DiamondDualConnectivityManager::getConnectivityToDiamondDualConnectivityDataMapper(const Connectivity<3>&);
diff --git a/src/mesh/DiamondDualConnectivityManager.hpp b/src/mesh/DiamondDualConnectivityManager.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..33dfbaf4cf7cc8e59eebfc095624d43c5bf038b3
--- /dev/null
+++ b/src/mesh/DiamondDualConnectivityManager.hpp
@@ -0,0 +1,94 @@
+#ifndef DIAMOND_DUAL_CONNECTIVITY_MANAGER_HPP
+#define DIAMOND_DUAL_CONNECTIVITY_MANAGER_HPP
+
+#include <mesh/IConnectivity.hpp>
+
+#include <memory>
+#include <unordered_map>
+
+template <size_t Dimension>
+class Connectivity;
+
+class IConnectivityToDiamondDualConnectivityDataMapper;
+
+template <size_t Dimension>
+class ConnectivityToDiamondDualConnectivityDataMapper;
+
+class DiamondDualConnectivityManager
+{
+ private:
+ class DiamondDualConnectivityInfo
+ {
+ private:
+ std::shared_ptr<const IConnectivity> m_diamond_dual_connectivity;
+ std::shared_ptr<const IConnectivityToDiamondDualConnectivityDataMapper>
+ m_connectivity_to_diamond_dual_connectivity_data_mapper;
+
+ public:
+ PUGS_INLINE
+ std::shared_ptr<const IConnectivity>
+ diamondDualConnectivity() const
+ {
+ return m_diamond_dual_connectivity;
+ }
+
+ PUGS_INLINE
+ std::shared_ptr<const IConnectivityToDiamondDualConnectivityDataMapper>
+ connectivityToDiamondDualConnectivityDataMapper()
+ {
+ return m_connectivity_to_diamond_dual_connectivity_data_mapper;
+ }
+
+ DiamondDualConnectivityInfo& operator=(const DiamondDualConnectivityInfo&) = default;
+ DiamondDualConnectivityInfo& operator=(DiamondDualConnectivityInfo&&) = default;
+
+ DiamondDualConnectivityInfo() = default;
+ DiamondDualConnectivityInfo(const DiamondDualConnectivityInfo&) = default;
+ DiamondDualConnectivityInfo(DiamondDualConnectivityInfo&&) = default;
+
+ DiamondDualConnectivityInfo(const std::shared_ptr<const IConnectivity>& diamond_dual_connectivity,
+ const std::shared_ptr<const IConnectivityToDiamondDualConnectivityDataMapper>&
+ connectivity_to_diamond_dual_connectivity_data_mapper)
+ : m_diamond_dual_connectivity{diamond_dual_connectivity},
+ m_connectivity_to_diamond_dual_connectivity_data_mapper{connectivity_to_diamond_dual_connectivity_data_mapper}
+ {}
+
+ ~DiamondDualConnectivityInfo() = default;
+ };
+
+ DiamondDualConnectivityInfo _getDiamondDualConnectivityInfo(const IConnectivity& connectivity);
+
+ std::unordered_map<const IConnectivity*, DiamondDualConnectivityInfo>
+ m_connectivity_to_diamond_dual_connectivity_info_map;
+
+ static DiamondDualConnectivityManager* m_instance;
+
+ DiamondDualConnectivityManager(const DiamondDualConnectivityManager&) = delete;
+ DiamondDualConnectivityManager(DiamondDualConnectivityManager&&) = delete;
+
+ DiamondDualConnectivityManager() = default;
+ ~DiamondDualConnectivityManager() = default;
+
+ public:
+ static void create();
+ static void destroy();
+
+ PUGS_INLINE
+ static DiamondDualConnectivityManager&
+ instance()
+ {
+ Assert(m_instance != nullptr, "DiamondDualConnectivityManager was not created!");
+ return *m_instance;
+ }
+
+ void deleteConnectivity(const IConnectivity*);
+
+ template <size_t Dimension>
+ std::shared_ptr<const Connectivity<Dimension>> getDiamondDualConnectivity(const Connectivity<Dimension>&);
+
+ template <size_t Dimension>
+ std::shared_ptr<const ConnectivityToDiamondDualConnectivityDataMapper<Dimension>>
+ getConnectivityToDiamondDualConnectivityDataMapper(const Connectivity<Dimension>&);
+};
+
+#endif // DIAMOND_DUAL_CONNECTIVITY_MANAGER_HPP
diff --git a/src/mesh/DiamondDualMeshBuilder.cpp b/src/mesh/DiamondDualMeshBuilder.cpp
index 2218541b338e196a913b3c0859b882ad668d252e..df4157c9c19b377f77b754d3f33ac54141539f2f 100644
--- a/src/mesh/DiamondDualMeshBuilder.cpp
+++ b/src/mesh/DiamondDualMeshBuilder.cpp
@@ -1,117 +1,60 @@
#include <mesh/DiamondDualMeshBuilder.hpp>
-#include <mesh/DiamondDualConnectivityBuilder.hpp>
-
#include <mesh/Connectivity.hpp>
-#include <mesh/ConnectivityDescriptor.hpp>
-#include <mesh/ConnectivityDispatcher.hpp>
+#include <mesh/ConnectivityToDiamondDualConnectivityDataMapper.hpp>
+#include <mesh/DiamondDualConnectivityBuilder.hpp>
+#include <mesh/DiamondDualConnectivityManager.hpp>
#include <mesh/ItemValueUtils.hpp>
#include <mesh/Mesh.hpp>
#include <mesh/MeshData.hpp>
#include <mesh/MeshDataManager.hpp>
-#include <mesh/RefId.hpp>
-#include <utils/Array.hpp>
-#include <utils/Messenger.hpp>
template <size_t Dimension>
void
-DiamondDualMeshBuilder::_buildDualDiamondMeshFrom(
- const Mesh<Connectivity<Dimension>>& primal_mesh,
- const std::shared_ptr<const Connectivity<Dimension>>& p_diamond_connectivity)
+DiamondDualMeshBuilder::_buildDualDiamondMeshFrom(const std::shared_ptr<const IMesh>& p_i_mesh)
{
using ConnectivityType = Connectivity<Dimension>;
using MeshType = Mesh<Connectivity<Dimension>>;
- const ConnectivityType& diamond_connectivity = *p_diamond_connectivity;
-
- NodeValue<TinyVector<Dimension>> diamond_xr{diamond_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();
-
- NodeId i_node = 0;
- for (; i_node < primal_mesh.numberOfNodes(); ++i_node) {
- diamond_xr[i_node] = primal_xr[i_node];
- }
+ std::shared_ptr p_primal_mesh = std::dynamic_pointer_cast<const MeshType>(p_i_mesh);
+ const MeshType& primal_mesh = *p_primal_mesh;
- for (CellId i_cell = 0; i_cell < primal_mesh.numberOfCells(); ++i_cell) {
- diamond_xr[i_node++] = primal_xj[i_cell];
- }
+ DiamondDualConnectivityManager& manager = DiamondDualConnectivityManager::instance();
- m_mesh = std::make_shared<MeshType>(p_diamond_connectivity, diamond_xr);
-}
-
-template <>
-void
-DiamondDualMeshBuilder::_buildDualDiamondMeshFrom(const Mesh<Connectivity<1>>& primal_mesh,
- const std::shared_ptr<const Connectivity<1>>& p_diamond_connectivity)
-{
- using ConnectivityType = Connectivity<1>;
- using MeshType = Mesh<Connectivity<1>>;
-
- const ConnectivityType& primal_connectivity = primal_mesh.connectivity();
- const auto& node_to_cell_matrix = primal_connectivity.nodeToCellMatrix();
+ std::shared_ptr<const ConnectivityType> p_diamond_connectivity =
+ manager.getDiamondDualConnectivity(primal_mesh.connectivity());
const ConnectivityType& diamond_connectivity = *p_diamond_connectivity;
- NodeValue<TinyVector<1>> diamond_xr{diamond_connectivity};
+ const NodeValue<const TinyVector<Dimension>> primal_xr = primal_mesh.xr();
- const NodeValue<const TinyVector<1>> primal_xr = primal_mesh.xr();
- MeshData<1>& primal_mesh_data = MeshDataManager::instance().getMeshData(primal_mesh);
- const CellValue<const TinyVector<1>> primal_xj = primal_mesh_data.xj();
+ MeshData<Dimension>& primal_mesh_data = MeshDataManager::instance().getMeshData(primal_mesh);
+ const CellValue<const TinyVector<Dimension>> primal_xj = primal_mesh_data.xj();
- NodeId next_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) {
- diamond_xr[next_node_id++] = primal_xr[primal_node_id];
- }
- }
+ std::shared_ptr connectivity_to_diamond_dual_connectivity_data_mapper =
+ manager.getConnectivityToDiamondDualConnectivityDataMapper(primal_mesh.connectivity());
- for (CellId i_cell = 0; i_cell < primal_mesh.numberOfCells(); ++i_cell) {
- diamond_xr[next_node_id++] = primal_xj[i_cell];
- }
+ NodeValue<TinyVector<Dimension>> diamond_xr{diamond_connectivity};
+ connectivity_to_diamond_dual_connectivity_data_mapper->toDualNode(primal_xr, primal_xj, diamond_xr);
m_mesh = std::make_shared<MeshType>(p_diamond_connectivity, diamond_xr);
}
DiamondDualMeshBuilder::DiamondDualMeshBuilder(const std::shared_ptr<const IMesh>& p_mesh)
{
+ std::cout << "building DiamondDualMesh\n";
+
switch (p_mesh->dimension()) {
case 1: {
- using ConnectivityType = Connectivity<1>;
- using MeshType = Mesh<ConnectivityType>;
-
- std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(p_mesh);
- DiamondDualConnectivityBuilder builder(mesh->connectivity());
-
- std::shared_ptr p_diamond_connectivity = std::dynamic_pointer_cast<const ConnectivityType>(builder.connectivity());
-
- this->_buildDualDiamondMeshFrom(*mesh, p_diamond_connectivity);
+ this->_buildDualDiamondMeshFrom<1>(p_mesh);
break;
}
case 2: {
- using ConnectivityType = Connectivity<2>;
- using MeshType = Mesh<ConnectivityType>;
-
- std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(p_mesh);
- DiamondDualConnectivityBuilder builder(mesh->connectivity());
-
- std::shared_ptr p_diamond_connectivity = std::dynamic_pointer_cast<const ConnectivityType>(builder.connectivity());
-
- this->_buildDualDiamondMeshFrom(*mesh, p_diamond_connectivity);
+ this->_buildDualDiamondMeshFrom<2>(p_mesh);
break;
}
case 3: {
- using ConnectivityType = Connectivity<3>;
- using MeshType = Mesh<ConnectivityType>;
-
- std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(p_mesh);
- DiamondDualConnectivityBuilder builder(mesh->connectivity());
-
- std::shared_ptr p_diamond_connectivity = std::dynamic_pointer_cast<const ConnectivityType>(builder.connectivity());
-
- this->_buildDualDiamondMeshFrom(*mesh, p_diamond_connectivity);
+ this->_buildDualDiamondMeshFrom<3>(p_mesh);
break;
}
default: {
diff --git a/src/mesh/DiamondDualMeshBuilder.hpp b/src/mesh/DiamondDualMeshBuilder.hpp
index ded3cfb0f40ca4d1b1bd9dc06f7981391f9401c6..2e47cfe21b4790a63e0e4aeba2fcfb9c59699286 100644
--- a/src/mesh/DiamondDualMeshBuilder.hpp
+++ b/src/mesh/DiamondDualMeshBuilder.hpp
@@ -5,21 +5,16 @@
#include <memory>
-template <size_t>
-class Connectivity;
-
-template <typename ConnectivityType>
-class Mesh;
-
class DiamondDualMeshBuilder : public MeshBuilderBase
{
private:
template <size_t Dimension>
- void _buildDualDiamondMeshFrom(const Mesh<Connectivity<Dimension>>&,
- const std::shared_ptr<const Connectivity<Dimension>>&);
+ void _buildDualDiamondMeshFrom(const std::shared_ptr<const IMesh>&);
- public:
+ friend class DiamondDualMeshManager;
DiamondDualMeshBuilder(const std::shared_ptr<const IMesh>&);
+
+ public:
~DiamondDualMeshBuilder() = default;
};
diff --git a/src/mesh/DiamondDualMeshManager.cpp b/src/mesh/DiamondDualMeshManager.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..0619736753f2b49a4229c60492ea6fd4c0de7212
--- /dev/null
+++ b/src/mesh/DiamondDualMeshManager.cpp
@@ -0,0 +1,65 @@
+#include <mesh/DiamondDualMeshManager.hpp>
+
+#include <mesh/Connectivity.hpp>
+#include <mesh/DiamondDualMeshBuilder.hpp>
+#include <mesh/Mesh.hpp>
+#include <utils/Exceptions.hpp>
+#include <utils/PugsAssert.hpp>
+
+#include <sstream>
+
+DiamondDualMeshManager* DiamondDualMeshManager::m_instance{nullptr};
+
+void
+DiamondDualMeshManager::create()
+{
+ Assert(m_instance == nullptr, "DiamondDualMeshManager is already created");
+ m_instance = new DiamondDualMeshManager;
+}
+
+void
+DiamondDualMeshManager::destroy()
+{
+ Assert(m_instance != nullptr, "DiamondDualMeshManager was not created!");
+
+ if (m_instance->m_mesh_to_diamond_dual_mesh_map.size() > 0) {
+ std::stringstream error;
+ error << ": some meshes are still registered\n";
+ for (const auto& i_mesh_data : m_instance->m_mesh_to_diamond_dual_mesh_map) {
+ error << " - mesh " << rang::fgB::magenta << i_mesh_data.first << rang::style::reset << '\n';
+ }
+ throw UnexpectedError(error.str());
+ }
+ delete m_instance;
+ m_instance = nullptr;
+}
+
+void
+DiamondDualMeshManager::deleteMesh(const IMesh* p_mesh)
+{
+ m_mesh_to_diamond_dual_mesh_map.erase(p_mesh);
+}
+
+template <size_t Dimension>
+std::shared_ptr<const Mesh<Connectivity<Dimension>>>
+DiamondDualMeshManager::getDiamondDualMesh(std::shared_ptr<const Mesh<Connectivity<Dimension>>> mesh)
+{
+ const IMesh* p_mesh = mesh.get();
+
+ if (auto i_mesh_data = m_mesh_to_diamond_dual_mesh_map.find(p_mesh);
+ i_mesh_data != m_mesh_to_diamond_dual_mesh_map.end()) {
+ return std::dynamic_pointer_cast<const Mesh<Connectivity<Dimension>>>(i_mesh_data->second);
+ } else {
+ DiamondDualMeshBuilder builder{mesh};
+
+ m_mesh_to_diamond_dual_mesh_map[p_mesh] = builder.mesh();
+ return std::dynamic_pointer_cast<const Mesh<Connectivity<Dimension>>>(builder.mesh());
+ }
+}
+
+template std::shared_ptr<const Mesh<Connectivity<1>>> DiamondDualMeshManager::getDiamondDualMesh(
+ std::shared_ptr<const Mesh<Connectivity<1>>>);
+template std::shared_ptr<const Mesh<Connectivity<2>>> DiamondDualMeshManager::getDiamondDualMesh(
+ std::shared_ptr<const Mesh<Connectivity<2>>>);
+template std::shared_ptr<const Mesh<Connectivity<3>>> DiamondDualMeshManager::getDiamondDualMesh(
+ std::shared_ptr<const Mesh<Connectivity<3>>>);
diff --git a/src/mesh/DiamondDualMeshManager.hpp b/src/mesh/DiamondDualMeshManager.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..9f802eb96d13704235f50e82ad75e84ba9be8a64
--- /dev/null
+++ b/src/mesh/DiamondDualMeshManager.hpp
@@ -0,0 +1,49 @@
+#ifndef DIAMOND_DUAL_MESH_MANAGER_HPP
+#define DIAMOND_DUAL_MESH_MANAGER_HPP
+
+#include <mesh/IMesh.hpp>
+#include <utils/PugsAssert.hpp>
+#include <utils/PugsMacros.hpp>
+
+#include <memory>
+#include <unordered_map>
+
+template <size_t Dimension>
+class Connectivity;
+
+template <typename ConnectivityType>
+class Mesh;
+
+class DiamondDualMeshManager
+{
+ private:
+ std::unordered_map<const IMesh*, std::shared_ptr<const IMesh>> m_mesh_to_diamond_dual_mesh_map;
+
+ static DiamondDualMeshManager* m_instance;
+
+ DiamondDualMeshManager(const DiamondDualMeshManager&) = delete;
+ DiamondDualMeshManager(DiamondDualMeshManager&&) = delete;
+
+ DiamondDualMeshManager() = default;
+ ~DiamondDualMeshManager() = default;
+
+ public:
+ static void create();
+ static void destroy();
+
+ PUGS_INLINE
+ static DiamondDualMeshManager&
+ instance()
+ {
+ Assert(m_instance != nullptr, "DiamondDualMeshManager was not created!");
+ return *m_instance;
+ }
+
+ void deleteMesh(const IMesh*);
+
+ template <size_t Dimension>
+ std::shared_ptr<const Mesh<Connectivity<Dimension>>> getDiamondDualMesh(
+ std::shared_ptr<const Mesh<Connectivity<Dimension>>>);
+};
+
+#endif // DIAMOND_DUAL_MESH_MANAGER_HPP
diff --git a/src/mesh/IConnectivity.cpp b/src/mesh/IConnectivity.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..48cf57d55ebceb7818f755fd9214807a247f6575
--- /dev/null
+++ b/src/mesh/IConnectivity.cpp
@@ -0,0 +1,10 @@
+#include <mesh/IConnectivity.hpp>
+
+#include <mesh/DiamondDualConnectivityManager.hpp>
+#include <mesh/SynchronizerManager.hpp>
+
+IConnectivity::~IConnectivity()
+{
+ SynchronizerManager::instance().deleteConnectivitySynchronizer(this);
+ DiamondDualConnectivityManager::instance().deleteConnectivity(this);
+}
diff --git a/src/mesh/IConnectivity.hpp b/src/mesh/IConnectivity.hpp
index 27be707cbcc4976b17e9db1e7bf6f37195c9d63a..9bfb8e5043088a18880440235cd72b1677c1f519 100644
--- a/src/mesh/IConnectivity.hpp
+++ b/src/mesh/IConnectivity.hpp
@@ -35,7 +35,7 @@ class IConnectivity : public std::enable_shared_from_this<IConnectivity>
IConnectivity() = default;
IConnectivity(IConnectivity&&) = delete;
IConnectivity(const IConnectivity&) = delete;
- ~IConnectivity() = default;
+ virtual ~IConnectivity();
};
template <>
diff --git a/src/mesh/IMesh.cpp b/src/mesh/IMesh.cpp
index cb2accc141ed09c778237587373494797fee726d..f315f4b74b01eb371c42eb4bf029f81e93ae6016 100644
--- a/src/mesh/IMesh.cpp
+++ b/src/mesh/IMesh.cpp
@@ -1,8 +1,10 @@
#include <mesh/IMesh.hpp>
+#include <mesh/DiamondDualMeshManager.hpp>
#include <mesh/MeshDataManager.hpp>
IMesh::~IMesh()
{
- MeshDataManager::instance().deleteMeshData(*this);
+ MeshDataManager::instance().deleteMeshData(this);
+ DiamondDualMeshManager::instance().deleteMesh(this);
}
diff --git a/src/mesh/ItemIdToItemIdMap.hpp b/src/mesh/ItemIdToItemIdMap.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..0e92eb13ce3943c4e2487081ee2e48fc3deac8ee
--- /dev/null
+++ b/src/mesh/ItemIdToItemIdMap.hpp
@@ -0,0 +1,30 @@
+#ifndef ITEM_ID_TO_ITEM_ID_MAP_HPP
+#define ITEM_ID_TO_ITEM_ID_MAP_HPP
+
+#include <mesh/ItemId.hpp>
+#include <utils/Array.hpp>
+
+template <ItemType type1, ItemType type2>
+using ItemIdToItemIdMap = Array<std::pair<ItemIdT<type1>, ItemIdT<type2>>>;
+
+using NodeIdToNodeIdMap = ItemIdToItemIdMap<ItemType::node, ItemType::node>;
+using NodeIdToEdgeIdMap = ItemIdToItemIdMap<ItemType::node, ItemType::edge>;
+using NodeIdToFaceIdMap = ItemIdToItemIdMap<ItemType::node, ItemType::face>;
+using NodeIdToCellIdMap = ItemIdToItemIdMap<ItemType::node, ItemType::cell>;
+
+using EdgeIdToNodeIdMap = ItemIdToItemIdMap<ItemType::edge, ItemType::node>;
+using EdgeIdToEdgeIdMap = ItemIdToItemIdMap<ItemType::edge, ItemType::edge>;
+using EdgeIdToFaceIdMap = ItemIdToItemIdMap<ItemType::edge, ItemType::face>;
+using EdgeIdToCellIdMap = ItemIdToItemIdMap<ItemType::edge, ItemType::cell>;
+
+using FaceIdToNodeIdMap = ItemIdToItemIdMap<ItemType::face, ItemType::node>;
+using FaceIdToEdgeIdMap = ItemIdToItemIdMap<ItemType::face, ItemType::edge>;
+using FaceIdToFaceIdMap = ItemIdToItemIdMap<ItemType::face, ItemType::face>;
+using FaceIdToCellIdMap = ItemIdToItemIdMap<ItemType::face, ItemType::cell>;
+
+using CellIdToNodeIdMap = ItemIdToItemIdMap<ItemType::cell, ItemType::node>;
+using CellIdToEdgeIdMap = ItemIdToItemIdMap<ItemType::cell, ItemType::edge>;
+using CellIdToFaceIdMap = ItemIdToItemIdMap<ItemType::cell, ItemType::face>;
+using CellIdToCellIdMap = ItemIdToItemIdMap<ItemType::cell, ItemType::cell>;
+
+#endif // ITEM_ID_TO_ITEM_ID_MAP_HPP
diff --git a/src/mesh/MeshDataManager.cpp b/src/mesh/MeshDataManager.cpp
index 281e7517aa704e6c4352769e88f854b86cbc7366..8ca184a6013c083f5e8df7412d19f80b8f50f95a 100644
--- a/src/mesh/MeshDataManager.cpp
+++ b/src/mesh/MeshDataManager.cpp
@@ -35,9 +35,9 @@ MeshDataManager::destroy()
}
void
-MeshDataManager::deleteMeshData(const IMesh& mesh)
+MeshDataManager::deleteMeshData(const IMesh* p_mesh)
{
- m_mesh_mesh_data_map.erase(&mesh);
+ m_mesh_mesh_data_map.erase(p_mesh);
}
template <size_t Dimension>
diff --git a/src/mesh/MeshDataManager.hpp b/src/mesh/MeshDataManager.hpp
index 7b127148b889bc585082d626b79c435293fa4cd7..fd04101dea29a2d15e0feb2193f22f785ff65dc3 100644
--- a/src/mesh/MeshDataManager.hpp
+++ b/src/mesh/MeshDataManager.hpp
@@ -5,8 +5,8 @@
#include <utils/PugsAssert.hpp>
#include <utils/PugsMacros.hpp>
-#include <map>
#include <memory>
+#include <unordered_map>
class IMesh;
@@ -22,7 +22,7 @@ class MeshData;
class MeshDataManager
{
private:
- std::map<const IMesh*, std::shared_ptr<IMeshData>> m_mesh_mesh_data_map;
+ std::unordered_map<const IMesh*, std::shared_ptr<IMeshData>> m_mesh_mesh_data_map;
static MeshDataManager* m_instance;
@@ -44,7 +44,7 @@ class MeshDataManager
return *m_instance;
}
- void deleteMeshData(const IMesh&);
+ void deleteMeshData(const IMesh*);
template <size_t Dimension>
MeshData<Dimension>& getMeshData(const Mesh<Connectivity<Dimension>>&);
diff --git a/src/mesh/SynchronizerManager.hpp b/src/mesh/SynchronizerManager.hpp
index 79cb3db0a3473bdbe9f88d4d20836a835aca9260..8bcb441f05c4e841a71c103f9a252bb52c59dcb2 100644
--- a/src/mesh/SynchronizerManager.hpp
+++ b/src/mesh/SynchronizerManager.hpp
@@ -4,8 +4,8 @@
#include <utils/PugsAssert.hpp>
#include <utils/PugsMacros.hpp>
-#include <map>
#include <memory>
+#include <unordered_map>
class IConnectivity;
class Synchronizer;
@@ -13,7 +13,7 @@ class Synchronizer;
class SynchronizerManager
{
private:
- std::map<const IConnectivity*, std::shared_ptr<Synchronizer>> m_connectivity_synchronizer_map;
+ std::unordered_map<const IConnectivity*, std::shared_ptr<Synchronizer>> m_connectivity_synchronizer_map;
static SynchronizerManager* m_instance;
diff --git a/src/output/VTKWriter.hpp b/src/output/VTKWriter.hpp
index 5d85b3630e2afc77ef7f2a3422b2f7d207abe5b7..d2459353e4bffa282eac8ff27e21c63d93273664 100644
--- a/src/output/VTKWriter.hpp
+++ b/src/output/VTKWriter.hpp
@@ -8,12 +8,14 @@
#include <output/OutputNamedItemValueSet.hpp>
#include <utils/Exceptions.hpp>
#include <utils/Messenger.hpp>
+#include <utils/PugsAssert.hpp>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <sstream>
#include <string>
+#include <unordered_map>
class VTKWriter
{
@@ -234,6 +236,10 @@ class VTKWriter
std::visit([&, name = name](auto&& item_value) { return this->_write_cell_pvtu(fout, name, item_value); },
item_value_variant);
}
+ if constexpr (MeshType::Dimension == 3) {
+ fout << "<PDataArray type=\"Int64\" IdType=\"1\" Name=\"faces\"/>\n";
+ fout << "<PDataArray type=\"Int64\" IdType=\"1\" Name=\"faceoffsets\"/>\n";
+ }
fout << "</PCells>\n";
fout << "<PPointData>\n";
@@ -347,6 +353,10 @@ class VTKWriter
types[j] = 12;
break;
}
+ case CellType::Diamond: {
+ types[j] = 41;
+ break;
+ }
default: {
std::ostringstream os;
os << __FILE__ << ':' << __LINE__ << ": unknown cell type";
@@ -356,7 +366,62 @@ class VTKWriter
});
_write_array(fout, "types", types);
}
+ if constexpr (MeshType::Dimension == 3) {
+ const auto& cell_to_face_matrix = mesh->connectivity().cellToFaceMatrix();
+ const auto& cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix();
+ const auto& face_to_node_matrix = mesh->connectivity().faceToNodeMatrix();
+ const auto& cell_face_is_reversed = mesh->connectivity().cellFaceIsReversed();
+
+ Array<size_t> faces_offsets(mesh->numberOfCells());
+ size_t next_offset = 0;
+ fout << "<DataArray type=\"Int64\" IdType=\"1\" Name=\"faces\">\n";
+ for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+ const auto& cell_nodes = cell_to_node_matrix[cell_id];
+ std::unordered_map<size_t, size_t> node_id_to_local_number_in_cell;
+ for (size_t i_cell_node = 0; i_cell_node < cell_nodes.size(); ++i_cell_node) {
+ node_id_to_local_number_in_cell[cell_nodes[i_cell_node]] = i_cell_node;
+ }
+
+ const auto& cell_faces = cell_to_face_matrix[cell_id];
+ fout << cell_faces.size() << '\n';
+ next_offset++;
+ for (size_t i_cell_face = 0; i_cell_face < cell_faces.size(); ++i_cell_face) {
+ const FaceId& face_id = cell_faces[i_cell_face];
+ const auto& face_nodes = face_to_node_matrix[face_id];
+ fout << face_nodes.size();
+ next_offset++;
+ Array<size_t> face_node_in_cell(face_nodes.size());
+
+ for (size_t i_face_node = 0; i_face_node < face_nodes.size(); ++i_face_node) {
+ const NodeId& node_id = face_nodes[i_face_node];
+ auto i_node_id_to_local_number_in_cell = node_id_to_local_number_in_cell.find(node_id);
+ Assert(i_node_id_to_local_number_in_cell != node_id_to_local_number_in_cell.end());
+ face_node_in_cell[i_face_node] = i_node_id_to_local_number_in_cell->second;
+ }
+
+ if (cell_face_is_reversed(cell_id, i_cell_face)) {
+ for (size_t i = 0; i < face_node_in_cell.size(); ++i) {
+ fout << ' ' << face_node_in_cell[face_node_in_cell.size() - 1 - i];
+ }
+ } else {
+ for (size_t i = 0; i < face_node_in_cell.size(); ++i) {
+ fout << ' ' << face_node_in_cell[i];
+ }
+ }
+
+ next_offset += face_nodes.size();
+ fout << '\n';
+ }
+ faces_offsets[cell_id] = next_offset;
+ }
+ fout << "</DataArray>\n";
+ fout << "<DataArray type=\"Int64\" IdType=\"1\" Name=\"faceoffsets\">\n";
+ for (size_t i_face_offsets = 0; i_face_offsets < faces_offsets.size(); ++i_face_offsets) {
+ fout << faces_offsets[i_face_offsets] << '\n';
+ }
+ fout << "</DataArray>\n";
+ }
fout << "</Cells>\n";
fout << "</Piece>\n";
fout << "</UnstructuredGrid>\n";