Skip to content
Snippets Groups Projects
Commit b9fe8d4f authored by Stéphane Del Pino's avatar Stéphane Del Pino
Browse files

Add dual diamond mesh in 1d

- code needs a lot of clean-up
parent 183b788d
No related branches found
No related tags found
1 merge request!37Feature/language
...@@ -16,12 +16,41 @@ DiamondDualMeshBuilder::_buildDiamondConnectivityDescriptor(const Connectivity<D ...@@ -16,12 +16,41 @@ DiamondDualMeshBuilder::_buildDiamondConnectivityDescriptor(const Connectivity<D
ConnectivityDescriptor& diamond_descriptor) ConnectivityDescriptor& diamond_descriptor)
{ {
const size_t primal_number_of_nodes = primal_connectivity.numberOfNodes(); 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 primal_number_of_cells = primal_connectivity.numberOfCells();
const size_t diamond_number_of_nodes = primal_number_of_cells + primal_number_of_nodes; const size_t diamond_number_of_nodes = [&]() {
if constexpr (Dimension == 1) {
return 2 * primal_number_of_nodes - primal_number_of_cells;
} else {
return primal_number_of_cells + primal_number_of_nodes;
}
}();
diamond_descriptor.node_number_vector.resize(diamond_number_of_nodes); diamond_descriptor.node_number_vector.resize(diamond_number_of_nodes);
if constexpr (Dimension == 1) {
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];
if (primal_node_cell_list.size() == 1) {
diamond_descriptor.node_number_vector[next_kept_node_id++] = primal_node_number[node_id];
}
}
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));
}
} else {
const auto& primal_node_number = primal_connectivity.nodeNumber(); const auto& primal_node_number = primal_connectivity.nodeNumber();
for (NodeId primal_node_id = 0; primal_node_id < primal_connectivity.numberOfNodes(); ++primal_node_id) { for (NodeId primal_node_id = 0; primal_node_id < primal_connectivity.numberOfNodes(); ++primal_node_id) {
...@@ -36,16 +65,21 @@ DiamondDualMeshBuilder::_buildDiamondConnectivityDescriptor(const Connectivity<D ...@@ -36,16 +65,21 @@ DiamondDualMeshBuilder::_buildDiamondConnectivityDescriptor(const Connectivity<D
diamond_descriptor.node_number_vector[primal_number_of_nodes + primal_cell_id] = diamond_descriptor.node_number_vector[primal_number_of_nodes + primal_cell_id] =
primal_cell_number[primal_cell_id] + max_node_number; primal_cell_number[primal_cell_id] + max_node_number;
} }
}
const size_t diamond_number_of_cells = primal_connectivity.numberOfFaces(); const size_t diamond_number_of_cells = primal_number_of_faces;
diamond_descriptor.cell_number_vector.resize(diamond_number_of_cells); diamond_descriptor.cell_number_vector.resize(diamond_number_of_cells);
if constexpr (Dimension == 1) {
const auto& primal_node_number = primal_connectivity.nodeNumber();
for (NodeId primal_node_id = 0; primal_node_id < primal_number_of_nodes; ++primal_node_id) {
diamond_descriptor.cell_number_vector[primal_node_id] = primal_node_number[primal_node_id];
}
} else {
const auto& primal_face_number = primal_connectivity.faceNumber(); const auto& primal_face_number = primal_connectivity.faceNumber();
for (FaceId i_primal_face = 0; i_primal_face < primal_number_of_faces; ++i_primal_face) {
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_number_vector[i_primal_face] = primal_face_number[i_primal_face];
} }
}
if constexpr (Dimension == 3) { if constexpr (Dimension == 3) {
const size_t number_of_edges = diamond_descriptor.edge_to_node_vector.size(); const size_t number_of_edges = diamond_descriptor.edge_to_node_vector.size();
diamond_descriptor.edge_number_vector.resize(number_of_edges); diamond_descriptor.edge_number_vector.resize(number_of_edges);
...@@ -57,8 +91,6 @@ DiamondDualMeshBuilder::_buildDiamondConnectivityDescriptor(const Connectivity<D ...@@ -57,8 +91,6 @@ DiamondDualMeshBuilder::_buildDiamondConnectivityDescriptor(const Connectivity<D
} }
} }
diamond_descriptor.cell_to_node_vector.resize(diamond_number_of_cells);
diamond_descriptor.cell_type_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(); const auto& primal_face_to_cell_matrix = primal_connectivity.faceToCellMatrix();
...@@ -68,7 +100,7 @@ DiamondDualMeshBuilder::_buildDiamondConnectivityDescriptor(const Connectivity<D ...@@ -68,7 +100,7 @@ DiamondDualMeshBuilder::_buildDiamondConnectivityDescriptor(const Connectivity<D
const auto& primal_face_cell_list = primal_face_to_cell_matrix[i_face]; const auto& primal_face_cell_list = primal_face_to_cell_matrix[i_face];
if constexpr (Dimension == 1) { if constexpr (Dimension == 1) {
throw NotImplementedError("dimension 1"); diamond_descriptor.cell_type_vector[i_cell] = CellType::Line;
} else if constexpr (Dimension == 2) { } else if constexpr (Dimension == 2) {
if (primal_face_cell_list.size() == 1) { if (primal_face_cell_list.size() == 1) {
diamond_descriptor.cell_type_vector[i_cell] = CellType::Triangle; diamond_descriptor.cell_type_vector[i_cell] = CellType::Triangle;
...@@ -90,6 +122,32 @@ DiamondDualMeshBuilder::_buildDiamondConnectivityDescriptor(const Connectivity<D ...@@ -90,6 +122,32 @@ DiamondDualMeshBuilder::_buildDiamondConnectivityDescriptor(const Connectivity<D
diamond_descriptor.cell_to_node_vector.resize(diamond_number_of_cells); diamond_descriptor.cell_to_node_vector.resize(diamond_number_of_cells);
if constexpr (Dimension == 1) {
const auto& primal_node_to_cell_matrix = primal_connectivity.nodeToCellMatrix();
const size_t number_of_kept_nodes = 2 * (diamond_number_of_nodes - diamond_number_of_cells);
size_t next_kept_node_id = 0;
const auto& primal_node_local_number_in_their_cells = primal_connectivity.nodeLocalNumbersInTheirCells();
for (NodeId i_node = 0; i_node < primal_connectivity.numberOfNodes(); ++i_node) {
const size_t& i_diamond_cell = i_node;
const auto& primal_node_cell_list = primal_node_to_cell_matrix[i_node];
diamond_descriptor.cell_to_node_vector[i_diamond_cell].resize(2);
if (primal_node_cell_list.size() == 1) {
const CellId cell_id = primal_node_cell_list[0];
const auto i_node_in_cell = primal_node_local_number_in_their_cells(i_node, 0);
diamond_descriptor.cell_to_node_vector[i_diamond_cell][i_node_in_cell] = next_kept_node_id++;
diamond_descriptor.cell_to_node_vector[i_diamond_cell][1 - i_node_in_cell] =
number_of_kept_nodes + primal_node_cell_list[0];
} else {
diamond_descriptor.cell_to_node_vector[i_diamond_cell][0] = number_of_kept_nodes + primal_node_cell_list[0];
diamond_descriptor.cell_to_node_vector[i_diamond_cell][1] = number_of_kept_nodes + primal_node_cell_list[1];
}
}
} else {
static_assert(Dimension > 1, "invalid dimension");
const auto& primal_face_to_node_matrix = primal_connectivity.faceToNodeMatrix(); const auto& primal_face_to_node_matrix = primal_connectivity.faceToNodeMatrix();
const auto& primal_face_local_number_in_their_cells = primal_connectivity.faceLocalNumbersInTheirCells(); const auto& primal_face_local_number_in_their_cells = primal_connectivity.faceLocalNumbersInTheirCells();
const auto& cell_face_is_reversed = primal_connectivity.cellFaceIsReversed(); const auto& cell_face_is_reversed = primal_connectivity.cellFaceIsReversed();
...@@ -128,7 +186,7 @@ DiamondDualMeshBuilder::_buildDiamondConnectivityDescriptor(const Connectivity<D ...@@ -128,7 +186,7 @@ DiamondDualMeshBuilder::_buildDiamondConnectivityDescriptor(const Connectivity<D
const CellId cell0_id = primal_face_cell_list[0]; const CellId cell0_id = primal_face_cell_list[0];
const CellId cell1_id = primal_face_cell_list[1]; 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); const auto i_face_in_cell0 = primal_face_local_number_in_their_cells(i_face, 0);
if constexpr (Dimension == 2) { if constexpr (Dimension == 2) {
Assert(primal_face_node_list.size() == 2); Assert(primal_face_node_list.size() == 2);
...@@ -137,7 +195,7 @@ DiamondDualMeshBuilder::_buildDiamondConnectivityDescriptor(const Connectivity<D ...@@ -137,7 +195,7 @@ DiamondDualMeshBuilder::_buildDiamondConnectivityDescriptor(const Connectivity<D
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][2] = primal_number_of_nodes + cell1_id;
diamond_descriptor.cell_to_node_vector[i_diamond_cell][3] = primal_face_node_list[1]; 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)) { if (cell_face_is_reversed(cell0_id, i_face_in_cell0)) {
std::swap(diamond_descriptor.cell_to_node_vector[i_diamond_cell][1], std::swap(diamond_descriptor.cell_to_node_vector[i_diamond_cell][1],
diamond_descriptor.cell_to_node_vector[i_diamond_cell][3]); diamond_descriptor.cell_to_node_vector[i_diamond_cell][3]);
} }
...@@ -149,7 +207,7 @@ DiamondDualMeshBuilder::_buildDiamondConnectivityDescriptor(const Connectivity<D ...@@ -149,7 +207,7 @@ DiamondDualMeshBuilder::_buildDiamondConnectivityDescriptor(const Connectivity<D
diamond_descriptor.cell_to_node_vector[i_diamond_cell][primal_face_node_list.size() + 1] = diamond_descriptor.cell_to_node_vector[i_diamond_cell][primal_face_node_list.size() + 1] =
primal_number_of_nodes + cell1_id; primal_number_of_nodes + cell1_id;
if (cell_face_is_reversed(cell0_id, i_face_in_cell)) { if (cell_face_is_reversed(cell0_id, i_face_in_cell0)) {
std::swap(diamond_descriptor.cell_to_node_vector[i_diamond_cell][0], 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]); diamond_descriptor.cell_to_node_vector[i_diamond_cell][primal_face_node_list.size() + 1]);
} }
...@@ -157,6 +215,7 @@ DiamondDualMeshBuilder::_buildDiamondConnectivityDescriptor(const Connectivity<D ...@@ -157,6 +215,7 @@ DiamondDualMeshBuilder::_buildDiamondConnectivityDescriptor(const Connectivity<D
} }
} }
} }
}
template <size_t Dimension> template <size_t Dimension>
void void
...@@ -174,10 +233,12 @@ DiamondDualMeshBuilder::_buildDiamondMeshFrom(const std::shared_ptr<const IMesh> ...@@ -174,10 +233,12 @@ DiamondDualMeshBuilder::_buildDiamondMeshFrom(const std::shared_ptr<const IMesh>
this->_buildDiamondConnectivityDescriptor(primal_connectivity, diamond_descriptor); this->_buildDiamondConnectivityDescriptor(primal_connectivity, diamond_descriptor);
if constexpr (Dimension > 1) {
MeshBuilderBase::_computeCellFaceAndFaceNodeConnectivities<Dimension>(diamond_descriptor); MeshBuilderBase::_computeCellFaceAndFaceNodeConnectivities<Dimension>(diamond_descriptor);
if constexpr (Dimension == 3) { if constexpr (Dimension > 2) {
MeshBuilderBase::_computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities<Dimension>(diamond_descriptor); MeshBuilderBase::_computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities<Dimension>(diamond_descriptor);
} }
}
{ {
const std::unordered_map<unsigned int, NodeId> node_to_id_map = [&] { const std::unordered_map<unsigned int, NodeId> node_to_id_map = [&] {
...@@ -329,6 +390,21 @@ DiamondDualMeshBuilder::_buildDiamondMeshFrom(const std::shared_ptr<const IMesh> ...@@ -329,6 +390,21 @@ DiamondDualMeshBuilder::_buildDiamondMeshFrom(const std::shared_ptr<const IMesh>
diamond_descriptor.node_owner_vector.resize(diamond_descriptor.node_number_vector.size()); diamond_descriptor.node_owner_vector.resize(diamond_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) {
diamond_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) {
diamond_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(); const auto& primal_node_owner = primal_connectivity.nodeOwner();
for (NodeId primal_node_id = 0; primal_node_id < primal_connectivity.numberOfNodes(); ++primal_node_id) { 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]; diamond_descriptor.node_owner_vector[primal_node_id] = primal_node_owner[primal_node_id];
...@@ -337,12 +413,22 @@ DiamondDualMeshBuilder::_buildDiamondMeshFrom(const std::shared_ptr<const IMesh> ...@@ -337,12 +413,22 @@ DiamondDualMeshBuilder::_buildDiamondMeshFrom(const std::shared_ptr<const IMesh>
for (CellId primal_cell_id = 0; primal_cell_id < primal_number_of_cells; ++primal_cell_id) { 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.node_owner_vector[primal_number_of_nodes + primal_cell_id] = primal_cell_owner[primal_cell_id];
} }
}
if constexpr (Dimension == 1) {
diamond_descriptor.cell_owner_vector.resize(diamond_descriptor.cell_number_vector.size()); diamond_descriptor.cell_owner_vector.resize(diamond_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) {
diamond_descriptor.cell_owner_vector[primal_node_id] = primal_node_owner[primal_node_id];
}
} else {
diamond_descriptor.cell_owner_vector.resize(diamond_descriptor.cell_number_vector.size());
const size_t primal_number_of_faces = primal_connectivity.numberOfFaces();
const auto& primal_face_owner = primal_connectivity.faceOwner(); const auto& primal_face_owner = primal_connectivity.faceOwner();
for (FaceId primal_face_id = 0; primal_face_id < primal_number_of_cells; ++primal_face_id) { for (FaceId primal_face_id = 0; primal_face_id < primal_number_of_faces; ++primal_face_id) {
diamond_descriptor.cell_owner_vector[primal_face_id] = primal_face_owner[primal_face_id]; diamond_descriptor.cell_owner_vector[primal_face_id] = primal_face_owner[primal_face_id];
} }
}
{ {
std::vector<int> face_cell_owner(diamond_descriptor.face_number_vector.size()); std::vector<int> face_cell_owner(diamond_descriptor.face_number_vector.size());
...@@ -391,6 +477,20 @@ DiamondDualMeshBuilder::_buildDiamondMeshFrom(const std::shared_ptr<const IMesh> ...@@ -391,6 +477,20 @@ DiamondDualMeshBuilder::_buildDiamondMeshFrom(const std::shared_ptr<const IMesh>
{ {
#warning define transfer functions #warning define transfer functions
if constexpr (Dimension == 1) {
const auto& node_to_cell_matrix = primal_connectivity.nodeToCellMatrix();
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];
}
}
for (CellId i_cell = 0; i_cell < primal_number_of_cells; ++i_cell) {
diamond_xr[next_node_id++] = primal_xj[i_cell];
}
} else {
NodeId i_node = 0; NodeId i_node = 0;
for (; i_node < primal_number_of_nodes; ++i_node) { for (; i_node < primal_number_of_nodes; ++i_node) {
diamond_xr[i_node] = primal_xr[i_node]; diamond_xr[i_node] = primal_xr[i_node];
...@@ -400,6 +500,7 @@ DiamondDualMeshBuilder::_buildDiamondMeshFrom(const std::shared_ptr<const IMesh> ...@@ -400,6 +500,7 @@ DiamondDualMeshBuilder::_buildDiamondMeshFrom(const std::shared_ptr<const IMesh>
diamond_xr[i_node++] = primal_xj[i_cell]; diamond_xr[i_node++] = primal_xj[i_cell];
} }
} }
}
std::shared_ptr p_diamond_mesh = std::make_shared<MeshType>(p_diamond_connectivity, diamond_xr); std::shared_ptr p_diamond_mesh = std::make_shared<MeshType>(p_diamond_connectivity, diamond_xr);
...@@ -418,13 +519,6 @@ DiamondDualMeshBuilder::_buildDiamondMeshFrom(const std::shared_ptr<const IMesh> ...@@ -418,13 +519,6 @@ DiamondDualMeshBuilder::_buildDiamondMeshFrom(const std::shared_ptr<const IMesh>
m_mesh = std::make_shared<MeshType>(p_diamond_connectivity, diamond_xr); m_mesh = std::make_shared<MeshType>(p_diamond_connectivity, diamond_xr);
} }
template <>
[[deprecated]] void
DiamondDualMeshBuilder::_buildDiamondMeshFrom<1>(const std::shared_ptr<const IMesh>&)
{
m_mesh = 0;
}
DiamondDualMeshBuilder::DiamondDualMeshBuilder(const std::shared_ptr<const IMesh>& p_mesh) DiamondDualMeshBuilder::DiamondDualMeshBuilder(const std::shared_ptr<const IMesh>& p_mesh)
{ {
switch (p_mesh->dimension()) { switch (p_mesh->dimension()) {
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment