Skip to content

Commit

Permalink
Apply suggestions from code review
Browse files Browse the repository at this point in the history
Co-authored-by: Wolfgang Bangerth <bangerth@colostate.edu>
  • Loading branch information
luca-heltai and bangerth committed Aug 6, 2023
1 parent fda43db commit 16d6a4a
Show file tree
Hide file tree
Showing 2 changed files with 86 additions and 36 deletions.
119 changes: 84 additions & 35 deletions include/deal.II/fe/fe_values_coupling.h
Original file line number Diff line number Diff line change
Expand Up @@ -40,14 +40,15 @@ namespace FEValuesViews
* A class that provides a renumbered view to a given FEValuesViews object.
*
* In general, the order of the degrees of freedom and quadrature points
* follows the one of the FEValues object. However, in some cases, it is
* convenient to group together degrees of freedom and quadrature points in a
* different order. This class provides a view to a given FEValuesViews
* object, where the degrees of freedom and quadrature points are renumbered
* according to the given renumbering vectors.
* follows the one of the FEValues object, which itself uses the numbering
* provided by the FiniteElement and Quadrature objects it uses. However, in
* some cases, it is convenient to group together degrees of freedom and
* quadrature points in a different order. This class provides a view to a
* given FEValuesViews object, where the degrees of freedom and quadrature
* points are renumbered according to the given renumbering vectors.
*/
template <typename ViewType>
class ReorderedView : public Subscriptor
class ReorderedView
{
public:
/**
Expand All @@ -65,15 +66,15 @@ namespace FEValuesViews
* of the underlying view type.
*/
template <typename Number>
using solution_value_type = typename ProductType<Number, value_type>::type;
using solution_value_type = typename ViewType::solution_value_type<Number>;

/**
* An alias for the data type of the product of a @p Number and the gradients
* of the underlying view type.
*/
template <typename Number>
using solution_gradient_type =
typename ProductType<Number, gradient_type>::type;
typename ViewType::solution_gradient_type<Number>;

/**
* Construct a new ReorderedView object.
Expand All @@ -83,15 +84,51 @@ namespace FEValuesViews
* quadrature points. An empty renumbering vector simply means that no
* renumbering is performed.
*
* @note The renumbering vectors may *not* match the dimensions of the
* underlying view, i.e., you could use this class to only run over half of
* the degrees of freedom, and integrating only over half of the current
* cell, if you wish to do so. The important part is that every entry of
* each renumbering vector is a legal value for the underlying view object.
* We allow the dof renumbering vector to contain
* `numbers::invalid_unsigned_int` values, which are ignored, and produce
* zero values, gradients, hessians, etc. for the corresponding shape
* function.
* @note The renumbering vectors are *not* required to match the dimensions
* of the underlying view, i.e., you could use this class to only run over
* half of the degrees of freedom, and integrating only over half of the
* current cell by selecting a subset of quadrature points, if you wish to
* do so. The important part is that every entry of each renumbering vector
* is a legal value for the underlying view object. We allow the dof
* renumbering vector to contain `numbers::invalid_unsigned_int` values,
* which are ignored, and produce zero values, gradients, hessians, etc. for
* the corresponding shape function.
*
* An admittedly not very useful example usage of this class is the
* following. Assume that you have an FE_Q(1) finite element space, and that
* you are manually taking care of an additional degree of freedom (say, a
* custom shape function in the middle of cell), which, for your own
* convenience, should be numbered locally in the middle of the others. You
* would like your loops on degrees of freedom to run over all degrees of
* freedom, and to have the same structure as if the additional degree of
* freedom was not there. This is where this class comes in handy:
* @code
* ...
* FEValues<dim> fe_values(fe, quadrature_formula, update_values); const
* auto n_dofs_per_cell = fe_values.dofs_per_cell; // say this is 4 const
* const std::vector<unsigned int> dof_renumbering{{0, 1,
* numbers::invalid_unsigned_int, 2, 3}};
*
* ...
* // inside the loop over cells fe_values.reinit(cell); const auto
* &fe_values_reordered =
* FEValuesViews::ReorderedView<FEValuesViews::Scalar<dim>>(fe_values,dof_renumbering);
*
* // loop over all indices in [0,5) // do something with the values: auto
* v_i_q = fe_values_reordered.value(i,q);
* @endcode
*
* Notice that in the code above, for `i=2`, the value `v_i_q` will be zero.
*
* Users will typically not use this class directly, but rather pass an
* extractor from the FEValuesExtractors namespace to the FECouplingValues
* class, which returns an object of this type. This is the same mechanism
* used in the FEValues classes, where passing an extractor returns an
* FEValuesViews object, and the user rarely instantiates an object of this
* type.
*
* @note The renumbering vectors are stored by reference in this class, so
* they must remain valid for the lifetime of the object.
*
* @param view The underlying FEValuesViews object.
* @param dof_renumbering A renumbering vector for the degrees of freedom.
Expand All @@ -118,7 +155,7 @@ namespace FEValuesViews
*
* @dealiiRequiresUpdateFlags{update_values}
*/
typename ViewType::value_type
value_type
value(const unsigned int shape_function, const unsigned int q_point) const;

/**
Expand All @@ -137,7 +174,7 @@ namespace FEValuesViews
*
* @dealiiRequiresUpdateFlags{update_gradients}
*/
typename ViewType::gradient_type
gradient_type
gradient(const unsigned int shape_function,
const unsigned int q_point) const;

Expand Down Expand Up @@ -235,6 +272,18 @@ namespace FEValuesViews
/**
* Helper function that constructs a unique name for a container, based on a
* string prefix, on its size, and on the type stored in the container.
*
* When the renumbering vectors are non empty, this class may need to
* construct temporary vectors to store the values of the solution and/or of
* its gradients with the sizes given by the underlying view object.
* Unfortunately, we don't know before hand the Number types with which the
* vectors will be instantiated, so we cannot use a simple cache internally,
* and we use a GeneralDataStorage object to avoid allocating memory at each
* call.
*
* This function constructs a unique name for temporary containers that will
* be stored upon the first request in the internal GeneralDataStorage
* object, to access the
*/
template <typename Number>
std::string
Expand Down Expand Up @@ -288,7 +337,7 @@ namespace FEValuesViews
/**
* Number of dofs in the underlying view (before any renumbering).
*/
const unsigned int inner_n_dofs;
const unsigned int n_inner_dofs;

/**
* Number of dofs in the renumbered view.
Expand All @@ -299,7 +348,7 @@ namespace FEValuesViews
* Number of quadrature points in the underlying view (before any
* renumbering).
*/
const unsigned int inner_n_q_points;
const unsigned int n_inner_q_points;

/**
* Number of quadrature points in the renumbered view.
Expand Down Expand Up @@ -329,10 +378,10 @@ namespace FEValuesViews
: view(view)
, dof_renumbering(dof_renumbering)
, quadrature_renumbering(quadrature_renumbering)
, inner_n_dofs(view.fe_values->dofs_per_cell)
, n_dofs(dof_renumbering.empty() ? inner_n_dofs : dof_renumbering.size())
, inner_n_q_points(view.fe_values->n_quadrature_points)
, n_q_points(quadrature_renumbering.empty() ? inner_n_q_points :
, n_inner_dofs(view.fe_values->dofs_per_cell)
, n_dofs(dof_renumbering.empty() ? n_inner_dofs : dof_renumbering.size())
, n_inner_q_points(view.fe_values->n_quadrature_points)
, n_q_points(quadrature_renumbering.empty() ? n_inner_q_points :
quadrature_renumbering.size())
{}

Expand Down Expand Up @@ -370,8 +419,8 @@ namespace FEValuesViews
return value_type(0);
else
{
AssertIndexRange(inner_shape_function, inner_n_dofs);
AssertIndexRange(inner_q_point, inner_n_q_points);
AssertIndexRange(inner_shape_function, n_inner_dofs);
AssertIndexRange(inner_q_point, n_inner_q_points);
return view.value(inner_shape_function, inner_q_point);
}
}
Expand All @@ -396,8 +445,8 @@ namespace FEValuesViews
return gradient_type();
else
{
AssertIndexRange(inner_shape_function, inner_n_dofs);
AssertIndexRange(inner_q_point, inner_n_q_points);
AssertIndexRange(inner_shape_function, n_inner_dofs);
AssertIndexRange(inner_q_point, n_inner_q_points);
return view.gradient(inner_shape_function, inner_q_point);
}
}
Expand All @@ -419,12 +468,12 @@ namespace FEValuesViews
{
const auto name =
get_unique_container_name("ReorderedView::outer_to_inner_values",
inner_n_q_points,
n_inner_q_points,
outer_values[0]);
auto &inner_values =
internal_data_storage
.template get_or_add_object_with_name<std::vector<ValueType>>(
name, inner_n_q_points);
name, n_inner_q_points);
return inner_values;
}
}
Expand All @@ -445,18 +494,18 @@ namespace FEValuesViews
{
const auto name =
get_unique_container_name("ReorderedView::outer_to_inner_dofs",
inner_n_dofs,
n_inner_dofs,
outer_dofs[0]);
auto &inner_dofs =
internal_data_storage
.template get_or_add_object_with_name<VectorType>(name,
inner_n_dofs);
n_inner_dofs);
for (unsigned int i = 0; i < n_dofs; ++i)
{
const auto inner_i = dof_renumbering[i];
if (inner_i != numbers::invalid_unsigned_int)
{
AssertIndexRange(inner_i, inner_n_dofs);
AssertIndexRange(inner_i, n_inner_dofs);
inner_dofs[inner_i] = outer_dofs[i];
}
}
Expand All @@ -479,7 +528,7 @@ namespace FEValuesViews
return;
}
AssertDimension(outer_values.size(), n_q_points);
AssertDimension(inner_values.size(), inner_n_q_points);
AssertDimension(inner_values.size(), n_inner_q_points);
for (unsigned int i = 0; i < quadrature_renumbering.size(); ++i)
outer_values[i] = inner_values[quadrature_renumbering[i]];
}
Expand Down
3 changes: 2 additions & 1 deletion tests/fe/fe_values_views_reordered_03.cc
Original file line number Diff line number Diff line change
Expand Up @@ -34,9 +34,10 @@
#include <fstream>
#include <iostream>

#include "../test_grids.h"
#include "../tests.h"

#include "../test_grids.h"

template <int dim>
void
test()
Expand Down

0 comments on commit 16d6a4a

Please sign in to comment.