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));
116 closed = point_value_history.
closed;
117 cleared = point_value_history.
cleared;
123 n_indep = point_value_history.
n_indep;
127 if (have_dof_handler)
130 dof_handler->get_triangulation().signals.any_change.connect(
131 [
this]() { this->tria_change_listener(); });
149 closed = point_value_history.
closed;
150 cleared = point_value_history.
cleared;
156 n_indep = point_value_history.
n_indep;
160 if (have_dof_handler)
163 dof_handler->get_triangulation().signals.any_change.connect(
164 [
this]() { this->tria_change_listener(); });
175 if (have_dof_handler)
177 tria_listener.disconnect();
191 AssertThrow(have_dof_handler, ExcDoFHandlerRequired());
192 AssertThrow(!triangulation_changed, ExcDoFHandlerChanged());
204 dof_handler->get_fe().get_unit_support_points());
206 support_point_quadrature,
208 unsigned int n_support_points =
209 dof_handler->get_fe().get_unit_support_points().size();
210 unsigned int n_components = dof_handler->get_fe(0).n_components();
223 std::vector<unsigned int> current_fe_index(
n_components,
227 for (
unsigned int support_point = 0; support_point < n_support_points;
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;
248 for (; cell != endc; cell++)
252 for (
unsigned int support_point = 0; support_point < n_support_points;
255 unsigned int component = dof_handler->get_fe()
256 .system_to_component_index(support_point)
258 Point<dim> test_point = fe_values.quadrature_point(support_point);
261 location.
distance(current_points[component]))
264 current_points[component] = test_point;
266 current_fe_index[component] = support_point;
272 std::vector<types::global_dof_index> local_dof_indices(
273 dof_handler->get_fe().dofs_per_cell);
274 std::vector<types::global_dof_index> new_solution_indices;
275 current_cell->get_dof_indices(local_dof_indices);
295 for (
unsigned int component = 0;
299 new_solution_indices.push_back(
300 local_dof_indices[current_fe_index[component]]);
304 new_point_geometry_data(location, current_points, new_solution_indices);
305 point_geometry_data.push_back(new_point_geometry_data);
307 for (
auto &data_entry : data_store)
311 (component_mask.find(data_entry.first))->
second;
313 data_entry.second.resize(data_entry.second.size() + n_stored);
331 AssertThrow(have_dof_handler, ExcDoFHandlerRequired());
332 AssertThrow(!triangulation_changed, ExcDoFHandlerChanged());
345 dof_handler->get_fe().get_unit_support_points());
347 support_point_quadrature,
349 unsigned int n_support_points =
350 dof_handler->get_fe().get_unit_support_points().size();
351 unsigned int n_components = dof_handler->get_fe(0).n_components();
366 std::vector<typename DoFHandler<dim>::active_cell_iterator> current_cell(
367 locations.size(), cell);
371 std::vector<unsigned int> temp_fe_index(
n_components, 0);
372 for (
unsigned int support_point = 0; support_point < n_support_points;
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;
382 std::vector<std::vector<Point<dim>>> current_points(
383 locations.size(), temp_points);
384 std::vector<std::vector<unsigned int>> current_fe_index(locations.size(),
396 for (; cell != endc; cell++)
399 for (
unsigned int support_point = 0; support_point < n_support_points;
402 unsigned int component = dof_handler->get_fe()
403 .system_to_component_index(support_point)
405 Point<dim> test_point = fe_values.quadrature_point(support_point);
409 if (locations[
point].distance(test_point) <
410 locations[
point].distance(current_points[
point][component]))
413 current_points[
point][component] = test_point;
414 current_cell[
point] = cell;
415 current_fe_index[
point][component] = support_point;
421 std::vector<types::global_dof_index> local_dof_indices(
422 dof_handler->get_fe().dofs_per_cell);
425 current_cell[
point]->get_dof_indices(local_dof_indices);
426 std::vector<types::global_dof_index> new_solution_indices;
428 for (
unsigned int component = 0;
432 new_solution_indices.push_back(
433 local_dof_indices[current_fe_index[
point][component]]);
437 new_point_geometry_data(locations[
point],
438 current_points[
point],
439 new_solution_indices);
441 point_geometry_data.push_back(new_point_geometry_data);
443 for (
auto &data_entry : data_store)
447 (component_mask.find(data_entry.first))->
second;
449 data_entry.second.resize(data_entry.second.size() + n_stored);
465 AssertThrow(have_dof_handler, ExcDoFHandlerRequired());
466 AssertThrow(!triangulation_changed, ExcDoFHandlerChanged());
470 component_mask.insert(std::make_pair(vector_name, mask));
472 component_mask.insert(
473 std::make_pair(vector_name,
475 dof_handler->get_fe(0).n_components(),
true))));
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);
486 std::pair<std::string, std::vector<std::vector<double>>> pair_data;
487 pair_data.first = vector_name;
488 const unsigned int n_stored =
491 dof_handler->get_fe(0).n_components());
494 point_geometry_data.size() * n_stored;
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);
508 add_field_name(vector_name, temp_mask);
515 const std::string & vector_name,
516 const std::vector<std::string> &component_names)
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(),
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();
528 Assert(component_names.size() == n_stored,
531 names->second = component_names;
538 const std::vector<std::string> &independent_names)
540 Assert(independent_names.size() == n_indep,
543 indep_names = independent_names;
561 dof_handler =
nullptr;
562 have_dof_handler =
false;
578 template <
typename VectorType>
587 AssertThrow(have_dof_handler, ExcDoFHandlerRequired());
588 AssertThrow(!triangulation_changed, ExcDoFHandlerChanged());
593 Assert(std::abs(
static_cast<int>(dataset_key.size()) -
594 static_cast<int>(independent_values[0].size())) < 2,
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(),
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"));
611 unsigned int n_stored =
612 mask->second.n_selected_components(dof_handler->get_fe(0).n_components());
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)
625 for (
unsigned int store_index = 0, comp = 0;
629 if (mask->second[comp])
631 unsigned int solution_index =
point->solution_indices[comp];
633 ->second[data_store_index * n_stored + store_index]
646 template <
typename VectorType>
649 const std::vector<std::string> &vector_names,
658 AssertThrow(have_dof_handler, ExcDoFHandlerRequired());
662 Assert(std::abs(
static_cast<int>(dataset_key.size()) -
663 static_cast<int>(independent_values[0].size())) < 2,
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();
679 unsigned int n_output_variables = data_postprocessor.
get_names().size();
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);
691 std::vector<Vector<typename VectorType::value_type>> vector_solution_values(
694 std::vector<std::vector<Tensor<1, dim, typename VectorType::value_type>>>
695 vector_solution_gradients(
700 std::vector<std::vector<Tensor<2, dim, typename VectorType::value_type>>>
701 vector_solution_hessians(
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)
723 std::vector<Vector<double>> computed_quantities(
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++)
736 selected_point = q_point;
760 fe_values.get_function_values(solution, scalar_solution_values);
762 std::vector<double>(1, scalar_solution_values[selected_point]);
766 fe_values.get_function_gradients(solution,
767 scalar_solution_gradients);
769 std::vector<Tensor<1, dim>>(
770 1, scalar_solution_gradients[selected_point]);
774 fe_values.get_function_hessians(solution,
775 scalar_solution_hessians);
777 std::vector<Tensor<2, dim>>(
778 1, scalar_solution_hessians[selected_point]);
787 computed_quantities);
796 fe_values.get_function_values(solution, vector_solution_values);
800 vector_solution_values[selected_point].
end(),
805 fe_values.get_function_gradients(solution,
806 vector_solution_gradients);
810 vector_solution_gradients[selected_point].
end(),
815 fe_values.get_function_hessians(solution,
816 vector_solution_hessians);
820 vector_solution_hessians[selected_point].
end(),
828 computed_quantities);
834 typename std::vector<std::string>::const_iterator name =
835 vector_names.begin();
836 for (; name != vector_names.end(); ++name)
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(),
844 typename std::map<std::string, ComponentMask>::iterator mask =
845 component_mask.find(*name);
846 Assert(mask != component_mask.end(),
849 unsigned int n_stored =
850 mask->second.n_selected_components(n_output_variables);
854 for (
unsigned int store_index = 0, comp = 0;
855 comp < n_output_variables;
858 if (mask->second[comp])
861 ->second[data_store_index * n_stored + store_index]
862 .push_back(computed_quantities[0](comp));
873 template <
typename VectorType>
876 const std::string & vector_name,
881 std::vector<std::string> vector_names;
882 vector_names.push_back(vector_name);
883 evaluate_field(vector_names, solution, data_postprocessor, quadrature);
889 template <
typename VectorType>
892 const std::string &vector_name,
895 using number =
typename VectorType::value_type;
900 AssertThrow(have_dof_handler, ExcDoFHandlerRequired());
905 Assert(std::abs(
static_cast<int>(dataset_key.size()) -
906 static_cast<int>(independent_values[0].size())) < 2,
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(),
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"));
923 unsigned int n_stored =
924 mask->second.n_selected_components(dof_handler->get_fe(0).n_components());
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)
938 point->requested_location,
943 for (
unsigned int store_index = 0, comp = 0; comp < mask->second.size();
946 if (mask->second[comp])
949 ->second[data_store_index * n_stored + store_index]
950 .push_back(
value(comp));
966 Assert(deep_check(
false), ExcDataLostSync());
968 dataset_key.push_back(key);
976 const std::vector<double> &indep_values)
982 Assert(indep_values.size() == n_indep,
984 Assert(n_indep != 0, ExcNoIndependent());
985 Assert(std::abs(
static_cast<int>(dataset_key.size()) -
986 static_cast<int>(independent_values[0].size())) < 2,
989 for (
unsigned int component = 0; component < n_indep; component++)
990 independent_values[component].push_back(indep_values[component]);
998 const std::string & base_name,
999 const std::vector<
Point<dim>> &postprocessor_locations)
1008 std::string filename = base_name +
"_indep.gpl";
1009 std::ofstream to_gnuplot(filename.c_str());
1011 to_gnuplot <<
"# Data independent of mesh location\n";
1014 to_gnuplot <<
"# <Key> ";
1016 if (indep_names.size() > 0)
1018 for (
const auto &indep_name : indep_names)
1020 to_gnuplot <<
"<" << indep_name <<
"> ";
1026 for (
unsigned int component = 0; component < n_indep; component++)
1028 to_gnuplot <<
"<Indep_" << component <<
"> ";
1033 for (
unsigned int key = 0; key < dataset_key.size(); key++)
1035 to_gnuplot << dataset_key[key];
1037 for (
unsigned int component = 0; component < n_indep; component++)
1039 to_gnuplot <<
" " << independent_values[component][key];
1050 if (have_dof_handler)
1052 AssertThrow(have_dof_handler, ExcDoFHandlerRequired());
1053 AssertThrow(postprocessor_locations.size() == 0 ||
1054 postprocessor_locations.size() ==
1055 point_geometry_data.size(),
1057 point_geometry_data.size()));
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)
1081 std::string filename = base_name +
"_" +
1088 std::ofstream to_gnuplot(filename.c_str());
1093 to_gnuplot <<
"# Requested location: " <<
point->requested_location
1095 to_gnuplot <<
"# DoF_index : Support location (for each component)\n";
1096 for (
unsigned int component = 0;
1100 to_gnuplot <<
"# " <<
point->solution_indices[component] <<
" : "
1101 <<
point->support_point_locations[component] <<
"\n";
1103 if (triangulation_changed)
1105 <<
"# (Original components and locations, may be invalidated by mesh change.)\n";
1107 if (postprocessor_locations.size() != 0)
1109 to_gnuplot <<
"# Postprocessor location: "
1110 << postprocessor_locations[data_store_index];
1111 if (triangulation_changed)
1112 to_gnuplot <<
" (may be approximate)\n";
1114 to_gnuplot <<
"#\n";
1118 to_gnuplot <<
"# <Key> ";
1120 if (indep_names.size() > 0)
1122 for (
const auto &indep_name : indep_names)
1124 to_gnuplot <<
"<" << indep_name <<
"> ";
1129 for (
unsigned int component = 0; component < n_indep; component++)
1131 to_gnuplot <<
"<Indep_" << component <<
"> ";
1135 for (
const auto &data_entry : data_store)
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;
1143 if (names.size() > 0)
1147 for (
const auto &name : names)
1149 to_gnuplot <<
"<" << name <<
"> ";
1154 for (
unsigned int component = 0; component < n_stored;
1157 to_gnuplot <<
"<" << data_entry.first <<
"_" << component
1165 for (
unsigned int key = 0; key < dataset_key.size(); key++)
1167 to_gnuplot << dataset_key[key];
1169 for (
unsigned int component = 0; component < n_indep; component++)
1171 to_gnuplot <<
" " << independent_values[component][key];
1174 for (
const auto &data_entry : data_store)
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();
1180 for (
unsigned int component = 0; component < n_stored;
1185 << (data_entry.second)[data_store_index * n_stored +
1206 AssertThrow(have_dof_handler, ExcDoFHandlerRequired());
1207 AssertThrow(!triangulation_changed, ExcDoFHandlerChanged());
1211 typename std::vector<
1213 point = point_geometry_data.begin();
1214 for (;
point != point_geometry_data.end(); ++
point)
1216 for (
unsigned int component = 0;
1220 dof_vector(
point->solution_indices[component]) = 1;
1230 std::vector<std::vector<
Point<dim>>> &locations)
1233 AssertThrow(have_dof_handler, ExcDoFHandlerRequired());
1234 AssertThrow(!triangulation_changed, ExcDoFHandlerChanged());
1236 std::vector<std::vector<Point<dim>>> actual_points;
1237 typename std::vector<
1239 point = point_geometry_data.begin();
1241 for (;
point != point_geometry_data.end(); ++
point)
1243 actual_points.push_back(
point->support_point_locations);
1245 locations = actual_points;
1252 std::vector<std::vector<
Point<dim>>> &locations)
1254 get_support_locations(locations);
1265 AssertThrow(have_dof_handler, ExcDoFHandlerRequired());
1267 locations = std::vector<Point<dim>>();
1272 unsigned int n_quadrature_points = quadrature.
size();
1273 std::vector<Point<dim>> evaluation_points;
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)
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";
1319 out <<
"Triangulation_changed: " << triangulation_changed <<
"\n";
1320 out <<
"Have_dof_handler: " << have_dof_handler <<
"\n";
1321 out <<
"Geometric Data"
1324 typename std::vector<
1326 point = point_geometry_data.begin();
1327 if (
point == point_geometry_data.end())
1329 out <<
"No points stored currently\n";
1335 for (;
point != point_geometry_data.end(); ++
point)
1337 out <<
"# Requested location: " <<
point->requested_location
1339 out <<
"# DoF_index : Support location (for each component)\n";
1340 for (
unsigned int component = 0;
1344 out <<
point->solution_indices[component] <<
" : "
1345 <<
point->support_point_locations[component] <<
"\n";
1352 out <<
"#Cannot access DoF_indices once cleared\n";
1357 if (independent_values.size() != 0)
1359 out <<
"Independent value(s): " << independent_values.size() <<
" : "
1360 << independent_values[0].size() <<
"\n";
1361 if (indep_names.size() > 0)
1364 for (
const auto &indep_name : indep_names)
1366 out <<
"<" << indep_name <<
"> ";
1373 out <<
"No independent values stored\n";
1376 if (data_store.begin() != data_store.end())
1379 <<
"Mnemonic: data set size (mask size, n true components) : n data sets\n";
1381 for (
const auto &data_entry : data_store)
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(),
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(),
1394 if (data_entry.second.size() != 0)
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";
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"
1410 if (component_names->second.size() > 0)
1412 for (
const auto &name : component_names->second)
1414 out <<
"<" << name <<
"> ";
1420 out <<
"***end of status output***\n\n";
1435 if (dataset_key.size() != independent_values[0].size())
1440 if (have_dof_handler)
1442 for (
const auto &data_entry : data_store)
1445 if ((data_entry.second)[0].size() != dataset_key.size())
1459 if (std::abs(
static_cast<int>(dataset_key.size()) -
1460 static_cast<int>(independent_values[0].size())) >= 2)
1466 if (have_dof_handler)
1468 for (
const auto &data_entry : data_store)
1472 if (std::abs(
static_cast<int>((data_entry.second)[0].size()) -
1473 static_cast<int>(dataset_key.size())) >= 2)
1502 triangulation_changed =
true;
1507 #include "point_value_history.inst"