Reference documentation for deal.II version Git fabdb8cac6 2021-10-20 09:44:48 +0200
\(\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 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 
114  template <int n_components, int dim, int spacedim, typename VectorType>
115  std::vector<typename FEPointEvaluation<n_components, dim>::value_type>
116  point_values(
117  const Mapping<dim> & mapping,
118  const DoFHandler<dim, spacedim> & dof_handler,
119  const VectorType & vector,
120  const std::vector<Point<spacedim>> & evaluation_points,
123 
136  template <int n_components, int dim, int spacedim, typename VectorType>
137  std::vector<typename FEPointEvaluation<n_components, dim>::value_type>
138  point_values(
140  const DoFHandler<dim, spacedim> & dof_handler,
141  const VectorType & vector,
143 
154  template <int n_components, int dim, int spacedim, typename VectorType>
155  std::vector<typename FEPointEvaluation<n_components, dim>::gradient_type>
157  const Mapping<dim> & mapping,
158  const DoFHandler<dim, spacedim> & dof_handler,
159  const VectorType & vector,
160  const std::vector<Point<spacedim>> & evaluation_points,
163 
175  template <int n_components, int dim, int spacedim, typename VectorType>
176  std::vector<typename FEPointEvaluation<n_components, dim>::gradient_type>
179  const DoFHandler<dim, spacedim> & dof_handler,
180  const VectorType & vector,
182 
183 
184 
185  // inlined functions
186 
187 
188 #ifndef DOXYGEN
189  template <int n_components, int dim, int spacedim, typename VectorType>
190  inline std::vector<typename FEPointEvaluation<n_components, dim>::value_type>
191  point_values(const Mapping<dim> & mapping,
192  const DoFHandler<dim, spacedim> & dof_handler,
193  const VectorType & vector,
194  const std::vector<Point<spacedim>> &evaluation_points,
197  {
198  cache.reinit(evaluation_points, dof_handler.get_triangulation(), mapping);
199 
200  return point_values<n_components>(cache, dof_handler, vector, flags);
201  }
202 
203 
204 
205  template <int n_components, int dim, int spacedim, typename VectorType>
206  inline std::vector<
208  point_gradients(const Mapping<dim> & mapping,
209  const DoFHandler<dim, spacedim> & dof_handler,
210  const VectorType & vector,
211  const std::vector<Point<spacedim>> &evaluation_points,
214  {
215  cache.reinit(evaluation_points, dof_handler.get_triangulation(), mapping);
216 
217  return point_gradients<n_components>(cache, dof_handler, vector, flags);
218  }
219 
220 
221 
222  namespace internal
223  {
227  template <typename T>
228  T
230  const ArrayView<const T> & values)
231  {
232  switch (flags)
233  {
235  {
236  return std::accumulate(values.begin(), values.end(), T{}) /
237  (T(1.0) * values.size());
238  }
240  return *std::max_element(values.begin(), values.end());
242  return *std::min_element(values.begin(), values.end());
244  return values[0];
245  default:
246  Assert(false, ExcNotImplemented());
247  return values[0];
248  }
249  }
250 
251 
252 
256  template <int rank, int dim, typename Number>
259  const ArrayView<const Tensor<rank, dim, Number>> &values)
260  {
261  switch (flags)
262  {
264  {
265  return std::accumulate(values.begin(),
266  values.end(),
268  (Number(1.0) * values.size());
269  }
271  return values[0];
272  default:
273  Assert(false, ExcNotImplemented());
274  return values[0];
275  }
276  }
277 
278 
279 
280  template <int n_components,
281  int dim,
282  int spacedim,
283  typename VectorType,
284  typename value_type>
285  inline std::vector<value_type>
286  evaluate_at_points(
288  const DoFHandler<dim, spacedim> & dof_handler,
289  const VectorType & vector,
291  const UpdateFlags update_flags,
293  const std::function<
294  value_type(const FEPointEvaluation<n_components, dim> &,
295  const unsigned int &)> process_quadrature_point)
296  {
297  Assert(cache.is_ready(),
298  ExcMessage(
299  "Utilities::MPI::RemotePointEvaluation is not ready yet! "
300  "Please call Utilities::MPI::RemotePointEvaluation::reinit() "
301  "yourself or another function that does this for you."));
302 
303  Assert(
304  &dof_handler.get_triangulation() == &cache.get_triangulation(),
305  ExcMessage(
306  "The provided Utilities::MPI::RemotePointEvaluation and DoFHandler "
307  "object have been set up with different Triangulation objects, "
308  "a scenario not supported!"));
309 
310  // evaluate values at points if possible
311  const auto evaluation_point_results = [&]() {
312  // helper function for accessing the global vector and interpolating
313  // the results onto the points
314  const auto evaluation_function = [&](auto & values,
315  const auto &cell_data) {
316  std::vector<typename VectorType::value_type> solution_values;
317 
318  std::vector<std::unique_ptr<FEPointEvaluation<n_components, dim>>>
319  evaluators(dof_handler.get_fe_collection().size());
320 
321  for (unsigned int i = 0; i < cell_data.cells.size(); ++i)
322  {
323  typename DoFHandler<dim>::active_cell_iterator cell = {
324  &cache.get_triangulation(),
325  cell_data.cells[i].first,
326  cell_data.cells[i].second,
327  &dof_handler};
328 
329  const ArrayView<const Point<dim>> unit_points(
330  cell_data.reference_point_values.data() +
331  cell_data.reference_point_ptrs[i],
332  cell_data.reference_point_ptrs[i + 1] -
333  cell_data.reference_point_ptrs[i]);
334 
335  solution_values.resize(
336  dof_handler.get_fe(cell->active_fe_index()).n_dofs_per_cell());
337  cell->get_dof_values(vector,
338  solution_values.begin(),
339  solution_values.end());
340 
341  if (evaluators[cell->active_fe_index()] == nullptr)
342  evaluators[cell->active_fe_index()] =
343  std::make_unique<FEPointEvaluation<n_components, dim>>(
344  cache.get_mapping(), cell->get_fe(), update_flags);
345  auto &evaluator = *evaluators[cell->active_fe_index()];
346 
347  evaluator.reinit(cell, unit_points);
348  evaluator.evaluate(solution_values, evaluation_flags);
349 
350  for (unsigned int q = 0; q < unit_points.size(); ++q)
351  values[q + cell_data.reference_point_ptrs[i]] =
352  process_quadrature_point(evaluator, q);
353  }
354  };
355 
356  std::vector<value_type> evaluation_point_results;
357  std::vector<value_type> buffer;
358 
359  cache.template evaluate_and_process<value_type>(
360  evaluation_point_results, buffer, evaluation_function);
361 
362  return evaluation_point_results;
363  }();
364 
365  if (cache.is_map_unique())
366  {
367  // each point has exactly one result (unique map)
368  return evaluation_point_results;
369  }
370  else
371  {
372  // map is not unique (multiple or no results): postprocessing is
373  // needed
374  std::vector<value_type> unique_evaluation_point_results(
375  cache.get_point_ptrs().size() - 1);
376 
377  const auto &ptr = cache.get_point_ptrs();
378 
379  for (unsigned int i = 0; i < ptr.size() - 1; ++i)
380  {
381  const auto n_entries = ptr[i + 1] - ptr[i];
382  if (n_entries == 0)
383  continue;
384 
385  unique_evaluation_point_results[i] =
386  reduce(flags,
388  evaluation_point_results.data() + ptr[i], n_entries));
389  }
390 
391  return unique_evaluation_point_results;
392  }
393  }
394  } // namespace internal
395 
396  template <int n_components, int dim, int spacedim, typename VectorType>
397  inline std::vector<typename FEPointEvaluation<n_components, dim>::value_type>
398  point_values(
400  const DoFHandler<dim, spacedim> & dof_handler,
401  const VectorType & vector,
403  {
404  return internal::evaluate_at_points<
405  n_components,
406  dim,
407  spacedim,
408  VectorType,
410  cache,
411  dof_handler,
412  vector,
413  flags,
416  [](const auto &evaluator, const auto &q) {
417  return evaluator.get_value(q);
418  });
419  }
420 
421  template <int n_components, int dim, int spacedim, typename VectorType>
422  inline std::vector<
426  const DoFHandler<dim, spacedim> & dof_handler,
427  const VectorType & vector,
429  {
430  return internal::evaluate_at_points<
431  n_components,
432  dim,
433  spacedim,
434  VectorType,
436  cache,
437  dof_handler,
438  vector,
439  flags,
442  [](const auto &evaluator, const auto &q) {
443  return evaluator.get_gradient(q);
444  });
445  }
446 
447 #endif
448 } // namespace VectorTools
449 
451 
452 #endif // dealii_vector_tools_boundary_h
Shape function values.
The namespace for the EvaluationFlags enum.
const Triangulation< dim, spacedim > & get_triangulation() const
void reinit(const std::vector< Point< spacedim >> &points, const Triangulation< dim, spacedim > &tria, const Mapping< dim, spacedim > &mapping)
iterator end() const
Definition: array_view.h:593
const std::vector< unsigned int > & get_point_ptrs() const
const FiniteElement< dim, spacedim > & get_fe(const unsigned int index=0) const
std::vector< typename FEPointEvaluation< n_components, dim >::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)
std::size_t size() const
Definition: array_view.h:575
static ::ExceptionBase & ExcMessage(std::string arg1)
static const char T
#define Assert(cond, exc)
Definition: exceptions.h:1461
UpdateFlags
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:401
std::vector< typename FEPointEvaluation< n_components, dim >::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)
const hp::FECollection< dim, spacedim > & get_fe_collection() const
const Mapping< dim, spacedim > & get_mapping() const
typename internal::FEPointEvaluation::EvaluatorTypeTraits< dim, n_components, Number >::value_type value_type
void reinit(const Triangulation< dim, spacedim > &tria)
const Triangulation< dim, spacedim > & get_triangulation() const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:400
Shape function gradients.
static ::ExceptionBase & ExcNotImplemented()
typename internal::FEPointEvaluation::EvaluatorTypeTraits< dim, n_components, Number >::gradient_type gradient_type
iterator begin() const
Definition: array_view.h:584
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)