Skip to content
Snippets Groups Projects

Feature/diamond dual mesh manager

1 file
+ 67
52
Compare changes
  • Side-by-side
  • Inline
+ 67
52
@@ -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";
Loading