Reference documentation for deal.II version 9.6.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\}}\)
Loading...
Searching...
No Matches
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 for (unsigned int component = 0;
294 component < dof_handler->get_fe(0).n_components();
295 component++)
296 {
297 new_solution_indices.push_back(
298 local_dof_indices[current_fe_index[component]]);
299 }
300
302 new_point_geometry_data(location, current_points, new_solution_indices);
303 point_geometry_data.push_back(new_point_geometry_data);
304
305 for (auto &data_entry : data_store)
306 {
307 // add an extra row to each vector entry
308 const ComponentMask &current_mask =
309 (component_mask.find(data_entry.first))->second;
310 unsigned int n_stored = current_mask.n_selected_components();
311 data_entry.second.resize(data_entry.second.size() + n_stored);
312 }
313}
314
315
316
317template <int dim>
318void
319PointValueHistory<dim>::add_points(const std::vector<Point<dim>> &locations)
320{
321 // This algorithm adds points in the same
322 // order as they appear in the vector
323 // locations and users may depend on this
324 // so do not change order added!
325
326 // can't be closed to add additional points or vectors
327 AssertThrow(!closed, ExcInvalidState());
328 AssertThrow(!cleared, ExcInvalidState());
329 AssertThrow(have_dof_handler, ExcDoFHandlerRequired());
330 AssertThrow(!triangulation_changed, ExcDoFHandlerChanged());
331
332
333 // Implementation assumes that support
334 // points locations are dofs locations
335 AssertThrow(dof_handler->get_fe().has_support_points(), ExcNotImplemented());
336
337 // While in general quadrature points seems
338 // to refer to Gauss quadrature points, in
339 // this case the quadrature points are
340 // forced to be the support points of the
341 // FE.
342 Quadrature<dim> support_point_quadrature(
343 dof_handler->get_fe().get_unit_support_points());
344 FEValues<dim> fe_values(dof_handler->get_fe(),
345 support_point_quadrature,
347 unsigned int n_support_points =
348 dof_handler->get_fe().get_unit_support_points().size();
349 unsigned int n_components = dof_handler->get_fe(0).n_components();
350
351 // set up a loop over all the cells in the
352 // DoFHandler
354 dof_handler->begin_active();
355 typename DoFHandler<dim>::active_cell_iterator endc = dof_handler->end();
356
357 // default values to be replaced as closer
358 // points are found however they need to be
359 // consistent in case they are actually
360 // chosen vector <vector>s defined where
361 // previously single vectors were used
362
363 // need to store one value per point per component
364 std::vector<typename DoFHandler<dim>::active_cell_iterator> current_cell(
365 locations.size(), cell);
366
367 fe_values.reinit(cell);
368 std::vector<Point<dim>> temp_points(n_components, Point<dim>());
369 std::vector<unsigned int> temp_fe_index(n_components, 0);
370 for (unsigned int support_point = 0; support_point < n_support_points;
371 support_point++)
372 {
373 // set up valid data in the empty vectors
374 unsigned int component =
375 dof_handler->get_fe().system_to_component_index(support_point).first;
376 temp_points[component] = fe_values.quadrature_point(support_point);
377 temp_fe_index[component] = support_point;
378 }
379 std::vector<std::vector<Point<dim>>> current_points(
380 locations.size(), temp_points); // give a valid start point
381 std::vector<std::vector<unsigned int>> current_fe_index(locations.size(),
382 temp_fe_index);
383
384 // check each cell to find suitable support
385 // points
386 // GridTools::find_active_cell_around_point
387 // is an alternative. That method is not
388 // used here mostly because of the history
389 // of the class. The algorithm used here
390 // may be slightly more
391 // efficient than find_active_cell_around_point
392 // because it operates on a set of points.
393 for (; cell != endc; ++cell)
394 {
395 fe_values.reinit(cell);
396 for (unsigned int support_point = 0; support_point < n_support_points;
397 support_point++)
398 {
399 unsigned int component = dof_handler->get_fe()
400 .system_to_component_index(support_point)
401 .first;
402 const Point<dim> &test_point =
403 fe_values.quadrature_point(support_point);
404
405 for (unsigned int point = 0; point < locations.size(); ++point)
406 {
407 if (locations[point].distance(test_point) <
408 locations[point].distance(current_points[point][component]))
409 {
410 // save the data
411 current_points[point][component] = test_point;
412 current_cell[point] = cell;
413 current_fe_index[point][component] = support_point;
414 }
415 }
416 }
417 }
418
419 std::vector<types::global_dof_index> local_dof_indices(
420 dof_handler->get_fe().n_dofs_per_cell());
421 for (unsigned int point = 0; point < locations.size(); ++point)
422 {
423 current_cell[point]->get_dof_indices(local_dof_indices);
424 std::vector<types::global_dof_index> new_solution_indices;
425
426 for (unsigned int component = 0;
427 component < dof_handler->get_fe(0).n_components();
428 component++)
429 {
430 new_solution_indices.push_back(
431 local_dof_indices[current_fe_index[point][component]]);
432 }
433
435 new_point_geometry_data(locations[point],
436 current_points[point],
437 new_solution_indices);
438
439 point_geometry_data.push_back(new_point_geometry_data);
440
441 for (auto &data_entry : data_store)
442 {
443 // add an extra row to each vector entry
444 const ComponentMask current_mask =
445 (component_mask.find(data_entry.first))->second;
446 unsigned int n_stored = current_mask.n_selected_components();
447 data_entry.second.resize(data_entry.second.size() + n_stored);
448 }
449 }
450}
451
452
453
454template <int dim>
455void
456PointValueHistory<dim>::add_field_name(const std::string &vector_name,
457 const ComponentMask &mask)
458{
459 // can't be closed to add additional points
460 // or vectors
461 AssertThrow(!closed, ExcInvalidState());
462 AssertThrow(!cleared, ExcInvalidState());
463 AssertThrow(have_dof_handler, ExcDoFHandlerRequired());
464 AssertThrow(!triangulation_changed, ExcDoFHandlerChanged());
465
466 // insert a component mask that is always of the right size
467 if (mask.represents_the_all_selected_mask() == false)
468 component_mask.insert(std::make_pair(vector_name, mask));
469 else
470 component_mask.insert(
471 std::make_pair(vector_name,
472 ComponentMask(std::vector<bool>(
473 dof_handler->get_fe(0).n_components(), true))));
474
475 // insert an empty vector of strings
476 // to ensure each field has an entry
477 // in the map
478 std::pair<std::string, std::vector<std::string>> empty_names(
479 vector_name, std::vector<std::string>());
480 component_names_map.insert(empty_names);
481
482 // make and add a new vector
483 // point_geometry_data.size() long
484 std::pair<std::string, std::vector<std::vector<double>>> pair_data;
485 pair_data.first = vector_name;
486 const unsigned int n_stored =
487 (mask.represents_the_all_selected_mask() == false ?
488 mask.n_selected_components() :
489 dof_handler->get_fe(0).n_components());
490
491 int n_datastreams =
492 point_geometry_data.size() * n_stored; // each point has n_stored sub parts
493 std::vector<std::vector<double>> vector_size(n_datastreams,
494 std::vector<double>(0));
495 pair_data.second = vector_size;
496 data_store.insert(pair_data);
497}
498
499
500template <int dim>
501void
502PointValueHistory<dim>::add_field_name(const std::string &vector_name,
503 const unsigned int n_components)
504{
505 ComponentMask temp_mask(std::vector<bool>(n_components, true));
506 add_field_name(vector_name, temp_mask);
507}
508
509
510template <int dim>
511void
513 const std::string &vector_name,
514 const std::vector<std::string> &component_names)
515{
516 typename std::map<std::string, std::vector<std::string>>::iterator names =
517 component_names_map.find(vector_name);
518 Assert(names != component_names_map.end(),
519 ExcMessage("vector_name not in class"));
520
521 typename std::map<std::string, ComponentMask>::iterator mask =
522 component_mask.find(vector_name);
523 Assert(mask != component_mask.end(), ExcMessage("vector_name not in class"));
524 unsigned int n_stored = mask->second.n_selected_components();
525 (void)n_stored;
526 Assert(component_names.size() == n_stored,
527 ExcDimensionMismatch(component_names.size(), n_stored));
528
529 names->second = component_names;
530}
531
532
533template <int dim>
534void
536 const std::vector<std::string> &independent_names)
537{
538 Assert(independent_names.size() == n_indep,
539 ExcDimensionMismatch(independent_names.size(), n_indep));
540
541 indep_names = independent_names;
542}
543
544
545template <int dim>
546void
548{
549 closed = true;
550}
551
552
553
554template <int dim>
555void
557{
558 cleared = true;
559 dof_handler = nullptr;
560 have_dof_handler = false;
561}
562
563// Need to test that the internal data has a full and complete dataset for
564// each key. That is that the data has not got 'out of sync'. Testing that
565// dataset_key is within 1 of independent_values is cheap and is done in all
566// three methods. Evaluate_field will check that its vector_name is within 1
567// of dataset_key. However this leaves the possibility that the user has
568// neglected to call evaluate_field on one vector_name consistently. To catch
569// this case start_new_dataset will call bool deap_check () which will test
570// all vector_names and return a bool. This can be called from an Assert
571// statement.
572
573
574
575template <int dim>
576template <typename VectorType>
577void
578PointValueHistory<dim>::evaluate_field(const std::string &vector_name,
579 const VectorType &solution)
580{
581 // must be closed to add data to internal
582 // members.
583 Assert(closed, ExcInvalidState());
584 Assert(!cleared, ExcInvalidState());
585 AssertThrow(have_dof_handler, ExcDoFHandlerRequired());
586 AssertThrow(!triangulation_changed, ExcDoFHandlerChanged());
587
588 if (n_indep != 0) // hopefully this will get optimized, can't test
589 // independent_values[0] unless n_indep > 0
590 {
591 Assert(std::abs(static_cast<int>(dataset_key.size()) -
592 static_cast<int>(independent_values[0].size())) < 2,
593 ExcDataLostSync());
594 }
595 // Look up the field name and get an
596 // iterator for the map. Doing this
597 // up front means that it only needs
598 // to be done once and also allows us
599 // to check vector_name is in the map.
600 typename std::map<std::string, std::vector<std::vector<double>>>::iterator
601 data_store_field = data_store.find(vector_name);
602 Assert(data_store_field != data_store.end(),
603 ExcMessage("vector_name not in class"));
604 // Repeat for component_mask
605 typename std::map<std::string, ComponentMask>::iterator mask =
606 component_mask.find(vector_name);
607 Assert(mask != component_mask.end(), ExcMessage("vector_name not in class"));
608
609 unsigned int n_stored =
610 mask->second.n_selected_components(dof_handler->get_fe(0).n_components());
611
612 typename std::vector<
614 point = point_geometry_data.begin();
615 for (unsigned int data_store_index = 0; point != point_geometry_data.end();
616 ++point, ++data_store_index)
617 {
618 // Look up the components to add
619 // in the component_mask, and
620 // access the data associated with
621 // those components
622
623 for (unsigned int store_index = 0, comp = 0;
624 comp < dof_handler->get_fe(0).n_components();
625 comp++)
626 {
627 if (mask->second[comp])
628 {
629 unsigned int solution_index = point->solution_indices[comp];
630 data_store_field
631 ->second[data_store_index * n_stored + store_index]
632 .push_back(
634 solution_index));
635 ++store_index;
636 }
637 }
638 }
639}
640
641
642
643template <int dim>
644template <typename VectorType>
645void
647 const std::vector<std::string> &vector_names,
648 const VectorType &solution,
649 const DataPostprocessor<dim> &data_postprocessor,
650 const Quadrature<dim> &quadrature)
651{
652 // must be closed to add data to internal
653 // members.
654 Assert(closed, ExcInvalidState());
655 Assert(!cleared, ExcInvalidState());
656 AssertThrow(have_dof_handler, ExcDoFHandlerRequired());
657 if (n_indep != 0) // hopefully this will get optimized, can't test
658 // independent_values[0] unless n_indep > 0
659 {
660 Assert(std::abs(static_cast<int>(dataset_key.size()) -
661 static_cast<int>(independent_values[0].size())) < 2,
662 ExcDataLostSync());
663 }
664
665 // Make an FEValues object
666 const UpdateFlags update_flags =
668 Assert(
669 !(update_flags & update_normal_vectors),
671 "The update of normal vectors may not be requested for evaluation of "
672 "data on cells via DataPostprocessor."));
673 FEValues<dim> fe_values(dof_handler->get_fe(), quadrature, update_flags);
674 unsigned int n_components = dof_handler->get_fe(0).n_components();
675 unsigned int n_quadrature_points = quadrature.size();
676
677 unsigned int n_output_variables = data_postprocessor.get_names().size();
678
679 // declare some temp objects for evaluating the solution at quadrature
680 // points. we will either need the scalar or vector version in the code
681 // below
682 std::vector<typename VectorType::value_type> scalar_solution_values(
683 n_quadrature_points);
684 std::vector<Tensor<1, dim, typename VectorType::value_type>>
685 scalar_solution_gradients(n_quadrature_points);
686 std::vector<Tensor<2, dim, typename VectorType::value_type>>
687 scalar_solution_hessians(n_quadrature_points);
688
689 std::vector<Vector<typename VectorType::value_type>> vector_solution_values(
690 n_quadrature_points, Vector<typename VectorType::value_type>(n_components));
691
692 std::vector<std::vector<Tensor<1, dim, typename VectorType::value_type>>>
693 vector_solution_gradients(
694 n_quadrature_points,
697
698 std::vector<std::vector<Tensor<2, dim, typename VectorType::value_type>>>
699 vector_solution_hessians(
700 n_quadrature_points,
703
704 // Loop over points and find correct cell
705 typename std::vector<
707 point = point_geometry_data.begin();
708 Assert(!dof_handler->get_triangulation().is_mixed_mesh(),
710 const auto reference_cell =
711 dof_handler->get_triangulation().get_reference_cells()[0];
712 for (unsigned int data_store_index = 0; point != point_geometry_data.end();
713 ++point, ++data_store_index)
714 {
715 // we now have a point to query, need to know what cell it is in
716 const Point<dim> requested_location = point->requested_location;
717 const typename DoFHandler<dim>::active_cell_iterator cell =
719 reference_cell.template get_default_linear_mapping<dim, dim>(),
720 *dof_handler,
721 requested_location)
722 .first;
723
724
725 fe_values.reinit(cell);
726 std::vector<Vector<double>> computed_quantities(
727 1, Vector<double>(n_output_variables)); // just one point needed
728
729 // find the closest quadrature point
730 std::vector<Point<dim>> quadrature_points =
731 fe_values.get_quadrature_points();
732 double distance = cell->diameter();
733 unsigned int selected_point = 0;
734 for (unsigned int q_point = 0; q_point < n_quadrature_points; ++q_point)
735 {
736 if (requested_location.distance(quadrature_points[q_point]) <
737 distance)
738 {
739 selected_point = q_point;
740 distance =
741 requested_location.distance(quadrature_points[q_point]);
742 }
743 }
744
745
746 // The case of a scalar FE
747 if (n_components == 1)
748 {
749 // Extract data for the DataPostprocessor object
750 DataPostprocessorInputs::Scalar<dim> postprocessor_input;
751
752 // for each quantity that is requested (values, gradients, hessians),
753 // first get them at all quadrature points, then restrict to the
754 // one value on the quadrature point closest to the evaluation
755 // point in question
756 //
757 // we need temporary objects because the underlying scalar
758 // type of the solution vector may be different from 'double',
759 // but the DataPostprocessorInputs only allow for 'double'
760 // data types
761 if (update_flags & update_values)
762 {
763 fe_values.get_function_values(solution, scalar_solution_values);
764 postprocessor_input.solution_values =
765 std::vector<double>(1, scalar_solution_values[selected_point]);
766 }
767 if (update_flags & update_gradients)
768 {
769 fe_values.get_function_gradients(solution,
770 scalar_solution_gradients);
771 postprocessor_input.solution_gradients =
772 std::vector<Tensor<1, dim>>(
773 1, scalar_solution_gradients[selected_point]);
774 }
775 if (update_flags & update_hessians)
776 {
777 fe_values.get_function_hessians(solution,
778 scalar_solution_hessians);
779 postprocessor_input.solution_hessians =
780 std::vector<Tensor<2, dim>>(
781 1, scalar_solution_hessians[selected_point]);
782 }
783
784 // then also set the single evaluation point
785 postprocessor_input.evaluation_points =
786 std::vector<Point<dim>>(1, quadrature_points[selected_point]);
787
788 // and finally do the postprocessing
789 data_postprocessor.evaluate_scalar_field(postprocessor_input,
790 computed_quantities);
791 }
792 else // The case of a vector FE
793 {
794 // exact same idea as above
795 DataPostprocessorInputs::Vector<dim> postprocessor_input;
796
797 if (update_flags & update_values)
798 {
799 fe_values.get_function_values(solution, vector_solution_values);
800 postprocessor_input.solution_values.resize(
801 1, Vector<double>(n_components));
802 std::copy(vector_solution_values[selected_point].begin(),
803 vector_solution_values[selected_point].end(),
804 postprocessor_input.solution_values[0].begin());
805 }
806 if (update_flags & update_gradients)
807 {
808 fe_values.get_function_gradients(solution,
809 vector_solution_gradients);
810 postprocessor_input.solution_gradients.resize(
811 1, std::vector<Tensor<1, dim>>(n_components));
812 std::copy(vector_solution_gradients[selected_point].begin(),
813 vector_solution_gradients[selected_point].end(),
814 postprocessor_input.solution_gradients[0].begin());
815 }
816 if (update_flags & update_hessians)
817 {
818 fe_values.get_function_hessians(solution,
819 vector_solution_hessians);
820 postprocessor_input.solution_hessians.resize(
821 1, std::vector<Tensor<2, dim>>(n_components));
822 std::copy(vector_solution_hessians[selected_point].begin(),
823 vector_solution_hessians[selected_point].end(),
824 postprocessor_input.solution_hessians[0].begin());
825 }
826
827 postprocessor_input.evaluation_points =
828 std::vector<Point<dim>>(1, quadrature_points[selected_point]);
829
830 data_postprocessor.evaluate_vector_field(postprocessor_input,
831 computed_quantities);
832 }
833
834
835 // we now have the data and need to save it
836 // loop over data names
837 typename std::vector<std::string>::const_iterator name =
838 vector_names.begin();
839 for (; name != vector_names.end(); ++name)
840 {
841 typename std::map<std::string,
842 std::vector<std::vector<double>>>::iterator
843 data_store_field = data_store.find(*name);
844 Assert(data_store_field != data_store.end(),
845 ExcMessage("vector_name not in class"));
846 // Repeat for component_mask
847 typename std::map<std::string, ComponentMask>::iterator mask =
848 component_mask.find(*name);
849 Assert(mask != component_mask.end(),
850 ExcMessage("vector_name not in class"));
851
852 unsigned int n_stored =
853 mask->second.n_selected_components(n_output_variables);
854
855 // Push back computed quantities according
856 // to the component_mask.
857 for (unsigned int store_index = 0, comp = 0;
858 comp < n_output_variables;
859 comp++)
860 {
861 if (mask->second[comp])
862 {
863 data_store_field
864 ->second[data_store_index * n_stored + store_index]
865 .push_back(computed_quantities[0](comp));
866 ++store_index;
867 }
868 }
869 }
870 } // end of loop over points
871}
872
873
874
875template <int dim>
876template <typename VectorType>
877void
879 const std::string &vector_name,
880 const VectorType &solution,
881 const DataPostprocessor<dim> &data_postprocessor,
882 const Quadrature<dim> &quadrature)
883{
884 std::vector<std::string> vector_names;
885 vector_names.push_back(vector_name);
886 evaluate_field(vector_names, solution, data_postprocessor, quadrature);
887}
888
889
890
891template <int dim>
892template <typename VectorType>
893void
895 const std::string &vector_name,
896 const VectorType &solution)
897{
898 using number = typename VectorType::value_type;
899 // must be closed to add data to internal
900 // members.
901 Assert(closed, ExcInvalidState());
902 Assert(!cleared, ExcInvalidState());
903 AssertThrow(have_dof_handler, ExcDoFHandlerRequired());
904
905 if (n_indep != 0) // hopefully this will get optimized, can't test
906 // independent_values[0] unless n_indep > 0
907 {
908 Assert(std::abs(static_cast<int>(dataset_key.size()) -
909 static_cast<int>(independent_values[0].size())) < 2,
910 ExcDataLostSync());
911 }
912 // Look up the field name and get an
913 // iterator for the map. Doing this
914 // up front means that it only needs
915 // to be done once and also allows us
916 // to check vector_name is in the map.
917 typename std::map<std::string, std::vector<std::vector<double>>>::iterator
918 data_store_field = data_store.find(vector_name);
919 Assert(data_store_field != data_store.end(),
920 ExcMessage("vector_name not in class"));
921 // Repeat for component_mask
922 typename std::map<std::string, ComponentMask>::iterator mask =
923 component_mask.find(vector_name);
924 Assert(mask != component_mask.end(), ExcMessage("vector_name not in class"));
925
926 unsigned int n_stored =
927 mask->second.n_selected_components(dof_handler->get_fe(0).n_components());
928
929 typename std::vector<
931 point = point_geometry_data.begin();
932 Vector<number> value(dof_handler->get_fe(0).n_components());
933 for (unsigned int data_store_index = 0; point != point_geometry_data.end();
934 ++point, ++data_store_index)
935 {
936 // Make a Vector <double> for the value
937 // at the point. It will have as many
938 // components as there are in the fe.
939 VectorTools::point_value(*dof_handler,
940 solution,
941 point->requested_location,
942 value);
943
944 // Look up the component_mask and add
945 // in components according to that mask
946 for (unsigned int store_index = 0, comp = 0; comp < mask->second.size();
947 comp++)
948 {
949 if (mask->second[comp])
950 {
951 data_store_field
952 ->second[data_store_index * n_stored + store_index]
953 .push_back(value(comp));
954 ++store_index;
955 }
956 }
957 }
958}
959
960
961template <int dim>
962void
964{
965 // must be closed to add data to internal
966 // members.
967 Assert(closed, ExcInvalidState());
968 Assert(!cleared, ExcInvalidState());
969 Assert(deep_check(false), ExcDataLostSync());
970
971 dataset_key.push_back(key);
972}
973
974
975
976template <int dim>
977void
979 const std::vector<double> &indep_values)
980{
981 // must be closed to add data to internal
982 // members.
983 Assert(closed, ExcInvalidState());
984 Assert(!cleared, ExcInvalidState());
985 Assert(indep_values.size() == n_indep,
986 ExcDimensionMismatch(indep_values.size(), n_indep));
987 Assert(n_indep != 0, ExcNoIndependent());
988 Assert(std::abs(static_cast<int>(dataset_key.size()) -
989 static_cast<int>(independent_values[0].size())) < 2,
990 ExcDataLostSync());
991
992 for (unsigned int component = 0; component < n_indep; ++component)
993 independent_values[component].push_back(indep_values[component]);
994}
995
996
997
998template <int dim>
999void
1001 const std::string &base_name,
1002 const std::vector<Point<dim>> &postprocessor_locations)
1003{
1004 AssertThrow(closed, ExcInvalidState());
1005 AssertThrow(!cleared, ExcInvalidState());
1006 AssertThrow(deep_check(true), ExcDataLostSync());
1007
1008 // write inputs to a file
1009 if (n_indep != 0)
1010 {
1011 std::string filename = base_name + "_indep.gpl";
1012 std::ofstream to_gnuplot(filename);
1013
1014 to_gnuplot << "# Data independent of mesh location\n";
1015
1016 // write column headings
1017 to_gnuplot << "# <Key> ";
1018
1019 if (indep_names.size() > 0)
1020 {
1021 for (const auto &indep_name : indep_names)
1022 {
1023 to_gnuplot << "<" << indep_name << "> ";
1024 }
1025 to_gnuplot << '\n';
1026 }
1027 else
1028 {
1029 for (unsigned int component = 0; component < n_indep; ++component)
1030 {
1031 to_gnuplot << "<Indep_" << component << "> ";
1032 }
1033 to_gnuplot << '\n';
1034 }
1035 // write general data stored
1036 for (unsigned int key = 0; key < dataset_key.size(); ++key)
1037 {
1038 to_gnuplot << dataset_key[key];
1039
1040 for (unsigned int component = 0; component < n_indep; ++component)
1041 {
1042 to_gnuplot << " " << independent_values[component][key];
1043 }
1044 to_gnuplot << '\n';
1045 }
1046
1047 to_gnuplot.close();
1048 }
1049
1050
1051
1052 // write points to a file
1053 if (have_dof_handler)
1054 {
1055 AssertThrow(have_dof_handler, ExcDoFHandlerRequired());
1056 AssertThrow(postprocessor_locations.empty() ||
1057 postprocessor_locations.size() ==
1058 point_geometry_data.size(),
1059 ExcDimensionMismatch(postprocessor_locations.size(),
1060 point_geometry_data.size()));
1061 // We previously required the
1062 // number of dofs to remain the
1063 // same to provide some sort of
1064 // test on the relevance of the
1065 // support point indices stored.
1066 // We now relax that to allow
1067 // adaptive refinement strategies
1068 // to make use of the
1069 // evaluate_field_requested_locations
1070 // method. Note that the support point
1071 // information is not meaningful if
1072 // the number of dofs has changed.
1073 // AssertThrow (!triangulation_changed, ExcDoFHandlerChanged ());
1074
1075 typename std::vector<internal::PointValueHistoryImplementation::
1076 PointGeometryData<dim>>::iterator point =
1077 point_geometry_data.begin();
1078 for (unsigned int data_store_index = 0;
1079 point != point_geometry_data.end();
1080 ++point, ++data_store_index)
1081 {
1082 // for each point, open a file to
1083 // be written to
1084 std::string filename = base_name + "_" +
1085 Utilities::int_to_string(data_store_index, 2) +
1086 ".gpl"; // store by order pushed back
1087 // due to
1088 // Utilities::int_to_string(data_store_index,
1089 // 2) call, can handle up to 100
1090 // points
1091 std::ofstream to_gnuplot(filename);
1092
1093 // put helpful info about the
1094 // support point into the file as
1095 // comments
1096 to_gnuplot << "# Requested location: " << point->requested_location
1097 << '\n';
1098 to_gnuplot << "# DoF_index : Support location (for each component)\n";
1099 for (unsigned int component = 0;
1100 component < dof_handler->get_fe(0).n_components();
1101 component++)
1102 {
1103 to_gnuplot << "# " << point->solution_indices[component] << " : "
1104 << point->support_point_locations[component] << '\n';
1105 }
1106 if (triangulation_changed)
1107 to_gnuplot
1108 << "# (Original components and locations, may be invalidated by mesh change.)\n";
1109
1110 if (postprocessor_locations.size() != 0)
1111 {
1112 to_gnuplot << "# Postprocessor location: "
1113 << postprocessor_locations[data_store_index];
1114 if (triangulation_changed)
1115 to_gnuplot << " (may be approximate)\n";
1116 }
1117 to_gnuplot << "#\n";
1118
1119
1120 // write column headings
1121 to_gnuplot << "# <Key> ";
1122
1123 if (indep_names.size() > 0)
1124 {
1125 for (const auto &indep_name : indep_names)
1126 {
1127 to_gnuplot << "<" << indep_name << "> ";
1128 }
1129 }
1130 else
1131 {
1132 for (unsigned int component = 0; component < n_indep; ++component)
1133 {
1134 to_gnuplot << "<Indep_" << component << "> ";
1135 }
1136 }
1137
1138 for (const auto &data_entry : data_store)
1139 {
1140 typename std::map<std::string, ComponentMask>::iterator mask =
1141 component_mask.find(data_entry.first);
1142 unsigned int n_stored = mask->second.n_selected_components();
1143 std::vector<std::string> names =
1144 (component_names_map.find(data_entry.first))->second;
1145
1146 if (names.size() > 0)
1147 {
1148 AssertThrow(names.size() == n_stored,
1149 ExcDimensionMismatch(names.size(), n_stored));
1150 for (const auto &name : names)
1151 {
1152 to_gnuplot << "<" << name << "> ";
1153 }
1154 }
1155 else
1156 {
1157 for (unsigned int component = 0; component < n_stored;
1158 component++)
1159 {
1160 to_gnuplot << "<" << data_entry.first << "_" << component
1161 << "> ";
1162 }
1163 }
1164 }
1165 to_gnuplot << '\n';
1166
1167 // write data stored for the point
1168 for (unsigned int key = 0; key < dataset_key.size(); ++key)
1169 {
1170 to_gnuplot << dataset_key[key];
1171
1172 for (unsigned int component = 0; component < n_indep; ++component)
1173 {
1174 to_gnuplot << " " << independent_values[component][key];
1175 }
1176
1177 for (const auto &data_entry : data_store)
1178 {
1179 typename std::map<std::string, ComponentMask>::iterator mask =
1180 component_mask.find(data_entry.first);
1181 unsigned int n_stored = mask->second.n_selected_components();
1182
1183 for (unsigned int component = 0; component < n_stored;
1184 component++)
1185 {
1186 to_gnuplot
1187 << " "
1188 << (data_entry.second)[data_store_index * n_stored +
1189 component][key];
1190 }
1191 }
1192 to_gnuplot << '\n';
1193 }
1194
1195 to_gnuplot.close();
1196 }
1197 }
1198}
1199
1200
1201
1202template <int dim>
1205{
1206 // a method to put a one at each point on
1207 // the grid where a location is defined
1208 AssertThrow(!cleared, ExcInvalidState());
1209 AssertThrow(have_dof_handler, ExcDoFHandlerRequired());
1210 AssertThrow(!triangulation_changed, ExcDoFHandlerChanged());
1211
1212 Vector<double> dof_vector(dof_handler->n_dofs());
1213
1214 typename std::vector<
1216 point = point_geometry_data.begin();
1217 for (; point != point_geometry_data.end(); ++point)
1218 {
1219 for (unsigned int component = 0;
1220 component < dof_handler->get_fe(0).n_components();
1221 component++)
1222 {
1223 dof_vector(point->solution_indices[component]) = 1;
1224 }
1225 }
1226 return dof_vector;
1227}
1228
1229
1230template <int dim>
1231void
1233 std::vector<std::vector<Point<dim>>> &locations)
1234{
1235 AssertThrow(!cleared, ExcInvalidState());
1236 AssertThrow(have_dof_handler, ExcDoFHandlerRequired());
1237 AssertThrow(!triangulation_changed, ExcDoFHandlerChanged());
1238
1239 std::vector<std::vector<Point<dim>>> actual_points;
1240 typename std::vector<
1242 point = point_geometry_data.begin();
1243
1244 for (; point != point_geometry_data.end(); ++point)
1245 {
1246 actual_points.push_back(point->support_point_locations);
1247 }
1248 locations = actual_points;
1249}
1250
1251
1252
1253template <int dim>
1254void
1256 const Quadrature<dim> &quadrature,
1257 std::vector<Point<dim>> &locations)
1258{
1259 Assert(!cleared, ExcInvalidState());
1260 AssertThrow(have_dof_handler, ExcDoFHandlerRequired());
1261
1262 locations = std::vector<Point<dim>>();
1263
1264 FEValues<dim> fe_values(dof_handler->get_fe(),
1265 quadrature,
1267 unsigned int n_quadrature_points = quadrature.size();
1268 std::vector<Point<dim>> evaluation_points;
1269
1270 // Loop over points and find correct cell
1271 typename std::vector<
1273 point = point_geometry_data.begin();
1274 Assert(!dof_handler->get_triangulation().is_mixed_mesh(),
1276 const auto reference_cell =
1277 dof_handler->get_triangulation().get_reference_cells()[0];
1278 for (unsigned int data_store_index = 0; point != point_geometry_data.end();
1279 ++point, ++data_store_index)
1280 {
1281 // we now have a point to query,
1282 // need to know what cell it is in
1283 Point<dim> requested_location = point->requested_location;
1286 reference_cell.template get_default_linear_mapping<dim, dim>(),
1287 *dof_handler,
1288 requested_location)
1289 .first;
1290 fe_values.reinit(cell);
1291
1292 evaluation_points = fe_values.get_quadrature_points();
1293 double distance = cell->diameter();
1294 unsigned int selected_point = 0;
1295
1296 for (unsigned int q_point = 0; q_point < n_quadrature_points; ++q_point)
1297 {
1298 if (requested_location.distance(evaluation_points[q_point]) <
1299 distance)
1300 {
1301 selected_point = q_point;
1302 distance =
1303 requested_location.distance(evaluation_points[q_point]);
1304 }
1305 }
1306
1307 locations.push_back(evaluation_points[selected_point]);
1308 }
1309}
1310
1311
1312template <int dim>
1313void
1315{
1316 out << "***PointValueHistory status output***\n\n";
1317 out << "Closed: " << closed << '\n';
1318 out << "Cleared: " << cleared << '\n';
1319 out << "Triangulation_changed: " << triangulation_changed << '\n';
1320 out << "Have_dof_handler: " << have_dof_handler << '\n';
1321 out << "Geometric Data" << '\n';
1322
1323 typename std::vector<
1325 point = point_geometry_data.begin();
1326 if (point == point_geometry_data.end())
1327 {
1328 out << "No points stored currently\n";
1329 }
1330 else
1331 {
1332 if (!cleared)
1333 {
1334 for (; point != point_geometry_data.end(); ++point)
1335 {
1336 out << "# Requested location: " << point->requested_location
1337 << '\n';
1338 out << "# DoF_index : Support location (for each component)\n";
1339 for (unsigned int component = 0;
1340 component < dof_handler->get_fe(0).n_components();
1341 component++)
1342 {
1343 out << point->solution_indices[component] << " : "
1344 << point->support_point_locations[component] << '\n';
1345 }
1346 out << '\n';
1347 }
1348 }
1349 else
1350 {
1351 out << "#Cannot access DoF_indices once cleared\n";
1352 }
1353 }
1354 out << '\n';
1355
1356 if (independent_values.size() != 0)
1357 {
1358 out << "Independent value(s): " << independent_values.size() << " : "
1359 << independent_values[0].size() << '\n';
1360 if (indep_names.size() > 0)
1361 {
1362 out << "Names: ";
1363 for (const auto &indep_name : indep_names)
1364 {
1365 out << "<" << indep_name << "> ";
1366 }
1367 out << '\n';
1368 }
1369 }
1370 else
1371 {
1372 out << "No independent values stored\n";
1373 }
1374
1375 if (data_store.begin() != data_store.end())
1376 {
1377 out
1378 << "Mnemonic: data set size (mask size, n true components) : n data sets\n";
1379 }
1380 for (const auto &data_entry : data_store)
1381 {
1382 // Find field mnemonic
1383 std::string vector_name = data_entry.first;
1384 typename std::map<std::string, ComponentMask>::iterator mask =
1385 component_mask.find(vector_name);
1386 Assert(mask != component_mask.end(),
1387 ExcMessage("vector_name not in class"));
1388 typename std::map<std::string, std::vector<std::string>>::iterator
1389 component_names = component_names_map.find(vector_name);
1390 Assert(component_names != component_names_map.end(),
1391 ExcMessage("vector_name not in class"));
1392
1393 if (data_entry.second.size() != 0)
1394 {
1395 out << data_entry.first << ": " << data_entry.second.size() << " (";
1396 out << mask->second.size() << ", "
1397 << mask->second.n_selected_components() << ") : ";
1398 out << (data_entry.second)[0].size() << '\n';
1399 }
1400 else
1401 {
1402 out << data_entry.first << ": " << data_entry.second.size() << " (";
1403 out << mask->second.size() << ", "
1404 << mask->second.n_selected_components() << ") : ";
1405 out << "No points added" << '\n';
1406 }
1407 // add names, if available
1408 if (component_names->second.size() > 0)
1409 {
1410 for (const auto &name : component_names->second)
1411 {
1412 out << "<" << name << "> ";
1413 }
1414 out << '\n';
1415 }
1416 }
1417 out << '\n';
1418 out << "***end of status output***\n\n";
1419}
1420
1421
1422
1423template <int dim>
1424bool
1426{
1427 // test ways that it can fail, if control
1428 // reaches last statement return true
1429 if (strict)
1430 {
1431 if (n_indep != 0)
1432 {
1433 if (dataset_key.size() != independent_values[0].size())
1434 {
1435 return false;
1436 }
1437 }
1438 if (have_dof_handler)
1439 {
1440 for (const auto &data_entry : data_store)
1441 {
1442 Assert(data_entry.second.size() > 0, ExcInternalError());
1443 if ((data_entry.second)[0].size() != dataset_key.size())
1444 return false;
1445 // this loop only tests one
1446 // member for each name,
1447 // i.e. checks the user it will
1448 // not catch internal errors
1449 // which do not update all
1450 // fields for a name.
1451 }
1452 }
1453 return true;
1454 }
1455 if (n_indep != 0)
1456 {
1457 if (std::abs(static_cast<int>(dataset_key.size()) -
1458 static_cast<int>(independent_values[0].size())) >= 2)
1459 {
1460 return false;
1461 }
1462 }
1463
1464 if (have_dof_handler)
1465 {
1466 for (const auto &data_entry : data_store)
1467 {
1468 Assert(data_entry.second.size() > 0, ExcInternalError());
1469
1470 if (std::abs(static_cast<int>((data_entry.second)[0].size()) -
1471 static_cast<int>(dataset_key.size())) >= 2)
1472 return false;
1473 // this loop only tests one member
1474 // for each name, i.e. checks the
1475 // user it will not catch internal
1476 // errors which do not update all
1477 // fields for a name.
1478 }
1479 }
1480 return true;
1481}
1482
1483
1484
1485template <int dim>
1486void
1488{
1489 // this function is called by the
1490 // Triangulation whenever something
1491 // changes, by virtue of having
1492 // attached the function to the
1493 // signal handler in the
1494 // triangulation object
1495
1496 // we record the fact that the mesh
1497 // has changed. we need to take
1498 // this into account next time we
1499 // evaluate the solution
1500 triangulation_changed = true;
1501}
1502
1503
1504// explicit instantiations
1505#include "point_value_history.inst"
1506
1507
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)
void evaluate_field(const std::string &name, const VectorType &solution)
SmartPointer< const DoFHandler< dim >, PointValueHistory< dim > > dof_handler
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:111
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:503
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:504
Point< 2 > second
Definition grid_out.cc:4624
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.
const Mapping< dim, spacedim > & get_default_linear_mapping(const Triangulation< dim, spacedim > &triangulation)
Definition mapping.cc:294
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:470
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