deal.II version GIT relicensing-6809-ge913b9bb34 2026-09-25 17:20:01+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
vector_tools_evaluate.h
Go to the documentation of this file.
1// -----------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
4// Copyright (C) 2021 - 2025 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Detailed license information governing the source code and contributions
9// can be found in LICENSE.md and CONTRIBUTING.md at the top level directory.
10//
11// -----------------------------------------------------------------------------
12
13
14#ifndef dealii_vector_tools_evaluation_h
15#define dealii_vector_tools_evaluation_h
16
17#include <deal.II/base/config.h>
18
20
22
24
26#include <deal.II/lac/vector.h>
27
29
31
32namespace VectorTools
33{
38 {
43 {
47 avg = 0,
53 max = 1,
59 min = 2,
63 insert = 3
64 };
65 } // namespace EvaluationFlags
66
137 template <int n_components,
138 template <int, int>
139 class MeshType,
140 int dim,
141 int spacedim,
142 typename VectorType>
145 concepts::is_triangulation_or_dof_handler<MeshType<dim, spacedim>>))
146 std::vector<
147 typename FEPointEvaluation<n_components,
148 dim,
149 spacedim,
150 typename VectorType::value_type>::
151 value_type> point_values(const Mapping<dim> &mapping,
152 const MeshType<dim, spacedim> &mesh,
153 const VectorType &vector,
154 const std::vector<Point<spacedim>>
155 &evaluation_points,
156 Utilities::MPI::RemotePointEvaluation<dim,
157 spacedim>
158 &cache,
159 const EvaluationFlags::EvaluationFlags flags =
160 EvaluationFlags::avg,
161 const unsigned int first_selected_component = 0);
162
178 template <int n_components,
179 template <int, int>
180 class MeshType,
181 int dim,
182 int spacedim,
183 typename VectorType>
185 (concepts::is_dealii_vector_type<VectorType> &&
186 concepts::is_triangulation_or_dof_handler<MeshType<dim, spacedim>>))
187 std::vector<
188 typename FEPointEvaluation<n_components,
189 dim,
190 spacedim,
191 typename VectorType::value_type>::
192 value_type> point_values(const Utilities::MPI::
193 RemotePointEvaluation<dim, spacedim> &cache,
194 const MeshType<dim, spacedim> &mesh,
195 const VectorType &vector,
196 const EvaluationFlags::EvaluationFlags flags =
197 EvaluationFlags::avg,
198 const unsigned int first_selected_component = 0);
199
213 template <int n_components,
214 template <int, int>
215 class MeshType,
216 int dim,
217 int spacedim,
218 typename VectorType>
220 (concepts::is_dealii_vector_type<VectorType> &&
221 concepts::is_triangulation_or_dof_handler<MeshType<dim, spacedim>>))
222 std::vector<
223 typename FEPointEvaluation<n_components,
224 dim,
225 spacedim,
226 typename VectorType::value_type>::
227 gradient_type> point_gradients(const Mapping<dim> &mapping,
228 const MeshType<dim, spacedim> &mesh,
229 const VectorType &vector,
230 const std::vector<Point<spacedim>>
231 &evaluation_points,
232 Utilities::MPI::RemotePointEvaluation<
233 dim,
234 spacedim> &cache,
236 flags = EvaluationFlags::avg,
237 const unsigned int
238 first_selected_component = 0);
239
254 template <int n_components,
255 template <int, int>
256 class MeshType,
257 int dim,
258 int spacedim,
259 typename VectorType>
261 (concepts::is_dealii_vector_type<VectorType> &&
262 concepts::is_triangulation_or_dof_handler<MeshType<dim, spacedim>>))
263 std::vector<
264 typename FEPointEvaluation<n_components,
265 dim,
266 spacedim,
267 typename VectorType::value_type>::
268 gradient_type> point_gradients(const Utilities::MPI::
269 RemotePointEvaluation<dim, spacedim>
270 &cache,
271 const MeshType<dim, spacedim> &mesh,
272 const VectorType &vector,
274 flags = EvaluationFlags::avg,
275 const unsigned int
276 first_selected_component = 0);
277
278
279
280 // inlined functions
281
282
283#ifndef DOXYGEN
284 template <int n_components,
285 template <int, int>
286 class MeshType,
287 int dim,
288 int spacedim,
289 typename VectorType>
292 concepts::is_triangulation_or_dof_handler<MeshType<dim, spacedim>>))
293 inline std::vector<
294 typename FEPointEvaluation<n_components,
295 dim,
296 spacedim,
297 typename VectorType::value_type>::
298 value_type> point_values(const Mapping<dim> &mapping,
299 const MeshType<dim, spacedim> &mesh,
300 const VectorType &vector,
301 const std::vector<Point<spacedim>>
302 &evaluation_points,
303 Utilities::MPI::RemotePointEvaluation<dim,
304 spacedim>
305 &cache,
306 const EvaluationFlags::EvaluationFlags flags,
307 const unsigned int first_selected_component)
308 {
309 cache.reinit(evaluation_points, mesh.get_triangulation(), mapping);
310
311 return point_values<n_components>(
312 cache, mesh, vector, flags, first_selected_component);
313 }
314
315
316
317 template <int n_components,
318 template <int, int>
319 class MeshType,
320 int dim,
321 int spacedim,
322 typename VectorType>
325 concepts::is_triangulation_or_dof_handler<MeshType<dim, spacedim>>))
326 inline std::vector<
327 typename FEPointEvaluation<n_components,
328 dim,
329 spacedim,
330 typename VectorType::value_type>::
331 gradient_type> point_gradients(const Mapping<dim> &mapping,
332 const MeshType<dim, spacedim> &mesh,
333 const VectorType &vector,
334 const std::vector<Point<spacedim>>
335 &evaluation_points,
336 Utilities::MPI::RemotePointEvaluation<
337 dim,
338 spacedim> &cache,
340 flags,
341 const unsigned int
342 first_selected_component)
343 {
344 cache.reinit(evaluation_points, mesh.get_triangulation(), mapping);
345
346 return point_gradients<n_components>(
347 cache, mesh, vector, flags, first_selected_component);
348 }
349
350
351
352 namespace internal
353 {
357 template <typename T>
358 T
360 const ArrayView<const T> &values)
361 {
362 switch (flags)
363 {
364 case EvaluationFlags::avg:
365 {
366 return std::accumulate(values.begin(), values.end(), T{}) /
367 (T(1.0) * values.size());
368 }
369 case EvaluationFlags::max:
370 return *std::max_element(values.begin(), values.end());
371 case EvaluationFlags::min:
372 return *std::min_element(values.begin(), values.end());
373 case EvaluationFlags::insert:
374 return values[0];
375 default:
377 return values[0];
378 }
379 }
380
381
382
386 template <int rank, int dim, typename Number>
389 const ArrayView<const Tensor<rank, dim, Number>> &values)
390 {
391 switch (flags)
392 {
393 case EvaluationFlags::avg:
394 {
395 return std::accumulate(values.begin(),
396 values.end(),
398 (Number(1.0) * values.size());
399 }
400 case EvaluationFlags::insert:
401 return values[0];
402 default:
404 return values[0];
405 }
406 }
407
408
409
414 template <int n_components, int rank, int dim, typename Number>
416 reduce(
418 const ArrayView<const Tensor<1, n_components, Tensor<rank, dim, Number>>>
419 &values)
420 {
421 switch (flags)
422 {
423 case EvaluationFlags::avg:
424 {
426
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];
430
431 for (unsigned int i = 0; i < n_components; ++i)
432 temp[i] /= Number(values.size());
433
434 return temp;
435 }
436 case EvaluationFlags::insert:
437 return values[0];
438 default:
440 return values[0];
441 }
442 }
443
444
445
446 template <int n_components,
447 int dim,
448 int spacedim,
449 typename VectorType,
450 typename value_type>
451 void
452 process_cell(
453 const unsigned int i,
454 const typename Utilities::MPI::RemotePointEvaluation<dim,
455 spacedim>::CellData
456 &cell_data,
458 const DoFHandler<dim, spacedim> &dof_handler,
459 const VectorType &vector,
460 const UpdateFlags update_flags,
461 const ::EvaluationFlags::EvaluationFlags evaluation_flags,
462 const unsigned int first_selected_component,
463 const std::function<
464 value_type(const FEPointEvaluation<n_components,
465 dim,
466 spacedim,
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,
471 std::vector<
472 std::unique_ptr<FEPointEvaluation<n_components,
473 dim,
474 spacedim,
475 typename VectorType::value_type>>>
476 &evaluators)
477 {
478 if (evaluators.empty())
479 evaluators.resize(dof_handler.get_fe_collection().size());
480
482 &cache.get_triangulation(),
483 cell_data.cells[i].first,
484 cell_data.cells[i].second,
485 &dof_handler};
486
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]);
492
493 solution_values.resize(
494 dof_handler.get_fe(cell->active_fe_index()).n_dofs_per_cell());
495 cell->get_dof_values(vector,
496 solution_values.begin(),
497 solution_values.end());
498
499 if (evaluators[cell->active_fe_index()] == nullptr)
500 evaluators[cell->active_fe_index()] =
501 std::make_unique<FEPointEvaluation<n_components,
502 dim,
503 spacedim,
504 typename VectorType::value_type>>(
505 cache.get_mapping(),
506 cell->get_fe(),
507 update_flags,
508 first_selected_component);
509 auto &evaluator = *evaluators[cell->active_fe_index()];
510
511 evaluator.reinit(cell, unit_points);
512 evaluator.evaluate(solution_values, evaluation_flags);
513
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);
517 }
518
519
520
521 template <int dim, int spacedim, typename Number>
522 Number
523 get_value(
525 const Vector<Number> &vector,
527 {
528 AssertDimension(tria.n_active_cells(), vector.size());
529 return vector[cell->active_cell_index()];
530 }
531
532
533
534 template <int dim, int spacedim, typename Number>
535 Number
536 get_value(
540 {
541 const auto distributed_tria =
542 dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(&tria);
543
544 const bool use_distributed_path =
545 (distributed_tria == nullptr) ?
546 false :
547 (vector.get_partitioner().get() ==
548 distributed_tria->global_active_cell_index_partitioner()
549 .lock()
550 .get());
551
552 if (use_distributed_path)
553 {
554 return vector[cell->global_active_cell_index()];
555 }
556 else
557 {
559 return vector[cell->active_cell_index()];
560 }
561 }
562
563
564
565 template <typename Number, typename Number2>
566 void
567 set_value(Number &dst, const Number2 &src)
568 {
569 dst = src;
570 }
571
572
573
574 template <typename Number, int rank, int dim, typename Number2>
575 void
576 set_value(Tensor<rank, dim, Number> &, const Number2 &)
577 {
578 Assert(false,
580 "A cell-data vector can only have a single component."));
581 }
582
583
584
585 template <int n_components,
586 int dim,
587 int spacedim,
588 typename VectorType,
589 typename value_type>
590 void
591 process_cell(
592 const unsigned int i,
593 const typename Utilities::MPI::RemotePointEvaluation<dim,
594 spacedim>::CellData
595 &cell_data,
597 const Triangulation<dim, spacedim> &triangulation,
598 const VectorType &vector,
599 const UpdateFlags,
600 const ::EvaluationFlags::EvaluationFlags evaluation_flags,
601 const unsigned int first_selected_component,
602 const std::function<
603 value_type(const FEPointEvaluation<n_components,
604 dim,
605 spacedim,
606 typename VectorType::value_type> &,
607 const unsigned int &)>,
608 const ArrayView<value_type> &values,
609 std::vector<typename VectorType::value_type> &,
610 std::vector<
611 std::unique_ptr<FEPointEvaluation<n_components,
612 dim,
613 spacedim,
614 typename VectorType::value_type>>> &)
615 {
616 Assert(n_components == 1 && first_selected_component == 0,
618 "A cell-data vector can only have a single component."));
619
620 Assert(evaluation_flags ==
622 ExcMessage("For cell-data vectors, only values can be queried."));
623
625 &triangulation, cell_data.cells[i].first, cell_data.cells[i].second};
626
627 const auto value = get_value(triangulation, vector, cell);
628
629 for (unsigned int q = cell_data.reference_point_ptrs[i];
630 q < cell_data.reference_point_ptrs[i + 1];
631 ++q)
632 set_value(values[q], value);
633 }
634
635
636
637 template <int n_components,
638 int dim,
639 int spacedim,
640 typename MeshType,
641 typename VectorType,
642 typename value_type>
646 inline std::vector<value_type> evaluate_at_points(
648 const MeshType &mesh,
649 const VectorType &vector,
651 const unsigned int first_selected_component,
652 const UpdateFlags update_flags,
653 const ::EvaluationFlags::EvaluationFlags evaluation_flags,
654 const std::function<
655 value_type(const FEPointEvaluation<n_components,
656 dim,
657 spacedim,
658 typename VectorType::value_type> &,
659 const unsigned int &)> process_quadrature_point)
660 {
661 Assert(cache.is_ready(),
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."));
666
667 Assert(
668 &mesh.get_triangulation() == &cache.get_triangulation(),
670 "The provided Utilities::MPI::RemotePointEvaluation and DoFHandler "
671 "object have been set up with different Triangulation objects, "
672 "a scenario not supported!"));
673
674 // evaluate values at points if possible
675 const auto evaluation_point_results = [&]() {
676 // helper function for accessing the global vector and interpolating
677 // the results onto the points
678 const auto evaluation_function = [&](auto &values,
679 const auto &cell_data) {
680 std::vector<typename VectorType::value_type> solution_values;
681
682 std::vector<
683 std::unique_ptr<FEPointEvaluation<n_components,
684 dim,
685 spacedim,
686 typename VectorType::value_type>>>
687 evaluators;
688
689 for (unsigned int i = 0; i < cell_data.cells.size(); ++i)
690 process_cell<n_components, dim, spacedim, VectorType, value_type>(
691 i,
692 cell_data,
693 cache,
694 mesh,
695 vector,
696 update_flags,
697 evaluation_flags,
698 first_selected_component,
699 process_quadrature_point,
700 values,
701 solution_values,
702 evaluators);
703 };
704
705 std::vector<value_type> evaluation_point_results;
706 std::vector<value_type> buffer;
707
708 cache.template evaluate_and_process<value_type>(
709 evaluation_point_results, buffer, evaluation_function);
710
711 return evaluation_point_results;
712 }();
713
714 if (cache.is_map_unique())
715 {
716 // each point has exactly one result (unique map)
717 return evaluation_point_results;
718 }
719 else
720 {
721 // map is not unique (multiple or no results): postprocessing is
722 // needed
723 std::vector<value_type> unique_evaluation_point_results(
724 cache.get_point_ptrs().size() - 1);
725
726 const auto &ptr = cache.get_point_ptrs();
727
728 for (unsigned int i = 0; i < ptr.size() - 1; ++i)
729 {
730 const auto n_entries = ptr[i + 1] - ptr[i];
731 if (n_entries == 0)
732 continue;
733
734 unique_evaluation_point_results[i] =
735 reduce(flags,
737 evaluation_point_results.data() + ptr[i], n_entries));
738 }
739
740 return unique_evaluation_point_results;
741 }
742 }
743 } // namespace internal
744
745 template <int n_components,
746 template <int, int>
747 class MeshType,
748 int dim,
749 int spacedim,
750 typename VectorType>
753 concepts::is_triangulation_or_dof_handler<MeshType<dim, spacedim>>))
754 inline std::vector<
755 typename FEPointEvaluation<n_components,
756 dim,
757 spacedim,
758 typename VectorType::value_type>::
760 RemotePointEvaluation<dim, spacedim> &cache,
761 const MeshType<dim, spacedim> &mesh,
762 const VectorType &vector,
763 const EvaluationFlags::EvaluationFlags flags,
764 const unsigned int first_selected_component)
765 {
766 return internal::evaluate_at_points<
767 n_components,
768 dim,
769 spacedim,
770 MeshType<dim, spacedim>,
772 typename FEPointEvaluation<n_components,
773 dim,
774 spacedim,
775 typename VectorType::value_type>::value_type>(
776 cache,
777 mesh,
778 vector,
779 flags,
780 first_selected_component,
783 [](const auto &evaluator, const auto &q) {
784 return evaluator.get_value(q);
785 });
786 }
787
788
789
790 template <int n_components,
791 template <int, int>
792 class MeshType,
793 int dim,
794 int spacedim,
795 typename VectorType>
798 concepts::is_triangulation_or_dof_handler<MeshType<dim, spacedim>>))
799 inline std::vector<
800 typename FEPointEvaluation<n_components,
801 dim,
802 spacedim,
803 typename VectorType::value_type>::
804 gradient_type> point_gradients(const Utilities::MPI::
805 RemotePointEvaluation<dim, spacedim>
806 &cache,
807 const MeshType<dim, spacedim> &mesh,
808 const VectorType &vector,
810 flags,
811 const unsigned int
812 first_selected_component)
813 {
814 return internal::evaluate_at_points<
815 n_components,
816 dim,
817 spacedim,
818 MeshType<dim, spacedim>,
820 typename FEPointEvaluation<
821 n_components,
822 dim,
823 spacedim,
824 typename VectorType::value_type>::gradient_type>(
825 cache,
826 mesh,
827 vector,
828 flags,
829 first_selected_component,
832 [](const auto &evaluator, const unsigned &q) {
833 return evaluator.get_gradient(q);
834 });
835 }
836
837#endif
838} // namespace VectorTools
839
841
842#endif // dealii_vector_tools_boundary_h
const hp::FECollection< dim, spacedim > & get_fe_collection() const
const FiniteElement< dim, spacedim > & get_fe(const types::fe_index index=0) 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
Abstract base class for mapping classes.
Definition mapping.h:318
Definition point.h:111
unsigned int n_active_cells() const
Communicate values between a mesh and arbitrary points.
const Triangulation< dim, spacedim > & get_triangulation() const
const std::vector< unsigned int > & get_point_ptrs() const
const Mapping< dim, spacedim > & get_mapping() const
virtual size_type size() const override
unsigned int size() const
Definition collection.h:314
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:38
#define DEAL_II_CXX20_REQUIRES(condition)
Definition config.h:249
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:39
#define DEAL_II_NOT_IMPLEMENTED()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcMessage(std::string arg1)
typename ActiveSelector::active_cell_iterator active_cell_iterator
UpdateFlags
@ update_values
Shape function values.
@ update_gradients
Shape function gradients.
The namespace for the EvaluationFlags enum.
EvaluationFlags
The EvaluationFlags enum.
constexpr char T
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)
std::vector< typename FEPointEvaluation< n_components, dim, spacedim, typename VectorType::value_type >::value_type > point_values(const Mapping< dim > &mapping, 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=EvaluationFlags::avg, const unsigned int first_selected_component=0)
std::vector< typename FEPointEvaluation< n_components, dim, spacedim, typename VectorType::value_type >::gradient_type > point_gradients(const Mapping< dim > &mapping, 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=EvaluationFlags::avg, const unsigned int first_selected_component=0)
STL namespace.