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