]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add function estimating Function bounds over a BoundingBox 12062/head
authorSimon Sticko <simon@sticko.se>
Tue, 20 Apr 2021 06:28:23 +0000 (08:28 +0200)
committerSimon Sticko <simon@sticko.se>
Sun, 25 Apr 2021 08:34:05 +0000 (10:34 +0200)
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 [new file with mode: 0644]
include/deal.II/base/function_tools.h [new file with mode: 0644]
source/base/CMakeLists.txt
source/base/function_tools.cc [new file with mode: 0644]
source/base/function_tools.inst.in [new file with mode: 0644]
tests/base/taylor_estimate_function_bounds.cc [new file with mode: 0644]
tests/base/taylor_estimate_function_bounds.output [new file with mode: 0644]

diff --git a/doc/news/changes/minor/20210425SimonSticko b/doc/news/changes/minor/20210425SimonSticko
new file mode 100644 (file)
index 0000000..cafb3b2
--- /dev/null
@@ -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<dim> over a
+BoundingBox<dim>, by approximating the function by a 2nd order Taylor
+polynomial.
+<br>
+(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 (file)
index 0000000..9f96cb1
--- /dev/null
@@ -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 <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 */
index 6cc778da4bb17374fcb610032d2c4d639103aff0..f8ce1e3469615bd86e6b2ef4b39ac3d3923c870f 100644 (file)
@@ -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 (file)
index 0000000..bc7b271
--- /dev/null
@@ -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 <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
diff --git a/source/base/function_tools.inst.in b/source/base/function_tools.inst.in
new file mode 100644 (file)
index 0000000..e8977eb
--- /dev/null
@@ -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<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);
+  }
diff --git a/tests/base/taylor_estimate_function_bounds.cc b/tests/base/taylor_estimate_function_bounds.cc
new file mode 100644 (file)
index 0000000..070475d
--- /dev/null
@@ -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 <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>();
+}
diff --git a/tests/base/taylor_estimate_function_bounds.output b/tests/base/taylor_estimate_function_bounds.output
new file mode 100644 (file)
index 0000000..b761881
--- /dev/null
@@ -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::

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.