Skip to content

Commit

Permalink
Merge pull request dealii#14066 from kronbichler/speed_up_access_dof_…
Browse files Browse the repository at this point in the history 
…indices

Speed up DoFAccessorImplementation::get_dof_indices
  • Loading branch information
@drwells
drwells committed Jun 29, 2022
2 parents c1f8b19 + 9acc046 commit 0a13f11
Show file tree
Hide file tree
Showing 3 changed files with 321 additions and 130 deletions.
249 changes: 122 additions & 127 deletions include/deal.II/dofs/dof_accessor.templates.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -997,11 +997,15 @@ namespace internal
const dealii::DoFAccessor<structdim, dim, spacedim, level_dof_access>
& accessor,
const DoFIndicesType &const_dof_indices,
const unsigned int fe_index,
const unsigned int fe_index_,
const DoFOperation & dof_operation,
const DoFProcessor & dof_processor,
const bool count_level_dofs)
{
const unsigned int fe_index =
internal::DoFAccessorImplementation::get_fe_index_or_default(
accessor, fe_index_);

// we cannot rely on the template parameter level_dof_access here, since
// the function get_mg_dof_indices()/set_mg_dof_indices() can be called
// even if level_dof_access==false.
Expand All @@ -1023,8 +1027,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
Expand All @@ -1034,27 +1041,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<int, 1>(),
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<int, 1>(),
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
Expand All @@ -1069,16 +1090,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<int, 2>(),
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,
Expand All @@ -1087,8 +1115,8 @@ namespace internal
accessor.face_flip(quad),
accessor.face_rotation(quad));
},
std::integral_constant<int, 2>(),
dof_indices_ptr,
fe_index,
dof_processor);
}

Expand All @@ -1100,10 +1128,13 @@ namespace internal
fe.max_dofs_per_quad() :
fe.template n_dofs_per_object<structdim>()) > 0)
dof_operation.process_dofs(
accessor,
accessor.get_dof_handler(),
accessor.level(),
accessor.index(),
fe_index,
[&](const auto d) { return d; },
std::integral_constant<int, structdim>(),
dof_indices_ptr,
fe_index,
dof_processor);

AssertDimension(n_dof_indices(accessor, fe_index, count_level_dofs),
Expand All @@ -1126,30 +1157,24 @@ namespace internal
* An internal struct encapsulating the task of getting (vertex)
* DoF indices.
*/
template <int dim, int spacedim, bool level_dof_access>
template <int dim, int spacedim>
struct DoFIndexProcessor
{
/**
* Return vertex DoF indices.
*/
template <int structdim, typename DoFProcessor>
template <typename DoFProcessor>
DEAL_II_ALWAYS_INLINE void
process_vertex_dofs(
const dealii::DoFAccessor<structdim, dim, spacedim, level_dof_access>
& 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<dim, spacedim> &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());
Expand All @@ -1165,41 +1190,24 @@ namespace internal
*/
template <int structdim, typename DoFMapping, typename DoFProcessor>
DEAL_II_ALWAYS_INLINE void
process_dofs(
const dealii::DoFAccessor<structdim, dim, spacedim, level_dof_access>
& 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<int, structdim>(),
dof_indices_ptr,
dof_processor);
}

/**
* Fallback for DoFInvalidAccessor.
*/
template <int structdim, typename DoFMapping, typename DoFProcessor>
DEAL_II_ALWAYS_INLINE void
process_dofs(
const dealii::DoFInvalidAccessor<structdim, dim, spacedim> &,
const DoFMapping &,
types::global_dof_index *&,
const unsigned int,
const DoFProcessor &) const
process_dofs(const DoFHandler<dim, spacedim> &dof_handler,
const unsigned int obj_level,
const unsigned int obj_index,
const unsigned int fe_index,
const DoFMapping & mapping,
const std::integral_constant<int, structdim>,
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<int, std::min(structdim, dim)>(),
dof_indices_ptr,
dof_processor);
}
};

Expand All @@ -1222,71 +1230,53 @@ namespace internal
/**
* Return vertex DoF indices.
*/
template <int structdim, typename DoFProcessor, bool level_dof_access>
template <typename DoFProcessor>
DEAL_II_ALWAYS_INLINE void
process_vertex_dofs(
const dealii::DoFAccessor<structdim, dim, spacedim, level_dof_access>
& 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<dim, spacedim> &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);
}

/**
* Return DoF indices for lines, quads, and inner degrees of freedom.
*/
template <int structdim,
typename DoFMapping,
typename DoFProcessor,
bool level_dof_access>
template <int structdim, typename DoFMapping, typename DoFProcessor>
DEAL_II_ALWAYS_INLINE void
process_dofs(
const dealii::DoFAccessor<structdim, dim, spacedim, level_dof_access>
& 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<dim, spacedim> &dof_handler,
const unsigned int,
const unsigned int obj_index,
const unsigned int fe_index,
const DoFMapping & mapping,
const std::integral_constant<int, structdim>,
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<structdim>();
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<int, structdim>());
dof_handler.get_fe(0).template n_dofs_per_object<structdim>();
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<int, std::min(structdim, dim)>());
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 <int structdim, typename DoFMapping, typename DoFProcessor>
DEAL_II_ALWAYS_INLINE void
process_dofs(
const dealii::DoFInvalidAccessor<structdim, dim, spacedim> &,
const DoFMapping &,
types::global_dof_index *&,
const unsigned int,
const DoFProcessor &) const
{
Assert(false, ExcInternalError());
}

private:
const unsigned int level;
};
Expand All @@ -1305,7 +1295,7 @@ namespace internal
accessor,
dof_indices,
fe_index,
DoFIndexProcessor<dim, spacedim, level_dof_access>(),
DoFIndexProcessor<dim, spacedim>(),
[](auto stored_index, auto dof_ptr) { *dof_ptr = stored_index; },
false);
}
Expand Down Expand Up @@ -1335,7 +1325,7 @@ namespace internal
accessor,
dof_indices,
fe_index,
DoFIndexProcessor<dim, spacedim, level_dof_access>(),
DoFIndexProcessor<dim, spacedim>(),
[](auto &stored_index, auto dof_ptr) { stored_index = *dof_ptr; },
false);
}
Expand All @@ -1351,6 +1341,8 @@ namespace internal
std::vector<types::global_dof_index> &dof_indices,
const unsigned int fe_index)
{
Assert((fe_index == DoFHandler<dim, spacedim>::default_fe_index),
ExcMessage("MG DoF indices cannot be queried in hp case"));
process_dof_indices(
accessor,
dof_indices,
Expand All @@ -1371,6 +1363,9 @@ namespace internal
const std::vector<types::global_dof_index> &dof_indices,
const unsigned int fe_index)
{
Assert((fe_index == DoFHandler<dim, spacedim>::default_fe_index),
ExcMessage("MG DoF indices cannot be queried in hp case"));

// Note: this function is as general as `get_mg_dof_indices()`. This
// assert is placed here since it is currently only used by the
// function DoFCellAccessor::set_mg_dof_indices(), which is called by
Expand Down

0 comments on commit 0a13f11

Please sign in to comment.