From 92dfa3d9c60f45995bd607bc0824068c342ffd51 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Tue, 28 Jun 2022 07:23:30 +0200 Subject: [PATCH] Use fast access to line indices to speed up get_dof_indices --- include/deal.II/dofs/dof_accessor.templates.h | 238 +++++++++--------- source/dofs/dof_handler_policy.cc | 6 +- 2 files changed, 115 insertions(+), 129 deletions(-) diff --git a/include/deal.II/dofs/dof_accessor.templates.h b/include/deal.II/dofs/dof_accessor.templates.h index 635f6e3e9b..fdbb395ddf 100644 --- a/include/deal.II/dofs/dof_accessor.templates.h +++ b/include/deal.II/dofs/dof_accessor.templates.h @@ -1023,8 +1023,11 @@ namespace internal // DoFs on them if (fe.n_dofs_per_vertex() > 0) for (const auto vertex : accessor.vertex_indices()) - dof_operation.process_vertex_dofs( - accessor, vertex, dof_indices_ptr, fe_index, dof_processor); + dof_operation.process_vertex_dofs(*accessor.dof_handler, + accessor.vertex_index(vertex), + fe_index, + dof_indices_ptr, + dof_processor); // 2) copy dof numbers from the LINE. for lines with the wrong // orientation (which might occur in 3d), we have already made sure that @@ -1034,27 +1037,41 @@ namespace internal // adjust the shape function indices that we see to correspond to the // correct (face/cell-local) ordering. if ((structdim == 2 || structdim == 3) && fe.n_dofs_per_line() > 0) - for (const auto line : accessor.line_indices()) - { - const bool line_orientation = accessor.line_orientation(line); - if (line_orientation) - dof_operation.process_dofs( - *accessor.line(line), - [](const auto d) { return d; }, - dof_indices_ptr, - fe_index, - dof_processor); - else - dof_operation.process_dofs( - *accessor.line(line), - [&fe, &line_orientation](const auto d) { - return fe.adjust_line_dof_index_for_line_orientation( - d, line_orientation); - }, - dof_indices_ptr, - fe_index, - dof_processor); - } + { + const auto line_indices = internal::TriaAccessorImplementation:: + Implementation::get_line_indices_of_cell(accessor); + const auto line_orientations = + internal::TriaAccessorImplementation::Implementation:: + get_line_orientations_of_cell(accessor); + + for (const auto line : accessor.line_indices()) + { + const bool line_orientation = line_orientations[line]; + if (line_orientation) + dof_operation.process_dofs( + accessor.get_dof_handler(), + 0, + line_indices[line], + fe_index, + [](const auto d) { return d; }, + std::integral_constant(), + dof_indices_ptr, + dof_processor); + else + dof_operation.process_dofs( + accessor.get_dof_handler(), + 0, + line_indices[line], + fe_index, + [&fe, line_orientation](const auto d) { + return fe.adjust_line_dof_index_for_line_orientation( + d, line_orientation); + }, + std::integral_constant(), + dof_indices_ptr, + dof_processor); + } + } // 3) copy dof numbers from the FACE. for faces with the wrong // orientation, we have already made sure that we're ok by picking the @@ -1069,16 +1086,23 @@ namespace internal { const unsigned int raw_orientation = TriaAccessorImplementation:: Implementation::face_orientation_raw(accessor, quad); + const unsigned int quad_index = accessor.quad_index(quad); if (raw_orientation == 1) dof_operation.process_dofs( - *accessor.quad(quad), + accessor.get_dof_handler(), + 0, + quad_index, + fe_index, [](const auto d) { return d; }, + std::integral_constant(), dof_indices_ptr, - fe_index, dof_processor); else dof_operation.process_dofs( - *accessor.quad(quad), + accessor.get_dof_handler(), + 0, + quad_index, + fe_index, [&](const auto d) { return fe.adjust_quad_dof_index_for_face_orientation( d, @@ -1087,8 +1111,8 @@ namespace internal accessor.face_flip(quad), accessor.face_rotation(quad)); }, + std::integral_constant(), dof_indices_ptr, - fe_index, dof_processor); } @@ -1100,10 +1124,13 @@ namespace internal fe.max_dofs_per_quad() : fe.template n_dofs_per_object()) > 0) dof_operation.process_dofs( - accessor, + accessor.get_dof_handler(), + accessor.level(), + accessor.index(), + fe_index, [&](const auto d) { return d; }, + std::integral_constant(), dof_indices_ptr, - fe_index, dof_processor); AssertDimension(n_dof_indices(accessor, fe_index, count_level_dofs), @@ -1126,30 +1153,24 @@ namespace internal * An internal struct encapsulating the task of getting (vertex) * DoF indices. */ - template + template struct DoFIndexProcessor { /** * Return vertex DoF indices. */ - template + template DEAL_II_ALWAYS_INLINE void - process_vertex_dofs( - const dealii::DoFAccessor - & accessor, - const unsigned int vertex, - types::global_dof_index *&dof_indices_ptr, - const unsigned int fe_index_, - const DoFProcessor & dof_processor) const + process_vertex_dofs(DoFHandler &dof_handler, + const unsigned int vertex_index, + const unsigned int fe_index, + types::global_dof_index *& dof_indices_ptr, + const DoFProcessor & dof_processor) const { - const auto fe_index = - internal::DoFAccessorImplementation::get_fe_index_or_default( - accessor, fe_index_); - process_object( - accessor.get_dof_handler(), + dof_handler, 0, - accessor.vertex_index(vertex), + vertex_index, fe_index, [](const auto d) { Assert(false, ExcInternalError()); @@ -1165,41 +1186,24 @@ namespace internal */ template DEAL_II_ALWAYS_INLINE void - process_dofs( - const dealii::DoFAccessor - & accessor, - const DoFMapping & mapping, - types::global_dof_index *&dof_indices_ptr, - const unsigned int fe_index_, - const DoFProcessor & dof_processor) const - { - const auto fe_index = - internal::DoFAccessorImplementation::get_fe_index_or_default( - accessor, fe_index_); - - process_object(accessor.get_dof_handler(), - accessor.level(), - accessor.index(), - fe_index, - mapping, - std::integral_constant(), - dof_indices_ptr, - dof_processor); - } - - /** - * Fallback for DoFInvalidAccessor. - */ - template - DEAL_II_ALWAYS_INLINE void - process_dofs( - const dealii::DoFInvalidAccessor &, - const DoFMapping &, - types::global_dof_index *&, - const unsigned int, - const DoFProcessor &) const + process_dofs(const DoFHandler &dof_handler, + const unsigned int obj_level, + const unsigned int obj_index, + const unsigned int fe_index, + const DoFMapping & mapping, + const std::integral_constant, + types::global_dof_index *&dof_indices_ptr, + const DoFProcessor & dof_processor) const { - Assert(false, ExcInternalError()); + process_object( + dof_handler, + obj_level, + obj_index, + fe_index, + mapping, + std::integral_constant(), + dof_indices_ptr, + dof_processor); } }; @@ -1222,21 +1226,20 @@ namespace internal /** * Return vertex DoF indices. */ - template + template DEAL_II_ALWAYS_INLINE void - process_vertex_dofs( - const dealii::DoFAccessor - & accessor, - const unsigned int vertex, - types::global_dof_index *&dof_indices_ptr, - const unsigned int fe_index, - const DoFProcessor & dof_processor) const + process_vertex_dofs(DoFHandler &dof_handler, + const unsigned int vertex_index, + const unsigned int, + types::global_dof_index *&dof_indices_ptr, + const DoFProcessor & dof_processor) const { const unsigned int n_indices = - accessor.get_fe(fe_index).template n_dofs_per_object<0>(); + dof_handler.get_fe(0).template n_dofs_per_object<0>(); types::global_dof_index *stored_indices = - &accessor.dof_handler->mg_vertex_dofs[accessor.vertex_index(vertex)] - .access_index(level, 0, n_indices); + &dof_handler.mg_vertex_dofs[vertex_index].access_index(level, + 0, + n_indices); for (unsigned int d = 0; d < n_indices; ++d, ++dof_indices_ptr) dof_processor(stored_indices[d], dof_indices_ptr); } @@ -1244,49 +1247,32 @@ namespace internal /** * Return DoF indices for lines, quads, and inner degrees of freedom. */ - template + template DEAL_II_ALWAYS_INLINE void - process_dofs( - const dealii::DoFAccessor - & accessor, - const DoFMapping & mapping, - types::global_dof_index *&dof_indices_ptr, - const unsigned int fe_index, - const DoFProcessor & dof_processor) const + process_dofs(const DoFHandler &dof_handler, + const unsigned int, + const unsigned int obj_index, + const unsigned int fe_index, + const DoFMapping & mapping, + const std::integral_constant, + types::global_dof_index *&dof_indices_ptr, + const DoFProcessor & dof_processor) const { const unsigned int n_indices = - accessor.get_fe(fe_index).template n_dofs_per_object(); - types::global_dof_index *stored_indices = - &get_mg_dof_index(*accessor.dof_handler, - accessor.dof_handler->mg_levels[level], - accessor.dof_handler->mg_faces, - accessor.index(), - fe_index, - 0, - std::integral_constant()); + dof_handler.get_fe(0).template n_dofs_per_object(); + types::global_dof_index *stored_indices = &get_mg_dof_index( + dof_handler, + dof_handler.mg_levels[level], + dof_handler.mg_faces, + obj_index, + fe_index, + 0, + std::integral_constant()); for (unsigned int d = 0; d < n_indices; ++d, ++dof_indices_ptr) dof_processor(stored_indices[structdim < dim ? mapping(d) : d], dof_indices_ptr); } - /** - * Fallback for DoFInvalidAccessor. - */ - template - DEAL_II_ALWAYS_INLINE void - process_dofs( - const dealii::DoFInvalidAccessor &, - const DoFMapping &, - types::global_dof_index *&, - const unsigned int, - const DoFProcessor &) const - { - Assert(false, ExcInternalError()); - } - private: const unsigned int level; }; @@ -1305,7 +1291,7 @@ namespace internal accessor, dof_indices, fe_index, - DoFIndexProcessor(), + DoFIndexProcessor(), [](auto stored_index, auto dof_ptr) { *dof_ptr = stored_index; }, false); } @@ -1335,7 +1321,7 @@ namespace internal accessor, dof_indices, fe_index, - DoFIndexProcessor(), + DoFIndexProcessor(), [](auto &stored_index, auto dof_ptr) { stored_index = *dof_ptr; }, false); } diff --git a/source/dofs/dof_handler_policy.cc b/source/dofs/dof_handler_policy.cc index bc2f88e7fd..90f57d392a 100644 --- a/source/dofs/dof_handler_policy.cc +++ b/source/dofs/dof_handler_policy.cc @@ -1725,7 +1725,7 @@ namespace internal std::make_tuple(), cell->active_fe_index(), DoFAccessorImplementation::Implementation:: - DoFIndexProcessor(), + DoFIndexProcessor(), [&next_free_dof](auto &stored_index, auto) { if (stored_index == numbers::invalid_dof_index) { @@ -3575,7 +3575,7 @@ namespace internal dofs, cell->active_fe_index(), DoFAccessorImplementation::Implementation:: - DoFIndexProcessor(), + DoFIndexProcessor(), [&complete](auto &stored_index, auto received_index) { if (*received_index != numbers::invalid_dof_index) { @@ -4105,7 +4105,7 @@ namespace internal std::make_tuple(), cell->active_fe_index(), DoFAccessorImplementation::Implementation:: - DoFIndexProcessor(), + DoFIndexProcessor(), [&owned_dofs](auto &stored_index, auto) { // delete a DoF index if it has not already been // deleted (e.g., by visiting a neighboring cell, if -- 2.39.5