Reference documentation for deal.II version GIT 574c7c8486 2023-09-22 10:20:02+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\}}\)
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 
35 namespace VectorTools
36 {
40  namespace EvaluationFlags
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>
144  (concepts::is_dealii_vector_type<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 =
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 =
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>
291  (concepts::is_dealii_vector_type<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>
324  (concepts::is_dealii_vector_type<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  {
365  {
366  return std::accumulate(values.begin(), values.end(), T{}) /
367  (T(1.0) * values.size());
368  }
370  return *std::max_element(values.begin(), values.end());
372  return *std::min_element(values.begin(), values.end());
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>
390  {
391  switch (flags)
392  {
394  {
395  return std::accumulate(values.begin(),
396  values.end(),
398  (Number(1.0) * values.size());
399  }
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  {
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  }
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,
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,
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 
481  typename DoFHandler<dim>::active_cell_iterator cell = {
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,
580  ExcMessage(
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,
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 &)>,
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,
621  ExcMessage(
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>
648  concepts::is_dealii_vector_type<VectorType>
649  &&concepts::is_triangulation_or_dof_handler<MeshType>)
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,
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(),
666  ExcMessage(
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(),
673  ExcMessage(
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>
756  (concepts::is_dealii_vector_type<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>
801  (concepts::is_dealii_vector_type<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 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.
Definition: mapping.h:317
Definition: point.h:112
Definition: tensor.h:516
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
Definition: vector.h:110
virtual size_type size() const override
unsigned int size() const
Definition: collection.h:265
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_CXX20_REQUIRES(condition)
Definition: config.h:166
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
#define Assert(cond, exc)
Definition: exceptions.h:1616
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1789
static ::ExceptionBase & ExcMessage(std::string arg1)
typename ActiveSelector::active_cell_iterator active_cell_iterator
Definition: dof_handler.h:441
UpdateFlags
@ update_values
Shape function values.
@ update_gradients
Shape function gradients.
The namespace for the EvaluationFlags enum.
EvaluationFlags
The EvaluationFlags enum.
static const char T
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 >::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)
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)
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