Reference documentation for deal.II version GIT c99bf6dfd7 2022-05-28 07: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 - 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 
28 #include <map>
29 
31 
32 namespace VectorTools
33 {
37  namespace EvaluationFlags
38  {
43  {
47  avg = 0,
53  max = 1,
59  min = 2,
63  insert = 3
64  };
65  } // namespace EvaluationFlags
66 
119  template <int n_components, int dim, int spacedim, typename VectorType>
120  std::vector<
121  typename FEPointEvaluation<n_components,
122  dim,
123  spacedim,
124  typename VectorType::value_type>::value_type>
126  const Mapping<dim> & mapping,
127  const DoFHandler<dim, spacedim> & dof_handler,
128  const VectorType & vector,
129  const std::vector<Point<spacedim>> & evaluation_points,
132 
145  template <int n_components, int dim, int spacedim, typename VectorType>
146  std::vector<
147  typename FEPointEvaluation<n_components,
148  dim,
149  spacedim,
150  typename VectorType::value_type>::value_type>
153  const DoFHandler<dim, spacedim> & dof_handler,
154  const VectorType & vector,
156 
167  template <int n_components, int dim, int spacedim, typename VectorType>
168  std::vector<
169  typename FEPointEvaluation<n_components,
170  dim,
171  spacedim,
172  typename VectorType::value_type>::gradient_type>
174  const Mapping<dim> & mapping,
175  const DoFHandler<dim, spacedim> & dof_handler,
176  const VectorType & vector,
177  const std::vector<Point<spacedim>> & evaluation_points,
180 
192  template <int n_components, int dim, int spacedim, typename VectorType>
193  std::vector<
194  typename FEPointEvaluation<n_components,
195  dim,
196  spacedim,
197  typename VectorType::value_type>::gradient_type>
200  const DoFHandler<dim, spacedim> & dof_handler,
201  const VectorType & vector,
203 
204 
205 
206  // inlined functions
207 
208 
209 #ifndef DOXYGEN
210  template <int n_components, int dim, int spacedim, typename VectorType>
211  inline std::vector<
212  typename FEPointEvaluation<n_components,
213  dim,
214  spacedim,
215  typename VectorType::value_type>::value_type>
216  point_values(const Mapping<dim> & mapping,
217  const DoFHandler<dim, spacedim> & dof_handler,
218  const VectorType & vector,
219  const std::vector<Point<spacedim>> &evaluation_points,
222  {
223  cache.reinit(evaluation_points, dof_handler.get_triangulation(), mapping);
224 
225  return point_values<n_components>(cache, dof_handler, vector, flags);
226  }
227 
228 
229 
230  template <int n_components, int dim, int spacedim, typename VectorType>
231  inline std::vector<
232  typename FEPointEvaluation<n_components,
233  dim,
234  spacedim,
235  typename VectorType::value_type>::gradient_type>
236  point_gradients(const Mapping<dim> & mapping,
237  const DoFHandler<dim, spacedim> & dof_handler,
238  const VectorType & vector,
239  const std::vector<Point<spacedim>> &evaluation_points,
242  {
243  cache.reinit(evaluation_points, dof_handler.get_triangulation(), mapping);
244 
245  return point_gradients<n_components>(cache, dof_handler, vector, flags);
246  }
247 
248 
249 
250  namespace internal
251  {
255  template <typename T>
256  T
258  const ArrayView<const T> & values)
259  {
260  switch (flags)
261  {
263  {
264  return std::accumulate(values.begin(), values.end(), T{}) /
265  (T(1.0) * values.size());
266  }
268  return *std::max_element(values.begin(), values.end());
270  return *std::min_element(values.begin(), values.end());
272  return values[0];
273  default:
274  Assert(false, ExcNotImplemented());
275  return values[0];
276  }
277  }
278 
279 
280 
284  template <int rank, int dim, typename Number>
288  {
289  switch (flags)
290  {
292  {
293  return std::accumulate(values.begin(),
294  values.end(),
296  (Number(1.0) * values.size());
297  }
299  return values[0];
300  default:
301  Assert(false, ExcNotImplemented());
302  return values[0];
303  }
304  }
305 
306 
307 
308  template <int n_components,
309  int dim,
310  int spacedim,
311  typename VectorType,
312  typename value_type>
313  inline std::vector<value_type>
314  evaluate_at_points(
316  const DoFHandler<dim, spacedim> & dof_handler,
317  const VectorType & vector,
319  const UpdateFlags update_flags,
321  const std::function<
322  value_type(const FEPointEvaluation<n_components,
323  dim,
324  spacedim,
325  typename VectorType::value_type> &,
326  const unsigned int &)> process_quadrature_point)
327  {
328  Assert(cache.is_ready(),
329  ExcMessage(
330  "Utilities::MPI::RemotePointEvaluation is not ready yet! "
331  "Please call Utilities::MPI::RemotePointEvaluation::reinit() "
332  "yourself or another function that does this for you."));
333 
334  Assert(
335  &dof_handler.get_triangulation() == &cache.get_triangulation(),
336  ExcMessage(
337  "The provided Utilities::MPI::RemotePointEvaluation and DoFHandler "
338  "object have been set up with different Triangulation objects, "
339  "a scenario not supported!"));
340 
341  // evaluate values at points if possible
342  const auto evaluation_point_results = [&]() {
343  // helper function for accessing the global vector and interpolating
344  // the results onto the points
345  const auto evaluation_function = [&](auto & values,
346  const auto &cell_data) {
347  std::vector<typename VectorType::value_type> solution_values;
348 
349  std::vector<
350  std::unique_ptr<FEPointEvaluation<n_components,
351  dim,
352  spacedim,
353  typename VectorType::value_type>>>
354  evaluators(dof_handler.get_fe_collection().size());
355 
356  for (unsigned int i = 0; i < cell_data.cells.size(); ++i)
357  {
358  typename DoFHandler<dim>::active_cell_iterator cell = {
359  &cache.get_triangulation(),
360  cell_data.cells[i].first,
361  cell_data.cells[i].second,
362  &dof_handler};
363 
364  const ArrayView<const Point<dim>> unit_points(
365  cell_data.reference_point_values.data() +
366  cell_data.reference_point_ptrs[i],
367  cell_data.reference_point_ptrs[i + 1] -
368  cell_data.reference_point_ptrs[i]);
369 
370  solution_values.resize(
371  dof_handler.get_fe(cell->active_fe_index()).n_dofs_per_cell());
372  cell->get_dof_values(vector,
373  solution_values.begin(),
374  solution_values.end());
375 
376  if (evaluators[cell->active_fe_index()] == nullptr)
377  evaluators[cell->active_fe_index()] = std::make_unique<
378  FEPointEvaluation<n_components,
379  dim,
380  spacedim,
381  typename VectorType::value_type>>(
382  cache.get_mapping(), cell->get_fe(), update_flags);
383  auto &evaluator = *evaluators[cell->active_fe_index()];
384 
385  evaluator.reinit(cell, unit_points);
386  evaluator.evaluate(solution_values, evaluation_flags);
387 
388  for (unsigned int q = 0; q < unit_points.size(); ++q)
389  values[q + cell_data.reference_point_ptrs[i]] =
390  process_quadrature_point(evaluator, q);
391  }
392  };
393 
394  std::vector<value_type> evaluation_point_results;
395  std::vector<value_type> buffer;
396 
397  cache.template evaluate_and_process<value_type>(
398  evaluation_point_results, buffer, evaluation_function);
399 
400  return evaluation_point_results;
401  }();
402 
403  if (cache.is_map_unique())
404  {
405  // each point has exactly one result (unique map)
406  return evaluation_point_results;
407  }
408  else
409  {
410  // map is not unique (multiple or no results): postprocessing is
411  // needed
412  std::vector<value_type> unique_evaluation_point_results(
413  cache.get_point_ptrs().size() - 1);
414 
415  const auto &ptr = cache.get_point_ptrs();
416 
417  for (unsigned int i = 0; i < ptr.size() - 1; ++i)
418  {
419  const auto n_entries = ptr[i + 1] - ptr[i];
420  if (n_entries == 0)
421  continue;
422 
423  unique_evaluation_point_results[i] =
424  reduce(flags,
426  evaluation_point_results.data() + ptr[i], n_entries));
427  }
428 
429  return unique_evaluation_point_results;
430  }
431  }
432  } // namespace internal
433 
434  template <int n_components, int dim, int spacedim, typename VectorType>
435  inline std::vector<
436  typename FEPointEvaluation<n_components,
437  dim,
438  spacedim,
439  typename VectorType::value_type>::value_type>
440  point_values(
442  const DoFHandler<dim, spacedim> & dof_handler,
443  const VectorType & vector,
445  {
446  return internal::evaluate_at_points<
447  n_components,
448  dim,
449  spacedim,
450  VectorType,
451  typename FEPointEvaluation<n_components,
452  dim,
453  spacedim,
454  typename VectorType::value_type>::value_type>(
455  cache,
456  dof_handler,
457  vector,
458  flags,
461  [](const auto &evaluator, const auto &q) {
462  return evaluator.get_value(q);
463  });
464  }
465 
466  template <int n_components, int dim, int spacedim, typename VectorType>
467  inline std::vector<
468  typename FEPointEvaluation<n_components,
469  dim,
470  spacedim,
471  typename VectorType::value_type>::gradient_type>
474  const DoFHandler<dim, spacedim> & dof_handler,
475  const VectorType & vector,
477  {
478  return internal::evaluate_at_points<
479  n_components,
480  dim,
481  spacedim,
482  VectorType,
483  typename FEPointEvaluation<
484  n_components,
485  dim,
486  spacedim,
487  typename VectorType::value_type>::gradient_type>(
488  cache,
489  dof_handler,
490  vector,
491  flags,
494  [](const auto &evaluator, const auto &q) {
495  return evaluator.get_gradient(q);
496  });
497  }
498 
499 #endif
500 } // namespace VectorTools
501 
503 
504 #endif // dealii_vector_tools_boundary_h
const Triangulation< dim, spacedim > & get_triangulation() const
const FiniteElement< dim, spacedim > & get_fe(const unsigned int 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
Abstract base class for mapping classes.
Definition: mapping.h:311
Definition: tensor.h:503
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)
unsigned int size() const
Definition: collection.h:263
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:416
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:417
UpdateFlags
@ update_values
Shape function values.
@ update_gradients
Shape function gradients.
#define Assert(cond, exc)
Definition: exceptions.h:1473
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcMessage(std::string arg1)
typename ActiveSelector::active_cell_iterator active_cell_iterator
Definition: dof_handler.h:438
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 DoFHandler< dim, spacedim > &dof_handler, const VectorType &vector, const std::vector< Point< spacedim >> &evaluation_points, Utilities::MPI::RemotePointEvaluation< dim, spacedim > &cache, const EvaluationFlags::EvaluationFlags flags=EvaluationFlags::avg)
std::vector< typename FEPointEvaluation< n_components, dim, spacedim, typename VectorType::value_type >::value_type > point_values(const Mapping< dim > &mapping, const DoFHandler< dim, spacedim > &dof_handler, const VectorType &vector, const std::vector< Point< spacedim >> &evaluation_points, Utilities::MPI::RemotePointEvaluation< dim, spacedim > &cache, const EvaluationFlags::EvaluationFlags flags=EvaluationFlags::avg)