Add test case

diff --git a/tests/matrix_free/tensor_product_evaluate_05.cc b/tests/matrix_free/tensor_product_evaluate_05.cc
new file mode 100644 (file)
index 0000000..745013c
--- /dev/null
+++ b/tests/matrix_free/tensor_product_evaluate_05.cc
@@ -0,0 +1,102 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// Tests point-wise evaluation of functions with
+// evaluate_tensor_product_hessian for a scalar function on FE_Q
+
+#include <deal.II/base/point.h>
+#include <deal.II/base/polynomial.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/symmetric_tensor.h>
+#include <deal.II/base/vectorization.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_tools.h>
+
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/matrix_free/tensor_product_kernels.h>
+
+#include "../tests.h"
+
+template <int dim>
+void
+test(const unsigned int degree)
+{
+  FE_Q<dim> fe(degree);
+
+  // choose a symmetric matrix and then construct f(x) = x^T A x
+  SymmetricTensor<2, dim> matrix = unit_symmetric_tensor<dim>();
+  for (unsigned int d = 0; d < dim; ++d)
+    for (unsigned int e = 0; e < dim; ++e)
+      matrix[d][e] += 0.1 * (d + 1 + 2 * e);
+
+  std::vector<double> coefficients(fe.dofs_per_cell);
+  for (unsigned int i = 0; i < fe.dofs_per_cell; ++i)
+    coefficients[i] = fe.get_unit_support_points()[i] *
+                      (matrix * fe.get_unit_support_points()[i]);
+
+  const std::vector<Polynomials::Polynomial<double>> polynomials =
+    Polynomials::generate_complete_Lagrange_basis(
+      QGaussLobatto<1>(degree + 1).get_points());
+
+  const std::vector<unsigned int> renumbering =
+    FETools::lexicographic_to_hierarchic_numbering<dim>(degree);
+
+  const std::vector<Point<dim>> evaluation_points =
+    dim == 3 ? QGauss<dim>(2).get_points() :
+               QIterated<dim>(QTrapez<1>(), 3).get_points();
+
+  deallog << "Evaluate in " << dim << "d with polynomial degree " << degree
+          << std::endl;
+  for (const auto &p : evaluation_points)
+    {
+      Point<dim, VectorizedArray<double>> p_vec;
+      for (unsigned int v = 0; v < VectorizedArray<double>::size(); ++v)
+        for (unsigned int d = 0; d < dim; ++d)
+          p_vec[d][v] = p[d] + 0.01 * v;
+
+      const auto hess = internal::evaluate_tensor_product_hessian(polynomials,
+                                                                  coefficients,
+                                                                  p_vec,
+                                                                  renumbering);
+
+      double error = 0;
+      for (unsigned int v = 0; v < VectorizedArray<double>::size(); ++v)
+        for (unsigned int d = 0; d < dim; ++d)
+          for (unsigned int e = 0; e < dim; ++e)
+            error += std::abs(hess[d][e][v] - 2. * matrix[e][d]);
+
+      deallog << "Hessian error " << error << std::endl;
+    }
+}
+
+
+
+int
+main()
+{
+  initlog();
+  deallog << std::setprecision(9);
+
+  for (unsigned int degree = 1; degree < 5; ++degree)
+    test<1>(degree);
+  for (unsigned int degree = 1; degree < 5; ++degree)
+    test<2>(degree);
+  for (unsigned int degree = 1; degree < 5; ++degree)
+    test<3>(degree);
+}

diff --git a/tests/matrix_free/tensor_product_evaluate_05.output b/tests/matrix_free/tensor_product_evaluate_05.output
new file mode 100644 (file)
index 0000000..731cec2
--- /dev/null
+++ b/tests/matrix_free/tensor_product_evaluate_05.output
@@ -0,0 +1,125 @@
+
+DEAL::Evaluate in 1d with polynomial degree 1
+DEAL::Hessian error 8.80000000
+DEAL::Hessian error 8.80000000
+DEAL::Hessian error 8.80000000
+DEAL::Hessian error 8.80000000
+DEAL::Evaluate in 1d with polynomial degree 2
+DEAL::Hessian error 0.00000000
+DEAL::Hessian error 0.00000000
+DEAL::Hessian error 0.00000000
+DEAL::Hessian error 0.00000000
+DEAL::Evaluate in 1d with polynomial degree 3
+DEAL::Hessian error 7.99360578e-15
+DEAL::Hessian error 2.22044605e-15
+DEAL::Hessian error 7.99360578e-15
+DEAL::Hessian error 1.46549439e-14
+DEAL::Evaluate in 1d with polynomial degree 4
+DEAL::Hessian error 1.99840144e-14
+DEAL::Hessian error 3.55271368e-15
+DEAL::Hessian error 7.10542736e-15
+DEAL::Hessian error 2.70894418e-14
+DEAL::Evaluate in 2d with polynomial degree 1
+DEAL::Hessian error 20.0000000
+DEAL::Hessian error 20.0000000
+DEAL::Hessian error 20.0000000
+DEAL::Hessian error 20.0000000
+DEAL::Hessian error 20.0000000
+DEAL::Hessian error 20.0000000
+DEAL::Hessian error 20.0000000
+DEAL::Hessian error 20.0000000
+DEAL::Hessian error 20.0000000
+DEAL::Hessian error 20.0000000
+DEAL::Hessian error 20.0000000
+DEAL::Hessian error 20.0000000
+DEAL::Hessian error 20.0000000
+DEAL::Hessian error 20.0000000
+DEAL::Hessian error 20.0000000
+DEAL::Hessian error 20.0000000
+DEAL::Evaluate in 2d with polynomial degree 2
+DEAL::Hessian error 4.88498131e-15
+DEAL::Hessian error 5.77315973e-15
+DEAL::Hessian error 7.10542736e-15
+DEAL::Hessian error 1.55431223e-14
+DEAL::Hessian error 5.77315973e-15
+DEAL::Hessian error 6.21724894e-15
+DEAL::Hessian error 9.32587341e-15
+DEAL::Hessian error 1.33226763e-14
+DEAL::Hessian error 1.22124533e-14
+DEAL::Hessian error 5.77315973e-15
+DEAL::Hessian error 1.17683641e-14
+DEAL::Hessian error 1.15463195e-14
+DEAL::Hessian error 1.82076576e-14
+DEAL::Hessian error 8.88178420e-15
+DEAL::Hessian error 1.53210777e-14
+DEAL::Hessian error 1.68753900e-14
+DEAL::Evaluate in 2d with polynomial degree 3
+DEAL::Hessian error 3.01980663e-14
+DEAL::Hessian error 4.61852778e-14
+DEAL::Hessian error 7.10542736e-14
+DEAL::Hessian error 9.59232693e-14
+DEAL::Hessian error 3.33066907e-14
+DEAL::Hessian error 1.75415238e-14
+DEAL::Hessian error 2.53130850e-14
+DEAL::Hessian error 8.83737528e-14
+DEAL::Hessian error 8.03801470e-14
+DEAL::Hessian error 3.28626015e-14
+DEAL::Hessian error 7.26085858e-14
+DEAL::Hessian error 1.62314606e-13
+DEAL::Hessian error 9.68114477e-14
+DEAL::Hessian error 1.28785871e-13
+DEAL::Hessian error 1.90514271e-13
+DEAL::Hessian error 2.10942375e-13
+DEAL::Evaluate in 2d with polynomial degree 4
+DEAL::Hessian error 1.00364161e-13
+DEAL::Hessian error 1.23678845e-13
+DEAL::Hessian error 1.19904087e-13
+DEAL::Hessian error 4.53193039e-13
+DEAL::Hessian error 1.16795462e-13
+DEAL::Hessian error 3.90798505e-14
+DEAL::Hessian error 5.50670620e-14
+DEAL::Hessian error 3.40172335e-13
+DEAL::Hessian error 1.88737914e-13
+DEAL::Hessian error 4.64073224e-14
+DEAL::Hessian error 6.01740879e-14
+DEAL::Hessian error 3.16191517e-13
+DEAL::Hessian error 3.14415161e-13
+DEAL::Hessian error 2.36477504e-13
+DEAL::Hessian error 2.56905608e-13
+DEAL::Hessian error 1.26387789e-12
+DEAL::Evaluate in 3d with polynomial degree 1
+DEAL::Hessian error 33.6000000
+DEAL::Hessian error 33.6000000
+DEAL::Hessian error 33.6000000
+DEAL::Hessian error 33.6000000
+DEAL::Hessian error 33.6000000
+DEAL::Hessian error 33.6000000
+DEAL::Hessian error 33.6000000
+DEAL::Hessian error 33.6000000
+DEAL::Evaluate in 3d with polynomial degree 2
+DEAL::Hessian error 4.84057239e-14
+DEAL::Hessian error 3.81916720e-14
+DEAL::Hessian error 6.26165786e-14
+DEAL::Hessian error 8.94839758e-14
+DEAL::Hessian error 4.77395901e-14
+DEAL::Hessian error 5.12923037e-14
+DEAL::Hessian error 4.88498131e-14
+DEAL::Hessian error 9.52571355e-14
+DEAL::Evaluate in 3d with polynomial degree 3
+DEAL::Hessian error 8.70414851e-14
+DEAL::Hessian error 1.37667655e-13
+DEAL::Hessian error 1.63646874e-13
+DEAL::Hessian error 3.19744231e-13
+DEAL::Hessian error 2.10276241e-13
+DEAL::Hessian error 3.64153152e-13
+DEAL::Hessian error 3.00204306e-13
+DEAL::Hessian error 5.62661029e-13
+DEAL::Evaluate in 3d with polynomial degree 4
+DEAL::Hessian error 8.45989945e-14
+DEAL::Hessian error 2.23820962e-13
+DEAL::Hessian error 2.29816166e-13
+DEAL::Hessian error 3.69038133e-13
+DEAL::Hessian error 2.57571742e-13
+DEAL::Hessian error 3.51718654e-13
+DEAL::Hessian error 3.60378394e-13
+DEAL::Hessian error 3.74811293e-13

diff --git a/tests/matrix_free/tensor_product_evaluate_06.cc b/tests/matrix_free/tensor_product_evaluate_06.cc
new file mode 100644 (file)
index 0000000..883736f
--- /dev/null
+++ b/tests/matrix_free/tensor_product_evaluate_06.cc
@@ -0,0 +1,99 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// Tests point-wise evaluation of functions with
+// evaluate_tensor_product_hessian for a scalar function on FE_DGQ
+
+#include <deal.II/base/point.h>
+#include <deal.II/base/polynomial.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/symmetric_tensor.h>
+#include <deal.II/base/vectorization.h>
+
+#include <deal.II/fe/fe_dgq.h>
+
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/matrix_free/tensor_product_kernels.h>
+
+#include "../tests.h"
+
+template <int dim>
+void
+test(const unsigned int degree)
+{
+  FE_DGQ<dim> fe(degree);
+
+  // choose a symmetric matrix and then construct f(x) = x^T A x
+  SymmetricTensor<2, dim> matrix = unit_symmetric_tensor<dim>();
+  for (unsigned int d = 0; d < dim; ++d)
+    for (unsigned int e = 0; e < dim; ++e)
+      matrix[d][e] += 0.1 * (d + 1 + 2 * e);
+
+  std::vector<double> coefficients(fe.dofs_per_cell);
+  for (unsigned int i = 0; i < fe.dofs_per_cell; ++i)
+    coefficients[i] = fe.get_unit_support_points()[i] *
+                      (matrix * fe.get_unit_support_points()[i]);
+
+  const std::vector<Polynomials::Polynomial<double>> polynomials =
+    Polynomials::generate_complete_Lagrange_basis(
+      QGaussLobatto<1>(degree + 1).get_points());
+
+  const std::vector<Point<dim>> evaluation_points =
+    dim == 3 ? QGauss<dim>(2).get_points() :
+               QIterated<dim>(QTrapez<1>(), 3).get_points();
+
+  deallog << "Evaluate in " << dim << "d with polynomial degree " << degree
+          << std::endl;
+  for (const auto &p : evaluation_points)
+    {
+      Point<dim, VectorizedArray<double>> p_vec;
+      for (unsigned int v = 0; v < VectorizedArray<double>::size(); ++v)
+        for (unsigned int d = 0; d < dim; ++d)
+          p_vec[d][v] = p[d] + 0.01 * v;
+
+      const auto hess = internal::evaluate_tensor_product_hessian(polynomials,
+                                                                  coefficients,
+                                                                  p_vec);
+
+      std::cout << hess << "    " << matrix << std::endl;
+
+      double error = 0;
+      for (unsigned int v = 0; v < VectorizedArray<double>::size(); ++v)
+        for (unsigned int d = 0; d < dim; ++d)
+          for (unsigned int e = 0; e < dim; ++e)
+            error += std::abs(hess[d][e][v] - 2. * matrix[e][d]);
+
+      deallog << "Hessian error " << error << std::endl;
+    }
+}
+
+
+
+int
+main()
+{
+  initlog();
+  deallog << std::setprecision(9);
+
+  for (unsigned int degree = 1; degree < 5; ++degree)
+    test<1>(degree);
+  for (unsigned int degree = 1; degree < 5; ++degree)
+    test<2>(degree);
+  for (unsigned int degree = 1; degree < 5; ++degree)
+    test<3>(degree);
+}

diff --git a/tests/matrix_free/tensor_product_evaluate_06.output b/tests/matrix_free/tensor_product_evaluate_06.output
new file mode 100644 (file)
index 0000000..731cec2
--- /dev/null
+++ b/tests/matrix_free/tensor_product_evaluate_06.output
@@ -0,0 +1,125 @@
+
+DEAL::Evaluate in 1d with polynomial degree 1
+DEAL::Hessian error 8.80000000
+DEAL::Hessian error 8.80000000
+DEAL::Hessian error 8.80000000
+DEAL::Hessian error 8.80000000
+DEAL::Evaluate in 1d with polynomial degree 2
+DEAL::Hessian error 0.00000000
+DEAL::Hessian error 0.00000000
+DEAL::Hessian error 0.00000000
+DEAL::Hessian error 0.00000000
+DEAL::Evaluate in 1d with polynomial degree 3
+DEAL::Hessian error 7.99360578e-15
+DEAL::Hessian error 2.22044605e-15
+DEAL::Hessian error 7.99360578e-15
+DEAL::Hessian error 1.46549439e-14
+DEAL::Evaluate in 1d with polynomial degree 4
+DEAL::Hessian error 1.99840144e-14
+DEAL::Hessian error 3.55271368e-15
+DEAL::Hessian error 7.10542736e-15
+DEAL::Hessian error 2.70894418e-14
+DEAL::Evaluate in 2d with polynomial degree 1
+DEAL::Hessian error 20.0000000
+DEAL::Hessian error 20.0000000
+DEAL::Hessian error 20.0000000
+DEAL::Hessian error 20.0000000
+DEAL::Hessian error 20.0000000
+DEAL::Hessian error 20.0000000
+DEAL::Hessian error 20.0000000
+DEAL::Hessian error 20.0000000
+DEAL::Hessian error 20.0000000
+DEAL::Hessian error 20.0000000
+DEAL::Hessian error 20.0000000
+DEAL::Hessian error 20.0000000
+DEAL::Hessian error 20.0000000
+DEAL::Hessian error 20.0000000
+DEAL::Hessian error 20.0000000
+DEAL::Hessian error 20.0000000
+DEAL::Evaluate in 2d with polynomial degree 2
+DEAL::Hessian error 4.88498131e-15
+DEAL::Hessian error 5.77315973e-15
+DEAL::Hessian error 7.10542736e-15
+DEAL::Hessian error 1.55431223e-14
+DEAL::Hessian error 5.77315973e-15
+DEAL::Hessian error 6.21724894e-15
+DEAL::Hessian error 9.32587341e-15
+DEAL::Hessian error 1.33226763e-14
+DEAL::Hessian error 1.22124533e-14
+DEAL::Hessian error 5.77315973e-15
+DEAL::Hessian error 1.17683641e-14
+DEAL::Hessian error 1.15463195e-14
+DEAL::Hessian error 1.82076576e-14
+DEAL::Hessian error 8.88178420e-15
+DEAL::Hessian error 1.53210777e-14
+DEAL::Hessian error 1.68753900e-14
+DEAL::Evaluate in 2d with polynomial degree 3
+DEAL::Hessian error 3.01980663e-14
+DEAL::Hessian error 4.61852778e-14
+DEAL::Hessian error 7.10542736e-14
+DEAL::Hessian error 9.59232693e-14
+DEAL::Hessian error 3.33066907e-14
+DEAL::Hessian error 1.75415238e-14
+DEAL::Hessian error 2.53130850e-14
+DEAL::Hessian error 8.83737528e-14
+DEAL::Hessian error 8.03801470e-14
+DEAL::Hessian error 3.28626015e-14
+DEAL::Hessian error 7.26085858e-14
+DEAL::Hessian error 1.62314606e-13
+DEAL::Hessian error 9.68114477e-14
+DEAL::Hessian error 1.28785871e-13
+DEAL::Hessian error 1.90514271e-13
+DEAL::Hessian error 2.10942375e-13
+DEAL::Evaluate in 2d with polynomial degree 4
+DEAL::Hessian error 1.00364161e-13
+DEAL::Hessian error 1.23678845e-13
+DEAL::Hessian error 1.19904087e-13
+DEAL::Hessian error 4.53193039e-13
+DEAL::Hessian error 1.16795462e-13
+DEAL::Hessian error 3.90798505e-14
+DEAL::Hessian error 5.50670620e-14
+DEAL::Hessian error 3.40172335e-13
+DEAL::Hessian error 1.88737914e-13
+DEAL::Hessian error 4.64073224e-14
+DEAL::Hessian error 6.01740879e-14
+DEAL::Hessian error 3.16191517e-13
+DEAL::Hessian error 3.14415161e-13
+DEAL::Hessian error 2.36477504e-13
+DEAL::Hessian error 2.56905608e-13
+DEAL::Hessian error 1.26387789e-12
+DEAL::Evaluate in 3d with polynomial degree 1
+DEAL::Hessian error 33.6000000
+DEAL::Hessian error 33.6000000
+DEAL::Hessian error 33.6000000
+DEAL::Hessian error 33.6000000
+DEAL::Hessian error 33.6000000
+DEAL::Hessian error 33.6000000
+DEAL::Hessian error 33.6000000
+DEAL::Hessian error 33.6000000
+DEAL::Evaluate in 3d with polynomial degree 2
+DEAL::Hessian error 4.84057239e-14
+DEAL::Hessian error 3.81916720e-14
+DEAL::Hessian error 6.26165786e-14
+DEAL::Hessian error 8.94839758e-14
+DEAL::Hessian error 4.77395901e-14
+DEAL::Hessian error 5.12923037e-14
+DEAL::Hessian error 4.88498131e-14
+DEAL::Hessian error 9.52571355e-14
+DEAL::Evaluate in 3d with polynomial degree 3
+DEAL::Hessian error 8.70414851e-14
+DEAL::Hessian error 1.37667655e-13
+DEAL::Hessian error 1.63646874e-13
+DEAL::Hessian error 3.19744231e-13
+DEAL::Hessian error 2.10276241e-13
+DEAL::Hessian error 3.64153152e-13
+DEAL::Hessian error 3.00204306e-13
+DEAL::Hessian error 5.62661029e-13
+DEAL::Evaluate in 3d with polynomial degree 4
+DEAL::Hessian error 8.45989945e-14
+DEAL::Hessian error 2.23820962e-13
+DEAL::Hessian error 2.29816166e-13
+DEAL::Hessian error 3.69038133e-13
+DEAL::Hessian error 2.57571742e-13
+DEAL::Hessian error 3.51718654e-13
+DEAL::Hessian error 3.60378394e-13
+DEAL::Hessian error 3.74811293e-13