From 241153c716bd07074566a3821ee807dd51e07520 Mon Sep 17 00:00:00 2001 From: Simon Sticko Date: Mon, 5 Apr 2021 13:34:10 +0200 Subject: [PATCH] Add classes restricting a function to a lower dimension Add a class CoordinateRestriction that corresponds to restricting one coordinate of a dealii::Function, f, to a given value. That is, if Z is a fixed value we create the new function g = g(x, y) = f(x, y, Z). Correspondingly, add a class PointRestriction that takes a point (Y, Z) and restricts a dealii::Function to the point, i.e, we create the new function g = g(x) = f(x, Y, Z). --- include/deal.II/base/function_restriction.h | 174 +++++++++++++++++ source/base/CMakeLists.txt | 2 + source/base/function_restriction.cc | 204 ++++++++++++++++++++ source/base/function_restriction.inst.in | 11 ++ tests/base/coordinate_restriction.cc | 115 +++++++++++ tests/base/coordinate_restriction.output | 28 +++ tests/base/create_higher_dim_point.cc | 113 +++++++++++ tests/base/create_higher_dim_point.output | 10 + tests/base/point_restriction.cc | 124 ++++++++++++ tests/base/point_restriction.output | 28 +++ 10 files changed, 809 insertions(+) create mode 100644 include/deal.II/base/function_restriction.h create mode 100644 source/base/function_restriction.cc create mode 100644 source/base/function_restriction.inst.in create mode 100644 tests/base/coordinate_restriction.cc create mode 100644 tests/base/coordinate_restriction.output create mode 100644 tests/base/create_higher_dim_point.cc create mode 100644 tests/base/create_higher_dim_point.output create mode 100644 tests/base/point_restriction.cc create mode 100644 tests/base/point_restriction.output diff --git a/include/deal.II/base/function_restriction.h b/include/deal.II/base/function_restriction.h new file mode 100644 index 0000000000..e96c4a0f07 --- /dev/null +++ b/include/deal.II/base/function_restriction.h @@ -0,0 +1,174 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2010 - 2019 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_function_restriction_h +#define dealii_function_restriction_h + +#include +#include +#include +#include + +DEAL_II_NAMESPACE_OPEN + +namespace Functions +{ + /** + * This class takes a function in `dim + 1` dimensions and creates a new + * function in one dimension lower by restricting one of the coordinates to a + * given value. Mathematically this corresponds to taking a function + * $f = f(x, y, z)$, + * a fixed value, $Z$, and defining a new function (the restriction) + * $g = g(x, y) = f(x, y, Z)$. + * Using this class, this translates to + * @code + * Function<3> & function = ... + * double z = ... + * unsigned int restricted_direction = 2; + * CoordinateRestriction<2> restriction(function, restricted_direction, z); + * @endcode + * + * The `dim`-dimensional coordinates on the restriction are ordered starting + * from the restricted (`dim + 1`)-coordinate. In particular, this means that + * if the $y$-coordinate is locked to $Y$ in 3D, the coordinates are ordered + * as $(z, x)$ on the restriction: + * $g = g(z, x) = f(x, Y, z)$. + * This is the same convention as in BoundingBox::cross_section. + */ + template + class CoordinateRestriction : public Function + { + public: + /** + * Constructor, takes the (`dim + 1`)-coordinate direction and the value + * that the incoming function should be restricted to. + * + * A pointer to the incoming function is stored internally, so the function + * must have a longer lifetime than the created restriction. + */ + CoordinateRestriction(const Function &function, + const unsigned int direction, + const double coordinate_value); + + double + value(const Point &point, const unsigned int component) const override; + + Tensor<1, dim> + gradient(const Point & point, + const unsigned int component) const override; + + SymmetricTensor<2, dim> + hessian(const Point & point, + const unsigned int component) const override; + + private: + // The higher-dimensional function that has been restricted. + const SmartPointer> function; + + // The (`dim + 1`)-coordinate direction that has been restricted. + const unsigned int restricted_direction; + + // Value of the restriced coordinate. + const double coordinate_value; + }; + + + + /** + * This class creates a 1-dimensional function from a `dim + 1` dimensional + * function by restricting `dim` of the coordinate values to a given point. + * Mathematically this corresponds to taking a function, $f = f(x, y, z)$, and + * a point $(Y, Z)$, and defining a new function $g = g(x) = f(x, Y, Z)$. + * Using this class, this translates to + * @code + * Function<3> & function = ... + * Point<2> point(y, z); + * unsigned int open_direction = 0; + * PointRestriction<2> restriction(function, open_direction, point); + * @endcode + * + * The coordinates of the point will be expanded in the higher-dimensional + * functions coordinates starting from the open-direction (and wrapping + * around). In particular, if we restrict to a point $(Z, X)$ and choose to + * keep the y-direction open, the restriction that is created is the function + * $g(y) = f(X, y, Z)$. + * This is consistent with the convention in BoundingBox::cross_section. + */ + template + class PointRestriction : public Function<1> + { + public: + /** + * Constructor, takes the point that the incoming function should be + * restricted to and which (`dim + 1`)-dimensional coordinate direction + * should be kept "open". + * + * A pointer to the incoming function is stored internally, so the function + * must have a longer lifetime than the created restriction. + */ + PointRestriction(const Function &function, + const unsigned int open_direction, + const Point & point); + + double + value(const Point<1> &point, const unsigned int component) const override; + + Tensor<1, 1> + gradient(const Point<1> & point, + const unsigned int component) const override; + + SymmetricTensor<2, 1> + hessian(const Point<1> &point, const unsigned int component) const override; + + private: + // The higher-dimensional function that has been restricted. + const SmartPointer> function; + + // The (`dim + 1`)-coordinate direction that is kept "open" + const unsigned int open_direction; + + // The point that we have restricted the above function to. + const Point point; + }; + +} // namespace Functions + + +namespace internal +{ + /** + * Creates a (`dim + 1`)-dimensional point by copying over the coordinates of + * the incoming `dim`-dimensional point and setting the "missing" + * (`dim + 1`)-dimensional component to the incoming coordinate value. + * + * For example, given the input + * $\{(x, y), 2, z \}$ this function creates the point $(x, y, z)$. + * + * The coordinates of the `dim`-dimensional point are written to the + * coordinates of the (`dim + 1`)-dimensional point in the order of the + * convention given by the function coordinate_to_one_dim_higher. Thus, the + * order of coordinates on the lower-dimensional point are not preserved: + * $\{(z, x), 1, y \}$ creates the point $(x, y, z)$. + */ + template + Point + create_higher_dim_point(const Point & point, + const unsigned int component_in_dim_plus_1, + const double coordinate_value); +} // namespace internal + +DEAL_II_NAMESPACE_CLOSE + +#endif /* dealii_function_restriction_h */ diff --git a/source/base/CMakeLists.txt b/source/base/CMakeLists.txt index 024d2a9b02..15551efd42 100644 --- a/source/base/CMakeLists.txt +++ b/source/base/CMakeLists.txt @@ -36,6 +36,7 @@ SET(_unity_include_src function_lib.cc function_lib_cutoff.cc function_parser.cc + function_restriction.cc function_spherical.cc function_time.cc function_tools.cc @@ -131,6 +132,7 @@ SET(_inst data_out_base.inst.in function.inst.in function_level_set.inst.in + function_restriction.inst.in function_time.inst.in function_tools.inst.in incremental_function.inst.in diff --git a/source/base/function_restriction.cc b/source/base/function_restriction.cc new file mode 100644 index 0000000000..71d9bf6590 --- /dev/null +++ b/source/base/function_restriction.cc @@ -0,0 +1,204 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2010 - 2019 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. +// +// --------------------------------------------------------------------- + +#include +#include + +DEAL_II_NAMESPACE_OPEN +namespace internal +{ + template + Point + create_higher_dim_point(const Point & point, + const unsigned int component_in_dim_plus_1, + const double coordinate_value) + { + AssertIndexRange(component_in_dim_plus_1, dim + 1); + + Point output; + output(component_in_dim_plus_1) = coordinate_value; + for (unsigned int d = 0; d < dim; ++d) + { + const unsigned int component_to_write_to = + dealii::internal::coordinate_to_one_dim_higher( + component_in_dim_plus_1, d); + output(component_to_write_to) = point(d); + } + + return output; + } +} // namespace internal + + + +namespace Functions +{ + template + CoordinateRestriction::CoordinateRestriction( + const Function &function, + const unsigned int direction, + const double coordinate_value) + : function(&function) + , restricted_direction(direction) + , coordinate_value(coordinate_value) + { + AssertIndexRange(restricted_direction, dim + 1); + } + + + + template + double + CoordinateRestriction::value(const Point & point, + const unsigned int component) const + { + const Point full_point = + internal::create_higher_dim_point(point, + restricted_direction, + coordinate_value); + + return function->value(full_point, component); + } + + + + template + Tensor<1, dim> + CoordinateRestriction::gradient(const Point & point, + const unsigned int component) const + { + const Point full_point = + internal::create_higher_dim_point(point, + restricted_direction, + coordinate_value); + + const Tensor<1, dim + 1> full_gradient = + function->gradient(full_point, component); + + // CoordinateRestriction is constant in restricted direction. Through away + // the derivatives with respect to this direction and copy the other + // values. + Tensor<1, dim> grad; + for (unsigned int d = 0; d < dim; ++d) + { + const unsigned int index_to_write_from = + internal::coordinate_to_one_dim_higher(restricted_direction, d); + grad[d] = full_gradient[index_to_write_from]; + } + return grad; + } + + + + template + SymmetricTensor<2, dim> + CoordinateRestriction::hessian(const Point & point, + const unsigned int component) const + { + const Point full_point = + internal::create_higher_dim_point(point, + restricted_direction, + coordinate_value); + + const Tensor<2, dim + 1> full_hessian = + function->hessian(full_point, component); + + // CoordinateRestriction is constant in restricted direction. Through away + // the derivatives with respect to this direction and copy the other + // values. + SymmetricTensor<2, dim> hess; + for (unsigned int i = 0; i < dim; ++i) + { + const unsigned int i_to_write_from = + internal::coordinate_to_one_dim_higher(restricted_direction, i); + for (unsigned int j = 0; j < dim; ++j) + { + const unsigned int j_to_write_from = + internal::coordinate_to_one_dim_higher(restricted_direction, + j); + hess[i][j] = full_hessian[i_to_write_from][j_to_write_from]; + } + } + return hess; + } + + + + template + PointRestriction::PointRestriction(const Function &function, + const unsigned int open_direction, + const Point & point) + : function(&function) + , open_direction(open_direction) + , point(point) + { + AssertIndexRange(open_direction, dim + 1); + } + + + + template + double + PointRestriction::value(const Point<1> & point_1D, + const unsigned int component) const + { + const Point full_point = + internal::create_higher_dim_point(point, open_direction, point_1D(0)); + return function->value(full_point, component); + } + + + + template + Tensor<1, 1> + PointRestriction::gradient(const Point<1> & point_1D, + const unsigned int component) const + { + const Point full_point = + internal::create_higher_dim_point(point, open_direction, point_1D(0)); + const Tensor<1, dim + 1> full_gradient = + function->gradient(full_point, component); + + // The PointRestrictions is constant in all but the open direction. Throw + // away the derivatives in all but this direction. + Tensor<1, 1> grad; + grad[0] = full_gradient[open_direction]; + return grad; + } + + + + template + SymmetricTensor<2, 1> + PointRestriction::hessian(const Point<1> & point_1D, + const unsigned int component) const + { + const Point full_point = + internal::create_higher_dim_point(point, open_direction, point_1D(0)); + const Tensor<2, dim + 1> full_hessian = + function->hessian(full_point, component); + + // The PointRestrictions is constant in all but the open direction. Throw + // away the derivatives in all but this direction. + SymmetricTensor<2, 1> hess; + hess[0][0] = full_hessian[open_direction][open_direction]; + return hess; + } + + + +} // namespace Functions +#include "function_restriction.inst" +DEAL_II_NAMESPACE_CLOSE diff --git a/source/base/function_restriction.inst.in b/source/base/function_restriction.inst.in new file mode 100644 index 0000000000..7ce0681563 --- /dev/null +++ b/source/base/function_restriction.inst.in @@ -0,0 +1,11 @@ +for (deal_II_dimension : DIMENSIONS) + { + template class Functions::PointRestriction; + + template Point internal::create_higher_dim_point( + const Point &, const unsigned int, const double); + +#if deal_II_dimension != 3 + template class Functions::CoordinateRestriction; +#endif + } diff --git a/tests/base/coordinate_restriction.cc b/tests/base/coordinate_restriction.cc new file mode 100644 index 0000000000..2bc454cee9 --- /dev/null +++ b/tests/base/coordinate_restriction.cc @@ -0,0 +1,115 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2010 - 2019 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 that CoordinateRestriction returns the derivatives of its underlying + * function truncated correctly. + */ + + +#include + +#include "../tests.h" + +namespace +{ + using namespace dealii; + + /* + * A simple test-function which is constant in space, but varies + * in each component. + */ + template + class TestFunction : public Function + { + public: + double + value(const Point &, const unsigned int) const override + { + return 1; + } + + Tensor<1, dim> + gradient(const Point &, const unsigned int) const override + { + Tensor<1, dim> tensor; + for (unsigned int i = 0; i < dim; ++i) + tensor[i] = i + 1; + return tensor; + } + + SymmetricTensor<2, dim> + hessian(const Point &, const unsigned int) const override + { + SymmetricTensor<2, dim> tensor; + for (unsigned int i = 0; i < dim; ++i) + for (unsigned int j = 0; j < dim; ++j) + tensor[i][j] = 1 + i + j; + + return tensor; + } + }; + + + /* + * Create a CoordinateRestriction that has the incoming direction of the + * above TestFunction restricted. Write value, gradient & hessian + * to deallog to check that they are truncated as expected. + */ + template + void + test_truncates_derivatives_correctly(const unsigned int restricted_direction) + { + deallog << "Restricted direction = " << restricted_direction << std::endl; + + const TestFunction function; + const double coordinate_value = 0; + + const Functions::CoordinateRestriction restriction( + function, restricted_direction, coordinate_value); + + const Point point; + const unsigned int component = 0; + + deallog << "value = " << restriction.value(point, component) << std::endl; + deallog << "gradient = " << restriction.gradient(point, component) + << std::endl; + deallog << "Hessian = " << restriction.hessian(point, component) + << std::endl; + + deallog << std::endl; + } + + + + template + void + test_all_coordinate_directions() + { + deallog << "dim = " << dim << std::endl; + for (unsigned int direction = 0; direction < dim; direction++) + test_truncates_derivatives_correctly(direction); + } + +} // namespace + +int +main() +{ + initlog(); + test_all_coordinate_directions<2>(); + test_all_coordinate_directions<3>(); +} diff --git a/tests/base/coordinate_restriction.output b/tests/base/coordinate_restriction.output new file mode 100644 index 0000000000..98470a126e --- /dev/null +++ b/tests/base/coordinate_restriction.output @@ -0,0 +1,28 @@ + +DEAL::dim = 2 +DEAL::Restricted direction = 0 +DEAL::value = 1.00000 +DEAL::gradient = 2.00000 +DEAL::Hessian = 3.00000 +DEAL:: +DEAL::Restricted direction = 1 +DEAL::value = 1.00000 +DEAL::gradient = 1.00000 +DEAL::Hessian = 1.00000 +DEAL:: +DEAL::dim = 3 +DEAL::Restricted direction = 0 +DEAL::value = 1.00000 +DEAL::gradient = 2.00000 3.00000 +DEAL::Hessian = 3.00000 4.00000 4.00000 5.00000 +DEAL:: +DEAL::Restricted direction = 1 +DEAL::value = 1.00000 +DEAL::gradient = 3.00000 1.00000 +DEAL::Hessian = 5.00000 3.00000 3.00000 1.00000 +DEAL:: +DEAL::Restricted direction = 2 +DEAL::value = 1.00000 +DEAL::gradient = 1.00000 2.00000 +DEAL::Hessian = 1.00000 2.00000 2.00000 3.00000 +DEAL:: diff --git a/tests/base/create_higher_dim_point.cc b/tests/base/create_higher_dim_point.cc new file mode 100644 index 0000000000..08461914d4 --- /dev/null +++ b/tests/base/create_higher_dim_point.cc @@ -0,0 +1,113 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2010 - 2019 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. +// +// --------------------------------------------------------------------- + + +#include + +#include "../tests.h" + +namespace +{ + using namespace dealii; + + // Test the function create_higher_dim_point by making sure that we + // can construct the point (x, y, z), + const double x = 1, y = 2, z = 3; + // for each possible component. That is, we test that we get the following + // results: + // {point, component, coordinate} -> result + // {(y, z), 0, x} -> (x, y, z) + // {(z, x), 1, y} -> (x, y, z) + // {(y, z), 2, x} -> (x, y, z) + // + // The same thing is done for the 2D case. + void + test3D() + { + const int dim = 3; + deallog << "dim=" << dim << std::endl; + { + const unsigned int component = 0; + const Point low_dim_point(y, z); + const Point point = + internal::create_higher_dim_point(low_dim_point, component, x); + deallog << "Point: " << point << std::endl; + } + { + const unsigned int component = 1; + // The order (z,x) here follows convention in + // coordinate_to_one_dim_higher. + const Point low_dim_point(z, x); + const Point point = + internal::create_higher_dim_point(low_dim_point, component, y); + deallog << "Point: " << point << std::endl; + } + { + const unsigned int component = 2; + Point low_dim_point(x, y); + Point point = + internal::create_higher_dim_point(low_dim_point, component, z); + deallog << "Point: " << point << std::endl; + } + } + + + + void + test2D() + { + const int dim = 2; + deallog << "dim=" << dim << std::endl; + { + const unsigned int component = 0; + const Point low_dim_point(y); + const Point point = + internal::create_higher_dim_point(low_dim_point, component, x); + deallog << "Point: " << point << std::endl; + } + { + const unsigned int component = 1; + const Point low_dim_point(x); + const Point point = + internal::create_higher_dim_point(low_dim_point, component, y); + deallog << "Point: " << point << std::endl; + } + } + + + + void + test1D() + { + const int dim = 1; + deallog << "dim=" << dim << std::endl; + const unsigned int component = 0; + const Point low_dim_point; + const Point point = + internal::create_higher_dim_point(low_dim_point, component, x); + deallog << "Point: " << point << std::endl; + } +} // namespace + + + +int +main() +{ + initlog(); + test1D(); + test2D(); + test3D(); +} diff --git a/tests/base/create_higher_dim_point.output b/tests/base/create_higher_dim_point.output new file mode 100644 index 0000000000..4f2d84f659 --- /dev/null +++ b/tests/base/create_higher_dim_point.output @@ -0,0 +1,10 @@ + +DEAL::dim=1 +DEAL::Point: 1.00000 +DEAL::dim=2 +DEAL::Point: 1.00000 2.00000 +DEAL::Point: 1.00000 2.00000 +DEAL::dim=3 +DEAL::Point: 1.00000 2.00000 3.00000 +DEAL::Point: 1.00000 2.00000 3.00000 +DEAL::Point: 1.00000 2.00000 3.00000 diff --git a/tests/base/point_restriction.cc b/tests/base/point_restriction.cc new file mode 100644 index 0000000000..abe76ee75c --- /dev/null +++ b/tests/base/point_restriction.cc @@ -0,0 +1,124 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2010 - 2019 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 that PointRestriction returns the derivatives of its underlying function + * truncated correctly. + */ + + +#include + +#include "../tests.h" + +namespace +{ + using namespace dealii; + + /* + * A simple test-function which is constant in space, but varies + * in each gradient component, and on the diagonal of the Hessian. + */ + template + class TestFunction : public Function + { + public: + double + value(const Point &, const unsigned int) const override + { + return 1; + } + + Tensor<1, dim> + gradient(const Point &, const unsigned int) const override + { + Tensor<1, dim> tensor; + for (unsigned int i = 0; i < dim; ++i) + tensor[i] = i + 1; + return tensor; + } + + SymmetricTensor<2, dim> + hessian(const Point &, const unsigned int) const override + { + SymmetricTensor<2, dim> tensor; + for (unsigned int i = 0; i < dim; ++i) + tensor[i][i] = i + 1; + + return tensor; + } + }; + + + /* + * Create a PointRestriction that has the incoming direction open and has + * the above TestFunction as underlying function. Check that the values that + * + * PointRestriction::values + * PointRestriction::gradient + * PointRestriction::hessian + * + * return are + * + * value + * gradient[open_direction] + * hessian[open_direction][open_direction] + * + * where value, gradient and hessian are the values from TestFunction. + */ + template + void + test_truncates_derivatives_correctly(const unsigned int open_direction) + { + deallog << "open direction = " << open_direction << std::endl; + const Point locked_point; + const TestFunction function; + + const Functions::PointRestriction restriction(function, + open_direction, + locked_point); + + const Point<1> point; + const unsigned int component = 0; + + deallog << restriction.value(point, component) << std::endl; + deallog << restriction.gradient(point, component) << std::endl; + deallog << restriction.hessian(point, component) << std::endl; + + deallog << std::endl; + } + + + + template + void + test_all_coordinate_directions() + { + deallog << "dim = " << dim << std::endl; + + for (unsigned int direction = 0; direction < dim; direction++) + test_truncates_derivatives_correctly(direction); + } + +} // namespace + +int +main() +{ + initlog(); + test_all_coordinate_directions<2>(); + test_all_coordinate_directions<3>(); +} diff --git a/tests/base/point_restriction.output b/tests/base/point_restriction.output new file mode 100644 index 0000000000..4c0411b666 --- /dev/null +++ b/tests/base/point_restriction.output @@ -0,0 +1,28 @@ + +DEAL::dim = 2 +DEAL::open direction = 0 +DEAL::1.00000 +DEAL::1.00000 +DEAL::1.00000 +DEAL:: +DEAL::open direction = 1 +DEAL::1.00000 +DEAL::2.00000 +DEAL::2.00000 +DEAL:: +DEAL::dim = 3 +DEAL::open direction = 0 +DEAL::1.00000 +DEAL::1.00000 +DEAL::1.00000 +DEAL:: +DEAL::open direction = 1 +DEAL::1.00000 +DEAL::2.00000 +DEAL::2.00000 +DEAL:: +DEAL::open direction = 2 +DEAL::1.00000 +DEAL::3.00000 +DEAL::3.00000 +DEAL:: -- 2.39.5