14#ifndef dealii_vector_tools_evaluation_h
15#define dealii_vector_tools_evaluation_h
137 template <
int n_components,
152 const MeshType<dim, spacedim> &mesh,
153 const VectorType &vector,
161 const
unsigned int first_selected_component = 0);
178 template <
int n_components,
185 (
concepts::is_dealii_vector_type<VectorType> &&
186 concepts::is_triangulation_or_dof_handler<MeshType<dim, spacedim>>))
191 typename VectorType::value_type>::
193 RemotePointEvaluation<dim, spacedim> &cache,
194 const MeshType<dim, spacedim> &mesh,
195 const VectorType &vector,
198 const
unsigned int first_selected_component = 0);
213 template <
int n_components,
220 (
concepts::is_dealii_vector_type<VectorType> &&
221 concepts::is_triangulation_or_dof_handler<MeshType<dim, spacedim>>))
226 typename VectorType::value_type>::
228 const MeshType<dim, spacedim> &mesh,
229 const VectorType &vector,
238 first_selected_component = 0);
254 template <
int n_components,
261 (
concepts::is_dealii_vector_type<VectorType> &&
262 concepts::is_triangulation_or_dof_handler<MeshType<dim, spacedim>>))
267 typename VectorType::value_type>::
269 RemotePointEvaluation<dim, spacedim>
271 const MeshType<dim, spacedim> &mesh,
272 const VectorType &vector,
276 first_selected_component = 0);
284 template <
int n_components,
297 typename VectorType::value_type>::
299 const MeshType<dim, spacedim> &mesh,
300 const VectorType &vector,
307 const
unsigned int first_selected_component)
309 cache.reinit(evaluation_points, mesh.get_triangulation(),
mapping);
312 cache, mesh, vector, flags, first_selected_component);
317 template <
int n_components,
332 const MeshType<dim, spacedim> &mesh,
333 const VectorType &vector,
342 first_selected_component)
344 cache.reinit(evaluation_points, mesh.get_triangulation(),
mapping);
347 cache, mesh, vector, flags, first_selected_component);
357 template <
typename T>
360 const ArrayView<const T> &values)
364 case EvaluationFlags::avg:
366 return std::accumulate(
values.begin(),
values.end(), T{}) /
369 case EvaluationFlags::max:
371 case EvaluationFlags::min:
373 case EvaluationFlags::insert:
386 template <
int rank,
int dim,
typename Number>
387 Tensor<rank, dim, Number>
389 const ArrayView<
const Tensor<rank, dim, Number>> &values)
393 case EvaluationFlags::avg:
395 return std::accumulate(
values.begin(),
397 Tensor<rank, dim, Number>{}) /
400 case EvaluationFlags::insert:
414 template <
int n_components,
int rank,
int dim,
typename Number>
415 Tensor<1, n_components, Tensor<rank, dim, Number>>
418 const ArrayView<
const Tensor<1, n_components, Tensor<rank, dim, Number>>>
423 case EvaluationFlags::avg:
425 Tensor<1, n_components, Tensor<rank, dim, Number>> temp;
427 for (
unsigned int j = 0; j <
values.size(); ++j)
428 for (
unsigned int i = 0; i < n_components; ++i)
429 temp[i] = temp[i] + values[j][i];
431 for (
unsigned int i = 0; i < n_components; ++i)
436 case EvaluationFlags::insert:
446 template <
int n_components,
453 const unsigned int i,
454 const typename Utilities::MPI::RemotePointEvaluation<dim,
457 const Utilities::MPI::RemotePointEvaluation<dim, spacedim> &cache,
458 const DoFHandler<dim, spacedim> &dof_handler,
459 const VectorType &vector,
461 const ::EvaluationFlags::EvaluationFlags evaluation_flags,
462 const unsigned int first_selected_component,
464 value_type(
const FEPointEvaluation<n_components,
467 typename VectorType::value_type> &,
468 const unsigned int &)> process_quadrature_point,
469 const ArrayView<value_type> &values,
470 std::vector<typename VectorType::value_type> &solution_values,
472 std::unique_ptr<FEPointEvaluation<n_components,
475 typename VectorType::value_type>>>
478 if (evaluators.empty())
483 cell_data.cells[i].first,
484 cell_data.cells[i].second,
487 const ArrayView<const Point<dim>> unit_points(
488 cell_data.reference_point_values.data() +
489 cell_data.reference_point_ptrs[i],
490 cell_data.reference_point_ptrs[i + 1] -
491 cell_data.reference_point_ptrs[i]);
493 solution_values.resize(
495 cell->get_dof_values(vector,
496 solution_values.begin(),
497 solution_values.end());
499 if (evaluators[cell->active_fe_index()] ==
nullptr)
500 evaluators[cell->active_fe_index()] =
501 std::make_unique<FEPointEvaluation<n_components,
504 typename VectorType::value_type>>(
508 first_selected_component);
509 auto &evaluator = *evaluators[cell->active_fe_index()];
511 evaluator.reinit(cell, unit_points);
512 evaluator.evaluate(solution_values, evaluation_flags);
514 for (
unsigned int q = 0; q < unit_points.size(); ++q)
515 values[q + cell_data.reference_point_ptrs[i]] =
516 process_quadrature_point(evaluator, q);
521 template <
int dim,
int spacedim,
typename Number>
524 const Triangulation<dim, spacedim> &tria,
525 const Vector<Number> &vector,
529 return vector[cell->active_cell_index()];
534 template <
int dim,
int spacedim,
typename Number>
537 const Triangulation<dim, spacedim> &tria,
538 const LinearAlgebra::distributed::Vector<Number> &vector,
541 const auto distributed_tria =
542 dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *
>(&tria);
544 const bool use_distributed_path =
545 (distributed_tria ==
nullptr) ?
547 (vector.get_partitioner().get() ==
548 distributed_tria->global_active_cell_index_partitioner()
552 if (use_distributed_path)
554 return vector[cell->global_active_cell_index()];
559 return vector[cell->active_cell_index()];
565 template <
typename Number,
typename Number2>
567 set_value(Number &dst,
const Number2 &src)
574 template <
typename Number,
int rank,
int dim,
typename Number2>
576 set_value(Tensor<rank, dim, Number> &,
const Number2 &)
580 "A cell-data vector can only have a single component."));
585 template <
int n_components,
592 const unsigned int i,
593 const typename Utilities::MPI::RemotePointEvaluation<dim,
596 const Utilities::MPI::RemotePointEvaluation<dim, spacedim> &,
597 const Triangulation<dim, spacedim> &triangulation,
598 const VectorType &vector,
600 const ::EvaluationFlags::EvaluationFlags evaluation_flags,
601 const unsigned int first_selected_component,
603 value_type(
const FEPointEvaluation<n_components,
606 typename VectorType::value_type> &,
607 const unsigned int &)>,
608 const ArrayView<value_type> &values,
609 std::vector<typename VectorType::value_type> &,
611 std::unique_ptr<FEPointEvaluation<n_components,
614 typename VectorType::value_type>>> &)
616 Assert(n_components == 1 && first_selected_component == 0,
618 "A cell-data vector can only have a single component."));
620 Assert(evaluation_flags ==
622 ExcMessage(
"For cell-data vectors, only values can be queried."));
625 &triangulation, cell_data.cells[i].first, cell_data.cells[i].second};
629 for (
unsigned int q = cell_data.reference_point_ptrs[i];
630 q < cell_data.reference_point_ptrs[i + 1];
637 template <
int n_components,
644 concepts::is_dealii_vector_type<VectorType>
645 &&concepts::is_triangulation_or_dof_handler<MeshType>)
646 inline std::vector<value_type> evaluate_at_points(
647 const Utilities::MPI::RemotePointEvaluation<dim, spacedim> &cache,
648 const MeshType &mesh,
649 const VectorType &vector,
651 const unsigned int first_selected_component,
653 const ::EvaluationFlags::EvaluationFlags evaluation_flags,
655 value_type(
const FEPointEvaluation<n_components,
658 typename VectorType::value_type> &,
659 const unsigned int &)> process_quadrature_point)
663 "Utilities::MPI::RemotePointEvaluation is not ready yet! "
664 "Please call Utilities::MPI::RemotePointEvaluation::reinit() "
665 "yourself or another function that does this for you."));
670 "The provided Utilities::MPI::RemotePointEvaluation and DoFHandler "
671 "object have been set up with different Triangulation objects, "
672 "a scenario not supported!"));
675 const auto evaluation_point_results = [&]() {
678 const auto evaluation_function = [&](
auto &
values,
679 const auto &cell_data) {
680 std::vector<typename VectorType::value_type> solution_values;
683 std::unique_ptr<FEPointEvaluation<n_components,
686 typename VectorType::value_type>>>
689 for (
unsigned int i = 0; i < cell_data.cells.size(); ++i)
690 process_cell<n_components, dim, spacedim, VectorType, value_type>(
698 first_selected_component,
699 process_quadrature_point,
705 std::vector<value_type> evaluation_point_results;
706 std::vector<value_type> buffer;
708 cache.template evaluate_and_process<value_type>(
709 evaluation_point_results, buffer, evaluation_function);
711 return evaluation_point_results;
717 return evaluation_point_results;
723 std::vector<value_type> unique_evaluation_point_results(
728 for (
unsigned int i = 0; i < ptr.size() - 1; ++i)
730 const auto n_entries = ptr[i + 1] - ptr[i];
734 unique_evaluation_point_results[i] =
736 ArrayView<const value_type>(
737 evaluation_point_results.data() + ptr[i], n_entries));
740 return unique_evaluation_point_results;
745 template <
int n_components,
752 (concepts::is_dealii_vector_type<VectorType> &&
753 concepts::is_triangulation_or_dof_handler<MeshType<dim, spacedim>>))
755 typename FEPointEvaluation<n_components,
760 RemotePointEvaluation<dim, spacedim> &cache,
761 const MeshType<dim, spacedim> &mesh,
762 const VectorType &vector,
764 const
unsigned int first_selected_component)
766 return internal::evaluate_at_points<
770 MeshType<dim, spacedim>,
772 typename FEPointEvaluation<n_components,
780 first_selected_component,
783 [](
const auto &evaluator,
const auto &q) {
784 return evaluator.get_value(q);
790 template <
int n_components,
797 (concepts::is_dealii_vector_type<VectorType> &&
798 concepts::is_triangulation_or_dof_handler<MeshType<dim, spacedim>>))
800 typename FEPointEvaluation<n_components,
805 RemotePointEvaluation<dim, spacedim>
807 const MeshType<dim, spacedim> &mesh,
808 const VectorType &vector,
812 first_selected_component)
814 return internal::evaluate_at_points<
818 MeshType<dim, spacedim>,
820 typename FEPointEvaluation<
824 typename VectorType::value_type>::gradient_type>(
829 first_selected_component,
832 [](
const auto &evaluator,
const unsigned &q) {
833 return evaluator.get_gradient(q);
const hp::FECollection< dim, spacedim > & get_fe_collection() const
const FiniteElement< dim, spacedim > & get_fe(const types::fe_index index=0) const
unsigned int n_dofs_per_cell() const
size_type locally_owned_size() 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
virtual size_type size() const override
unsigned int size() const
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_CXX20_REQUIRES(condition)
#define DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NOT_IMPLEMENTED()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcMessage(std::string arg1)
TriaActiveIterator< CellAccessor< dim, spacedim > > active_cell_iterator
typename ActiveSelector::active_cell_iterator active_cell_iterator
const bool IsBlockVector< VectorType >::value
@ update_values
Shape function values.
@ update_gradients
Shape function gradients.
MappingQ< dim, spacedim > StaticMappingQ1< dim, spacedim >::mapping
The namespace for the EvaluationFlags enum.
EvaluationFlags
The EvaluationFlags enum.
Tpetra::Vector< Number, LO, GO, NodeType< MemorySpace > > VectorType
* * * * ValueType TimeRateRequest< ValueType, dim, Number > get_value() const
* * * RotationFunction< dim, Number >::RotationFunction Number(dim)
* * * * void TimeRateRequest< ValueType, dim, Number > set_value(const ValueType &value)
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)