diff --git a/src/mesh/DiamondDualConnectivityBuilder.cpp b/src/mesh/DiamondDualConnectivityBuilder.cpp
index 76adb60de4b047cd5d456662db96d10c6e1677e7..31d54be8997e9571958399f17c44dbcede1c35e4 100644
--- a/src/mesh/DiamondDualConnectivityBuilder.cpp
+++ b/src/mesh/DiamondDualConnectivityBuilder.cpp
@@ -44,19 +44,19 @@ DiamondDualConnectivityBuilder::_buildDiamondConnectivityDescriptor(const Connec
diamond_descriptor.node_number_vector.resize(diamond_number_of_nodes);
- for (size_t i = 0; i < m_primal_node_to_dual_node_map.size(); ++i) {
+ parallel_for(m_primal_node_to_dual_node_map.size(), [&](size_t i) {
const auto [primal_node_id, diamond_dual_node_id] = m_primal_node_to_dual_node_map[i];
diamond_descriptor.node_number_vector[diamond_dual_node_id] = primal_node_number[primal_node_id];
- }
+ });
const size_t cell_number_shift = max(primal_node_number) + 1;
- for (size_t i = 0; i < primal_number_of_cells; ++i) {
+ parallel_for(primal_number_of_cells, [&](size_t i) {
const auto [primal_cell_id, diamond_dual_node_id] = m_primal_cell_to_dual_node_map[i];
diamond_descriptor.node_number_vector[diamond_dual_node_id] =
primal_cell_number[primal_cell_id] + cell_number_shift;
- }
+ });
{
m_primal_face_to_dual_cell_map = FaceIdToCellIdMap{primal_number_of_faces};
@@ -68,10 +68,10 @@ DiamondDualConnectivityBuilder::_buildDiamondConnectivityDescriptor(const Connec
diamond_descriptor.cell_number_vector.resize(diamond_number_of_cells);
const auto& primal_face_number = primal_connectivity.faceNumber();
- for (size_t i = 0; i < diamond_number_of_cells; ++i) {
+ parallel_for(diamond_number_of_cells, [&](size_t 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) {
const size_t number_of_edges = diamond_descriptor.edge_to_node_vector.size();
@@ -88,9 +88,9 @@ DiamondDualConnectivityBuilder::_buildDiamondConnectivityDescriptor(const Connec
const auto& primal_face_to_cell_matrix = primal_connectivity.faceToCellMatrix();
- 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];
+ parallel_for(primal_number_of_faces, [&](FaceId face_id) {
+ const size_t i_cell = face_id;
+ const auto& primal_face_cell_list = primal_face_to_cell_matrix[face_id];
if (primal_face_cell_list.size() == 1) {
diamond_descriptor.cell_type_vector[i_cell] = CellType::Triangle;
@@ -107,22 +107,22 @@ DiamondDualConnectivityBuilder::_buildDiamondConnectivityDescriptor(const Connec
diamond_descriptor.cell_type_vector[i_cell] = CellType::Diamond;
}
}
- }
+ });
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_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];
+ parallel_for(primal_number_of_faces, [&](FaceId face_id) {
+ const size_t& i_diamond_cell = face_id;
+ const auto& primal_face_cell_list = primal_face_to_cell_matrix[face_id];
+ const auto& primal_face_node_list = primal_face_to_node_matrix[face_id];
if (primal_face_cell_list.size() == 1) {
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);
+ const auto i_face_in_cell = primal_face_local_number_in_their_cells(face_id, 0);
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];
@@ -149,7 +149,7 @@ DiamondDualConnectivityBuilder::_buildDiamondConnectivityDescriptor(const Connec
const CellId cell0_id = primal_face_cell_list[0];
const CellId cell1_id = primal_face_cell_list[1];
- const auto i_face_in_cell0 = primal_face_local_number_in_their_cells(i_face, 0);
+ const auto i_face_in_cell0 = primal_face_local_number_in_their_cells(face_id, 0);
if constexpr (Dimension == 2) {
Assert(primal_face_node_list.size() == 2);
@@ -176,7 +176,7 @@ DiamondDualConnectivityBuilder::_buildDiamondConnectivityDescriptor(const Connec
}
}
}
- }
+ });
}
template <>
@@ -368,7 +368,6 @@ DiamondDualConnectivityBuilder::_buildDiamondConnectivityFrom(const IConnectivit
face_array[i] = diamond_face_list[i];
}
diamond_descriptor.addRefItemList(RefFaceList{primal_ref_face_list.refId(), face_array});
- std::cout << "stored " << primal_ref_face_list.refId() << '\n';
}
}
}
@@ -420,7 +419,6 @@ DiamondDualConnectivityBuilder::_buildDiamondConnectivityFrom(const IConnectivit
edge_array[i] = diamond_edge_list[i];
}
diamond_descriptor.addRefItemList(RefEdgeList{primal_ref_edge_list.refId(), edge_array});
- std::cout << "stored " << primal_ref_edge_list.refId() << '\n';
}
}
}
@@ -508,40 +506,38 @@ DiamondDualConnectivityBuilder::_buildDiamondConnectivityFrom(const IConnectivit
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);
- }();
+ if constexpr (Dimension == 1) {
+ const auto& node_to_cell_matrix = primal_connectivity.nodeToCellMatrix();
- 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++);
+ 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 primal_cell_to_dual_node_map;
- }();
+ }
+ return convert_to_array(primal_node_to_dual_node_vector);
+ }();
- 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;
+ 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;
+ }();
- primal_face_to_dual_cell_map[id] = std::make_pair(primal_face_id, dual_cell_id);
- }
- return primal_face_to_dual_cell_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 =
diff --git a/src/output/VTKWriter.hpp b/src/output/VTKWriter.hpp
index d2459353e4bffa282eac8ff27e21c63d93273664..f3c202ffc51633d291d36150d662fe950023ef58 100644
--- a/src/output/VTKWriter.hpp
+++ b/src/output/VTKWriter.hpp
@@ -321,7 +321,9 @@ class VTKWriter
{
Array<int8_t> types(mesh->numberOfCells());
- const auto& cell_type = mesh->connectivity().cellType();
+ const auto& cell_type = mesh->connectivity().cellType();
+ const auto& cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix();
+
parallel_for(
mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
switch (cell_type[j]) {
@@ -342,7 +344,11 @@ class VTKWriter
break;
}
case CellType::Pyramid: {
- types[j] = 14;
+ if (cell_to_node_matrix[j].size() == 5) {
+ types[j] = 14; // quadrangle basis
+ } else {
+ types[j] = 41; // polygonal basis
+ }
break;
}
case CellType::Prism: {
@@ -365,63 +371,72 @@ 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 constexpr (MeshType::Dimension == 3) {
+ const bool has_general_polyhedron = [&] {
+ for (size_t i = 0; i < types.size(); ++i) {
+ if (types[i] == 41)
+ return true;
}
-
- 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];
+ return false;
+ }();
+ if (has_general_polyhedron) {
+ const auto& cell_to_face_matrix = mesh->connectivity().cellToFaceMatrix();
+ 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;
}
- } else {
- for (size_t i = 0; i < face_node_in_cell.size(); ++i) {
- fout << ' ' << face_node_in_cell[i];
+
+ 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;
}
-
- next_offset += face_nodes.size();
-
- fout << '\n';
+ 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";
}
- 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";