17 #ifndef dealii_vector_tools_evaluation_h
18 #define dealii_vector_tools_evaluation_h
134 template <
int n_components,
141 (concepts::is_dealii_vector_type<VectorType> &&
147 typename VectorType::value_type>::
149 const MeshType<dim, spacedim> &mesh,
150 const VectorType & vector,
151 const std::vector<
Point<spacedim>>
153 Utilities::MPI::RemotePointEvaluation<dim,
158 const
unsigned int first_selected_component = 0);
172 template <
int n_components,
185 typename VectorType::value_type>::
187 RemotePointEvaluation<dim, spacedim> &cache,
188 const MeshType<dim, spacedim> & mesh,
189 const VectorType & vector,
192 const
unsigned int first_selected_component = 0);
204 template <
int n_components,
217 typename VectorType::value_type>::
219 const MeshType<dim, spacedim> &mesh,
220 const VectorType & vector,
221 const std::vector<
Point<spacedim>>
229 first_selected_component = 0);
242 template <
int n_components,
255 typename VectorType::value_type>::
257 RemotePointEvaluation<dim, spacedim>
259 const MeshType<dim, spacedim> &mesh,
260 const VectorType & vector,
264 first_selected_component = 0);
272 template <
int n_components,
279 (concepts::is_dealii_vector_type<VectorType> &&
285 typename VectorType::value_type>::
287 const MeshType<dim, spacedim> &mesh,
288 const VectorType & vector,
289 const std::vector<
Point<spacedim>>
291 Utilities::MPI::RemotePointEvaluation<dim,
295 const
unsigned int first_selected_component)
297 cache.reinit(evaluation_points, mesh.get_triangulation(), mapping);
299 return point_values<n_components>(
300 cache, mesh, vector, flags, first_selected_component);
305 template <
int n_components,
312 (concepts::is_dealii_vector_type<VectorType> &&
318 typename VectorType::value_type>::
320 const MeshType<dim, spacedim> &mesh,
321 const VectorType & vector,
322 const std::vector<
Point<spacedim>>
330 first_selected_component)
332 cache.reinit(evaluation_points, mesh.get_triangulation(), mapping);
334 return point_gradients<n_components>(
335 cache, mesh, vector, flags, first_selected_component);
345 template <
typename T>
374 template <
int rank,
int dim,
typename Number>
383 return std::accumulate(
values.begin(),
386 (Number(1.0) *
values.size());
402 template <
int n_components,
int rank,
int dim,
typename Number>
415 for (
unsigned int j = 0; j <
values.size(); ++j)
416 for (
unsigned int i = 0; i < n_components; ++i)
417 temp[i] = temp[i] +
values[j][i];
419 for (
unsigned int i = 0; i < n_components; ++i)
420 temp[i] /= Number(
values.size());
434 template <
int n_components,
441 const unsigned int i,
447 const VectorType & vector,
450 const unsigned int first_selected_component,
455 typename VectorType::value_type> &,
456 const unsigned int &)> process_quadrature_point,
458 std::vector<typename VectorType::value_type> &solution_values,
463 typename VectorType::value_type>>>
466 if (evaluators.size() == 0)
471 cell_data.cells[i].first,
472 cell_data.cells[i].second,
476 cell_data.reference_point_values.data() +
477 cell_data.reference_point_ptrs[i],
478 cell_data.reference_point_ptrs[i + 1] -
479 cell_data.reference_point_ptrs[i]);
481 solution_values.resize(
483 cell->get_dof_values(vector,
484 solution_values.begin(),
485 solution_values.end());
487 if (evaluators[cell->active_fe_index()] ==
nullptr)
488 evaluators[cell->active_fe_index()] =
492 typename VectorType::value_type>>(
496 first_selected_component);
497 auto &evaluator = *evaluators[cell->active_fe_index()];
499 evaluator.
reinit(cell, unit_points);
500 evaluator.evaluate(solution_values, evaluation_flags);
502 for (
unsigned int q = 0; q < unit_points.size(); ++q)
503 values[q + cell_data.reference_point_ptrs[i]] =
504 process_quadrature_point(evaluator, q);
509 template <
int dim,
int spacedim,
typename Number>
518 return vector[cell->active_cell_index()];
523 template <
int dim,
int spacedim,
typename Number>
530 const auto distributed_tria =
533 const bool use_distributed_path =
534 (distributed_tria ==
nullptr) ?
537 distributed_tria->global_active_cell_index_partitioner()
541 if (use_distributed_path)
543 return vector[cell->global_active_cell_index()];
548 return vector[cell->active_cell_index()];
554 template <
typename Number,
typename Number2>
556 set_value(Number &dst,
const Number2 &src)
563 template <
typename Number,
int rank,
int dim,
typename Number2>
569 "A cell-data vector can only have a single component."));
574 template <
int n_components,
581 const unsigned int i,
587 const VectorType & vector,
590 const unsigned int first_selected_component,
595 typename VectorType::value_type> &,
596 const unsigned int &)>,
598 std::vector<typename VectorType::value_type> &,
603 typename VectorType::value_type>>> &)
605 (void)evaluation_flags;
606 (void)first_selected_component;
608 Assert(n_components == 1 && first_selected_component == 0,
610 "A cell-data vector can only have a single component."));
612 Assert(evaluation_flags ==
614 ExcMessage(
"For cell-data vectors, only values can be queried."));
617 &
triangulation, cell_data.cells[i].first, cell_data.cells[i].second};
621 for (
unsigned int q = cell_data.reference_point_ptrs[i];
622 q < cell_data.reference_point_ptrs[i + 1];
624 set_value(
values[q], value);
629 template <
int n_components,
636 concepts::is_dealii_vector_type<VectorType>
637 &&concepts::is_triangulation_or_dof_handler<MeshType>)
638 inline std::vector<value_type> evaluate_at_points(
640 const MeshType & mesh,
641 const VectorType & vector,
643 const unsigned int first_selected_component,
650 typename VectorType::value_type> &,
651 const unsigned int &)> process_quadrature_point)
655 "Utilities::MPI::RemotePointEvaluation is not ready yet! "
656 "Please call Utilities::MPI::RemotePointEvaluation::reinit() "
657 "yourself or another function that does this for you."));
662 "The provided Utilities::MPI::RemotePointEvaluation and DoFHandler "
663 "object have been set up with different Triangulation objects, "
664 "a scenario not supported!"));
667 const auto evaluation_point_results = [&]() {
670 const auto evaluation_function = [&](
auto &
values,
671 const auto &cell_data) {
672 std::vector<typename VectorType::value_type> solution_values;
678 typename VectorType::value_type>>>
681 for (
unsigned int i = 0; i < cell_data.cells.size(); ++i)
682 process_cell<n_components, dim, spacedim, VectorType, value_type>(
690 first_selected_component,
691 process_quadrature_point,
697 std::vector<value_type> evaluation_point_results;
698 std::vector<value_type> buffer;
700 cache.template evaluate_and_process<value_type>(
701 evaluation_point_results, buffer, evaluation_function);
703 return evaluation_point_results;
709 return evaluation_point_results;
715 std::vector<value_type> unique_evaluation_point_results(
720 for (
unsigned int i = 0; i < ptr.size() - 1; ++i)
722 const auto n_entries = ptr[i + 1] - ptr[i];
726 unique_evaluation_point_results[i] =
729 evaluation_point_results.data() + ptr[i], n_entries));
732 return unique_evaluation_point_results;
737 template <
int n_components,
744 (concepts::is_dealii_vector_type<VectorType> &&
750 typename VectorType::value_type>::
752 RemotePointEvaluation<dim, spacedim> &cache,
753 const MeshType<dim, spacedim> & mesh,
754 const VectorType & vector,
756 const
unsigned int first_selected_component)
758 return internal::evaluate_at_points<
762 MeshType<dim, spacedim>,
767 typename VectorType::value_type>::value_type>(
772 first_selected_component,
775 [](
const auto &evaluator,
const auto &q) {
776 return evaluator.get_value(q);
782 template <
int n_components,
789 (concepts::is_dealii_vector_type<VectorType> &&
795 typename VectorType::value_type>::
797 RemotePointEvaluation<dim, spacedim>
799 const MeshType<dim, spacedim> &mesh,
800 const VectorType & vector,
804 first_selected_component)
806 return internal::evaluate_at_points<
810 MeshType<dim, spacedim>,
816 typename VectorType::value_type>::gradient_type>(
821 first_selected_component,
824 [](
const auto &evaluator,
const unsigned &q) {
825 return evaluator.get_gradient(q);
const FiniteElement< dim, spacedim > & get_fe(const types::fe_index index=0) const
const hp::FECollection< dim, spacedim > & get_fe_collection() const
void reinit(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const ArrayView< const Point< dim >> &unit_points)
unsigned int n_dofs_per_cell() const
size_type locally_owned_size() const
const std::shared_ptr< const Utilities::MPI::Partitioner > & get_partitioner() const
Abstract base class for mapping classes.
unsigned int n_active_cells() const
const Triangulation< dim, spacedim > & get_triangulation() const
const std::vector< unsigned int > & get_point_ptrs() const
bool is_map_unique() const
const Mapping< dim, spacedim > & get_mapping() const
unsigned int size() const
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_CXX20_REQUIRES(condition)
#define DEAL_II_NAMESPACE_CLOSE
#define Assert(cond, exc)
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcMessage(std::string arg1)
typename ActiveSelector::active_cell_iterator active_cell_iterator
@ update_values
Shape function values.
@ update_gradients
Shape function gradients.
The namespace for the EvaluationFlags enum.
EvaluationFlags
The EvaluationFlags enum.
T reduce(const T &local_value, const MPI_Comm &comm, const std::function< T(const T &, const T &)> &combiner, const unsigned int root_process=0)
constexpr bool is_triangulation_or_dof_handler
constexpr bool is_dealii_vector_type
concept is_triangulation_or_dof_handler
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
const ::Triangulation< dim, spacedim > & tria