From f423646c9da6079f99b270d74c44f7db6febe8f7 Mon Sep 17 00:00:00 2001 From: Simon Sticko Date: Tue, 20 Apr 2021 08:28:23 +0200 Subject: [PATCH] Add function estimating Function bounds over a BoundingBox Add a function taylor_estimate_function_bounds, which estimates the range of a Function over a BoundingBox by approximating the function by a 2nd order Taylor polynomial. Both the range of the function value and the range of each gradient component is estimated. --- doc/news/changes/minor/20210425SimonSticko | 6 + include/deal.II/base/function_tools.h | 71 +++++++++ source/base/CMakeLists.txt | 2 + source/base/function_tools.cc | 73 +++++++++ source/base/function_tools.inst.in | 25 +++ tests/base/taylor_estimate_function_bounds.cc | 144 ++++++++++++++++++ .../taylor_estimate_function_bounds.output | 16 ++ 7 files changed, 337 insertions(+) create mode 100644 doc/news/changes/minor/20210425SimonSticko create mode 100644 include/deal.II/base/function_tools.h create mode 100644 source/base/function_tools.cc create mode 100644 source/base/function_tools.inst.in create mode 100644 tests/base/taylor_estimate_function_bounds.cc create mode 100644 tests/base/taylor_estimate_function_bounds.output diff --git a/doc/news/changes/minor/20210425SimonSticko b/doc/news/changes/minor/20210425SimonSticko new file mode 100644 index 0000000000..cafb3b23f7 --- /dev/null +++ b/doc/news/changes/minor/20210425SimonSticko @@ -0,0 +1,6 @@ +New: Add a function taylor_estimate_function_bounds, which estimates the range +of the value and gradient components of a Function over a +BoundingBox, by approximating the function by a 2nd order Taylor +polynomial. +
+(Simon Sticko, 2021/04/25) diff --git a/include/deal.II/base/function_tools.h b/include/deal.II/base/function_tools.h new file mode 100644 index 0000000000..9f96cb1df4 --- /dev/null +++ b/include/deal.II/base/function_tools.h @@ -0,0 +1,71 @@ +// --------------------------------------------------------------------- +// +// 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_tools_h +#define dealii_function_tools_h + + +#include + +#include +#include + +DEAL_II_NAMESPACE_OPEN + +namespace FunctionTools +{ + /** + * Estimate bounds on the value and bounds on each gradient component of a + * Function, $f$, over a BoundingBox, by approximating it by a 2nd order + * Taylor polynomial starting from the box center. + * + * Each lower and upper bound is returned as a + * std::pair, such that the first entry is the + * lower bound, $L$, and the second is the upper bound, $U$, i.e. + * $f(x) \in [L, U]$. + * + * The function value, gradient, and Hessian are computed at the box center. + * The bounds on the value of the function are then estimated as + * + * $f(x) \in [f(x_c) - F, f(x_c) + F]$, + * where + * $F = \sum_i |\partial_i f(x_c)| h_i + * + 1/2 \sum_i \sum_j |\partial_i \partial_j f(x_c)| h_i h_j$. + * + * Here, $h_i$ is half the side length of the box in the $i$th coordinate + * direction, which is the distance we extrapolate. The bounds on the gradient + * components are estimated similarly as + * + * $\partial_i f \in [\partial_i f(x_c) - G_i, \partial_i f(x_c) + G_i]$, + * where + * $G_i = \sum_j |\partial_i \partial_j f(x_c)| h_j$. + * + * If the function has more than 1 component the @p component parameter can + * be used to specify which function component the bounds should be computed + * for. + */ + template + void + taylor_estimate_function_bounds( + const Function & function, + const BoundingBox & box, + std::pair & value_bounds, + std::array, dim> &gradient_bounds, + const unsigned int component = 0); + +} // namespace FunctionTools +DEAL_II_NAMESPACE_CLOSE + +#endif /* dealii_function_tools_h */ diff --git a/source/base/CMakeLists.txt b/source/base/CMakeLists.txt index 6cc778da4b..f8ce1e3469 100644 --- a/source/base/CMakeLists.txt +++ b/source/base/CMakeLists.txt @@ -37,6 +37,7 @@ SET(_unity_include_src function_parser.cc function_spherical.cc function_time.cc + function_tools.cc geometry_info.cc geometric_utilities.cc graph_coloring.cc @@ -129,6 +130,7 @@ SET(_inst data_out_base.inst.in function.inst.in function_time.inst.in + function_tools.inst.in incremental_function.inst.in geometric_utilities.inst.in mpi.inst.in diff --git a/source/base/function_tools.cc b/source/base/function_tools.cc new file mode 100644 index 0000000000..bc7b2711b9 --- /dev/null +++ b/source/base/function_tools.cc @@ -0,0 +1,73 @@ +// --------------------------------------------------------------------- +// +// 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 FunctionTools +{ + template + void + taylor_estimate_function_bounds( + const Function & function, + const BoundingBox & box, + std::pair & value_bounds, + std::array, dim> &gradient_bounds, + const unsigned int component) + { + const Point center = box.center(); + const double value = function.value(center, component); + const Tensor<1, dim> gradient = function.gradient(center, component); + const SymmetricTensor<2, dim> hessian = function.hessian(center, component); + + // Deviation from function value at the center, based on the + // Taylor-expansion: |f'| * dx + 1/2 * |f''| * dx^2, (in 1D). dx is half + // the side-length of the box. + double taylor_bound_f = 0; + + for (unsigned int i = 0; i < dim; ++i) + { + const double dx_i = .5 * box.side_length(i); + + taylor_bound_f += std::abs(gradient[i]) * dx_i; + + // Deviation from value of df/dx_i at the center, + // |f''| * dx, (in 1D). + double taylor_bound_dfdxi = 0; + + for (unsigned int j = 0; j < dim; ++j) + { + const double dx_j = .5 * box.side_length(j); + + taylor_bound_dfdxi += std::abs(hessian[i][j]) * dx_j; + taylor_bound_f += .5 * std::abs(hessian[i][j]) * dx_i * dx_j; + } + + gradient_bounds[i].first = gradient[i] - taylor_bound_dfdxi; + gradient_bounds[i].second = gradient[i] + taylor_bound_dfdxi; + } + + value_bounds.first = value - taylor_bound_f; + value_bounds.second = value + taylor_bound_f; + } + +} // namespace FunctionTools + +#include "function_tools.inst" + +DEAL_II_NAMESPACE_CLOSE diff --git a/source/base/function_tools.inst.in b/source/base/function_tools.inst.in new file mode 100644 index 0000000000..e8977eb5c4 --- /dev/null +++ b/source/base/function_tools.inst.in @@ -0,0 +1,25 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2015 - 2018 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. +// +// --------------------------------------------------------------------- + +for (deal_II_dimension : DIMENSIONS) + { + template void + FunctionTools::taylor_estimate_function_bounds( + const Function &, + const BoundingBox &, + std::pair &, + std::array, deal_II_dimension> &, + const unsigned int); + } diff --git a/tests/base/taylor_estimate_function_bounds.cc b/tests/base/taylor_estimate_function_bounds.cc new file mode 100644 index 0000000000..070475d25b --- /dev/null +++ b/tests/base/taylor_estimate_function_bounds.cc @@ -0,0 +1,144 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2017 - 2020 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" + +/* + * Test taylor_estimate_function_bounds, by calling it with a function and a box + * such that the function bounds are simple to compute. + */ + +namespace +{ + using namespace dealii; + + /** + * Returns a bounding box with different side lengths, 2*h_i, in all + * directions. Here, h_i = 1/(1 + i), i.e. the box is + * [-1, 1] x [-1/2, 1/2] x [-1/3, 1/3] in 3D. + */ + template + BoundingBox + create_box() + { + std::pair, Point> lower_upper_corner; + for (unsigned int i = 0; i < dim; i++) + { + lower_upper_corner.first[i] = -1.0 / (i + 1); + lower_upper_corner.second[i] = 1.0 / (i + 1); + } + + return BoundingBox(lower_upper_corner); + } + + + + /* + * The Taylor expansion multiplies the gradient entries with h_i, + * and the Hessian entries with h_i * h_j, where h_i is the half the + * side length of the of the box in the ith coordinate direction. Choose the + * gradient and Hessian as + * + * grad[i] = 1 / h_i, + * hessian[i][j] = 1 / (h_i * h_j), + * + * to get function bounds that are simple to compute. + */ + 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> grad; + for (unsigned int i = 0; i < dim; i++) + grad[i] = i + 1; + + return grad; + } + + SymmetricTensor<2, dim> + hessian(const Point &, const unsigned int) const + { + SymmetricTensor<2, dim> hess; + for (unsigned int i = 0; i < dim; i++) + for (unsigned int j = 0; j < dim; j++) + hess[i][j] = (i + 1) * (j + 1); + + return hess; + } + }; + + + + void + print_bounds(const std::pair &bounds) + { + deallog << "[" << bounds.first << ", " << bounds.second << "]" << std::endl; + } + + + + /** + * Call taylor_estimate_function_bounds with the above TestFunction and + * the box returned by create_box. Print the computed bounds to deallog. + */ + template + void + run_test() + { + deallog << dim << "D" << std::endl; + + const TestFunction function; + const BoundingBox box = create_box(); + + std::pair value_bounds; + std::array, dim> gradient_bounds; + FunctionTools::taylor_estimate_function_bounds(function, + box, + value_bounds, + gradient_bounds); + + deallog << "value: "; + print_bounds(value_bounds); + + for (unsigned int i = 0; i < dim; i++) + { + deallog << "gradient[" << i << "]: "; + print_bounds(gradient_bounds[i]); + } + deallog << std::endl; + } +} // namespace + + + +int +main() +{ + initlog(); + run_test<1>(); + run_test<2>(); + run_test<3>(); +} diff --git a/tests/base/taylor_estimate_function_bounds.output b/tests/base/taylor_estimate_function_bounds.output new file mode 100644 index 0000000000..b7618811ea --- /dev/null +++ b/tests/base/taylor_estimate_function_bounds.output @@ -0,0 +1,16 @@ + +DEAL::1D +DEAL::value: [-0.500000, 2.50000] +DEAL::gradient[0]: [0.00000, 2.00000] +DEAL:: +DEAL::2D +DEAL::value: [-3.00000, 5.00000] +DEAL::gradient[0]: [-1.00000, 3.00000] +DEAL::gradient[1]: [-2.00000, 6.00000] +DEAL:: +DEAL::3D +DEAL::value: [-6.50000, 8.50000] +DEAL::gradient[0]: [-2.00000, 4.00000] +DEAL::gradient[1]: [-4.00000, 8.00000] +DEAL::gradient[2]: [-6.00000, 12.0000] +DEAL:: -- 2.39.5