Line data Source code
1 : // rbOOmit: An implementation of the Certified Reduced Basis method.
2 : // Copyright (C) 2009, 2010 David J. Knezevic
3 :
4 : // This file is part of rbOOmit.
5 :
6 : // rbOOmit is free software; you can redistribute it and/or
7 : // modify it under the terms of the GNU Lesser General Public
8 : // License as published by the Free Software Foundation; either
9 : // version 2.1 of the License, or (at your option) any later version.
10 :
11 : // rbOOmit is distributed in the hope that it will be useful,
12 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 : // Lesser General Public License for more details.
15 :
16 : // You should have received a copy of the GNU Lesser General Public
17 : // License along with this library; if not, write to the Free Software
18 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 :
20 : // libmesh includes
21 : #include "libmesh/rb_parametrized_function.h"
22 : #include "libmesh/int_range.h"
23 : #include "libmesh/point.h"
24 : #include "libmesh/libmesh_logging.h"
25 : #include "libmesh/utility.h"
26 : #include "libmesh/rb_parameters.h"
27 : #include "libmesh/system.h"
28 : #include "libmesh/elem.h"
29 : #include "libmesh/fem_context.h"
30 : #include "libmesh/quadrature.h"
31 : #include "timpi/parallel_implementation.h"
32 :
33 : namespace libMesh
34 : {
35 :
36 0 : void VectorizedEvalInput::clear()
37 : {
38 0 : all_xyz.clear();
39 0 : elem_ids.clear();
40 0 : qps.clear();
41 0 : sbd_ids.clear();
42 0 : all_xyz_perturb.clear();
43 0 : phi_i_qp.clear();
44 0 : side_indices.clear();
45 0 : boundary_ids.clear();
46 0 : node_ids.clear();
47 0 : elem_types.clear();
48 :
49 0 : elem_id_to_local_index.clear();
50 0 : JxW_all_qp.clear();
51 0 : phi_i_all_qp.clear();
52 0 : dxyzdxi_elem_center.clear();
53 0 : dxyzdeta_elem_center.clear();
54 0 : qrule_orders.clear();
55 :
56 0 : rb_property_map.clear();
57 0 : }
58 :
59 0 : RBParametrizedFunction::RBParametrizedFunction()
60 : :
61 0 : requires_xyz_perturbations(false),
62 0 : requires_all_elem_qp_data(false),
63 0 : requires_all_elem_center_data(false),
64 0 : is_lookup_table(false),
65 0 : fd_delta(1.e-6),
66 0 : _is_nodal_boundary(false)
67 0 : {}
68 :
69 0 : RBParametrizedFunction::~RBParametrizedFunction() = default;
70 :
71 : Number
72 0 : RBParametrizedFunction::evaluate_comp(const RBParameters & mu,
73 : unsigned int comp,
74 : const Point & xyz,
75 : dof_id_type elem_id,
76 : unsigned int qp,
77 : subdomain_id_type subdomain_id,
78 : const std::vector<Point> & xyz_perturb,
79 : const std::vector<Real> & phi_i_qp)
80 : {
81 0 : std::vector<Number> values = evaluate(mu, xyz, elem_id, qp, subdomain_id, xyz_perturb, phi_i_qp);
82 :
83 0 : libmesh_error_msg_if(comp >= values.size(), "Error: Invalid value of comp");
84 :
85 0 : return values[comp];
86 : }
87 :
88 : Number
89 0 : RBParametrizedFunction::side_evaluate_comp(const RBParameters & mu,
90 : unsigned int comp,
91 : const Point & xyz,
92 : dof_id_type elem_id,
93 : unsigned int side_index,
94 : unsigned int qp,
95 : subdomain_id_type subdomain_id,
96 : boundary_id_type boundary_id,
97 : const std::vector<Point> & xyz_perturb,
98 : const std::vector<Real> & phi_i_qp)
99 : {
100 0 : std::vector<Number> values = side_evaluate(mu, xyz, elem_id, side_index, qp, subdomain_id, boundary_id, xyz_perturb, phi_i_qp);
101 :
102 0 : libmesh_error_msg_if(comp >= values.size(), "Error: Invalid value of comp");
103 :
104 0 : return values[comp];
105 : }
106 :
107 : Number
108 0 : RBParametrizedFunction::node_evaluate_comp(const RBParameters & mu,
109 : unsigned int comp,
110 : const Point & xyz,
111 : dof_id_type node_id,
112 : boundary_id_type boundary_id)
113 : {
114 0 : std::vector<Number> values = node_evaluate(mu, xyz, node_id, boundary_id);
115 :
116 0 : libmesh_error_msg_if(comp >= values.size(), "Error: Invalid value of comp");
117 :
118 0 : return values[comp];
119 : }
120 :
121 : std::vector<Number>
122 0 : RBParametrizedFunction::evaluate(const RBParameters & /*mu*/,
123 : const Point & /*xyz*/,
124 : dof_id_type /*elem_id*/,
125 : unsigned int /*qp*/,
126 : subdomain_id_type /*subdomain_id*/,
127 : const std::vector<Point> & /*xyz_perturb*/,
128 : const std::vector<Real> & /*phi_i_qp*/)
129 : {
130 : // This method should be overridden in subclasses, so we just give a not implemented error message here
131 0 : libmesh_not_implemented();
132 :
133 : return std::vector<Number>();
134 : }
135 :
136 : std::vector<Number>
137 0 : RBParametrizedFunction::side_evaluate(const RBParameters & /*mu*/,
138 : const Point & /*xyz*/,
139 : dof_id_type /*elem_id*/,
140 : unsigned int /*side_index*/,
141 : unsigned int /*qp*/,
142 : subdomain_id_type /*subdomain_id*/,
143 : boundary_id_type /*boundary_id*/,
144 : const std::vector<Point> & /*xyz_perturb*/,
145 : const std::vector<Real> & /*phi_i_qp*/)
146 : {
147 : // This method should be overridden in subclasses, so we just give a not implemented error message here
148 0 : libmesh_not_implemented();
149 :
150 : return std::vector<Number>();
151 : }
152 :
153 : std::vector<Number>
154 0 : RBParametrizedFunction::node_evaluate(const RBParameters & /*mu*/,
155 : const Point & /*xyz*/,
156 : dof_id_type /*node_id*/,
157 : boundary_id_type /*boundary_id*/)
158 : {
159 : // This method should be overridden in subclasses, so we just give a not implemented error message here
160 0 : libmesh_not_implemented();
161 :
162 : return std::vector<Number>();
163 : }
164 :
165 0 : void RBParametrizedFunction::vectorized_evaluate(const std::vector<RBParameters> & mus,
166 : const VectorizedEvalInput & v,
167 : std::vector<std::vector<std::vector<Number>>> & output)
168 : {
169 0 : LOG_SCOPE("vectorized_evaluate()", "RBParametrizedFunction");
170 :
171 0 : output.clear();
172 0 : unsigned int n_points = v.all_xyz.size();
173 :
174 0 : libmesh_error_msg_if(v.sbd_ids.size() != n_points, "Error: invalid vector sizes");
175 0 : libmesh_error_msg_if(requires_xyz_perturbations && (v.all_xyz_perturb.size() != n_points), "Error: invalid vector sizes");
176 :
177 : // Dummy vector to be used when xyz perturbations are not required
178 0 : std::vector<Point> empty_perturbs;
179 :
180 : // The number of components returned by this RBParametrizedFunction
181 0 : auto n_components = this->get_n_components();
182 :
183 : // We first loop over all mus and all n_points, calling evaluate()
184 : // for each and storing the results. It is easier to first
185 : // pre-compute all the values before filling output, since, in the
186 : // case of multi-sample RBParameters, the ordering of the loops is a
187 : // bit complicated otherwise.
188 0 : std::vector<std::vector<std::vector<Number>>> all_evals(mus.size());
189 0 : for (auto mu_index : index_range(mus))
190 : {
191 : // Allocate enough space to store all points for the current mu
192 0 : all_evals[mu_index].resize(n_points);
193 0 : for (auto point_index : index_range(all_evals[mu_index]))
194 : {
195 : // Evaluate all samples for the current mu at the current interpolation point
196 0 : all_evals[mu_index][point_index] =
197 0 : this->evaluate(mus[mu_index],
198 0 : v.all_xyz[point_index],
199 0 : v.elem_ids[point_index],
200 0 : v.qps[point_index],
201 0 : v.sbd_ids[point_index],
202 0 : requires_xyz_perturbations ? v.all_xyz_perturb[point_index] : empty_perturbs,
203 0 : v.phi_i_qp[point_index]);
204 :
205 : // The vector returned by evaluate() should contain:
206 : // n_components * mus[mu_index].n_samples()
207 : // entries. That is, for multi-sample RBParameters objects,
208 : // the vector will be packed with entries as follows:
209 : // [sample0_component0, sample0_component1, ..., sample0_componentN,
210 : // sample1_component0, sample1_component1, ..., sample1_componentN,
211 : // ...
212 : // sampleM_component0, sampleM_component1, ..., sampleM_componentN]
213 0 : auto n_samples = mus[mu_index].n_samples();
214 0 : auto received_data = all_evals[mu_index][point_index].size();
215 0 : libmesh_error_msg_if(received_data != n_components * n_samples,
216 : "Recieved " << received_data <<
217 : " evaluated values but expected to receive " << n_components * n_samples);
218 : }
219 : }
220 :
221 : // TODO: move this code for computing the total number of samples
222 : // represented by a std::vector of RBParameters objects to a helper
223 : // function.
224 0 : unsigned int output_size = 0;
225 0 : for (const auto & mu : mus)
226 0 : output_size += mu.n_samples();
227 :
228 0 : output.resize(output_size);
229 :
230 : // We use traditional for-loops here (rather than range-based) so that we can declare and
231 : // increment multiple loop counters all within the local scope of the for-loop.
232 0 : for (auto [mu_index, output_index] = std::make_tuple(0u, 0u); mu_index < mus.size(); ++mu_index)
233 : {
234 0 : auto n_samples = mus[mu_index].n_samples();
235 0 : for (auto mu_sample_idx = 0u; mu_sample_idx < n_samples; ++mu_sample_idx, ++output_index)
236 : {
237 0 : output[output_index].resize(n_points);
238 0 : for (auto point_index : make_range(n_points))
239 : {
240 0 : output[output_index][point_index].resize(n_components);
241 :
242 0 : for (auto comp : make_range(n_components))
243 0 : output[output_index][point_index][comp] = all_evals[mu_index][point_index][n_components*mu_sample_idx + comp];
244 : }
245 : }
246 : }
247 0 : }
248 :
249 0 : void RBParametrizedFunction::side_vectorized_evaluate(const std::vector<RBParameters> & mus,
250 : const VectorizedEvalInput & v,
251 : std::vector<std::vector<std::vector<Number>>> & output)
252 : {
253 0 : LOG_SCOPE("side_vectorized_evaluate()", "RBParametrizedFunction");
254 :
255 0 : output.clear();
256 0 : unsigned int n_points = v.all_xyz.size();
257 :
258 0 : libmesh_error_msg_if(v.sbd_ids.size() != n_points, "Error: invalid vector sizes");
259 0 : libmesh_error_msg_if(requires_xyz_perturbations && (v.all_xyz_perturb.size() != n_points), "Error: invalid vector sizes");
260 :
261 : // Dummy vector to be used when xyz perturbations are not required
262 0 : std::vector<Point> empty_perturbs;
263 :
264 : // The number of components returned by this RBParametrizedFunction
265 0 : auto n_components = this->get_n_components();
266 :
267 : // We first loop over all mus and all n_points, calling side_evaluate()
268 : // for each and storing the results. It is easier to first
269 : // pre-compute all the values before filling output, since, in the
270 : // case of multi-sample RBParameters, the ordering of the loops is a
271 : // bit complicated otherwise.
272 0 : std::vector<std::vector<std::vector<Number>>> all_evals(mus.size());
273 0 : for (auto mu_index : index_range(mus))
274 : {
275 : // Allocate enough space to store all points for the current mu
276 0 : all_evals[mu_index].resize(n_points);
277 0 : for (auto point_index : index_range(all_evals[mu_index]))
278 : {
279 : // Evaluate all samples for the current mu at the current interpolation point
280 0 : all_evals[mu_index][point_index] =
281 0 : this->side_evaluate(mus[mu_index],
282 0 : v.all_xyz[point_index],
283 0 : v.elem_ids[point_index],
284 0 : v.side_indices[point_index],
285 0 : v.qps[point_index],
286 0 : v.sbd_ids[point_index],
287 0 : v.boundary_ids[point_index],
288 0 : requires_xyz_perturbations ? v.all_xyz_perturb[point_index] : empty_perturbs,
289 0 : v.phi_i_qp[point_index]);
290 :
291 : // The vector returned by side_evaluate() should contain:
292 : // n_components * mus[mu_index].n_samples()
293 : // entries. That is, for multi-sample RBParameters objects,
294 : // the vector will be packed with entries as follows:
295 : // [sample0_component0, sample0_component1, ..., sample0_componentN,
296 : // sample1_component0, sample1_component1, ..., sample1_componentN,
297 : // ...
298 : // sampleM_component0, sampleM_component1, ..., sampleM_componentN]
299 0 : auto n_samples = mus[mu_index].n_samples();
300 0 : auto received_data = all_evals[mu_index][point_index].size();
301 0 : libmesh_error_msg_if(received_data != n_components * n_samples,
302 : "Recieved " << received_data <<
303 : " evaluated values but expected to receive " << n_components * n_samples);
304 : }
305 : }
306 :
307 : // TODO: move this code for computing the total number of samples
308 : // represented by a std::vector of RBParameters objects to a helper
309 : // function.
310 0 : unsigned int output_size = 0;
311 0 : for (const auto & mu : mus)
312 0 : output_size += mu.n_samples();
313 :
314 0 : output.resize(output_size);
315 :
316 : // We use traditional for-loops here (rather than range-based) so that we can declare and
317 : // increment multiple loop counters all within the local scope of the for-loop.
318 0 : for (auto [mu_index, output_index] = std::make_tuple(0u, 0u); mu_index < mus.size(); ++mu_index)
319 : {
320 0 : auto n_samples = mus[mu_index].n_samples();
321 0 : for (auto mu_sample_idx = 0u; mu_sample_idx < n_samples; ++mu_sample_idx, ++output_index)
322 : {
323 0 : output[output_index].resize(n_points);
324 0 : for (auto point_index : make_range(n_points))
325 : {
326 0 : output[output_index][point_index].resize(n_components);
327 :
328 0 : for (auto comp : make_range(n_components))
329 0 : output[output_index][point_index][comp] = all_evals[mu_index][point_index][n_components*mu_sample_idx + comp];
330 : }
331 : }
332 : }
333 0 : }
334 :
335 0 : void RBParametrizedFunction::node_vectorized_evaluate(const std::vector<RBParameters> & mus,
336 : const VectorizedEvalInput & v,
337 : std::vector<std::vector<std::vector<Number>>> & output)
338 : {
339 0 : LOG_SCOPE("node_vectorized_evaluate()", "RBParametrizedFunction");
340 :
341 0 : output.clear();
342 0 : unsigned int n_points = v.all_xyz.size();
343 :
344 : // The number of components returned by this RBParametrizedFunction
345 0 : auto n_components = this->get_n_components();
346 :
347 : // We first loop over all mus and all n_points, calling node_evaluate()
348 : // for each and storing the results. It is easier to first
349 : // pre-compute all the values before filling output, since, in the
350 : // case of multi-sample RBParameters, the ordering of the loops is a
351 : // bit complicated otherwise.
352 0 : std::vector<std::vector<std::vector<Number>>> all_evals(mus.size());
353 0 : for (auto mu_index : index_range(mus))
354 : {
355 : // Allocate enough space to store all points for the current mu
356 0 : all_evals[mu_index].resize(n_points);
357 0 : for (auto point_index : index_range(all_evals[mu_index]))
358 : {
359 : // Evaluate all samples for the current mu at the current interpolation point
360 0 : all_evals[mu_index][point_index] =
361 0 : this->node_evaluate(mus[mu_index],
362 0 : v.all_xyz[point_index],
363 0 : v.node_ids[point_index],
364 0 : v.boundary_ids[point_index]);
365 :
366 : // The vector returned by node_evaluate() should contain:
367 : // n_components * mus[mu_index].n_samples()
368 : // entries. That is, for multi-sample RBParameters objects,
369 : // the vector will be packed with entries as follows:
370 : // [sample0_component0, sample0_component1, ..., sample0_componentN,
371 : // sample1_component0, sample1_component1, ..., sample1_componentN,
372 : // ...
373 : // sampleM_component0, sampleM_component1, ..., sampleM_componentN]
374 0 : auto n_samples = mus[mu_index].n_samples();
375 0 : auto received_data = all_evals[mu_index][point_index].size();
376 0 : libmesh_error_msg_if(received_data != n_components * n_samples,
377 : "Recieved " << received_data <<
378 : " evaluated values but expected to receive " << n_components * n_samples);
379 : }
380 : }
381 :
382 : // TODO: move this code for computing the total number of samples
383 : // represented by a std::vector of RBParameters objects to a helper
384 : // function.
385 0 : unsigned int output_size = 0;
386 0 : for (const auto & mu : mus)
387 0 : output_size += mu.n_samples();
388 :
389 0 : output.resize(output_size);
390 :
391 : // We use traditional for-loops here (rather than range-based) so that we can declare and
392 : // increment multiple loop counters all within the local scope of the for-loop.
393 0 : for (auto [mu_index, output_index] = std::make_tuple(0u, 0u); mu_index < mus.size(); ++mu_index)
394 : {
395 0 : auto n_samples = mus[mu_index].n_samples();
396 0 : for (auto mu_sample_idx = 0u; mu_sample_idx < n_samples; ++mu_sample_idx, ++output_index)
397 : {
398 0 : output[output_index].resize(n_points);
399 0 : for (auto point_index : make_range(n_points))
400 : {
401 0 : output[output_index][point_index].resize(n_components);
402 :
403 0 : for (auto comp : make_range(n_components))
404 0 : output[output_index][point_index][comp] = all_evals[mu_index][point_index][n_components*mu_sample_idx + comp];
405 : }
406 : }
407 : }
408 0 : }
409 :
410 0 : void RBParametrizedFunction::preevaluate_parametrized_function_on_mesh(const RBParameters & mu,
411 : const std::unordered_map<dof_id_type, std::vector<Point>> & all_xyz,
412 : const std::unordered_map<dof_id_type, subdomain_id_type> & sbd_ids,
413 : const std::unordered_map<dof_id_type, std::vector<std::vector<Point>> > & all_xyz_perturb,
414 : const System & sys)
415 : {
416 0 : mesh_to_preevaluated_values_map.clear();
417 :
418 0 : unsigned int n_elems = all_xyz.size();
419 0 : unsigned int n_points = 0;
420 0 : for (const auto & xyz_pair : all_xyz)
421 : {
422 0 : const std::vector<Point> & xyz_vec = xyz_pair.second;
423 0 : n_points += xyz_vec.size();
424 : }
425 :
426 0 : VectorizedEvalInput v;
427 0 : v.all_xyz.resize(n_points);
428 0 : v.elem_ids.resize(n_points);
429 0 : v.qps.resize(n_points);
430 0 : v.sbd_ids.resize(n_points);
431 0 : v.all_xyz_perturb.resize(n_points);
432 0 : v.phi_i_qp.resize(n_points);
433 0 : v.elem_types.resize(n_points);
434 :
435 0 : if (requires_all_elem_qp_data)
436 : {
437 0 : v.elem_id_to_local_index.clear();
438 0 : v.JxW_all_qp.resize(n_elems);
439 0 : v.phi_i_all_qp.resize(n_elems);
440 : }
441 0 : if (requires_all_elem_center_data)
442 : {
443 0 : v.dxyzdxi_elem_center.resize(n_elems);
444 0 : v.dxyzdeta_elem_center.resize(n_elems);
445 : // At the time of writing, the quadrature order is only used in conjunction with element
446 : // center data so we should not compute it elsewhere for now.
447 0 : v.qrule_orders.resize(n_elems);
448 : }
449 :
450 : // Empty vector to be used when xyz perturbations are not required
451 0 : std::vector<Point> empty_perturbs;
452 :
453 : // In order to compute phi_i_qp, we initialize a FEMContext
454 0 : FEMContext con(sys);
455 0 : for (auto dim : con.elem_dimensions())
456 : {
457 0 : auto fe = con.get_element_fe(/*var=*/0, dim);
458 0 : fe->get_JxW();
459 0 : fe->get_phi();
460 0 : fe->get_dxyzdxi();
461 0 : fe->get_dxyzdeta();
462 : }
463 :
464 0 : unsigned int counter = 0;
465 0 : unsigned int elem_counter = 0;
466 0 : for (const auto & [elem_id, xyz_vec] : all_xyz)
467 : {
468 0 : subdomain_id_type subdomain_id = libmesh_map_find(sbd_ids, elem_id);
469 :
470 : // The amount of data to be stored for each component
471 0 : auto n_qp = xyz_vec.size();
472 0 : mesh_to_preevaluated_values_map[elem_id].resize(n_qp);
473 :
474 : // Also initialize phi in order to compute phi_i_qp
475 0 : const Elem & elem_ref = sys.get_mesh().elem_ref(elem_id);
476 0 : con.pre_fe_reinit(sys, &elem_ref);
477 :
478 0 : auto elem_fe = con.get_element_fe(/*var=*/0, elem_ref.dim());
479 0 : const std::vector<std::vector<Real>> & phi = elem_fe->get_phi();
480 0 : const std::vector<Real> & JxW = elem_fe->get_JxW();
481 0 : const auto & dxyzdxi = elem_fe->get_dxyzdxi();
482 0 : const auto & dxyzdeta = elem_fe->get_dxyzdeta();
483 :
484 0 : elem_fe->reinit(&elem_ref);
485 :
486 0 : for (auto qp : index_range(xyz_vec))
487 : {
488 0 : mesh_to_preevaluated_values_map[elem_id][qp] = counter;
489 :
490 0 : v.all_xyz[counter] = xyz_vec[qp];
491 0 : v.elem_ids[counter] = elem_id;
492 0 : v.qps[counter] = qp;
493 0 : v.sbd_ids[counter] = subdomain_id;
494 0 : v.elem_types[counter] = elem_ref.type();
495 :
496 0 : v.phi_i_qp[counter].resize(phi.size());
497 0 : for(auto i : index_range(phi))
498 0 : v.phi_i_qp[counter][i] = phi[i][qp];
499 :
500 0 : if (requires_xyz_perturbations)
501 : {
502 : const auto & qps_and_perturbs =
503 0 : libmesh_map_find(all_xyz_perturb, elem_id);
504 0 : libmesh_error_msg_if(qp >= qps_and_perturbs.size(), "Error: Invalid qp");
505 :
506 0 : v.all_xyz_perturb[counter] = qps_and_perturbs[qp];
507 : }
508 : else
509 : {
510 0 : v.all_xyz_perturb[counter] = empty_perturbs;
511 : }
512 :
513 0 : if (requires_all_elem_qp_data)
514 : {
515 0 : if (v.elem_id_to_local_index.count(elem_id) == 0)
516 : {
517 : // In this case we store data for all qps on this element
518 : // at each point.
519 0 : v.JxW_all_qp[elem_counter].resize(JxW.size());
520 0 : for(auto i : index_range(JxW))
521 0 : v.JxW_all_qp[elem_counter][i] = JxW[i];
522 :
523 0 : v.phi_i_all_qp[elem_counter].resize(phi.size());
524 0 : for(auto i : index_range(phi))
525 : {
526 0 : v.phi_i_all_qp[elem_counter][i].resize(phi[i].size());
527 0 : for(auto j : index_range(phi[i]))
528 0 : v.phi_i_all_qp[elem_counter][i][j] = phi[i][j];
529 : }
530 0 : v.elem_id_to_local_index[elem_id] = elem_counter;
531 : }
532 : }
533 :
534 0 : counter++;
535 : }
536 :
537 : // Here we presume that if requires_all_elem_center_data is set to true
538 : // then requires_all_elem_qp_data is also set to true so that elem_id_to_local_index entry
539 : // for this elem has already been set.
540 0 : if (requires_all_elem_center_data)
541 : {
542 : // Get data derivatives at vertex average
543 0 : std::vector<Point> nodes = { elem_ref.reference_elem()->vertex_average() };
544 0 : elem_fe->reinit (&elem_ref, &nodes);
545 : // Set qrule_order here to prevent calling getter multiple times.
546 0 : Order qrule_order = con.get_element_qrule().get_order();
547 :
548 : // We add qrule_order in this loop as it is used in conjunction with elem center
549 : // quantities for now.
550 0 : v.qrule_orders[elem_counter] = qrule_order;
551 0 : Point dxyzdxi_pt, dxyzdeta_pt;
552 0 : if (con.get_elem_dim()>0)
553 0 : dxyzdxi_pt = dxyzdxi[0];
554 0 : if (con.get_elem_dim()>1)
555 0 : dxyzdeta_pt = dxyzdeta[0];
556 : // Here we do an implicit conversion from RealGradient which is a VectorValue<Real>
557 : // which in turn is a TypeVector<T> to a Point which is a TypeVector<Real>.
558 : // They are essentially the same thing. This helps us limiting the number of includes
559 : // in serialization and deserialization as RealGradient is a typedef and we cannot
560 : // forward declare typedefs. As a result we leverage the fact that point.h is already
561 : // included in most places we need RealGradient.
562 0 : v.dxyzdxi_elem_center[elem_counter] = dxyzdxi_pt;
563 0 : v.dxyzdeta_elem_center[elem_counter] = dxyzdeta_pt;
564 : }
565 0 : elem_counter++;
566 : }
567 :
568 0 : std::vector<RBParameters> mus {mu};
569 0 : vectorized_evaluate(mus, v, preevaluated_values);
570 :
571 0 : preevaluate_parametrized_function_cleanup();
572 0 : }
573 :
574 0 : void RBParametrizedFunction::preevaluate_parametrized_function_on_mesh_sides(const RBParameters & mu,
575 : const std::map<std::pair<dof_id_type,unsigned int>, std::vector<Point>> & side_all_xyz,
576 : const std::map<std::pair<dof_id_type,unsigned int>, subdomain_id_type> & sbd_ids,
577 : const std::map<std::pair<dof_id_type,unsigned int>, boundary_id_type> & side_boundary_ids,
578 : const std::map<std::pair<dof_id_type,unsigned int>, unsigned int> & side_types,
579 : const std::map<std::pair<dof_id_type,unsigned int>, std::vector<std::vector<Point>> > & side_all_xyz_perturb,
580 : const System & sys)
581 : {
582 0 : mesh_to_preevaluated_side_values_map.clear();
583 :
584 0 : unsigned int n_points = 0;
585 0 : for (const auto & xyz_pair : side_all_xyz)
586 : {
587 0 : const std::vector<Point> & xyz_vec = xyz_pair.second;
588 0 : n_points += xyz_vec.size();
589 : }
590 :
591 0 : VectorizedEvalInput v;
592 0 : v.all_xyz.resize(n_points);
593 0 : v.elem_ids.resize(n_points);
594 0 : v.side_indices.resize(n_points);
595 0 : v.qps.resize(n_points);
596 0 : v.sbd_ids.resize(n_points);
597 0 : v.boundary_ids.resize(n_points);
598 0 : v.all_xyz_perturb.resize(n_points);
599 0 : v.phi_i_qp.resize(n_points);
600 :
601 : // Empty vector to be used when xyz perturbations are not required
602 0 : std::vector<Point> empty_perturbs;
603 :
604 : // In order to compute phi_i_qp, we initialize a FEMContext
605 0 : FEMContext con(sys);
606 0 : for (auto dim : con.elem_dimensions())
607 : {
608 0 : auto fe = con.get_element_fe(/*var=*/0, dim);
609 0 : fe->get_phi();
610 :
611 0 : auto side_fe = con.get_side_fe(/*var=*/0, dim);
612 0 : side_fe->get_phi();
613 : }
614 :
615 0 : unsigned int counter = 0;
616 0 : for (const auto & xyz_pair : side_all_xyz)
617 : {
618 0 : auto elem_side_pair = xyz_pair.first;
619 0 : dof_id_type elem_id = elem_side_pair.first;
620 0 : unsigned int side_index = elem_side_pair.second;
621 :
622 0 : const std::vector<Point> & xyz_vec = xyz_pair.second;
623 :
624 0 : subdomain_id_type subdomain_id = libmesh_map_find(sbd_ids, elem_side_pair);
625 0 : boundary_id_type boundary_id = libmesh_map_find(side_boundary_ids, elem_side_pair);
626 0 : unsigned int side_type = libmesh_map_find(side_types, elem_side_pair);
627 :
628 : // The amount of data to be stored for each component
629 0 : auto n_qp = xyz_vec.size();
630 0 : mesh_to_preevaluated_side_values_map[elem_side_pair].resize(n_qp);
631 :
632 0 : const Elem & elem_ref = sys.get_mesh().elem_ref(elem_id);
633 0 : con.pre_fe_reinit(sys, &elem_ref);
634 :
635 : // side_type == 0 --> standard side
636 : // side_type == 1 --> shellface
637 0 : if (side_type == 0)
638 : {
639 0 : std::unique_ptr<const Elem> elem_side;
640 0 : elem_ref.build_side_ptr(elem_side, side_index);
641 :
642 0 : auto side_fe = con.get_side_fe(/*var=*/0, elem_ref.dim());
643 0 : side_fe->reinit(&elem_ref, side_index);
644 :
645 0 : const std::vector<std::vector<Real>> & phi = side_fe->get_phi();
646 0 : for (auto qp : index_range(xyz_vec))
647 : {
648 0 : mesh_to_preevaluated_side_values_map[elem_side_pair][qp] = counter;
649 :
650 0 : v.all_xyz[counter] = xyz_vec[qp];
651 0 : v.elem_ids[counter] = elem_side_pair.first;
652 0 : v.side_indices[counter] = elem_side_pair.second;
653 0 : v.qps[counter] = qp;
654 0 : v.sbd_ids[counter] = subdomain_id;
655 0 : v.boundary_ids[counter] = boundary_id;
656 :
657 0 : v.phi_i_qp[counter].resize(phi.size());
658 0 : for(auto i : index_range(phi))
659 0 : v.phi_i_qp[counter][i] = phi[i][qp];
660 :
661 0 : if (requires_xyz_perturbations)
662 : {
663 : const auto & qps_and_perturbs =
664 0 : libmesh_map_find(side_all_xyz_perturb, elem_side_pair);
665 0 : libmesh_error_msg_if(qp >= qps_and_perturbs.size(), "Error: Invalid qp");
666 :
667 0 : v.all_xyz_perturb[counter] = qps_and_perturbs[qp];
668 : }
669 : else
670 : {
671 0 : v.all_xyz_perturb[counter] = empty_perturbs;
672 : }
673 0 : counter++;
674 : }
675 0 : }
676 0 : else if (side_type == 1)
677 : {
678 0 : auto elem_fe = con.get_element_fe(/*var=*/0, elem_ref.dim());
679 0 : const std::vector<std::vector<Real>> & phi = elem_fe->get_phi();
680 :
681 0 : elem_fe->reinit(&elem_ref);
682 :
683 0 : for (auto qp : index_range(xyz_vec))
684 : {
685 0 : mesh_to_preevaluated_side_values_map[elem_side_pair][qp] = counter;
686 :
687 0 : v.all_xyz[counter] = xyz_vec[qp];
688 0 : v.elem_ids[counter] = elem_side_pair.first;
689 0 : v.side_indices[counter] = elem_side_pair.second;
690 0 : v.qps[counter] = qp;
691 0 : v.sbd_ids[counter] = subdomain_id;
692 0 : v.boundary_ids[counter] = boundary_id;
693 :
694 0 : v.phi_i_qp[counter].resize(phi.size());
695 0 : for(auto i : index_range(phi))
696 0 : v.phi_i_qp[counter][i] = phi[i][qp];
697 :
698 0 : if (requires_xyz_perturbations)
699 : {
700 : const auto & qps_and_perturbs =
701 0 : libmesh_map_find(side_all_xyz_perturb, elem_side_pair);
702 0 : libmesh_error_msg_if(qp >= qps_and_perturbs.size(), "Error: Invalid qp");
703 :
704 0 : v.all_xyz_perturb[counter] = qps_and_perturbs[qp];
705 : }
706 : else
707 : {
708 0 : v.all_xyz_perturb[counter] = empty_perturbs;
709 : }
710 0 : counter++;
711 : }
712 : }
713 : else
714 0 : libmesh_error_msg ("Unrecognized side_type: " << side_type);
715 : }
716 :
717 0 : std::vector<RBParameters> mus {mu};
718 0 : side_vectorized_evaluate(mus, v, preevaluated_values);
719 :
720 0 : preevaluate_parametrized_function_cleanup();
721 0 : }
722 :
723 0 : void RBParametrizedFunction::preevaluate_parametrized_function_on_mesh_nodes(const RBParameters & mu,
724 : const std::unordered_map<dof_id_type, Point> & all_xyz,
725 : const std::unordered_map<dof_id_type, boundary_id_type> & node_boundary_ids,
726 : const System & /*sys*/)
727 : {
728 0 : mesh_to_preevaluated_node_values_map.clear();
729 :
730 0 : unsigned int n_points = all_xyz.size();
731 :
732 0 : VectorizedEvalInput v;
733 0 : v.all_xyz.resize(n_points);
734 0 : v.node_ids.resize(n_points);
735 0 : v.boundary_ids.resize(n_points);
736 :
737 : // Empty vector to be used when xyz perturbations are not required
738 0 : std::vector<Point> empty_perturbs;
739 :
740 0 : unsigned int counter = 0;
741 0 : for (const auto & [node_id, p] : all_xyz)
742 : {
743 0 : boundary_id_type boundary_id = libmesh_map_find(node_boundary_ids, node_id);
744 :
745 0 : mesh_to_preevaluated_node_values_map[node_id] = counter;
746 :
747 0 : v.all_xyz[counter] = p;
748 0 : v.node_ids[counter] = node_id;
749 0 : v.boundary_ids[counter] = boundary_id;
750 :
751 0 : counter++;
752 : }
753 :
754 0 : std::vector<RBParameters> mus {mu};
755 0 : node_vectorized_evaluate(mus, v, preevaluated_values);
756 :
757 0 : preevaluate_parametrized_function_cleanup();
758 0 : }
759 :
760 0 : Number RBParametrizedFunction::lookup_preevaluated_value_on_mesh(unsigned int comp,
761 : dof_id_type elem_id,
762 : unsigned int qp) const
763 : {
764 : const std::vector<unsigned int> & indices_at_qps =
765 0 : libmesh_map_find(mesh_to_preevaluated_values_map, elem_id);
766 :
767 0 : libmesh_error_msg_if(qp >= indices_at_qps.size(), "Error: invalid qp");
768 :
769 0 : unsigned int index = indices_at_qps[qp];
770 0 : libmesh_error_msg_if(preevaluated_values.size() != 1, "Error: we expect only one parameter index");
771 0 : libmesh_error_msg_if(index >= preevaluated_values[0].size(), "Error: invalid index");
772 :
773 0 : return preevaluated_values[0][index][comp];
774 : }
775 :
776 0 : Number RBParametrizedFunction::lookup_preevaluated_side_value_on_mesh(unsigned int comp,
777 : dof_id_type elem_id,
778 : unsigned int side_index,
779 : unsigned int qp) const
780 : {
781 : const std::vector<unsigned int> & indices_at_qps =
782 0 : libmesh_map_find(mesh_to_preevaluated_side_values_map, std::make_pair(elem_id,side_index));
783 :
784 0 : libmesh_error_msg_if(qp >= indices_at_qps.size(), "Error: invalid qp");
785 :
786 0 : unsigned int index = indices_at_qps[qp];
787 0 : libmesh_error_msg_if(preevaluated_values.size() != 1, "Error: we expect only one parameter index");
788 0 : libmesh_error_msg_if(index >= preevaluated_values[0].size(), "Error: invalid index");
789 :
790 0 : return preevaluated_values[0][index][comp];
791 : }
792 :
793 0 : Number RBParametrizedFunction::lookup_preevaluated_node_value_on_mesh(unsigned int comp,
794 : dof_id_type node_id) const
795 : {
796 : unsigned int index =
797 0 : libmesh_map_find(mesh_to_preevaluated_node_values_map, node_id);
798 :
799 0 : libmesh_error_msg_if(preevaluated_values.size() != 1, "Error: we expect only one parameter index");
800 0 : libmesh_error_msg_if(index >= preevaluated_values[0].size(), "Error: invalid index");
801 :
802 0 : return preevaluated_values[0][index][comp];
803 : }
804 :
805 0 : void RBParametrizedFunction::initialize_lookup_table()
806 : {
807 : // No-op by default, override in subclasses as needed
808 0 : }
809 :
810 0 : Number RBParametrizedFunction::get_parameter_independent_data(const std::string & property_name,
811 : subdomain_id_type sbd_id) const
812 : {
813 0 : return libmesh_map_find(libmesh_map_find(_parameter_independent_data, property_name), sbd_id);
814 : }
815 :
816 0 : void RBParametrizedFunction::get_spatial_indices(std::vector<std::vector<unsigned int>> & /*spatial_indices*/,
817 : const VectorizedEvalInput & /*v*/)
818 : {
819 : // No-op by default
820 0 : }
821 :
822 0 : void RBParametrizedFunction::initialize_spatial_indices(const std::vector<std::vector<unsigned int>> & /*spatial_indices*/,
823 : const VectorizedEvalInput & /*v*/)
824 : {
825 : // No-op by default
826 0 : }
827 :
828 0 : void RBParametrizedFunction::preevaluate_parametrized_function_cleanup()
829 : {
830 : // No-op by default
831 0 : }
832 :
833 0 : const std::unordered_map<std::string, std::set<dof_id_type>> & RBParametrizedFunction::get_rb_property_map() const
834 : {
835 0 : return _rb_property_map;
836 : }
837 :
838 0 : void RBParametrizedFunction::add_rb_property_map_entry(std::string & property_name, std::set<dof_id_type> & entity_ids)
839 : {
840 0 : bool insert_succeed = _rb_property_map.insert({property_name, entity_ids}).second;
841 0 : libmesh_error_msg_if(!insert_succeed, "Entry already added, duplicate detected.");
842 0 : }
843 :
844 0 : void RBParametrizedFunction::add_interpolation_data_to_rb_property_map(
845 : const Parallel::Communicator & /*comm*/,
846 : std::unordered_map<std::string, std::set<dof_id_type>> & /*rb_property_map*/,
847 : dof_id_type /*elem_id*/)
848 : {
849 : // No-op by default
850 0 : }
851 :
852 0 : const std::set<boundary_id_type> & RBParametrizedFunction::get_parametrized_function_boundary_ids() const
853 : {
854 0 : return _parametrized_function_boundary_ids;
855 : }
856 :
857 0 : void RBParametrizedFunction::set_parametrized_function_boundary_ids(const std::set<boundary_id_type> & boundary_ids, bool is_nodal_boundary)
858 : {
859 0 : _parametrized_function_boundary_ids = boundary_ids;
860 0 : _is_nodal_boundary = is_nodal_boundary;
861 0 : }
862 :
863 0 : bool RBParametrizedFunction::on_mesh_sides() const
864 : {
865 0 : return !get_parametrized_function_boundary_ids().empty() && !_is_nodal_boundary;
866 : }
867 :
868 0 : bool RBParametrizedFunction::on_mesh_nodes() const
869 : {
870 0 : return !get_parametrized_function_boundary_ids().empty() && _is_nodal_boundary;
871 : }
872 :
873 : }
|
