Loading [MathJax]/extensions/TeX/AMSsymbols.js
 deal.II version GIT relicensing-3104-g8c3bc2e695 2025-04-21 18:40:00+00:00
\(\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\}}\)
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages Concepts
point_value_history.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2009 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
16
24#include <deal.II/lac/vector.h>
26
29
30#include <algorithm>
31
32
34
35
36namespace internal
37{
38 namespace PointValueHistoryImplementation
39 {
40 template <int dim>
42 const Point<dim> &new_requested_location,
43 const std::vector<Point<dim>> &new_locations,
44 const std::vector<types::global_dof_index> &new_sol_indices)
45 {
46 requested_location = new_requested_location;
47 support_point_locations = new_locations;
48 solution_indices = new_sol_indices;
49 }
50 } // namespace PointValueHistoryImplementation
51} // namespace internal
52
53
54
55template <int dim>
57 const unsigned int n_independent_variables)
58 : n_indep(n_independent_variables)
59{
60 closed = false;
61 cleared = false;
63 have_dof_handler = false;
64
65 // make a vector for keys
66 dataset_key = std::vector<double>(); // initialize the std::vector
67
68 // make a vector of independent values
70 std::vector<std::vector<double>>(n_indep, std::vector<double>(0));
71 indep_names = std::vector<std::string>();
72}
73
74
75
76template <int dim>
78 const DoFHandler<dim> &dof_handler,
79 const unsigned int n_independent_variables)
80 : dof_handler(&dof_handler)
81 , n_indep(n_independent_variables)
82{
83 closed = false;
84 cleared = false;
86 have_dof_handler = true;
87
88 // make a vector to store keys
89 dataset_key = std::vector<double>(); // initialize the std::vector
90
91 // make a vector for the independent values
93 std::vector<std::vector<double>>(n_indep, std::vector<double>(0));
94 indep_names = std::vector<std::string>();
95
96 tria_listener = dof_handler.get_triangulation().signals.any_change.connect(
97 [this]() { this->tria_change_listener(); });
98}
99
100
101
102template <int dim>
104 const PointValueHistory &point_value_history)
105{
106 dataset_key = point_value_history.dataset_key;
107 independent_values = point_value_history.independent_values;
108 indep_names = point_value_history.indep_names;
109 data_store = point_value_history.data_store;
110 component_mask = point_value_history.component_mask;
111 component_names_map = point_value_history.component_names_map;
112 point_geometry_data = point_value_history.point_geometry_data;
113
114 closed = point_value_history.closed;
115 cleared = point_value_history.cleared;
116
117 dof_handler = point_value_history.dof_handler;
118
119 triangulation_changed = point_value_history.triangulation_changed;
120 have_dof_handler = point_value_history.have_dof_handler;
121 n_indep = point_value_history.n_indep;
122
123 // What to do with tria_listener?
124 // Presume subscribe new instance?
125 if (have_dof_handler)
126 {
127 tria_listener =
128 dof_handler->get_triangulation().signals.any_change.connect(
129 [this]() { this->tria_change_listener(); });
130 }
131}
132
133
134
135template <int dim>
138{
139 dataset_key = point_value_history.dataset_key;
140 independent_values = point_value_history.independent_values;
141 indep_names = point_value_history.indep_names;
142 data_store = point_value_history.data_store;
143 component_mask = point_value_history.component_mask;
144 component_names_map = point_value_history.component_names_map;
145 point_geometry_data = point_value_history.point_geometry_data;
146
147 closed = point_value_history.closed;
148 cleared = point_value_history.cleared;
149
150 dof_handler = point_value_history.dof_handler;
151
152 triangulation_changed = point_value_history.triangulation_changed;
153 have_dof_handler = point_value_history.have_dof_handler;
154 n_indep = point_value_history.n_indep;
155
156 // What to do with tria_listener?
157 // Presume subscribe new instance?
158 if (have_dof_handler)
159 {
160 tria_listener =
161 dof_handler->get_triangulation().signals.any_change.connect(
162 [this]() { this->tria_change_listener(); });
163 }
164
165 return *this;
166}
167
168
169
170template <int dim>
172{
173 if (have_dof_handler)
174 {
175 tria_listener.disconnect();
176 }
177}
178
179
180
181template <int dim>
182void
184{
185 // can't be closed to add additional points
186 // or vectors
187 AssertThrow(!closed, ExcInvalidState());
188 AssertThrow(!cleared, ExcInvalidState());
189 AssertThrow(have_dof_handler, ExcDoFHandlerRequired());
190 AssertThrow(!triangulation_changed, ExcDoFHandlerChanged());
191
192 // Implementation assumes that support
193 // points locations are dofs locations
194 AssertThrow(dof_handler->get_fe().has_support_points(), ExcNotImplemented());
195
196 // While in general quadrature points seems
197 // to refer to Gauss quadrature points, in
198 // this case the quadrature points are
199 // forced to be the support points of the
200 // FE.
201 Quadrature<dim> support_point_quadrature(
202 dof_handler->get_fe().get_unit_support_points());
203 FEValues<dim> fe_values(dof_handler->get_fe(),
204 support_point_quadrature,
206 unsigned int n_support_points =
207 dof_handler->get_fe().get_unit_support_points().size();
208 unsigned int n_components = dof_handler->get_fe(0).n_components();
209
210 // set up a loop over all the cells in the
211 // DoFHandler
213 dof_handler->begin_active();
214 typename DoFHandler<dim>::active_cell_iterator endc = dof_handler->end();
215
216 // default values to be replaced as closer
217 // points are found however they need to be
218 // consistent in case they are actually
219 // chosen
220 typename DoFHandler<dim>::active_cell_iterator current_cell = cell;
221 std::vector<unsigned int> current_fe_index(n_components,
222 0); // need one index per component
223 fe_values.reinit(cell);
224 std::vector<Point<dim>> current_points(n_components, Point<dim>());
225 for (unsigned int support_point = 0; support_point < n_support_points;
226 support_point++)
227 {
228 // set up valid data in the empty vectors
229 unsigned int component =
230 dof_handler->get_fe().system_to_component_index(support_point).first;
231 current_points[component] = fe_values.quadrature_point(support_point);
232 current_fe_index[component] = support_point;
233 }
234
235 // check each cell to find a suitable
236 // support points
237 // GridTools::find_active_cell_around_point
238 // is an alternative. That method is not
239 // used here mostly because of the history
240 // of the class. The algorithm used in
241 // add_points below may be slightly more
242 // efficient than find_active_cell_around_point
243 // because it operates on a set of points.
244
245 for (; cell != endc; ++cell)
246 {
247 fe_values.reinit(cell);
248
249 for (unsigned int support_point = 0; support_point < n_support_points;
250 support_point++)
251 {
252 unsigned int component = dof_handler->get_fe()
253 .system_to_component_index(support_point)
254 .first;
255 const Point<dim> &test_point =
256 fe_values.quadrature_point(support_point);
257
258 if (location.distance(test_point) <
259 location.distance(current_points[component]))
260 {
261 // save the data
262 current_points[component] = test_point;
263 current_cell = cell;
264 current_fe_index[component] = support_point;
265 }
266 }
267 }
268
269
270 std::vector<types::global_dof_index> local_dof_indices(
271 dof_handler->get_fe().n_dofs_per_cell());
272 std::vector<types::global_dof_index> new_solution_indices;
273 current_cell->get_dof_indices(local_dof_indices);
274 // there is an implicit assumption here
275 // that all the closest support point to
276 // the requested point for all finite
277 // element components lie in the same cell.
278 // this could possibly be violated if
279 // components use different FE orders,
280 // requested points are on the edge or
281 // vertex of a cell and we are unlucky with
282 // floating point rounding. Worst case
283 // scenario however is that the point
284 // selected isn't the closest possible, it
285 // will still lie within one cell distance.
286 // calling
287 // GridTools::find_active_cell_around_point
288 // to obtain a cell to search is an
289 // option for these methods, but currently
290 // the GridTools function does not cater for
291 // a vector of points, and does not seem to
292 // be intrinsicly faster than this method.
293 new_solution_indices.reserve(dof_handler->get_fe(0).n_components());
294 for (unsigned int component = 0;
295 component < dof_handler->get_fe(0).n_components();
296 component++)
297 {
298 new_solution_indices.push_back(
299 local_dof_indices[current_fe_index[component]]);
300 }
301
303 new_point_geometry_data(location, current_points, new_solution_indices);
304 point_geometry_data.push_back(new_point_geometry_data);
305
306 for (auto &data_entry : data_store)
307 {
308 // add an extra row to each vector entry
309 const ComponentMask &current_mask =
310 (component_mask.find(data_entry.first))->second;
311 unsigned int n_stored = current_mask.n_selected_components();
312 data_entry.second.resize(data_entry.second.size() + n_stored);
313 }
314}
315
316
317
318template <int dim>
319void
320PointValueHistory<dim>::add_points(const std::vector<Point<dim>> &locations)
321{
322 // This algorithm adds points in the same
323 // order as they appear in the vector
324 // locations and users may depend on this
325 // so do not change order added!
326
327 // can't be closed to add additional points or vectors
328 AssertThrow(!closed, ExcInvalidState());
329 AssertThrow(!cleared, ExcInvalidState());
330 AssertThrow(have_dof_handler, ExcDoFHandlerRequired());
331 AssertThrow(!triangulation_changed, ExcDoFHandlerChanged());
332
333
334 // Implementation assumes that support
335 // points locations are dofs locations
336 AssertThrow(dof_handler->get_fe().has_support_points(), ExcNotImplemented());
337
338 // While in general quadrature points seems
339 // to refer to Gauss quadrature points, in
340 // this case the quadrature points are
341 // forced to be the support points of the
342 // FE.
343 Quadrature<dim> support_point_quadrature(
344 dof_handler->get_fe().get_unit_support_points());
345 FEValues<dim> fe_values(dof_handler->get_fe(),
346 support_point_quadrature,
348 unsigned int n_support_points =
349 dof_handler->get_fe().get_unit_support_points().size();
350 unsigned int n_components = dof_handler->get_fe(0).n_components();
351
352 // set up a loop over all the cells in the
353 // DoFHandler
355 dof_handler->begin_active();
356 typename DoFHandler<dim>::active_cell_iterator endc = dof_handler->end();
357
358 // default values to be replaced as closer
359 // points are found however they need to be
360 // consistent in case they are actually
361 // chosen vector <vector>s defined where
362 // previously single vectors were used
363
364 // need to store one value per point per component
365 std::vector<typename DoFHandler<dim>::active_cell_iterator> current_cell(
366 locations.size(), cell);
367
368 fe_values.reinit(cell);
369 std::vector<Point<dim>> temp_points(n_components, Point<dim>());
370 std::vector<unsigned int> temp_fe_index(n_components, 0);
371 for (unsigned int support_point = 0; support_point < n_support_points;
372 support_point++)
373 {
374 // set up valid data in the empty vectors
375 unsigned int component =
376 dof_handler->get_fe().system_to_component_index(support_point).first;
377 temp_points[component] = fe_values.quadrature_point(support_point);
378 temp_fe_index[component] = support_point;
379 }
380 std::vector<std::vector<Point<dim>>> current_points(
381 locations.size(), temp_points); // give a valid start point
382 std::vector<std::vector<unsigned int>> current_fe_index(locations.size(),
383 temp_fe_index);
384
385 // check each cell to find suitable support
386 // points
387 // GridTools::find_active_cell_around_point
388 // is an alternative. That method is not
389 // used here mostly because of the history
390 // of the class. The algorithm used here
391 // may be slightly more
392 // efficient than find_active_cell_around_point
393 // because it operates on a set of points.
394 for (; cell != endc; ++cell)
395 {
396 fe_values.reinit(cell);
397 for (unsigned int support_point = 0; support_point < n_support_points;
398 support_point++)
399 {
400 unsigned int component = dof_handler->get_fe()
401 .system_to_component_index(support_point)
402 .first;
403 const Point<dim> &test_point =
404 fe_values.quadrature_point(support_point);
405
406 for (unsigned int point = 0; point < locations.size(); ++point)
407 {
408 if (locations[point].distance(test_point) <
409 locations[point].distance(current_points[point][component]))
410 {
411 // save the data
412 current_points[point][component] = test_point;
413 current_cell[point] = cell;
414 current_fe_index[point][component] = support_point;
415 }
416 }
417 }
418 }
419
420 std::vector<types::global_dof_index> local_dof_indices(
421 dof_handler->get_fe().n_dofs_per_cell());
422 for (unsigned int point = 0; point < locations.size(); ++point)
423 {
424 current_cell[point]->get_dof_indices(local_dof_indices);
425 std::vector<types::global_dof_index> new_solution_indices;
426
427 new_solution_indices.reserve(dof_handler->get_fe(0).n_components());
428 for (unsigned int component = 0;
429 component < dof_handler->get_fe(0).n_components();
430 component++)
431 {
432 new_solution_indices.push_back(
433 local_dof_indices[current_fe_index[point][component]]);
434 }
435
437 new_point_geometry_data(locations[point],
438 current_points[point],
439 new_solution_indices);
440
441 point_geometry_data.push_back(new_point_geometry_data);
442
443 for (auto &data_entry : data_store)
444 {
445 // add an extra row to each vector entry
446 const ComponentMask current_mask =
447 (component_mask.find(data_entry.first))->second;
448 unsigned int n_stored = current_mask.n_selected_components();
449 data_entry.second.resize(data_entry.second.size() + n_stored);
450 }
451 }
452}
453
454
455
456template <int dim>
457void
458PointValueHistory<dim>::add_field_name(const std::string &vector_name,
459 const ComponentMask &mask)
460{
461 // can't be closed to add additional points
462 // or vectors
463 AssertThrow(!closed, ExcInvalidState());
464 AssertThrow(!cleared, ExcInvalidState());
465 AssertThrow(have_dof_handler, ExcDoFHandlerRequired());
466 AssertThrow(!triangulation_changed, ExcDoFHandlerChanged());
467
468 // insert a component mask that is always of the right size
469 if (mask.represents_the_all_selected_mask() == false)
470 component_mask.insert(std::make_pair(vector_name, mask));
471 else
472 component_mask.insert(
473 std::make_pair(vector_name,
474 ComponentMask(std::vector<bool>(
475 dof_handler->get_fe(0).n_components(), true))));
476
477 // insert an empty vector of strings
478 // to ensure each field has an entry
479 // in the map
480 std::pair<std::string, std::vector<std::string>> empty_names(
481 vector_name, std::vector<std::string>());
482 component_names_map.insert(empty_names);
483
484 // make and add a new vector
485 // point_geometry_data.size() long
486 std::pair<std::string, std::vector<std::vector<double>>> pair_data;
487 pair_data.first = vector_name;
488 const unsigned int n_stored =
489 (mask.represents_the_all_selected_mask() == false ?
490 mask.n_selected_components() :
491 dof_handler->get_fe(0).n_components());
492
493 int n_datastreams =
494 point_geometry_data.size() * n_stored; // each point has n_stored sub parts
495 std::vector<std::vector<double>> vector_size(n_datastreams,
496 std::vector<double>(0));
497 pair_data.second = vector_size;
498 data_store.insert(pair_data);
499}
500
501
502template <int dim>
503void
504PointValueHistory<dim>::add_field_name(const std::string &vector_name,
505 const unsigned int n_components)
506{
507 ComponentMask temp_mask(std::vector<bool>(n_components, true));
508 add_field_name(vector_name, temp_mask);
509}
510
511
512template <int dim>
513void
515 const std::string &vector_name,
516 const std::vector<std::string> &component_names)
517{
518 typename std::map<std::string, std::vector<std::string>>::iterator names =
519 component_names_map.find(vector_name);
520 Assert(names != component_names_map.end(),
521 ExcMessage("vector_name not in class"));
522
523 typename std::map<std::string, ComponentMask>::iterator mask =
524 component_mask.find(vector_name);
525 Assert(mask != component_mask.end(), ExcMessage("vector_name not in class"));
526 unsigned int n_stored = mask->second.n_selected_components();
527 (void)n_stored;
528 Assert(component_names.size() == n_stored,
529 ExcDimensionMismatch(component_names.size(), n_stored));
530
531 names->second = component_names;
532}
533
534
535template <int dim>
536void
538 const std::vector<std::string> &independent_names)
539{
540 Assert(independent_names.size() == n_indep,
541 ExcDimensionMismatch(independent_names.size(), n_indep));
542
543 indep_names = independent_names;
544}
545
546
547template <int dim>
548void
550{
551 closed = true;
552}
553
554
555
556template <int dim>
557void
559{
560 cleared = true;
561 dof_handler = nullptr;
562 have_dof_handler = false;
563}
564
565// Need to test that the internal data has a full and complete dataset for
566// each key. That is that the data has not got 'out of sync'. Testing that
567// dataset_key is within 1 of independent_values is cheap and is done in all
568// three methods. Evaluate_field will check that its vector_name is within 1
569// of dataset_key. However this leaves the possibility that the user has
570// neglected to call evaluate_field on one vector_name consistently. To catch
571// this case start_new_dataset will call bool deap_check () which will test
572// all vector_names and return a bool. This can be called from an Assert
573// statement.
574
575
576
577template <int dim>
578template <typename VectorType>
579void
580PointValueHistory<dim>::evaluate_field(const std::string &vector_name,
581 const VectorType &solution)
582{
583 // must be closed to add data to internal
584 // members.
585 Assert(closed, ExcInvalidState());
586 Assert(!cleared, ExcInvalidState());
587 AssertThrow(have_dof_handler, ExcDoFHandlerRequired());
588 AssertThrow(!triangulation_changed, ExcDoFHandlerChanged());
589
590 if (n_indep != 0) // hopefully this will get optimized, can't test
591 // independent_values[0] unless n_indep > 0
592 {
593 Assert(std::abs(static_cast<int>(dataset_key.size()) -
594 static_cast<int>(independent_values[0].size())) < 2,
595 ExcDataLostSync());
596 }
597 // Look up the field name and get an
598 // iterator for the map. Doing this
599 // up front means that it only needs
600 // to be done once and also allows us
601 // to check vector_name is in the map.
602 typename std::map<std::string, std::vector<std::vector<double>>>::iterator
603 data_store_field = data_store.find(vector_name);
604 Assert(data_store_field != data_store.end(),
605 ExcMessage("vector_name not in class"));
606 // Repeat for component_mask
607 typename std::map<std::string, ComponentMask>::iterator mask =
608 component_mask.find(vector_name);
609 Assert(mask != component_mask.end(), ExcMessage("vector_name not in class"));
610
611 unsigned int n_stored =
612 mask->second.n_selected_components(dof_handler->get_fe(0).n_components());
613
614 typename std::vector<
616 point = point_geometry_data.begin();
617 for (unsigned int data_store_index = 0; point != point_geometry_data.end();
618 ++point, ++data_store_index)
619 {
620 // Look up the components to add
621 // in the component_mask, and
622 // access the data associated with
623 // those components
624
625 for (unsigned int store_index = 0, comp = 0;
626 comp < dof_handler->get_fe(0).n_components();
627 comp++)
628 {
629 if (mask->second[comp])
630 {
631 unsigned int solution_index = point->solution_indices[comp];
632 data_store_field
633 ->second[data_store_index * n_stored + store_index]
634 .push_back(
636 solution_index));
637 ++store_index;
638 }
639 }
640 }
641}
642
643
644
645template <int dim>
646template <typename VectorType>
647void
649 const std::vector<std::string> &vector_names,
650 const VectorType &solution,
651 const DataPostprocessor<dim> &data_postprocessor,
652 const Quadrature<dim> &quadrature)
653{
654 // must be closed to add data to internal
655 // members.
656 Assert(closed, ExcInvalidState());
657 Assert(!cleared, ExcInvalidState());
658 AssertThrow(have_dof_handler, ExcDoFHandlerRequired());
659 if (n_indep != 0) // hopefully this will get optimized, can't test
660 // independent_values[0] unless n_indep > 0
661 {
662 Assert(std::abs(static_cast<int>(dataset_key.size()) -
663 static_cast<int>(independent_values[0].size())) < 2,
664 ExcDataLostSync());
665 }
666
667 // Make an FEValues object
668 const UpdateFlags update_flags =
670 Assert(
671 !(update_flags & update_normal_vectors),
673 "The update of normal vectors may not be requested for evaluation of "
674 "data on cells via DataPostprocessor."));
675 FEValues<dim> fe_values(dof_handler->get_fe(), quadrature, update_flags);
676 unsigned int n_components = dof_handler->get_fe(0).n_components();
677 unsigned int n_quadrature_points = quadrature.size();
678
679 unsigned int n_output_variables = data_postprocessor.get_names().size();
680
681 // declare some temp objects for evaluating the solution at quadrature
682 // points. we will either need the scalar or vector version in the code
683 // below
684 std::vector<typename VectorType::value_type> scalar_solution_values(
685 n_quadrature_points);
686 std::vector<Tensor<1, dim, typename VectorType::value_type>>
687 scalar_solution_gradients(n_quadrature_points);
688 std::vector<Tensor<2, dim, typename VectorType::value_type>>
689 scalar_solution_hessians(n_quadrature_points);
690
691 std::vector<Vector<typename VectorType::value_type>> vector_solution_values(
692 n_quadrature_points, Vector<typename VectorType::value_type>(n_components));
693
694 std::vector<std::vector<Tensor<1, dim, typename VectorType::value_type>>>
695 vector_solution_gradients(
696 n_quadrature_points,
699
700 std::vector<std::vector<Tensor<2, dim, typename VectorType::value_type>>>
701 vector_solution_hessians(
702 n_quadrature_points,
705
706 // Loop over points and find correct cell
707 typename std::vector<
709 point = point_geometry_data.begin();
710 Assert(!dof_handler->get_triangulation().is_mixed_mesh(),
712 const auto reference_cell =
713 dof_handler->get_triangulation().get_reference_cells()[0];
714 for (unsigned int data_store_index = 0; point != point_geometry_data.end();
715 ++point, ++data_store_index)
716 {
717 // we now have a point to query, need to know what cell it is in
718 const Point<dim> requested_location = point->requested_location;
719 const typename DoFHandler<dim>::active_cell_iterator cell =
721 reference_cell.template get_default_linear_mapping<dim, dim>(),
722 *dof_handler,
723 requested_location)
724 .first;
725
726
727 fe_values.reinit(cell);
728 std::vector<Vector<double>> computed_quantities(
729 1, Vector<double>(n_output_variables)); // just one point needed
730
731 // find the closest quadrature point
732 std::vector<Point<dim>> quadrature_points =
733 fe_values.get_quadrature_points();
734 double distance = cell->diameter();
735 unsigned int selected_point = 0;
736 for (unsigned int q_point = 0; q_point < n_quadrature_points; ++q_point)
737 {
738 if (requested_location.distance(quadrature_points[q_point]) <
739 distance)
740 {
741 selected_point = q_point;
742 distance =
743 requested_location.distance(quadrature_points[q_point]);
744 }
745 }
746
747
748 // The case of a scalar FE
749 if (n_components == 1)
750 {
751 // Extract data for the DataPostprocessor object
752 DataPostprocessorInputs::Scalar<dim> postprocessor_input;
753
754 // for each quantity that is requested (values, gradients, hessians),
755 // first get them at all quadrature points, then restrict to the
756 // one value on the quadrature point closest to the evaluation
757 // point in question
758 //
759 // we need temporary objects because the underlying scalar
760 // type of the solution vector may be different from 'double',
761 // but the DataPostprocessorInputs only allow for 'double'
762 // data types
763 if (update_flags & update_values)
764 {
765 fe_values.get_function_values(solution, scalar_solution_values);
766 postprocessor_input.solution_values =
767 std::vector<double>(1, scalar_solution_values[selected_point]);
768 }
769 if (update_flags & update_gradients)
770 {
771 fe_values.get_function_gradients(solution,
772 scalar_solution_gradients);
773 postprocessor_input.solution_gradients =
774 std::vector<Tensor<1, dim>>(
775 1, scalar_solution_gradients[selected_point]);
776 }
777 if (update_flags & update_hessians)
778 {
779 fe_values.get_function_hessians(solution,
780 scalar_solution_hessians);
781 postprocessor_input.solution_hessians =
782 std::vector<Tensor<2, dim>>(
783 1, scalar_solution_hessians[selected_point]);
784 }
785
786 // then also set the single evaluation point
787 postprocessor_input.evaluation_points =
788 std::vector<Point<dim>>(1, quadrature_points[selected_point]);
789
790 // and finally do the postprocessing
791 data_postprocessor.evaluate_scalar_field(postprocessor_input,
792 computed_quantities);
793 }
794 else // The case of a vector FE
795 {
796 // exact same idea as above
797 DataPostprocessorInputs::Vector<dim> postprocessor_input;
798
799 if (update_flags & update_values)
800 {
801 fe_values.get_function_values(solution, vector_solution_values);
802 postprocessor_input.solution_values.resize(
803 1, Vector<double>(n_components));
804 std::copy(vector_solution_values[selected_point].begin(),
805 vector_solution_values[selected_point].end(),
806 postprocessor_input.solution_values[0].begin());
807 }
808 if (update_flags & update_gradients)
809 {
810 fe_values.get_function_gradients(solution,
811 vector_solution_gradients);
812 postprocessor_input.solution_gradients.resize(
813 1, std::vector<Tensor<1, dim>>(n_components));
814 std::copy(vector_solution_gradients[selected_point].begin(),
815 vector_solution_gradients[selected_point].end(),
816 postprocessor_input.solution_gradients[0].begin());
817 }
818 if (update_flags & update_hessians)
819 {
820 fe_values.get_function_hessians(solution,
821 vector_solution_hessians);
822 postprocessor_input.solution_hessians.resize(
823 1, std::vector<Tensor<2, dim>>(n_components));
824 std::copy(vector_solution_hessians[selected_point].begin(),
825 vector_solution_hessians[selected_point].end(),
826 postprocessor_input.solution_hessians[0].begin());
827 }
828
829 postprocessor_input.evaluation_points =
830 std::vector<Point<dim>>(1, quadrature_points[selected_point]);
831
832 data_postprocessor.evaluate_vector_field(postprocessor_input,
833 computed_quantities);
834 }
835
836
837 // we now have the data and need to save it
838 // loop over data names
839 typename std::vector<std::string>::const_iterator name =
840 vector_names.begin();
841 for (; name != vector_names.end(); ++name)
842 {
843 typename std::map<std::string,
844 std::vector<std::vector<double>>>::iterator
845 data_store_field = data_store.find(*name);
846 Assert(data_store_field != data_store.end(),
847 ExcMessage("vector_name not in class"));
848 // Repeat for component_mask
849 typename std::map<std::string, ComponentMask>::iterator mask =
850 component_mask.find(*name);
851 Assert(mask != component_mask.end(),
852 ExcMessage("vector_name not in class"));
853
854 unsigned int n_stored =
855 mask->second.n_selected_components(n_output_variables);
856
857 // Push back computed quantities according
858 // to the component_mask.
859 for (unsigned int store_index = 0, comp = 0;
860 comp < n_output_variables;
861 comp++)
862 {
863 if (mask->second[comp])
864 {
865 data_store_field
866 ->second[data_store_index * n_stored + store_index]
867 .push_back(computed_quantities[0](comp));
868 ++store_index;
869 }
870 }
871 }
872 } // end of loop over points
873}
874
875
876
877template <int dim>
878template <typename VectorType>
879void
881 const std::string &vector_name,
882 const VectorType &solution,
883 const DataPostprocessor<dim> &data_postprocessor,
884 const Quadrature<dim> &quadrature)
885{
886 std::vector<std::string> vector_names;
887 vector_names.push_back(vector_name);
888 evaluate_field(vector_names, solution, data_postprocessor, quadrature);
889}
890
891
892
893template <int dim>
894template <typename VectorType>
895void
897 const std::string &vector_name,
898 const VectorType &solution)
899{
900 using number = typename VectorType::value_type;
901 // must be closed to add data to internal
902 // members.
903 Assert(closed, ExcInvalidState());
904 Assert(!cleared, ExcInvalidState());
905 AssertThrow(have_dof_handler, ExcDoFHandlerRequired());
906
907 if (n_indep != 0) // hopefully this will get optimized, can't test
908 // independent_values[0] unless n_indep > 0
909 {
910 Assert(std::abs(static_cast<int>(dataset_key.size()) -
911 static_cast<int>(independent_values[0].size())) < 2,
912 ExcDataLostSync());
913 }
914 // Look up the field name and get an
915 // iterator for the map. Doing this
916 // up front means that it only needs
917 // to be done once and also allows us
918 // to check vector_name is in the map.
919 typename std::map<std::string, std::vector<std::vector<double>>>::iterator
920 data_store_field = data_store.find(vector_name);
921 Assert(data_store_field != data_store.end(),
922 ExcMessage("vector_name not in class"));
923 // Repeat for component_mask
924 typename std::map<std::string, ComponentMask>::iterator mask =
925 component_mask.find(vector_name);
926 Assert(mask != component_mask.end(), ExcMessage("vector_name not in class"));
927
928 unsigned int n_stored =
929 mask->second.n_selected_components(dof_handler->get_fe(0).n_components());
930
931 typename std::vector<
933 point = point_geometry_data.begin();
934 Vector<number> value(dof_handler->get_fe(0).n_components());
935 for (unsigned int data_store_index = 0; point != point_geometry_data.end();
936 ++point, ++data_store_index)
937 {
938 // Make a Vector <double> for the value
939 // at the point. It will have as many
940 // components as there are in the fe.
941 VectorTools::point_value(*dof_handler,
942 solution,
943 point->requested_location,
944 value);
945
946 // Look up the component_mask and add
947 // in components according to that mask
948 for (unsigned int store_index = 0, comp = 0; comp < mask->second.size();
949 comp++)
950 {
951 if (mask->second[comp])
952 {
953 data_store_field
954 ->second[data_store_index * n_stored + store_index]
955 .push_back(value(comp));
956 ++store_index;
957 }
958 }
959 }
960}
961
962
963template <int dim>
964void
966{
967 // must be closed to add data to internal
968 // members.
969 Assert(closed, ExcInvalidState());
970 Assert(!cleared, ExcInvalidState());
971 Assert(deep_check(false), ExcDataLostSync());
972
973 dataset_key.push_back(key);
974}
975
976
977
978template <int dim>
979void
981 const std::vector<double> &indep_values)
982{
983 // must be closed to add data to internal
984 // members.
985 Assert(closed, ExcInvalidState());
986 Assert(!cleared, ExcInvalidState());
987 Assert(indep_values.size() == n_indep,
988 ExcDimensionMismatch(indep_values.size(), n_indep));
989 Assert(n_indep != 0, ExcNoIndependent());
990 Assert(std::abs(static_cast<int>(dataset_key.size()) -
991 static_cast<int>(independent_values[0].size())) < 2,
992 ExcDataLostSync());
993
994 for (unsigned int component = 0; component < n_indep; ++component)
995 independent_values[component].push_back(indep_values[component]);
996}
997
998
999
1000template <int dim>
1001void
1003 const std::string &base_name,
1004 const std::vector<Point<dim>> &postprocessor_locations)
1005{
1006 AssertThrow(closed, ExcInvalidState());
1007 AssertThrow(!cleared, ExcInvalidState());
1008 AssertThrow(deep_check(true), ExcDataLostSync());
1009
1010 // write inputs to a file
1011 if (n_indep != 0)
1012 {
1013 std::string filename = base_name + "_indep.gpl";
1014 std::ofstream to_gnuplot(filename);
1015
1016 to_gnuplot << "# Data independent of mesh location\n";
1017
1018 // write column headings
1019 to_gnuplot << "# <Key> ";
1020
1021 if (indep_names.size() > 0)
1022 {
1023 for (const auto &indep_name : indep_names)
1024 {
1025 to_gnuplot << "<" << indep_name << "> ";
1026 }
1027 to_gnuplot << '\n';
1028 }
1029 else
1030 {
1031 for (unsigned int component = 0; component < n_indep; ++component)
1032 {
1033 to_gnuplot << "<Indep_" << component << "> ";
1034 }
1035 to_gnuplot << '\n';
1036 }
1037 // write general data stored
1038 for (unsigned int key = 0; key < dataset_key.size(); ++key)
1039 {
1040 to_gnuplot << dataset_key[key];
1041
1042 for (unsigned int component = 0; component < n_indep; ++component)
1043 {
1044 to_gnuplot << " " << independent_values[component][key];
1045 }
1046 to_gnuplot << '\n';
1047 }
1048
1049 to_gnuplot.close();
1050 }
1051
1052
1053
1054 // write points to a file
1055 if (have_dof_handler)
1056 {
1057 AssertThrow(have_dof_handler, ExcDoFHandlerRequired());
1058 AssertThrow(postprocessor_locations.empty() ||
1059 postprocessor_locations.size() ==
1060 point_geometry_data.size(),
1061 ExcDimensionMismatch(postprocessor_locations.size(),
1062 point_geometry_data.size()));
1063 // We previously required the
1064 // number of dofs to remain the
1065 // same to provide some sort of
1066 // test on the relevance of the
1067 // support point indices stored.
1068 // We now relax that to allow
1069 // adaptive refinement strategies
1070 // to make use of the
1071 // evaluate_field_requested_locations
1072 // method. Note that the support point
1073 // information is not meaningful if
1074 // the number of dofs has changed.
1075 // AssertThrow (!triangulation_changed, ExcDoFHandlerChanged ());
1076
1077 typename std::vector<internal::PointValueHistoryImplementation::
1078 PointGeometryData<dim>>::iterator point =
1079 point_geometry_data.begin();
1080 for (unsigned int data_store_index = 0;
1081 point != point_geometry_data.end();
1082 ++point, ++data_store_index)
1083 {
1084 // for each point, open a file to
1085 // be written to
1086 std::string filename = base_name + "_" +
1087 Utilities::int_to_string(data_store_index, 2) +
1088 ".gpl"; // store by order pushed back
1089 // due to
1090 // Utilities::int_to_string(data_store_index,
1091 // 2) call, can handle up to 100
1092 // points
1093 std::ofstream to_gnuplot(filename);
1094
1095 // put helpful info about the
1096 // support point into the file as
1097 // comments
1098 to_gnuplot << "# Requested location: " << point->requested_location
1099 << '\n';
1100 to_gnuplot << "# DoF_index : Support location (for each component)\n";
1101 for (unsigned int component = 0;
1102 component < dof_handler->get_fe(0).n_components();
1103 component++)
1104 {
1105 to_gnuplot << "# " << point->solution_indices[component] << " : "
1106 << point->support_point_locations[component] << '\n';
1107 }
1108 if (triangulation_changed)
1109 to_gnuplot
1110 << "# (Original components and locations, may be invalidated by mesh change.)\n";
1111
1112 if (postprocessor_locations.size() != 0)
1113 {
1114 to_gnuplot << "# Postprocessor location: "
1115 << postprocessor_locations[data_store_index];
1116 if (triangulation_changed)
1117 to_gnuplot << " (may be approximate)\n";
1118 }
1119 to_gnuplot << "#\n";
1120
1121
1122 // write column headings
1123 to_gnuplot << "# <Key> ";
1124
1125 if (indep_names.size() > 0)
1126 {
1127 for (const auto &indep_name : indep_names)
1128 {
1129 to_gnuplot << "<" << indep_name << "> ";
1130 }
1131 }
1132 else
1133 {
1134 for (unsigned int component = 0; component < n_indep; ++component)
1135 {
1136 to_gnuplot << "<Indep_" << component << "> ";
1137 }
1138 }
1139
1140 for (const auto &data_entry : data_store)
1141 {
1142 typename std::map<std::string, ComponentMask>::iterator mask =
1143 component_mask.find(data_entry.first);
1144 unsigned int n_stored = mask->second.n_selected_components();
1145 std::vector<std::string> names =
1146 (component_names_map.find(data_entry.first))->second;
1147
1148 if (names.size() > 0)
1149 {
1150 AssertThrow(names.size() == n_stored,
1151 ExcDimensionMismatch(names.size(), n_stored));
1152 for (const auto &name : names)
1153 {
1154 to_gnuplot << "<" << name << "> ";
1155 }
1156 }
1157 else
1158 {
1159 for (unsigned int component = 0; component < n_stored;
1160 component++)
1161 {
1162 to_gnuplot << "<" << data_entry.first << "_" << component
1163 << "> ";
1164 }
1165 }
1166 }
1167 to_gnuplot << '\n';
1168
1169 // write data stored for the point
1170 for (unsigned int key = 0; key < dataset_key.size(); ++key)
1171 {
1172 to_gnuplot << dataset_key[key];
1173
1174 for (unsigned int component = 0; component < n_indep; ++component)
1175 {
1176 to_gnuplot << " " << independent_values[component][key];
1177 }
1178
1179 for (const auto &data_entry : data_store)
1180 {
1181 typename std::map<std::string, ComponentMask>::iterator mask =
1182 component_mask.find(data_entry.first);
1183 unsigned int n_stored = mask->second.n_selected_components();
1184
1185 for (unsigned int component = 0; component < n_stored;
1186 component++)
1187 {
1188 to_gnuplot
1189 << " "
1190 << (data_entry.second)[data_store_index * n_stored +
1191 component][key];
1192 }
1193 }
1194 to_gnuplot << '\n';
1195 }
1196
1197 to_gnuplot.close();
1198 }
1199 }
1200}
1201
1202
1203
1204template <int dim>
1207{
1208 // a method to put a one at each point on
1209 // the grid where a location is defined
1210 AssertThrow(!cleared, ExcInvalidState());
1211 AssertThrow(have_dof_handler, ExcDoFHandlerRequired());
1212 AssertThrow(!triangulation_changed, ExcDoFHandlerChanged());
1213
1214 Vector<double> dof_vector(dof_handler->n_dofs());
1215
1216 typename std::vector<
1218 point = point_geometry_data.begin();
1219 for (; point != point_geometry_data.end(); ++point)
1220 {
1221 for (unsigned int component = 0;
1222 component < dof_handler->get_fe(0).n_components();
1223 component++)
1224 {
1225 dof_vector(point->solution_indices[component]) = 1;
1226 }
1227 }
1228 return dof_vector;
1229}
1230
1231
1232template <int dim>
1233void
1235 std::vector<std::vector<Point<dim>>> &locations)
1236{
1237 AssertThrow(!cleared, ExcInvalidState());
1238 AssertThrow(have_dof_handler, ExcDoFHandlerRequired());
1239 AssertThrow(!triangulation_changed, ExcDoFHandlerChanged());
1240
1241 std::vector<std::vector<Point<dim>>> actual_points;
1242 typename std::vector<
1244 point = point_geometry_data.begin();
1245
1246 for (; point != point_geometry_data.end(); ++point)
1247 {
1248 actual_points.push_back(point->support_point_locations);
1249 }
1250 locations = actual_points;
1251}
1252
1253
1254
1255template <int dim>
1256void
1258 const Quadrature<dim> &quadrature,
1259 std::vector<Point<dim>> &locations)
1260{
1261 Assert(!cleared, ExcInvalidState());
1262 AssertThrow(have_dof_handler, ExcDoFHandlerRequired());
1263
1264 locations = std::vector<Point<dim>>();
1265
1266 FEValues<dim> fe_values(dof_handler->get_fe(),
1267 quadrature,
1269 unsigned int n_quadrature_points = quadrature.size();
1270 std::vector<Point<dim>> evaluation_points;
1271
1272 // Loop over points and find correct cell
1273 typename std::vector<
1275 point = point_geometry_data.begin();
1276 Assert(!dof_handler->get_triangulation().is_mixed_mesh(),
1278 const auto reference_cell =
1279 dof_handler->get_triangulation().get_reference_cells()[0];
1280 for (unsigned int data_store_index = 0; point != point_geometry_data.end();
1281 ++point, ++data_store_index)
1282 {
1283 // we now have a point to query,
1284 // need to know what cell it is in
1285 Point<dim> requested_location = point->requested_location;
1288 reference_cell.template get_default_linear_mapping<dim, dim>(),
1289 *dof_handler,
1290 requested_location)
1291 .first;
1292 fe_values.reinit(cell);
1293
1294 evaluation_points = fe_values.get_quadrature_points();
1295 double distance = cell->diameter();
1296 unsigned int selected_point = 0;
1297
1298 for (unsigned int q_point = 0; q_point < n_quadrature_points; ++q_point)
1299 {
1300 if (requested_location.distance(evaluation_points[q_point]) <
1301 distance)
1302 {
1303 selected_point = q_point;
1304 distance =
1305 requested_location.distance(evaluation_points[q_point]);
1306 }
1307 }
1308
1309 locations.push_back(evaluation_points[selected_point]);
1310 }
1311}
1312
1313
1314template <int dim>
1315void
1317{
1318 out << "***PointValueHistory status output***\n\n";
1319 out << "Closed: " << closed << '\n';
1320 out << "Cleared: " << cleared << '\n';
1321 out << "Triangulation_changed: " << triangulation_changed << '\n';
1322 out << "Have_dof_handler: " << have_dof_handler << '\n';
1323 out << "Geometric Data" << '\n';
1324
1325 typename std::vector<
1327 point = point_geometry_data.begin();
1328 if (point == point_geometry_data.end())
1329 {
1330 out << "No points stored currently\n";
1331 }
1332 else
1333 {
1334 if (!cleared)
1335 {
1336 for (; point != point_geometry_data.end(); ++point)
1337 {
1338 out << "# Requested location: " << point->requested_location
1339 << '\n';
1340 out << "# DoF_index : Support location (for each component)\n";
1341 for (unsigned int component = 0;
1342 component < dof_handler->get_fe(0).n_components();
1343 component++)
1344 {
1345 out << point->solution_indices[component] << " : "
1346 << point->support_point_locations[component] << '\n';
1347 }
1348 out << '\n';
1349 }
1350 }
1351 else
1352 {
1353 out << "#Cannot access DoF_indices once cleared\n";
1354 }
1355 }
1356 out << '\n';
1357
1358 if (independent_values.size() != 0)
1359 {
1360 out << "Independent value(s): " << independent_values.size() << " : "
1361 << independent_values[0].size() << '\n';
1362 if (indep_names.size() > 0)
1363 {
1364 out << "Names: ";
1365 for (const auto &indep_name : indep_names)
1366 {
1367 out << "<" << indep_name << "> ";
1368 }
1369 out << '\n';
1370 }
1371 }
1372 else
1373 {
1374 out << "No independent values stored\n";
1375 }
1376
1377 if (data_store.begin() != data_store.end())
1378 {
1379 out
1380 << "Mnemonic: data set size (mask size, n true components) : n data sets\n";
1381 }
1382 for (const auto &data_entry : data_store)
1383 {
1384 // Find field mnemonic
1385 std::string vector_name = data_entry.first;
1386 typename std::map<std::string, ComponentMask>::iterator mask =
1387 component_mask.find(vector_name);
1388 Assert(mask != component_mask.end(),
1389 ExcMessage("vector_name not in class"));
1390 typename std::map<std::string, std::vector<std::string>>::iterator
1391 component_names = component_names_map.find(vector_name);
1392 Assert(component_names != component_names_map.end(),
1393 ExcMessage("vector_name not in class"));
1394
1395 if (data_entry.second.size() != 0)
1396 {
1397 out << data_entry.first << ": " << data_entry.second.size() << " (";
1398 out << mask->second.size() << ", "
1399 << mask->second.n_selected_components() << ") : ";
1400 out << (data_entry.second)[0].size() << '\n';
1401 }
1402 else
1403 {
1404 out << data_entry.first << ": " << data_entry.second.size() << " (";
1405 out << mask->second.size() << ", "
1406 << mask->second.n_selected_components() << ") : ";
1407 out << "No points added" << '\n';
1408 }
1409 // add names, if available
1410 if (component_names->second.size() > 0)
1411 {
1412 for (const auto &name : component_names->second)
1413 {
1414 out << "<" << name << "> ";
1415 }
1416 out << '\n';
1417 }
1418 }
1419 out << '\n';
1420 out << "***end of status output***\n\n";
1421}
1422
1423
1424
1425template <int dim>
1426bool
1428{
1429 // test ways that it can fail, if control
1430 // reaches last statement return true
1431 if (strict)
1432 {
1433 if (n_indep != 0)
1434 {
1435 if (dataset_key.size() != independent_values[0].size())
1436 {
1437 return false;
1438 }
1439 }
1440 if (have_dof_handler)
1441 {
1442 for (const auto &data_entry : data_store)
1443 {
1444 Assert(data_entry.second.size() > 0, ExcInternalError());
1445 if ((data_entry.second)[0].size() != dataset_key.size())
1446 return false;
1447 // this loop only tests one
1448 // member for each name,
1449 // i.e. checks the user it will
1450 // not catch internal errors
1451 // which do not update all
1452 // fields for a name.
1453 }
1454 }
1455 return true;
1456 }
1457 if (n_indep != 0)
1458 {
1459 if (std::abs(static_cast<int>(dataset_key.size()) -
1460 static_cast<int>(independent_values[0].size())) >= 2)
1461 {
1462 return false;
1463 }
1464 }
1465
1466 if (have_dof_handler)
1467 {
1468 for (const auto &data_entry : data_store)
1469 {
1470 Assert(data_entry.second.size() > 0, ExcInternalError());
1471
1472 if (std::abs(static_cast<int>((data_entry.second)[0].size()) -
1473 static_cast<int>(dataset_key.size())) >= 2)
1474 return false;
1475 // this loop only tests one member
1476 // for each name, i.e. checks the
1477 // user it will not catch internal
1478 // errors which do not update all
1479 // fields for a name.
1480 }
1481 }
1482 return true;
1483}
1484
1485
1486
1487template <int dim>
1488void
1490{
1491 // this function is called by the
1492 // Triangulation whenever something
1493 // changes, by virtue of having
1494 // attached the function to the
1495 // signal handler in the
1496 // triangulation object
1497
1498 // we record the fact that the mesh
1499 // has changed. we need to take
1500 // this into account next time we
1501 // evaluate the solution
1502 triangulation_changed = true;
1503}
1504
1505
1506// explicit instantiations
1507#include "numerics/point_value_history.inst"
1508
1509
unsigned int n_selected_components(const unsigned int overall_number_of_components=numbers::invalid_unsigned_int) const
virtual UpdateFlags get_needed_update_flags() const =0
virtual void evaluate_vector_field(const DataPostprocessorInputs::Vector< dim > &input_data, std::vector< Vector< double > > &computed_quantities) const
virtual void evaluate_scalar_field(const DataPostprocessorInputs::Scalar< dim > &input_data, std::vector< Vector< double > > &computed_quantities) const
virtual std::vector< std::string > get_names() const =0
void get_function_values(const ReadVector< Number > &fe_function, std::vector< Number > &values) const
const std::vector< Point< spacedim > > & get_quadrature_points() const
void get_function_hessians(const ReadVector< Number > &fe_function, std::vector< Tensor< 2, spacedim, Number > > &hessians) const
const Point< spacedim > & quadrature_point(const unsigned int q_point) const
void get_function_gradients(const ReadVector< Number > &fe_function, std::vector< Tensor< 1, spacedim, Number > > &gradients) const
void reinit(const TriaIterator< DoFCellAccessor< dim, spacedim, level_dof_access > > &cell)
void evaluate_field_at_requested_location(const std::string &name, const VectorType &solution)
void add_field_name(const std::string &vector_name, const ComponentMask &component_mask={})
boost::signals2::connection tria_listener
void add_point(const Point< dim > &location)
std::map< std::string, ComponentMask > component_mask
std::vector< internal::PointValueHistoryImplementation::PointGeometryData< dim > > point_geometry_data
std::map< std::string, std::vector< std::string > > component_names_map
PointValueHistory & operator=(const PointValueHistory &point_value_history)
ObserverPointer< const DoFHandler< dim >, PointValueHistory< dim > > dof_handler
void evaluate_field(const std::string &name, const VectorType &solution)
Vector< double > mark_support_locations()
std::vector< std::string > indep_names
std::vector< double > dataset_key
void write_gnuplot(const std::string &base_name, const std::vector< Point< dim > > &postprocessor_locations=std::vector< Point< dim > >())
void status(std::ostream &out)
void add_independent_names(const std::vector< std::string > &independent_names)
PointValueHistory(const unsigned int n_independent_variables=0)
void get_support_locations(std::vector< std::vector< Point< dim > > > &locations)
std::map< std::string, std::vector< std::vector< double > > > data_store
void add_component_names(const std::string &vector_name, const std::vector< std::string > &component_names)
void start_new_dataset(const double key)
void add_points(const std::vector< Point< dim > > &locations)
std::vector< std::vector< double > > independent_values
void push_back_independent(const std::vector< double > &independent_values)
void get_postprocessor_locations(const Quadrature< dim > &quadrature, std::vector< Point< dim > > &locations)
bool deep_check(const bool strict)
Definition point.h:113
numbers::NumberTraits< Number >::real_type distance(const Point< dim, Number > &p) const
unsigned int size() const
PointGeometryData(const Point< dim > &new_requested_location, const std::vector< Point< dim > > &new_locations, const std::vector< types::global_dof_index > &new_sol_indices)
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:35
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:36
Point< 2 > second
Definition grid_out.cc:4633
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcInvalidState()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
typename ActiveSelector::active_cell_iterator active_cell_iterator
UpdateFlags
@ update_hessians
Second derivatives of shape functions.
@ update_values
Shape function values.
@ update_normal_vectors
Normal vectors.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
std::size_t size
Definition mpi.cc:745
std::pair< typename MeshType< dim, spacedim >::active_cell_iterator, Point< dim > > find_active_cell_around_point(const Mapping< dim, spacedim > &mapping, const MeshType< dim, spacedim > &mesh, const Point< spacedim > &p, const std::vector< bool > &marked_vertices={}, const double tolerance=1.e-10)
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition utilities.cc:466
void point_value(const DoFHandler< dim, spacedim > &dof, const VectorType &fe_function, const Point< spacedim, double > &point, Vector< typename VectorType::value_type > &value)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
std::vector< Point< spacedim > > evaluation_points
std::vector< double > solution_values
std::vector< Tensor< 2, spacedim > > solution_hessians
std::vector< Tensor< 1, spacedim > > solution_gradients
std::vector< std::vector< Tensor< 2, spacedim > > > solution_hessians
std::vector<::Vector< double > > solution_values
std::vector< std::vector< Tensor< 1, spacedim > > > solution_gradients