Reference documentation for deal.II version 9.4.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 - 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
32namespace VectorTools
33{
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
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(),
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(),
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 {
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>
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 FiniteElement< dim, spacedim > & get_fe(const unsigned int index=0) const
const hp::FECollection< dim, spacedim > & get_fe_collection() const
const Triangulation< dim, spacedim > & get_triangulation() 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:264
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
UpdateFlags
@ update_values
Shape function values.
@ update_gradients
Shape function gradients.
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
Definition: exceptions.h:1473
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)