From 54b4ebaf2ae38019366b24b27abf9ba9fff3506c Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Wed, 8 Jun 2022 12:03:37 +0200 Subject: [PATCH] Speed up computation of cell vertex indices cache --- source/grid/tria.cc | 37 ++++++++++++++++++++++++++++++++++--- 1 file changed, 34 insertions(+), 3 deletions(-) diff --git a/source/grid/tria.cc b/source/grid/tria.cc index 472f3979c5..f2d99de001 100644 --- a/source/grid/tria.cc +++ b/source/grid/tria.cc @@ -14561,9 +14561,40 @@ Triangulation::reset_cell_vertex_indices_cache() for (const auto &cell : cell_iterators_on_level(l)) { const unsigned int my_index = cell->index() * max_vertices_per_cell; - for (const unsigned int i : cell->vertex_indices()) - cache[my_index + i] = internal::TriaAccessorImplementation:: - Implementation::vertex_index(*cell, i); + + // to reduce the cost of this function when passing down into quads, + // then lines, then vertices, we use a more low-level access method + // for hexahedral cells, where we can streamline most of the logic + const ReferenceCell ref_cell = cell->reference_cell(); + if (ref_cell == ReferenceCells::Hexahedron) + for (unsigned int face = 4; face < 6; ++face) + { + const auto face_iter = cell->face(face); + const std::array line_orientations{ + {face_iter->line_orientation(0), + face_iter->line_orientation(1)}}; + std::array raw_vertex_indices{ + {face_iter->line(0)->vertex_index(1 - line_orientations[0]), + face_iter->line(1)->vertex_index(1 - line_orientations[1]), + face_iter->line(0)->vertex_index(line_orientations[0]), + face_iter->line(1)->vertex_index(line_orientations[1])}}; + + const unsigned char orientate = + levels[l]->face_orientations[cell->index() * 6 + face]; + std::array vertex_order{ + {ref_cell.standard_to_real_face_vertex(0, face, orientate), + ref_cell.standard_to_real_face_vertex(1, face, orientate), + ref_cell.standard_to_real_face_vertex(2, face, orientate), + ref_cell.standard_to_real_face_vertex(3, face, orientate)}}; + + const unsigned int index = my_index + 4 * (face - 4); + for (unsigned int i = 0; i < 4; ++i) + cache[index + i] = raw_vertex_indices[vertex_order[i]]; + } + else + for (const unsigned int i : cell->vertex_indices()) + cache[my_index + i] = internal::TriaAccessorImplementation:: + Implementation::vertex_index(*cell, i); } } } -- 2.39.5