From ff1aac08840a4201174ef93f53cf66627fe0ca51 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Tue, 10 Jan 2023 10:30:16 -0700 Subject: [PATCH] Make special casing in VTU output more explicit. This also allows adding an assertion for the non-special case. --- source/base/data_out_base.cc | 51 +++++++++++++++++++++++++++++++++--- 1 file changed, 48 insertions(+), 3 deletions(-) diff --git a/source/base/data_out_base.cc b/source/base/data_out_base.cc index fb6ccc6adf..f835cc354c 100644 --- a/source/base/data_out_base.cc +++ b/source/base/data_out_base.cc @@ -730,6 +730,15 @@ namespace /** * Return the tuple (vtk cell type, number of cells, number of nodes) * for a patch. + * + * The logic used here is as follows: + * - If a cell is not subdivided or we don't use higher order cells, + * then we use linear cells + * - For hypercubes, we support subdividing cells into sub-cells, + * which are then treated as each being linear + * - For triangles and tetrahedra, we special-case the situation of + * n_subdivisions==2, in which case we treat the cell as a single + * quadratic cell (i.e., higher order) */ template std::array @@ -753,12 +762,14 @@ namespace else if (patch.reference_cell == ReferenceCells::Triangle && patch.data.n_cols() == 6) { + Assert(patch.n_subdivisions == 2, ExcInternalError()); vtk_cell_id[0] = patch.reference_cell.vtk_quadratic_type(); vtk_cell_id[2] = patch.data.n_cols(); } else if (patch.reference_cell == ReferenceCells::Tetrahedron && patch.data.n_cols() == 10) { + Assert(patch.n_subdivisions == 2, ExcInternalError()); vtk_cell_id[0] = patch.reference_cell.vtk_quadratic_type(); vtk_cell_id[2] = patch.data.n_cols(); } @@ -5867,10 +5878,44 @@ namespace DataOutBase for (const auto &patch : patches) { - // First treat non-hypercubes since they can currently - // not be subdivided (into sub-cells, or into higher-order cells): - if (patch.reference_cell != ReferenceCells::get_hypercube()) + // First treat a slight oddball case: For triangles and tetrahedra, + // the case with n_subdivisions==2 is treated as if the cell was + // output as a single, quadratic, cell rather than as one would + // expect as 4 sub-cells (for triangles; and the corresponding + // number of sub-cells for tetrahedra). This is courtesy of some + // special-casing in the function extract_vtk_patch_info(). + if ((dim >= 2) && + (patch.reference_cell == ReferenceCells::get_simplex()) && + (patch.n_subdivisions == 2)) { + const unsigned int n_points = patch.data.n_cols(); + Assert((dim == 2 && n_points == 6) || + (dim == 3 && n_points == 10), + ExcInternalError()); + + if (deal_ii_with_zlib && + (flags.compression_level != + DataOutBase::CompressionLevel::plain_text)) + { + for (unsigned int i = 0; i < n_points; ++i) + cells.push_back(first_vertex_of_patch + i); + } + else + { + for (unsigned int i = 0; i < n_points; ++i) + o << '\t' << first_vertex_of_patch + i; + o << '\n'; + } + + first_vertex_of_patch += n_points; + } + // Then treat all of the other non-hypercube cases since they can + // currently not be subdivided (into sub-cells, or into higher-order + // cells): + else if (patch.reference_cell != ReferenceCells::get_hypercube()) + { + Assert(patch.n_subdivisions == 1, ExcNotImplemented()); + const unsigned int n_points = patch.data.n_cols(); static const std::array pyramid_index_translation_table = {{0, 1, 3, 2, 4}}; -- 2.39.5