#include <deal.II/base/mpi_remote_point_evaluation.h>
+#include <deal.II/distributed/tria_base.h>
+
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/matrix_free/fe_point_evaluation.h>
* 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()`
* @warning This is a collective call that needs to be executed by all
* processors in the communicator.
*/
- template <int n_components, int dim, int spacedim, typename VectorType>
+ template <int n_components,
+ template <int, int>
+ class MeshType,
+ int dim,
+ int spacedim,
+ typename VectorType>
std::vector<
typename FEPointEvaluation<n_components,
dim,
typename VectorType::value_type>::value_type>
point_values(
const Mapping<dim> & mapping,
- const DoFHandler<dim, spacedim> & dof_handler,
+ const MeshType<dim, spacedim> & mesh,
const VectorType & vector,
const std::vector<Point<spacedim>> & evaluation_points,
Utilities::MPI::RemotePointEvaluation<dim, spacedim> &cache,
* @warning This is a collective call that needs to be executed by all
* processors in the communicator.
*/
- template <int n_components, int dim, int spacedim, typename VectorType>
+ template <int n_components,
+ template <int, int>
+ class MeshType,
+ int dim,
+ int spacedim,
+ typename VectorType>
std::vector<
typename FEPointEvaluation<n_components,
dim,
typename VectorType::value_type>::value_type>
point_values(
const Utilities::MPI::RemotePointEvaluation<dim, spacedim> &cache,
- const DoFHandler<dim, spacedim> & dof_handler,
+ const MeshType<dim, spacedim> & mesh,
const VectorType & vector,
const EvaluationFlags::EvaluationFlags flags = EvaluationFlags::avg);
* @warning This is a collective call that needs to be executed by all
* processors in the communicator.
*/
- template <int n_components, int dim, int spacedim, typename VectorType>
+ template <int n_components,
+ template <int, int>
+ class MeshType,
+ int dim,
+ int spacedim,
+ typename VectorType>
std::vector<
typename FEPointEvaluation<n_components,
dim,
typename VectorType::value_type>::gradient_type>
point_gradients(
const Mapping<dim> & mapping,
- const DoFHandler<dim, spacedim> & dof_handler,
+ const MeshType<dim, spacedim> & mesh,
const VectorType & vector,
const std::vector<Point<spacedim>> & evaluation_points,
Utilities::MPI::RemotePointEvaluation<dim, spacedim> &cache,
* @warning This is a collective call that needs to be executed by all
* processors in the communicator.
*/
- template <int n_components, int dim, int spacedim, typename VectorType>
+ template <int n_components,
+ template <int, int>
+ class MeshType,
+ int dim,
+ int spacedim,
+ typename VectorType>
std::vector<
typename FEPointEvaluation<n_components,
dim,
typename VectorType::value_type>::gradient_type>
point_gradients(
const Utilities::MPI::RemotePointEvaluation<dim, spacedim> &cache,
- const DoFHandler<dim, spacedim> & dof_handler,
+ const MeshType<dim, spacedim> & mesh,
const VectorType & vector,
const EvaluationFlags::EvaluationFlags flags = EvaluationFlags::avg);
#ifndef DOXYGEN
- template <int n_components, int dim, int spacedim, typename VectorType>
+ template <int n_components,
+ template <int, int>
+ class MeshType,
+ int dim,
+ int spacedim,
+ typename VectorType>
inline std::vector<
typename FEPointEvaluation<n_components,
dim,
spacedim,
typename VectorType::value_type>::value_type>
point_values(const Mapping<dim> & mapping,
- const DoFHandler<dim, spacedim> & dof_handler,
+ const MeshType<dim, spacedim> & mesh,
const VectorType & vector,
const std::vector<Point<spacedim>> &evaluation_points,
Utilities::MPI::RemotePointEvaluation<dim, spacedim> &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<n_components>(cache, dof_handler, vector, flags);
+ return point_values<n_components>(cache, mesh, vector, flags);
}
- template <int n_components, int dim, int spacedim, typename VectorType>
+ template <int n_components,
+ template <int, int>
+ class MeshType,
+ int dim,
+ int spacedim,
+ typename VectorType>
inline std::vector<
typename FEPointEvaluation<n_components,
dim,
spacedim,
typename VectorType::value_type>::gradient_type>
point_gradients(const Mapping<dim> & mapping,
- const DoFHandler<dim, spacedim> & dof_handler,
+ const MeshType<dim, spacedim> & mesh,
const VectorType & vector,
const std::vector<Point<spacedim>> &evaluation_points,
Utilities::MPI::RemotePointEvaluation<dim, spacedim> &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<n_components>(cache, dof_handler, vector, flags);
+ return point_gradients<n_components>(cache, mesh, vector, flags);
}
int spacedim,
typename VectorType,
typename value_type>
+ void
+ process_cell(
+ const unsigned int i,
+ const typename Utilities::MPI::RemotePointEvaluation<dim,
+ spacedim>::CellData
+ & cell_data,
+ const Utilities::MPI::RemotePointEvaluation<dim, spacedim> &cache,
+ const DoFHandler<dim, spacedim> & dof_handler,
+ const VectorType & vector,
+ const UpdateFlags update_flags,
+ const dealii::EvaluationFlags::EvaluationFlags evaluation_flags,
+ const std::function<
+ value_type(const FEPointEvaluation<n_components,
+ dim,
+ spacedim,
+ typename VectorType::value_type> &,
+ const unsigned int &)> process_quadrature_point,
+ const ArrayView<value_type> & values,
+ std::vector<typename VectorType::value_type> &solution_values,
+ std::vector<
+ std::unique_ptr<FEPointEvaluation<n_components,
+ dim,
+ spacedim,
+ typename VectorType::value_type>>>
+ &evaluators)
+ {
+ if (evaluators.size() == 0)
+ evaluators.resize(dof_handler.get_fe_collection().size());
+
+ typename DoFHandler<dim>::active_cell_iterator cell = {
+ &cache.get_triangulation(),
+ cell_data.cells[i].first,
+ cell_data.cells[i].second,
+ &dof_handler};
+
+ const ArrayView<const Point<dim>> 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<n_components,
+ dim,
+ spacedim,
+ typename VectorType::value_type>>(
+ 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 <int dim, int spacedim, typename Number>
+ Number
+ get_value(
+ const Triangulation<dim, spacedim> & tria,
+ const Vector<Number> & vector,
+ const typename Triangulation<dim, spacedim>::active_cell_iterator &cell)
+ {
+ (void)tria;
+ AssertDimension(tria.n_active_cells(), vector.size());
+ return vector[cell->active_cell_index()];
+ }
+
+
+
+ template <int dim, int spacedim, typename Number>
+ Number
+ get_value(
+ const Triangulation<dim, spacedim> & tria,
+ const LinearAlgebra::distributed::Vector<Number> & vector,
+ const typename Triangulation<dim, spacedim>::active_cell_iterator &cell)
+ {
+ const auto distributed_tria =
+ dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(&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 <int n_components,
+ int dim,
+ int spacedim,
+ typename VectorType,
+ typename value_type>
+ void
+ process_cell(
+ const unsigned int i,
+ const typename Utilities::MPI::RemotePointEvaluation<dim,
+ spacedim>::CellData
+ &cell_data,
+ const Utilities::MPI::RemotePointEvaluation<dim, spacedim> &,
+ const Triangulation<dim, spacedim> &triangulation,
+ const VectorType & vector,
+ const UpdateFlags,
+ const dealii::EvaluationFlags::EvaluationFlags evaluation_flags,
+ const std::function<
+ value_type(const FEPointEvaluation<n_components,
+ dim,
+ spacedim,
+ typename VectorType::value_type> &,
+ const unsigned int &)>,
+ const ArrayView<value_type> &values,
+ std::vector<typename VectorType::value_type> &,
+ std::vector<
+ std::unique_ptr<FEPointEvaluation<n_components,
+ dim,
+ spacedim,
+ typename VectorType::value_type>>> &)
+ {
+ (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<dim>::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 <int n_components,
+ int dim,
+ int spacedim,
+ typename MeshType,
+ typename VectorType,
+ typename value_type>
inline std::vector<value_type>
evaluate_at_points(
const Utilities::MPI::RemotePointEvaluation<dim, spacedim> &cache,
- const DoFHandler<dim, spacedim> & dof_handler,
+ const MeshType & mesh,
const VectorType & vector,
const EvaluationFlags::EvaluationFlags flags,
const UpdateFlags update_flags,
"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, "
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<dim>::active_cell_iterator cell = {
- &cache.get_triangulation(),
- cell_data.cells[i].first,
- cell_data.cells[i].second,
- &dof_handler};
-
- const ArrayView<const Point<dim>> 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<n_components,
- dim,
- spacedim,
- typename VectorType::value_type>>(
- 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<n_components, dim, spacedim, VectorType, value_type>(
+ i,
+ cell_data,
+ cache,
+ mesh,
+ vector,
+ update_flags,
+ evaluation_flags,
+ process_quadrature_point,
+ values,
+ solution_values,
+ evaluators);
};
std::vector<value_type> evaluation_point_results;
}
} // namespace internal
- template <int n_components, int dim, int spacedim, typename VectorType>
+ template <int n_components,
+ template <int, int>
+ class MeshType,
+ int dim,
+ int spacedim,
+ typename VectorType>
inline std::vector<
typename FEPointEvaluation<n_components,
dim,
typename VectorType::value_type>::value_type>
point_values(
const Utilities::MPI::RemotePointEvaluation<dim, spacedim> &cache,
- const DoFHandler<dim, spacedim> & dof_handler,
+ const MeshType<dim, spacedim> & mesh,
const VectorType & vector,
const EvaluationFlags::EvaluationFlags flags)
{
n_components,
dim,
spacedim,
+ MeshType<dim, spacedim>,
VectorType,
typename FEPointEvaluation<n_components,
dim,
spacedim,
typename VectorType::value_type>::value_type>(
cache,
- dof_handler,
+ mesh,
vector,
flags,
update_values,
});
}
- template <int n_components, int dim, int spacedim, typename VectorType>
+ template <int n_components,
+ template <int, int>
+ class MeshType,
+ int dim,
+ int spacedim,
+ typename VectorType>
inline std::vector<
typename FEPointEvaluation<n_components,
dim,
typename VectorType::value_type>::gradient_type>
point_gradients(
const Utilities::MPI::RemotePointEvaluation<dim, spacedim> &cache,
- const DoFHandler<dim, spacedim> & dof_handler,
+ const MeshType<dim, spacedim> & mesh,
const VectorType & vector,
const EvaluationFlags::EvaluationFlags flags)
{
n_components,
dim,
spacedim,
+ MeshType<dim, spacedim>,
VectorType,
typename FEPointEvaluation<
n_components,
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);
});
}
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/mapping_fe.h>
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/lac/la_parallel_vector.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/solver_control.h>
+
+#include <deal.II/matrix_free/fe_point_evaluation.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/numerics/vector_tools_evaluate.h>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+template <int dim, int spacedim>
+std::shared_ptr<const Utilities::MPI::Partitioner>
+create_partitioner(const DoFHandler<dim, spacedim> &dof_handler)
+{
+ IndexSet locally_relevant_dofs;
+
+ DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs);
+
+ return std::make_shared<const Utilities::MPI::Partitioner>(
+ 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<dim> tria_1(MPI_COMM_WORLD);
+ GridGenerator::hyper_cube(tria_1, -1.0, +1.0);
+ tria_1.refine_global(n_refinements_1);
+ DoFHandler<dim> dof_handler_1(tria_1);
+
+ FE_DGQ<dim> fe_1(0);
+ dof_handler_1.distribute_dofs(fe_1);
+
+ // approach 1:
+ LinearAlgebra::distributed::Vector<double> vector_1(
+ create_partitioner(dof_handler_1));
+
+ // approach 2:
+ LinearAlgebra::distributed::Vector<double> vector_2(
+ tria_1.global_active_cell_index_partitioner().lock());
+
+ // approach 3:
+ Vector<double> 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<double> 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<dim> mapping_1;
+
+ std::vector<Point<dim>> 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<dim> 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();
+}
projects / dealii.git / commitdiff