Reference documentation for deal.II version GIT 741a3088b8 2023-04-01 07:05: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 - 2022 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 
134  template <int n_components,
135  template <int, int>
136  class MeshType,
137  int dim,
138  int spacedim,
139  typename VectorType>
141  (concepts::is_dealii_vector_type<VectorType> &&
142  concepts::is_triangulation_or_dof_handler<MeshType<dim, spacedim>>))
143  std::vector<
144  typename FEPointEvaluation<n_components,
145  dim,
146  spacedim,
147  typename VectorType::value_type>::
148  value_type> point_values(const Mapping<dim> & mapping,
149  const MeshType<dim, spacedim> &mesh,
150  const VectorType & vector,
151  const std::vector<Point<spacedim>>
152  &evaluation_points,
153  Utilities::MPI::RemotePointEvaluation<dim,
154  spacedim>
155  & cache,
156  const EvaluationFlags::EvaluationFlags flags =
158  const unsigned int first_selected_component = 0);
159 
172  template <int n_components,
173  template <int, int>
174  class MeshType,
175  int dim,
176  int spacedim,
177  typename VectorType>
179  (concepts::is_dealii_vector_type<VectorType> &&
180  concepts::is_triangulation_or_dof_handler<MeshType<dim, spacedim>>))
181  std::vector<
182  typename FEPointEvaluation<n_components,
183  dim,
184  spacedim,
185  typename VectorType::value_type>::
186  value_type> point_values(const Utilities::MPI::
187  RemotePointEvaluation<dim, spacedim> &cache,
188  const MeshType<dim, spacedim> & mesh,
189  const VectorType & vector,
190  const EvaluationFlags::EvaluationFlags flags =
192  const unsigned int first_selected_component = 0);
193 
204  template <int n_components,
205  template <int, int>
206  class MeshType,
207  int dim,
208  int spacedim,
209  typename VectorType>
211  (concepts::is_dealii_vector_type<VectorType> &&
212  concepts::is_triangulation_or_dof_handler<MeshType<dim, spacedim>>))
213  std::vector<
214  typename FEPointEvaluation<n_components,
215  dim,
216  spacedim,
217  typename VectorType::value_type>::
218  gradient_type> point_gradients(const Mapping<dim> & mapping,
219  const MeshType<dim, spacedim> &mesh,
220  const VectorType & vector,
221  const std::vector<Point<spacedim>>
222  &evaluation_points,
223  Utilities::MPI::RemotePointEvaluation<
224  dim,
225  spacedim> &cache,
227  flags = EvaluationFlags::avg,
228  const unsigned int
229  first_selected_component = 0);
230 
242  template <int n_components,
243  template <int, int>
244  class MeshType,
245  int dim,
246  int spacedim,
247  typename VectorType>
249  (concepts::is_dealii_vector_type<VectorType> &&
250  concepts::is_triangulation_or_dof_handler<MeshType<dim, spacedim>>))
251  std::vector<
252  typename FEPointEvaluation<n_components,
253  dim,
254  spacedim,
255  typename VectorType::value_type>::
256  gradient_type> point_gradients(const Utilities::MPI::
257  RemotePointEvaluation<dim, spacedim>
258  & cache,
259  const MeshType<dim, spacedim> &mesh,
260  const VectorType & vector,
262  flags = EvaluationFlags::avg,
263  const unsigned int
264  first_selected_component = 0);
265 
266 
267 
268  // inlined functions
269 
270 
271 #ifndef DOXYGEN
272  template <int n_components,
273  template <int, int>
274  class MeshType,
275  int dim,
276  int spacedim,
277  typename VectorType>
279  (concepts::is_dealii_vector_type<VectorType> &&
280  concepts::is_triangulation_or_dof_handler<MeshType<dim, spacedim>>))
281  inline std::vector<
282  typename FEPointEvaluation<n_components,
283  dim,
284  spacedim,
285  typename VectorType::value_type>::
286  value_type> point_values(const Mapping<dim> & mapping,
287  const MeshType<dim, spacedim> &mesh,
288  const VectorType & vector,
289  const std::vector<Point<spacedim>>
290  &evaluation_points,
291  Utilities::MPI::RemotePointEvaluation<dim,
292  spacedim>
293  & cache,
294  const EvaluationFlags::EvaluationFlags flags,
295  const unsigned int first_selected_component)
296  {
297  cache.reinit(evaluation_points, mesh.get_triangulation(), mapping);
298 
299  return point_values<n_components>(
300  cache, mesh, vector, flags, first_selected_component);
301  }
302 
303 
304 
305  template <int n_components,
306  template <int, int>
307  class MeshType,
308  int dim,
309  int spacedim,
310  typename VectorType>
312  (concepts::is_dealii_vector_type<VectorType> &&
313  concepts::is_triangulation_or_dof_handler<MeshType<dim, spacedim>>))
314  inline std::vector<
315  typename FEPointEvaluation<n_components,
316  dim,
317  spacedim,
318  typename VectorType::value_type>::
319  gradient_type> point_gradients(const Mapping<dim> & mapping,
320  const MeshType<dim, spacedim> &mesh,
321  const VectorType & vector,
322  const std::vector<Point<spacedim>>
323  &evaluation_points,
324  Utilities::MPI::RemotePointEvaluation<
325  dim,
326  spacedim> &cache,
328  flags,
329  const unsigned int
330  first_selected_component)
331  {
332  cache.reinit(evaluation_points, mesh.get_triangulation(), mapping);
333 
334  return point_gradients<n_components>(
335  cache, mesh, vector, flags, first_selected_component);
336  }
337 
338 
339 
340  namespace internal
341  {
345  template <typename T>
346  T
348  const ArrayView<const T> & values)
349  {
350  switch (flags)
351  {
353  {
354  return std::accumulate(values.begin(), values.end(), T{}) /
355  (T(1.0) * values.size());
356  }
358  return *std::max_element(values.begin(), values.end());
360  return *std::min_element(values.begin(), values.end());
362  return values[0];
363  default:
364  Assert(false, ExcNotImplemented());
365  return values[0];
366  }
367  }
368 
369 
370 
374  template <int rank, int dim, typename Number>
378  {
379  switch (flags)
380  {
382  {
383  return std::accumulate(values.begin(),
384  values.end(),
386  (Number(1.0) * values.size());
387  }
389  return values[0];
390  default:
391  Assert(false, ExcNotImplemented());
392  return values[0];
393  }
394  }
395 
396 
397 
402  template <int n_components, int rank, int dim, typename Number>
404  reduce(
406  const ArrayView<const Tensor<1, n_components, Tensor<rank, dim, Number>>>
407  &values)
408  {
409  switch (flags)
410  {
412  {
414 
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];
418 
419  for (unsigned int i = 0; i < n_components; ++i)
420  temp[i] /= Number(values.size());
421 
422  return temp;
423  }
425  return values[0];
426  default:
427  Assert(false, ExcNotImplemented());
428  return values[0];
429  }
430  }
431 
432 
433 
434  template <int n_components,
435  int dim,
436  int spacedim,
437  typename VectorType,
438  typename value_type>
439  void
440  process_cell(
441  const unsigned int i,
442  const typename Utilities::MPI::RemotePointEvaluation<dim,
443  spacedim>::CellData
444  & cell_data,
446  const DoFHandler<dim, spacedim> & dof_handler,
447  const VectorType & vector,
448  const UpdateFlags update_flags,
450  const unsigned int first_selected_component,
451  const std::function<
452  value_type(const FEPointEvaluation<n_components,
453  dim,
454  spacedim,
455  typename VectorType::value_type> &,
456  const unsigned int &)> process_quadrature_point,
458  std::vector<typename VectorType::value_type> &solution_values,
459  std::vector<
460  std::unique_ptr<FEPointEvaluation<n_components,
461  dim,
462  spacedim,
463  typename VectorType::value_type>>>
464  &evaluators)
465  {
466  if (evaluators.size() == 0)
467  evaluators.resize(dof_handler.get_fe_collection().size());
468 
469  typename DoFHandler<dim>::active_cell_iterator cell = {
470  &cache.get_triangulation(),
471  cell_data.cells[i].first,
472  cell_data.cells[i].second,
473  &dof_handler};
474 
475  const ArrayView<const Point<dim>> unit_points(
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]);
480 
481  solution_values.resize(
482  dof_handler.get_fe(cell->active_fe_index()).n_dofs_per_cell());
483  cell->get_dof_values(vector,
484  solution_values.begin(),
485  solution_values.end());
486 
487  if (evaluators[cell->active_fe_index()] == nullptr)
488  evaluators[cell->active_fe_index()] =
489  std::make_unique<FEPointEvaluation<n_components,
490  dim,
491  spacedim,
492  typename VectorType::value_type>>(
493  cache.get_mapping(),
494  cell->get_fe(),
495  update_flags,
496  first_selected_component);
497  auto &evaluator = *evaluators[cell->active_fe_index()];
498 
499  evaluator.reinit(cell, unit_points);
500  evaluator.evaluate(solution_values, evaluation_flags);
501 
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);
505  }
506 
507 
508 
509  template <int dim, int spacedim, typename Number>
510  Number
511  get_value(
513  const Vector<Number> & vector,
515  {
516  (void)tria;
518  return vector[cell->active_cell_index()];
519  }
520 
521 
522 
523  template <int dim, int spacedim, typename Number>
524  Number
525  get_value(
529  {
530  const auto distributed_tria =
531  dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(&tria);
532 
533  const bool use_distributed_path =
534  (distributed_tria == nullptr) ?
535  false :
536  (vector.get_partitioner().get() ==
537  distributed_tria->global_active_cell_index_partitioner()
538  .lock()
539  .get());
540 
541  if (use_distributed_path)
542  {
543  return vector[cell->global_active_cell_index()];
544  }
545  else
546  {
548  return vector[cell->active_cell_index()];
549  }
550  }
551 
552 
553 
554  template <typename Number, typename Number2>
555  void
556  set_value(Number &dst, const Number2 &src)
557  {
558  dst = src;
559  }
560 
561 
562 
563  template <typename Number, int rank, int dim, typename Number2>
564  void
565  set_value(Tensor<rank, dim, Number> &, const Number2 &)
566  {
567  Assert(false,
568  ExcMessage(
569  "A cell-data vector can only have a single component."));
570  }
571 
572 
573 
574  template <int n_components,
575  int dim,
576  int spacedim,
577  typename VectorType,
578  typename value_type>
579  void
580  process_cell(
581  const unsigned int i,
582  const typename Utilities::MPI::RemotePointEvaluation<dim,
583  spacedim>::CellData
584  &cell_data,
587  const VectorType & vector,
588  const UpdateFlags,
590  const unsigned int first_selected_component,
591  const std::function<
592  value_type(const FEPointEvaluation<n_components,
593  dim,
594  spacedim,
595  typename VectorType::value_type> &,
596  const unsigned int &)>,
598  std::vector<typename VectorType::value_type> &,
599  std::vector<
600  std::unique_ptr<FEPointEvaluation<n_components,
601  dim,
602  spacedim,
603  typename VectorType::value_type>>> &)
604  {
605  (void)evaluation_flags;
606  (void)first_selected_component;
607 
608  Assert(n_components == 1 && first_selected_component == 0,
609  ExcMessage(
610  "A cell-data vector can only have a single component."));
611 
612  Assert(evaluation_flags ==
614  ExcMessage("For cell-data vectors, only values can be queried."));
615 
617  &triangulation, cell_data.cells[i].first, cell_data.cells[i].second};
618 
619  const auto value = get_value(triangulation, vector, cell);
620 
621  for (unsigned int q = cell_data.reference_point_ptrs[i];
622  q < cell_data.reference_point_ptrs[i + 1];
623  ++q)
624  set_value(values[q], value);
625  }
626 
627 
628 
629  template <int n_components,
630  int dim,
631  int spacedim,
632  typename MeshType,
633  typename VectorType,
634  typename value_type>
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,
644  const UpdateFlags update_flags,
646  const std::function<
647  value_type(const FEPointEvaluation<n_components,
648  dim,
649  spacedim,
650  typename VectorType::value_type> &,
651  const unsigned int &)> process_quadrature_point)
652  {
653  Assert(cache.is_ready(),
654  ExcMessage(
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."));
658 
659  Assert(
660  &mesh.get_triangulation() == &cache.get_triangulation(),
661  ExcMessage(
662  "The provided Utilities::MPI::RemotePointEvaluation and DoFHandler "
663  "object have been set up with different Triangulation objects, "
664  "a scenario not supported!"));
665 
666  // evaluate values at points if possible
667  const auto evaluation_point_results = [&]() {
668  // helper function for accessing the global vector and interpolating
669  // the results onto the points
670  const auto evaluation_function = [&](auto & values,
671  const auto &cell_data) {
672  std::vector<typename VectorType::value_type> solution_values;
673 
674  std::vector<
675  std::unique_ptr<FEPointEvaluation<n_components,
676  dim,
677  spacedim,
678  typename VectorType::value_type>>>
679  evaluators;
680 
681  for (unsigned int i = 0; i < cell_data.cells.size(); ++i)
682  process_cell<n_components, dim, spacedim, VectorType, value_type>(
683  i,
684  cell_data,
685  cache,
686  mesh,
687  vector,
688  update_flags,
689  evaluation_flags,
690  first_selected_component,
691  process_quadrature_point,
692  values,
693  solution_values,
694  evaluators);
695  };
696 
697  std::vector<value_type> evaluation_point_results;
698  std::vector<value_type> buffer;
699 
700  cache.template evaluate_and_process<value_type>(
701  evaluation_point_results, buffer, evaluation_function);
702 
703  return evaluation_point_results;
704  }();
705 
706  if (cache.is_map_unique())
707  {
708  // each point has exactly one result (unique map)
709  return evaluation_point_results;
710  }
711  else
712  {
713  // map is not unique (multiple or no results): postprocessing is
714  // needed
715  std::vector<value_type> unique_evaluation_point_results(
716  cache.get_point_ptrs().size() - 1);
717 
718  const auto &ptr = cache.get_point_ptrs();
719 
720  for (unsigned int i = 0; i < ptr.size() - 1; ++i)
721  {
722  const auto n_entries = ptr[i + 1] - ptr[i];
723  if (n_entries == 0)
724  continue;
725 
726  unique_evaluation_point_results[i] =
727  reduce(flags,
729  evaluation_point_results.data() + ptr[i], n_entries));
730  }
731 
732  return unique_evaluation_point_results;
733  }
734  }
735  } // namespace internal
736 
737  template <int n_components,
738  template <int, int>
739  class MeshType,
740  int dim,
741  int spacedim,
742  typename VectorType>
744  (concepts::is_dealii_vector_type<VectorType> &&
745  concepts::is_triangulation_or_dof_handler<MeshType<dim, spacedim>>))
746  inline std::vector<
747  typename FEPointEvaluation<n_components,
748  dim,
749  spacedim,
750  typename VectorType::value_type>::
751  value_type> point_values(const Utilities::MPI::
752  RemotePointEvaluation<dim, spacedim> &cache,
753  const MeshType<dim, spacedim> & mesh,
754  const VectorType & vector,
755  const EvaluationFlags::EvaluationFlags flags,
756  const unsigned int first_selected_component)
757  {
758  return internal::evaluate_at_points<
759  n_components,
760  dim,
761  spacedim,
762  MeshType<dim, spacedim>,
763  VectorType,
764  typename FEPointEvaluation<n_components,
765  dim,
766  spacedim,
767  typename VectorType::value_type>::value_type>(
768  cache,
769  mesh,
770  vector,
771  flags,
772  first_selected_component,
775  [](const auto &evaluator, const auto &q) {
776  return evaluator.get_value(q);
777  });
778  }
779 
780 
781 
782  template <int n_components,
783  template <int, int>
784  class MeshType,
785  int dim,
786  int spacedim,
787  typename VectorType>
789  (concepts::is_dealii_vector_type<VectorType> &&
790  concepts::is_triangulation_or_dof_handler<MeshType<dim, spacedim>>))
791  inline std::vector<
792  typename FEPointEvaluation<n_components,
793  dim,
794  spacedim,
795  typename VectorType::value_type>::
796  gradient_type> point_gradients(const Utilities::MPI::
797  RemotePointEvaluation<dim, spacedim>
798  & cache,
799  const MeshType<dim, spacedim> &mesh,
800  const VectorType & vector,
802  flags,
803  const unsigned int
804  first_selected_component)
805  {
806  return internal::evaluate_at_points<
807  n_components,
808  dim,
809  spacedim,
810  MeshType<dim, spacedim>,
811  VectorType,
812  typename FEPointEvaluation<
813  n_components,
814  dim,
815  spacedim,
816  typename VectorType::value_type>::gradient_type>(
817  cache,
818  mesh,
819  vector,
820  flags,
821  first_selected_component,
824  [](const auto &evaluator, const unsigned &q) {
825  return evaluator.get_gradient(q);
826  });
827  }
828 
829 #endif
830 } // namespace VectorTools
831 
833 
834 #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:314
Definition: point.h:110
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:109
size_type size() const
unsigned int size() const
Definition: collection.h:264
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:474
#define DEAL_II_CXX20_REQUIRES(condition)
Definition: config.h:162
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:475
#define Assert(cond, exc)
Definition: exceptions.h:1586
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1759
static ::ExceptionBase & ExcMessage(std::string arg1)
typename ActiveSelector::active_cell_iterator active_cell_iterator
Definition: dof_handler.h:439
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