From d557770c5ad49587e9b607750a3b7f6a3bcc1a9e Mon Sep 17 00:00:00 2001 From: Luca Heltai Date: Wed, 2 Aug 2023 02:22:34 +0200 Subject: [PATCH] FEValuesViews::ReorderedView --- doc/news/changes/minor/20230802LucaHeltai | 4 + include/deal.II/fe/fe_coupling_values.h | 667 ++++++++++++++++++ tests/fe/fe_values_views_renumbered_01.cc | 148 ++++ tests/fe/fe_values_views_renumbered_01.output | 25 + tests/fe/fe_values_views_renumbered_02.cc | 148 ++++ tests/fe/fe_values_views_renumbered_02.output | 9 + tests/fe/fe_values_views_renumbered_03.cc | 149 ++++ tests/fe/fe_values_views_renumbered_03.output | 9 + 8 files changed, 1159 insertions(+) create mode 100644 doc/news/changes/minor/20230802LucaHeltai create mode 100644 include/deal.II/fe/fe_coupling_values.h create mode 100644 tests/fe/fe_values_views_renumbered_01.cc create mode 100644 tests/fe/fe_values_views_renumbered_01.output create mode 100644 tests/fe/fe_values_views_renumbered_02.cc create mode 100644 tests/fe/fe_values_views_renumbered_02.output create mode 100644 tests/fe/fe_values_views_renumbered_03.cc create mode 100644 tests/fe/fe_values_views_renumbered_03.output diff --git a/doc/news/changes/minor/20230802LucaHeltai b/doc/news/changes/minor/20230802LucaHeltai new file mode 100644 index 0000000000..428bf6c29e --- /dev/null +++ b/doc/news/changes/minor/20230802LucaHeltai @@ -0,0 +1,4 @@ +New: The new class FEValuesViews::RenumberedView allows one to filter an existing +FEValuesViews object via two renumbering vectors, one acting on the degrees of +freedom, and the other acting on the quadrature points. +
(Luca Heltai, 2023/08/02) diff --git a/include/deal.II/fe/fe_coupling_values.h b/include/deal.II/fe/fe_coupling_values.h new file mode 100644 index 0000000000..c6ab804fec --- /dev/null +++ b/include/deal.II/fe/fe_coupling_values.h @@ -0,0 +1,667 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2023 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +#ifndef dealii_fe_coupling_values_h +#define dealii_fe_coupling_values_h + +#include + +#include + +#include +#include +#include + +#include + +DEAL_II_NAMESPACE_OPEN + +// Forward declaration +#ifndef DOXYGEN +template +class FEValuesBase; +#endif + +namespace FEValuesViews +{ + /** + * A class that stores the renumbering of degrees of freedom and quadrature + * points for the RenumberedView class. This class is common to all possible + * renumbered views, and therefore it can be shared between all views that use + * the same set of renumbering vectors. + * + * The renumbering is stored in two vectors, one for the degrees of freedom + * and one for the quadrature points, and it is used by the + * FEValuesViews::RenumberedView class. + */ + struct RenumberingData + { + /** + * Construct a new renumbering data object. + * + * The @p dof_renumbering vector is used to renumber the degrees of freedom, + * while the @p quadrature_renumbering vector is used to renumber the + * quadrature points. An empty renumbering vector simply means that no + * renumbering is performed. + * + * @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, or to run over some of the degrees of freedom more than once with + * regard to the underlying view. The important part is that every entry of + * each renumbering vector is a legal index within 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 example of the renumbering ordering is the following. Assume that you + * have an FE_Q(1) finite element in dimension two (i.e., four degrees of + * freedom), and that you are manually taking care of an additional degree + * of freedom (say, a custom shape function in the middle of the 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. Then the renumbering data + * object should look like this: + * @code + * ... + * RenumberingData data({{0, 1, numbers::invalid_unsigned_int, 2, 3}}); + * @endcode + * + * Using this RenumberingData object, the RenumberedView will return zero + * when we ask for values, gradients, etc., of the shape function with index + * two. + */ + RenumberingData( + const unsigned int n_inner_dofs = numbers::invalid_unsigned_int, + const unsigned int n_inner_quadrature_points = + numbers::invalid_unsigned_int, + const std::vector &dof_renumbering = {}, + const std::vector &quadrature_renumbering = {}); + + /** + * Make sure we never copy this object, for performance reasons. + */ + RenumberingData(const RenumberingData &other) = delete; + + /** + * The number of dofs in the underlying view (before any renumbering). + */ + const unsigned int n_inner_dofs; + + /** + * The number of dofs in the renumbered view. + */ + const unsigned int n_dofs; + + /** + * The number of quadrature points in the underlying view (before any + * renumbering). + */ + const unsigned int n_inner_quadrature_points; + + /** + * The number of quadrature points in the renumbered view. + */ + const unsigned int n_quadrature_points; + + /** + * The renumbering of degrees of freedom. + */ + const std::vector dof_renumbering; + + /** + * The renumbering of quadrature points. + */ + const std::vector quadrature_renumbering; + + /** + * General data storage to store temporary vectors. + * + * When the renumbering vectors are non empty, the RenumberedView 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 beforehand 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. + */ + mutable Threads::ThreadLocalStorage data_storage; + }; + + /** + * 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, 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, to select only a subset of the + * degrees of freedom, or to combine two different sets of degrees of freedom + * together. This class provides a view to a given FEValuesViews object, where + * the degrees of freedom and quadrature points are renumbered according to + * the given RenumberingData object (see there for a documentation of how the + * renumbering is interpreted). + * + * 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. + */ + template + class RenumberedView + { + public: + /** + * An alias for the data type of values of the underlying view. + */ + using value_type = typename ViewType::value_type; + + /** + * An alias for the data type of gradients of the underlying view. + */ + using gradient_type = typename ViewType::gradient_type; + + /** + * An alias for the data type of the product of a @p Number and the values + * of the underlying view type. + */ + template + using solution_value_type = + typename ViewType::template solution_value_type; + + /** + * An alias for the data type of the product of a @p Number and the gradients + * of the underlying view type. + */ + template + using solution_gradient_type = + typename ViewType::template solution_gradient_type; + + /** + * Construct a new RenumberedView object. + * + * The renumbering information is taken from the class RenumberingData (see + * there for a documentation of how the renumbering is interpreted). + * + * @note The renumbering information is stored as a reference in this + * object, so you have to make sure that the RenumberingData object lives + * longer than this object. + * + * @param view The underlying FEValuesViews object. + * @param data A RenumberingData object, containing renumbering information. + */ + RenumberedView(const ViewType &view, const RenumberingData &data); + + /** + * Return the value of the underlying view, for the shape function and + * quadrature point selected by the arguments. + * + * @param shape_function Number of the shape function to be evaluated. Note + * that this number runs from zero to size of the renumbering vector + * provided at construction time, or to `dofs_per_cell`, if the renumbering + * vector is empty. + * + * @param q_point Number of the quadrature point at which the function is to + * be evaluated. Note that this number runs from zero to the size of the + * renumbering vector provided at construction time, or to + * `n_quadrature_points`, if the renumbering vector is empty. + * + * @dealiiRequiresUpdateFlags{update_values} + */ + value_type + value(const unsigned int shape_function, const unsigned int q_point) const; + + /** + * Return the gradient of the underlying view, for the shape function and + * quadrature point selected by the arguments. + * + * @param shape_function Number of the shape function to be evaluated. This + * number runs from zero to size of the renumbering vector provided at + * construction time, or to `dofs_per_cell`, if the renumbering vector is + * empty. + * + * @param q_point Number of the quadrature point at which the function is to + * be evaluated. This number runs from zero to the size of the renumbering + * vector provided at construction time, or to `n_quadrature_points`, if the + * renumbering vector is empty. + * + * @dealiiRequiresUpdateFlags{update_gradients} + */ + gradient_type + gradient(const unsigned int shape_function, + const unsigned int q_point) const; + + /** + * Return the values of the underlying view characterized by + * fe_function at the renumbered quadrature points. + * + * This function is the equivalent of the FEValuesBase::get_function_values + * function but it only works on the selected view. + * + * The data type stored by the output vector must be what you get when you + * multiply the values of shape functions (i.e., @p value_type) times the + * type used to store the values of the unknowns $U_j$ of your finite + * element vector $U$ (represented by the @p fe_function argument). + * + * @dealiiRequiresUpdateFlags{update_values} + */ + template + void + get_function_values(const ReadVector &fe_function, + std::vector> &values) const; + + /** + * Same as above, but using a vector of renumbered local degree-of-freedom + * values. In other words, instead of extracting the nodal values of the + * degrees of freedom located on the current cell from a global vector + * associated with a DoFHandler object (as the function above does), this + * function instead takes these local nodal values through its first + * argument. + * + * @param[in] dof_values A vector of local nodal values. This vector must + * have a length equal to the size of the dof renumbering vector. + * + * @param[out] values A vector of values of the given finite element field, + * at the renumbered quadrature points on the current object. + * + * @tparam InputVector The @p InputVector type must allow creation of an + * ArrayView object from it; this is satisfied by the `std::vector` class, + * among others. + * + * @dealiiRequiresUpdateFlags{update_values} + */ + template + void + get_function_values_from_local_dof_values( + const InputVector &dof_values, + std::vector> + &values) const; + + /** + * Return the gradients of the underlying view characterized by + * fe_function at the renumbered quadrature points. + * + * This function is the equivalent of the + * FEValuesBase::get_function_gradients function but it only works on the + * selected view. + * + * The data type stored by the output vector must be what you get when you + * multiply the gradients of shape functions (i.e., @p value_type) times the + * type used to store the gradients of the unknowns $U_j$ of your finite + * element vector $U$ (represented by the @p fe_function argument). + * + * @dealiiRequiresUpdateFlags{update_gradients} + */ + template + void + get_function_gradients( + const ReadVector &fe_function, + std::vector> &gradients) const; + + /** + * Same as above, but using a vector of renumbered local degree-of-freedom + * gradients. In other words, instead of extracting the nodal values of the + * degrees of freedom located on the current cell from a global vector + * associated with a DoFHandler object (as the function above does), this + * function instead takes these local nodal values through its first + * argument. + * + * @param[in] dof_values A vector of local nodal values. This vector must + * have a length equal to the size of the dof renumbering vector. + * + * @param[out] gradients A vector of gradients of the given finite element + * field, at the renumbered quadrature points on the current object. + * + * @tparam InputVector The @p InputVector type must allow creation of an + * ArrayView object from it; this is satisfied by the `std::vector` class, + * among others. + * + * @dealiiRequiresUpdateFlags{update_gradients} + */ + template + void + get_function_gradients_from_local_dof_values( + const InputVector &dof_values, + std::vector> + &gradients) const; + + private: + /** + * The data structure that stores the renumbering of the degrees of freedom + * and of the quadrature points. + */ + const RenumberingData &data; + + + /** + * 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 + std::string + get_unique_container_name(const std::string &prefix, + const unsigned int size, + const Number &exemplar_number) const; + + /** + * Produce an inner vector compatible with the inner view, after copying the + * values with the correct numbering from the outer vector. + */ + template + const InputVector & + outer_to_inner_dofs(const InputVector &outer_vector) const; + + /** + * Produce an inner vector compatible with the inner view, and zero out its + * entries if necessary. + */ + template + std::vector & + outer_to_inner_values(std::vector &outer_values) const; + + /** + * Return the outer argument renumbered according to the quadrature + * renumbering. The values in the inner values are copied to the right + * position in the outer vector. + */ + template + void + inner_to_outer_values(const std::vector &inner_values, + std::vector &outer_values) const; + + /** + * Store a reference to the underlying view. + */ + const ViewType &view; + }; +} // namespace FEValuesViews + + +#ifndef DOXYGEN + + +/*------------------------ Inline functions: namespace FEValuesViews --------*/ + +namespace FEValuesViews +{ + RenumberingData::RenumberingData( + const unsigned int n_inner_dofs, + const unsigned int n_inner_quadrature_points, + const std::vector &dof_renumbering, + const std::vector &quadrature_renumbering) + : n_inner_dofs(n_inner_dofs) + , n_dofs(dof_renumbering.empty() ? n_inner_dofs : dof_renumbering.size()) + , n_inner_quadrature_points(n_inner_quadrature_points) + , n_quadrature_points(quadrature_renumbering.empty() ? + n_inner_quadrature_points : + quadrature_renumbering.size()) + , dof_renumbering(dof_renumbering) + , quadrature_renumbering(quadrature_renumbering) + { +// Check that the renumbering vectors are valid. +# ifdef DEBUG + // While for dofs we admit invalid values, this is not the case for + // quadrature points. + for (const auto i : dof_renumbering) + Assert(i < n_inner_dofs || i == numbers::invalid_unsigned_int, + ExcIndexRange(i, 0, n_inner_dofs)); + + for (const auto q : quadrature_renumbering) + AssertIndexRange(q, n_inner_quadrature_points); +# endif + } + + + + template + RenumberedView::RenumberedView(const ViewType &view, + const RenumberingData &data) + : view(view) + , data(data) + {} + + + + template + template + inline std::string + RenumberedView::get_unique_container_name( + const std::string &prefix, + const unsigned int size, + const Number &exemplar_number) const + { + return prefix + "_" + Utilities::int_to_string(size) + "_" + + Utilities::type_to_string(exemplar_number); + } + + + + template + typename RenumberedView::value_type + RenumberedView::value(const unsigned int shape_function, + const unsigned int q_point) const + { + AssertIndexRange(shape_function, data.n_dofs); + AssertIndexRange(q_point, data.n_quadrature_points); + + const auto inner_shape_function = data.dof_renumbering.empty() ? + shape_function : + data.dof_renumbering[shape_function]; + const auto inner_q_point = data.quadrature_renumbering.empty() ? + q_point : + data.quadrature_renumbering[q_point]; + if (inner_shape_function == numbers::invalid_unsigned_int) + return value_type(0); + else + { + AssertIndexRange(inner_shape_function, data.n_inner_dofs); + AssertIndexRange(inner_q_point, data.n_inner_quadrature_points); + return view.value(inner_shape_function, inner_q_point); + } + } + + + + template + typename RenumberedView::gradient_type + RenumberedView::gradient(const unsigned int shape_function, + const unsigned int q_point) const + { + AssertIndexRange(shape_function, data.n_dofs); + AssertIndexRange(q_point, data.n_quadrature_points); + + const auto inner_shape_function = data.dof_renumbering.empty() ? + shape_function : + data.dof_renumbering[shape_function]; + const auto inner_q_point = data.quadrature_renumbering.empty() ? + q_point : + data.quadrature_renumbering[q_point]; + if (inner_shape_function == numbers::invalid_unsigned_int) + return gradient_type(); + else + return view.gradient(inner_shape_function, inner_q_point); + } + + + + template + template + std::vector & + RenumberedView::outer_to_inner_values( + std::vector &outer_values) const + { + AssertDimension(outer_values.size(), data.n_quadrature_points); + if (data.quadrature_renumbering.empty()) + { + return outer_values; + } + else + { + const auto name = + get_unique_container_name("RenumberedView::outer_to_inner_values", + data.n_inner_quadrature_points, + outer_values[0]); + auto &inner_values = + data.data_storage.get() + .template get_or_add_object_with_name>( + name, data.n_inner_quadrature_points); + return inner_values; + } + } + + + + template + template + const VectorType & + RenumberedView::outer_to_inner_dofs( + const VectorType &outer_dofs) const + { + AssertDimension(outer_dofs.size(), data.n_dofs); + if (data.dof_renumbering.empty()) + { + return outer_dofs; + } + else + { + const auto name = + get_unique_container_name("RenumberedView::outer_to_inner_dofs", + data.n_inner_dofs, + outer_dofs[0]); + + auto &inner_dofs = data.data_storage.get() + .template get_or_add_object_with_name( + name, data.n_inner_dofs); + for (unsigned int i = 0; i < data.n_dofs; ++i) + { + const auto inner_i = data.dof_renumbering[i]; + if (inner_i != numbers::invalid_unsigned_int) + { + AssertIndexRange(inner_i, data.n_inner_dofs); + inner_dofs[inner_i] = outer_dofs[i]; + } + } + return inner_dofs; + } + } + + + + template + template + void + RenumberedView::inner_to_outer_values( + const std::vector &inner_values, + std::vector &outer_values) const + { + AssertDimension(outer_values.size(), data.n_quadrature_points); + AssertDimension(inner_values.size(), data.n_inner_quadrature_points); + if (data.quadrature_renumbering.empty()) + { + Assert(&inner_values == &outer_values, ExcInternalError()); + return; + } + for (unsigned int i = 0; i < data.quadrature_renumbering.size(); ++i) + { + outer_values[i] = inner_values[data.quadrature_renumbering[i]]; + } + } + + + + template + template + void + RenumberedView::get_function_values( + const ReadVector &fe_function, + std::vector> &values) const + { + auto &inner_values = outer_to_inner_values(values); + view.get_function_values(fe_function, inner_values); + inner_to_outer_values(inner_values, values); + } + + + + template + template + void + RenumberedView::get_function_values_from_local_dof_values( + const InputVector &dof_values, + std::vector> &values) + const + { + const auto &inner_dof_values = outer_to_inner_dofs(dof_values); + auto &inner_values = outer_to_inner_values(values); + + view.get_function_values_from_local_dof_values(inner_dof_values, + inner_values); + inner_to_outer_values(inner_values, values); + } + + + template + template + void + RenumberedView::get_function_gradients( + const ReadVector &fe_function, + std::vector> &gradients) const + { + auto &inner_gradients = outer_to_inner_values(gradients); + view.get_function_gradients(fe_function, inner_gradients); + inner_to_outer_values(inner_gradients, gradients); + } + + + + template + template + void + RenumberedView::get_function_gradients_from_local_dof_values( + const InputVector &dof_values, + std::vector> + &gradients) const + { + const auto &inner_dof_values = outer_to_inner_dofs(dof_values); + auto &inner_gradients = outer_to_inner_values(gradients); + + view.get_function_gradients_from_local_dof_values(inner_dof_values, + inner_gradients); + inner_to_outer_values(inner_gradients, gradients); + } +} // namespace FEValuesViews + +#endif // DOXYGEN + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/tests/fe/fe_values_views_renumbered_01.cc b/tests/fe/fe_values_views_renumbered_01.cc new file mode 100644 index 0000000000..f4e69e7cc4 --- /dev/null +++ b/tests/fe/fe_values_views_renumbered_01.cc @@ -0,0 +1,148 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2023 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + +// Test renumbering of FEValuesViews::Scalar using +// FEValuesViews::RenumberedView + +#include + +#include +#include +#include +#include + +#include +#include +#include + +#include +#include + +#include "../tests.h" + +#include "../test_grids.h" + +template +void +test() +{ + Triangulation tria; + TestGrids::hyper_line(tria, 2); + + DoFHandler dofh(tria); + FE_Q fe(1); + dofh.distribute_dofs(fe); + + MappingQ mapping(1); + UpdateFlags update_flags = + update_values | update_quadrature_points | update_JxW_values; + + FEValues fv1(mapping, fe, QGauss(fe.degree + 1), update_flags); + fv1.reinit(dofh.begin_active()); + + FEValuesViews::Scalar fv1_scalar(fv1, 0); + + using iota = std_cxx20::ranges::iota_view; + { + // No renumbering + for (const auto &i : fv1.dof_indices()) + { + deallog << "i = " << i << ", fv1_scalar = " << fv1_scalar.value(i, 0) + << std::endl; + } + } + { + // Trivial renumbering + auto id = iota(0, fv1.dofs_per_cell); + std::vector renumber(id.begin(), id.end()); + FEValuesViews::RenumberingData data(fv1.dofs_per_cell, + fv1.n_quadrature_points, + renumber); + FEValuesViews::RenumberedView> + fv1_scalar_reordered(fv1_scalar, data); + + for (const auto &i : fv1.dof_indices()) + { + deallog << "i = " << i << ", fv1_scalar_trivial_renumber = " + << fv1_scalar_reordered.value(i, 0) << std::endl; + } + } + { + // Inverse renumbering + auto id = iota(0, fv1.dofs_per_cell); + std::vector renumber(id.begin(), id.end()); + std::reverse(renumber.begin(), renumber.end()); + FEValuesViews::RenumberingData data(fv1.dofs_per_cell, + fv1.n_quadrature_points, + renumber); + FEValuesViews::RenumberedView> + fv1_scalar_reordered(fv1_scalar, data); + for (const auto &i : fv1.dof_indices()) + { + deallog << "i = " << i << ", fv1_scalar_inverse_renumber = " + << fv1_scalar_reordered.value(i, 0) << std::endl; + } + } + // Now test quadrature renumbering + { + // No renumbering + for (const auto &q : fv1.quadrature_point_indices()) + { + deallog << "q = " << q + << ", fv1_scalar_0_q = " << fv1_scalar.value(0, q) << std::endl; + } + } + { + // Trivial renumbering + auto id = iota(0, fv1.n_quadrature_points); + std::vector renumber(id.begin(), id.end()); + FEValuesViews::RenumberingData data(fv1.dofs_per_cell, + fv1.n_quadrature_points, + {}, + renumber); + FEValuesViews::RenumberedView> + fv1_scalar_reordered(fv1_scalar, data); + for (const auto &q : fv1.quadrature_point_indices()) + { + deallog << "q = " << q << ", fv1_scalar_0_q_trivial_renumbering = " + << fv1_scalar_reordered.value(0, q) << std::endl; + } + } + { + // Inverse renumbering + auto id = iota(0, fv1.n_quadrature_points); + std::vector renumber(id.begin(), id.end()); + std::reverse(renumber.begin(), renumber.end()); + FEValuesViews::RenumberingData data(fv1.dofs_per_cell, + fv1.n_quadrature_points, + {}, + renumber); + FEValuesViews::RenumberedView> + fv1_scalar_reordered(fv1_scalar, data); + for (const auto &q : fv1.quadrature_point_indices()) + { + deallog << "q = " << q << ", fv1_scalar_0_q_inverse_renumbering = " + << fv1_scalar_reordered.value(0, q) << std::endl; + } + } +} + +int +main() +{ + initlog(); + test<2>(); +} diff --git a/tests/fe/fe_values_views_renumbered_01.output b/tests/fe/fe_values_views_renumbered_01.output new file mode 100644 index 0000000000..86d8fb06ac --- /dev/null +++ b/tests/fe/fe_values_views_renumbered_01.output @@ -0,0 +1,25 @@ + +DEAL::i = 0, fv1_scalar = 0.622008 +DEAL::i = 1, fv1_scalar = 0.166667 +DEAL::i = 2, fv1_scalar = 0.166667 +DEAL::i = 3, fv1_scalar = 0.0446582 +DEAL::i = 0, fv1_scalar_trivial_renumber = 0.622008 +DEAL::i = 1, fv1_scalar_trivial_renumber = 0.166667 +DEAL::i = 2, fv1_scalar_trivial_renumber = 0.166667 +DEAL::i = 3, fv1_scalar_trivial_renumber = 0.0446582 +DEAL::i = 0, fv1_scalar_inverse_renumber = 0.0446582 +DEAL::i = 1, fv1_scalar_inverse_renumber = 0.166667 +DEAL::i = 2, fv1_scalar_inverse_renumber = 0.166667 +DEAL::i = 3, fv1_scalar_inverse_renumber = 0.622008 +DEAL::q = 0, fv1_scalar_0_q = 0.622008 +DEAL::q = 1, fv1_scalar_0_q = 0.166667 +DEAL::q = 2, fv1_scalar_0_q = 0.166667 +DEAL::q = 3, fv1_scalar_0_q = 0.0446582 +DEAL::q = 0, fv1_scalar_0_q_trivial_renumbering = 0.622008 +DEAL::q = 1, fv1_scalar_0_q_trivial_renumbering = 0.166667 +DEAL::q = 2, fv1_scalar_0_q_trivial_renumbering = 0.166667 +DEAL::q = 3, fv1_scalar_0_q_trivial_renumbering = 0.0446582 +DEAL::q = 0, fv1_scalar_0_q_inverse_renumbering = 0.0446582 +DEAL::q = 1, fv1_scalar_0_q_inverse_renumbering = 0.166667 +DEAL::q = 2, fv1_scalar_0_q_inverse_renumbering = 0.166667 +DEAL::q = 3, fv1_scalar_0_q_inverse_renumbering = 0.622008 diff --git a/tests/fe/fe_values_views_renumbered_02.cc b/tests/fe/fe_values_views_renumbered_02.cc new file mode 100644 index 0000000000..4129ffc4f0 --- /dev/null +++ b/tests/fe/fe_values_views_renumbered_02.cc @@ -0,0 +1,148 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2023 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + +// Test renumbering of FEValuesViews::Scalar using +// FEValuesViews::RenumberedView. Test get_function_values and +// get_function_values_from_local_dof_values + +#include +#include + +#include +#include +#include + +#include +#include +#include + +#include + +#include +#include + +#include "../tests.h" + +#include "../test_grids.h" + +template +void +test() +{ + Triangulation tria; + TestGrids::hyper_line(tria, 2); + + DoFHandler dofh(tria); + FE_Q fe(1); + dofh.distribute_dofs(fe); + + // Create a non trivial function of x and y + FunctionParser f("x+10*y"); + + Vector vec(dofh.n_dofs()); + VectorTools::interpolate(dofh, f, vec); + + MappingQ mapping(1); + UpdateFlags update_flags = + update_values | update_quadrature_points | update_JxW_values; + + FEValues fe_v(mapping, fe, QGauss(fe.degree + 1), update_flags); + const auto &cell = dofh.begin_active(); + fe_v.reinit(cell); + + const auto &fv1_scalar = fe_v[FEValuesExtractors::Scalar(0)]; + + // Small print function + auto print = [](const auto &str, const auto &v) { + deallog << str << " = "; + for (const auto &i : v) + deallog << i << " "; + deallog << std::endl; + }; + + // Now build an artificial renumbering vector for dofs and one for quadratures + // - The first renumbering, simply shifts the dof indices by 1, as if there + // was an additional dof at the beginning, which is ignored by this view + // - The second renumbering, repeats each quadrature point twice, so that the + // returned values are actually repeated. + using iota = std_cxx20::ranges::iota_view; + + std::vector dof_renumbering(fe_v.dofs_per_cell + 1); + std::vector quad_renumbering(fe_v.n_quadrature_points * 2); + { + auto id = iota(0, fe_v.dofs_per_cell); + dof_renumbering[0] = numbers::invalid_unsigned_int; + std::copy(id.begin(), id.end(), dof_renumbering.begin() + 1); + } + { + auto id = iota(0, fe_v.n_quadrature_points); + std::copy(id.begin(), id.end(), quad_renumbering.begin()); + std::copy(id.begin(), id.end(), quad_renumbering.begin() + id.size()); + } + + deallog << "dof_renumbering = "; + for (const auto &i : dof_renumbering) + deallog << (int)i << " "; + deallog << std::endl; + print("quad_renumbering", quad_renumbering); + + std::vector function_values(fe_v.n_quadrature_points); + std::vector function_values_renumbered(fe_v.n_quadrature_points * 2); + + fv1_scalar.get_function_values(vec, function_values); + print("function_values", function_values); + + FEValuesViews::RenumberingData data(fe_v.dofs_per_cell, + fe_v.n_quadrature_points, + dof_renumbering, + quad_renumbering); + FEValuesViews::RenumberedView> + fv1_scalar_renumbered(fv1_scalar, data); + + fv1_scalar_renumbered.get_function_values(vec, function_values_renumbered); + print("function_values_renumbered", function_values_renumbered); + + std::vector dof_values(fe_v.dofs_per_cell); + std::vector dof_values_renumbered(fe_v.dofs_per_cell + 1); + + cell->get_dof_values(vec, dof_values.begin(), dof_values.end()); + + // Now fill the first dof with a dummy value + dof_values_renumbered[0] = 42; + std::copy(dof_values.begin(), + dof_values.end(), + dof_values_renumbered.begin() + 1); + + print("dof_values", dof_values); + print("dof_values_renumbered", dof_values_renumbered); + + fv1_scalar.get_function_values_from_local_dof_values(dof_values, + function_values); + print("function_values from dof_values", function_values); + + fv1_scalar_renumbered.get_function_values_from_local_dof_values( + dof_values_renumbered, function_values_renumbered); + + print("function_values_renumbered from dof_values_renumbered", + function_values_renumbered); +} + +int +main() +{ + initlog(); + test<2>(); +} diff --git a/tests/fe/fe_values_views_renumbered_02.output b/tests/fe/fe_values_views_renumbered_02.output new file mode 100644 index 0000000000..78fc11cdd4 --- /dev/null +++ b/tests/fe/fe_values_views_renumbered_02.output @@ -0,0 +1,9 @@ + +DEAL::dof_renumbering = -1 0 1 2 3 +DEAL::quad_renumbering = 0 1 2 3 0 1 2 3 +DEAL::function_values = 2.32457 2.90192 8.09808 8.67543 +DEAL::function_values_renumbered = 2.32457 2.90192 8.09808 8.67543 2.32457 2.90192 8.09808 8.67543 +DEAL::dof_values = 0.00000 1.00000 10.0000 11.0000 +DEAL::dof_values_renumbered = 42.0000 0.00000 1.00000 10.0000 11.0000 +DEAL::function_values from dof_values = 2.32457 2.90192 8.09808 8.67543 +DEAL::function_values_renumbered from dof_values_renumbered = 2.32457 2.90192 8.09808 8.67543 2.32457 2.90192 8.09808 8.67543 diff --git a/tests/fe/fe_values_views_renumbered_03.cc b/tests/fe/fe_values_views_renumbered_03.cc new file mode 100644 index 0000000000..0b58974373 --- /dev/null +++ b/tests/fe/fe_values_views_renumbered_03.cc @@ -0,0 +1,149 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2023 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + +// Test renumbering of FEValuesViews::Scalar using +// FEValuesViews::RenumberedView. Test get_function_gradients and +// get_function_gradients_from_local_dof_values + +#include +#include + +#include +#include +#include + +#include +#include +#include + +#include + +#include +#include + +#include "../tests.h" + +#include "../test_grids.h" + +template +void +test() +{ + Triangulation tria; + TestGrids::hyper_line(tria, 2); + + DoFHandler dofh(tria); + FE_Q fe(2); + dofh.distribute_dofs(fe); + + // Create a non trivial function of x and y + FunctionParser f("(x^2)/2+(10*y^2)/2"); + + Vector vec(dofh.n_dofs()); + VectorTools::interpolate(dofh, f, vec); + + MappingQ mapping(1); + UpdateFlags update_flags = update_gradients; + + FEValues fe_v(mapping, fe, QGauss(1), update_flags); + const auto &cell = dofh.begin_active(); + fe_v.reinit(cell); + + const auto &fv1_scalar = fe_v[FEValuesExtractors::Scalar(0)]; + + // Small print function + auto print = [](const auto &str, const auto &v) { + deallog << str << " = "; + for (const auto &i : v) + deallog << i << "; "; + deallog << std::endl; + }; + + // Now build an artificial renumbering vector for dofs and one for quadratures + // - The first renumbering, simply shifts the dof indices by 1, as if there + // was an additional dof at the beginning, which is ignored by this view + // - The second renumbering, repeats each quadrature point twice, so that the + // returned values are actually repeated. + using iota = std_cxx20::ranges::iota_view; + + std::vector dof_renumbering(fe_v.dofs_per_cell + 1); + std::vector quad_renumbering(fe_v.n_quadrature_points * 2); + { + auto id = iota(0, fe_v.dofs_per_cell); + dof_renumbering[0] = numbers::invalid_unsigned_int; + std::copy(id.begin(), id.end(), dof_renumbering.begin() + 1); + } + { + auto id = iota(0, fe_v.n_quadrature_points); + std::copy(id.begin(), id.end(), quad_renumbering.begin()); + std::copy(id.begin(), id.end(), quad_renumbering.begin() + id.size()); + } + + deallog << "dof_renumbering = "; + for (const auto &i : dof_renumbering) + deallog << (int)i << " "; + deallog << std::endl; + print("quad_renumbering", quad_renumbering); + + std::vector> function_gradients(fe_v.n_quadrature_points); + std::vector> function_gradients_renumbered( + fe_v.n_quadrature_points * 2); + + fv1_scalar.get_function_gradients(vec, function_gradients); + print("function_gradients", function_gradients); + + FEValuesViews::RenumberingData data(fe_v.dofs_per_cell, + fe_v.n_quadrature_points, + dof_renumbering, + quad_renumbering); + FEValuesViews::RenumberedView> + fv1_scalar_renumbered(fv1_scalar, data); + + fv1_scalar_renumbered.get_function_gradients(vec, + function_gradients_renumbered); + print("function_gradients_renumbered", function_gradients_renumbered); + + std::vector dof_values(fe_v.dofs_per_cell); + std::vector dof_values_renumbered(fe_v.dofs_per_cell + 1); + + cell->get_dof_values(vec, dof_values.begin(), dof_values.end()); + + // Now fill the first dof with a dummy value + dof_values_renumbered[0] = 42; + std::copy(dof_values.begin(), + dof_values.end(), + dof_values_renumbered.begin() + 1); + + print("dof_values", dof_values); + print("dof_values_renumbered", dof_values_renumbered); + + fv1_scalar.get_function_gradients_from_local_dof_values(dof_values, + function_gradients); + print("function_gradients from dof_values", function_gradients); + + fv1_scalar_renumbered.get_function_gradients_from_local_dof_values( + dof_values_renumbered, function_gradients_renumbered); + + print("function_gradients_renumbered from dof_values_renumbered", + function_gradients_renumbered); +} + +int +main() +{ + initlog(); + test<2>(); +} diff --git a/tests/fe/fe_values_views_renumbered_03.output b/tests/fe/fe_values_views_renumbered_03.output new file mode 100644 index 0000000000..92fce61ad6 --- /dev/null +++ b/tests/fe/fe_values_views_renumbered_03.output @@ -0,0 +1,9 @@ + +DEAL::dof_renumbering = -1 0 1 2 3 4 5 6 7 8 +DEAL::quad_renumbering = 0; 0; +DEAL::function_gradients = 0.500000 5.00000; +DEAL::function_gradients_renumbered = 0.500000 5.00000; 0.500000 5.00000; +DEAL::dof_values = 0.00000; 0.500000; 5.00000; 5.50000; 1.25000; 1.75000; 0.125000; 5.12500; 1.37500; +DEAL::dof_values_renumbered = 42.0000; 0.00000; 0.500000; 5.00000; 5.50000; 1.25000; 1.75000; 0.125000; 5.12500; 1.37500; +DEAL::function_gradients from dof_values = 0.500000 5.00000; +DEAL::function_gradients_renumbered from dof_values_renumbered = 0.500000 5.00000; 0.500000 5.00000; -- 2.39.5