From 1fa30e55f82d02b08469a8bc2305c0e75eb2012e Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Wed, 18 Jun 2025 18:06:18 -0600 Subject: [PATCH] Add a function that converts a lambda function into a vector Function object. --- include/deal.II/base/function.h | 118 +++++++++++++++++++++- include/deal.II/base/function.templates.h | 104 +++++++++++++++++++ source/base/function.inst.in | 2 + 3 files changed, 221 insertions(+), 3 deletions(-) diff --git a/include/deal.II/base/function.h b/include/deal.II/base/function.h index 2f555f2ff3..acb9cc0072 100644 --- a/include/deal.II/base/function.h +++ b/include/deal.II/base/function.h @@ -1126,12 +1126,12 @@ public: * By default, create a Vector object of the same size as * tensor_function returns, i.e., with dim components. * - * @param tensor_function The TensorFunction that will form one component of + * @param tensor_function The TensorFunction that will form `dim` components of * the resulting Vector Function object. * @param n_components The total number of vector components of the - * resulting TensorFunction object. + * resulting TensorFunction object. This clearly has to be at least `dim`. * @param selected_component The first component that should be filled by - * the first argument. This should be such that the entire tensor_function + * the first argument. This should be such that the entire `tensor_function` * fits inside the n_component length return vector. */ explicit VectorFunctionFromTensorFunction( @@ -1243,6 +1243,118 @@ private: }; +/** + * This class is built as a means of translating the Tensor<1,dim, + * RangeNumberType> values produced by function objects that + * for a given point return a tensor, + * and returning them as a multiple component version of the same thing as a + * Vector for use in, for example, the VectorTools::interpolate or the many + * other functions taking Function objects. It allows the user to place the + * desired components into an n_components long vector starting at + * the selected_component location in that vector and have all other + * components be 0. + * + * For example: Say you created a function object that returns the gravity + * (here, a radially inward pointing vector of magnitude 9.81): + * @code + * const auto gravity + * = [](const Point &p) -> Tensor<1,dim> { return -9.81 * (p / + * p.norm()); } + * @endcode + * and you want to interpolate this onto your mesh using the + * VectorTools::interpolate function, with a finite element for the + * DoFHandler object has 3 copies of a finite element with dim + * components, for a total of 3*dim components. To interpolate onto that + * DoFHandler, you need an object of type Function that has 3*dim vector + * components. Creating such an object from the existing gravity + * object is done using this piece of code: + * @code + * VectorFunctionFromTensorFunctionObject + * gravity_vector_function(gravity, 0, 3*dim); + * @endcode + * where the dim components of the `gravity` function are placed + * into the first dim components of the function object. + * + * @ingroup functions + */ +template +class VectorFunctionFromTensorFunctionObject + : public Function +{ +public: + /** + * Given a function object that takes a Point and returns a + * Tensor<1,dim, RangeNumberType> value, convert this into an object + * that matches the Function@ interface. + * + * By default, create a Vector object of the same size as + * tensor_function returns, i.e., with dim components. + * + * @param tensor_function_object The TensorFunction that will form `dim` components of + * the resulting Vector Function object. + * @param n_components The total number of vector components of the + * resulting TensorFunction object. This clearly has to be at least `dim`. + * @param selected_component The first component that should be filled by + * the first argument. This should be such that the entire tensor_function + * fits inside the n_component length return vector. + */ + explicit VectorFunctionFromTensorFunctionObject( + const std::function(const Point &)> + &tensor_function_object, + const unsigned int selected_component = 0, + const unsigned int n_components = dim); + + /** + * This destructor is defined as virtual so as to coincide with all other + * aspects of class. + */ + virtual ~VectorFunctionFromTensorFunctionObject() override = default; + + /** + * Return a single component of a vector-valued function at a given point. + */ + virtual RangeNumberType + value(const Point &p, const unsigned int component = 0) const override; + + /** + * Return all components of a vector-valued function at a given point. + * + * values shall have the right size beforehand, i.e. #n_components. + */ + virtual void + vector_value(const Point &p, + Vector &values) const override; + + /** + * Return all components of a vector-valued function at a list of points. + * + * value_list shall be the same size as points and each + * element of the vector will be passed to vector_value() to evaluate the + * function + */ + virtual void + vector_value_list( + const std::vector> &points, + std::vector> &value_list) const override; + +private: + /** + * The TensorFunction object which we call when this class's vector_value() + * or vector_value_list() functions are called. + */ + const std::function(const Point &)> + tensor_function_object; + + /** + * The first vector component whose value is to be filled by the given + * TensorFunction. The values will be placed in components + * selected_component to selected_component+dim-1 for a + * TensorFunction<1,dim, RangeNumberType> object. + */ + const unsigned int selected_component; +}; + + #ifndef DOXYGEN // icc 2018 complains about an undefined reference // if we put this in the templates.h file diff --git a/include/deal.II/base/function.templates.h b/include/deal.II/base/function.templates.h index d263832b0b..371946e42f 100644 --- a/include/deal.II/base/function.templates.h +++ b/include/deal.II/base/function.templates.h @@ -867,6 +867,8 @@ VectorFunctionFromTensorFunction::vector_value_list( points[p], value_list[p]); } + + template inline Tensor<1, dim, RangeNumberType> VectorFunctionFromTensorFunction::gradient( @@ -891,6 +893,8 @@ VectorFunctionFromTensorFunction::gradient( return tensor_gradient[component - selected_component]; } + + template void VectorFunctionFromTensorFunction::vector_gradient( @@ -904,6 +908,8 @@ VectorFunctionFromTensorFunction::vector_gradient( gradients[i] = gradient(p, i); } + + template void VectorFunctionFromTensorFunction::gradient_list( @@ -939,6 +945,8 @@ VectorFunctionFromTensorFunction::vector_gradient_list( } } + + template void VectorFunctionFromTensorFunction::vector_gradients( @@ -952,6 +960,7 @@ VectorFunctionFromTensorFunction::vector_gradients( } + template FunctionFromFunctionObjects::FunctionFromFunctionObjects( const unsigned int n_components, @@ -1058,6 +1067,101 @@ FunctionFromFunctionObjects::set_function_gradients( }; } + + +template +VectorFunctionFromTensorFunctionObject:: + VectorFunctionFromTensorFunctionObject( + const std::function(const Point &)> + &tensor_function_object, + const unsigned int selected_component, + const unsigned int n_components) + : Function(n_components) + , tensor_function_object(tensor_function_object) + , selected_component(selected_component) +{ + // Verify that the Tensor<1,dim,RangeNumberType> will fit in the given length + // selected_components and not hang over the end of the vector. + AssertIndexRange(selected_component + dim - 1, this->n_components); +} + + + +template +inline RangeNumberType +VectorFunctionFromTensorFunctionObject::value( + const Point &p, + const unsigned int component) const +{ + AssertIndexRange(component, this->n_components); + + // if the requested component is out of the range selected, then we can + // return early + if ((component < selected_component) || + (component >= selected_component + dim)) + return 0; + + // otherwise retrieve the values from the tensor_function to be + // placed at the selected_component to + // selected_component + dim - 1 elements of the Vector + // values and pick the correct one + const Tensor<1, dim, RangeNumberType> tensor_value = + tensor_function_object(p); + + return tensor_value[component - selected_component]; +} + + +template +inline void +VectorFunctionFromTensorFunctionObject::vector_value( + const Point &p, + Vector &values) const +{ + Assert(values.size() == this->n_components, + ExcDimensionMismatch(values.size(), this->n_components)); + + // Retrieve the values from the tensor_function to be placed at + // the selected_component to + // selected_component + dim - 1 elements of the Vector + // values. + const Tensor<1, dim, RangeNumberType> tensor_value = + tensor_function_object(p); + + // First we make all elements of values = 0 + values = 0; + + // Second we adjust the desired components to take on the values in + // tensor_value. + for (unsigned int i = 0; i < dim; ++i) + values(i + selected_component) = tensor_value[i]; +} + + +/** + * Member function vector_value_list is the interface for giving a + * list of points (vector >) of which to evaluate + * using the vector_value member function. Again, this function is + * written so as to not replicate the function definition but passes each + * point on to vector_value to be evaluated. + */ +template +void +VectorFunctionFromTensorFunctionObject::vector_value_list( + const std::vector> &points, + std::vector> &value_list) const +{ + Assert(value_list.size() == points.size(), + ExcDimensionMismatch(value_list.size(), points.size())); + + const unsigned int n_points = points.size(); + + for (unsigned int p = 0; p < n_points; ++p) + VectorFunctionFromTensorFunctionObject::vector_value( + points[p], value_list[p]); +} + + DEAL_II_NAMESPACE_CLOSE #endif /* dealii_function_templates_h */ diff --git a/source/base/function.inst.in b/source/base/function.inst.in index 483d1d3ff5..442476435f 100644 --- a/source/base/function.inst.in +++ b/source/base/function.inst.in @@ -24,6 +24,7 @@ for (S : REAL_SCALARS; dim : SPACE_DIMENSIONS) template class VectorFunctionFromScalarFunctionObject; template class FunctionFromFunctionObjects; template class VectorFunctionFromTensorFunction; + template class VectorFunctionFromTensorFunctionObject; } for (S : COMPLEX_SCALARS; dim : SPACE_DIMENSIONS) @@ -37,4 +38,5 @@ for (S : COMPLEX_SCALARS; dim : SPACE_DIMENSIONS) template class VectorFunctionFromScalarFunctionObject; template class FunctionFromFunctionObjects; template class VectorFunctionFromTensorFunction; + template class VectorFunctionFromTensorFunctionObject; } -- 2.39.5