13 #include <unordered_set>
19 namespace ComponentTransport
22 std::vector<std::size_t>
const& vec1, std::vector<std::size_t>
const& vec2)
24 std::unordered_set<std::size_t>
set(vec1.begin(), vec1.end());
25 std::vector<std::size_t> vec;
26 for (
auto const a : vec2)
42 WARN(
"The interpolation point is below the lower bound.");
43 auto const nx = std::next(it);
44 return std::make_pair(*it, *nx);
48 WARN(
"The interpolation point is above the upper bound.");
52 auto const pv = std::prev(it);
53 return std::make_pair(*pv, *it);
57 std::vector<GlobalVector*>
const& x_prev,
58 std::size_t
const n_nodes)
const
60 using EntryInput = std::vector<double>;
62 for (std::size_t node_id = 0; node_id < n_nodes; ++node_id)
64 std::vector<InterpolationPoint> interpolation_points;
65 EntryInput base_entry_input;
70 auto const process_id = input_field.variable_id;
71 auto const& variable_name = input_field.name;
72 double input_field_value =
73 variable_name.find(
"_prev") == std::string::npos
74 ? x[process_id]->get(node_id)
75 : x_prev[process_id]->get(node_id);
77 (std::abs(input_field_value) + input_field_value) / 2;
78 auto bounding_seed_points =
79 input_field.getBoundingSeedPoints(input_field_value);
83 interpolation_points.push_back(interpolation_point);
85 base_entry_input.push_back(bounding_seed_points.first);
91 EntryInput bounding_entry_input{base_entry_input};
92 std::vector<std::size_t> bounding_entry_ids;
93 for (std::size_t i = 0; i < interpolation_points.size(); ++i)
95 bounding_entry_input[i] =
96 interpolation_points[i].bounding_points.second;
98 bounding_entry_input[i] =
99 interpolation_points[i].bounding_points.first;
104 if (input_field.name.find(
"_prev") != std::string::npos)
109 auto const output_field_name = input_field.name +
"_new";
110 auto const base_value =
112 auto new_value = base_value;
115 for (std::size_t i = 0; i < interpolation_points.size(); ++i)
117 auto const interpolation_point_value =
118 tabular_data.at(output_field_name)[bounding_entry_ids[i]];
120 (interpolation_point_value - base_value) /
121 (interpolation_points[i].bounding_points.second -
122 interpolation_points[i].bounding_points.first);
125 slope * (interpolation_points[i].value -
126 interpolation_points[i].bounding_points.first);
129 x[input_field.variable_id]->set(node_id, new_value);
135 std::vector<double>
const& entry_input)
const
137 std::vector<std::size_t> temp_vec;
139 std::vector<std::size_t> intersected_vec =
150 std::vector<std::size_t>
const vec =
156 std::swap(intersected_vec, temp_vec);
160 return intersected_vec[0];
void WARN(char const *fmt, Args const &... args)
std::size_t findIndex(Container const &container, typename Container::value_type const &element)
void set(PETScVector &x, double const a)
static std::vector< std::size_t > intersection(std::vector< std::size_t > const &vec1, std::vector< std::size_t > const &vec2)
std::pair< double, double > getBoundingSeedPoints(double const value) const
std::vector< double > const seed_points
void lookup(std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, std::size_t const n_nodes) const
std::size_t getTableEntryID(std::vector< double > const &entry_input) const
std::map< std::string, std::vector< double > > const tabular_data
std::vector< Field > const input_fields