From: Maximilian Bergbauer Date: Mon, 25 Sep 2023 08:20:57 +0000 (+0200) Subject: Implement integrate sum into values option X-Git-Tag: relicensing~431^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=bce0707290ef19f7f8795fffb9e36873170187bf;p=dealii.git Implement integrate sum into values option --- diff --git a/doc/news/changes/minor/20231004Bergbauer_2 b/doc/news/changes/minor/20231004Bergbauer_2 new file mode 100644 index 0000000000..a4e04403b5 --- /dev/null +++ b/doc/news/changes/minor/20231004Bergbauer_2 @@ -0,0 +1,4 @@ +Added: An option to sum into the solution values vector was added +to FEPointEvaluation integrate() and test_and_sum(). +
+(Maximilian Bergbauer, 2023/10/04) diff --git a/include/deal.II/matrix_free/fe_point_evaluation.h b/include/deal.II/matrix_free/fe_point_evaluation.h index c1e2c28216..15baefe1f0 100644 --- a/include/deal.II/matrix_free/fe_point_evaluation.h +++ b/include/deal.II/matrix_free/fe_point_evaluation.h @@ -885,12 +885,15 @@ public: * @param[in] integration_flags Flags specifying which quantities should be * integrated at the points. * + * @param[in] sum_into_values Flag specifying if the integrated values + * should be summed into the solution values. Defaults to false. + * */ template void integrate(const StridedArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &integration_flags); - + const EvaluationFlags::EvaluationFlags &integration_flags, + const bool sum_into_values = false); /** * This function multiplies the quantities passed in by previous @@ -909,15 +912,14 @@ public: * @param[in] integration_flags Flags specifying which quantities should be * integrated at the points. * + * @param[in] sum_into_values Flag specifying if the integrated values + * should be summed into the solution values. Defaults to false. + * */ void integrate(const ArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &integration_flags) - { - integrate(StridedArrayView(solution_values.data(), - solution_values.size()), - integration_flags); - } + const EvaluationFlags::EvaluationFlags &integration_flags, + const bool sum_into_values = false); /** * This function multiplies the quantities passed in by previous @@ -940,12 +942,16 @@ public: * @param[in] integration_flags Flags specifying which quantities should be * integrated at the points. * + * @param[in] sum_into_values Flag specifying if the integrated values + * should be summed into the solution values. Defaults to false. + * */ template void test_and_sum( const StridedArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &integration_flags); + const EvaluationFlags::EvaluationFlags &integration_flags, + const bool sum_into_values = false); /** * This function multiplies the quantities passed in by previous @@ -968,15 +974,14 @@ public: * @param[in] integration_flags Flags specifying which quantities should be * integrated at the points. * + * @param[in] sum_into_values Flag specifying if the integrated values + * should be summed into the solution values. Defaults to false. + * */ void test_and_sum(const ArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &integration_flags) - { - test_and_sum(StridedArrayView(solution_values.data(), - solution_values.size()), - integration_flags); - } + const EvaluationFlags::EvaluationFlags &integration_flags, + const bool sum_into_values = false); /** * Return the value at quadrature point number @p point_index after a call to @@ -1172,7 +1177,8 @@ private: finish_integrate_fast( const StridedArrayView &solution_values, const EvaluationFlags::EvaluationFlags &integration_flags, - vectorized_value_type *solution_values_vectorized_linear); + vectorized_value_type *solution_values_vectorized_linear, + const bool sum_into_values); /** * Fast path of the integrate function. @@ -1184,7 +1190,8 @@ private: void integrate_fast( const StridedArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &integration_flags); + const EvaluationFlags::EvaluationFlags &integration_flags, + const bool sum_into_values); /** * Slow path of the integrate function using FEValues. @@ -1193,7 +1200,8 @@ private: void integrate_slow( const StridedArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &integration_flags); + const EvaluationFlags::EvaluationFlags &integration_flags, + const bool sum_into_values); /** * Implementation of the integrate/test_and_sum function. @@ -1202,7 +1210,8 @@ private: void do_integrate( const StridedArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &integration_flags); + const EvaluationFlags::EvaluationFlags &integration_flags, + const bool sum_into_values); /** * Number of quadrature batches of the current cell/face. @@ -2393,9 +2402,10 @@ inline void FEPointEvaluation::finish_integrate_fast( const StridedArrayView &solution_values, const EvaluationFlags::EvaluationFlags &integration_flags, - vectorized_value_type *solution_values_vectorized_linear) + vectorized_value_type *solution_values_vectorized_linear, + const bool sum_into_values) { - if (fe->n_components() > n_components) + if (!sum_into_values && fe->n_components() > n_components) for (unsigned int i = 0; i < solution_values.size(); ++i) solution_values[i] = 0; @@ -2434,13 +2444,19 @@ FEPointEvaluation::finish_integrate_fast( if (is_linear || renumber.empty()) { for (unsigned int i = 0; i < dofs_per_comp; ++i) - solution_values[i + offset] = output[i]; + if (sum_into_values) + solution_values[i + offset] += output[i]; + else + solution_values[i + offset] = output[i]; } else { const unsigned int *renumber_ptr = renumber.data() + offset; for (unsigned int i = 0; i < dofs_per_comp; ++i) - solution_values[renumber_ptr[i]] = output[i]; + if (sum_into_values) + solution_values[renumber_ptr[i]] += output[i]; + else + solution_values[renumber_ptr[i]] = output[i]; } } else @@ -2448,18 +2464,29 @@ FEPointEvaluation::finish_integrate_fast( if (is_linear || renumber.empty()) { for (unsigned int i = 0; i < dofs_per_comp; ++i) - solution_values[i + offset] = - ETT::sum_value(comp, - is_linear ? - *(solution_values_vectorized_linear + i) : - solution_renumbered_vectorized[i]); + if (sum_into_values) + solution_values[i + offset] += + ETT::sum_value(comp, + is_linear ? + *(solution_values_vectorized_linear + i) : + solution_renumbered_vectorized[i]); + else + solution_values[i + offset] = + ETT::sum_value(comp, + is_linear ? + *(solution_values_vectorized_linear + i) : + solution_renumbered_vectorized[i]); } else { const unsigned int *renumber_ptr = renumber.data() + offset; for (unsigned int i = 0; i < dofs_per_comp; ++i) - solution_values[renumber_ptr[i]] = - ETT::sum_value(comp, solution_renumbered_vectorized[i]); + if (sum_into_values) + solution_values[renumber_ptr[i]] += + ETT::sum_value(comp, solution_renumbered_vectorized[i]); + else + solution_values[renumber_ptr[i]] = + ETT::sum_value(comp, solution_renumbered_vectorized[i]); } } } @@ -2475,7 +2502,8 @@ template ::integrate_fast( const StridedArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &integration_flags) + const EvaluationFlags::EvaluationFlags &integration_flags, + const bool sum_into_values) { // zero out lanes of incomplete last quadrature point batch if constexpr (stride == 1) @@ -2547,7 +2575,8 @@ FEPointEvaluation::integrate_fast( finish_integrate_fast( solution_values, integration_flags, - solution_values_vectorized_linear.data()); + solution_values_vectorized_linear.data(), + sum_into_values); } @@ -2557,14 +2586,16 @@ template inline void FEPointEvaluation::integrate_slow( const StridedArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &integration_flags) + const EvaluationFlags::EvaluationFlags &integration_flags, + const bool sum_into_values) { // slow path with FEValues Assert(fe_values.get() != nullptr, ExcMessage( "Not initialized. Please call FEPointEvaluation::reinit()!")); - for (unsigned int i = 0; i < solution_values.size(); ++i) - solution_values[i] = 0; + if (!sum_into_values) + for (unsigned int i = 0; i < solution_values.size(); ++i) + solution_values[i] = 0; const std::size_t n_points = fe_values->get_quadrature().size(); @@ -2634,7 +2665,8 @@ template void FEPointEvaluation::do_integrate( const StridedArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &integration_flags) + const EvaluationFlags::EvaluationFlags &integration_flags, + const bool sum_into_values) { if (!is_reinitialized) reinit(); @@ -2646,8 +2678,9 @@ FEPointEvaluation::do_integrate( (integration_flags & EvaluationFlags::gradients))) // no integration flags { - for (unsigned int i = 0; i < solution_values.size(); ++i) - solution_values[i] = 0; + if (!sum_into_values) + for (unsigned int i = 0; i < solution_values.size(); ++i) + solution_values[i] = 0; return; } @@ -2663,23 +2696,27 @@ FEPointEvaluation::do_integrate( { if (use_linear_path) integrate_fast(solution_values, - integration_flags); + integration_flags, + sum_into_values); else integrate_fast(solution_values, - integration_flags); + integration_flags, + sum_into_values); } else { if (use_linear_path) integrate_fast(solution_values, - integration_flags); + integration_flags, + sum_into_values); else integrate_fast(solution_values, - integration_flags); + integration_flags, + sum_into_values); } } else - integrate_slow(solution_values, integration_flags); + integrate_slow(solution_values, integration_flags, sum_into_values); } @@ -2689,9 +2726,25 @@ template void FEPointEvaluation::integrate( const StridedArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &integration_flags) + const EvaluationFlags::EvaluationFlags &integration_flags, + const bool sum_into_values) +{ + do_integrate(solution_values, integration_flags, sum_into_values); +} + + + +template +void +FEPointEvaluation::integrate( + const ArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &integration_flags, + const bool sum_into_values) { - do_integrate(solution_values, integration_flags); + integrate(StridedArrayView(solution_values.data(), + solution_values.size()), + integration_flags, + sum_into_values); } @@ -2701,9 +2754,25 @@ template void FEPointEvaluation::test_and_sum( const StridedArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &integration_flags) + const EvaluationFlags::EvaluationFlags &integration_flags, + const bool sum_into_values) +{ + do_integrate(solution_values, integration_flags, sum_into_values); +} + + + +template +void +FEPointEvaluation::test_and_sum( + const ArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &integration_flags, + const bool sum_into_values) { - do_integrate(solution_values, integration_flags); + test_and_sum(StridedArrayView(solution_values.data(), + solution_values.size()), + integration_flags, + sum_into_values); } diff --git a/tests/matrix_free/point_evaluation_28.cc b/tests/matrix_free/point_evaluation_28.cc new file mode 100644 index 0000000000..e91a8b3aa9 --- /dev/null +++ b/tests/matrix_free/point_evaluation_28.cc @@ -0,0 +1,142 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2020 - 2023 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 FEPointEvaluation integrate_add() + +#include + +#include + +#include +#include + +#include +#include + +#include + +#include + +#include + +#include + +#include "../tests.h" + + + +template +void +test(const unsigned int degree) +{ + using namespace dealii; + Triangulation tria; + + if (dim > 1) + GridGenerator::hyper_shell(tria, Point(), 0.5, 1, 6); + else + GridGenerator::subdivided_hyper_cube(tria, 2, 0, 1); + + MappingQ mapping(degree); + deallog << "Mapping of degree " << degree << std::endl; + + std::vector> unit_points; + for (unsigned int i = 0; i < Utilities::pow(2, dim); ++i) + { + Point p; + for (unsigned int d = 0; d < dim; ++d) + p[d] = static_cast(i) / 17. + 0.015625 * d; + unit_points.push_back(p); + } + + QGauss q_gauss(2); + + Quadrature quad(unit_points, q_gauss.get_weights()); + + FE_Q fe(degree); + + DoFHandler dof_handler(tria); + dof_handler.distribute_dofs(fe); + Vector vector(dof_handler.n_dofs()); + + NonMatching::MappingInfo mapping_info(mapping, + update_values | update_gradients | + update_JxW_values); + + std::vector> quad_vec_cell; + quad_vec_cell.reserve(tria.n_cells()); + + for (const auto &cell : tria.active_cell_iterators()) + quad_vec_cell.push_back(quad); + + mapping_info.reinit_cells(tria.active_cell_iterators(), quad_vec_cell); + + FEPointEvaluation<1, dim, dim, Number> evaluator(mapping_info, fe); + + Tensor<1, dim, Number> exponents; + exponents[0] = 1.; + VectorTools::interpolate(mapping, + dof_handler, + Functions::Monomial(exponents), + vector); + + std::vector solution_values(fe.dofs_per_cell); + + for (const auto &cell : dof_handler.active_cell_iterators()) + { + cell->get_dof_values(vector, + solution_values.begin(), + solution_values.end()); + + evaluator.reinit(cell->active_cell_index()); + evaluator.evaluate(solution_values, + EvaluationFlags::values | EvaluationFlags::gradients); + + for (const unsigned int i : evaluator.quadrature_point_indices()) + { + evaluator.submit_value(evaluator.get_value(i), i); + evaluator.submit_gradient(evaluator.get_gradient(i), i); + } + + evaluator.integrate(solution_values, + EvaluationFlags::values | EvaluationFlags::gradients); + + deallog << "Cell with center " << cell->center(true) << std::endl; + for (const auto i : solution_values) + deallog << i << ' '; + deallog << std::endl; + + evaluator.integrate(solution_values, + EvaluationFlags::values | EvaluationFlags::gradients, + true); + + for (const auto i : solution_values) + deallog << i << ' '; + deallog << std::endl; + } +} + + + +int +main() +{ + initlog(); + deallog << std::setprecision(10); + + test<2>(1); + test<2>(2); +} diff --git a/tests/matrix_free/point_evaluation_28.output b/tests/matrix_free/point_evaluation_28.output new file mode 100644 index 0000000000..ff41259658 --- /dev/null +++ b/tests/matrix_free/point_evaluation_28.output @@ -0,0 +1,39 @@ + +DEAL::Mapping of degree 1 +DEAL::Cell with center 0.6495190528 0.3750000000 +DEAL::1.094128055 0.06405871628 -0.7129529194 -0.07186179469 +DEAL::2.188256110 0.1281174326 -1.425905839 -0.1437235894 +DEAL::Cell with center 4.592425497e-17 0.7500000000 +DEAL::0.5256716676 -0.3778726922 0.05767485942 -0.04324345223 +DEAL::1.051343335 -0.7557453844 0.1153497188 -0.08648690447 +DEAL::Cell with center -0.6495190528 0.3750000000 +DEAL::-0.5684563874 -0.4419314085 0.7706277789 0.02861834246 +DEAL::-1.136912775 -0.8838628170 1.541255558 0.05723668491 +DEAL::Cell with center -0.6495190528 -0.3750000000 +DEAL::-1.094128055 -0.06405871628 0.7129529194 0.07186179469 +DEAL::-2.188256110 -0.1281174326 1.425905839 0.1437235894 +DEAL::Cell with center -1.377727649e-16 -0.7500000000 +DEAL::-0.5256716676 0.3778726922 -0.05767485942 0.04324345223 +DEAL::-1.051343335 0.7557453844 -0.1153497188 0.08648690447 +DEAL::Cell with center 0.6495190528 -0.3750000000 +DEAL::0.5684563874 0.4419314085 -0.7706277789 -0.02861834246 +DEAL::1.136912775 0.8838628170 -1.541255558 -0.05723668491 +DEAL::Mapping of degree 2 +DEAL::Cell with center 0.6495190528 0.3750000000 +DEAL::2.482631076 -0.1564725804 0.4910617067 -0.02558423728 -2.474880528 0.1748480233 0.7176813542 0.1215956751 -0.8420572937 +DEAL::4.965262152 -0.3129451607 0.9821234133 -0.05116847457 -4.949761056 0.3496960465 1.435362708 0.2431913503 -1.684114587 +DEAL::Cell with center 4.592425497e-17 0.7500000000 +DEAL::1.903653500 0.1627020517 0.1365654459 -0.02631613009 -0.6782922627 0.1422251499 -0.8006550044 0.1439386250 -0.7778032374 +DEAL::3.807306999 0.3254041034 0.2731308918 -0.05263226019 -1.356584525 0.2844502999 -1.601310009 0.2878772501 -1.555606475 +DEAL::Cell with center -0.6495190528 0.3750000000 +DEAL::-0.5789775763 0.3191746321 -0.3544962608 -0.0007318928105 1.796588265 -0.03262287333 -1.518336359 0.02234294989 0.06425405628 +DEAL::-1.157955153 0.6383492641 -0.7089925215 -0.001463785621 3.593176530 -0.06524574667 -3.036672717 0.04468589979 0.1285081126 +DEAL::Cell with center -0.6495190528 -0.3750000000 +DEAL::-2.482631076 0.1564725804 -0.4910617067 0.02558423728 2.474880528 -0.1748480233 -0.7176813542 -0.1215956751 0.8420572937 +DEAL::-4.965262152 0.3129451607 -0.9821234133 0.05116847457 4.949761056 -0.3496960465 -1.435362708 -0.2431913503 1.684114587 +DEAL::Cell with center -1.377727649e-16 -0.7500000000 +DEAL::-1.903653500 -0.1627020517 -0.1365654459 0.02631613009 0.6782922627 -0.1422251499 0.8006550044 -0.1439386250 0.7778032374 +DEAL::-3.807306999 -0.3254041034 -0.2731308918 0.05263226019 1.356584525 -0.2844502999 1.601310009 -0.2878772501 1.555606475 +DEAL::Cell with center 0.6495190528 -0.3750000000 +DEAL::0.5789775763 -0.3191746321 0.3544962608 0.0007318928105 -1.796588265 0.03262287333 1.518336359 -0.02234294989 -0.06425405628 +DEAL::1.157955153 -0.6383492641 0.7089925215 0.001463785621 -3.593176530 0.06524574667 3.036672717 -0.04468589979 -0.1285081126