From f2be671f59790e80b5d759415fb657c7bebbeb75 Mon Sep 17 00:00:00 2001
From: Stephane Del Pino <stephane.delpino44@gmail.com>
Date: Mon, 29 Jun 2020 17:34:58 +0200
Subject: [PATCH] Generalize diamond dual mesh to 3D
---
src/mesh/DiamondDualMeshBuilder.cpp | 210 ++++++++++++++++++----------
src/mesh/DiamondDualMeshBuilder.hpp | 8 ++
2 files changed, 148 insertions(+), 70 deletions(-)
diff --git a/src/mesh/DiamondDualMeshBuilder.cpp b/src/mesh/DiamondDualMeshBuilder.cpp
index 480693993..faa7f2779 100644
--- a/src/mesh/DiamondDualMeshBuilder.cpp
+++ b/src/mesh/DiamondDualMeshBuilder.cpp
@@ -12,32 +12,9 @@
template <size_t Dimension>
void
-DiamondDualMeshBuilder::_buildDiamondMeshFrom(const std::shared_ptr<const IMesh>&)
+DiamondDualMeshBuilder::_buildDiamondConnectivityDescriptor(const Connectivity<Dimension>& primal_connectivity,
+ ConnectivityDescriptor& diamond_descriptor)
{
- static_assert(Dimension <= 3, "invalid mesh dimension");
-}
-
-template <>
-void
-DiamondDualMeshBuilder::_buildDiamondMeshFrom<1>(const std::shared_ptr<const IMesh>& p_mesh)
-{
- m_mesh = p_mesh;
-}
-
-template <>
-void
-DiamondDualMeshBuilder::_buildDiamondMeshFrom<2>(const std::shared_ptr<const IMesh>& p_mesh)
-{
- m_mesh = p_mesh;
-
- using MeshType = Mesh<Connectivity2D>;
-
- const MeshType& primal_mesh = dynamic_cast<const MeshType&>(*p_mesh);
-
- const Connectivity2D& primal_connectivity = primal_mesh.connectivity();
-
- ConnectivityDescriptor diamond_descriptor;
-
const size_t primal_number_of_nodes = primal_connectivity.numberOfNodes();
const size_t primal_number_of_cells = primal_connectivity.numberOfCells();
@@ -45,18 +22,19 @@ DiamondDualMeshBuilder::_buildDiamondMeshFrom<2>(const std::shared_ptr<const IMe
diamond_descriptor.node_number_vector.resize(diamond_number_of_nodes);
- // const auto& primal_node_number = primal_connectivity.nodeNumber();
- // const auto& primal_cell_number = primal_connectivity.cellNumber();
+ const auto& primal_node_number = primal_connectivity.nodeNumber();
+
+ 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];
+ }
- // for (size_t i_primal_node = 0; i_primal_node < primal_connectivity.numberOfNodes(); ++i_primal_node) {
- // diamond_descriptor.node_number_vector[i_primal_node] = primal_node_number[i_primal_node];
- // }
+ const auto& primal_cell_number = primal_connectivity.cellNumber();
-#warning fix max and compute diamond node numbers correctly in parallel
- // const size_t max_node_number = max(primal_node_number);
+ const size_t max_node_number = max(primal_node_number);
- for (size_t i = 0; i < diamond_number_of_nodes; ++i) {
- diamond_descriptor.node_number_vector[i] = i;
+ 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;
}
const size_t diamond_number_of_cells = primal_connectivity.numberOfFaces();
@@ -78,11 +56,24 @@ DiamondDualMeshBuilder::_buildDiamondMeshFrom<2>(const std::shared_ptr<const IMe
const size_t i_cell = i_face;
const auto& primal_face_cell_list = primal_face_to_cell_matrix[i_face];
- if (primal_face_cell_list.size() == 1) {
- diamond_descriptor.cell_type_vector[i_cell] = CellType::Triangle;
+ if constexpr (Dimension == 1) {
+ throw NotImplementedError("dimension 1");
+ } else if constexpr (Dimension == 2) {
+ if (primal_face_cell_list.size() == 1) {
+ diamond_descriptor.cell_type_vector[i_cell] = CellType::Triangle;
+ } else {
+ Assert(primal_face_cell_list.size() == 2);
+ diamond_descriptor.cell_type_vector[i_cell] = CellType::Quadrangle;
+ }
} else {
- Assert(primal_face_cell_list.size() == 2);
- diamond_descriptor.cell_type_vector[i_cell] = CellType::Quadrangle;
+ static_assert(Dimension == 3, "unexpected dimension");
+
+ if (primal_face_cell_list.size() == 1) {
+ diamond_descriptor.cell_type_vector[i_cell] = CellType::Pyramid;
+ } else {
+ Assert(primal_face_cell_list.size() == 2);
+ diamond_descriptor.cell_type_vector[i_cell] = CellType::Diamond;
+ }
}
}
@@ -96,43 +87,88 @@ DiamondDualMeshBuilder::_buildDiamondMeshFrom<2>(const std::shared_ptr<const IMe
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];
if (primal_face_cell_list.size() == 1) {
- diamond_descriptor.cell_to_node_vector[i_diamond_cell].resize(3);
+ diamond_descriptor.cell_to_node_vector[i_diamond_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(i_face, 0);
- diamond_descriptor.cell_to_node_vector[i_diamond_cell][0] = primal_face_node_list[0];
- diamond_descriptor.cell_to_node_vector[i_diamond_cell][1] = primal_face_node_list[1];
- diamond_descriptor.cell_to_node_vector[i_diamond_cell][2] = primal_number_of_nodes + cell_id;
+ for (size_t i_node = 0; i_node < primal_face_node_list.size(); ++i_node) {
+ diamond_descriptor.cell_to_node_vector[i_diamond_cell][i_node] = primal_face_node_list[i_node];
+ }
+ diamond_descriptor.cell_to_node_vector[i_diamond_cell][primal_face_node_list.size()] =
+ primal_number_of_nodes + cell_id;
if (cell_face_is_reversed(cell_id, i_face_in_cell)) {
- std::swap(diamond_descriptor.cell_to_node_vector[i_diamond_cell][0],
- diamond_descriptor.cell_to_node_vector[i_diamond_cell][1]);
+ if constexpr (Dimension == 2) {
+ std::swap(diamond_descriptor.cell_to_node_vector[i_diamond_cell][0],
+ diamond_descriptor.cell_to_node_vector[i_diamond_cell][1]);
+
+ } else {
+ for (size_t i_node = 0; i_node < primal_face_node_list.size() / 2; ++i_node) {
+ std::swap(diamond_descriptor.cell_to_node_vector[i_diamond_cell][i_node],
+ diamond_descriptor
+ .cell_to_node_vector[i_diamond_cell][primal_face_node_list.size() - 1 - i_node]);
+ }
+ }
}
} else {
Assert(primal_face_cell_list.size() == 2);
- diamond_descriptor.cell_to_node_vector[i_diamond_cell].resize(4);
+ diamond_descriptor.cell_to_node_vector[i_diamond_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_cell = primal_face_local_number_in_their_cells(i_face, 0);
- diamond_descriptor.cell_to_node_vector[i_diamond_cell][0] = primal_number_of_nodes + cell0_id;
- diamond_descriptor.cell_to_node_vector[i_diamond_cell][1] = primal_face_node_list[0];
- diamond_descriptor.cell_to_node_vector[i_diamond_cell][2] = primal_number_of_nodes + cell1_id;
- diamond_descriptor.cell_to_node_vector[i_diamond_cell][3] = primal_face_node_list[1];
+ if constexpr (Dimension == 2) {
+ Assert(primal_face_node_list.size() == 2);
+ diamond_descriptor.cell_to_node_vector[i_diamond_cell][0] = primal_number_of_nodes + cell0_id;
+ diamond_descriptor.cell_to_node_vector[i_diamond_cell][1] = primal_face_node_list[0];
+ diamond_descriptor.cell_to_node_vector[i_diamond_cell][2] = primal_number_of_nodes + cell1_id;
+ diamond_descriptor.cell_to_node_vector[i_diamond_cell][3] = primal_face_node_list[1];
+
+ if (cell_face_is_reversed(cell0_id, i_face_in_cell)) {
+ std::swap(diamond_descriptor.cell_to_node_vector[i_diamond_cell][1],
+ diamond_descriptor.cell_to_node_vector[i_diamond_cell][3]);
+ }
+ } else {
+ diamond_descriptor.cell_to_node_vector[i_diamond_cell][0] = primal_number_of_nodes + cell0_id;
+ for (size_t i_node = 0; i_node < primal_face_node_list.size(); ++i_node) {
+ diamond_descriptor.cell_to_node_vector[i_diamond_cell][i_node + 1] = primal_face_node_list[i_node];
+ }
+ diamond_descriptor.cell_to_node_vector[i_diamond_cell][primal_face_node_list.size() + 1] =
+ primal_number_of_nodes + cell1_id;
- if (cell_face_is_reversed(cell0_id, i_face_in_cell)) {
- std::swap(diamond_descriptor.cell_to_node_vector[i_diamond_cell][1],
- diamond_descriptor.cell_to_node_vector[i_diamond_cell][3]);
+ if (cell_face_is_reversed(cell0_id, i_face_in_cell)) {
+ std::swap(diamond_descriptor.cell_to_node_vector[i_diamond_cell][0],
+ diamond_descriptor.cell_to_node_vector[i_diamond_cell][primal_face_node_list.size() + 1]);
+ }
}
}
}
+}
- MeshBuilderBase::_computeCellFaceAndFaceNodeConnectivities<2>(diamond_descriptor);
+template <size_t Dimension>
+void
+DiamondDualMeshBuilder::_buildDiamondMeshFrom(const std::shared_ptr<const IMesh>& p_mesh)
+{
+ static_assert(Dimension <= 3, "invalid mesh dimension");
+ using ConnectivityType = Connectivity<Dimension>;
+ using MeshType = Mesh<ConnectivityType>;
+
+ const MeshType& primal_mesh = dynamic_cast<const MeshType&>(*p_mesh);
+
+ const ConnectivityType& primal_connectivity = primal_mesh.connectivity();
+
+ ConnectivityDescriptor diamond_descriptor;
+
+ this->_buildDiamondConnectivityDescriptor(primal_connectivity, diamond_descriptor);
+
+ MeshBuilderBase::_computeCellFaceAndFaceNodeConnectivities<Dimension>(diamond_descriptor);
{
- using Face = ConnectivityFace<2>;
+ 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;
@@ -144,9 +180,9 @@ DiamondDualMeshBuilder::_buildDiamondMeshFrom<2>(const std::shared_ptr<const IMe
return face_to_id_map;
}();
- for (size_t i_face_list = 0; i_face_list < primal_connectivity.numberOfRefItemList<ItemType::face>();
+ 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.refItemList<ItemType::face>(i_face_list);
+ const auto& primal_ref_face_list = primal_connectivity.template refItemList<ItemType::face>(i_face_list);
std::cout << "treating " << primal_ref_face_list.refId() << '\n';
const auto& primal_face_list = primal_ref_face_list.list();
@@ -183,24 +219,55 @@ DiamondDualMeshBuilder::_buildDiamondMeshFrom<2>(const std::shared_ptr<const IMe
}
}
-#warning completly wrong
+ const size_t primal_number_of_nodes = primal_connectivity.numberOfNodes();
+ const size_t primal_number_of_cells = primal_connectivity.numberOfCells();
+
diamond_descriptor.node_owner_vector.resize(diamond_descriptor.node_number_vector.size());
- std::fill(diamond_descriptor.node_owner_vector.begin(), diamond_descriptor.node_owner_vector.end(), parallel::rank());
+
+ const auto& primal_node_owner = primal_connectivity.nodeOwner();
+ for (NodeId primal_node_id = 0; primal_node_id < primal_connectivity.numberOfNodes(); ++primal_node_id) {
+ diamond_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) {
+ diamond_descriptor.node_owner_vector[primal_number_of_nodes + primal_cell_id] = primal_cell_owner[primal_cell_id];
+ }
diamond_descriptor.cell_owner_vector.resize(diamond_descriptor.cell_number_vector.size());
- std::fill(diamond_descriptor.cell_owner_vector.begin(), diamond_descriptor.cell_owner_vector.end(), parallel::rank());
+ const auto& primal_face_owner = primal_connectivity.faceOwner();
+ for (FaceId primal_face_id = 0; primal_face_id < primal_number_of_cells; ++primal_face_id) {
+ diamond_descriptor.cell_owner_vector[primal_face_id] = primal_face_owner[primal_face_id];
+ }
- std::shared_ptr p_diamond_connectivity = Connectivity2D::build(diamond_descriptor);
- Connectivity2D& diamond_connectivity = *p_diamond_connectivity;
+ {
+ std::vector<int> face_cell_owner(diamond_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 < diamond_descriptor.cell_to_face_vector.size(); ++i_cell) {
+ const auto& cell_face_list = diamond_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], diamond_descriptor.cell_number_vector[i_cell]);
+ }
+ }
- NodeValue<TinyVector<2>> diamond_xr{diamond_connectivity};
+ diamond_descriptor.face_owner_vector.resize(face_cell_owner.size());
+ for (size_t i_face = 0; i_face < face_cell_owner.size(); ++i_face) {
+ diamond_descriptor.face_owner_vector[i_face] = diamond_descriptor.cell_owner_vector[face_cell_owner[i_face]];
+ }
+ }
+
+ std::shared_ptr p_diamond_connectivity = ConnectivityType::build(diamond_descriptor);
+ ConnectivityType& diamond_connectivity = *p_diamond_connectivity;
+
+ NodeValue<TinyVector<Dimension>> diamond_xr{diamond_connectivity};
const auto primal_xr = primal_mesh.xr();
- MeshData<Mesh<Connectivity2D>> primal_mesh_data{primal_mesh};
+ MeshData<MeshType> primal_mesh_data{primal_mesh};
const auto primal_xj = primal_mesh_data.xj();
{
-#warning to recode
+#warning define transfer functions
NodeId i_node = 0;
for (; i_node < primal_number_of_nodes; ++i_node) {
diamond_xr[i_node] = primal_xr[i_node];
@@ -211,9 +278,11 @@ DiamondDualMeshBuilder::_buildDiamondMeshFrom<2>(const std::shared_ptr<const IMe
}
}
- std::shared_ptr p_diamond_mesh = std::make_shared<Mesh<Connectivity2D>>(p_diamond_connectivity, diamond_xr);
+ std::shared_ptr p_diamond_mesh = std::make_shared<MeshType>(p_diamond_connectivity, diamond_xr);
- MeshData<Mesh<Connectivity2D>> dual_mesh_data{*p_diamond_mesh};
+#warning USELESS TEST
+ // -->>
+ MeshData<MeshType> dual_mesh_data{*p_diamond_mesh};
double sum = 0;
for (CellId cell_id = 0; cell_id < p_diamond_mesh->numberOfCells(); ++cell_id) {
@@ -221,15 +290,16 @@ DiamondDualMeshBuilder::_buildDiamondMeshFrom<2>(const std::shared_ptr<const IMe
}
std::cout << "volume = " << sum << '\n';
+ // <<--
- m_mesh = std::make_shared<Mesh<Connectivity2D>>(p_diamond_connectivity, diamond_xr);
+ m_mesh = std::make_shared<MeshType>(p_diamond_connectivity, diamond_xr);
}
template <>
-void
-DiamondDualMeshBuilder::_buildDiamondMeshFrom<3>(const std::shared_ptr<const IMesh>& p_mesh)
+[[deprecated]] void
+DiamondDualMeshBuilder::_buildDiamondMeshFrom<1>(const std::shared_ptr<const IMesh>&)
{
- m_mesh = p_mesh;
+ m_mesh = 0;
}
DiamondDualMeshBuilder::DiamondDualMeshBuilder(const std::shared_ptr<const IMesh>& p_mesh)
diff --git a/src/mesh/DiamondDualMeshBuilder.hpp b/src/mesh/DiamondDualMeshBuilder.hpp
index 7efbafa24..a6af1b830 100644
--- a/src/mesh/DiamondDualMeshBuilder.hpp
+++ b/src/mesh/DiamondDualMeshBuilder.hpp
@@ -6,9 +6,17 @@
#include <memory>
+template <size_t>
+class Connectivity;
+
+class ConnectivityDescriptor;
+
class DiamondDualMeshBuilder : public MeshBuilderBase
{
private:
+ template <size_t Dimension>
+ void _buildDiamondConnectivityDescriptor(const Connectivity<Dimension>&, ConnectivityDescriptor&);
+
template <size_t Dimension>
void _buildDiamondMeshFrom(const std::shared_ptr<const IMesh>&);
--
GitLab