Reference documentation for deal.II version GIT 6dd81a2340 2022-11-29 05:55: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>
140  std::vector<
141  typename FEPointEvaluation<n_components,
142  dim,
143  spacedim,
144  typename VectorType::value_type>::value_type>
146  const Mapping<dim> & mapping,
147  const MeshType<dim, spacedim> & mesh,
148  const VectorType & vector,
149  const std::vector<Point<spacedim>> & evaluation_points,
152  const unsigned int first_selected_component = 0);
153 
166  template <int n_components,
167  template <int, int>
168  class MeshType,
169  int dim,
170  int spacedim,
171  typename VectorType>
172  std::vector<
173  typename FEPointEvaluation<n_components,
174  dim,
175  spacedim,
176  typename VectorType::value_type>::value_type>
179  const MeshType<dim, spacedim> & mesh,
180  const VectorType & vector,
182  const unsigned int first_selected_component = 0);
183 
194  template <int n_components,
195  template <int, int>
196  class MeshType,
197  int dim,
198  int spacedim,
199  typename VectorType>
200  std::vector<
201  typename FEPointEvaluation<n_components,
202  dim,
203  spacedim,
204  typename VectorType::value_type>::gradient_type>
206  const Mapping<dim> & mapping,
207  const MeshType<dim, spacedim> & mesh,
208  const VectorType & vector,
209  const std::vector<Point<spacedim>> & evaluation_points,
212  const unsigned int first_selected_component = 0);
213 
225  template <int n_components,
226  template <int, int>
227  class MeshType,
228  int dim,
229  int spacedim,
230  typename VectorType>
231  std::vector<
232  typename FEPointEvaluation<n_components,
233  dim,
234  spacedim,
235  typename VectorType::value_type>::gradient_type>
238  const MeshType<dim, spacedim> & mesh,
239  const VectorType & vector,
241  const unsigned int first_selected_component = 0);
242 
243 
244 
245  // inlined functions
246 
247 
248 #ifndef DOXYGEN
249  template <int n_components,
250  template <int, int>
251  class MeshType,
252  int dim,
253  int spacedim,
254  typename VectorType>
255  inline std::vector<
256  typename FEPointEvaluation<n_components,
257  dim,
258  spacedim,
259  typename VectorType::value_type>::value_type>
260  point_values(const Mapping<dim> & mapping,
261  const MeshType<dim, spacedim> & mesh,
262  const VectorType & vector,
263  const std::vector<Point<spacedim>> &evaluation_points,
266  const unsigned int first_selected_component)
267  {
268  cache.reinit(evaluation_points, mesh.get_triangulation(), mapping);
269 
270  return point_values<n_components>(
271  cache, mesh, vector, flags, first_selected_component);
272  }
273 
274 
275 
276  template <int n_components,
277  template <int, int>
278  class MeshType,
279  int dim,
280  int spacedim,
281  typename VectorType>
282  inline std::vector<
283  typename FEPointEvaluation<n_components,
284  dim,
285  spacedim,
286  typename VectorType::value_type>::gradient_type>
287  point_gradients(const Mapping<dim> & mapping,
288  const MeshType<dim, spacedim> & mesh,
289  const VectorType & vector,
290  const std::vector<Point<spacedim>> &evaluation_points,
293  const unsigned int first_selected_component)
294  {
295  cache.reinit(evaluation_points, mesh.get_triangulation(), mapping);
296 
297  return point_gradients<n_components>(
298  cache, mesh, vector, flags, first_selected_component);
299  }
300 
301 
302 
303  namespace internal
304  {
308  template <typename T>
309  T
311  const ArrayView<const T> & values)
312  {
313  switch (flags)
314  {
316  {
317  return std::accumulate(values.begin(), values.end(), T{}) /
318  (T(1.0) * values.size());
319  }
321  return *std::max_element(values.begin(), values.end());
323  return *std::min_element(values.begin(), values.end());
325  return values[0];
326  default:
327  Assert(false, ExcNotImplemented());
328  return values[0];
329  }
330  }
331 
332 
333 
337  template <int rank, int dim, typename Number>
341  {
342  switch (flags)
343  {
345  {
346  return std::accumulate(values.begin(),
347  values.end(),
349  (Number(1.0) * values.size());
350  }
352  return values[0];
353  default:
354  Assert(false, ExcNotImplemented());
355  return values[0];
356  }
357  }
358 
359 
360 
365  template <int n_components, int rank, int dim, typename Number>
367  reduce(
369  const ArrayView<const Tensor<1, n_components, Tensor<rank, dim, Number>>>
370  &values)
371  {
372  switch (flags)
373  {
375  {
377 
378  for (unsigned int j = 0; j < values.size(); ++j)
379  for (unsigned int i = 0; i < n_components; ++i)
380  temp[i] = temp[i] + values[j][i];
381 
382  for (unsigned int i = 0; i < n_components; ++i)
383  temp[i] /= Number(values.size());
384 
385  return temp;
386  }
388  return values[0];
389  default:
390  Assert(false, ExcNotImplemented());
391  return values[0];
392  }
393  }
394 
395 
396 
397  template <int n_components,
398  int dim,
399  int spacedim,
400  typename VectorType,
401  typename value_type>
402  void
403  process_cell(
404  const unsigned int i,
405  const typename Utilities::MPI::RemotePointEvaluation<dim,
406  spacedim>::CellData
407  & cell_data,
409  const DoFHandler<dim, spacedim> & dof_handler,
410  const VectorType & vector,
411  const UpdateFlags update_flags,
413  const unsigned int first_selected_component,
414  const std::function<
415  value_type(const FEPointEvaluation<n_components,
416  dim,
417  spacedim,
418  typename VectorType::value_type> &,
419  const unsigned int &)> process_quadrature_point,
421  std::vector<typename VectorType::value_type> &solution_values,
422  std::vector<
423  std::unique_ptr<FEPointEvaluation<n_components,
424  dim,
425  spacedim,
426  typename VectorType::value_type>>>
427  &evaluators)
428  {
429  if (evaluators.size() == 0)
430  evaluators.resize(dof_handler.get_fe_collection().size());
431 
432  typename DoFHandler<dim>::active_cell_iterator cell = {
433  &cache.get_triangulation(),
434  cell_data.cells[i].first,
435  cell_data.cells[i].second,
436  &dof_handler};
437 
438  const ArrayView<const Point<dim>> unit_points(
439  cell_data.reference_point_values.data() +
440  cell_data.reference_point_ptrs[i],
441  cell_data.reference_point_ptrs[i + 1] -
442  cell_data.reference_point_ptrs[i]);
443 
444  solution_values.resize(
445  dof_handler.get_fe(cell->active_fe_index()).n_dofs_per_cell());
446  cell->get_dof_values(vector,
447  solution_values.begin(),
448  solution_values.end());
449 
450  if (evaluators[cell->active_fe_index()] == nullptr)
451  evaluators[cell->active_fe_index()] =
452  std::make_unique<FEPointEvaluation<n_components,
453  dim,
454  spacedim,
455  typename VectorType::value_type>>(
456  cache.get_mapping(),
457  cell->get_fe(),
458  update_flags,
459  first_selected_component);
460  auto &evaluator = *evaluators[cell->active_fe_index()];
461 
462  evaluator.reinit(cell, unit_points);
463  evaluator.evaluate(solution_values, evaluation_flags);
464 
465  for (unsigned int q = 0; q < unit_points.size(); ++q)
466  values[q + cell_data.reference_point_ptrs[i]] =
467  process_quadrature_point(evaluator, q);
468  }
469 
470 
471 
472  template <int dim, int spacedim, typename Number>
473  Number
474  get_value(
476  const Vector<Number> & vector,
478  {
479  (void)tria;
481  return vector[cell->active_cell_index()];
482  }
483 
484 
485 
486  template <int dim, int spacedim, typename Number>
487  Number
488  get_value(
492  {
493  const auto distributed_tria =
494  dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(&tria);
495 
496  const bool use_distributed_path =
497  (distributed_tria == nullptr) ?
498  false :
499  (vector.get_partitioner().get() ==
500  distributed_tria->global_active_cell_index_partitioner()
501  .lock()
502  .get());
503 
504  if (use_distributed_path)
505  {
506  return vector[cell->global_active_cell_index()];
507  }
508  else
509  {
511  return vector[cell->active_cell_index()];
512  }
513  }
514 
515 
516 
517  template <typename Number, typename Number2>
518  void
519  set_value(Number &dst, const Number2 &src)
520  {
521  dst = src;
522  }
523 
524 
525 
526  template <typename Number, int rank, int dim, typename Number2>
527  void
528  set_value(Tensor<rank, dim, Number> &, const Number2 &)
529  {
530  Assert(false,
531  ExcMessage(
532  "A cell-data vector can only have a single component."));
533  }
534 
535 
536 
537  template <int n_components,
538  int dim,
539  int spacedim,
540  typename VectorType,
541  typename value_type>
542  void
543  process_cell(
544  const unsigned int i,
545  const typename Utilities::MPI::RemotePointEvaluation<dim,
546  spacedim>::CellData
547  &cell_data,
550  const VectorType & vector,
551  const UpdateFlags,
553  const unsigned int first_selected_component,
554  const std::function<
555  value_type(const FEPointEvaluation<n_components,
556  dim,
557  spacedim,
558  typename VectorType::value_type> &,
559  const unsigned int &)>,
561  std::vector<typename VectorType::value_type> &,
562  std::vector<
563  std::unique_ptr<FEPointEvaluation<n_components,
564  dim,
565  spacedim,
566  typename VectorType::value_type>>> &)
567  {
568  (void)evaluation_flags;
569  (void)first_selected_component;
570 
571  Assert(n_components == 1 && first_selected_component == 0,
572  ExcMessage(
573  "A cell-data vector can only have a single component."));
574 
575  Assert(evaluation_flags ==
577  ExcMessage("For cell-data vectors, only values can be queried."));
578 
580  &triangulation, cell_data.cells[i].first, cell_data.cells[i].second};
581 
582  const auto value = get_value(triangulation, vector, cell);
583 
584  for (unsigned int q = cell_data.reference_point_ptrs[i];
585  q < cell_data.reference_point_ptrs[i + 1];
586  ++q)
587  set_value(values[q], value);
588  }
589 
590 
591 
592  template <int n_components,
593  int dim,
594  int spacedim,
595  typename MeshType,
596  typename VectorType,
597  typename value_type>
598  inline std::vector<value_type>
599  evaluate_at_points(
601  const MeshType & mesh,
602  const VectorType & vector,
604  const unsigned int first_selected_component,
605  const UpdateFlags update_flags,
607  const std::function<
608  value_type(const FEPointEvaluation<n_components,
609  dim,
610  spacedim,
611  typename VectorType::value_type> &,
612  const unsigned int &)> process_quadrature_point)
613  {
614  Assert(cache.is_ready(),
615  ExcMessage(
616  "Utilities::MPI::RemotePointEvaluation is not ready yet! "
617  "Please call Utilities::MPI::RemotePointEvaluation::reinit() "
618  "yourself or another function that does this for you."));
619 
620  Assert(
621  &mesh.get_triangulation() == &cache.get_triangulation(),
622  ExcMessage(
623  "The provided Utilities::MPI::RemotePointEvaluation and DoFHandler "
624  "object have been set up with different Triangulation objects, "
625  "a scenario not supported!"));
626 
627  // evaluate values at points if possible
628  const auto evaluation_point_results = [&]() {
629  // helper function for accessing the global vector and interpolating
630  // the results onto the points
631  const auto evaluation_function = [&](auto & values,
632  const auto &cell_data) {
633  std::vector<typename VectorType::value_type> solution_values;
634 
635  std::vector<
636  std::unique_ptr<FEPointEvaluation<n_components,
637  dim,
638  spacedim,
639  typename VectorType::value_type>>>
640  evaluators;
641 
642  for (unsigned int i = 0; i < cell_data.cells.size(); ++i)
643  process_cell<n_components, dim, spacedim, VectorType, value_type>(
644  i,
645  cell_data,
646  cache,
647  mesh,
648  vector,
649  update_flags,
650  evaluation_flags,
651  first_selected_component,
652  process_quadrature_point,
653  values,
654  solution_values,
655  evaluators);
656  };
657 
658  std::vector<value_type> evaluation_point_results;
659  std::vector<value_type> buffer;
660 
661  cache.template evaluate_and_process<value_type>(
662  evaluation_point_results, buffer, evaluation_function);
663 
664  return evaluation_point_results;
665  }();
666 
667  if (cache.is_map_unique())
668  {
669  // each point has exactly one result (unique map)
670  return evaluation_point_results;
671  }
672  else
673  {
674  // map is not unique (multiple or no results): postprocessing is
675  // needed
676  std::vector<value_type> unique_evaluation_point_results(
677  cache.get_point_ptrs().size() - 1);
678 
679  const auto &ptr = cache.get_point_ptrs();
680 
681  for (unsigned int i = 0; i < ptr.size() - 1; ++i)
682  {
683  const auto n_entries = ptr[i + 1] - ptr[i];
684  if (n_entries == 0)
685  continue;
686 
687  unique_evaluation_point_results[i] =
688  reduce(flags,
690  evaluation_point_results.data() + ptr[i], n_entries));
691  }
692 
693  return unique_evaluation_point_results;
694  }
695  }
696  } // namespace internal
697 
698  template <int n_components,
699  template <int, int>
700  class MeshType,
701  int dim,
702  int spacedim,
703  typename VectorType>
704  inline std::vector<
705  typename FEPointEvaluation<n_components,
706  dim,
707  spacedim,
708  typename VectorType::value_type>::value_type>
709  point_values(
711  const MeshType<dim, spacedim> & mesh,
712  const VectorType & vector,
714  const unsigned int first_selected_component)
715  {
716  return internal::evaluate_at_points<
717  n_components,
718  dim,
719  spacedim,
720  MeshType<dim, spacedim>,
721  VectorType,
722  typename FEPointEvaluation<n_components,
723  dim,
724  spacedim,
725  typename VectorType::value_type>::value_type>(
726  cache,
727  mesh,
728  vector,
729  flags,
730  first_selected_component,
733  [](const auto &evaluator, const auto &q) {
734  return evaluator.get_value(q);
735  });
736  }
737 
738  template <int n_components,
739  template <int, int>
740  class MeshType,
741  int dim,
742  int spacedim,
743  typename VectorType>
744  inline std::vector<
745  typename FEPointEvaluation<n_components,
746  dim,
747  spacedim,
748  typename VectorType::value_type>::gradient_type>
751  const MeshType<dim, spacedim> & mesh,
752  const VectorType & vector,
754  const unsigned int first_selected_component)
755  {
756  return internal::evaluate_at_points<
757  n_components,
758  dim,
759  spacedim,
760  MeshType<dim, spacedim>,
761  VectorType,
762  typename FEPointEvaluation<
763  n_components,
764  dim,
765  spacedim,
766  typename VectorType::value_type>::gradient_type>(
767  cache,
768  mesh,
769  vector,
770  flags,
771  first_selected_component,
774  [](const auto &evaluator, const unsigned &q) {
775  return evaluator.get_gradient(q);
776  });
777  }
778 
779 #endif
780 } // namespace VectorTools
781 
783 
784 #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:311
Definition: tensor.h:503
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
void reinit(const std::vector< Point< spacedim >> &points, const Triangulation< dim, spacedim > &tria, const Mapping< dim, spacedim > &mapping)
Definition: vector.h:109
size_type size() const
unsigned int size() const
Definition: collection.h:264
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:458
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:459
#define Assert(cond, exc)
Definition: exceptions.h:1501
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1695
static ::ExceptionBase & ExcMessage(std::string arg1)
typename ActiveSelector::active_cell_iterator active_cell_iterator
Definition: dof_handler.h:437
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)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
const ::Triangulation< dim, spacedim > & tria