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

Use std::variant instead of a hand-rolled version. #16355

Merged
merged 1 commit into from
Dec 17, 2023
Merged

Use std::variant instead of a hand-rolled version. #16355

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
40 changes: 21 additions & 19 deletions include/deal.II/fe/fe_values_base.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -48,6 +48,7 @@
#include <memory>
#include <optional>
#include <type_traits>
#include <variant>

DEAL_II_NAMESPACE_OPEN

Expand Down Expand Up @@ -1571,23 +1572,29 @@ class FEValuesBase : public Subscriptor
"FEValues::reinit() with an iterator type that allows to extract "
"degrees of freedom, such as DoFHandler<dim,spacedim>::cell_iterator.");

/**
* Constructor. Creates an unusable object that is not associated with
* any cell at all.
*/
CellIteratorContainer() = default;

/**
* Constructor.
*/
CellIteratorContainer();
CellIteratorContainer(
const typename Triangulation<dim, spacedim>::cell_iterator &cell);

/**
* Constructor.
*/
template <bool lda>
CellIteratorContainer(
const TriaIterator<DoFCellAccessor<dim, spacedim, lda>> &cell);
const typename DoFHandler<dim, spacedim>::cell_iterator &cell);

/**
* Constructor.
*/
CellIteratorContainer(
const typename Triangulation<dim, spacedim>::cell_iterator &cell);
const typename DoFHandler<dim, spacedim>::level_cell_iterator &cell);

/**
* Indicate whether FEValues::reinit() was called.
Expand Down Expand Up @@ -1628,9 +1635,16 @@ class FEValuesBase : public Subscriptor
Vector<IndexSet::value_type> &out) const;

private:
std::optional<typename Triangulation<dim, spacedim>::cell_iterator> cell;
const DoFHandler<dim, spacedim> *dof_handler;
bool level_dof_access;
/**
* The cell in question, if one has been assigned to this object. The
* concrete data type can either be a Triangulation cell iterator, a
* DoFHandler cell iterator, or a DoFHandler level cell iterator.
*/
std::optional<
std::variant<typename Triangulation<dim, spacedim>::cell_iterator,
typename DoFHandler<dim, spacedim>::cell_iterator,
typename DoFHandler<dim, spacedim>::level_cell_iterator>>
cell;
Comment on lines -1631 to +1647
Copy link
Member Author

Choose a reason for hiding this comment

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

This is the key part of the patch. We really do store a variant already, we just break the data stored into the base object plus potentially the members of a derived class. Using the std::variant is, in my view, a clear win.

};

/**
Expand Down Expand Up @@ -1784,18 +1798,6 @@ class FEValuesBase : public Subscriptor

/*---------------------- Inline functions: FEValuesBase ---------------------*/

template <int dim, int spacedim>
template <bool lda>
inline FEValuesBase<dim, spacedim>::CellIteratorContainer::
CellIteratorContainer(
const TriaIterator<DoFCellAccessor<dim, spacedim, lda>> &cell)
: cell(cell)
, dof_handler(&cell->get_dof_handler())
, level_dof_access(lda)
{}



template <int dim, int spacedim>
inline const FEValuesViews::Scalar<dim, spacedim> &
FEValuesBase<dim, spacedim>::operator[](
Expand Down
105 changes: 69 additions & 36 deletions source/fe/fe_values_base.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -107,21 +107,27 @@ namespace internal

/* ------------ FEValuesBase<dim,spacedim>::CellIteratorContainer ----------- */


template <int dim, int spacedim>
FEValuesBase<dim, spacedim>::CellIteratorContainer::CellIteratorContainer()
: cell()
, dof_handler(nullptr)
, level_dof_access(false)
FEValuesBase<dim, spacedim>::CellIteratorContainer::CellIteratorContainer(
const typename Triangulation<dim, spacedim>::cell_iterator &cell)
: cell(cell)
{}



template <int dim, int spacedim>
FEValuesBase<dim, spacedim>::CellIteratorContainer::CellIteratorContainer(
const typename Triangulation<dim, spacedim>::cell_iterator &cell)
const typename DoFHandler<dim, spacedim>::cell_iterator &cell)
: cell(cell)
{}



template <int dim, int spacedim>
FEValuesBase<dim, spacedim>::CellIteratorContainer::CellIteratorContainer(
const typename DoFHandler<dim, spacedim>::level_cell_iterator &cell)
: cell(cell)
, dof_handler(nullptr)
, level_dof_access(false)
{}


Expand All @@ -141,7 +147,14 @@ operator typename Triangulation<dim, spacedim>::cell_iterator() const
{
Assert(is_initialized(), ExcNotReinited());

return cell.value();
// We can always convert to a tria iterator, regardless of which of
// the three types of cell we store.
return std::visit(
[](auto &cell_iterator) ->
typename Triangulation<dim, spacedim>::cell_iterator {
return cell_iterator;
},
cell.value());
}


Expand All @@ -152,9 +165,17 @@ FEValuesBase<dim, spacedim>::CellIteratorContainer::n_dofs_for_dof_handler()
const
{
Assert(is_initialized(), ExcNotReinited());
Assert(dof_handler != nullptr, ExcNeedsDoFHandler());

return dof_handler->n_dofs();
switch (cell.value().index())
{
case 1:
return std::get<1>(cell.value())->get_dof_handler().n_dofs();
case 2:
return std::get<2>(cell.value())->get_dof_handler().n_dofs();
Copy link
Member

Choose a reason for hiding this comment

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

Unrelated to this patch, but shouldn't we return

Suggested change
return std::get<2>(cell.value())->get_dof_handler().n_dofs();
return std::get<2>(cell.value())->get_dof_handler().n_dofs(std::get<2>(cell.value())->level());

i.e., the number of level dofs rather than the number of dofs globally?

Copy link
Member Author

Choose a reason for hiding this comment

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

That's not what the code originally did, nor what is apparently necessary to pass the test suite. This function is used in contexts such as this:

template <int dim, int spacedim>
template <typename Number>
void
FEValuesBase<dim, spacedim>::get_function_values(
  const ReadVector<Number> &fe_function,
  std::vector<Number>      &values) const
{
  Assert(this->update_flags & update_values,
         ExcAccessToUninitializedField("update_values"));
  AssertDimension(fe->n_components(), 1);
  Assert(present_cell.is_initialized(), ExcNotReinited());
  AssertDimension(fe_function.size(), present_cell.n_dofs_for_dof_handler());

  // get function values of dofs on this cell
  Vector<Number> dof_values(dofs_per_cell);
  present_cell.get_interpolated_dof_values(fe_function, dof_values);

We pass in a global vector and interpolate it onto the current cell (active or not).

Copy link
Member Author

Choose a reason for hiding this comment

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

In other words, I believe that the change I made is correct.

Copy link
Member

Choose a reason for hiding this comment

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

Thanks for checking, this is not what I expected to happen, but I do see it now.

default:
Assert(false, ExcNeedsDoFHandler());
return numbers::invalid_dof_index;
}
}


Expand All @@ -167,20 +188,21 @@ FEValuesBase<dim, spacedim>::CellIteratorContainer::get_interpolated_dof_values(
Vector<Number> &out) const
{
Assert(is_initialized(), ExcNotReinited());
Assert(dof_handler != nullptr, ExcNeedsDoFHandler());

if (level_dof_access)
DoFCellAccessor<dim, spacedim, true>(&cell.value()->get_triangulation(),
cell.value()->level(),
cell.value()->index(),
dof_handler)
.get_interpolated_dof_values(in, out);
else
DoFCellAccessor<dim, spacedim, false>(&cell.value()->get_triangulation(),
cell.value()->level(),
cell.value()->index(),
dof_handler)
.get_interpolated_dof_values(in, out);

switch (cell.value().index())
{
case 1:
std::get<1>(cell.value())->get_interpolated_dof_values(in, out);
break;

case 2:
std::get<2>(cell.value())->get_interpolated_dof_values(in, out);
break;

default:
Assert(false, ExcNeedsDoFHandler());
break;
}
}


Expand All @@ -192,21 +214,32 @@ FEValuesBase<dim, spacedim>::CellIteratorContainer::get_interpolated_dof_values(
Vector<IndexSet::value_type> &out) const
{
Assert(is_initialized(), ExcNotReinited());
Assert(dof_handler != nullptr, ExcNeedsDoFHandler());
Assert(level_dof_access == false, ExcNotImplemented());

const DoFCellAccessor<dim, spacedim, false> cell_dofs(
&cell.value()->get_triangulation(),
cell.value()->level(),
cell.value()->index(),
dof_handler);
switch (cell.value().index())
{
case 1:
{
const typename DoFHandler<dim, spacedim>::cell_iterator cell =
std::get<1>(this->cell.value());

std::vector<types::global_dof_index> dof_indices(
cell->get_fe().n_dofs_per_cell());

cell->get_dof_indices(dof_indices);

std::vector<types::global_dof_index> dof_indices(
cell_dofs.get_fe().n_dofs_per_cell());
cell_dofs.get_dof_indices(dof_indices);
for (unsigned int i = 0; i < cell->get_fe().n_dofs_per_cell(); ++i)
out[i] = (in.is_element(dof_indices[i]) ? 1 : 0);

for (unsigned int i = 0; i < cell_dofs.get_fe().n_dofs_per_cell(); ++i)
out[i] = (in.is_element(dof_indices[i]) ? 1 : 0);
break;
}
case 2:
Assert(false, ExcNotImplemented());
break;

default:
Assert(false, ExcNeedsDoFHandler());
break;
}
}


Expand Down