--- /dev/null
+New: Add a function taylor_estimate_function_bounds, which estimates the range
+of the value and gradient components of a Function<dim> over a
+BoundingBox<dim>, by approximating the function by a 2nd order Taylor
+polynomial.
+<br>
+(Simon Sticko, 2021/04/25)
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/config.h>
+
+#include <deal.II/base/bounding_box.h>
+#include <deal.II/base/function.h>
+
+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
+ * <code>std::pair<double, double></code>, 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 <int dim>
+ void
+ taylor_estimate_function_bounds(
+ const Function<dim> & function,
+ const BoundingBox<dim> & box,
+ std::pair<double, double> & value_bounds,
+ std::array<std::pair<double, double>, dim> &gradient_bounds,
+ const unsigned int component = 0);
+
+} // namespace FunctionTools
+DEAL_II_NAMESPACE_CLOSE
+
+#endif /* dealii_function_tools_h */
function_parser.cc
function_spherical.cc
function_time.cc
+ function_tools.cc
geometry_info.cc
geometric_utilities.cc
graph_coloring.cc
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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/function_tools.h>
+
+#include <cmath>
+
+
+DEAL_II_NAMESPACE_OPEN
+namespace FunctionTools
+{
+ template <int dim>
+ void
+ taylor_estimate_function_bounds(
+ const Function<dim> & function,
+ const BoundingBox<dim> & box,
+ std::pair<double, double> & value_bounds,
+ std::array<std::pair<double, double>, dim> &gradient_bounds,
+ const unsigned int component)
+ {
+ const Point<dim> 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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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<deal_II_dimension>(
+ const Function<deal_II_dimension> &,
+ const BoundingBox<deal_II_dimension> &,
+ std::pair<double, double> &,
+ std::array<std::pair<double, double>, deal_II_dimension> &,
+ const unsigned int);
+ }
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/function_tools.h>
+
+#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 <int dim>
+ BoundingBox<dim>
+ create_box()
+ {
+ std::pair<Point<dim>, Point<dim>> 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<dim>(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 <int dim>
+ class TestFunction : public Function<dim>
+ {
+ public:
+ double
+ value(const Point<dim> &, const unsigned int) const override
+ {
+ return 1;
+ }
+
+ Tensor<1, dim>
+ gradient(const Point<dim> &, 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<dim> &, 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<double, double> &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 <int dim>
+ void
+ run_test()
+ {
+ deallog << dim << "D" << std::endl;
+
+ const TestFunction<dim> function;
+ const BoundingBox<dim> box = create_box<dim>();
+
+ std::pair<double, double> value_bounds;
+ std::array<std::pair<double, double>, dim> gradient_bounds;
+ FunctionTools::taylor_estimate_function_bounds<dim>(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>();
+}
--- /dev/null
+
+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::
projects / dealii.git / commitdiff