44 const std::vector<types::global_dof_index> &new_sol_indices)
57 const unsigned int n_independent_variables)
58 :
n_indep(n_independent_variables)
70 std::vector<std::vector<double>>(
n_indep, std::vector<double>(0));
79 const unsigned int n_independent_variables)
81 ,
n_indep(n_independent_variables)
93 std::vector<std::vector<double>>(
n_indep, std::vector<double>(0));
128 dof_handler->get_triangulation().signals.any_change.connect(
161 dof_handler->get_triangulation().signals.any_change.connect(
204 support_point_quadrature,
206 unsigned int n_support_points =
207 dof_handler->get_fe().get_unit_support_points().size();
208 unsigned int n_components =
dof_handler->get_fe(0).n_components();
221 std::vector<unsigned int> current_fe_index(n_components,
224 std::vector<Point<dim>> current_points(n_components,
Point<dim>());
225 for (
unsigned int support_point = 0; support_point < n_support_points;
229 unsigned int component =
230 dof_handler->get_fe().system_to_component_index(support_point).first;
232 current_fe_index[component] = support_point;
245 for (; cell != endc; ++cell)
249 for (
unsigned int support_point = 0; support_point < n_support_points;
253 .system_to_component_index(support_point)
259 location.
distance(current_points[component]))
262 current_points[component] = test_point;
264 current_fe_index[component] = support_point;
270 std::vector<types::global_dof_index> local_dof_indices(
272 std::vector<types::global_dof_index> new_solution_indices;
273 current_cell->get_dof_indices(local_dof_indices);
293 for (
unsigned int component = 0;
297 new_solution_indices.push_back(
298 local_dof_indices[current_fe_index[component]]);
302 new_point_geometry_data(location, current_points, new_solution_indices);
311 data_entry.second.resize(data_entry.second.size() + n_stored);
345 support_point_quadrature,
347 unsigned int n_support_points =
348 dof_handler->get_fe().get_unit_support_points().size();
349 unsigned int n_components =
dof_handler->get_fe(0).n_components();
364 std::vector<typename DoFHandler<dim>::active_cell_iterator> current_cell(
365 locations.size(), cell);
368 std::vector<Point<dim>> temp_points(n_components,
Point<dim>());
369 std::vector<unsigned int> temp_fe_index(n_components, 0);
370 for (
unsigned int support_point = 0; support_point < n_support_points;
374 unsigned int component =
375 dof_handler->get_fe().system_to_component_index(support_point).first;
377 temp_fe_index[component] = support_point;
379 std::vector<std::vector<Point<dim>>> current_points(
380 locations.size(), temp_points);
381 std::vector<std::vector<unsigned int>> current_fe_index(locations.size(),
393 for (; cell != endc; ++cell)
396 for (
unsigned int support_point = 0; support_point < n_support_points;
400 .system_to_component_index(support_point)
405 for (
unsigned int point = 0; point < locations.size(); ++point)
407 if (locations[point].distance(test_point) <
408 locations[point].distance(current_points[point][component]))
411 current_points[point][component] = test_point;
412 current_cell[point] = cell;
413 current_fe_index[point][component] = support_point;
419 std::vector<types::global_dof_index> local_dof_indices(
421 for (
unsigned int point = 0; point < locations.size(); ++point)
423 current_cell[point]->get_dof_indices(local_dof_indices);
424 std::vector<types::global_dof_index> new_solution_indices;
426 for (
unsigned int component = 0;
430 new_solution_indices.push_back(
431 local_dof_indices[current_fe_index[point][component]]);
435 new_point_geometry_data(locations[point],
436 current_points[point],
437 new_solution_indices);
447 data_entry.second.resize(data_entry.second.size() + n_stored);
467 if (mask.represents_the_all_selected_mask() ==
false)
471 std::make_pair(vector_name,
478 std::pair<std::string, std::vector<std::string>> empty_names(
479 vector_name, std::vector<std::string>());
484 std::pair<std::string, std::vector<std::vector<double>>> pair_data;
485 pair_data.first = vector_name;
486 const unsigned int n_stored =
487 (mask.represents_the_all_selected_mask() ==
false ?
488 mask.n_selected_components() :
493 std::vector<std::vector<double>> vector_size(n_datastreams,
494 std::vector<double>(0));
495 pair_data.second = vector_size;
503 const unsigned int n_components)
505 ComponentMask temp_mask(std::vector<bool>(n_components,
true));
513 const std::string &vector_name,
514 const std::vector<std::string> &component_names)
516 typename std::map<std::string, std::vector<std::string>>::iterator names =
521 typename std::map<std::string, ComponentMask>::iterator mask =
524 unsigned int n_stored = mask->second.n_selected_components();
526 Assert(component_names.size() == n_stored,
529 names->second = component_names;
536 const std::vector<std::string> &independent_names)
576template <
typename VectorType>
579 const VectorType &solution)
600 typename std::map<std::string, std::vector<std::vector<double>>>::iterator
601 data_store_field =
data_store.find(vector_name);
605 typename std::map<std::string, ComponentMask>::iterator mask =
609 unsigned int n_stored =
610 mask->second.n_selected_components(
dof_handler->get_fe(0).n_components());
612 typename std::vector<
616 ++point, ++data_store_index)
623 for (
unsigned int store_index = 0, comp = 0;
627 if (mask->second[comp])
629 unsigned int solution_index = point->solution_indices[comp];
631 ->second[data_store_index * n_stored + store_index]
644template <
typename VectorType>
647 const std::vector<std::string> &vector_names,
648 const VectorType &solution,
666 const UpdateFlags update_flags =
671 "The update of normal vectors may not be requested for evaluation of "
672 "data on cells via DataPostprocessor."));
674 unsigned int n_components =
dof_handler->get_fe(0).n_components();
675 unsigned int n_quadrature_points = quadrature.
size();
677 unsigned int n_output_variables = data_postprocessor.
get_names().size();
682 std::vector<typename VectorType::value_type> scalar_solution_values(
683 n_quadrature_points);
684 std::vector<Tensor<1, dim, typename VectorType::value_type>>
685 scalar_solution_gradients(n_quadrature_points);
686 std::vector<Tensor<2, dim, typename VectorType::value_type>>
687 scalar_solution_hessians(n_quadrature_points);
689 std::vector<Vector<typename VectorType::value_type>> vector_solution_values(
692 std::vector<std::vector<Tensor<1, dim, typename VectorType::value_type>>>
693 vector_solution_gradients(
698 std::vector<std::vector<Tensor<2, dim, typename VectorType::value_type>>>
699 vector_solution_hessians(
705 typename std::vector<
710 const auto reference_cell =
711 dof_handler->get_triangulation().get_reference_cells()[0];
713 ++point, ++data_store_index)
716 const Point<dim> requested_location = point->requested_location;
726 std::vector<Vector<double>> computed_quantities(
730 std::vector<Point<dim>> 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)
736 if (requested_location.
distance(quadrature_points[q_point]) <
739 selected_point = q_point;
741 requested_location.
distance(quadrature_points[q_point]);
747 if (n_components == 1)
765 std::vector<double>(1, scalar_solution_values[selected_point]);
770 scalar_solution_gradients);
772 std::vector<Tensor<1, dim>>(
773 1, scalar_solution_gradients[selected_point]);
778 scalar_solution_hessians);
780 std::vector<Tensor<2, dim>>(
781 1, scalar_solution_hessians[selected_point]);
786 std::vector<Point<dim>>(1, quadrature_points[selected_point]);
790 computed_quantities);
802 std::copy(vector_solution_values[selected_point].begin(),
803 vector_solution_values[selected_point].end(),
809 vector_solution_gradients);
812 std::copy(vector_solution_gradients[selected_point].begin(),
813 vector_solution_gradients[selected_point].end(),
819 vector_solution_hessians);
822 std::copy(vector_solution_hessians[selected_point].begin(),
823 vector_solution_hessians[selected_point].end(),
828 std::vector<Point<dim>>(1, quadrature_points[selected_point]);
831 computed_quantities);
837 typename std::vector<std::string>::const_iterator name =
838 vector_names.begin();
839 for (; name != vector_names.end(); ++name)
841 typename std::map<std::string,
842 std::vector<std::vector<double>>>::iterator
847 typename std::map<std::string, ComponentMask>::iterator mask =
852 unsigned int n_stored =
853 mask->second.n_selected_components(n_output_variables);
857 for (
unsigned int store_index = 0, comp = 0;
858 comp < n_output_variables;
861 if (mask->second[comp])
864 ->second[data_store_index * n_stored + store_index]
865 .push_back(computed_quantities[0](comp));
876template <
typename VectorType>
879 const std::string &vector_name,
880 const VectorType &solution,
884 std::vector<std::string> vector_names;
885 vector_names.push_back(vector_name);
886 evaluate_field(vector_names, solution, data_postprocessor, quadrature);
892template <
typename VectorType>
895 const std::string &vector_name,
896 const VectorType &solution)
898 using number =
typename VectorType::value_type;
917 typename std::map<std::string, std::vector<std::vector<double>>>::iterator
918 data_store_field =
data_store.find(vector_name);
922 typename std::map<std::string, ComponentMask>::iterator mask =
926 unsigned int n_stored =
927 mask->second.n_selected_components(
dof_handler->get_fe(0).n_components());
929 typename std::vector<
934 ++point, ++data_store_index)
941 point->requested_location,
946 for (
unsigned int store_index = 0, comp = 0; comp < mask->second.size();
949 if (mask->second[comp])
952 ->second[data_store_index * n_stored + store_index]
953 .push_back(
value(comp));
979 const std::vector<double> &indep_values)
992 for (
unsigned int component = 0; component <
n_indep; ++component)
1001 const std::string &base_name,
1002 const std::vector<
Point<dim>> &postprocessor_locations)
1011 std::string filename = base_name +
"_indep.gpl";
1012 std::ofstream to_gnuplot(filename);
1014 to_gnuplot <<
"# Data independent of mesh location\n";
1017 to_gnuplot <<
"# <Key> ";
1023 to_gnuplot <<
"<" << indep_name <<
"> ";
1029 for (
unsigned int component = 0; component <
n_indep; ++component)
1031 to_gnuplot <<
"<Indep_" << component <<
"> ";
1036 for (
unsigned int key = 0; key <
dataset_key.size(); ++key)
1040 for (
unsigned int component = 0; component <
n_indep; ++component)
1057 postprocessor_locations.size() ==
1075 typename std::vector<internal::PointValueHistoryImplementation::
1076 PointGeometryData<dim>>::iterator point =
1078 for (
unsigned int data_store_index = 0;
1080 ++point, ++data_store_index)
1084 std::string filename = base_name +
"_" +
1091 std::ofstream to_gnuplot(filename);
1096 to_gnuplot <<
"# Requested location: " << point->requested_location
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();
1103 to_gnuplot <<
"# " << point->solution_indices[component] <<
" : "
1104 << point->support_point_locations[component] <<
'\n';
1108 <<
"# (Original components and locations, may be invalidated by mesh change.)\n";
1110 if (postprocessor_locations.size() != 0)
1112 to_gnuplot <<
"# Postprocessor location: "
1113 << postprocessor_locations[data_store_index];
1115 to_gnuplot <<
" (may be approximate)\n";
1117 to_gnuplot <<
"#\n";
1121 to_gnuplot <<
"# <Key> ";
1127 to_gnuplot <<
"<" << indep_name <<
"> ";
1132 for (
unsigned int component = 0; component <
n_indep; ++component)
1134 to_gnuplot <<
"<Indep_" << component <<
"> ";
1140 typename std::map<std::string, ComponentMask>::iterator mask =
1142 unsigned int n_stored = mask->second.n_selected_components();
1143 std::vector<std::string> names =
1146 if (names.size() > 0)
1150 for (
const auto &name : names)
1152 to_gnuplot <<
"<" << name <<
"> ";
1157 for (
unsigned int component = 0; component < n_stored;
1160 to_gnuplot <<
"<" << data_entry.first <<
"_" << component
1168 for (
unsigned int key = 0; key <
dataset_key.size(); ++key)
1172 for (
unsigned int component = 0; component <
n_indep; ++component)
1179 typename std::map<std::string, ComponentMask>::iterator mask =
1181 unsigned int n_stored = mask->second.n_selected_components();
1183 for (
unsigned int component = 0; component < n_stored;
1188 << (data_entry.second)[data_store_index * n_stored +
1214 typename std::vector<
1219 for (
unsigned int component = 0;
1220 component <
dof_handler->get_fe(0).n_components();
1223 dof_vector(point->solution_indices[component]) = 1;
1233 std::vector<std::vector<
Point<dim>>> &locations)
1239 std::vector<std::vector<Point<dim>>> actual_points;
1240 typename std::vector<
1246 actual_points.push_back(point->support_point_locations);
1248 locations = actual_points;
1262 locations = std::vector<Point<dim>>();
1267 unsigned int n_quadrature_points = quadrature.
size();
1268 std::vector<Point<dim>> evaluation_points;
1271 typename std::vector<
1276 const auto reference_cell =
1277 dof_handler->get_triangulation().get_reference_cells()[0];
1279 ++point, ++data_store_index)
1283 Point<dim> requested_location = point->requested_location;
1290 fe_values.reinit(cell);
1292 evaluation_points = fe_values.get_quadrature_points();
1293 double distance = cell->diameter();
1294 unsigned int selected_point = 0;
1296 for (
unsigned int q_point = 0; q_point < n_quadrature_points; ++q_point)
1298 if (requested_location.
distance(evaluation_points[q_point]) <
1301 selected_point = q_point;
1303 requested_location.
distance(evaluation_points[q_point]);
1307 locations.push_back(evaluation_points[selected_point]);
1316 out <<
"***PointValueHistory status output***\n\n";
1317 out <<
"Closed: " <<
closed <<
'\n';
1318 out <<
"Cleared: " <<
cleared <<
'\n';
1321 out <<
"Geometric Data" <<
'\n';
1323 typename std::vector<
1328 out <<
"No points stored currently\n";
1336 out <<
"# Requested location: " << point->requested_location
1338 out <<
"# DoF_index : Support location (for each component)\n";
1339 for (
unsigned int component = 0;
1340 component <
dof_handler->get_fe(0).n_components();
1343 out << point->solution_indices[component] <<
" : "
1344 << point->support_point_locations[component] <<
'\n';
1351 out <<
"#Cannot access DoF_indices once cleared\n";
1365 out <<
"<" << indep_name <<
"> ";
1372 out <<
"No independent values stored\n";
1378 <<
"Mnemonic: data set size (mask size, n true components) : n data sets\n";
1383 std::string vector_name = data_entry.first;
1384 typename std::map<std::string, ComponentMask>::iterator mask =
1388 typename std::map<std::string, std::vector<std::string>>::iterator
1393 if (data_entry.second.size() != 0)
1395 out << data_entry.first <<
": " << data_entry.second.size() <<
" (";
1396 out << mask->second.size() <<
", "
1397 << mask->second.n_selected_components() <<
") : ";
1398 out << (data_entry.second)[0].size() <<
'\n';
1402 out << data_entry.first <<
": " << data_entry.second.size() <<
" (";
1403 out << mask->second.size() <<
", "
1404 << mask->second.n_selected_components() <<
") : ";
1405 out <<
"No points added" <<
'\n';
1408 if (component_names->second.size() > 0)
1410 for (
const auto &name : component_names->second)
1412 out <<
"<" << name <<
"> ";
1418 out <<
"***end of status output***\n\n";
1443 if ((data_entry.second)[0].size() !=
dataset_key.size())
1470 if (
std::abs(
static_cast<int>((data_entry.second)[0].size()) -
1505#include "point_value_history.inst"
unsigned int n_selected_components(const unsigned int overall_number_of_components=numbers::invalid_unsigned_int) const
virtual UpdateFlags get_needed_update_flags() const =0
virtual void evaluate_vector_field(const DataPostprocessorInputs::Vector< dim > &input_data, std::vector< Vector< double > > &computed_quantities) const
virtual void evaluate_scalar_field(const DataPostprocessorInputs::Scalar< dim > &input_data, std::vector< Vector< double > > &computed_quantities) const
virtual std::vector< std::string > get_names() const =0
void get_function_values(const ReadVector< Number > &fe_function, std::vector< Number > &values) const
const std::vector< Point< spacedim > > & get_quadrature_points() const
void get_function_hessians(const ReadVector< Number > &fe_function, std::vector< Tensor< 2, spacedim, Number > > &hessians) const
const Point< spacedim > & quadrature_point(const unsigned int q_point) const
void get_function_gradients(const ReadVector< Number > &fe_function, std::vector< Tensor< 1, spacedim, Number > > &gradients) const
void reinit(const TriaIterator< DoFCellAccessor< dim, spacedim, level_dof_access > > &cell)
void evaluate_field_at_requested_location(const std::string &name, const VectorType &solution)
void add_field_name(const std::string &vector_name, const ComponentMask &component_mask={})
boost::signals2::connection tria_listener
void add_point(const Point< dim > &location)
std::map< std::string, ComponentMask > component_mask
std::vector< internal::PointValueHistoryImplementation::PointGeometryData< dim > > point_geometry_data
std::map< std::string, std::vector< std::string > > component_names_map
PointValueHistory & operator=(const PointValueHistory &point_value_history)
void evaluate_field(const std::string &name, const VectorType &solution)
SmartPointer< const DoFHandler< dim >, PointValueHistory< dim > > dof_handler
Vector< double > mark_support_locations()
std::vector< std::string > indep_names
bool triangulation_changed
std::vector< double > dataset_key
void write_gnuplot(const std::string &base_name, const std::vector< Point< dim > > &postprocessor_locations=std::vector< Point< dim > >())
void status(std::ostream &out)
void add_independent_names(const std::vector< std::string > &independent_names)
PointValueHistory(const unsigned int n_independent_variables=0)
void get_support_locations(std::vector< std::vector< Point< dim > > > &locations)
std::map< std::string, std::vector< std::vector< double > > > data_store
void add_component_names(const std::string &vector_name, const std::vector< std::string > &component_names)
void start_new_dataset(const double key)
void add_points(const std::vector< Point< dim > > &locations)
std::vector< std::vector< double > > independent_values
void push_back_independent(const std::vector< double > &independent_values)
void get_postprocessor_locations(const Quadrature< dim > &quadrature, std::vector< Point< dim > > &locations)
void tria_change_listener()
bool deep_check(const bool strict)
numbers::NumberTraits< Number >::real_type distance(const Point< dim, Number > &p) const
unsigned int size() const
Point< dim > requested_location
std::vector< types::global_dof_index > solution_indices
std::vector< Point< dim > > support_point_locations
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
#define DEAL_II_NAMESPACE_CLOSE
static ::ExceptionBase & ExcNoIndependent()
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
static ::ExceptionBase & ExcDataLostSync()
static ::ExceptionBase & ExcDoFHandlerRequired()
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDoFHandlerChanged()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcInvalidState()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
typename ActiveSelector::active_cell_iterator active_cell_iterator
const bool IsBlockVector< VectorType >::value
@ update_hessians
Second derivatives of shape functions.
@ update_values
Shape function values.
@ update_normal_vectors
Normal vectors.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
const Mapping< dim, spacedim > & get_default_linear_mapping(const Triangulation< dim, spacedim > &triangulation)
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
std::vector< Point< spacedim > > evaluation_points
std::vector< double > solution_values
std::vector< Tensor< 2, spacedim > > solution_hessians
std::vector< Tensor< 1, spacedim > > solution_gradients
std::vector< std::vector< Tensor< 2, spacedim > > > solution_hessians
std::vector<::Vector< double > > solution_values
std::vector< std::vector< Tensor< 1, spacedim > > > solution_gradients
static VectorType::value_type get(const VectorType &V, const types::global_dof_index i)