17 #include <deal.II/lac/vector.h> 18 #include <deal.II/lac/block_vector.h> 19 #include <deal.II/lac/la_vector.h> 20 #include <deal.II/lac/la_parallel_vector.h> 21 #include <deal.II/lac/la_parallel_block_vector.h> 22 #include <deal.II/lac/petsc_vector.h> 23 #include <deal.II/lac/petsc_block_vector.h> 24 #include <deal.II/lac/trilinos_vector.h> 25 #include <deal.II/lac/trilinos_block_vector.h> 27 #include <deal.II/numerics/vector_tools.h> 29 #include <deal.II/numerics/point_value_history.h> 34 DEAL_II_NAMESPACE_OPEN
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;
60 n_indep (n_independent_variables)
64 triangulation_changed =
false;
65 have_dof_handler =
false;
68 dataset_key = std::vector <double> ();
72 = std::vector<std::vector <double> > (n_indep, std::vector <double> (0));
73 indep_names = std::vector <std::string> ();
80 const unsigned int n_independent_variables) :
81 dof_handler (&dof_handler),
82 n_indep (n_independent_variables)
94 = std::vector<std::vector <double> > (
n_indep, std::vector <double> (0));
98 std_cxx11::ref(*
this)));
128 std_cxx11::ref(*
this)));
160 std_cxx11::ref(*
this)));
203 support_point_quadrature,
205 unsigned int n_support_points
207 unsigned int n_components
222 std::vector <unsigned int> current_fe_index (n_components, 0);
223 fe_values.reinit (cell);
224 std::vector <Point <dim> > current_points (n_components,
Point <dim > ());
225 for (
unsigned int support_point = 0;
226 support_point < n_support_points; support_point++)
230 unsigned int component
232 current_points [component] = fe_values.quadrature_point (support_point);
233 current_fe_index [component] = support_point;
246 for (; cell != endc; cell++)
248 fe_values.reinit (cell);
250 for (
unsigned int support_point = 0;
251 support_point < n_support_points; support_point++)
253 unsigned int component
256 = fe_values.quadrature_point (support_point);
258 if (location.
distance (test_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>
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 .push_back (local_dof_indices[current_fe_index [component]]);
301 new_point_geometry_data (location, current_points, new_solution_indices);
304 std::map <std::string, std::vector <std::vector <double> > >::iterator
306 for (; data_store_begin !=
data_store.end (); ++data_store_begin)
312 for (
unsigned int component = 0; component < n_stored; component++)
314 data_store_begin->second.push_back (std::vector<double> (0));
363 std::vector <typename DoFHandler<dim>::active_cell_iterator > current_cell (locations.size (), cell);
365 fe_values.reinit (cell);
366 std::vector <Point <dim> > temp_points (n_components,
Point <dim > ());
367 std::vector <unsigned int> temp_fe_index (n_components, 0);
368 for (
unsigned int support_point = 0; support_point < n_support_points; support_point++)
373 temp_points [component] = fe_values.quadrature_point (support_point);
374 temp_fe_index [component] = support_point;
376 std::vector <std::vector <Point <dim> > > current_points (locations.size (), temp_points);
377 std::vector <std::vector <unsigned int> > current_fe_index (locations.size (), temp_fe_index);
388 for (; cell != endc; cell++)
390 fe_values.reinit (cell);
391 for (
unsigned int support_point = 0; support_point < n_support_points; support_point++)
394 Point<dim> test_point = fe_values.quadrature_point (support_point);
396 for (
unsigned int point = 0; point < locations.size (); point++)
398 if (locations[point].distance (test_point) < locations[point].distance (current_points[point][component]))
401 current_points[point][component] = test_point;
402 current_cell[point] = cell;
403 current_fe_index[point][component] = support_point;
410 for (
unsigned int point = 0; point < locations.size (); point++)
412 current_cell[point]->get_dof_indices (local_dof_indices);
413 std::vector<types::global_dof_index> new_solution_indices;
417 new_solution_indices.push_back (local_dof_indices[current_fe_index[point][component]]);
424 std::map <std::string, std::vector <std::vector <double> > >::iterator
426 for (; data_store_begin !=
data_store.end (); ++data_store_begin)
432 for (
unsigned int component = 0; component < n_stored; component++)
434 data_store_begin->second.push_back (std::vector<double> (0));
467 std::pair <std::string, std::vector <std::string> >
468 empty_names (vector_name, std::vector <std::string> ());
473 std::pair<std::string, std::vector <std::vector <double> > > pair_data;
474 pair_data.first = vector_name;
482 std::vector < std::vector <double> > vector_size (n_datastreams,
483 std::vector <double> (0));
484 pair_data.second = vector_size;
493 std::vector <bool> temp_mask (n_components,
true);
501 const std::vector <std::string> &component_names)
503 typename std::map <std::string, std::vector <std::string> >::iterator names =
component_names_map.find(vector_name);
506 typename std::map <std::string, ComponentMask>::iterator mask =
component_mask.find(vector_name);
508 unsigned int n_stored = mask->second.n_selected_components();
512 names->second = component_names;
557 template <
typename VectorType>
577 typename std::map <std::string, std::vector <std::vector <double> > >::iterator data_store_field =
data_store.find(vector_name);
580 typename std::map <std::string, ComponentMask>::iterator mask =
component_mask.find(vector_name);
585 typename std::vector <internal::PointValueHistory::PointGeometryData <dim> >::iterator point =
point_geometry_data.begin ();
586 for (
unsigned int data_store_index = 0; point !=
point_geometry_data.end (); ++point, ++data_store_index)
595 if (mask->second[comp])
597 unsigned int solution_index = point->solution_indices[comp];
598 data_store_field->second[data_store_index * n_stored + store_index].push_back (solution (solution_index));
610 template <
typename VectorType>
613 const VectorType &solution,
630 ExcMessage(
"The update of normal vectors may not be requested for evaluation of " 631 "data on cells via DataPostprocessor."));
634 unsigned int n_quadrature_points = quadrature.
size();
636 unsigned int n_output_variables = data_postprocessor.
get_names().size();
639 typename std::vector <internal::PointValueHistory::PointGeometryData <dim> >::iterator point =
point_geometry_data.begin ();
640 for (
unsigned int data_store_index = 0; point !=
point_geometry_data.end (); ++point, ++data_store_index)
644 Point <dim> requested_location = point->requested_location;
648 fe_values.reinit (cell);
649 std::vector< Vector< double > > computed_quantities (1,
Vector <double> (n_output_variables));
652 if (n_components == 1)
656 std::vector< typename VectorType::value_type > solution_values (n_quadrature_points, 0.0);
659 std::vector<Point<dim> > dummy_normals (1,
Point<dim> ());
660 std::vector<Point<dim> > evaluation_points;
665 fe_values.get_function_values (solution,
668 fe_values.get_function_gradients (solution,
671 fe_values.get_function_hessians (solution,
675 evaluation_points = fe_values.get_quadrature_points();
676 double distance = cell->diameter ();
677 unsigned int selected_point = 0;
678 for (
unsigned int q_point = 0; q_point < n_quadrature_points; q_point++)
680 if (requested_location.
distance (evaluation_points[q_point]) < distance)
682 selected_point = q_point;
683 distance = requested_location.
distance (evaluation_points[q_point]);
691 compute_derived_quantities_scalar(std::vector< double > (1, solution_values[selected_point]),
695 std::vector<
Point<dim> > (1, evaluation_points[selected_point]),
696 computed_quantities);
705 std::vector<Point<dim> > dummy_normals (1,
Point<dim> ());
706 std::vector<Point<dim> > evaluation_points;
712 fe_values.get_function_values (solution,
715 fe_values.get_function_gradients (solution,
718 fe_values.get_function_hessians (solution,
722 evaluation_points = fe_values.get_quadrature_points();
723 double distance = cell->diameter ();
724 unsigned int selected_point = 0;
725 for (
unsigned int q_point = 0; q_point < n_quadrature_points; q_point++)
727 if (requested_location.
distance (evaluation_points[q_point]) < distance)
729 selected_point = q_point;
730 distance = requested_location.
distance (evaluation_points[q_point]);
738 const std::vector< Tensor< 1, dim, typename VectorType::value_type > > &solution_gradients_s = solution_gradients[selected_point];
739 const std::vector< Tensor< 2, dim, typename VectorType::value_type > > &solution_hessians_s = solution_hessians[selected_point];
740 std::vector< Tensor< 1, dim > > tmp_d (solution_gradients_s.size());
741 for (
unsigned int i = 0; i < solution_gradients_s.size(); i++)
742 tmp_d[i] = solution_gradients_s[i];
744 std::vector< Tensor< 2, dim > > tmp_dd (solution_hessians_s.size());
745 for (
unsigned int i = 0; i < solution_hessians_s.size(); i++)
746 tmp_dd[i] = solution_hessians_s[i];
749 for (
unsigned int i = 0; i < solution_values_s.
size(); i++)
750 tmp[i] = solution_values_s[i];
758 std::vector<
Point<dim> > (1, evaluation_points[selected_point]),
759 computed_quantities);
765 typename std::vector<std::string>::const_iterator name = vector_names.begin();
766 for (; name != vector_names.end(); ++name)
768 typename std::map <std::string, std::vector <std::vector <double> > >::iterator data_store_field =
data_store.find(*name);
771 typename std::map <std::string, ComponentMask>::iterator mask =
component_mask.find(*name);
774 unsigned int n_stored = mask->second.n_selected_components(n_output_variables);
778 for (
unsigned int store_index = 0, comp = 0; comp < n_output_variables; comp++)
780 if (mask->second[comp])
782 data_store_field->second[data_store_index * n_stored + store_index].push_back (computed_quantities[0](comp));
792 template <
typename VectorType>
795 const VectorType &solution,
799 std::vector <std::string> vector_names;
800 vector_names.push_back (vector_name);
801 evaluate_field (vector_names, solution, data_postprocessor, quadrature);
807 template <
typename VectorType>
810 const VectorType &solution)
812 typedef typename VectorType::value_type number;
828 typename std::map <std::string, std::vector <std::vector <double> > >::iterator data_store_field =
data_store.find(vector_name);
831 typename std::map <std::string, ComponentMask>::iterator mask =
component_mask.find(vector_name);
836 typename std::vector <internal::PointValueHistory::PointGeometryData <dim> >::iterator point =
point_geometry_data.begin ();
838 for (
unsigned int data_store_index = 0; point !=
point_geometry_data.end (); ++point, ++data_store_index)
847 for (
unsigned int store_index = 0, comp = 0; comp < mask->second.size(); comp++)
849 if (mask->second[comp])
851 data_store_field->second[data_store_index * n_stored + store_index].push_back (value (comp));
886 for (
unsigned int component = 0; component <
n_indep; component++)
895 const std::vector <
Point <dim> > &postprocessor_locations)
904 std::string filename = base_name +
"_indep.gpl";
905 std::ofstream to_gnuplot (filename.c_str ());
907 to_gnuplot <<
"# Data independent of mesh location\n";
910 to_gnuplot <<
"# <Key> ";
914 for (
unsigned int name = 0; name <
indep_names.size(); name++)
922 for (
unsigned int component = 0; component <
n_indep; component++)
924 to_gnuplot <<
"<Indep_" << component <<
"> ";
929 for (
unsigned int key = 0; key <
dataset_key.size (); key++)
933 for (
unsigned int component = 0; component <
n_indep; component++)
964 typename std::vector <internal::PointValueHistory::PointGeometryData <dim> >::iterator point =
point_geometry_data.begin ();
965 for (
unsigned int data_store_index = 0; point !=
point_geometry_data.end (); ++point, ++data_store_index)
974 std::ofstream to_gnuplot (filename.c_str ());
979 to_gnuplot <<
"# Requested location: " << point->requested_location <<
"\n";
980 to_gnuplot <<
"# DoF_index : Support location (for each component)\n";
983 to_gnuplot <<
"# " << point->solution_indices[component] <<
" : " << point->support_point_locations [component] <<
"\n";
986 to_gnuplot <<
"# (Original components and locations, may be invalidated by mesh change.)\n";
988 if (postprocessor_locations.size() != 0)
990 to_gnuplot <<
"# Postprocessor location: " << postprocessor_locations[data_store_index];
992 to_gnuplot <<
" (may be approximate)\n";
998 to_gnuplot <<
"# <Key> ";
1002 for (
unsigned int name = 0; name <
indep_names.size(); name++)
1009 for (
unsigned int component = 0; component <
n_indep; component++)
1011 to_gnuplot <<
"<Indep_" << component <<
"> ";
1015 for (std::map <std::string, std::vector <std::vector <double> > >::iterator
1016 data_store_begin =
data_store.begin (); data_store_begin !=
data_store.end (); ++data_store_begin)
1018 typename std::map <std::string, ComponentMask>::iterator mask =
component_mask.find(data_store_begin->first);
1019 unsigned int n_stored = mask->second.n_selected_components();
1020 std::vector <std::string> names = (
component_names_map.find (data_store_begin->first))->second;
1022 if (names.size() > 0)
1025 for (
unsigned int component = 0; component < names.size(); component++)
1027 to_gnuplot <<
"<" << names[component] <<
"> ";
1032 for (
unsigned int component = 0; component < n_stored; component++)
1034 to_gnuplot <<
"<" << data_store_begin->first <<
"_" << component <<
"> ";
1041 for (
unsigned int key = 0; key <
dataset_key.size (); key++)
1045 for (
unsigned int component = 0; component <
n_indep; component++)
1050 for (std::map <std::string, std::vector <std::vector <double> > >::iterator
1052 data_store_begin !=
data_store.end (); ++data_store_begin)
1054 typename std::map <std::string, ComponentMask>::iterator mask =
component_mask.find(data_store_begin->first);
1055 unsigned int n_stored = mask->second.n_selected_components();
1057 for (
unsigned int component = 0; component < n_stored; component++)
1059 to_gnuplot <<
" " << (data_store_begin->second)[data_store_index * n_stored + component][key];
1065 to_gnuplot.close ();
1084 typename std::vector <internal::PointValueHistory::PointGeometryData <dim> >::iterator point =
point_geometry_data.begin ();
1089 dof_vector (point->solution_indices[component]) = 1;
1104 std::vector <std::vector <Point <dim> > > actual_points;
1105 typename std::vector <internal::PointValueHistory::PointGeometryData <dim> >::iterator point =
point_geometry_data.begin ();
1109 actual_points.push_back (point->support_point_locations);
1111 locations = actual_points;
1130 locations = std::vector<Point <dim> > ();
1133 unsigned int n_quadrature_points = quadrature.
size();
1134 std::vector<Point<dim> > evaluation_points;
1137 typename std::vector <internal::PointValueHistory::PointGeometryData <dim> >::iterator point =
point_geometry_data.begin ();
1138 for (
unsigned int data_store_index = 0; point !=
point_geometry_data.end (); ++point, ++data_store_index)
1142 Point <dim> requested_location = point->requested_location;
1144 fe_values.reinit (cell);
1146 evaluation_points = fe_values.get_quadrature_points();
1147 double distance = cell->diameter ();
1148 unsigned int selected_point = 0;
1150 for (
unsigned int q_point = 0; q_point < n_quadrature_points; q_point++)
1152 if (requested_location.
distance (evaluation_points[q_point]) < distance)
1154 selected_point = q_point;
1155 distance = requested_location.
distance (evaluation_points[q_point]);
1159 locations.push_back (evaluation_points[selected_point]);
1168 out <<
"***PointValueHistory status output***\n\n";
1169 out <<
"Closed: " <<
closed <<
"\n";
1170 out <<
"Cleared: " <<
cleared <<
"\n";
1173 out <<
"Geometric Data" <<
"\n";
1175 typename std::vector <internal::PointValueHistory::PointGeometryData <dim> >::iterator point =
point_geometry_data.begin ();
1178 out <<
"No points stored currently\n";
1186 out <<
"# Requested location: " << point->requested_location <<
"\n";
1187 out <<
"# DoF_index : Support location (for each component)\n";
1190 out << point->solution_indices[component] <<
" : " << point->support_point_locations [component] <<
"\n";
1197 out <<
"#Cannot access DoF_indices once cleared\n";
1208 for (
unsigned int name = 0; name <
indep_names.size(); name++)
1217 out <<
"No independent values stored\n";
1220 std::map <std::string, std::vector <std::vector <double> > >::iterator
1224 out <<
"Mnemonic: data set size (mask size, n true components) : n data sets\n";
1226 for (; data_store_begin !=
data_store.end (); ++data_store_begin)
1229 std::string vector_name = data_store_begin->first;
1230 typename std::map <std::string, ComponentMask>::iterator mask =
component_mask.find(vector_name);
1232 typename std::map <std::string, std::vector <std::string> >::iterator component_names =
component_names_map.find(vector_name);
1235 if (data_store_begin->second.size () != 0)
1237 out << data_store_begin->first <<
": " << data_store_begin->second.size () <<
" (";
1238 out << mask->second.size() <<
", " << mask->second.n_selected_components() <<
") : ";
1239 out << (data_store_begin->second)[0].size () <<
"\n";
1243 out << data_store_begin->first <<
": " << data_store_begin->second.size () <<
" (";
1244 out << mask->second.size() <<
", " << mask->second.n_selected_components() <<
") : ";
1245 out <<
"No points added" <<
"\n";
1248 if (component_names->second.size() > 0)
1250 for (
unsigned int name = 0; name < component_names->second.size(); name++)
1252 out <<
"<" << component_names->second[name] <<
"> ";
1258 out <<
"***end of status output***\n\n";
1278 std::map <std::string, std::vector <std::vector <double> > >::iterator
1282 for (; data_store_begin !=
data_store.end (); ++data_store_begin)
1284 Assert (data_store_begin->second.size() > 0,
1286 if ((data_store_begin->second)[0].size () !=
dataset_key.size ())
1308 std::map <std::string, std::vector <std::vector <double> > >::iterator
1310 for (; data_store_begin !=
data_store.end (); ++data_store_begin)
1312 Assert (data_store_begin->second.size() > 0,
1315 if (std::abs ((
int) (data_store_begin->second)[0].size () - (int)
dataset_key.size ()) >= 2)
1349 #include "point_value_history.inst" 1352 DEAL_II_NAMESPACE_CLOSE
std::map< std::string, std::vector< std::vector< double > > > data_store
const std::vector< Point< dim > > & get_unit_support_points() const
active_cell_iterator begin_active(const unsigned int level=0) const
static::ExceptionBase & ExcDataLostSync()
void get_postprocessor_locations(const Quadrature< dim > &quadrature, std::vector< Point< dim > > &locations)
PointValueHistory(const unsigned int n_independent_variables=0)
unsigned int n_selected_components(const unsigned int overall_number_of_components=numbers::invalid_unsigned_int) const
const Triangulation< dim, spacedim > & get_triangulation() const
static::ExceptionBase & ExcDoFHandlerChanged()
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)
std::vector< double > dataset_key
Transformed quadrature points.
numbers::NumberTraits< Number >::real_type distance(const Point< dim, Number > &p) const
#define AssertThrow(cond, exc)
std::vector< std::string > indep_names
void add_points(const std::vector< Point< dim > > &locations)
std::vector< internal::PointValueHistory::PointGeometryData< dim > > point_geometry_data
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
static::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
void evaluate_field_at_requested_location(const std::string &name, const VectorType &solution)
std::vector< std::vector< double > > independent_values
void evaluate_field(const std::string &name, const VectorType &solution)
void start_new_dataset(const double key)
bool represents_the_all_selected_mask() const
SmartPointer< const DoFHandler< dim >, PointValueHistory< dim > > dof_handler
PointValueHistory & operator=(const PointValueHistory &point_value_history)
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
void write_gnuplot(const std::string &base_name, const std::vector< Point< dim > > &postprocessor_locations=std::vector< Point< dim > >())
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)
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 status(std::ostream &out)
static::ExceptionBase & ExcNotImplemented()
void add_point(const Point< dim > &location)
void get_points(std::vector< std::vector< Point< dim > > > &locations)
bool has_support_points() const
std::map< std::string, std::vector< std::string > > component_names_map
void get_support_locations(std::vector< std::vector< Point< dim > > > &locations)
virtual std::vector< std::string > get_names() const =0
static::ExceptionBase & ExcDoFHandlerRequired()
const FiniteElement< dim, spacedim > & get_fe() const
static::ExceptionBase & ExcInternalError()
virtual UpdateFlags get_needed_update_flags() const =0
void push_back_independent(const std::vector< double > &independent_values)