Reference documentation for deal.II version 9.5.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// Copyright (C) 2009 - 2022 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
17
26#include <deal.II/lac/vector.h>
28
31
32#include <algorithm>
33
34
36
37
38namespace internal
39{
40 namespace PointValueHistoryImplementation
41 {
42 template <int dim>
44 const Point<dim> & new_requested_location,
45 const std::vector<Point<dim>> & new_locations,
46 const std::vector<types::global_dof_index> &new_sol_indices)
47 {
48 requested_location = new_requested_location;
49 support_point_locations = new_locations;
50 solution_indices = new_sol_indices;
51 }
52 } // namespace PointValueHistoryImplementation
53} // namespace internal
54
55
56
57template <int dim>
59 const unsigned int n_independent_variables)
60 : n_indep(n_independent_variables)
61{
62 closed = false;
63 cleared = false;
65 have_dof_handler = false;
66
67 // make a vector for keys
68 dataset_key = std::vector<double>(); // initialize the std::vector
69
70 // make a vector of independent values
72 std::vector<std::vector<double>>(n_indep, std::vector<double>(0));
73 indep_names = std::vector<std::string>();
74}
75
76
77
78template <int dim>
80 const DoFHandler<dim> &dof_handler,
81 const unsigned int n_independent_variables)
82 : dof_handler(&dof_handler)
83 , n_indep(n_independent_variables)
84{
85 closed = false;
86 cleared = false;
88 have_dof_handler = true;
89
90 // make a vector to store keys
91 dataset_key = std::vector<double>(); // initialize the std::vector
92
93 // make a vector for the independent values
95 std::vector<std::vector<double>>(n_indep, std::vector<double>(0));
96 indep_names = std::vector<std::string>();
97
98 tria_listener = dof_handler.get_triangulation().signals.any_change.connect(
99 [this]() { this->tria_change_listener(); });
100}
101
102
103
104template <int dim>
106 const PointValueHistory &point_value_history)
107{
108 dataset_key = point_value_history.dataset_key;
109 independent_values = point_value_history.independent_values;
110 indep_names = point_value_history.indep_names;
111 data_store = point_value_history.data_store;
112 component_mask = point_value_history.component_mask;
113 component_names_map = point_value_history.component_names_map;
114 point_geometry_data = point_value_history.point_geometry_data;
115
116 closed = point_value_history.closed;
117 cleared = point_value_history.cleared;
118
119 dof_handler = point_value_history.dof_handler;
120
121 triangulation_changed = point_value_history.triangulation_changed;
122 have_dof_handler = point_value_history.have_dof_handler;
123 n_indep = point_value_history.n_indep;
124
125 // What to do with tria_listener?
126 // Presume subscribe new instance?
127 if (have_dof_handler)
128 {
129 tria_listener =
130 dof_handler->get_triangulation().signals.any_change.connect(
131 [this]() { this->tria_change_listener(); });
132 }
133}
134
135
136
137template <int dim>
140{
141 dataset_key = point_value_history.dataset_key;
142 independent_values = point_value_history.independent_values;
143 indep_names = point_value_history.indep_names;
144 data_store = point_value_history.data_store;
145 component_mask = point_value_history.component_mask;
146 component_names_map = point_value_history.component_names_map;
147 point_geometry_data = point_value_history.point_geometry_data;
148
149 closed = point_value_history.closed;
150 cleared = point_value_history.cleared;
151
152 dof_handler = point_value_history.dof_handler;
153
154 triangulation_changed = point_value_history.triangulation_changed;
155 have_dof_handler = point_value_history.have_dof_handler;
156 n_indep = point_value_history.n_indep;
157
158 // What to do with tria_listener?
159 // Presume subscribe new instance?
160 if (have_dof_handler)
161 {
162 tria_listener =
163 dof_handler->get_triangulation().signals.any_change.connect(
164 [this]() { this->tria_change_listener(); });
165 }
166
167 return *this;
168}
169
170
171
172template <int dim>
174{
175 if (have_dof_handler)
176 {
177 tria_listener.disconnect();
178 }
179}
180
181
182
183template <int dim>
184void
186{
187 // can't be closed to add additional points
188 // or vectors
189 AssertThrow(!closed, ExcInvalidState());
190 AssertThrow(!cleared, ExcInvalidState());
191 AssertThrow(have_dof_handler, ExcDoFHandlerRequired());
192 AssertThrow(!triangulation_changed, ExcDoFHandlerChanged());
193
194 // Implementation assumes that support
195 // points locations are dofs locations
196 AssertThrow(dof_handler->get_fe().has_support_points(), ExcNotImplemented());
197
198 // While in general quadrature points seems
199 // to refer to Gauss quadrature points, in
200 // this case the quadrature points are
201 // forced to be the support points of the
202 // FE.
203 Quadrature<dim> support_point_quadrature(
204 dof_handler->get_fe().get_unit_support_points());
205 FEValues<dim> fe_values(dof_handler->get_fe(),
206 support_point_quadrature,
208 unsigned int n_support_points =
209 dof_handler->get_fe().get_unit_support_points().size();
210 unsigned int n_components = dof_handler->get_fe(0).n_components();
211
212 // set up a loop over all the cells in the
213 // DoFHandler
215 dof_handler->begin_active();
216 typename DoFHandler<dim>::active_cell_iterator endc = dof_handler->end();
217
218 // default values to be replaced as closer
219 // points are found however they need to be
220 // consistent in case they are actually
221 // chosen
222 typename DoFHandler<dim>::active_cell_iterator current_cell = cell;
223 std::vector<unsigned int> current_fe_index(n_components,
224 0); // need one index per component
225 fe_values.reinit(cell);
226 std::vector<Point<dim>> current_points(n_components, Point<dim>());
227 for (unsigned int support_point = 0; support_point < n_support_points;
228 support_point++)
229 {
230 // set up valid data in the empty vectors
231 unsigned int component =
232 dof_handler->get_fe().system_to_component_index(support_point).first;
233 current_points[component] = fe_values.quadrature_point(support_point);
234 current_fe_index[component] = support_point;
235 }
236
237 // check each cell to find a suitable
238 // support points
239 // GridTools::find_active_cell_around_point
240 // is an alternative. That method is not
241 // used here mostly because of the history
242 // of the class. The algorithm used in
243 // add_points below may be slightly more
244 // efficient than find_active_cell_around_point
245 // because it operates on a set of points.
246
247 for (; cell != endc; ++cell)
248 {
249 fe_values.reinit(cell);
250
251 for (unsigned int support_point = 0; support_point < n_support_points;
252 support_point++)
253 {
254 unsigned int component = dof_handler->get_fe()
255 .system_to_component_index(support_point)
256 .first;
257 const Point<dim> &test_point =
258 fe_values.quadrature_point(support_point);
259
260 if (location.distance(test_point) <
261 location.distance(current_points[component]))
262 {
263 // save the data
264 current_points[component] = test_point;
265 current_cell = cell;
266 current_fe_index[component] = support_point;
267 }
268 }
269 }
270
271
272 std::vector<types::global_dof_index> local_dof_indices(
273 dof_handler->get_fe().n_dofs_per_cell());
274 std::vector<types::global_dof_index> new_solution_indices;
275 current_cell->get_dof_indices(local_dof_indices);
276 // there is an implicit assumption here
277 // that all the closest support point to
278 // the requested point for all finite
279 // element components lie in the same cell.
280 // this could possibly be violated if
281 // components use different FE orders,
282 // requested points are on the edge or
283 // vertex of a cell and we are unlucky with
284 // floating point rounding. Worst case
285 // scenario however is that the point
286 // selected isn't the closest possible, it
287 // will still lie within one cell distance.
288 // calling
289 // GridTools::find_active_cell_around_point
290 // to obtain a cell to search is an
291 // option for these methods, but currently
292 // the GridTools function does not cater for
293 // a vector of points, and does not seem to
294 // be intrinsicly faster than this method.
295 for (unsigned int component = 0;
296 component < dof_handler->get_fe(0).n_components();
297 component++)
298 {
299 new_solution_indices.push_back(
300 local_dof_indices[current_fe_index[component]]);
301 }
302
304 new_point_geometry_data(location, current_points, new_solution_indices);
305 point_geometry_data.push_back(new_point_geometry_data);
306
307 for (auto &data_entry : data_store)
308 {
309 // add an extra row to each vector entry
310 const ComponentMask &current_mask =
311 (component_mask.find(data_entry.first))->second;
312 unsigned int n_stored = current_mask.n_selected_components();
313 data_entry.second.resize(data_entry.second.size() + n_stored);
314 }
315}
316
317
318
319template <int dim>
320void
321PointValueHistory<dim>::add_points(const std::vector<Point<dim>> &locations)
322{
323 // This algorithm adds points in the same
324 // order as they appear in the vector
325 // locations and users may depend on this
326 // so do not change order added!
327
328 // can't be closed to add additional points or vectors
329 AssertThrow(!closed, ExcInvalidState());
330 AssertThrow(!cleared, ExcInvalidState());
331 AssertThrow(have_dof_handler, ExcDoFHandlerRequired());
332 AssertThrow(!triangulation_changed, ExcDoFHandlerChanged());
333
334
335 // Implementation assumes that support
336 // points locations are dofs locations
337 AssertThrow(dof_handler->get_fe().has_support_points(), ExcNotImplemented());
338
339 // While in general quadrature points seems
340 // to refer to Gauss quadrature points, in
341 // this case the quadrature points are
342 // forced to be the support points of the
343 // FE.
344 Quadrature<dim> support_point_quadrature(
345 dof_handler->get_fe().get_unit_support_points());
346 FEValues<dim> fe_values(dof_handler->get_fe(),
347 support_point_quadrature,
349 unsigned int n_support_points =
350 dof_handler->get_fe().get_unit_support_points().size();
351 unsigned int n_components = dof_handler->get_fe(0).n_components();
352
353 // set up a loop over all the cells in the
354 // DoFHandler
356 dof_handler->begin_active();
357 typename DoFHandler<dim>::active_cell_iterator endc = dof_handler->end();
358
359 // default values to be replaced as closer
360 // points are found however they need to be
361 // consistent in case they are actually
362 // chosen vector <vector>s defined where
363 // previously single vectors were used
364
365 // need to store one value per point per component
366 std::vector<typename DoFHandler<dim>::active_cell_iterator> current_cell(
367 locations.size(), cell);
368
369 fe_values.reinit(cell);
370 std::vector<Point<dim>> temp_points(n_components, Point<dim>());
371 std::vector<unsigned int> temp_fe_index(n_components, 0);
372 for (unsigned int support_point = 0; support_point < n_support_points;
373 support_point++)
374 {
375 // set up valid data in the empty vectors
376 unsigned int component =
377 dof_handler->get_fe().system_to_component_index(support_point).first;
378 temp_points[component] = fe_values.quadrature_point(support_point);
379 temp_fe_index[component] = support_point;
380 }
381 std::vector<std::vector<Point<dim>>> current_points(
382 locations.size(), temp_points); // give a valid start point
383 std::vector<std::vector<unsigned int>> current_fe_index(locations.size(),
384 temp_fe_index);
385
386 // check each cell to find suitable support
387 // points
388 // GridTools::find_active_cell_around_point
389 // is an alternative. That method is not
390 // used here mostly because of the history
391 // of the class. The algorithm used here
392 // may be slightly more
393 // efficient than find_active_cell_around_point
394 // because it operates on a set of points.
395 for (; cell != endc; ++cell)
396 {
397 fe_values.reinit(cell);
398 for (unsigned int support_point = 0; support_point < n_support_points;
399 support_point++)
400 {
401 unsigned int component = dof_handler->get_fe()
402 .system_to_component_index(support_point)
403 .first;
404 const Point<dim> &test_point =
405 fe_values.quadrature_point(support_point);
406
407 for (unsigned int point = 0; point < locations.size(); ++point)
408 {
409 if (locations[point].distance(test_point) <
410 locations[point].distance(current_points[point][component]))
411 {
412 // save the data
413 current_points[point][component] = test_point;
414 current_cell[point] = cell;
415 current_fe_index[point][component] = support_point;
416 }
417 }
418 }
419 }
420
421 std::vector<types::global_dof_index> local_dof_indices(
422 dof_handler->get_fe().n_dofs_per_cell());
423 for (unsigned int point = 0; point < locations.size(); ++point)
424 {
425 current_cell[point]->get_dof_indices(local_dof_indices);
426 std::vector<types::global_dof_index> new_solution_indices;
427
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 std::vector<bool> temp_mask(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.c_str());
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.size() == 0 ||
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.c_str());
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 "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
const std::vector< Point< spacedim > > & get_quadrature_points() const
void get_function_values(const InputVector &fe_function, std::vector< typename InputVector::value_type > &values) const
const Point< spacedim > & quadrature_point(const unsigned int q_point) const
void get_function_gradients(const InputVector &fe_function, std::vector< Tensor< 1, spacedim, typename InputVector::value_type > > &gradients) const
void get_function_hessians(const InputVector &fe_function, std::vector< Tensor< 2, spacedim, typename InputVector::value_type > > &hessians) const
void reinit(const TriaIterator< DoFCellAccessor< dim, spacedim, level_dof_access > > &cell)
void add_field_name(const std::string &vector_name, const ComponentMask &component_mask=ComponentMask())
void evaluate_field_at_requested_location(const std::string &name, const VectorType &solution)
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:112
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:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
Point< 2 > second
Definition grid_out.cc:4616
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::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:471
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