New test cases

dealii · Apr 21, 2023 · ae3c6b0 · ae3c6b0
1 parent f94d31c
commit ae3c6b0
Show file tree

Hide file tree

Showing 4 changed files with 419 additions and 0 deletions.
diff --git a/tests/matrix_free/tensor_product_evaluate_07.cc b/tests/matrix_free/tensor_product_evaluate_07.cc
@@ -0,0 +1,114 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 - 2022 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 higher order derivatives with
+// evaluate_tensor_product_higher_derivatives for a scalar function on FE_Q
+
+#include <deal.II/base/function_lib.h>
+#include <deal.II/base/point.h>
+#include <deal.II/base/polynomial.h>
+#include <deal.II/base/quadrature_lib.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, int derivative_order>
+void
+test(const unsigned int degree)
+{
+  FE_Q<dim> fe(degree);
+
+  // go through all monomials of degree 'derivative_order + 1'
+  std::vector<std::array<int, 3>> exponents;
+  for (int i2 = 0; i2 <= (dim > 2 ? derivative_order : 0); ++i2)
+    for (int i1 = 0; i1 <= (dim > 1 ? derivative_order : 0); ++i1)
+      for (int i0 = 0; i0 <= derivative_order; ++i0)
+        if (i0 + i1 + i2 == derivative_order)
+          exponents.push_back(std::array<int, 3>{{i0, i1, i2}});
+
+  deallog << "Evaluate derivative " << derivative_order << " in " << dim
+          << "d with polynomial degree " << degree << std::endl;
+
+  const std::vector<unsigned int> renumbering =
+    FETools::lexicographic_to_hierarchic_numbering<dim>(degree);
+  const std::vector<Polynomials::Polynomial<double>> polynomials =
+    Polynomials::generate_complete_Lagrange_basis(
+      QGaussLobatto<1>(degree + 1).get_points());
+
+  Point<dim> p;
+  for (unsigned int d = 0; d < dim; ++d)
+    p[d] = 1.;
+
+  std::vector<double> function_values(fe.dofs_per_cell);
+
+  for (const std::array<int, 3> &exponent : exponents)
+    {
+      Tensor<1, dim> exp;
+      for (unsigned int d = 0; d < dim; ++d)
+        exp[d] = exponent[d];
+      Functions::Monomial<dim> monomial(exp);
+      for (unsigned int i = 0; i < fe.dofs_per_cell; ++i)
+        function_values[i] = monomial.value(fe.get_unit_support_points()[i]);
+
+      deallog << "Monomial [";
+      for (unsigned int d = 0; d < dim; ++d)
+        deallog << exponent[d] << (d == dim - 1 ? "" : " ");
+      deallog << "]: ";
+      const auto derivative =
+        internal::evaluate_tensor_product_higher_derivatives<derivative_order>(
+          polynomials, function_values, p, renumbering);
+
+      for (unsigned int d = 0; d < derivative.dimension; ++d)
+        deallog << (std::abs(derivative[d]) < 1e-11 ? 0. : derivative[d])
+                << " ";
+      deallog << std::endl;
+    }
+  deallog << std::endl;
+}
+
+
+
+int
+main()
+{
+  initlog();
+  deallog << std::setprecision(3);
+
+  // test 2nd derivatives
+  test<1, 2>(2);
+  test<2, 2>(2);
+  test<3, 2>(2);
+  test<3, 2>(3);
+
+  // test 3rd derivatives
+  test<1, 3>(3);
+  test<2, 3>(3);
+  test<3, 3>(3);
+  test<3, 3>(2);
+  test<3, 3>(4);
+
+  // test 4th derivatives
+  test<1, 4>(4);
+  test<2, 4>(4);
+  test<3, 4>(4);
+}
diff --git a/tests/matrix_free/tensor_product_evaluate_07.output b/tests/matrix_free/tensor_product_evaluate_07.output
@@ -0,0 +1,97 @@
+
+DEAL::Evaluate derivative 2 in 1d with polynomial degree 2
+DEAL::Monomial [2]: 2.00 
+DEAL::
+DEAL::Evaluate derivative 2 in 2d with polynomial degree 2
+DEAL::Monomial [2 0]: 2.00 0.00 0.00 
+DEAL::Monomial [1 1]: 0.00 1.00 0.00 
+DEAL::Monomial [0 2]: 0.00 0.00 2.00 
+DEAL::
+DEAL::Evaluate derivative 2 in 3d with polynomial degree 2
+DEAL::Monomial [2 0 0]: 2.00 0.00 0.00 0.00 0.00 0.00 
+DEAL::Monomial [1 1 0]: 0.00 1.00 0.00 0.00 0.00 0.00 
+DEAL::Monomial [0 2 0]: 0.00 0.00 2.00 0.00 0.00 0.00 
+DEAL::Monomial [1 0 1]: 0.00 0.00 0.00 1.00 0.00 0.00 
+DEAL::Monomial [0 1 1]: 0.00 0.00 0.00 0.00 1.00 0.00 
+DEAL::Monomial [0 0 2]: 0.00 0.00 0.00 0.00 0.00 2.00 
+DEAL::
+DEAL::Evaluate derivative 2 in 3d with polynomial degree 3
+DEAL::Monomial [2 0 0]: 2.00 0.00 0.00 0.00 0.00 0.00 
+DEAL::Monomial [1 1 0]: 0.00 1.00 0.00 0.00 0.00 0.00 
+DEAL::Monomial [0 2 0]: 0.00 0.00 2.00 0.00 0.00 0.00 
+DEAL::Monomial [1 0 1]: 0.00 0.00 0.00 1.00 0.00 0.00 
+DEAL::Monomial [0 1 1]: 0.00 0.00 0.00 0.00 1.00 0.00 
+DEAL::Monomial [0 0 2]: 0.00 0.00 0.00 0.00 0.00 2.00 
+DEAL::
+DEAL::Evaluate derivative 3 in 1d with polynomial degree 3
+DEAL::Monomial [3]: 6.00 
+DEAL::
+DEAL::Evaluate derivative 3 in 2d with polynomial degree 3
+DEAL::Monomial [3 0]: 6.00 0.00 0.00 0.00 
+DEAL::Monomial [2 1]: 0.00 2.00 0.00 0.00 
+DEAL::Monomial [1 2]: 0.00 0.00 2.00 0.00 
+DEAL::Monomial [0 3]: 0.00 0.00 0.00 6.00 
+DEAL::
+DEAL::Evaluate derivative 3 in 3d with polynomial degree 3
+DEAL::Monomial [3 0 0]: 6.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 
+DEAL::Monomial [2 1 0]: 0.00 2.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 
+DEAL::Monomial [1 2 0]: 0.00 0.00 2.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 
+DEAL::Monomial [0 3 0]: 0.00 0.00 0.00 6.00 0.00 0.00 0.00 0.00 0.00 0.00 
+DEAL::Monomial [2 0 1]: 0.00 0.00 0.00 0.00 2.00 0.00 0.00 0.00 0.00 0.00 
+DEAL::Monomial [1 1 1]: 0.00 0.00 0.00 0.00 0.00 1.00 0.00 0.00 0.00 0.00 
+DEAL::Monomial [0 2 1]: 0.00 0.00 0.00 0.00 0.00 0.00 2.00 0.00 0.00 0.00 
+DEAL::Monomial [1 0 2]: 0.00 0.00 0.00 0.00 0.00 0.00 0.00 2.00 0.00 0.00 
+DEAL::Monomial [0 1 2]: 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 2.00 0.00 
+DEAL::Monomial [0 0 3]: 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 6.00 
+DEAL::
+DEAL::Evaluate derivative 3 in 3d with polynomial degree 2
+DEAL::Monomial [3 0 0]: 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 
+DEAL::Monomial [2 1 0]: 0.00 2.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 
+DEAL::Monomial [1 2 0]: 0.00 0.00 2.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 
+DEAL::Monomial [0 3 0]: 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 
+DEAL::Monomial [2 0 1]: 0.00 0.00 0.00 0.00 2.00 0.00 0.00 0.00 0.00 0.00 
+DEAL::Monomial [1 1 1]: 0.00 0.00 0.00 0.00 0.00 1.00 0.00 0.00 0.00 0.00 
+DEAL::Monomial [0 2 1]: 0.00 0.00 0.00 0.00 0.00 0.00 2.00 0.00 0.00 0.00 
+DEAL::Monomial [1 0 2]: 0.00 0.00 0.00 0.00 0.00 0.00 0.00 2.00 0.00 0.00 
+DEAL::Monomial [0 1 2]: 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 2.00 0.00 
+DEAL::Monomial [0 0 3]: 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 
+DEAL::
+DEAL::Evaluate derivative 3 in 3d with polynomial degree 4
+DEAL::Monomial [3 0 0]: 6.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 
+DEAL::Monomial [2 1 0]: 0.00 2.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 
+DEAL::Monomial [1 2 0]: 0.00 0.00 2.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 
+DEAL::Monomial [0 3 0]: 0.00 0.00 0.00 6.00 0.00 0.00 0.00 0.00 0.00 0.00 
+DEAL::Monomial [2 0 1]: 0.00 0.00 0.00 0.00 2.00 0.00 0.00 0.00 0.00 0.00 
+DEAL::Monomial [1 1 1]: 0.00 0.00 0.00 0.00 0.00 1.00 0.00 0.00 0.00 0.00 
+DEAL::Monomial [0 2 1]: 0.00 0.00 0.00 0.00 0.00 0.00 2.00 0.00 0.00 0.00 
+DEAL::Monomial [1 0 2]: 0.00 0.00 0.00 0.00 0.00 0.00 0.00 2.00 0.00 0.00 
+DEAL::Monomial [0 1 2]: 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 2.00 0.00 
+DEAL::Monomial [0 0 3]: 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 6.00 
+DEAL::
+DEAL::Evaluate derivative 4 in 1d with polynomial degree 4
+DEAL::Monomial [4]: 24.0 
+DEAL::
+DEAL::Evaluate derivative 4 in 2d with polynomial degree 4
+DEAL::Monomial [4 0]: 24.0 0.00 0.00 0.00 0.00 
+DEAL::Monomial [3 1]: 0.00 6.00 0.00 0.00 0.00 
+DEAL::Monomial [2 2]: 0.00 0.00 4.00 0.00 0.00 
+DEAL::Monomial [1 3]: 0.00 0.00 0.00 6.00 0.00 
+DEAL::Monomial [0 4]: 0.00 0.00 0.00 0.00 24.0 
+DEAL::
+DEAL::Evaluate derivative 4 in 3d with polynomial degree 4
+DEAL::Monomial [4 0 0]: 24.0 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 
+DEAL::Monomial [3 1 0]: 0.00 6.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 
+DEAL::Monomial [2 2 0]: 0.00 0.00 4.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 
+DEAL::Monomial [1 3 0]: 0.00 0.00 0.00 6.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 
+DEAL::Monomial [0 4 0]: 0.00 0.00 0.00 0.00 24.0 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 
+DEAL::Monomial [3 0 1]: 0.00 0.00 0.00 0.00 0.00 6.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 
+DEAL::Monomial [2 1 1]: 0.00 0.00 0.00 0.00 0.00 0.00 2.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 
+DEAL::Monomial [1 2 1]: 0.00 0.00 0.00 0.00 0.00 0.00 0.00 2.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 
+DEAL::Monomial [0 3 1]: 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 6.00 0.00 0.00 0.00 0.00 0.00 0.00 
+DEAL::Monomial [2 0 2]: 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 4.00 0.00 0.00 0.00 0.00 0.00 
+DEAL::Monomial [1 1 2]: 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 2.00 0.00 0.00 0.00 0.00 
+DEAL::Monomial [0 2 2]: 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 4.00 0.00 0.00 0.00 
+DEAL::Monomial [1 0 3]: 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 6.00 0.00 0.00 
+DEAL::Monomial [0 1 3]: 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 6.00 0.00 
+DEAL::Monomial [0 0 4]: 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 24.0 
+DEAL::
diff --git a/tests/matrix_free/tensor_product_evaluate_08.cc b/tests/matrix_free/tensor_product_evaluate_08.cc
@@ -0,0 +1,111 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 - 2022 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 higher order derivatives with
+// evaluate_tensor_product_higher_derivatives for a scalar function on FE_DGQ
+
+#include <deal.II/base/function_lib.h>
+#include <deal.II/base/point.h>
+#include <deal.II/base/polynomial.h>
+#include <deal.II/base/quadrature_lib.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, int derivative_order>
+void
+test(const unsigned int degree)
+{
+  FE_DGQ<dim> fe(degree);
+
+  // go through all monomials of degree 'derivative_order + 1'
+  std::vector<std::array<int, 3>> exponents;
+  for (int i2 = 0; i2 <= (dim > 2 ? derivative_order : 0); ++i2)
+    for (int i1 = 0; i1 <= (dim > 1 ? derivative_order : 0); ++i1)
+      for (int i0 = 0; i0 <= derivative_order; ++i0)
+        if (i0 + i1 + i2 == derivative_order)
+          exponents.push_back(std::array<int, 3>{{i0, i1, i2}});
+
+  deallog << "Evaluate derivative " << derivative_order << " in " << dim
+          << "d with polynomial degree " << degree << std::endl;
+
+  const std::vector<Polynomials::Polynomial<double>> polynomials =
+    Polynomials::generate_complete_Lagrange_basis(
+      QGaussLobatto<1>(degree + 1).get_points());
+
+  Point<dim> p;
+  for (unsigned int d = 0; d < dim; ++d)
+    p[d] = 1.;
+
+  std::vector<double> function_values(fe.dofs_per_cell);
+
+  for (const std::array<int, 3> &exponent : exponents)
+    {
+      Tensor<1, dim> exp;
+      for (unsigned int d = 0; d < dim; ++d)
+        exp[d] = exponent[d];
+      Functions::Monomial<dim> monomial(exp);
+      for (unsigned int i = 0; i < fe.dofs_per_cell; ++i)
+        function_values[i] = monomial.value(fe.get_unit_support_points()[i]);
+
+      deallog << "Monomial [";
+      for (unsigned int d = 0; d < dim; ++d)
+        deallog << exponent[d] << (d == dim - 1 ? "" : " ");
+      deallog << "]: ";
+      const auto derivative =
+        internal::evaluate_tensor_product_higher_derivatives<derivative_order>(
+          polynomials, function_values, p);
+
+      for (unsigned int d = 0; d < derivative.dimension; ++d)
+        deallog << (std::abs(derivative[d]) < 1e-11 ? 0. : derivative[d])
+                << " ";
+      deallog << std::endl;
+    }
+  deallog << std::endl;
+}
+
+
+
+int
+main()
+{
+  initlog();
+  deallog << std::setprecision(3);
+
+  // test 2nd derivatives
+  test<1, 2>(2);
+  test<2, 2>(2);
+  test<3, 2>(2);
+  test<3, 2>(3);
+
+  // test 3rd derivatives
+  test<1, 3>(3);
+  test<2, 3>(3);
+  test<3, 3>(3);
+  test<3, 3>(2);
+  test<3, 3>(4);
+
+  // test 4th derivatives
+  test<1, 4>(4);
+  test<2, 4>(4);
+  test<3, 4>(4);
+}