From cbbc24e71fff1ad659bf3669ee0887f50a48d042 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Thu, 21 Mar 2024 10:15:43 +0100 Subject: [PATCH] Introduce FEPointEvaluation::integrate_value() --- include/deal.II/matrix_free/fe_evaluation.h | 5 +- .../deal.II/matrix_free/fe_point_evaluation.h | 84 +++++ tests/matrix_free/integrate_value_01.cc | 313 ++++++++++++++++++ tests/matrix_free/integrate_value_01.output | 8 + 4 files changed, 408 insertions(+), 2 deletions(-) create mode 100644 tests/matrix_free/integrate_value_01.cc create mode 100644 tests/matrix_free/integrate_value_01.output diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h index a8247d1a44..882d8839b1 100644 --- a/include/deal.II/matrix_free/fe_evaluation.h +++ b/include/deal.II/matrix_free/fe_evaluation.h @@ -558,8 +558,9 @@ public: #endif /** - * Takes values at quadrature points, multiplies by the Jacobian determinant - * and quadrature weights (JxW) and sums the values for all quadrature + * Take values collected at quadrature points via the submit_value() function, + * multiply by the Jacobian determinant + * and quadrature weights (JxW) and sum the values for all quadrature * points on the cell. The result is a scalar, representing the integral * over the function over the cell. If a vector-element is used, the * resulting components are still separated. Moreover, if vectorization is diff --git a/include/deal.II/matrix_free/fe_point_evaluation.h b/include/deal.II/matrix_free/fe_point_evaluation.h index e7f1426163..cae51b5b7c 100644 --- a/include/deal.II/matrix_free/fe_point_evaluation.h +++ b/include/deal.II/matrix_free/fe_point_evaluation.h @@ -85,6 +85,23 @@ namespace internal result[component] = vector_entry; } + static scalar_value_type + sum_value(const scalar_value_type &result) + { + return result; + } + + static scalar_value_type + sum_value(const vectorized_value_type &result) + { + scalar_value_type result_scalar = {}; + + for (unsigned int c = 0; c < n_components; ++c) + result_scalar[c] = result[c].sum(); + + return result_scalar; + } + static ScalarNumber sum_value(const unsigned int component, const vectorized_value_type &result) @@ -236,6 +253,18 @@ namespace internal result = vector_entry; } + static scalar_value_type + sum_value(const scalar_value_type &result) + { + return result; + } + + static scalar_value_type + sum_value(const vectorized_value_type &result) + { + return result.sum(); + } + static ScalarNumber sum_value(const unsigned int, const vectorized_value_type &result) { @@ -390,6 +419,23 @@ namespace internal result[component] = vector_entry; } + static scalar_value_type + sum_value(const scalar_value_type &result) + { + return result; + } + + static scalar_value_type + sum_value(const vectorized_value_type &result) + { + scalar_value_type result_scalar = {}; + + for (unsigned int c = 0; c < dim; ++c) + result_scalar[c] = result[c].sum(); + + return result_scalar; + } + static ScalarNumber sum_value(const unsigned int component, const vectorized_value_type &result) @@ -541,6 +587,18 @@ namespace internal result = vector_entry; } + static scalar_value_type + sum_value(const scalar_value_type &result) + { + return result; + } + + static scalar_value_type + sum_value(const vectorized_value_type &result) + { + return result.sum(); + } + static ScalarNumber sum_value(const unsigned int, const vectorized_value_type &result) { @@ -859,6 +917,16 @@ public: Point unit_point(const unsigned int point_index) const; + /** + * Take values collected at quadrature points via the submit_value() + * function, multiply by the Jacobian determinant + * and quadrature weights (JxW) and sum the values for all quadrature + * points on the cell. The result is a scalar, representing the integral + * of the function over the cell. + */ + scalar_value_type + integrate_value() const; + /** * Return an object that can be thought of as an array containing all indices * from zero to n_quadrature_points. This allows to write code using @@ -2466,6 +2534,22 @@ FEPointEvaluation::integrate( +template +inline typename FEPointEvaluationBase:: + scalar_value_type + FEPointEvaluationBase::integrate_value() + const +{ + value_type return_value = {}; + + for (const auto point_index : this->quadrature_point_indices()) + return_value += values[point_index] * this->JxW(point_index); + + return ETT::sum_value(return_value); +} + + + template template void diff --git a/tests/matrix_free/integrate_value_01.cc b/tests/matrix_free/integrate_value_01.cc new file mode 100644 index 0000000000..0dfee3258c --- /dev/null +++ b/tests/matrix_free/integrate_value_01.cc @@ -0,0 +1,313 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2024 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 FEEvaluation::integrate_value() and +// FEPointEvaluation::integrate_value(). + +#include + +#include +#include + +#include +#include +#include + +#include +#include +#include + +#include +#include + +#include +#include + +#include + +#include + +#include "../tests.h" + +template +class RightHandSideFunction : public Function +{ +public: + RightHandSideFunction(const unsigned int n_components) + : Function(n_components) + {} + + Number + value(const Point &p, const unsigned int component) const + { + if (component == 0) + return p[0]; + else if (component == 1) + return p[1] * p[1]; + else + return 0; + } + + VectorizedArray + value(const Point> &p, + const unsigned int component) const + { + if (component == 0) + return p[0]; + else if (component == 1) + return p[1] * p[1]; + else + return 0; + } +}; + + +template +Number +sum(const Number &value) +{ + return value; +} + + +template +Number +sum(const Tensor<1, dim, Number> &value) +{ + Number result = {}; + + for (unsigned int i = 0; i < dim; ++i) + result += value[i]; + + return result; +} + + +template +Number & +get(Number &value, const unsigned int) +{ + return value; +} + + +template +Number & +get(Tensor<1, dim, Number> &value, const unsigned int component) +{ + return value[component]; +} + + + +template +void +test() +{ + const unsigned int fe_degree = 2; + const unsigned int n_global_mesh_refinements = 2; + + Triangulation tria; + GridGenerator::hyper_cube(tria); + tria.refine_global(n_global_mesh_refinements); + + FESystem fe(FE_Q(fe_degree), n_components); + QGauss quad(fe_degree + 1); + + MappingQ1 mapping; + + DoFHandler dof_handler(tria); + dof_handler.distribute_dofs(fe); + + Vector vector(dof_handler.n_dofs()); + + RightHandSideFunction fu(n_components); + + Vector integrals(tria.n_active_cells()); + + // create reference solution with VectorTools::integrate_difference() + { + Functions::ConstantFunction zero(0.0, n_components); + VectorTools::integrate_difference(mapping, + dof_handler, + vector, + fu, + integrals, + quad, + VectorTools::NormType::L2_norm); + + for (unsigned int i = 0; i < integrals.size(); ++i) + integrals[i] = integrals[i] * integrals[i]; + } + + // create solution with FEEvaluation::integrate_value() + { + Vector my_integrals(tria.n_active_cells()); + + typename MatrixFree::AdditionalData ad; + ad.mapping_update_flags = + update_JxW_values | update_values | update_quadrature_points; + MatrixFree matrix_free; + + matrix_free.reinit( + mapping, dof_handler, AffineConstraints(), quad, ad); + + FEEvaluation fe_eval(matrix_free); + + for (unsigned int cell = 0; cell < matrix_free.n_cell_batches(); ++cell) + { + fe_eval.reinit(cell); + + for (const auto q : fe_eval.quadrature_point_indices()) + { + const auto point = fe_eval.quadrature_point(q); + + typename FEEvaluation::value_type value; + + for (unsigned int c = 0; c < n_components; ++c) + get(value, c) = fu.value(point, c) * fu.value(point, c); + + fe_eval.submit_value(value, q); + } + + const auto integrated_cell_value = sum(fe_eval.integrate_value()); + + for (unsigned int v = 0; + v < matrix_free.n_active_entries_per_cell_batch(cell); + ++v) + my_integrals[matrix_free.get_cell_iterator(cell, v) + ->active_cell_index()] = integrated_cell_value[v]; + } + + for (unsigned int i = 0; i < my_integrals.size(); ++i) + Assert(std::abs(my_integrals[i] - integrals[i]) < 1e-12, + ExcInternalError()); + } + + // create solution with FEEvaluation::integrate_value() (non-vectorized) + { + Vector my_integrals(tria.n_active_cells()); + + NonMatching::MappingInfo mapping_info(mapping, + update_JxW_values | + update_values); + + std::vector::cell_iterator> cells; + std::vector> quadratures; + + for (const auto &cell : tria.active_cell_iterators()) + { + cells.push_back(cell); + quadratures.push_back(quad); + } + + mapping_info.reinit_cells(cells, quadratures); + + FEPointEvaluation fe_eval(mapping_info, fe); + + for (unsigned int cell = 0; cell < cells.size(); ++cell) + { + fe_eval.reinit(cell); + + for (const auto q : fe_eval.quadrature_point_indices()) + { + const auto point = fe_eval.real_point(q); + + typename FEPointEvaluation:: + value_type value; + + for (unsigned int c = 0; c < n_components; ++c) + get(value, c) = fu.value(point, c) * fu.value(point, c); + + fe_eval.submit_value(value, q); + } + + my_integrals[cell] = sum(fe_eval.integrate_value()); + } + + for (unsigned int i = 0; i < my_integrals.size(); ++i) + Assert(std::abs(my_integrals[i] - integrals[i]) < 1e-12, + ExcInternalError()); + } + + // create solution with FEEvaluation::integrate_value() (vectorized) + { + Vector my_integrals(tria.n_active_cells()); + + NonMatching::MappingInfo> mapping_info( + mapping, update_JxW_values | update_values); + + std::vector::cell_iterator> cells; + std::vector> quadratures; + + for (const auto &cell : tria.active_cell_iterators()) + { + cells.push_back(cell); + quadratures.push_back(quad); + } + + mapping_info.reinit_cells(cells, quadratures); + + FEPointEvaluation> fe_eval( + mapping_info, fe); + + for (unsigned int cell = 0; cell < cells.size(); ++cell) + { + fe_eval.reinit(cell); + + for (const auto q : fe_eval.quadrature_point_indices()) + { + const auto point = fe_eval.real_point(q); + + typename FEPointEvaluation>::value_type + value; + + for (unsigned int c = 0; c < n_components; ++c) + get(value, c) = fu.value(point, c) * fu.value(point, c); + + fe_eval.submit_value(value, q); + } + + my_integrals[cell] = sum(fe_eval.integrate_value()); + } + + for (unsigned int i = 0; i < my_integrals.size(); ++i) + Assert(std::abs(my_integrals[i] - integrals[i]) < 1e-12, + ExcInternalError()); + } + + deallog << "OK" << std::endl; +} + + +int +main() +{ + initlog(); + + test<2, 1>(); + test<2, 2>(); + test<2, 3>(); + + test<3, 1>(); + test<3, 2>(); + test<3, 3>(); + test<3, 4>(); +} diff --git a/tests/matrix_free/integrate_value_01.output b/tests/matrix_free/integrate_value_01.output new file mode 100644 index 0000000000..5cf798a265 --- /dev/null +++ b/tests/matrix_free/integrate_value_01.output @@ -0,0 +1,8 @@ + +DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK -- 2.39.5