Reference documentation for deal.II version 9.3.0
\(\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 
71  template <int n_components, int dim, int spacedim, typename VectorType>
72  std::vector<typename FEPointEvaluation<n_components, dim>::value_type>
74  const Mapping<dim> & mapping,
75  const DoFHandler<dim, spacedim> & dof_handler,
76  const VectorType & vector,
77  const std::vector<Point<spacedim>> & evaluation_points,
80 
89  template <int n_components, int dim, int spacedim, typename VectorType>
90  std::vector<typename FEPointEvaluation<n_components, dim>::value_type>
93  const DoFHandler<dim, spacedim> & dof_handler,
94  const VectorType & vector,
96 
97 
98 
99  // inlined functions
100 
101 
102 #ifndef DOXYGEN
103  template <int n_components, int dim, int spacedim, typename VectorType>
104  inline std::vector<typename FEPointEvaluation<n_components, dim>::value_type>
105  point_values(const Mapping<dim> & mapping,
106  const DoFHandler<dim, spacedim> & dof_handler,
107  const VectorType & vector,
108  const std::vector<Point<spacedim>> &evaluation_points,
111  {
112  cache.reinit(evaluation_points, dof_handler.get_triangulation(), mapping);
113 
114  return point_values<n_components>(cache, dof_handler, vector, flags);
115  }
116 
117 
118 
119  namespace internal
120  {
124  template <typename T>
125  T
126  reduce(const EvaluationFlags::EvaluationFlags &flags,
127  const ArrayView<const T> & values)
128  {
129  switch (flags)
130  {
132  {
133  return std::accumulate(values.begin(), values.end(), T{}) /
134  (T(1.0) * values.size());
135  }
137  return *std::max_element(values.begin(), values.end());
139  return *std::min_element(values.begin(), values.end());
141  return values[0];
142  default:
143  Assert(false, ExcNotImplemented());
144  return values[0];
145  }
146  }
147 
148 
149 
153  template <int rank, int dim, typename Number>
155  reduce(const EvaluationFlags::EvaluationFlags & flags,
156  const ArrayView<const Tensor<rank, dim, Number>> &values)
157  {
158  switch (flags)
159  {
161  {
162  return std::accumulate(values.begin(),
163  values.end(),
165  (Number(1.0) * values.size());
166  }
168  return values[0];
169  default:
170  Assert(false, ExcNotImplemented());
171  return values[0];
172  }
173  }
174  } // namespace internal
175 
176 
177 
178  template <int n_components, int dim, int spacedim, typename VectorType>
179  inline std::vector<typename FEPointEvaluation<n_components, dim>::value_type>
180  point_values(
182  const DoFHandler<dim, spacedim> & dof_handler,
183  const VectorType & vector,
185  {
186  using value_type =
188 
189  Assert(cache.is_ready(),
190  ExcMessage(
191  "Utilities::MPI::RemotePointEvaluation is not ready yet! "
192  "Please call Utilities::MPI::RemotePointEvaluation::reinit() "
193  "yourself or the other point_values(), which does this for"
194  "you."));
195 
196  Assert(
197  &dof_handler.get_triangulation() == &cache.get_triangulation(),
198  ExcMessage(
199  "The provided Utilities::MPI::RemotePointEvaluation and DoFHandler "
200  "object have been set up with different Triangulation objects, "
201  "a scenario not supported!"));
202 
203  // evaluate values at points if possible
204  const auto evaluation_point_results = [&]() {
205  // helper function for accessing the global vector and interpolating
206  // the results onto the points
207  const auto evaluation_function = [&](auto & values,
208  const auto &cell_data) {
209  std::vector<typename VectorType::value_type> solution_values;
210 
211  std::vector<std::unique_ptr<FEPointEvaluation<n_components, dim>>>
212  evaluators(dof_handler.get_fe_collection().size());
213 
214  const auto get_evaluator = [&](const unsigned int active_fe_index)
216  if (evaluators[active_fe_index] == nullptr)
217  evaluators[active_fe_index] =
218  std::make_unique<FEPointEvaluation<n_components, dim>>(
219  cache.get_mapping(),
220  dof_handler.get_fe(active_fe_index),
221  update_values);
222 
223  return *evaluators[active_fe_index];
224  };
225 
226  for (unsigned int i = 0; i < cell_data.cells.size(); ++i)
227  {
228  typename DoFHandler<dim>::active_cell_iterator cell = {
229  &cache.get_triangulation(),
230  cell_data.cells[i].first,
231  cell_data.cells[i].second,
232  &dof_handler};
233 
234  const ArrayView<const Point<dim>> unit_points(
235  cell_data.reference_point_values.data() +
236  cell_data.reference_point_ptrs[i],
237  cell_data.reference_point_ptrs[i + 1] -
238  cell_data.reference_point_ptrs[i]);
239 
240  solution_values.resize(
241  dof_handler.get_fe(cell->active_fe_index()).n_dofs_per_cell());
242  cell->get_dof_values(vector,
243  solution_values.begin(),
244  solution_values.end());
245 
246  auto &evaluator = get_evaluator(cell->active_fe_index());
247 
248  evaluator.reinit(cell, unit_points);
249  evaluator.evaluate(solution_values,
251 
252  for (unsigned int q = 0; q < unit_points.size(); ++q)
253  values[q + cell_data.reference_point_ptrs[i]] =
254  evaluator.get_value(q);
255  }
256  };
257 
258  std::vector<value_type> evaluation_point_results;
259  std::vector<value_type> buffer;
260 
261  cache.template evaluate_and_process<value_type>(evaluation_point_results,
262  buffer,
263  evaluation_function);
264 
265  return evaluation_point_results;
266  }();
267 
268  if (cache.is_map_unique())
269  {
270  // each point has exactly one result (unique map)
271  return evaluation_point_results;
272  }
273  else
274  {
275  // map is not unique (multiple or no results): postprocessing is needed
276  std::vector<value_type> unique_evaluation_point_results(
277  cache.get_point_ptrs().size() - 1);
278 
279  const auto &ptr = cache.get_point_ptrs();
280 
281  for (unsigned int i = 0; i < ptr.size() - 1; ++i)
282  {
283  const auto n_entries = ptr[i + 1] - ptr[i];
284  if (n_entries == 0)
285  continue;
286 
287  unique_evaluation_point_results[i] =
288  internal::reduce(flags,
290  evaluation_point_results.data() + ptr[i],
291  n_entries));
292  }
293 
294  return unique_evaluation_point_results;
295  }
296  }
297 #endif
298 } // namespace VectorTools
299 
301 
302 #endif // dealii_vector_tools_boundary_h
Shape function values.
The namespace for the EvaluationFlags enum.
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
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
const hp::FECollection< dim, spacedim > & get_fe_collection() const
const FiniteElement< dim, spacedim > & get_fe(const unsigned int index=0) const
static ::ExceptionBase & ExcMessage(std::string arg1)
static const char T
#define Assert(cond, exc)
Definition: exceptions.h:1465
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:395
const Mapping< dim, spacedim > & get_mapping() const
typename internal::FEPointEvaluation::EvaluatorTypeTraits< dim, n_components, Number >::value_type value_type
const Triangulation< dim, spacedim > & get_triangulation() const
const Triangulation< dim, spacedim > & get_triangulation() const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:394
static ::ExceptionBase & ExcNotImplemented()
iterator begin() const
Definition: array_view.h:583