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