From e1d2cdae935f08507a329546de51dbca8b939983 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Wed, 6 Jul 2022 11:48:18 +0200 Subject: [PATCH] RPE: allow to evaluate cell-data vector --- doc/news/changes/minor/20220706Munch | 4 + .../deal.II/numerics/vector_tools_evaluate.h | 318 ++++++++++++++---- .../vector_tools_evaluate_at_points_06.cc | 145 ++++++++ ..._points_06.with_p4est=true.mpirun=4.output | 2 + 4 files changed, 408 insertions(+), 61 deletions(-) create mode 100644 doc/news/changes/minor/20220706Munch create mode 100644 tests/remote_point_evaluation/vector_tools_evaluate_at_points_06.cc create mode 100644 tests/remote_point_evaluation/vector_tools_evaluate_at_points_06.with_p4est=true.mpirun=4.output diff --git a/doc/news/changes/minor/20220706Munch b/doc/news/changes/minor/20220706Munch new file mode 100644 index 0000000000..76eb69ba48 --- /dev/null +++ b/doc/news/changes/minor/20220706Munch @@ -0,0 +1,4 @@ +Improved: The function VectorTools::point_values() can now handle +cell-data vectors. +
+(Peter Munch, 2022/07/06) diff --git a/include/deal.II/numerics/vector_tools_evaluate.h b/include/deal.II/numerics/vector_tools_evaluate.h index 15ee805097..8660498930 100644 --- a/include/deal.II/numerics/vector_tools_evaluate.h +++ b/include/deal.II/numerics/vector_tools_evaluate.h @@ -21,6 +21,8 @@ #include +#include + #include #include @@ -108,6 +110,12 @@ namespace VectorTools * cache, dof_handler_2, vector_2); * @endcode * + * The function can also be used to evaluate cell-data vectors. For this + * purpose, one passes in a Triangulation instead of a DoFHandler and a + * vector of size Trinagulation::n_active_cells() or a vector, which + * has been initialized with the partitioner returned by + * parallel::TriangulationBase::global_active_cell_index_partitioner(). + * * @note If a point cannot be found, the result for these points will * be undefined (most probably 0). If you want to be sure that all * points received a valid result, call `cache.all_points_found()` @@ -116,7 +124,12 @@ namespace VectorTools * @warning This is a collective call that needs to be executed by all * processors in the communicator. */ - template + template + class MeshType, + int dim, + int spacedim, + typename VectorType> std::vector< typename FEPointEvaluation::value_type> point_values( const Mapping & mapping, - const DoFHandler & dof_handler, + const MeshType & mesh, const VectorType & vector, const std::vector> & evaluation_points, Utilities::MPI::RemotePointEvaluation &cache, @@ -142,7 +155,12 @@ namespace VectorTools * @warning This is a collective call that needs to be executed by all * processors in the communicator. */ - template + template + class MeshType, + int dim, + int spacedim, + typename VectorType> std::vector< typename FEPointEvaluation::value_type> point_values( const Utilities::MPI::RemotePointEvaluation &cache, - const DoFHandler & dof_handler, + const MeshType & mesh, const VectorType & vector, const EvaluationFlags::EvaluationFlags flags = EvaluationFlags::avg); @@ -164,7 +182,12 @@ namespace VectorTools * @warning This is a collective call that needs to be executed by all * processors in the communicator. */ - template + template + class MeshType, + int dim, + int spacedim, + typename VectorType> std::vector< typename FEPointEvaluation::gradient_type> point_gradients( const Mapping & mapping, - const DoFHandler & dof_handler, + const MeshType & mesh, const VectorType & vector, const std::vector> & evaluation_points, Utilities::MPI::RemotePointEvaluation &cache, @@ -189,7 +212,12 @@ namespace VectorTools * @warning This is a collective call that needs to be executed by all * processors in the communicator. */ - template + template + class MeshType, + int dim, + int spacedim, + typename VectorType> std::vector< typename FEPointEvaluation::gradient_type> point_gradients( const Utilities::MPI::RemotePointEvaluation &cache, - const DoFHandler & dof_handler, + const MeshType & mesh, const VectorType & vector, const EvaluationFlags::EvaluationFlags flags = EvaluationFlags::avg); @@ -207,42 +235,52 @@ namespace VectorTools #ifndef DOXYGEN - template + template + class MeshType, + int dim, + int spacedim, + typename VectorType> inline std::vector< typename FEPointEvaluation::value_type> point_values(const Mapping & mapping, - const DoFHandler & dof_handler, + const MeshType & mesh, 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); + cache.reinit(evaluation_points, mesh.get_triangulation(), mapping); - return point_values(cache, dof_handler, vector, flags); + return point_values(cache, mesh, vector, flags); } - template + template + class MeshType, + int dim, + int spacedim, + typename VectorType> inline std::vector< typename FEPointEvaluation::gradient_type> point_gradients(const Mapping & mapping, - const DoFHandler & dof_handler, + const MeshType & mesh, 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); + cache.reinit(evaluation_points, mesh.get_triangulation(), mapping); - return point_gradients(cache, dof_handler, vector, flags); + return point_gradients(cache, mesh, vector, flags); } @@ -310,10 +348,179 @@ namespace VectorTools int spacedim, typename VectorType, typename value_type> + void + process_cell( + const unsigned int i, + const typename Utilities::MPI::RemotePointEvaluation::CellData + & cell_data, + const Utilities::MPI::RemotePointEvaluation &cache, + const DoFHandler & dof_handler, + const VectorType & vector, + const UpdateFlags update_flags, + const dealii::EvaluationFlags::EvaluationFlags evaluation_flags, + const std::function< + value_type(const FEPointEvaluation &, + const unsigned int &)> process_quadrature_point, + const ArrayView & values, + std::vector &solution_values, + std::vector< + std::unique_ptr>> + &evaluators) + { + if (evaluators.size() == 0) + evaluators.resize(dof_handler.get_fe_collection().size()); + + 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()); + + if (evaluators[cell->active_fe_index()] == nullptr) + evaluators[cell->active_fe_index()] = + std::make_unique>( + cache.get_mapping(), cell->get_fe(), update_flags); + auto &evaluator = *evaluators[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); + } + + + + template + Number + get_value( + const Triangulation & tria, + const Vector & vector, + const typename Triangulation::active_cell_iterator &cell) + { + (void)tria; + AssertDimension(tria.n_active_cells(), vector.size()); + return vector[cell->active_cell_index()]; + } + + + + template + Number + get_value( + const Triangulation & tria, + const LinearAlgebra::distributed::Vector & vector, + const typename Triangulation::active_cell_iterator &cell) + { + const auto distributed_tria = + dynamic_cast *>(&tria); + + const bool use_distributed_path = + (distributed_tria == nullptr) ? + false : + (vector.get_partitioner().get() == + distributed_tria->global_active_cell_index_partitioner() + .lock() + .get()); + + if (use_distributed_path) + { + return vector[cell->global_active_cell_index()]; + } + else + { + AssertDimension(tria.n_active_cells(), vector.locally_owned_size()); + return vector[cell->active_cell_index()]; + } + } + + + + template + void + process_cell( + const unsigned int i, + const typename Utilities::MPI::RemotePointEvaluation::CellData + &cell_data, + const Utilities::MPI::RemotePointEvaluation &, + const Triangulation &triangulation, + const VectorType & vector, + const UpdateFlags, + const dealii::EvaluationFlags::EvaluationFlags evaluation_flags, + const std::function< + value_type(const FEPointEvaluation &, + const unsigned int &)>, + const ArrayView &values, + std::vector &, + std::vector< + std::unique_ptr>> &) + { + (void)evaluation_flags; + + static_assert(n_components == 1, + "A cell-data vector can only have a single component."); + + Assert(evaluation_flags == + dealii::EvaluationFlags::EvaluationFlags::values, + ExcMessage("For cell-data vectors, only values can be queried.")); + + typename Triangulation::active_cell_iterator cell = { + &triangulation, cell_data.cells[i].first, cell_data.cells[i].second}; + + const auto value = get_value(triangulation, vector, cell); + + for (unsigned int q = cell_data.reference_point_ptrs[i]; + q < cell_data.reference_point_ptrs[i + 1]; + ++q) + values[q] = value; + } + + + + template inline std::vector evaluate_at_points( const Utilities::MPI::RemotePointEvaluation &cache, - const DoFHandler & dof_handler, + const MeshType & mesh, const VectorType & vector, const EvaluationFlags::EvaluationFlags flags, const UpdateFlags update_flags, @@ -332,7 +539,7 @@ namespace VectorTools "yourself or another function that does this for you.")); Assert( - &dof_handler.get_triangulation() == &cache.get_triangulation(), + &mesh.get_triangulation() == &cache.get_triangulation(), ExcMessage( "The provided Utilities::MPI::RemotePointEvaluation and DoFHandler " "object have been set up with different Triangulation objects, " @@ -351,44 +558,21 @@ namespace VectorTools dim, spacedim, typename VectorType::value_type>>> - evaluators(dof_handler.get_fe_collection().size()); + evaluators; 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()); - - if (evaluators[cell->active_fe_index()] == nullptr) - evaluators[cell->active_fe_index()] = std::make_unique< - FEPointEvaluation>( - cache.get_mapping(), cell->get_fe(), update_flags); - auto &evaluator = *evaluators[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); - } + process_cell( + i, + cell_data, + cache, + mesh, + vector, + update_flags, + evaluation_flags, + process_quadrature_point, + values, + solution_values, + evaluators); }; std::vector evaluation_point_results; @@ -431,7 +615,12 @@ namespace VectorTools } } // namespace internal - template + template + class MeshType, + int dim, + int spacedim, + typename VectorType> inline std::vector< typename FEPointEvaluation::value_type> point_values( const Utilities::MPI::RemotePointEvaluation &cache, - const DoFHandler & dof_handler, + const MeshType & mesh, const VectorType & vector, const EvaluationFlags::EvaluationFlags flags) { @@ -447,13 +636,14 @@ namespace VectorTools n_components, dim, spacedim, + MeshType, VectorType, typename FEPointEvaluation::value_type>( cache, - dof_handler, + mesh, vector, flags, update_values, @@ -463,7 +653,12 @@ namespace VectorTools }); } - template + template + class MeshType, + int dim, + int spacedim, + typename VectorType> inline std::vector< typename FEPointEvaluation::gradient_type> point_gradients( const Utilities::MPI::RemotePointEvaluation &cache, - const DoFHandler & dof_handler, + const MeshType & mesh, const VectorType & vector, const EvaluationFlags::EvaluationFlags flags) { @@ -479,6 +674,7 @@ namespace VectorTools n_components, dim, spacedim, + MeshType, VectorType, typename FEPointEvaluation< n_components, @@ -486,12 +682,12 @@ namespace VectorTools spacedim, typename VectorType::value_type>::gradient_type>( cache, - dof_handler, + mesh, vector, flags, update_gradients, dealii::EvaluationFlags::gradients, - [](const auto &evaluator, const auto &q) { + [](const auto &evaluator, const unsigned &q) { return evaluator.get_gradient(q); }); } diff --git a/tests/remote_point_evaluation/vector_tools_evaluate_at_points_06.cc b/tests/remote_point_evaluation/vector_tools_evaluate_at_points_06.cc new file mode 100644 index 0000000000..81563f4d6a --- /dev/null +++ b/tests/remote_point_evaluation/vector_tools_evaluate_at_points_06.cc @@ -0,0 +1,145 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2022 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 VectorTools::point_values() for cell-data vectors. + +#include + +#include + +#include +#include + +#include +#include +#include +#include + +#include + +#include +#include +#include +#include + +#include + +#include +#include +#include + +#include "../tests.h" + +using namespace dealii; + +template +std::shared_ptr +create_partitioner(const DoFHandler &dof_handler) +{ + IndexSet locally_relevant_dofs; + + DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs); + + return std::make_shared( + dof_handler.locally_owned_dofs(), + locally_relevant_dofs, + dof_handler.get_communicator()); +} + +void +test() +{ + const unsigned int dim = 2; + const unsigned int n_refinements_1 = 3; + + parallel::distributed::Triangulation tria_1(MPI_COMM_WORLD); + GridGenerator::hyper_cube(tria_1, -1.0, +1.0); + tria_1.refine_global(n_refinements_1); + DoFHandler dof_handler_1(tria_1); + + FE_DGQ fe_1(0); + dof_handler_1.distribute_dofs(fe_1); + + // approach 1: + LinearAlgebra::distributed::Vector vector_1( + create_partitioner(dof_handler_1)); + + // approach 2: + LinearAlgebra::distributed::Vector vector_2( + tria_1.global_active_cell_index_partitioner().lock()); + + // approach 3: + Vector vector_3(tria_1.n_active_cells()); + + for (const auto &cell : dof_handler_1.active_cell_iterators()) + if (cell->is_locally_owned()) + { + const auto value = cell->center()[0]; + + // approach 1: + Vector temp(cell->get_fe().n_dofs_per_cell()); + temp = value; + cell->set_dof_values(temp, vector_1); + + // approach 2: + vector_2[cell->global_active_cell_index()] = value; + + // approach 3: + vector_3[cell->active_cell_index()] = value; + } + + const MappingQ1 mapping_1; + + std::vector> evaluation_points; + + const unsigned int n_intervals = 100; + + if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0) + for (unsigned int i = 0; i < n_intervals; ++i) + { + const double rad = 2.0 * numbers::PI / n_intervals * i; + evaluation_points.emplace_back(0.75 * std::cos(rad), + 0.75 * std::sin(rad)); + } + + Utilities::MPI::RemotePointEvaluation evaluation_cache; + const auto evaluation_point_results_1 = VectorTools::point_values<1>( + mapping_1, dof_handler_1, vector_1, evaluation_points, evaluation_cache); + const auto evaluation_point_results_2 = + VectorTools::point_values<1>(evaluation_cache, tria_1, vector_2); + const auto evaluation_point_results_3 = + VectorTools::point_values<1>(evaluation_cache, tria_1, vector_3); + + if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0) + for (unsigned int i = 0; i < n_intervals; ++i) + { + if (std::abs(evaluation_point_results_1[i] - + evaluation_point_results_2[i]) > 1e-10 || + std::abs(evaluation_point_results_1[i] - + evaluation_point_results_3[i]) > 1e-10) + Assert(false, ExcNotImplemented()); + } + + deallog << "OK!" << std::endl; +} + +int +main(int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi(argc, argv, 1); + initlog(); + + test(); +} diff --git a/tests/remote_point_evaluation/vector_tools_evaluate_at_points_06.with_p4est=true.mpirun=4.output b/tests/remote_point_evaluation/vector_tools_evaluate_at_points_06.with_p4est=true.mpirun=4.output new file mode 100644 index 0000000000..5cfb783b8f --- /dev/null +++ b/tests/remote_point_evaluation/vector_tools_evaluate_at_points_06.with_p4est=true.mpirun=4.output @@ -0,0 +1,2 @@ + +DEAL::OK! -- 2.39.5