From 19c84e26271d0ec8295c7bd87c1cd3861cf6085f Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Tue, 28 Jun 2022 07:23:02 +0200 Subject: [PATCH] Provide fast access function to all line indices of cell --- .../deal.II/grid/tria_accessor.templates.h | 196 ++++++++++++++++++ 1 file changed, 196 insertions(+) diff --git a/include/deal.II/grid/tria_accessor.templates.h b/include/deal.II/grid/tria_accessor.templates.h index 3ff63d7c74..63c46968bf 100644 --- a/include/deal.II/grid/tria_accessor.templates.h +++ b/include/deal.II/grid/tria_accessor.templates.h @@ -1050,6 +1050,202 @@ namespace internal return accessor.quad(face_index)->vertex_index(vertex_index); } + + + + template + static std::array + get_line_indices_of_cell(const TriaAccessor<1, dim, spacedim> &) + { + Assert(false, ExcInternalError()); + return {}; + } + + + + template + static std::array + get_line_indices_of_cell(const TriaAccessor<2, dim, spacedim> &cell) + { + // For 2D cells the access cell->line_orientation() is already + // efficient + std::array line_indices; + for (unsigned int line : cell.line_indices()) + line_indices[line] = cell.line_index(line); + return line_indices; + } + + /** + * A helper function to provide faster access to cell->line_index() in + * 3D + */ + template + static std::array + get_line_indices_of_cell( + const TriaAccessor &cell) + { + std::array line_indices; + + // For hexahedra, the classical access via quads -> lines is too + // inefficient. Unroll this code here to allow the compiler to inline + // the necessary functions. + const auto ref_cell = cell.reference_cell(); + if (ref_cell == ReferenceCells::Hexahedron) + { + for (unsigned int f = 4; f < 6; ++f) + { + const unsigned char orientation = + cell.get_triangulation() + .levels[cell.level()] + ->face_orientations[cell.index() * 6 + f]; + + // It might seem superfluous to spell out the four indices + // that get later consumed by a for loop over these four + // elements; however, for the compiler it is easier to inline + // the statement of standard_to_real_face_line() when next to + // each other, as opposed to be interleaved with a + // line_index() call. + const std::array my_indices{ + {ref_cell.standard_to_real_face_line(0, f, orientation), + ref_cell.standard_to_real_face_line(1, f, orientation), + ref_cell.standard_to_real_face_line(2, f, orientation), + ref_cell.standard_to_real_face_line(3, f, orientation)}}; + const auto quad = cell.quad(f); + for (unsigned int l = 0; l < 4; ++l) + line_indices[4 * (f - 4) + l] = + quad->line_index(my_indices[l]); + } + for (unsigned int f = 0; f < 2; ++f) + { + const unsigned char orientation = + cell.get_triangulation() + .levels[cell.level()] + ->face_orientations[cell.index() * 6 + f]; + const std::array my_indices{ + {ref_cell.standard_to_real_face_line(0, f, orientation), + ref_cell.standard_to_real_face_line(1, f, orientation)}}; + const auto quad = cell.quad(f); + line_indices[8 + f] = quad->line_index(my_indices[0]); + line_indices[10 + f] = quad->line_index(my_indices[1]); + } + } + else + // For other shapes (tetrahedra, wedges, pyramids), we do not + // currently implement an optimized function. + for (unsigned int l = 0; l < std::min(12U, cell.n_lines()); ++l) + line_indices[l] = cell.line_index(l); + + return line_indices; + } + + + + /** + * A helper function to provide faster access to + * cell->line_orientation(), 1D specialization + */ + template + static std::array + get_line_orientations_of_cell(const TriaAccessor<1, dim, spacedim> &) + { + Assert(false, ExcInternalError()); + return {}; + } + + + + /** + * A helper function to provide faster access to + * cell->line_orientation(), 2D specialization + */ + template + static std::array + get_line_orientations_of_cell(const TriaAccessor<2, dim, spacedim> &cell) + { + // For 2D cells the access cell->line_orientation() is already + // efficient + std::array line_orientations; + for (unsigned int line : cell.line_indices()) + line_orientations[line] = cell.line_orientation(line); + return line_orientations; + } + + + + /** + * A helper function to provide faster access to + * cell->line_orientation(), 3D specialization + */ + template + static std::array + get_line_orientations_of_cell(const TriaAccessor<3, dim, spacedim> &cell) + { + std::array line_orientations; + + // For hexahedra, the classical access via quads -> lines is too + // inefficient. Unroll this code here to allow the compiler to inline + // the necessary functions. + const auto ref_cell = cell.reference_cell(); + if (ref_cell == ReferenceCells::Hexahedron) + { + for (unsigned int f = 4; f < 6; ++f) + { + const unsigned char orientation = + cell.get_triangulation() + .levels[cell.level()] + ->face_orientations[cell.index() * 6 + f]; + + // It might seem superfluous to spell out the four indices and + // orientations that get later consumed by a for loop over + // these four elements; however, for the compiler it is easier + // to inline the statement of standard_to_real_face_line() + // when next to each other, as opposed to be interleaved with + // a line_index() call. + const std::array my_indices{ + {ref_cell.standard_to_real_face_line(0, f, orientation), + ref_cell.standard_to_real_face_line(1, f, orientation), + ref_cell.standard_to_real_face_line(2, f, orientation), + ref_cell.standard_to_real_face_line(3, f, orientation)}}; + const auto quad = cell.quad(f); + const std::array my_orientations{ + {ref_cell.standard_vs_true_line_orientation( + 0, orientation, quad->line_orientation(my_indices[0])), + ref_cell.standard_vs_true_line_orientation( + 1, orientation, quad->line_orientation(my_indices[1])), + ref_cell.standard_vs_true_line_orientation( + 2, orientation, quad->line_orientation(my_indices[2])), + ref_cell.standard_vs_true_line_orientation( + 3, orientation, quad->line_orientation(my_indices[3]))}}; + for (unsigned int l = 0; l < 4; ++l) + line_orientations[4 * (f - 4) + l] = my_orientations[l]; + } + for (unsigned int f = 0; f < 2; ++f) + { + const unsigned char orientation = + cell.get_triangulation() + .levels[cell.level()] + ->face_orientations[cell.index() * 6 + f]; + const std::array my_indices{ + {ref_cell.standard_to_real_face_line(0, f, orientation), + ref_cell.standard_to_real_face_line(1, f, orientation)}}; + const auto quad = cell.quad(f); + const std::array my_orientations{ + {ref_cell.standard_vs_true_line_orientation( + 0, orientation, quad->line_orientation(my_indices[0])), + ref_cell.standard_vs_true_line_orientation( + 1, orientation, quad->line_orientation(my_indices[1]))}}; + line_orientations[8 + f] = my_orientations[0]; + line_orientations[10 + f] = my_orientations[1]; + } + } + else + // For other shapes (tetrahedra, wedges, pyramids), we do not + // currently implement an optimized function + for (unsigned int l = 0; l < std::min(12U, cell.n_lines()); ++l) + line_orientations[l] = cell.line_orientation(l); + + return line_orientations; + } }; } // namespace TriaAccessorImplementation } // namespace internal -- 2.39.5