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> 28 #include <deal.II/numerics/point_value_history.h> 29 #include <deal.II/numerics/vector_tools.h> 34 DEAL_II_NAMESPACE_OPEN
39 namespace PointValueHistoryImplementation
46 const std::vector<types::global_dof_index> &new_sol_indices)
48 requested_location = new_requested_location;
49 support_point_locations = new_locations;
50 solution_indices = new_sol_indices;
59 const unsigned int n_independent_variables)
60 : n_indep(n_independent_variables)
72 std::vector<std::vector<double>>(
n_indep, std::vector<double>(0));
81 const unsigned int n_independent_variables)
82 : dof_handler(&dof_handler)
83 ,
n_indep(n_independent_variables)
95 std::vector<std::vector<double>>(
n_indep, std::vector<double>(0));
208 support_point_quadrature,
210 unsigned int n_support_points =
225 std::vector<unsigned int> current_fe_index(n_components,
227 fe_values.reinit(cell);
228 std::vector<Point<dim>> current_points(n_components,
Point<dim>());
229 for (
unsigned int support_point = 0; support_point < n_support_points;
234 unsigned int component =
236 current_points[component] = fe_values.quadrature_point(support_point);
237 current_fe_index[component] = support_point;
250 for (; cell != endc; cell++)
252 fe_values.reinit(cell);
254 for (
unsigned int support_point = 0; support_point < n_support_points;
260 Point<dim> test_point = fe_values.quadrature_point(support_point);
263 location.
distance(current_points[component]))
266 current_points[component] = test_point;
268 current_fe_index[component] = support_point;
274 std::vector<types::global_dof_index> local_dof_indices(
276 std::vector<types::global_dof_index> new_solution_indices;
277 current_cell->get_dof_indices(local_dof_indices);
297 for (
unsigned int component = 0;
301 new_solution_indices.push_back(
302 local_dof_indices[current_fe_index[component]]);
306 new_point_geometry_data(location, current_points, new_solution_indices);
309 std::map<std::string, std::vector<std::vector<double>>>::iterator
311 for (; data_store_begin !=
data_store.end(); ++data_store_begin)
318 data_store_begin->second.resize(data_store_begin->second.size() +
353 support_point_quadrature,
355 unsigned int n_support_points =
372 std::vector<typename DoFHandler<dim>::active_cell_iterator> current_cell(
373 locations.size(), cell);
375 fe_values.reinit(cell);
376 std::vector<Point<dim>> temp_points(n_components,
Point<dim>());
377 std::vector<unsigned int> temp_fe_index(n_components, 0);
378 for (
unsigned int support_point = 0; support_point < n_support_points;
383 unsigned int component =
385 temp_points[component] = fe_values.quadrature_point(support_point);
386 temp_fe_index[component] = support_point;
388 std::vector<std::vector<Point<dim>>> current_points(
389 locations.size(), temp_points);
390 std::vector<std::vector<unsigned int>> current_fe_index(locations.size(),
402 for (; cell != endc; cell++)
404 fe_values.reinit(cell);
405 for (
unsigned int support_point = 0; support_point < n_support_points;
411 Point<dim> test_point = fe_values.quadrature_point(support_point);
413 for (
unsigned int point = 0; point < locations.size(); point++)
415 if (locations[point].distance(test_point) <
416 locations[point].distance(current_points[point][component]))
419 current_points[point][component] = test_point;
420 current_cell[point] = cell;
421 current_fe_index[point][component] = support_point;
427 std::vector<types::global_dof_index> local_dof_indices(
429 for (
unsigned int point = 0; point < locations.size(); point++)
431 current_cell[point]->get_dof_indices(local_dof_indices);
432 std::vector<types::global_dof_index> new_solution_indices;
434 for (
unsigned int component = 0;
438 new_solution_indices.push_back(
439 local_dof_indices[current_fe_index[point][component]]);
443 new_point_geometry_data(locations[point],
444 current_points[point],
445 new_solution_indices);
449 std::map<std::string, std::vector<std::vector<double>>>::iterator
451 for (; data_store_begin !=
data_store.end(); ++data_store_begin)
458 data_store_begin->second.resize(data_store_begin->second.size() +
483 std::make_pair(vector_name,
490 std::pair<std::string, std::vector<std::string>> empty_names(
491 vector_name, std::vector<std::string>());
496 std::pair<std::string, std::vector<std::vector<double>>> pair_data;
497 pair_data.first = vector_name;
498 const unsigned int n_stored =
505 std::vector<std::vector<double>> vector_size(n_datastreams,
506 std::vector<double>(0));
507 pair_data.second = vector_size;
515 const unsigned int n_components)
517 std::vector<bool> temp_mask(n_components,
true);
525 const std::string & vector_name,
526 const std::vector<std::string> &component_names)
528 typename std::map<std::string, std::vector<std::string>>::iterator names =
533 typename std::map<std::string, ComponentMask>::iterator mask =
536 unsigned int n_stored = mask->second.n_selected_components();
538 Assert(component_names.size() == n_stored,
541 names->second = component_names;
548 const std::vector<std::string> &independent_names)
588 template <
typename VectorType>
591 const VectorType & solution)
612 typename std::map<std::string, std::vector<std::vector<double>>>::iterator
613 data_store_field =
data_store.find(vector_name);
617 typename std::map<std::string, ComponentMask>::iterator mask =
621 unsigned int n_stored =
624 typename std::vector<
628 ++point, ++data_store_index)
635 for (
unsigned int store_index = 0, comp = 0;
639 if (mask->second[comp])
641 unsigned int solution_index = point->solution_indices[comp];
643 ->second[data_store_index * n_stored + store_index]
645 internal::ElementAccess<VectorType>::get(solution,
656 template <
typename VectorType>
659 const std::vector<std::string> &vector_names,
660 const VectorType & solution,
683 "The update of normal vectors may not be requested for evaluation of " 684 "data on cells via DataPostprocessor."));
687 unsigned int n_quadrature_points = quadrature.
size();
689 unsigned int n_output_variables = data_postprocessor.
get_names().size();
694 std::vector<typename VectorType::value_type> scalar_solution_values(
695 n_quadrature_points);
696 std::vector<Tensor<1, dim, typename VectorType::value_type>>
697 scalar_solution_gradients(n_quadrature_points);
698 std::vector<Tensor<2, dim, typename VectorType::value_type>>
699 scalar_solution_hessians(n_quadrature_points);
701 std::vector<Vector<typename VectorType::value_type>> vector_solution_values(
704 std::vector<std::vector<Tensor<1, dim, typename VectorType::value_type>>>
705 vector_solution_gradients(
710 std::vector<std::vector<Tensor<2, dim, typename VectorType::value_type>>>
711 vector_solution_hessians(
717 typename std::vector<
721 ++point, ++data_store_index)
724 const Point<dim> requested_location = point->requested_location;
732 fe_values.reinit(cell);
733 std::vector<Vector<double>> computed_quantities(
737 std::vector<Point<dim>> quadrature_points =
738 fe_values.get_quadrature_points();
739 double distance = cell->diameter();
740 unsigned int selected_point = 0;
741 for (
unsigned int q_point = 0; q_point < n_quadrature_points; q_point++)
743 if (requested_location.
distance(quadrature_points[q_point]) <
746 selected_point = q_point;
748 requested_location.
distance(quadrature_points[q_point]);
754 if (n_components == 1)
770 fe_values.get_function_values(solution, scalar_solution_values);
772 std::vector<double>(1, scalar_solution_values[selected_point]);
776 fe_values.get_function_gradients(solution,
777 scalar_solution_gradients);
779 std::vector<Tensor<1, dim>>(
780 1, scalar_solution_gradients[selected_point]);
784 fe_values.get_function_hessians(solution,
785 scalar_solution_hessians);
787 std::vector<Tensor<2, dim>>(
788 1, scalar_solution_hessians[selected_point]);
793 std::vector<Point<dim>>(1, quadrature_points[selected_point]);
797 computed_quantities);
806 fe_values.get_function_values(solution, vector_solution_values);
809 std::copy(vector_solution_values[selected_point].begin(),
810 vector_solution_values[selected_point].end(),
815 fe_values.get_function_gradients(solution,
816 vector_solution_gradients);
819 std::copy(vector_solution_gradients[selected_point].begin(),
820 vector_solution_gradients[selected_point].end(),
825 fe_values.get_function_hessians(solution,
826 vector_solution_hessians);
829 std::copy(vector_solution_hessians[selected_point].begin(),
830 vector_solution_hessians[selected_point].end(),
835 std::vector<Point<dim>>(1, quadrature_points[selected_point]);
838 computed_quantities);
844 typename std::vector<std::string>::const_iterator name =
845 vector_names.begin();
846 for (; name != vector_names.end(); ++name)
848 typename std::map<std::string,
849 std::vector<std::vector<double>>>::iterator
854 typename std::map<std::string, ComponentMask>::iterator mask =
859 unsigned int n_stored =
860 mask->second.n_selected_components(n_output_variables);
864 for (
unsigned int store_index = 0, comp = 0;
865 comp < n_output_variables;
868 if (mask->second[comp])
871 ->second[data_store_index * n_stored + store_index]
872 .push_back(computed_quantities[0](comp));
883 template <
typename VectorType>
886 const std::string & vector_name,
887 const VectorType & solution,
891 std::vector<std::string> vector_names;
892 vector_names.push_back(vector_name);
893 evaluate_field(vector_names, solution, data_postprocessor, quadrature);
899 template <
typename VectorType>
902 const std::string &vector_name,
903 const VectorType & solution)
905 using number =
typename VectorType::value_type;
924 typename std::map<std::string, std::vector<std::vector<double>>>::iterator
925 data_store_field =
data_store.find(vector_name);
929 typename std::map<std::string, ComponentMask>::iterator mask =
933 unsigned int n_stored =
936 typename std::vector<
941 ++point, ++data_store_index)
948 point->requested_location,
953 for (
unsigned int store_index = 0, comp = 0; comp < mask->second.size();
956 if (mask->second[comp])
959 ->second[data_store_index * n_stored + store_index]
960 .push_back(value(comp));
986 const std::vector<double> &indep_values)
999 for (
unsigned int component = 0; component <
n_indep; component++)
1008 const std::string & base_name,
1009 const std::vector<
Point<dim>> &postprocessor_locations)
1018 std::string filename = base_name +
"_indep.gpl";
1019 std::ofstream to_gnuplot(filename.c_str());
1021 to_gnuplot <<
"# Data independent of mesh location\n";
1024 to_gnuplot <<
"# <Key> ";
1028 for (
unsigned int name = 0; name <
indep_names.size(); name++)
1036 for (
unsigned int component = 0; component <
n_indep; component++)
1038 to_gnuplot <<
"<Indep_" << component <<
"> ";
1043 for (
unsigned int key = 0; key <
dataset_key.size(); key++)
1047 for (
unsigned int component = 0; component <
n_indep; component++)
1063 AssertThrow(postprocessor_locations.size() == 0 ||
1064 postprocessor_locations.size() ==
1085 for (
unsigned int data_store_index = 0;
1087 ++point, ++data_store_index)
1091 std::string filename = base_name +
"_" +
1098 std::ofstream to_gnuplot(filename.c_str());
1103 to_gnuplot <<
"# Requested location: " << point->requested_location
1105 to_gnuplot <<
"# DoF_index : Support location (for each component)\n";
1106 for (
unsigned int component = 0;
1110 to_gnuplot <<
"# " << point->solution_indices[component] <<
" : " 1111 << point->support_point_locations[component] <<
"\n";
1115 <<
"# (Original components and locations, may be invalidated by mesh change.)\n";
1117 if (postprocessor_locations.size() != 0)
1119 to_gnuplot <<
"# Postprocessor location: " 1120 << postprocessor_locations[data_store_index];
1122 to_gnuplot <<
" (may be approximate)\n";
1124 to_gnuplot <<
"#\n";
1128 to_gnuplot <<
"# <Key> ";
1132 for (
unsigned int name = 0; name <
indep_names.size(); name++)
1139 for (
unsigned int component = 0; component <
n_indep; component++)
1141 to_gnuplot <<
"<Indep_" << component <<
"> ";
1145 for (std::map<std::string, std::vector<std::vector<double>>>::iterator
1150 typename std::map<std::string, ComponentMask>::iterator mask =
1152 unsigned int n_stored = mask->second.n_selected_components();
1153 std::vector<std::string> names =
1156 if (names.size() > 0)
1160 for (
unsigned int component = 0; component < names.size();
1163 to_gnuplot <<
"<" << names[component] <<
"> ";
1168 for (
unsigned int component = 0; component < n_stored;
1171 to_gnuplot <<
"<" << data_store_begin->first <<
"_" 1172 << component <<
"> ";
1179 for (
unsigned int key = 0; key <
dataset_key.size(); key++)
1183 for (
unsigned int component = 0; component <
n_indep; component++)
1188 for (std::map<std::string,
1189 std::vector<std::vector<double>>>::iterator
1194 typename std::map<std::string, ComponentMask>::iterator mask =
1196 unsigned int n_stored = mask->second.n_selected_components();
1198 for (
unsigned int component = 0; component < n_stored;
1203 << (data_store_begin
1204 ->second)[data_store_index * n_stored + component]
1230 typename std::vector<
1235 for (
unsigned int component = 0;
1239 dof_vector(point->solution_indices[component]) = 1;
1249 std::vector<std::vector<
Point<dim>>> &locations)
1255 std::vector<std::vector<Point<dim>>> actual_points;
1256 typename std::vector<
1262 actual_points.push_back(point->support_point_locations);
1264 locations = actual_points;
1271 std::vector<std::vector<
Point<dim>>> &locations)
1286 locations = std::vector<Point<dim>>();
1291 unsigned int n_quadrature_points = quadrature.
size();
1292 std::vector<Point<dim>> evaluation_points;
1295 typename std::vector<
1299 ++point, ++data_store_index)
1303 Point<dim> requested_location = point->requested_location;
1309 fe_values.reinit(cell);
1311 evaluation_points = fe_values.get_quadrature_points();
1312 double distance = cell->diameter();
1313 unsigned int selected_point = 0;
1315 for (
unsigned int q_point = 0; q_point < n_quadrature_points; q_point++)
1317 if (requested_location.
distance(evaluation_points[q_point]) <
1320 selected_point = q_point;
1322 requested_location.
distance(evaluation_points[q_point]);
1326 locations.push_back(evaluation_points[selected_point]);
1335 out <<
"***PointValueHistory status output***\n\n";
1336 out <<
"Closed: " <<
closed <<
"\n";
1337 out <<
"Cleared: " <<
cleared <<
"\n";
1340 out <<
"Geometric Data" 1343 typename std::vector<
1348 out <<
"No points stored currently\n";
1356 out <<
"# Requested location: " << point->requested_location
1358 out <<
"# DoF_index : Support location (for each component)\n";
1359 for (
unsigned int component = 0;
1363 out << point->solution_indices[component] <<
" : " 1364 << point->support_point_locations[component] <<
"\n";
1371 out <<
"#Cannot access DoF_indices once cleared\n";
1383 for (
unsigned int name = 0; name <
indep_names.size(); name++)
1392 out <<
"No independent values stored\n";
1395 std::map<std::string, std::vector<std::vector<double>>>::iterator
1400 <<
"Mnemonic: data set size (mask size, n true components) : n data sets\n";
1402 for (; data_store_begin !=
data_store.end(); ++data_store_begin)
1405 std::string vector_name = data_store_begin->first;
1406 typename std::map<std::string, ComponentMask>::iterator mask =
1410 typename std::map<std::string, std::vector<std::string>>::iterator
1415 if (data_store_begin->second.size() != 0)
1417 out << data_store_begin->first <<
": " 1418 << data_store_begin->second.size() <<
" (";
1419 out << mask->second.size() <<
", " 1420 << mask->second.n_selected_components() <<
") : ";
1421 out << (data_store_begin->second)[0].size() <<
"\n";
1425 out << data_store_begin->first <<
": " 1426 << data_store_begin->second.size() <<
" (";
1427 out << mask->second.size() <<
", " 1428 << mask->second.n_selected_components() <<
") : ";
1429 out <<
"No points added" 1433 if (component_names->second.size() > 0)
1435 for (
unsigned int name = 0; name < component_names->second.size();
1438 out <<
"<" << component_names->second[name] <<
"> ";
1444 out <<
"***end of status output***\n\n";
1464 std::map<std::string, std::vector<std::vector<double>>>::iterator
1468 for (; data_store_begin !=
data_store.end(); ++data_store_begin)
1471 if ((data_store_begin->second)[0].size() !=
dataset_key.size())
1494 std::map<std::string, std::vector<std::vector<double>>>::iterator
1496 for (; data_store_begin !=
data_store.end(); ++data_store_begin)
1500 if (std::abs((
int)(data_store_begin->second)[0].size() -
1535 #include "point_value_history.inst" 1538 DEAL_II_NAMESPACE_CLOSE
std::vector< std::vector< double > > independent_values
virtual void evaluate_scalar_field(const DataPostprocessorInputs::Scalar< dim > &input_data, std::vector< Vector< double >> &computed_quantities) const
const std::vector< Point< dim > > & get_unit_support_points() const
active_cell_iterator begin_active(const unsigned int level=0) const
static::ExceptionBase & ExcDataLostSync()
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
unsigned int n_selected_components(const unsigned int overall_number_of_components=numbers::invalid_unsigned_int) const
std::vector< double > solution_values
void get_support_locations(std::vector< std::vector< Point< dim >>> &locations)
static::ExceptionBase & ExcDoFHandlerChanged()
std::vector< double > dataset_key
Transformed quadrature points.
numbers::NumberTraits< Number >::real_type distance(const Point< dim, Number > &p) const
#define AssertThrow(cond, exc)
const Triangulation< dim, spacedim > & get_triangulation() const
std::vector< std::string > indep_names
SmartPointer< const DoFHandler< dim >, PointValueHistory< dim > > dof_handler
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
unsigned int size() const
static::ExceptionBase & ExcInvalidState()
static::ExceptionBase & ExcMessage(std::string arg1)
#define Assert(cond, exc)
types::global_dof_index n_dofs() const
std::vector< std::vector< Tensor< 1, spacedim > > > solution_gradients
static::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
void evaluate_field_at_requested_location(const std::string &name, const VectorType &solution)
PointGeometryData(const Point< dim > &new_requested_location, const std::vector< Point< dim >> &new_locations, const std::vector< types::global_dof_index > &new_sol_indices)
Only a constructor needed for this class (a struct really)
void get_postprocessor_locations(const Quadrature< dim > &quadrature, std::vector< Point< dim >> &locations)
std::map< std::string, std::vector< std::vector< double > > > data_store
void evaluate_field(const std::string &name, const VectorType &solution)
void start_new_dataset(const double key)
bool represents_the_all_selected_mask() const
PointValueHistory & operator=(const PointValueHistory &point_value_history)
std::vector< Tensor< 2, spacedim > > solution_hessians
unsigned int n_components() const
Second derivatives of shape functions.
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
cell_iterator end() const
const unsigned int dofs_per_cell
virtual void evaluate_vector_field(const DataPostprocessorInputs::Vector< dim > &input_data, std::vector< Vector< double >> &computed_quantities) const
std::vector< Tensor< 1, spacedim > > solution_gradients
void add_independent_names(const std::vector< std::string > &independent_names)
void tria_change_listener()
Vector< double > mark_support_locations()
bool deep_check(const bool strict)
std::vector< Point< spacedim > > evaluation_points
const FiniteElement< dim, spacedim > & get_fe(const unsigned int index=0) const
bool triangulation_changed
std::pair< unsigned int, unsigned int > system_to_component_index(const unsigned int index) const
std::map< std::string, ComponentMask > component_mask
Shape function gradients.
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)
void status(std::ostream &out)
static::ExceptionBase & ExcNotImplemented()
void add_point(const Point< dim > &location)
bool has_support_points() const
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)