--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/config.h>
+
+#include <deal.II/algorithms/general_data_storage.h>
+
+#include <deal.II/base/subscriptor.h>
+#include <deal.II/base/thread_local_storage.h>
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/fe/fe_values_views.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+// Forward declaration
+#ifndef DOXYGEN
+template <int dim, int spacedim>
+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<unsigned int> &dof_renumbering = {},
+ const std::vector<unsigned int> &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<unsigned int> dof_renumbering;
+
+ /**
+ * The renumbering of quadrature points.
+ */
+ const std::vector<unsigned int> 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<GeneralDataStorage> 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 <typename ViewType>
+ 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 <typename Number>
+ using solution_value_type =
+ typename ViewType::template 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 ViewType::template solution_gradient_type<Number>;
+
+ /**
+ * 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
+ * <tt>fe_function</tt> 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 <typename Number>
+ void
+ get_function_values(const ReadVector<Number> &fe_function,
+ std::vector<solution_value_type<Number>> &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 <class InputVector>
+ void
+ get_function_values_from_local_dof_values(
+ const InputVector &dof_values,
+ std::vector<solution_value_type<typename InputVector::value_type>>
+ &values) const;
+
+ /**
+ * Return the gradients of the underlying view characterized by
+ * <tt>fe_function</tt> 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 <typename Number>
+ void
+ get_function_gradients(
+ const ReadVector<Number> &fe_function,
+ std::vector<solution_gradient_type<Number>> &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 <class InputVector>
+ void
+ get_function_gradients_from_local_dof_values(
+ const InputVector &dof_values,
+ std::vector<solution_gradient_type<typename InputVector::value_type>>
+ &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 <typename Number>
+ 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 <typename InputVector>
+ 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 <typename ValueType>
+ std::vector<ValueType> &
+ outer_to_inner_values(std::vector<ValueType> &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 <typename ValueType>
+ void
+ inner_to_outer_values(const std::vector<ValueType> &inner_values,
+ std::vector<ValueType> &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<unsigned int> &dof_renumbering,
+ const std::vector<unsigned int> &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 <typename ViewType>
+ RenumberedView<ViewType>::RenumberedView(const ViewType &view,
+ const RenumberingData &data)
+ : view(view)
+ , data(data)
+ {}
+
+
+
+ template <typename ViewType>
+ template <typename Number>
+ inline std::string
+ RenumberedView<ViewType>::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 ViewType>
+ typename RenumberedView<ViewType>::value_type
+ RenumberedView<ViewType>::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 ViewType>
+ typename RenumberedView<ViewType>::gradient_type
+ RenumberedView<ViewType>::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 <typename ViewType>
+ template <typename ValueType>
+ std::vector<ValueType> &
+ RenumberedView<ViewType>::outer_to_inner_values(
+ std::vector<ValueType> &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<std::vector<ValueType>>(
+ name, data.n_inner_quadrature_points);
+ return inner_values;
+ }
+ }
+
+
+
+ template <typename ViewType>
+ template <typename VectorType>
+ const VectorType &
+ RenumberedView<ViewType>::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<VectorType>(
+ 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 <typename ViewType>
+ template <typename ValueType>
+ void
+ RenumberedView<ViewType>::inner_to_outer_values(
+ const std::vector<ValueType> &inner_values,
+ std::vector<ValueType> &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 <typename ViewType>
+ template <typename Number>
+ void
+ RenumberedView<ViewType>::get_function_values(
+ const ReadVector<Number> &fe_function,
+ std::vector<solution_value_type<Number>> &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 <typename ViewType>
+ template <typename InputVector>
+ void
+ RenumberedView<ViewType>::get_function_values_from_local_dof_values(
+ const InputVector &dof_values,
+ std::vector<solution_value_type<typename InputVector::value_type>> &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 <typename ViewType>
+ template <typename Number>
+ void
+ RenumberedView<ViewType>::get_function_gradients(
+ const ReadVector<Number> &fe_function,
+ std::vector<solution_gradient_type<Number>> &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 <typename ViewType>
+ template <typename InputVector>
+ void
+ RenumberedView<ViewType>::get_function_gradients_from_local_dof_values(
+ const InputVector &dof_values,
+ std::vector<solution_gradient_type<typename InputVector::value_type>>
+ &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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/fe/fe_coupling_values.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/tria.h>
+
+#include <fstream>
+#include <iostream>
+
+#include "../tests.h"
+
+#include "../test_grids.h"
+
+template <int dim>
+void
+test()
+{
+ Triangulation<dim> tria;
+ TestGrids::hyper_line(tria, 2);
+
+ DoFHandler<dim> dofh(tria);
+ FE_Q<dim> fe(1);
+ dofh.distribute_dofs(fe);
+
+ MappingQ<dim> mapping(1);
+ UpdateFlags update_flags =
+ update_values | update_quadrature_points | update_JxW_values;
+
+ FEValues<dim> fv1(mapping, fe, QGauss<dim>(fe.degree + 1), update_flags);
+ fv1.reinit(dofh.begin_active());
+
+ FEValuesViews::Scalar<dim> fv1_scalar(fv1, 0);
+
+ using iota = std_cxx20::ranges::iota_view<unsigned int, unsigned int>;
+ {
+ // 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<unsigned int> renumber(id.begin(), id.end());
+ FEValuesViews::RenumberingData data(fv1.dofs_per_cell,
+ fv1.n_quadrature_points,
+ renumber);
+ FEValuesViews::RenumberedView<FEValuesViews::Scalar<dim>>
+ 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<unsigned int> 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<FEValuesViews::Scalar<dim>>
+ 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<unsigned int> renumber(id.begin(), id.end());
+ FEValuesViews::RenumberingData data(fv1.dofs_per_cell,
+ fv1.n_quadrature_points,
+ {},
+ renumber);
+ FEValuesViews::RenumberedView<FEValuesViews::Scalar<dim>>
+ 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<unsigned int> 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<FEValuesViews::Scalar<dim>>
+ 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>();
+}
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/function_parser.h>
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/fe/fe_coupling_values.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include <fstream>
+#include <iostream>
+
+#include "../tests.h"
+
+#include "../test_grids.h"
+
+template <int dim>
+void
+test()
+{
+ Triangulation<dim> tria;
+ TestGrids::hyper_line(tria, 2);
+
+ DoFHandler<dim> dofh(tria);
+ FE_Q<dim> fe(1);
+ dofh.distribute_dofs(fe);
+
+ // Create a non trivial function of x and y
+ FunctionParser<dim> f("x+10*y");
+
+ Vector<double> vec(dofh.n_dofs());
+ VectorTools::interpolate(dofh, f, vec);
+
+ MappingQ<dim> mapping(1);
+ UpdateFlags update_flags =
+ update_values | update_quadrature_points | update_JxW_values;
+
+ FEValues<dim> fe_v(mapping, fe, QGauss<dim>(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<unsigned int, unsigned int>;
+
+ std::vector<unsigned int> dof_renumbering(fe_v.dofs_per_cell + 1);
+ std::vector<unsigned int> 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<double> function_values(fe_v.n_quadrature_points);
+ std::vector<double> 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<FEValuesViews::Scalar<dim>>
+ 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<double> dof_values(fe_v.dofs_per_cell);
+ std::vector<double> 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>();
+}
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/function_parser.h>
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/fe/fe_coupling_values.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include <fstream>
+#include <iostream>
+
+#include "../tests.h"
+
+#include "../test_grids.h"
+
+template <int dim>
+void
+test()
+{
+ Triangulation<dim> tria;
+ TestGrids::hyper_line(tria, 2);
+
+ DoFHandler<dim> dofh(tria);
+ FE_Q<dim> fe(2);
+ dofh.distribute_dofs(fe);
+
+ // Create a non trivial function of x and y
+ FunctionParser<dim> f("(x^2)/2+(10*y^2)/2");
+
+ Vector<double> vec(dofh.n_dofs());
+ VectorTools::interpolate(dofh, f, vec);
+
+ MappingQ<dim> mapping(1);
+ UpdateFlags update_flags = update_gradients;
+
+ FEValues<dim> fe_v(mapping, fe, QGauss<dim>(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<unsigned int, unsigned int>;
+
+ std::vector<unsigned int> dof_renumbering(fe_v.dofs_per_cell + 1);
+ std::vector<unsigned int> 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<Tensor<1, dim>> function_gradients(fe_v.n_quadrature_points);
+ std::vector<Tensor<1, dim>> 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<FEValuesViews::Scalar<dim>>
+ 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<double> dof_values(fe_v.dofs_per_cell);
+ std::vector<double> 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>();
+}
projects / dealii.git / commitdiff