From f08b5f6aacf288733c62e5280c02a46480c2e243 Mon Sep 17 00:00:00 2001
From: David Wells <drwells@email.unc.edu>
Date: Tue, 16 Mar 2021 18:59:25 -0400
Subject: [PATCH] Add Functions::IdentityFunction.

---
 doc/news/changes/minor/20210316DavidWells |   3 +
 include/deal.II/base/function.h           |  45 +++++
 include/deal.II/base/function.templates.h |  53 ++++++
 source/base/function.inst.in              |   2 +
 tests/base/functions_14.cc                | 205 ++++++++++++++++++++
 tests/base/functions_14.output            | 218 ++++++++++++++++++++++
 6 files changed, 526 insertions(+)
 create mode 100644 doc/news/changes/minor/20210316DavidWells
 create mode 100644 tests/base/functions_14.cc
 create mode 100644 tests/base/functions_14.output

diff --git a/doc/news/changes/minor/20210316DavidWells b/doc/news/changes/minor/20210316DavidWells
new file mode 100644
index 0000000000..6984bb9f57
--- /dev/null
+++ b/doc/news/changes/minor/20210316DavidWells
@@ -0,0 +1,3 @@
+New: Added Functions::IdentityFunction.
+<br>
+(David Wells, 2021/03/16)
diff --git a/include/deal.II/base/function.h b/include/deal.II/base/function.h
index f415820400..e3ec0fc127 100644
--- a/include/deal.II/base/function.h
+++ b/include/deal.II/base/function.h
@@ -516,6 +516,51 @@ namespace Functions
     explicit ZeroFunction(const unsigned int n_components = 1);
   };
 
+  /**
+   * A function whose output is also its input. One possible application of this
+   * function is interpolating or projecting a finite element field that
+   * represents spatial coordinates: e.g., one can set up a finite element field
+   * to interpolate the vertices of a Triangulation with this function, which is
+   * useful when doing calculations in a Lagrangian reference frame.
+   *
+   * @ingroup functions
+   */
+  template <int dim, typename RangeNumberType = double>
+  class IdentityFunction : public Function<dim, RangeNumberType>
+  {
+  public:
+    /**
+     * Constructor. The number of components is set to dim.
+     */
+    IdentityFunction();
+
+    /**
+     * @copydoc Function::value()
+     */
+    virtual RangeNumberType
+    value(const Point<dim> &p, const unsigned int component = 0) const override;
+
+    /**
+     * @copydoc Function::gradient()
+     */
+    virtual Tensor<1, dim, RangeNumberType>
+    gradient(const Point<dim> & p,
+             const unsigned int component = 0) const override;
+
+    /**
+     * @copydoc Function::laplacian()
+     */
+    virtual RangeNumberType
+    laplacian(const Point<dim> & p,
+              const unsigned int component = 0) const override;
+
+    /**
+     * @copydoc Function::hessian()
+     */
+    virtual SymmetricTensor<2, dim, RangeNumberType>
+    hessian(const Point<dim> & p,
+            const unsigned int component = 0) const override;
+  };
 } // namespace Functions
 
 
diff --git a/include/deal.II/base/function.templates.h b/include/deal.II/base/function.templates.h
index 313ce62b63..46ded4f4db 100644
--- a/include/deal.II/base/function.templates.h
+++ b/include/deal.II/base/function.templates.h
@@ -533,6 +533,59 @@ namespace Functions
     const unsigned int n_components)
     : ConstantFunction<dim, RangeNumberType>(RangeNumberType(), n_components)
   {}
+
+
+
+  template <int dim, typename RangeNumberType>
+  IdentityFunction<dim, RangeNumberType>::IdentityFunction()
+    : Function<dim, RangeNumberType>(dim)
+  {}
+
+
+
+  template <int dim, typename RangeNumberType>
+  RangeNumberType
+  IdentityFunction<dim, RangeNumberType>::value(
+    const Point<dim> & p,
+    const unsigned int component) const
+  {
+    AssertIndexRange(component, this->n_components);
+    return p[component];
+  }
+
+
+
+  template <int dim, typename RangeNumberType>
+  Tensor<1, dim, RangeNumberType>
+  IdentityFunction<dim, RangeNumberType>::gradient(
+    const Point<dim> &,
+    const unsigned int component) const
+  {
+    AssertIndexRange(component, this->n_components);
+    Tensor<1, dim, RangeNumberType> result;
+    result[component] = RangeNumberType(1);
+    return result;
+  }
+
+
+
+  template <int dim, typename RangeNumberType>
+  SymmetricTensor<2, dim, RangeNumberType>
+  IdentityFunction<dim, RangeNumberType>::hessian(const Point<dim> &,
+                                                  const unsigned int) const
+  {
+    return SymmetricTensor<2, dim, RangeNumberType>();
+  }
+
+
+
+  template <int dim, typename RangeNumberType>
+  RangeNumberType
+  IdentityFunction<dim, RangeNumberType>::laplacian(const Point<dim> &,
+                                                    const unsigned int) const
+  {
+    return 0;
+  }
 } // namespace Functions
 
 //---------------------------------------------------------------------------
diff --git a/source/base/function.inst.in b/source/base/function.inst.in
index 4c1b540e24..ac59661111 100644
--- a/source/base/function.inst.in
+++ b/source/base/function.inst.in
@@ -19,6 +19,7 @@ for (S : REAL_SCALARS; dim : SPACE_DIMENSIONS)
     template class Function<dim, S>;
     template class Functions::ZeroFunction<dim, S>;
     template class Functions::ConstantFunction<dim, S>;
+    template class Functions::IdentityFunction<dim, S>;
     template class ComponentSelectFunction<dim, S>;
     template class ScalarFunctionFromFunctionObject<dim, S>;
     template class VectorFunctionFromScalarFunctionObject<dim, S>;
@@ -31,6 +32,7 @@ for (S : COMPLEX_SCALARS; dim : SPACE_DIMENSIONS)
     template class Function<dim, S>;
     template class Functions::ZeroFunction<dim, S>;
     template class Functions::ConstantFunction<dim, S>;
+    template class Functions::IdentityFunction<dim, S>;
     template class ComponentSelectFunction<dim, S>;
     template class ScalarFunctionFromFunctionObject<dim, S>;
     template class VectorFunctionFromScalarFunctionObject<dim, S>;
diff --git a/tests/base/functions_14.cc b/tests/base/functions_14.cc
new file mode 100644
index 0000000000..43edd767d1
--- /dev/null
+++ b/tests/base/functions_14.cc
@@ -0,0 +1,205 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2021 - 2021 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.
+//
+// ---------------------------------------------------------------------
+
+
+// Check IdentityFunction
+
+#include <deal.II/base/data_out_base.h>
+#include <deal.II/base/function_lib.h>
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/lac/vector.h>
+
+#include <string>
+#include <vector>
+
+#include "../tests.h"
+
+
+
+template <int dim>
+void
+check_value(const Function<dim> &f)
+{
+  for (unsigned int d = 0; d < dim; ++d)
+    {
+      deallog << " d = " << d << std::endl;
+      Point<dim> p;
+      for (unsigned int i = 0; i < dim; i++)
+        p[i] = i;
+
+      deallog << f.value(p, d) << std::endl;
+    }
+  deallog << " values checked" << std::endl;
+}
+
+template <int dim>
+void
+check_value_list(const Function<dim> &f)
+{
+  for (unsigned int d = 0; d < dim; ++d)
+    {
+      deallog << " d = " << d << std::endl;
+      const unsigned int      max_number_of_points = 5;
+      std::vector<Point<dim>> points(max_number_of_points);
+
+      for (unsigned int i = 0; i < max_number_of_points; ++i)
+        {
+          Point<dim> p;
+          for (unsigned int j = 0; j < dim; ++j)
+            p[j] = i + 1;
+
+          points[i] = p;
+        }
+      std::vector<double> values(max_number_of_points);
+      f.value_list(points, values, d);
+
+      for (unsigned int j = 0; j < max_number_of_points; ++j)
+        deallog << values[j] << std::endl;
+    }
+
+  deallog << " value_list checked" << std::endl;
+}
+
+
+template <int dim>
+void
+check_gradient(const Function<dim> &f)
+{
+  for (unsigned int d = 0; d < dim; ++d)
+    {
+      deallog << " d = " << d << std::endl;
+      Point<dim> p;
+      for (unsigned int i = 0; i < dim; i++)
+        p[i] = i;
+
+      deallog << f.gradient(p, d) << std::endl;
+    }
+  deallog << " gradients checked" << std::endl;
+}
+
+template <int dim>
+void
+check_gradient_list(const Function<dim> &f)
+{
+  const unsigned int      max_number_of_points = 5;
+  std::vector<Point<dim>> points(max_number_of_points);
+
+  for (unsigned int d = 0; d < dim; ++d)
+    {
+      deallog << " d = " << d << std::endl;
+      for (unsigned int i = 0; i < max_number_of_points; ++i)
+        {
+          Point<dim> p;
+          for (unsigned int j = 0; j < dim; ++j)
+            p[j] = i + 1;
+
+          points[i] = p;
+        }
+      std::vector<Tensor<1, dim>> tensors(max_number_of_points);
+      f.gradient_list(points, tensors, d);
+
+      for (unsigned int j = 0; j < max_number_of_points; ++j)
+        deallog << tensors[j] << std::endl;
+    }
+
+  deallog << " gradient_list checked" << std::endl;
+}
+
+
+template <int dim>
+void
+check_laplacian(const Function<dim> &f)
+{
+  for (unsigned int d = 0; d < dim; ++d)
+    {
+      deallog << " d = " << d << std::endl;
+      Point<dim> p;
+      for (unsigned int i = 0; i < dim; i++)
+        p[i] = i;
+
+      deallog << f.laplacian(p, d) << std::endl;
+    }
+  deallog << " laplacians checked" << std::endl;
+}
+
+
+template <int dim>
+void
+check_laplacian_list(const Function<dim> &f)
+{
+  for (unsigned int d = 0; d < dim; ++d)
+    {
+      deallog << " d = " << d << std::endl;
+      const unsigned int      max_number_of_points = 5;
+      std::vector<Point<dim>> points(max_number_of_points);
+
+      for (unsigned int i = 0; i < max_number_of_points; ++i)
+        {
+          Point<dim> p;
+          for (unsigned int j = 0; j < dim; ++j)
+            p[j] = i + 1;
+
+          points[i] = p;
+        }
+      std::vector<double> values(max_number_of_points);
+      f.laplacian_list(points, values, d);
+
+      for (unsigned int j = 0; j < max_number_of_points; ++j)
+        deallog << values[j] << std::endl;
+    }
+
+  deallog << " laplacian_list checked" << std::endl;
+}
+
+
+int
+main()
+{
+  initlog();
+
+  deallog << "Functions IdentityFunction" << std::endl;
+  check_value(Functions::IdentityFunction<1>());
+  check_value(Functions::IdentityFunction<2>());
+  check_value(Functions::IdentityFunction<3>());
+
+  check_value_list(Functions::IdentityFunction<1>());
+  check_value_list(Functions::IdentityFunction<2>());
+  check_value_list(Functions::IdentityFunction<3>());
+
+  check_gradient(Functions::IdentityFunction<1>());
+  check_gradient(Functions::IdentityFunction<2>());
+  check_gradient(Functions::IdentityFunction<3>());
+
+  check_gradient_list(Functions::IdentityFunction<1>());
+  check_gradient_list(Functions::IdentityFunction<2>());
+  check_gradient_list(Functions::IdentityFunction<3>());
+
+  check_gradient(Functions::IdentityFunction<1>());
+  check_gradient(Functions::IdentityFunction<2>());
+  check_gradient(Functions::IdentityFunction<3>());
+
+  check_gradient_list(Functions::IdentityFunction<1>());
+  check_gradient_list(Functions::IdentityFunction<2>());
+  check_gradient_list(Functions::IdentityFunction<3>());
+
+  check_laplacian(Functions::IdentityFunction<1>());
+  check_laplacian(Functions::IdentityFunction<2>());
+  check_laplacian(Functions::IdentityFunction<3>());
+
+  check_laplacian_list(Functions::IdentityFunction<1>());
+  check_laplacian_list(Functions::IdentityFunction<2>());
+  check_laplacian_list(Functions::IdentityFunction<3>());
+}
diff --git a/tests/base/functions_14.output b/tests/base/functions_14.output
new file mode 100644
index 0000000000..6eaa59d88f
--- /dev/null
+++ b/tests/base/functions_14.output
@@ -0,0 +1,218 @@
+
+DEAL::Functions IdentityFunction
+DEAL:: d = 0
+DEAL::0.00000
+DEAL:: values checked
+DEAL:: d = 0
+DEAL::0.00000
+DEAL:: d = 1
+DEAL::1.00000
+DEAL:: values checked
+DEAL:: d = 0
+DEAL::0.00000
+DEAL:: d = 1
+DEAL::1.00000
+DEAL:: d = 2
+DEAL::2.00000
+DEAL:: values checked
+DEAL:: d = 0
+DEAL::1.00000
+DEAL::2.00000
+DEAL::3.00000
+DEAL::4.00000
+DEAL::5.00000
+DEAL:: value_list checked
+DEAL:: d = 0
+DEAL::1.00000
+DEAL::2.00000
+DEAL::3.00000
+DEAL::4.00000
+DEAL::5.00000
+DEAL:: d = 1
+DEAL::1.00000
+DEAL::2.00000
+DEAL::3.00000
+DEAL::4.00000
+DEAL::5.00000
+DEAL:: value_list checked
+DEAL:: d = 0
+DEAL::1.00000
+DEAL::2.00000
+DEAL::3.00000
+DEAL::4.00000
+DEAL::5.00000
+DEAL:: d = 1
+DEAL::1.00000
+DEAL::2.00000
+DEAL::3.00000
+DEAL::4.00000
+DEAL::5.00000
+DEAL:: d = 2
+DEAL::1.00000
+DEAL::2.00000
+DEAL::3.00000
+DEAL::4.00000
+DEAL::5.00000
+DEAL:: value_list checked
+DEAL:: d = 0
+DEAL::1.00000
+DEAL:: gradients checked
+DEAL:: d = 0
+DEAL::1.00000 0.00000
+DEAL:: d = 1
+DEAL::0.00000 1.00000
+DEAL:: gradients checked
+DEAL:: d = 0
+DEAL::1.00000 0.00000 0.00000
+DEAL:: d = 1
+DEAL::0.00000 1.00000 0.00000
+DEAL:: d = 2
+DEAL::0.00000 0.00000 1.00000
+DEAL:: gradients checked
+DEAL:: d = 0
+DEAL::1.00000
+DEAL::1.00000
+DEAL::1.00000
+DEAL::1.00000
+DEAL::1.00000
+DEAL:: gradient_list checked
+DEAL:: d = 0
+DEAL::1.00000 0.00000
+DEAL::1.00000 0.00000
+DEAL::1.00000 0.00000
+DEAL::1.00000 0.00000
+DEAL::1.00000 0.00000
+DEAL:: d = 1
+DEAL::0.00000 1.00000
+DEAL::0.00000 1.00000
+DEAL::0.00000 1.00000
+DEAL::0.00000 1.00000
+DEAL::0.00000 1.00000
+DEAL:: gradient_list checked
+DEAL:: d = 0
+DEAL::1.00000 0.00000 0.00000
+DEAL::1.00000 0.00000 0.00000
+DEAL::1.00000 0.00000 0.00000
+DEAL::1.00000 0.00000 0.00000
+DEAL::1.00000 0.00000 0.00000
+DEAL:: d = 1
+DEAL::0.00000 1.00000 0.00000
+DEAL::0.00000 1.00000 0.00000
+DEAL::0.00000 1.00000 0.00000
+DEAL::0.00000 1.00000 0.00000
+DEAL::0.00000 1.00000 0.00000
+DEAL:: d = 2
+DEAL::0.00000 0.00000 1.00000
+DEAL::0.00000 0.00000 1.00000
+DEAL::0.00000 0.00000 1.00000
+DEAL::0.00000 0.00000 1.00000
+DEAL::0.00000 0.00000 1.00000
+DEAL:: gradient_list checked
+DEAL:: d = 0
+DEAL::1.00000
+DEAL:: gradients checked
+DEAL:: d = 0
+DEAL::1.00000 0.00000
+DEAL:: d = 1
+DEAL::0.00000 1.00000
+DEAL:: gradients checked
+DEAL:: d = 0
+DEAL::1.00000 0.00000 0.00000
+DEAL:: d = 1
+DEAL::0.00000 1.00000 0.00000
+DEAL:: d = 2
+DEAL::0.00000 0.00000 1.00000
+DEAL:: gradients checked
+DEAL:: d = 0
+DEAL::1.00000
+DEAL::1.00000
+DEAL::1.00000
+DEAL::1.00000
+DEAL::1.00000
+DEAL:: gradient_list checked
+DEAL:: d = 0
+DEAL::1.00000 0.00000
+DEAL::1.00000 0.00000
+DEAL::1.00000 0.00000
+DEAL::1.00000 0.00000
+DEAL::1.00000 0.00000
+DEAL:: d = 1
+DEAL::0.00000 1.00000
+DEAL::0.00000 1.00000
+DEAL::0.00000 1.00000
+DEAL::0.00000 1.00000
+DEAL::0.00000 1.00000
+DEAL:: gradient_list checked
+DEAL:: d = 0
+DEAL::1.00000 0.00000 0.00000
+DEAL::1.00000 0.00000 0.00000
+DEAL::1.00000 0.00000 0.00000
+DEAL::1.00000 0.00000 0.00000
+DEAL::1.00000 0.00000 0.00000
+DEAL:: d = 1
+DEAL::0.00000 1.00000 0.00000
+DEAL::0.00000 1.00000 0.00000
+DEAL::0.00000 1.00000 0.00000
+DEAL::0.00000 1.00000 0.00000
+DEAL::0.00000 1.00000 0.00000
+DEAL:: d = 2
+DEAL::0.00000 0.00000 1.00000
+DEAL::0.00000 0.00000 1.00000
+DEAL::0.00000 0.00000 1.00000
+DEAL::0.00000 0.00000 1.00000
+DEAL::0.00000 0.00000 1.00000
+DEAL:: gradient_list checked
+DEAL:: d = 0
+DEAL::0.00000
+DEAL:: laplacians checked
+DEAL:: d = 0
+DEAL::0.00000
+DEAL:: d = 1
+DEAL::0.00000
+DEAL:: laplacians checked
+DEAL:: d = 0
+DEAL::0.00000
+DEAL:: d = 1
+DEAL::0.00000
+DEAL:: d = 2
+DEAL::0.00000
+DEAL:: laplacians checked
+DEAL:: d = 0
+DEAL::0.00000
+DEAL::0.00000
+DEAL::0.00000
+DEAL::0.00000
+DEAL::0.00000
+DEAL:: laplacian_list checked
+DEAL:: d = 0
+DEAL::0.00000
+DEAL::0.00000
+DEAL::0.00000
+DEAL::0.00000
+DEAL::0.00000
+DEAL:: d = 1
+DEAL::0.00000
+DEAL::0.00000
+DEAL::0.00000
+DEAL::0.00000
+DEAL::0.00000
+DEAL:: laplacian_list checked
+DEAL:: d = 0
+DEAL::0.00000
+DEAL::0.00000
+DEAL::0.00000
+DEAL::0.00000
+DEAL::0.00000
+DEAL:: d = 1
+DEAL::0.00000
+DEAL::0.00000
+DEAL::0.00000
+DEAL::0.00000
+DEAL::0.00000
+DEAL:: d = 2
+DEAL::0.00000
+DEAL::0.00000
+DEAL::0.00000
+DEAL::0.00000
+DEAL::0.00000
+DEAL:: laplacian_list checked
-- 
2.39.5