From 39ef20c01defe09e8e6afc56e74eb09897b115e8 Mon Sep 17 00:00:00 2001
From: Stephane Del Pino <stephane.delpino44@gmail.com>
Date: Fri, 26 Jun 2020 19:20:08 +0200
Subject: [PATCH] Continue Diamond mesh construction
- ItemValueUtils: min/max/sum functions require improvement or fixes
- Diamond meshes are built correctly in 2D and sequential
- faces boundaries references lists are correctly translated
- node boundaries (ref lists) are missing
- 1D and 3D are to be done
---
src/mesh/DiamondDualMeshBuilder.cpp | 220 +++++++++++++++++++++++++++-
src/mesh/DiamondDualMeshBuilder.hpp | 4 +-
src/mesh/ItemValueUtils.hpp | 15 +-
3 files changed, 227 insertions(+), 12 deletions(-)
diff --git a/src/mesh/DiamondDualMeshBuilder.cpp b/src/mesh/DiamondDualMeshBuilder.cpp
index 2c5108775..480693993 100644
--- a/src/mesh/DiamondDualMeshBuilder.cpp
+++ b/src/mesh/DiamondDualMeshBuilder.cpp
@@ -3,16 +3,232 @@
#include <mesh/Connectivity.hpp>
#include <mesh/ConnectivityDescriptor.hpp>
#include <mesh/ConnectivityDispatcher.hpp>
+#include <mesh/ItemValueUtils.hpp>
+#include <mesh/Mesh.hpp>
+#include <mesh/MeshData.hpp>
#include <mesh/RefId.hpp>
#include <utils/Array.hpp>
#include <utils/Messenger.hpp>
template <size_t Dimension>
void
-DiamondDualMeshBuilder::_buildDiamondMeshFrom(const std::shared_ptr<const IMesh>& p_mesh)
+DiamondDualMeshBuilder::_buildDiamondMeshFrom(const std::shared_ptr<const IMesh>&)
{
- static_assert(Dimension <= 3, "invalid dimension");
+ 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();
+
+ const size_t diamond_number_of_nodes = primal_number_of_cells + primal_number_of_nodes;
+
+ 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();
+
+ // 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];
+ // }
+
+#warning fix max and compute diamond node numbers correctly in parallel
+ // 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;
+ }
+
+ const size_t diamond_number_of_cells = primal_connectivity.numberOfFaces();
+ 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_connectivity.numberOfFaces(); ++i_primal_face) {
+ diamond_descriptor.cell_number_vector[i_primal_face] = primal_face_number[i_primal_face];
+ }
+
+ diamond_descriptor.cell_to_node_vector.resize(diamond_number_of_cells);
+
+ diamond_descriptor.cell_type_vector.resize(diamond_number_of_cells);
+
+ const auto& primal_face_to_cell_matrix = primal_connectivity.faceToCellMatrix();
+
+ for (FaceId i_face = 0; i_face < primal_connectivity.numberOfFaces(); ++i_face) {
+ 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;
+ } else {
+ Assert(primal_face_cell_list.size() == 2);
+ diamond_descriptor.cell_type_vector[i_cell] = CellType::Quadrangle;
+ }
+ }
+
+ diamond_descriptor.cell_to_node_vector.resize(diamond_number_of_cells);
+
+ 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) {
+ 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];
+ if (primal_face_cell_list.size() == 1) {
+ diamond_descriptor.cell_to_node_vector[i_diamond_cell].resize(3);
+
+ 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;
+
+ 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]);
+ }
+ } else {
+ Assert(primal_face_cell_list.size() == 2);
+ diamond_descriptor.cell_to_node_vector[i_diamond_cell].resize(4);
+
+ 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 (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]);
+ }
+ }
+ }
+
+ MeshBuilderBase::_computeCellFaceAndFaceNodeConnectivities<2>(diamond_descriptor);
+
+ {
+ using Face = ConnectivityFace<2>;
+
+ 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 < diamond_descriptor.face_to_node_vector.size(); ++l) {
+ const auto& node_vector = diamond_descriptor.face_to_node_vector[l];
+
+ face_to_id_map[Face(node_vector, diamond_descriptor.node_number_vector)] = l;
+ }
+ return face_to_id_map;
+ }();
+
+ for (size_t i_face_list = 0; i_face_list < primal_connectivity.numberOfRefItemList<ItemType::face>();
+ ++i_face_list) {
+ const auto& primal_ref_face_list = primal_connectivity.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();
+
+ const std::vector<FaceId> diamond_face_list = [&]() {
+ std::vector<FaceId> diamond_face_list;
+ diamond_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_diamond_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, diamond_descriptor.node_number_vector));
+ }();
+
+ if (i_diamond_face != face_to_id_map.end()) {
+ diamond_face_list.push_back(i_diamond_face->second);
+ }
+ }
+ return diamond_face_list;
+ }();
+
+ if (parallel::allReduceOr(diamond_face_list.size() > 0)) {
+ Array<FaceId> face_array(diamond_face_list.size());
+ for (size_t i = 0; i < diamond_face_list.size(); ++i) {
+ face_array[i] = diamond_face_list[i];
+ }
+ diamond_descriptor.addRefItemList(RefFaceList{primal_ref_face_list.refId(), face_array});
+ }
+ }
+ }
+
+#warning completly wrong
+ 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());
+
+ 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());
+
+ std::shared_ptr p_diamond_connectivity = Connectivity2D::build(diamond_descriptor);
+ Connectivity2D& diamond_connectivity = *p_diamond_connectivity;
+
+ NodeValue<TinyVector<2>> diamond_xr{diamond_connectivity};
+
+ const auto primal_xr = primal_mesh.xr();
+ MeshData<Mesh<Connectivity2D>> primal_mesh_data{primal_mesh};
+ const auto primal_xj = primal_mesh_data.xj();
+
+ {
+#warning to recode
+ NodeId i_node = 0;
+ for (; i_node < primal_number_of_nodes; ++i_node) {
+ diamond_xr[i_node] = primal_xr[i_node];
+ }
+
+ for (CellId i_cell = 0; i_cell < primal_number_of_cells; ++i_cell) {
+ diamond_xr[i_node++] = primal_xj[i_cell];
+ }
+ }
+
+ std::shared_ptr p_diamond_mesh = std::make_shared<Mesh<Connectivity2D>>(p_diamond_connectivity, diamond_xr);
+
+ MeshData<Mesh<Connectivity2D>> dual_mesh_data{*p_diamond_mesh};
+
+ double sum = 0;
+ for (CellId cell_id = 0; cell_id < p_diamond_mesh->numberOfCells(); ++cell_id) {
+ sum += dual_mesh_data.Vj()[cell_id];
+ }
+
+ std::cout << "volume = " << sum << '\n';
+
+ m_mesh = std::make_shared<Mesh<Connectivity2D>>(p_diamond_connectivity, diamond_xr);
+}
+
+template <>
+void
+DiamondDualMeshBuilder::_buildDiamondMeshFrom<3>(const std::shared_ptr<const IMesh>& p_mesh)
+{
m_mesh = p_mesh;
}
diff --git a/src/mesh/DiamondDualMeshBuilder.hpp b/src/mesh/DiamondDualMeshBuilder.hpp
index 16b2a107d..7efbafa24 100644
--- a/src/mesh/DiamondDualMeshBuilder.hpp
+++ b/src/mesh/DiamondDualMeshBuilder.hpp
@@ -10,10 +10,10 @@ class DiamondDualMeshBuilder : public MeshBuilderBase
{
private:
template <size_t Dimension>
- void _buildDiamondMeshFrom(const std::shared_ptr<const IMesh>& mesh);
+ void _buildDiamondMeshFrom(const std::shared_ptr<const IMesh>&);
public:
- DiamondDualMeshBuilder(const std::shared_ptr<const IMesh>& mesh);
+ DiamondDualMeshBuilder(const std::shared_ptr<const IMesh>&);
~DiamondDualMeshBuilder() = default;
};
diff --git a/src/mesh/ItemValueUtils.hpp b/src/mesh/ItemValueUtils.hpp
index ca0c38ecd..6e73e0abe 100644
--- a/src/mesh/ItemValueUtils.hpp
+++ b/src/mesh/ItemValueUtils.hpp
@@ -10,11 +10,11 @@
#include <iostream>
-template <typename DataType, ItemType item_type>
+template <typename DataType, ItemType item_type, typename ConnectivityPtr>
std::remove_const_t<DataType>
-min(const ItemValue<DataType, item_type>& item_value)
+min(const ItemValue<DataType, item_type, ConnectivityPtr>& item_value)
{
- using ItemValueType = ItemValue<DataType, item_type>;
+ using ItemValueType = ItemValue<DataType, item_type, ConnectivityPtr>;
using data_type = std::remove_const_t<typename ItemValueType::data_type>;
using index_type = typename ItemValueType::index_type;
@@ -100,11 +100,11 @@ min(const ItemValue<DataType, item_type>& item_value)
return parallel::allReduceMin(local_min);
}
-template <typename DataType, ItemType item_type>
+template <typename DataType, ItemType item_type, typename ConnectivityPtr>
std::remove_const_t<DataType>
-max(const ItemValue<DataType, item_type>& item_value)
+max(const ItemValue<DataType, item_type, ConnectivityPtr>& item_value)
{
- using ItemValueType = ItemValue<DataType, item_type>;
+ using ItemValueType = ItemValue<DataType, item_type, ConnectivityPtr>;
using data_type = std::remove_const_t<typename ItemValueType::data_type>;
using index_type = typename ItemValueType::index_type;
@@ -268,8 +268,7 @@ sum(const ItemValue<DataType, item_type>& item_value)
}
PUGS_INLINE
- ItemValueSum(const ItemValueType& item_value, const IsOwnedType& is_owned)
- : m_item_value(item_value), m_is_owned(is_owned)
+ ItemValueSum(const ItemValueType& item_value) : m_item_value(item_value), m_is_owned{is_owned(item_value)}
{
;
}
--
GitLab