Reference documentation for deal.II version Git 42dd89c428 2021-07-27 06:40:55 +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 
74  template <int n_components, int dim, int spacedim, typename VectorType>
75  std::vector<typename FEPointEvaluation<n_components, dim>::value_type>
77  const Mapping<dim> & mapping,
78  const DoFHandler<dim, spacedim> & dof_handler,
79  const VectorType & vector,
80  const std::vector<Point<spacedim>> & evaluation_points,
83 
95  template <int n_components, int dim, int spacedim, typename VectorType>
96  std::vector<typename FEPointEvaluation<n_components, dim>::value_type>
99  const DoFHandler<dim, spacedim> & dof_handler,
100  const VectorType & vector,
102 
107  template <int n_components, int dim, int spacedim, typename VectorType>
108  std::vector<typename FEPointEvaluation<n_components, dim>::gradient_type>
110  const Mapping<dim> & mapping,
111  const DoFHandler<dim, spacedim> & dof_handler,
112  const VectorType & vector,
113  const std::vector<Point<spacedim>> & evaluation_points,
116 
125  template <int n_components, int dim, int spacedim, typename VectorType>
126  std::vector<typename FEPointEvaluation<n_components, dim>::gradient_type>
129  const DoFHandler<dim, spacedim> & dof_handler,
130  const VectorType & vector,
132 
133 
134 
135  // inlined functions
136 
137 
138 #ifndef DOXYGEN
139  template <int n_components, int dim, int spacedim, typename VectorType>
140  inline std::vector<typename FEPointEvaluation<n_components, dim>::value_type>
141  point_values(const Mapping<dim> & mapping,
142  const DoFHandler<dim, spacedim> & dof_handler,
143  const VectorType & vector,
144  const std::vector<Point<spacedim>> &evaluation_points,
147  {
148  cache.reinit(evaluation_points, dof_handler.get_triangulation(), mapping);
149 
150  return point_values<n_components>(cache, dof_handler, vector, flags);
151  }
152 
153 
154 
155  template <int n_components, int dim, int spacedim, typename VectorType>
156  inline std::vector<
158  point_gradients(const Mapping<dim> & mapping,
159  const DoFHandler<dim, spacedim> & dof_handler,
160  const VectorType & vector,
161  const std::vector<Point<spacedim>> &evaluation_points,
164  {
165  cache.reinit(evaluation_points, dof_handler.get_triangulation(), mapping);
166 
167  return point_gradients<n_components>(cache, dof_handler, vector, flags);
168  }
169 
170 
171 
172  namespace internal
173  {
177  template <typename T>
178  T
180  const ArrayView<const T> & values)
181  {
182  switch (flags)
183  {
185  {
186  return std::accumulate(values.begin(), values.end(), T{}) /
187  (T(1.0) * values.size());
188  }
190  return *std::max_element(values.begin(), values.end());
192  return *std::min_element(values.begin(), values.end());
194  return values[0];
195  default:
196  Assert(false, ExcNotImplemented());
197  return values[0];
198  }
199  }
200 
201 
202 
206  template <int rank, int dim, typename Number>
209  const ArrayView<const Tensor<rank, dim, Number>> &values)
210  {
211  switch (flags)
212  {
214  {
215  return std::accumulate(values.begin(),
216  values.end(),
218  (Number(1.0) * values.size());
219  }
221  return values[0];
222  default:
223  Assert(false, ExcNotImplemented());
224  return values[0];
225  }
226  }
227 
228 
229 
230  template <int n_components,
231  int dim,
232  int spacedim,
233  typename VectorType,
234  typename value_type>
235  inline std::vector<value_type>
236  evaluate_at_points(
238  const DoFHandler<dim, spacedim> & dof_handler,
239  const VectorType & vector,
241  const UpdateFlags update_flags,
243  const std::function<
244  value_type(const FEPointEvaluation<n_components, dim> &,
245  const unsigned int &)> process_quadrature_point)
246  {
247  Assert(cache.is_ready(),
248  ExcMessage(
249  "Utilities::MPI::RemotePointEvaluation is not ready yet! "
250  "Please call Utilities::MPI::RemotePointEvaluation::reinit() "
251  "yourself or another function that does this for you."));
252 
253  Assert(
254  &dof_handler.get_triangulation() == &cache.get_triangulation(),
255  ExcMessage(
256  "The provided Utilities::MPI::RemotePointEvaluation and DoFHandler "
257  "object have been set up with different Triangulation objects, "
258  "a scenario not supported!"));
259 
260  // evaluate values at points if possible
261  const auto evaluation_point_results = [&]() {
262  // helper function for accessing the global vector and interpolating
263  // the results onto the points
264  const auto evaluation_function = [&](auto & values,
265  const auto &cell_data) {
266  std::vector<typename VectorType::value_type> solution_values;
267 
268  std::vector<std::unique_ptr<FEPointEvaluation<n_components, dim>>>
269  evaluators(dof_handler.get_fe_collection().size());
270 
271  for (unsigned int i = 0; i < cell_data.cells.size(); ++i)
272  {
273  typename DoFHandler<dim>::active_cell_iterator cell = {
274  &cache.get_triangulation(),
275  cell_data.cells[i].first,
276  cell_data.cells[i].second,
277  &dof_handler};
278 
279  const ArrayView<const Point<dim>> unit_points(
280  cell_data.reference_point_values.data() +
281  cell_data.reference_point_ptrs[i],
282  cell_data.reference_point_ptrs[i + 1] -
283  cell_data.reference_point_ptrs[i]);
284 
285  solution_values.resize(
286  dof_handler.get_fe(cell->active_fe_index()).n_dofs_per_cell());
287  cell->get_dof_values(vector,
288  solution_values.begin(),
289  solution_values.end());
290 
291  if (evaluators[cell->active_fe_index()] == nullptr)
292  evaluators[cell->active_fe_index()] =
293  std::make_unique<FEPointEvaluation<n_components, dim>>(
294  cache.get_mapping(), cell->get_fe(), update_flags);
295  auto &evaluator = *evaluators[cell->active_fe_index()];
296 
297  evaluator.reinit(cell, unit_points);
298  evaluator.evaluate(solution_values, evaluation_flags);
299 
300  for (unsigned int q = 0; q < unit_points.size(); ++q)
301  values[q + cell_data.reference_point_ptrs[i]] =
302  process_quadrature_point(evaluator, q);
303  }
304  };
305 
306  std::vector<value_type> evaluation_point_results;
307  std::vector<value_type> buffer;
308 
309  cache.template evaluate_and_process<value_type>(
310  evaluation_point_results, buffer, evaluation_function);
311 
312  return evaluation_point_results;
313  }();
314 
315  if (cache.is_map_unique())
316  {
317  // each point has exactly one result (unique map)
318  return evaluation_point_results;
319  }
320  else
321  {
322  // map is not unique (multiple or no results): postprocessing is
323  // needed
324  std::vector<value_type> unique_evaluation_point_results(
325  cache.get_point_ptrs().size() - 1);
326 
327  const auto &ptr = cache.get_point_ptrs();
328 
329  for (unsigned int i = 0; i < ptr.size() - 1; ++i)
330  {
331  const auto n_entries = ptr[i + 1] - ptr[i];
332  if (n_entries == 0)
333  continue;
334 
335  unique_evaluation_point_results[i] =
336  reduce(flags,
338  evaluation_point_results.data() + ptr[i], n_entries));
339  }
340 
341  return unique_evaluation_point_results;
342  }
343  }
344  } // namespace internal
345 
346  template <int n_components, int dim, int spacedim, typename VectorType>
347  inline std::vector<typename FEPointEvaluation<n_components, dim>::value_type>
348  point_values(
350  const DoFHandler<dim, spacedim> & dof_handler,
351  const VectorType & vector,
353  {
354  return internal::evaluate_at_points<
355  n_components,
356  dim,
357  spacedim,
358  VectorType,
360  cache,
361  dof_handler,
362  vector,
363  flags,
366  [](const auto &evaluator, const auto &q) {
367  return evaluator.get_value(q);
368  });
369  }
370 
371  template <int n_components, int dim, int spacedim, typename VectorType>
372  inline std::vector<
376  const DoFHandler<dim, spacedim> & dof_handler,
377  const VectorType & vector,
379  {
380  return internal::evaluate_at_points<
381  n_components,
382  dim,
383  spacedim,
384  VectorType,
386  cache,
387  dof_handler,
388  vector,
389  flags,
392  [](const auto &evaluator, const auto &q) {
393  return evaluator.get_gradient(q);
394  });
395  }
396 
397 #endif
398 } // namespace VectorTools
399 
401 
402 #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:592
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:574
static ::ExceptionBase & ExcMessage(std::string arg1)
static const char T
#define Assert(cond, exc)
Definition: exceptions.h:1465
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:583
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)