Reference documentation for deal.II version 9.5.0
\(\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// Copyright (C) 2021 - 2023 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
16
17#ifndef dealii_vector_tools_evaluation_h
18#define dealii_vector_tools_evaluation_h
19
20#include <deal.II/base/config.h>
21
23
25
27
29#include <deal.II/lac/vector.h>
30
32
34
35namespace VectorTools
36{
41 {
46 {
50 avg = 0,
56 max = 1,
62 min = 2,
66 insert = 3
67 };
68 } // namespace EvaluationFlags
69
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:
376 Assert(false, ExcNotImplemented());
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:
403 Assert(false, ExcNotImplemented());
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:
439 Assert(false, ExcNotImplemented());
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.size() == 0)
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 (void)tria;
530 return vector[cell->active_cell_index()];
531 }
532
533
534
535 template <int dim, int spacedim, typename Number>
536 Number
537 get_value(
541 {
542 const auto distributed_tria =
543 dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(&tria);
544
545 const bool use_distributed_path =
546 (distributed_tria == nullptr) ?
547 false :
548 (vector.get_partitioner().get() ==
549 distributed_tria->global_active_cell_index_partitioner()
550 .lock()
551 .get());
552
553 if (use_distributed_path)
554 {
555 return vector[cell->global_active_cell_index()];
556 }
557 else
558 {
560 return vector[cell->active_cell_index()];
561 }
562 }
563
564
565
566 template <typename Number, typename Number2>
567 void
568 set_value(Number &dst, const Number2 &src)
569 {
570 dst = src;
571 }
572
573
574
575 template <typename Number, int rank, int dim, typename Number2>
576 void
577 set_value(Tensor<rank, dim, Number> &, const Number2 &)
578 {
579 Assert(false,
581 "A cell-data vector can only have a single component."));
582 }
583
584
585
586 template <int n_components,
587 int dim,
588 int spacedim,
589 typename VectorType,
590 typename value_type>
591 void
592 process_cell(
593 const unsigned int i,
594 const typename Utilities::MPI::RemotePointEvaluation<dim,
595 spacedim>::CellData
596 &cell_data,
599 const VectorType & vector,
600 const UpdateFlags,
601 const ::EvaluationFlags::EvaluationFlags evaluation_flags,
602 const unsigned int first_selected_component,
603 const std::function<
604 value_type(const FEPointEvaluation<n_components,
605 dim,
606 spacedim,
607 typename VectorType::value_type> &,
608 const unsigned int &)>,
609 const ArrayView<value_type> &values,
610 std::vector<typename VectorType::value_type> &,
611 std::vector<
612 std::unique_ptr<FEPointEvaluation<n_components,
613 dim,
614 spacedim,
615 typename VectorType::value_type>>> &)
616 {
617 (void)evaluation_flags;
618 (void)first_selected_component;
619
620 Assert(n_components == 1 && first_selected_component == 0,
622 "A cell-data vector can only have a single component."));
623
624 Assert(evaluation_flags ==
626 ExcMessage("For cell-data vectors, only values can be queried."));
627
629 &triangulation, cell_data.cells[i].first, cell_data.cells[i].second};
630
631 const auto value = get_value(triangulation, vector, cell);
632
633 for (unsigned int q = cell_data.reference_point_ptrs[i];
634 q < cell_data.reference_point_ptrs[i + 1];
635 ++q)
636 set_value(values[q], value);
637 }
638
639
640
641 template <int n_components,
642 int dim,
643 int spacedim,
644 typename MeshType,
645 typename VectorType,
646 typename value_type>
650 inline std::vector<value_type> evaluate_at_points(
652 const MeshType & mesh,
653 const VectorType & vector,
655 const unsigned int first_selected_component,
656 const UpdateFlags update_flags,
657 const ::EvaluationFlags::EvaluationFlags evaluation_flags,
658 const std::function<
659 value_type(const FEPointEvaluation<n_components,
660 dim,
661 spacedim,
662 typename VectorType::value_type> &,
663 const unsigned int &)> process_quadrature_point)
664 {
665 Assert(cache.is_ready(),
667 "Utilities::MPI::RemotePointEvaluation is not ready yet! "
668 "Please call Utilities::MPI::RemotePointEvaluation::reinit() "
669 "yourself or another function that does this for you."));
670
671 Assert(
672 &mesh.get_triangulation() == &cache.get_triangulation(),
674 "The provided Utilities::MPI::RemotePointEvaluation and DoFHandler "
675 "object have been set up with different Triangulation objects, "
676 "a scenario not supported!"));
677
678 // evaluate values at points if possible
679 const auto evaluation_point_results = [&]() {
680 // helper function for accessing the global vector and interpolating
681 // the results onto the points
682 const auto evaluation_function = [&](auto & values,
683 const auto &cell_data) {
684 std::vector<typename VectorType::value_type> solution_values;
685
686 std::vector<
687 std::unique_ptr<FEPointEvaluation<n_components,
688 dim,
689 spacedim,
690 typename VectorType::value_type>>>
691 evaluators;
692
693 for (unsigned int i = 0; i < cell_data.cells.size(); ++i)
694 process_cell<n_components, dim, spacedim, VectorType, value_type>(
695 i,
696 cell_data,
697 cache,
698 mesh,
699 vector,
700 update_flags,
701 evaluation_flags,
702 first_selected_component,
703 process_quadrature_point,
704 values,
705 solution_values,
706 evaluators);
707 };
708
709 std::vector<value_type> evaluation_point_results;
710 std::vector<value_type> buffer;
711
712 cache.template evaluate_and_process<value_type>(
713 evaluation_point_results, buffer, evaluation_function);
714
715 return evaluation_point_results;
716 }();
717
718 if (cache.is_map_unique())
719 {
720 // each point has exactly one result (unique map)
721 return evaluation_point_results;
722 }
723 else
724 {
725 // map is not unique (multiple or no results): postprocessing is
726 // needed
727 std::vector<value_type> unique_evaluation_point_results(
728 cache.get_point_ptrs().size() - 1);
729
730 const auto &ptr = cache.get_point_ptrs();
731
732 for (unsigned int i = 0; i < ptr.size() - 1; ++i)
733 {
734 const auto n_entries = ptr[i + 1] - ptr[i];
735 if (n_entries == 0)
736 continue;
737
738 unique_evaluation_point_results[i] =
739 reduce(flags,
741 evaluation_point_results.data() + ptr[i], n_entries));
742 }
743
744 return unique_evaluation_point_results;
745 }
746 }
747 } // namespace internal
748
749 template <int n_components,
750 template <int, int>
751 class MeshType,
752 int dim,
753 int spacedim,
754 typename VectorType>
757 concepts::is_triangulation_or_dof_handler<MeshType<dim, spacedim>>))
758 inline std::vector<
759 typename FEPointEvaluation<n_components,
760 dim,
761 spacedim,
762 typename VectorType::value_type>::
763 value_type> point_values(const Utilities::MPI::
764 RemotePointEvaluation<dim, spacedim> &cache,
765 const MeshType<dim, spacedim> & mesh,
766 const VectorType & vector,
767 const EvaluationFlags::EvaluationFlags flags,
768 const unsigned int first_selected_component)
769 {
770 return internal::evaluate_at_points<
771 n_components,
772 dim,
773 spacedim,
774 MeshType<dim, spacedim>,
775 VectorType,
776 typename FEPointEvaluation<n_components,
777 dim,
778 spacedim,
779 typename VectorType::value_type>::value_type>(
780 cache,
781 mesh,
782 vector,
783 flags,
784 first_selected_component,
787 [](const auto &evaluator, const auto &q) {
788 return evaluator.get_value(q);
789 });
790 }
791
792
793
794 template <int n_components,
795 template <int, int>
796 class MeshType,
797 int dim,
798 int spacedim,
799 typename VectorType>
802 concepts::is_triangulation_or_dof_handler<MeshType<dim, spacedim>>))
803 inline std::vector<
804 typename FEPointEvaluation<n_components,
805 dim,
806 spacedim,
807 typename VectorType::value_type>::
808 gradient_type> point_gradients(const Utilities::MPI::
809 RemotePointEvaluation<dim, spacedim>
810 & cache,
811 const MeshType<dim, spacedim> &mesh,
812 const VectorType & vector,
814 flags,
815 const unsigned int
816 first_selected_component)
817 {
818 return internal::evaluate_at_points<
819 n_components,
820 dim,
821 spacedim,
822 MeshType<dim, spacedim>,
823 VectorType,
824 typename FEPointEvaluation<
825 n_components,
826 dim,
827 spacedim,
828 typename VectorType::value_type>::gradient_type>(
829 cache,
830 mesh,
831 vector,
832 flags,
833 first_selected_component,
836 [](const auto &evaluator, const unsigned &q) {
837 return evaluator.get_gradient(q);
838 });
839 }
840
841#endif
842} // namespace VectorTools
843
845
846#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:317
Definition point.h:112
unsigned int n_active_cells() const
const Triangulation< dim, spacedim > & get_triangulation() const
const std::vector< unsigned int > & get_point_ptrs() const
const Mapping< dim, spacedim > & get_mapping() const
size_type size() const
unsigned int size() const
Definition collection.h:265
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_CXX20_REQUIRES(condition)
Definition config.h:160
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
static ::ExceptionBase & ExcNotImplemented()
#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.
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.
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
const ::Triangulation< dim, spacedim > & tria