From 2ccfde68d449c8583375f3a9c54b328bc4c7ae6d Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Fri, 28 Jan 2022 22:26:17 +0100 Subject: [PATCH] Add test case --- .../matrix_free/tensor_product_evaluate_05.cc | 102 ++++++++++++++ .../tensor_product_evaluate_05.output | 125 ++++++++++++++++++ .../matrix_free/tensor_product_evaluate_06.cc | 99 ++++++++++++++ .../tensor_product_evaluate_06.output | 125 ++++++++++++++++++ 4 files changed, 451 insertions(+) create mode 100644 tests/matrix_free/tensor_product_evaluate_05.cc create mode 100644 tests/matrix_free/tensor_product_evaluate_05.output create mode 100644 tests/matrix_free/tensor_product_evaluate_06.cc create mode 100644 tests/matrix_free/tensor_product_evaluate_06.output diff --git a/tests/matrix_free/tensor_product_evaluate_05.cc b/tests/matrix_free/tensor_product_evaluate_05.cc new file mode 100644 index 0000000000..745013c9a5 --- /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 +#include +#include +#include +#include + +#include +#include + +#include + +#include + +#include "../tests.h" + +template +void +test(const unsigned int degree) +{ + FE_Q fe(degree); + + // choose a symmetric matrix and then construct f(x) = x^T A x + SymmetricTensor<2, dim> matrix = unit_symmetric_tensor(); + 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 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 = + Polynomials::generate_complete_Lagrange_basis( + QGaussLobatto<1>(degree + 1).get_points()); + + const std::vector renumbering = + FETools::lexicographic_to_hierarchic_numbering(degree); + + const std::vector> evaluation_points = + dim == 3 ? QGauss(2).get_points() : + QIterated(QTrapez<1>(), 3).get_points(); + + deallog << "Evaluate in " << dim << "d with polynomial degree " << degree + << std::endl; + for (const auto &p : evaluation_points) + { + Point> p_vec; + for (unsigned int v = 0; v < VectorizedArray::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::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 index 0000000000..731cec2af8 --- /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 index 0000000000..883736fe75 --- /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 +#include +#include +#include +#include + +#include + +#include + +#include + +#include "../tests.h" + +template +void +test(const unsigned int degree) +{ + FE_DGQ fe(degree); + + // choose a symmetric matrix and then construct f(x) = x^T A x + SymmetricTensor<2, dim> matrix = unit_symmetric_tensor(); + 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 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 = + Polynomials::generate_complete_Lagrange_basis( + QGaussLobatto<1>(degree + 1).get_points()); + + const std::vector> evaluation_points = + dim == 3 ? QGauss(2).get_points() : + QIterated(QTrapez<1>(), 3).get_points(); + + deallog << "Evaluate in " << dim << "d with polynomial degree " << degree + << std::endl; + for (const auto &p : evaluation_points) + { + Point> p_vec; + for (unsigned int v = 0; v < VectorizedArray::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::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 index 0000000000..731cec2af8 --- /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 -- 2.39.5