From ac2ea78e76794fab2ca9abc83f58394fc79c3cfb Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Sun, 30 May 2021 08:52:56 +0200 Subject: [PATCH] Introduce VectorTools::point_gradients() --- .../deal.II/numerics/vector_tools_evaluate.h | 326 ++++++++++++------ .../vector_tools_evaluate_at_points_04.cc | 6 +- 2 files changed, 219 insertions(+), 113 deletions(-) diff --git a/include/deal.II/numerics/vector_tools_evaluate.h b/include/deal.II/numerics/vector_tools_evaluate.h index cde9f9508f..ec8aec1776 100644 --- a/include/deal.II/numerics/vector_tools_evaluate.h +++ b/include/deal.II/numerics/vector_tools_evaluate.h @@ -32,12 +32,12 @@ DEAL_II_NAMESPACE_OPEN namespace VectorTools { /** - * Namespace for the flags for point_values(). + * Namespace for the flags for point_values() and point_gradients(). */ namespace EvaluationFlags { /** - * Flags for point_values(). + * Flags for point_values() and point_gradients(). */ enum EvaluationFlags { @@ -100,6 +100,36 @@ namespace VectorTools const VectorType & vector, const EvaluationFlags::EvaluationFlags flags = EvaluationFlags::avg); + /** + * Given a (distributed) solution vector @p vector, evaluate the gradients at + * the (arbitrary and even remote) points specified by @p evaluation_points. + */ + template + std::vector::gradient_type> + point_gradients( + const Mapping & mapping, + const DoFHandler & dof_handler, + const VectorType & vector, + const std::vector> & evaluation_points, + Utilities::MPI::RemotePointEvaluation &cache, + const EvaluationFlags::EvaluationFlags flags = EvaluationFlags::avg); + + /** + * Given a (distributed) solution vector @p vector, evaluate the gradients at + * the points specified by @p cache which might have been set up by the + * above function. + * + * @note Refinement/coarsening/repartitioning leads to the invalidation of the + * cache so that the above function has to be called again. + */ + template + std::vector::gradient_type> + point_gradients( + const Utilities::MPI::RemotePointEvaluation &cache, + const DoFHandler & dof_handler, + const VectorType & vector, + const EvaluationFlags::EvaluationFlags flags = EvaluationFlags::avg); + // inlined functions @@ -122,6 +152,23 @@ namespace VectorTools + template + inline std::vector< + typename FEPointEvaluation::gradient_type> + point_gradients(const Mapping & mapping, + const DoFHandler & dof_handler, + const VectorType & vector, + const std::vector> &evaluation_points, + Utilities::MPI::RemotePointEvaluation &cache, + const EvaluationFlags::EvaluationFlags flags) + { + cache.reinit(evaluation_points, dof_handler.get_triangulation(), mapping); + + return point_gradients(cache, dof_handler, vector, flags); + } + + + namespace internal { /** @@ -177,10 +224,133 @@ namespace VectorTools return values[0]; } } - } // namespace internal + template + inline std::vector + evaluate_at_points( + const Utilities::MPI::RemotePointEvaluation &cache, + const DoFHandler & dof_handler, + const VectorType & vector, + const EvaluationFlags::EvaluationFlags flags, + const UpdateFlags update_flags, + const dealii::EvaluationFlags::EvaluationFlags evaluation_flags, + const std::function< + value_type(const FEPointEvaluation &, + const unsigned int &)> process_quadrature_point) + { + Assert(cache.is_ready(), + ExcMessage( + "Utilities::MPI::RemotePointEvaluation is not ready yet! " + "Please call Utilities::MPI::RemotePointEvaluation::reinit() " + "yourself or another function that does this for you.")); + + Assert( + &dof_handler.get_triangulation() == &cache.get_triangulation(), + ExcMessage( + "The provided Utilities::MPI::RemotePointEvaluation and DoFHandler " + "object have been set up with different Triangulation objects, " + "a scenario not supported!")); + + // evaluate values at points if possible + const auto evaluation_point_results = [&]() { + // helper function for accessing the global vector and interpolating + // the results onto the points + const auto evaluation_function = [&](auto & values, + const auto &cell_data) { + std::vector solution_values; + + std::vector>> + evaluators(dof_handler.get_fe_collection().size()); + + const auto get_evaluator = [&](const unsigned int active_fe_index) + -> FEPointEvaluation & { + if (evaluators[active_fe_index] == nullptr) + evaluators[active_fe_index] = + std::make_unique>( + cache.get_mapping(), + dof_handler.get_fe(active_fe_index), + update_flags); + + return *evaluators[active_fe_index]; + }; + + for (unsigned int i = 0; i < cell_data.cells.size(); ++i) + { + typename DoFHandler::active_cell_iterator cell = { + &cache.get_triangulation(), + cell_data.cells[i].first, + cell_data.cells[i].second, + &dof_handler}; + + const ArrayView> unit_points( + cell_data.reference_point_values.data() + + cell_data.reference_point_ptrs[i], + cell_data.reference_point_ptrs[i + 1] - + cell_data.reference_point_ptrs[i]); + + solution_values.resize( + dof_handler.get_fe(cell->active_fe_index()).n_dofs_per_cell()); + cell->get_dof_values(vector, + solution_values.begin(), + solution_values.end()); + + auto &evaluator = get_evaluator(cell->active_fe_index()); + + evaluator.reinit(cell, unit_points); + evaluator.evaluate(solution_values, evaluation_flags); + + for (unsigned int q = 0; q < unit_points.size(); ++q) + values[q + cell_data.reference_point_ptrs[i]] = + process_quadrature_point(evaluator, q); + } + }; + + std::vector evaluation_point_results; + std::vector buffer; + + cache.template evaluate_and_process( + evaluation_point_results, buffer, evaluation_function); + + return evaluation_point_results; + }(); + + if (cache.is_map_unique()) + { + // each point has exactly one result (unique map) + return evaluation_point_results; + } + else + { + // map is not unique (multiple or no results): postprocessing is + // needed + std::vector unique_evaluation_point_results( + cache.get_point_ptrs().size() - 1); + + const auto &ptr = cache.get_point_ptrs(); + + for (unsigned int i = 0; i < ptr.size() - 1; ++i) + { + const auto n_entries = ptr[i + 1] - ptr[i]; + if (n_entries == 0) + continue; + + unique_evaluation_point_results[i] = + reduce(flags, + ArrayView( + evaluation_point_results.data() + ptr[i], n_entries)); + } + + return unique_evaluation_point_results; + } + } + } // namespace internal + template inline std::vector::value_type> point_values( @@ -189,117 +359,49 @@ namespace VectorTools const VectorType & vector, const EvaluationFlags::EvaluationFlags flags) { - using value_type = - typename FEPointEvaluation::value_type; - - Assert(cache.is_ready(), - ExcMessage( - "Utilities::MPI::RemotePointEvaluation is not ready yet! " - "Please call Utilities::MPI::RemotePointEvaluation::reinit() " - "yourself or the other point_values(), which does this for" - "you.")); - - Assert( - &dof_handler.get_triangulation() == &cache.get_triangulation(), - ExcMessage( - "The provided Utilities::MPI::RemotePointEvaluation and DoFHandler " - "object have been set up with different Triangulation objects, " - "a scenario not supported!")); - - // evaluate values at points if possible - const auto evaluation_point_results = [&]() { - // helper function for accessing the global vector and interpolating - // the results onto the points - const auto evaluation_function = [&](auto & values, - const auto &cell_data) { - std::vector solution_values; - - std::vector>> - evaluators(dof_handler.get_fe_collection().size()); - - const auto get_evaluator = [&](const unsigned int active_fe_index) - -> FEPointEvaluation & { - if (evaluators[active_fe_index] == nullptr) - evaluators[active_fe_index] = - std::make_unique>( - cache.get_mapping(), - dof_handler.get_fe(active_fe_index), - update_values); - - return *evaluators[active_fe_index]; - }; + return internal::evaluate_at_points< + n_components, + dim, + spacedim, + VectorType, + typename FEPointEvaluation::value_type>( + cache, + dof_handler, + vector, + flags, + update_values, + dealii::EvaluationFlags::values, + [](const auto &evaluator, const auto &q) { + return evaluator.get_value(q); + }); + } - for (unsigned int i = 0; i < cell_data.cells.size(); ++i) - { - typename DoFHandler::active_cell_iterator cell = { - &cache.get_triangulation(), - cell_data.cells[i].first, - cell_data.cells[i].second, - &dof_handler}; - - const ArrayView> unit_points( - cell_data.reference_point_values.data() + - cell_data.reference_point_ptrs[i], - cell_data.reference_point_ptrs[i + 1] - - cell_data.reference_point_ptrs[i]); - - solution_values.resize( - dof_handler.get_fe(cell->active_fe_index()).n_dofs_per_cell()); - cell->get_dof_values(vector, - solution_values.begin(), - solution_values.end()); - - auto &evaluator = get_evaluator(cell->active_fe_index()); - - evaluator.reinit(cell, unit_points); - evaluator.evaluate(solution_values, - dealii::EvaluationFlags::values); - - for (unsigned int q = 0; q < unit_points.size(); ++q) - values[q + cell_data.reference_point_ptrs[i]] = - evaluator.get_value(q); - } - }; - - std::vector evaluation_point_results; - std::vector buffer; - - cache.template evaluate_and_process(evaluation_point_results, - buffer, - evaluation_function); - - return evaluation_point_results; - }(); - - if (cache.is_map_unique()) - { - // each point has exactly one result (unique map) - return evaluation_point_results; - } - else - { - // map is not unique (multiple or no results): postprocessing is needed - std::vector unique_evaluation_point_results( - cache.get_point_ptrs().size() - 1); - - const auto &ptr = cache.get_point_ptrs(); - - for (unsigned int i = 0; i < ptr.size() - 1; ++i) - { - const auto n_entries = ptr[i + 1] - ptr[i]; - if (n_entries == 0) - continue; - - unique_evaluation_point_results[i] = - internal::reduce(flags, - ArrayView( - evaluation_point_results.data() + ptr[i], - n_entries)); - } - - return unique_evaluation_point_results; - } + template + inline std::vector< + typename FEPointEvaluation::gradient_type> + point_gradients( + const Utilities::MPI::RemotePointEvaluation &cache, + const DoFHandler & dof_handler, + const VectorType & vector, + const EvaluationFlags::EvaluationFlags flags) + { + return internal::evaluate_at_points< + n_components, + dim, + spacedim, + VectorType, + typename FEPointEvaluation::gradient_type>( + cache, + dof_handler, + vector, + flags, + update_gradients, + dealii::EvaluationFlags::gradients, + [](const auto &evaluator, const auto &q) { + return evaluator.get_gradient(q); + }); } + #endif } // namespace VectorTools diff --git a/tests/remote_point_evaluation/vector_tools_evaluate_at_points_04.cc b/tests/remote_point_evaluation/vector_tools_evaluate_at_points_04.cc index 1e15ae2727..c30f0ec9bc 100644 --- a/tests/remote_point_evaluation/vector_tools_evaluate_at_points_04.cc +++ b/tests/remote_point_evaluation/vector_tools_evaluate_at_points_04.cc @@ -155,13 +155,17 @@ test() Utilities::MPI::RemotePointEvaluation evaluation_cache; const auto evaluation_point_results = VectorTools::point_values<1>( mapping_1, dof_handler_1, vector_1, evaluation_points, evaluation_cache); + const auto evaluation_point_gradient_results = + VectorTools::point_gradients<1>( + mapping_1, dof_handler_1, vector_1, evaluation_points, evaluation_cache); vector_1.zero_out_ghosts(); if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0) for (unsigned int i = 0; i <= n_intervals; ++i) { if (std::abs(evaluation_points[i][0] - evaluation_point_results[i]) > - 1e-10) + 1e-10 || + std::abs(evaluation_point_gradient_results[i][0] - 1.0) > 1e-10) Assert(false, ExcNotImplemented()); } -- 2.39.5