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