Reference documentation for deal.II version Git 191d06ed00 2021-05-11 21:15:49 -0400
\(\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>
106  const Mapping<dim> & mapping,
107  const DoFHandler<dim, spacedim> & dof_handler,
108  const VectorType & vector,
109  const std::vector<Point<spacedim>> & evaluation_points,
112  {
113  cache.reinit(evaluation_points, dof_handler.get_triangulation(), mapping);
114 
115  return evaluate_at_points<n_components>(cache, dof_handler, vector, flags);
116  }
117 
118 
119 
120  namespace internal
121  {
125  template <typename T>
126  T
127  reduce(const EvaluationFlags::EvaluationFlags &flags,
128  const ArrayView<const T> & values)
129  {
130  switch (flags)
131  {
133  {
134  return std::accumulate(values.begin(), values.end(), T{}) /
135  (T(1.0) * values.size());
136  }
138  return *std::max_element(values.begin(), values.end());
140  return *std::min_element(values.begin(), values.end());
142  return values[0];
143  default:
144  Assert(false, ExcNotImplemented());
145  return values[0];
146  }
147  }
148 
149 
150 
154  template <int rank, int dim, typename Number>
156  reduce(const EvaluationFlags::EvaluationFlags & flags,
157  const ArrayView<const Tensor<rank, dim, Number>> &values)
158  {
159  switch (flags)
160  {
162  {
163  return std::accumulate(values.begin(),
164  values.end(),
166  (Number(1.0) * values.size());
167  }
169  return values[0];
170  default:
171  Assert(false, ExcNotImplemented());
172  return values[0];
173  }
174  }
175  } // namespace internal
176 
177 
178 
179  template <int n_components, int dim, int spacedim, typename VectorType>
180  inline std::vector<typename FEPointEvaluation<n_components, dim>::value_type>
183  const DoFHandler<dim, spacedim> & dof_handler,
184  const VectorType & vector,
186  {
187  using value_type =
189 
190  Assert(cache.is_ready(),
191  ExcMessage(
192  "Utilties::MPI::RemotePointEvaluation is not ready yet! "
193  "Please call Utilties::MPI::RemotePointEvaluation::reinit() "
194  "yourself or the other evaluate_at_points(), which does this for"
195  "you."));
196 
197  Assert(&dof_handler.get_triangulation() == &cache.get_triangulation(),
198  ExcMessage(
199  "The provided Utilties::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(), dof_handler.get_fe(active_fe_index));
220 
221  return *evaluators[active_fe_index];
222  };
223 
224  for (unsigned int i = 0; i < cell_data.cells.size(); ++i)
225  {
226  typename DoFHandler<dim>::active_cell_iterator cell = {
227  &cache.get_triangulation(),
228  cell_data.cells[i].first,
229  cell_data.cells[i].second,
230  &dof_handler};
231 
232  const ArrayView<const Point<dim>> unit_points(
233  cell_data.reference_point_values.data() +
234  cell_data.reference_point_ptrs[i],
235  cell_data.reference_point_ptrs[i + 1] -
236  cell_data.reference_point_ptrs[i]);
237 
238  solution_values.resize(
239  dof_handler.get_fe(cell->active_fe_index()).n_dofs_per_cell());
240  cell->get_dof_values(vector,
241  solution_values.begin(),
242  solution_values.end());
243 
244  auto &evaluator = get_evaluator(cell->active_fe_index());
245 
246  evaluator.evaluate(cell,
247  unit_points,
248  solution_values,
250 
251  for (unsigned int q = 0; q < unit_points.size(); ++q)
252  values[q + cell_data.reference_point_ptrs[i]] =
253  evaluator.get_value(q);
254  }
255  };
256 
257  std::vector<value_type> evaluation_point_results;
258  std::vector<value_type> buffer;
259 
260  cache.template evaluate_and_process<value_type>(evaluation_point_results,
261  buffer,
262  evaluation_function);
263 
264  return evaluation_point_results;
265  }();
266 
267  if (cache.is_map_unique())
268  {
269  // each point has exactly one result (unique map)
270  return evaluation_point_results;
271  }
272  else
273  {
274  // map is not unique (multiple or no results): postprocessing is needed
275  std::vector<value_type> unique_evaluation_point_results(
276  cache.get_point_ptrs().size() - 1);
277 
278  const auto &ptr = cache.get_point_ptrs();
279 
280  for (unsigned int i = 0; i < ptr.size() - 1; ++i)
281  {
282  const auto n_entries = ptr[i + 1] - ptr[i];
283  if (n_entries == 0)
284  continue;
285 
286  unique_evaluation_point_results[i] =
287  internal::reduce(flags,
289  evaluation_point_results.data() + ptr[i],
290  n_entries));
291  }
292 
293  return unique_evaluation_point_results;
294  }
295  }
296 #endif
297 } // namespace VectorTools
298 
300 
301 #endif // dealii_vector_tools_boundary_h
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::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:396
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:395
static ::ExceptionBase & ExcNotImplemented()
iterator begin() const
Definition: array_view.h:583
std::vector< typename FEPointEvaluation< n_components, dim >::value_type > evaluate_at_points(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)