]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add classes restricting a function to a lower dimension 12015/head
authorSimon Sticko <simon@sticko.se>
Mon, 5 Apr 2021 11:34:10 +0000 (13:34 +0200)
committerSimon Sticko <simon@sticko.se>
Wed, 28 Apr 2021 05:57:30 +0000 (07:57 +0200)
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 [new file with mode: 0644]
source/base/CMakeLists.txt
source/base/function_restriction.cc [new file with mode: 0644]
source/base/function_restriction.inst.in [new file with mode: 0644]
tests/base/coordinate_restriction.cc [new file with mode: 0644]
tests/base/coordinate_restriction.output [new file with mode: 0644]
tests/base/create_higher_dim_point.cc [new file with mode: 0644]
tests/base/create_higher_dim_point.output [new file with mode: 0644]
tests/base/point_restriction.cc [new file with mode: 0644]
tests/base/point_restriction.output [new file with mode: 0644]

diff --git a/include/deal.II/base/function_restriction.h b/include/deal.II/base/function_restriction.h
new file mode 100644 (file)
index 0000000..e96c4a0
--- /dev/null
@@ -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 <deal.II/base/function.h>
+#include <deal.II/base/point.h>
+#include <deal.II/base/smartpointer.h>
+#include <deal.II/base/tensor.h>
+
+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 <int dim>
+  class CoordinateRestriction : public Function<dim>
+  {
+  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<dim + 1> &function,
+                          const unsigned int       direction,
+                          const double             coordinate_value);
+
+    double
+    value(const Point<dim> &point, const unsigned int component) const override;
+
+    Tensor<1, dim>
+    gradient(const Point<dim> & point,
+             const unsigned int component) const override;
+
+    SymmetricTensor<2, dim>
+    hessian(const Point<dim> & point,
+            const unsigned int component) const override;
+
+  private:
+    // The higher-dimensional function that has been restricted.
+    const SmartPointer<const Function<dim + 1>> 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 <int dim>
+  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<dim + 1> &function,
+                     const unsigned int       open_direction,
+                     const Point<dim> &       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<const Function<dim + 1>> 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<dim> 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 <int dim>
+  Point<dim + 1>
+  create_higher_dim_point(const Point<dim> & point,
+                          const unsigned int component_in_dim_plus_1,
+                          const double       coordinate_value);
+} // namespace internal
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif /* dealii_function_restriction_h */
index 024d2a9b023084c0bd1b6bf2484959773ddb16d5..15551efd42df28d2762419f19223e594253ea4fb 100644 (file)
@@ -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 (file)
index 0000000..71d9bf6
--- /dev/null
@@ -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 <deal.II/base/bounding_box.h>
+#include <deal.II/base/function_restriction.h>
+
+DEAL_II_NAMESPACE_OPEN
+namespace internal
+{
+  template <int dim>
+  Point<dim + 1>
+  create_higher_dim_point(const Point<dim> & point,
+                          const unsigned int component_in_dim_plus_1,
+                          const double       coordinate_value)
+  {
+    AssertIndexRange(component_in_dim_plus_1, dim + 1);
+
+    Point<dim + 1> 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<dim>(
+            component_in_dim_plus_1, d);
+        output(component_to_write_to) = point(d);
+      }
+
+    return output;
+  }
+} // namespace internal
+
+
+
+namespace Functions
+{
+  template <int dim>
+  CoordinateRestriction<dim>::CoordinateRestriction(
+    const Function<dim + 1> &function,
+    const unsigned int       direction,
+    const double             coordinate_value)
+    : function(&function)
+    , restricted_direction(direction)
+    , coordinate_value(coordinate_value)
+  {
+    AssertIndexRange(restricted_direction, dim + 1);
+  }
+
+
+
+  template <int dim>
+  double
+  CoordinateRestriction<dim>::value(const Point<dim> & point,
+                                    const unsigned int component) const
+  {
+    const Point<dim + 1> full_point =
+      internal::create_higher_dim_point(point,
+                                        restricted_direction,
+                                        coordinate_value);
+
+    return function->value(full_point, component);
+  }
+
+
+
+  template <int dim>
+  Tensor<1, dim>
+  CoordinateRestriction<dim>::gradient(const Point<dim> & point,
+                                       const unsigned int component) const
+  {
+    const Point<dim + 1> 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<dim>(restricted_direction, d);
+        grad[d] = full_gradient[index_to_write_from];
+      }
+    return grad;
+  }
+
+
+
+  template <int dim>
+  SymmetricTensor<2, dim>
+  CoordinateRestriction<dim>::hessian(const Point<dim> & point,
+                                      const unsigned int component) const
+  {
+    const Point<dim + 1> 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<dim>(restricted_direction, i);
+        for (unsigned int j = 0; j < dim; ++j)
+          {
+            const unsigned int j_to_write_from =
+              internal::coordinate_to_one_dim_higher<dim>(restricted_direction,
+                                                          j);
+            hess[i][j] = full_hessian[i_to_write_from][j_to_write_from];
+          }
+      }
+    return hess;
+  }
+
+
+
+  template <int dim>
+  PointRestriction<dim>::PointRestriction(const Function<dim + 1> &function,
+                                          const unsigned int open_direction,
+                                          const Point<dim> & point)
+    : function(&function)
+    , open_direction(open_direction)
+    , point(point)
+  {
+    AssertIndexRange(open_direction, dim + 1);
+  }
+
+
+
+  template <int dim>
+  double
+  PointRestriction<dim>::value(const Point<1> &   point_1D,
+                               const unsigned int component) const
+  {
+    const Point<dim + 1> full_point =
+      internal::create_higher_dim_point(point, open_direction, point_1D(0));
+    return function->value(full_point, component);
+  }
+
+
+
+  template <int dim>
+  Tensor<1, 1>
+  PointRestriction<dim>::gradient(const Point<1> &   point_1D,
+                                  const unsigned int component) const
+  {
+    const Point<dim + 1> 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 <int dim>
+  SymmetricTensor<2, 1>
+  PointRestriction<dim>::hessian(const Point<1> &   point_1D,
+                                 const unsigned int component) const
+  {
+    const Point<dim + 1> 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 (file)
index 0000000..7ce0681
--- /dev/null
@@ -0,0 +1,11 @@
+for (deal_II_dimension : DIMENSIONS)
+  {
+    template class Functions::PointRestriction<deal_II_dimension>;
+
+    template Point<deal_II_dimension> internal::create_higher_dim_point(
+      const Point<deal_II_dimension - 1> &, const unsigned int, const double);
+
+#if deal_II_dimension != 3
+    template class Functions::CoordinateRestriction<deal_II_dimension>;
+#endif
+  }
diff --git a/tests/base/coordinate_restriction.cc b/tests/base/coordinate_restriction.cc
new file mode 100644 (file)
index 0000000..2bc454c
--- /dev/null
@@ -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 <deal.II/base/function_restriction.h>
+
+#include "../tests.h"
+
+namespace
+{
+  using namespace dealii;
+
+  /*
+   * A simple test-function which is constant in space, but varies
+   * in each component.
+   */
+  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> tensor;
+      for (unsigned int i = 0; i < dim; ++i)
+        tensor[i] = i + 1;
+      return tensor;
+    }
+
+    SymmetricTensor<2, dim>
+    hessian(const Point<dim> &, 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 <int dim>
+  void
+  test_truncates_derivatives_correctly(const unsigned int restricted_direction)
+  {
+    deallog << "Restricted direction = " << restricted_direction << std::endl;
+
+    const TestFunction<dim> function;
+    const double            coordinate_value = 0;
+
+    const Functions::CoordinateRestriction<dim - 1> restriction(
+      function, restricted_direction, coordinate_value);
+
+    const Point<dim - 1> 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 <int dim>
+  void
+  test_all_coordinate_directions()
+  {
+    deallog << "dim = " << dim << std::endl;
+    for (unsigned int direction = 0; direction < dim; direction++)
+      test_truncates_derivatives_correctly<dim>(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 (file)
index 0000000..98470a1
--- /dev/null
@@ -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 (file)
index 0000000..0846191
--- /dev/null
@@ -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 <deal.II/base/function_restriction.h>
+
+#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<dim - 1> low_dim_point(y, z);
+      const Point<dim>     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<dim - 1> low_dim_point(z, x);
+      const Point<dim>     point =
+        internal::create_higher_dim_point(low_dim_point, component, y);
+      deallog << "Point: " << point << std::endl;
+    }
+    {
+      const unsigned int component = 2;
+      Point<dim - 1>     low_dim_point(x, y);
+      Point<dim>         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<dim - 1> low_dim_point(y);
+      const Point<dim>     point =
+        internal::create_higher_dim_point(low_dim_point, component, x);
+      deallog << "Point: " << point << std::endl;
+    }
+    {
+      const unsigned int   component = 1;
+      const Point<dim - 1> low_dim_point(x);
+      const Point<dim>     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<dim - 1> low_dim_point;
+    const Point<dim>     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 (file)
index 0000000..4f2d84f
--- /dev/null
@@ -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 (file)
index 0000000..abe76ee
--- /dev/null
@@ -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 <deal.II/base/function_restriction.h>
+
+#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 <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> tensor;
+      for (unsigned int i = 0; i < dim; ++i)
+        tensor[i] = i + 1;
+      return tensor;
+    }
+
+    SymmetricTensor<2, dim>
+    hessian(const Point<dim> &, 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 <int dim>
+  void
+  test_truncates_derivatives_correctly(const unsigned int open_direction)
+  {
+    deallog << "open direction = " << open_direction << std::endl;
+    const Point<dim - 1>    locked_point;
+    const TestFunction<dim> function;
+
+    const Functions::PointRestriction<dim - 1> 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 <int dim>
+  void
+  test_all_coordinate_directions()
+  {
+    deallog << "dim = " << dim << std::endl;
+
+    for (unsigned int direction = 0; direction < dim; direction++)
+      test_truncates_derivatives_correctly<dim>(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 (file)
index 0000000..4c0411b
--- /dev/null
@@ -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::

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.