Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Speed up DoFAccessorImplementation::get_dof_indices #14066

Merged
merged 3 commits into from
Jun 29, 2022
Merged

Speed up DoFAccessorImplementation::get_dof_indices #14066

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
19c84e2
Provide fast access function to all line indices of cell
kronbichler Jun 28, 2022
92dfa3d
Use fast access to line indices to speed up get_dof_indices
kronbichler Jun 28, 2022
4d99dd6
Add assertions for fe_index in MG case
kronbichler Jun 29, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
249 changes: 122 additions & 127 deletions include/deal.II/dofs/dof_accessor.templates.h
View file
Open in desktop
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>();
Comment on lines -1236 to +1242
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What happened with fe_index?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

OK. This is the MG path. Wouldn't it make sense to add an assert that fe_index is ineed 0 or default?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, that makes sense.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Note that I added the check in line 1344 below, not here, as I do not want to pick up the index here again.

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