diff --git a/src/mesh/DiamondDualMeshBuilder.cpp b/src/mesh/DiamondDualMeshBuilder.cpp
index faa7f27793f6a654700f7c3f143a500a816a4a98..52857868c67382d0f8f766ec243a2c39aeb95146 100644
--- a/src/mesh/DiamondDualMeshBuilder.cpp
+++ b/src/mesh/DiamondDualMeshBuilder.cpp
@@ -164,6 +164,9 @@ DiamondDualMeshBuilder::_buildDiamondMeshFrom(const std::shared_ptr<const IMesh>
this->_buildDiamondConnectivityDescriptor(primal_connectivity, diamond_descriptor);
MeshBuilderBase::_computeCellFaceAndFaceNodeConnectivities<Dimension>(diamond_descriptor);
+ if constexpr (Dimension == 3) {
+ MeshBuilderBase::_computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities<Dimension>(diamond_descriptor);
+ }
{
const auto& primal_face_to_node_matrix = primal_connectivity.faceToNodeMatrix();
@@ -257,6 +260,33 @@ DiamondDualMeshBuilder::_buildDiamondMeshFrom(const std::shared_ptr<const IMesh>
}
}
+ if constexpr (Dimension == 3) {
+ const size_t number_of_edges = diamond_descriptor.edge_to_node_vector.size();
+ diamond_descriptor.edge_number_vector.resize(number_of_edges);
+ for (size_t i_edge = 0; i_edge < number_of_edges; ++i_edge) {
+ diamond_descriptor.edge_number_vector[i_edge] = i_edge;
+ }
+ if (parallel::size() > 1) {
+ throw NotImplementedError("parallel edge numbering is undefined");
+ }
+
+ std::vector<int> edge_cell_owner(number_of_edges);
+ std::fill(std::begin(edge_cell_owner), std::end(edge_cell_owner), parallel::size());
+
+ for (size_t i_cell = 0; i_cell < diamond_descriptor.cell_to_face_vector.size(); ++i_cell) {
+ const auto& cell_edge_list = diamond_descriptor.cell_to_edge_vector[i_cell];
+ for (size_t i_edge = 0; i_edge < cell_edge_list.size(); ++i_edge) {
+ const size_t edge_id = cell_edge_list[i_edge];
+ edge_cell_owner[edge_id] = std::min(edge_cell_owner[edge_id], diamond_descriptor.cell_number_vector[i_cell]);
+ }
+ }
+
+ 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]];
+ }
+ }
+
std::shared_ptr p_diamond_connectivity = ConnectivityType::build(diamond_descriptor);
ConnectivityType& diamond_connectivity = *p_diamond_connectivity;
diff --git a/src/mesh/MeshBuilderBase.cpp b/src/mesh/MeshBuilderBase.cpp
index 33236a9e2619d546c3bd6db1e7dbd307afa5da23..7f7221b049872fac2d78b936093be706f03ea080 100644
--- a/src/mesh/MeshBuilderBase.cpp
+++ b/src/mesh/MeshBuilderBase.cpp
@@ -321,7 +321,7 @@ MeshBuilderBase::_computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities(Connectivi
}
{
- descriptor.cell_to_node_vector.reserve(descriptor.cell_to_node_vector.size());
+ descriptor.cell_to_edge_vector.reserve(descriptor.cell_to_node_vector.size());
for (CellId j = 0; j < descriptor.cell_to_node_vector.size(); ++j) {
const auto& cell_nodes = descriptor.cell_to_node_vector[j];
@@ -359,10 +359,78 @@ MeshBuilderBase::_computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities(Connectivi
descriptor.cell_to_edge_vector.emplace_back(cell_edge_vector);
break;
}
+ case CellType::Pyramid: {
+ const size_t number_of_edges = 2 * cell_nodes.size();
+ std::vector<unsigned int> base_nodes;
+ std::copy_n(cell_nodes.begin(), cell_nodes.size() - 1, std::back_inserter(base_nodes));
+
+ std::vector<unsigned int> cell_edge_vector;
+ cell_edge_vector.reserve(number_of_edges);
+ for (size_t i_edge = 0; i_edge < base_nodes.size(); ++i_edge) {
+ Edge edge{{base_nodes[i_edge], base_nodes[(i_edge + 1) % base_nodes.size()]}, node_number_vector};
+ auto i = edge_id_map.find(edge);
+ if (i == edge_id_map.end()) {
+ throw NormalError("could not find this edge");
+ }
+ cell_edge_vector.push_back(i->second);
+ }
+
+ const unsigned int top_vertex = cell_nodes[cell_nodes.size() - 1];
+ for (size_t i_edge = 0; i_edge < base_nodes.size(); ++i_edge) {
+ Edge edge{{base_nodes[i_edge], top_vertex}, node_number_vector};
+ auto i = edge_id_map.find(edge);
+ if (i == edge_id_map.end()) {
+ throw NormalError("could not find this edge");
+ }
+ cell_edge_vector.push_back(i->second);
+ }
+ descriptor.cell_to_edge_vector.emplace_back(cell_edge_vector);
+ break;
+ }
+ case CellType::Diamond: {
+ const size_t number_of_edges = 3 * cell_nodes.size();
+ std::vector<unsigned int> base_nodes;
+ std::copy_n(cell_nodes.begin() + 1, cell_nodes.size() - 2, std::back_inserter(base_nodes));
+
+ std::vector<unsigned int> cell_edge_vector;
+ cell_edge_vector.reserve(number_of_edges);
+ for (size_t i_edge = 0; i_edge < base_nodes.size(); ++i_edge) {
+ Edge edge{{base_nodes[i_edge], base_nodes[(i_edge + 1) % base_nodes.size()]}, node_number_vector};
+ auto i = edge_id_map.find(edge);
+ if (i == edge_id_map.end()) {
+ throw NormalError("could not find this edge");
+ }
+ cell_edge_vector.push_back(i->second);
+ }
+
+ const unsigned int top_vertex = cell_nodes[cell_nodes.size() - 1];
+ for (size_t i_edge = 0; i_edge < base_nodes.size(); ++i_edge) {
+ Edge edge{{base_nodes[i_edge], top_vertex}, node_number_vector};
+ auto i = edge_id_map.find(edge);
+ if (i == edge_id_map.end()) {
+ throw NormalError("could not find this edge");
+ }
+ cell_edge_vector.push_back(i->second);
+ }
+
+ const unsigned int bottom_vertex = cell_nodes[0];
+ for (size_t i_edge = 0; i_edge < base_nodes.size(); ++i_edge) {
+ Edge edge{{base_nodes[i_edge], bottom_vertex}, node_number_vector};
+ auto i = edge_id_map.find(edge);
+ if (i == edge_id_map.end()) {
+ throw NormalError("could not find this edge");
+ }
+ cell_edge_vector.push_back(i->second);
+ }
+
+ descriptor.cell_to_edge_vector.emplace_back(cell_edge_vector);
+
+ break;
+ }
default: {
std::stringstream error_msg;
error_msg << name(descriptor.cell_type_vector[j]) << ": unexpected cell type in dimension 3";
- throw UnexpectedError(error_msg.str());
+ throw NotImplementedError(error_msg.str());
}
}
}