From d928af7b50910311636a647fe0708b5257478167 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Tue, 26 Jan 2021 12:20:09 +0100 Subject: [PATCH] Update interface of FEPointEvaluation Conflicts: .github/workflows/linux.yml --- include/deal.II/fe/fe_point_evaluation.h | 118 ++++++++++++++++------- tests/fe/point_evaluation_01.cc | 15 +-- tests/fe/point_evaluation_02.cc | 15 +-- tests/fe/point_evaluation_03.cc | 16 ++- tests/fe/point_evaluation_04.cc | 26 ++--- tests/fe/point_evaluation_05.cc | 15 +-- tests/fe/point_evaluation_06.cc | 15 +-- tests/fe/point_evaluation_07.cc | 16 ++- tests/fe/point_evaluation_08.cc | 26 ++--- tests/fe/point_evaluation_09.cc | 17 ++-- 10 files changed, 141 insertions(+), 138 deletions(-) diff --git a/include/deal.II/fe/fe_point_evaluation.h b/include/deal.II/fe/fe_point_evaluation.h index 2a637cd27c..eeeae8efcc 100644 --- a/include/deal.II/fe/fe_point_evaluation.h +++ b/include/deal.II/fe/fe_point_evaluation.h @@ -27,6 +27,7 @@ #include +#include #include #include @@ -271,23 +272,36 @@ public: * values on the element as returned by `cell->get_dof_values(global_vector, * solution_values)`. * - * @param[out] values Array where the values of the evaluation of the - * solution interpolated to the given points is stored. The array must have - * the same size as the `unit_points` argument if values should be - * computed. If empty, values will not be computed. - * - * @param[out] gradients Array where the gradients of the evaluation of the - * solution interpolated to the given points is stored. The gradients are - * expressed in real space via the Mapping defined in the constructor. The - * array must have the same size as the `unit_points` argument if gradients - * should be computed. If empty, gradients will not be computed. + * @param[in] flags Flags specifying which quantities should be evaluated + * at the points. */ void evaluate(const typename Triangulation::cell_iterator &cell, - const ArrayView> &unit_points, - const ArrayView & solution_values, - const ArrayView & values, - const ArrayView & gradients = {}); + const ArrayView> & unit_points, + const ArrayView & solution_values, + const EvaluationFlags::EvaluationFlags &flags); + + /** + * Return the value at quadrature point number @p q_point after a call to + * FEPointEvaluation::evaluate() with EvaluationFlags::value set, or + * the value that has been stored there with a call to + * FEPointEvaluation::submit_value(). If the object is vector-valued, a + * vector-valued return argument is given. Note that when + * vectorization is enabled, values from several points are grouped together. + */ + const value_type & + get_value(const unsigned int q_point) const; + + /** + * Return the gradient at quadrature point number @p q_point after a call to + * FEPointEvaluation::evaluate() with EvaluationFlags::gradient set, or + * the gradient that has been stored there with a call to + * FEPointEvaluation::submit_gradient(). If the object is vector-valued, a + * vector-valued return argument is given. Note that when + * vectorization is enabled, values from several points are grouped together. + */ + const gradient_type & + get_gradient(const unsigned int q_point) const; private: /** @@ -326,6 +340,16 @@ private: */ std::vector solution_renumbered; + /** + * Temporary array to store the values at the points. + */ + std::vector values; + + /** + * Temporary array to store the gradients at the points. + */ + std::vector gradients; + /** * Number of unknowns per component, i.e., number of unique basis functions, * for the chosen FiniteElement (or base element). @@ -406,11 +430,12 @@ FEPointEvaluation::evaluate( const typename Triangulation::cell_iterator &cell, const ArrayView> & unit_points, const ArrayView & solution_values, - const ArrayView & values, - const ArrayView & gradients) + const EvaluationFlags::EvaluationFlags & evaluation_flag) { AssertDimension(solution_values.size(), fe->dofs_per_cell); - if ((values.size() > 0 || gradients.size() > 0) && !poly.empty()) + if (((evaluation_flag & EvaluationFlags::values) || + (evaluation_flag & EvaluationFlags::gradients)) && + !poly.empty()) { // fast path with tensor product evaluation if (solution_renumbered.size() != dofs_per_component) @@ -422,10 +447,10 @@ FEPointEvaluation::evaluate( comp, solution_renumbered[i]); - if (values.size() > 0) - AssertDimension(unit_points.size(), values.size()); - if (gradients.size() > 0) - AssertDimension(unit_points.size(), gradients.size()); + if (evaluation_flag & EvaluationFlags::values) + values.resize(unit_points.size()); + if (evaluation_flag & EvaluationFlags::gradients) + gradients.resize(unit_points.size()); const std::size_t n_points = unit_points.size(); const std::size_t n_lanes = VectorizedArray::size(); @@ -443,12 +468,12 @@ FEPointEvaluation::evaluate( poly, solution_renumbered, vectorized_points, poly.size() == 2); // convert back to standard format - if (values.size() > 0) + if (evaluation_flag & EvaluationFlags::values) for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j) internal::FEPointEvaluation::EvaluatorTypeTraits< dim, n_components>::set_value(val_and_grad.first, j, values[i + j]); - if (gradients.size() > 0) + if (evaluation_flag & EvaluationFlags::gradients) for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j) internal::FEPointEvaluation:: EvaluatorTypeTraits::set_gradient( @@ -456,23 +481,26 @@ FEPointEvaluation::evaluate( } // let mapping compute the transformation - if (gradients.size() > 0) + if (evaluation_flag & EvaluationFlags::gradients) { Assert(mapping_q_generic != nullptr, ExcInternalError()); mapping_q_generic->transform_variable( cell, mapping_covariant, unit_points, - ArrayView(gradients.begin(), gradients.size()), - gradients); + ArrayView(gradients.data(), gradients.size()), + ArrayView(gradients.data(), gradients.size())); } } - else if (values.size() > 0 || gradients.size() > 0) + else if ((evaluation_flag & EvaluationFlags::values) || + (evaluation_flag & EvaluationFlags::gradients)) { // slow path with FEValues const UpdateFlags flags = - (values.size() > 0 ? update_values : update_default) | - (gradients.size() > 0 ? update_gradients : update_default); + ((evaluation_flag & EvaluationFlags::values) ? update_values : + update_default) | + ((evaluation_flag & EvaluationFlags::gradients) ? update_gradients : + update_default); FEValues fe_values( *mapping, *fe, @@ -481,9 +509,9 @@ FEPointEvaluation::evaluate( flags); fe_values.reinit(cell); - if (values.size() > 0) + if (evaluation_flag & EvaluationFlags::values) { - AssertDimension(unit_points.size(), values.size()); + values.resize(unit_points.size()); std::fill(values.begin(), values.end(), value_type()); for (unsigned int i = 0; i < fe->n_dofs_per_cell(); ++i) { @@ -505,9 +533,9 @@ FEPointEvaluation::evaluate( } } - if (gradients.size() > 0) + if (evaluation_flag & EvaluationFlags::gradients) { - AssertDimension(unit_points.size(), gradients.size()); + gradients.resize(unit_points.size()); std::fill(gradients.begin(), gradients.end(), gradient_type()); for (unsigned int i = 0; i < fe->n_dofs_per_cell(); ++i) { @@ -530,6 +558,30 @@ FEPointEvaluation::evaluate( } } + + +template +inline const typename FEPointEvaluation::value_type + & + FEPointEvaluation::get_value( + const unsigned int q_point) const +{ + AssertIndexRange(q_point, values.size()); + return values[q_point]; +} + + + +template +inline const typename FEPointEvaluation:: + gradient_type & + FEPointEvaluation::get_gradient( + const unsigned int q_point) const +{ + AssertIndexRange(q_point, gradients.size()); + return gradients[q_point]; +} + DEAL_II_NAMESPACE_CLOSE #endif diff --git a/tests/fe/point_evaluation_01.cc b/tests/fe/point_evaluation_01.cc index 3fe84a69eb..92afd4347a 100644 --- a/tests/fe/point_evaluation_01.cc +++ b/tests/fe/point_evaluation_01.cc @@ -86,10 +86,6 @@ test(const unsigned int degree) std::vector function_values(unit_points.size()); std::vector> function_gradients(unit_points.size()); - std::vector dof_indices(fe.dofs_per_cell); - std::vector function_values_1(unit_points.size()); - std::vector> function_gradients_1(unit_points.size()); - for (const auto &cell : dof_handler.active_cell_iterators()) { fe_values.reinit(cell); @@ -103,15 +99,14 @@ test(const unsigned int degree) evaluator.evaluate(cell, unit_points, solution_values, - function_values_1, - function_gradients_1); + EvaluationFlags::values | EvaluationFlags::gradients); deallog << "Cell with center " << cell->center(true) << std::endl; - for (unsigned int i = 0; i < function_values_1.size(); ++i) + for (unsigned int i = 0; i < function_values.size(); ++i) deallog << mapping.transform_unit_to_real_cell(cell, unit_points[i]) - << ": " << function_values_1[i] << " error value " - << function_values[i] - function_values_1[i] << " error grad " - << (function_gradients_1[i] - function_gradients[i]).norm() + << ": " << evaluator.get_value(i) << " error value " + << function_values[i] - evaluator.get_value(i) << " error grad " + << (evaluator.get_gradient(i) - function_gradients[i]).norm() << std::endl; deallog << std::endl; } diff --git a/tests/fe/point_evaluation_02.cc b/tests/fe/point_evaluation_02.cc index 678983ff55..a6e8c07ba6 100644 --- a/tests/fe/point_evaluation_02.cc +++ b/tests/fe/point_evaluation_02.cc @@ -86,10 +86,6 @@ test(const unsigned int degree) std::vector function_values(unit_points.size()); std::vector> function_gradients(unit_points.size()); - std::vector dof_indices(fe.dofs_per_cell); - std::vector function_values_1(unit_points.size()); - std::vector> function_gradients_1(unit_points.size()); - for (const auto &cell : dof_handler.active_cell_iterators()) { fe_values.reinit(cell); @@ -103,15 +99,14 @@ test(const unsigned int degree) evaluator.evaluate(cell, unit_points, solution_values, - function_values_1, - function_gradients_1); + EvaluationFlags::values | EvaluationFlags::gradients); deallog << "Cell with center " << cell->center(true) << std::endl; - for (unsigned int i = 0; i < function_values_1.size(); ++i) + for (unsigned int i = 0; i < function_values.size(); ++i) deallog << mapping.transform_unit_to_real_cell(cell, unit_points[i]) - << ": " << function_values_1[i] << " error value " - << function_values[i] - function_values_1[i] << " error grad " - << (function_gradients_1[i] - function_gradients[i]).norm() + << ": " << evaluator.get_value(i) << " error value " + << function_values[i] - evaluator.get_value(i) << " error grad " + << (evaluator.get_gradient(i) - function_gradients[i]).norm() << std::endl; deallog << std::endl; } diff --git a/tests/fe/point_evaluation_03.cc b/tests/fe/point_evaluation_03.cc index a7bc6e1663..a51f158dd0 100644 --- a/tests/fe/point_evaluation_03.cc +++ b/tests/fe/point_evaluation_03.cc @@ -99,10 +99,7 @@ test(const unsigned int degree) std::vector> function_values(unit_points.size()); std::vector> function_gradients(unit_points.size()); - std::vector dof_indices(fe.dofs_per_cell); - std::vector> function_values_1(unit_points.size()); - std::vector> function_gradients_1(unit_points.size()); - FEValuesExtractors::Vector extractor(0); + FEValuesExtractors::Vector extractor(0); for (const auto &cell : dof_handler.active_cell_iterators()) { @@ -117,16 +114,15 @@ test(const unsigned int degree) evaluator.evaluate(cell, unit_points, solution_values, - function_values_1, - function_gradients_1); + EvaluationFlags::values | EvaluationFlags::gradients); deallog << "Cell with center " << cell->center(true) << std::endl; - for (unsigned int i = 0; i < function_values_1.size(); ++i) + for (unsigned int i = 0; i < function_values.size(); ++i) deallog << mapping.transform_unit_to_real_cell(cell, unit_points[i]) - << ": " << function_values_1[i] << " error value " - << (function_values[i] - function_values_1[i]).norm() + << ": " << evaluator.get_value(i) << " error value " + << (function_values[i] - evaluator.get_value(i)).norm() << " error grad " - << (function_gradients_1[i] - function_gradients[i]).norm() + << (evaluator.get_gradient(i) - function_gradients[i]).norm() << std::endl; deallog << std::endl; } diff --git a/tests/fe/point_evaluation_04.cc b/tests/fe/point_evaluation_04.cc index 5625ed676d..c9df96a372 100644 --- a/tests/fe/point_evaluation_04.cc +++ b/tests/fe/point_evaluation_04.cc @@ -105,13 +105,6 @@ test(const unsigned int degree) std::vector function_values_scalar(unit_points.size()); std::vector> function_gradients_scalar(unit_points.size()); - std::vector dof_indices(fe.dofs_per_cell); - - std::vector> function_values_1(unit_points.size()); - std::vector> function_gradients_1(unit_points.size()); - std::vector function_values_scalar_1(unit_points.size()); - std::vector> function_gradients_scalar_1(unit_points.size()); - FEValuesExtractors::Vector extractor(0); FEValuesExtractors::Scalar extractor_s(dim); @@ -132,26 +125,25 @@ test(const unsigned int degree) evaluator.evaluate(cell, unit_points, solution_values, - function_values_1, - function_gradients_1); + EvaluationFlags::values | EvaluationFlags::gradients); evaluator_scalar.evaluate(cell, unit_points, solution_values, - function_values_scalar_1, - function_gradients_scalar_1); + EvaluationFlags::values | + EvaluationFlags::gradients); deallog << "Cell with center " << cell->center(true) << std::endl; - for (unsigned int i = 0; i < function_values_1.size(); ++i) + for (unsigned int i = 0; i < function_values.size(); ++i) deallog << mapping.transform_unit_to_real_cell(cell, unit_points[i]) - << ": " << function_values_1[i] << " error value " - << (function_values[i] - function_values_1[i]).norm() + << ": " << evaluator.get_value(i) << " error value " + << (function_values[i] - evaluator.get_value(i)).norm() << " error grad " - << (function_gradients_1[i] - function_gradients[i]).norm() + << (evaluator.get_gradient(i) - function_gradients[i]).norm() << " error scalar value " - << function_values_scalar[i] - function_values_scalar_1[i] + << function_values_scalar[i] - evaluator_scalar.get_value(i) << " error scalar grad " << (function_gradients_scalar[i] - - function_gradients_scalar_1[i]) + evaluator_scalar.get_gradient(i)) .norm() << std::endl; deallog << std::endl; diff --git a/tests/fe/point_evaluation_05.cc b/tests/fe/point_evaluation_05.cc index a7728cf4cc..3f88d8100e 100644 --- a/tests/fe/point_evaluation_05.cc +++ b/tests/fe/point_evaluation_05.cc @@ -94,10 +94,6 @@ test(const unsigned int degree) std::vector function_values(unit_points.size()); std::vector> function_gradients(unit_points.size()); - std::vector dof_indices(fe.dofs_per_cell); - std::vector function_values_1(unit_points.size()); - std::vector> function_gradients_1(unit_points.size()); - for (const auto &cell : dof_handler.active_cell_iterators()) { fe_values.reinit(cell); @@ -111,15 +107,14 @@ test(const unsigned int degree) evaluator.evaluate(cell, unit_points, solution_values, - function_values_1, - function_gradients_1); + EvaluationFlags::values | EvaluationFlags::gradients); deallog << "Cell with center " << cell->center(true) << std::endl; - for (unsigned int i = 0; i < function_values_1.size(); ++i) + for (unsigned int i = 0; i < function_values.size(); ++i) deallog << mapping.transform_unit_to_real_cell(cell, unit_points[i]) - << ": " << function_values_1[i] << " error value " - << function_values[i] - function_values_1[i] << " error grad " - << (function_gradients_1[i] - function_gradients[i]).norm() + << ": " << evaluator.get_value(i) << " error value " + << function_values[i] - evaluator.get_value(i) << " error grad " + << (evaluator.get_gradient(i) - function_gradients[i]).norm() << std::endl; deallog << std::endl; } diff --git a/tests/fe/point_evaluation_06.cc b/tests/fe/point_evaluation_06.cc index b76a91971d..34516a6292 100644 --- a/tests/fe/point_evaluation_06.cc +++ b/tests/fe/point_evaluation_06.cc @@ -95,10 +95,6 @@ test(const unsigned int degree) std::vector function_values(unit_points.size()); std::vector> function_gradients(unit_points.size()); - std::vector dof_indices(fe.dofs_per_cell); - std::vector function_values_1(unit_points.size()); - std::vector> function_gradients_1(unit_points.size()); - for (const auto &cell : dof_handler.active_cell_iterators()) { fe_values.reinit(cell); @@ -112,15 +108,14 @@ test(const unsigned int degree) evaluator.evaluate(cell, unit_points, solution_values, - function_values_1, - function_gradients_1); + EvaluationFlags::values | EvaluationFlags::gradients); deallog << "Cell with center " << cell->center(true) << std::endl; - for (unsigned int i = 0; i < function_values_1.size(); ++i) + for (unsigned int i = 0; i < function_values.size(); ++i) deallog << mapping.transform_unit_to_real_cell(cell, unit_points[i]) - << ": " << function_values_1[i] << " error value " - << function_values[i] - function_values_1[i] << " error grad " - << (function_gradients_1[i] - function_gradients[i]).norm() + << ": " << evaluator.get_value(i) << " error value " + << function_values[i] - evaluator.get_value(i) << " error grad " + << (evaluator.get_gradient(i) - function_gradients[i]).norm() << std::endl; deallog << std::endl; } diff --git a/tests/fe/point_evaluation_07.cc b/tests/fe/point_evaluation_07.cc index b921dbb17b..cfceb95494 100644 --- a/tests/fe/point_evaluation_07.cc +++ b/tests/fe/point_evaluation_07.cc @@ -106,10 +106,7 @@ test(const unsigned int degree) std::vector> function_values(unit_points.size()); std::vector> function_gradients(unit_points.size()); - std::vector dof_indices(fe.dofs_per_cell); - std::vector> function_values_1(unit_points.size()); - std::vector> function_gradients_1(unit_points.size()); - FEValuesExtractors::Vector extractor(0); + FEValuesExtractors::Vector extractor(0); for (const auto &cell : dof_handler.active_cell_iterators()) { @@ -124,16 +121,15 @@ test(const unsigned int degree) evaluator.evaluate(cell, unit_points, solution_values, - function_values_1, - function_gradients_1); + EvaluationFlags::values | EvaluationFlags::gradients); deallog << "Cell with center " << cell->center(true) << std::endl; - for (unsigned int i = 0; i < function_values_1.size(); ++i) + for (unsigned int i = 0; i < function_values.size(); ++i) deallog << mapping.transform_unit_to_real_cell(cell, unit_points[i]) - << ": " << function_values_1[i] << " error value " - << (function_values[i] - function_values_1[i]).norm() + << ": " << evaluator.get_value(i) << " error value " + << (function_values[i] - evaluator.get_value(i)).norm() << " error grad " - << (function_gradients_1[i] - function_gradients[i]).norm() + << (evaluator.get_gradient(i) - function_gradients[i]).norm() << std::endl; deallog << std::endl; } diff --git a/tests/fe/point_evaluation_08.cc b/tests/fe/point_evaluation_08.cc index 5512e2da59..dbac85db7e 100644 --- a/tests/fe/point_evaluation_08.cc +++ b/tests/fe/point_evaluation_08.cc @@ -112,13 +112,6 @@ test(const unsigned int degree) std::vector function_values_scalar(unit_points.size()); std::vector> function_gradients_scalar(unit_points.size()); - std::vector dof_indices(fe.dofs_per_cell); - - std::vector> function_values_1(unit_points.size()); - std::vector> function_gradients_1(unit_points.size()); - std::vector function_values_scalar_1(unit_points.size()); - std::vector> function_gradients_scalar_1(unit_points.size()); - FEValuesExtractors::Vector extractor(0); FEValuesExtractors::Scalar extractor_s(dim); @@ -139,26 +132,25 @@ test(const unsigned int degree) evaluator.evaluate(cell, unit_points, solution_values, - function_values_1, - function_gradients_1); + EvaluationFlags::values | EvaluationFlags::gradients); evaluator_scalar.evaluate(cell, unit_points, solution_values, - function_values_scalar_1, - function_gradients_scalar_1); + EvaluationFlags::values | + EvaluationFlags::gradients); deallog << "Cell with center " << cell->center(true) << std::endl; - for (unsigned int i = 0; i < function_values_1.size(); ++i) + for (unsigned int i = 0; i < function_values.size(); ++i) deallog << mapping.transform_unit_to_real_cell(cell, unit_points[i]) - << ": " << function_values_1[i] << " error value " - << (function_values[i] - function_values_1[i]).norm() + << ": " << evaluator.get_value(i) << " error value " + << (function_values[i] - evaluator.get_value(i)).norm() << " error grad " - << (function_gradients_1[i] - function_gradients[i]).norm() + << (evaluator.get_gradient(i) - function_gradients[i]).norm() << " error scalar value " - << function_values_scalar[i] - function_values_scalar_1[i] + << function_values_scalar[i] - evaluator_scalar.get_value(i) << " error scalar grad " << (function_gradients_scalar[i] - - function_gradients_scalar_1[i]) + evaluator_scalar.get_gradient(i)) .norm() << std::endl; deallog << std::endl; diff --git a/tests/fe/point_evaluation_09.cc b/tests/fe/point_evaluation_09.cc index 2650e9dc0b..7df599f2a3 100644 --- a/tests/fe/point_evaluation_09.cc +++ b/tests/fe/point_evaluation_09.cc @@ -101,11 +101,6 @@ test(const unsigned int degree) std::vector>> function_gradients( unit_points.size(), std::vector>(dim + 1)); - std::vector dof_indices(fe.dofs_per_cell); - std::vector> function_values_1(unit_points.size()); - std::vector>> function_gradients_1( - unit_points.size()); - for (const auto &cell : dof_handler.active_cell_iterators()) { fe_values.reinit(cell); @@ -119,22 +114,22 @@ test(const unsigned int degree) evaluator.evaluate(cell, unit_points, solution_values, - function_values_1, - function_gradients_1); + EvaluationFlags::values | EvaluationFlags::gradients); deallog << "Cell with center " << cell->center(true) << std::endl; - for (unsigned int i = 0; i < function_values_1.size(); ++i) + for (unsigned int i = 0; i < function_values.size(); ++i) { deallog << mapping.transform_unit_to_real_cell(cell, unit_points[i]) - << ": " << function_values_1[i] << " error value "; + << ": " << evaluator.get_value(i) << " error value "; double error = 0; for (unsigned int d = 0; d < dim + 1; ++d) - error += std::abs(function_values[i][d] - function_values_1[i][d]); + error += + std::abs(function_values[i][d] - evaluator.get_value(i)[d]); deallog << error << " error grad "; error = 0; for (unsigned int d = 0; d < dim + 1; ++d) error += - (function_gradients[i][d] - function_gradients_1[i][d]).norm(); + (function_gradients[i][d] - evaluator.get_gradient(i)[d]).norm(); deallog << error << std::endl; } deallog << std::endl; -- 2.39.5