Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #include "SolutionUserObjectBase.h"
11 :
12 : // MOOSE includes
13 : #include "ConsoleUtils.h"
14 : #include "MooseError.h"
15 : #include "MooseMesh.h"
16 : #include "MooseUtils.h"
17 : #include "MooseVariableFE.h"
18 : #include "RotationMatrix.h"
19 : #include "Function.h"
20 :
21 : // libMesh includes
22 : #include "libmesh/equation_systems.h"
23 : #include "libmesh/mesh_function.h"
24 : #include "libmesh/numeric_vector.h"
25 : #include "libmesh/nonlinear_implicit_system.h"
26 : #include "libmesh/transient_system.h"
27 : #include "libmesh/parallel_mesh.h"
28 : #include "libmesh/serial_mesh.h"
29 : #include "libmesh/exodusII_io.h"
30 : #include "libmesh/exodusII_io_helper.h"
31 : #include "libmesh/nemesis_io.h"
32 : #include "libmesh/nemesis_io_helper.h"
33 : #include "libmesh/enum_xdr_mode.h"
34 : #include "libmesh/string_to_enum.h"
35 :
36 : InputParameters
37 4063 : SolutionUserObjectBase::validParams()
38 : {
39 : // Get the input parameters from the parent class
40 4063 : InputParameters params = GeneralUserObject::validParams();
41 :
42 : // Add required parameters
43 16252 : params.addRequiredParam<MeshFileName>(
44 : "mesh", "The name of the mesh file (must be xda/xdr, exodusII or nemesis file).");
45 12189 : params.addParam<std::vector<std::string>>(
46 : "system_variables",
47 8126 : std::vector<std::string>(),
48 : "The name of the nodal and elemental variables from the file you want to use for values");
49 12189 : params.addParam<bool>(
50 : "force_replicated_source_mesh",
51 8126 : false,
52 : "Whether to force the serialization of the source mesh. This can be useful for discontinuous "
53 : "variables or if the partitioning of the source mesh and the current mesh are different");
54 :
55 : // When using XDA/XDR files the following must be defined
56 16252 : params.addParam<FileName>(
57 : "es",
58 : "<not supplied>",
59 : "The name of the file holding the equation system info in xda/xdr format (xda/xdr only).");
60 16252 : params.addParam<std::string>(
61 : "system",
62 : "nl0",
63 : "The name of the system to pull values out of (xda/xdr only). The default name for the "
64 : "nonlinear system is 'nl0', auxiliary system is 'aux0'");
65 :
66 : // When using ExodusII a specific time is extracted
67 16252 : params.addParam<std::string>("timestep",
68 : "Index of the single timestep used or \"LATEST\" for "
69 : "the last timestep (exodusII only). If not supplied, "
70 : "time interpolation will occur.");
71 :
72 : // Choose the interpolation order for incoming nodal variables
73 12189 : params.addParam<MooseEnum>("nodal_variable_order",
74 20315 : MooseEnum("FIRST SECOND", "FIRST"),
75 : "Specifies the order of the nodal solution data.");
76 :
77 : // Add ability to perform coordinate transformation: scale, factor
78 12189 : params.addParam<std::vector<Real>>(
79 8126 : "scale", std::vector<Real>(LIBMESH_DIM, 1), "Scale factor for points in the simulation");
80 12189 : params.addParam<std::vector<Real>>("scale_multiplier",
81 8126 : std::vector<Real>(LIBMESH_DIM, 1),
82 : "Scale multiplying factor for points in the simulation");
83 12189 : params.addParam<std::vector<Real>>("translation",
84 8126 : std::vector<Real>(LIBMESH_DIM, 0),
85 : "Translation factors for x,y,z coordinates of the simulation");
86 12189 : params.addParam<RealVectorValue>("rotation0_vector",
87 8126 : RealVectorValue(0, 0, 1),
88 : "Vector about which to rotate points of the simulation.");
89 12189 : params.addParam<Real>(
90 : "rotation0_angle",
91 8126 : 0.0,
92 : "Anticlockwise rotation angle (in degrees) to use for rotation about rotation0_vector.");
93 12189 : params.addParam<RealVectorValue>("rotation1_vector",
94 8126 : RealVectorValue(0, 0, 1),
95 : "Vector about which to rotate points of the simulation.");
96 12189 : params.addParam<Real>(
97 : "rotation1_angle",
98 8126 : 0.0,
99 : "Anticlockwise rotation angle (in degrees) to use for rotation about rotation1_vector.");
100 :
101 : // following lines build the default_transformation_order
102 : MultiMooseEnum default_transformation_order(
103 8126 : "rotation0 translation scale rotation1 scale_multiplier", "translation scale");
104 16252 : params.addParam<MultiMooseEnum>(
105 : "transformation_order",
106 : default_transformation_order,
107 : "The order to perform the operations in. Define R0 to be the rotation matrix encoded by "
108 : "rotation0_vector and rotation0_angle. Similarly for R1. Denote the scale by s, the "
109 : "scale_multiplier by m, and the translation by t. Then, given a point x in the simulation, "
110 : "if transformation_order = 'rotation0 scale_multiplier translation scale rotation1' then "
111 : "form p = R1*(R0*x*m - t)/s. Then the values provided by the SolutionUserObjectBase at "
112 : "point x "
113 : "in the simulation are the variable values at point p in the mesh.");
114 16252 : params.addParamNamesToGroup("scale scale_multiplier translation rotation0_vector rotation0_angle "
115 : "rotation1_angle transformation_order",
116 : "Coordinate system transformation");
117 4063 : params.addClassDescription("Reads a variable from a mesh in one simulation to another");
118 8126 : return params;
119 4063 : }
120 :
121 : // Static mutex definition
122 : Threads::spin_mutex SolutionUserObjectBase::_solution_user_object_mutex;
123 :
124 501 : SolutionUserObjectBase::SolutionUserObjectBase(const InputParameters & parameters)
125 : : GeneralUserObject(parameters),
126 1503 : _file_type(MooseEnum("xda=0 exodusII=1 xdr=2 nemesis=3")),
127 1002 : _mesh_file(getParam<MeshFileName>("mesh")),
128 1002 : _es_file(getParam<FileName>("es")),
129 1002 : _system_name(getParam<std::string>("system")),
130 1002 : _system_variables(getParam<std::vector<std::string>>("system_variables")),
131 501 : _exodus_time_index(-1),
132 501 : _interpolate_times(false),
133 501 : _system(nullptr),
134 501 : _system2(nullptr),
135 501 : _interpolation_time(0.0),
136 501 : _interpolation_factor(0.0),
137 501 : _exodus_times(nullptr),
138 501 : _exodus_index1(0),
139 501 : _exodus_index2(0),
140 1002 : _nodal_variable_order(getParam<MooseEnum>("nodal_variable_order")),
141 1002 : _scale(getParam<std::vector<Real>>("scale")),
142 1002 : _scale_multiplier(getParam<std::vector<Real>>("scale_multiplier")),
143 1002 : _translation(getParam<std::vector<Real>>("translation")),
144 1002 : _rotation0_vector(getParam<RealVectorValue>("rotation0_vector")),
145 1002 : _rotation0_angle(getParam<Real>("rotation0_angle")),
146 501 : _r0(RealTensorValue()),
147 1002 : _rotation1_vector(getParam<RealVectorValue>("rotation1_vector")),
148 1002 : _rotation1_angle(getParam<Real>("rotation1_angle")),
149 501 : _r1(RealTensorValue()),
150 1002 : _transformation_order(getParam<MultiMooseEnum>("transformation_order")),
151 1002 : _force_replicated_source(getParam<bool>("force_replicated_source_mesh")),
152 1503 : _initialized(false)
153 : {
154 : // form rotation matrices with the specified angles
155 501 : Real halfPi = libMesh::pi / 2.0;
156 : Real a;
157 : Real b;
158 :
159 501 : a = std::cos(halfPi * -_rotation0_angle / 90.0);
160 501 : b = std::sin(halfPi * -_rotation0_angle / 90.0);
161 : // the following is an anticlockwise rotation about z
162 501 : RealTensorValue rot0_z(a, -b, 0, b, a, 0, 0, 0, 1);
163 : // form the rotation matrix that will take rotation0_vector to the z axis
164 501 : RealTensorValue vec0_to_z = RotationMatrix::rotVecToZ(_rotation0_vector);
165 : // _r0 is then: rotate points so vec0 lies along z; then rotate about angle0; then rotate points
166 : // back
167 501 : _r0 = vec0_to_z.transpose() * (rot0_z * vec0_to_z);
168 :
169 501 : a = std::cos(halfPi * -_rotation1_angle / 90.0);
170 501 : b = std::sin(halfPi * -_rotation1_angle / 90.0);
171 : // the following is an anticlockwise rotation about z
172 501 : RealTensorValue rot1_z(a, -b, 0, b, a, 0, 0, 0, 1);
173 : // form the rotation matrix that will take rotation1_vector to the z axis
174 501 : RealTensorValue vec1_to_z = RotationMatrix::rotVecToZ(_rotation1_vector);
175 : // _r1 is then: rotate points so vec1 lies along z; then rotate about angle1; then rotate points
176 : // back
177 501 : _r1 = vec1_to_z.transpose() * (rot1_z * vec1_to_z);
178 :
179 1905 : if (isParamValid("timestep") && getParam<std::string>("timestep") == "-1")
180 0 : mooseError("A \"timestep\" of -1 is no longer supported for interpolation. Instead simply "
181 : "remove this parameter altogether for interpolation");
182 :
183 : // check for possible inconsistencies between mesh and input data
184 501 : if (_fe_problem.mesh().hasSecondOrderElements())
185 24 : mooseInfo(
186 : "The mesh has second-order elements, be sure to set 'nodal_variable_order' if needed.");
187 501 : }
188 :
189 : void
190 110 : SolutionUserObjectBase::readXda()
191 : {
192 330 : if (!isParamSetByUser("es"))
193 0 : paramError("es", "Equation system file (.xda or .xdr) should have been specified");
194 :
195 : // Check that the required files exist
196 110 : MooseUtils::checkFileReadable(_es_file);
197 110 : MooseUtils::checkFileReadable(_mesh_file);
198 :
199 : // Read the libmesh::mesh from the xda file
200 110 : _mesh->read(_mesh_file);
201 :
202 : // Create the libmesh::EquationSystems
203 110 : _es = std::make_unique<EquationSystems>(*_mesh);
204 :
205 : // Use new read syntax (binary)
206 110 : if (_file_type == "xdr")
207 0 : _es->read(_es_file,
208 : libMesh::DECODE,
209 : EquationSystems::READ_HEADER | EquationSystems::READ_DATA |
210 : EquationSystems::READ_ADDITIONAL_DATA);
211 :
212 : // Use new read syntax
213 110 : else if (_file_type == "xda")
214 110 : _es->read(_es_file,
215 : libMesh::READ,
216 : EquationSystems::READ_HEADER | EquationSystems::READ_DATA |
217 : EquationSystems::READ_ADDITIONAL_DATA);
218 :
219 : // This should never occur, just in case produce an error
220 : else
221 0 : mooseError("Failed to determine proper read method for XDA/XDR equation system file: ",
222 0 : _es_file);
223 :
224 : // Update and store the EquationSystems name locally
225 110 : _es->update();
226 110 : _system = &_es->get_system(_system_name);
227 110 : }
228 :
229 : void
230 391 : SolutionUserObjectBase::readExodusIIOrNemesis()
231 : {
232 : // Define a default system name
233 391 : if (_system_name == "")
234 0 : _system_name = "SolutionUserObjectSystem";
235 :
236 : // Read the Exodus or Nemesis file
237 391 : if (_file_type == "exodusII")
238 : {
239 367 : _exodusII_io = std::make_unique<libMesh::ExodusII_IO>(*_mesh);
240 367 : _exodusII_io->read(_mesh_file);
241 367 : _exodus_times = &_exodusII_io->get_time_steps();
242 : }
243 : else
244 : {
245 24 : _nemesis_io = std::make_unique<libMesh::Nemesis_IO>(*_mesh);
246 24 : _nemesis_io->read(_mesh_file);
247 24 : _exodus_times = &_nemesis_io->get_time_steps();
248 : }
249 : // Same routine for Nemesis
250 391 : readBlockIdMapFromExodusII();
251 :
252 1173 : if (isParamValid("timestep"))
253 : {
254 402 : std::string s_timestep = getParam<std::string>("timestep");
255 201 : int n_steps = (_file_type == "exodusII") ? _exodusII_io->get_num_time_steps()
256 0 : : _nemesis_io->get_num_time_steps();
257 201 : if (s_timestep == "LATEST")
258 54 : _exodus_time_index = n_steps;
259 : else
260 : {
261 147 : std::istringstream ss(s_timestep);
262 147 : if (!((ss >> _exodus_time_index) && ss.eof()) || _exodus_time_index > n_steps)
263 0 : mooseError("Invalid value passed as \"timestep\". Expected \"LATEST\" or a valid integer "
264 : "less than ",
265 : n_steps,
266 : ", received ",
267 : s_timestep);
268 147 : }
269 201 : }
270 : else
271 : // Interpolate between times rather than using values from a set timestep
272 190 : _interpolate_times = true;
273 :
274 : // Check that the number of time steps is valid
275 391 : int num_exo_times = _exodus_times->size();
276 391 : if (num_exo_times == 0)
277 0 : mooseError("In SolutionUserObjectBase, exodus file contains no timesteps.");
278 :
279 : // Account for parallel mesh
280 395 : if (_fe_problem.mesh().isDistributedMesh() && (_file_type != "exodusII") &&
281 4 : !_force_replicated_source)
282 : {
283 2 : _mesh->allow_renumbering(true);
284 2 : _mesh->prepare_for_use();
285 : }
286 : else
287 : {
288 389 : _mesh->allow_renumbering(false);
289 389 : _mesh->prepare_for_use();
290 : }
291 :
292 : // Create EquationSystems object for solution
293 391 : _es = std::make_unique<EquationSystems>(*_mesh);
294 391 : _es->add_system<libMesh::ExplicitSystem>(_system_name);
295 391 : _system = &_es->get_system(_system_name);
296 :
297 : // Get the variable name lists as set; these need to be sets to perform set_intersection
298 391 : const std::vector<std::string> & all_nodal((_file_type == "exodusII")
299 391 : ? _exodusII_io->get_nodal_var_names()
300 24 : : _nemesis_io->get_nodal_var_names());
301 391 : const std::vector<std::string> & all_elemental((_file_type == "exodusII")
302 391 : ? _exodusII_io->get_elem_var_names()
303 24 : : _nemesis_io->get_elem_var_names());
304 391 : const std::vector<std::string> & all_scalar((_file_type == "exodusII")
305 391 : ? _exodusII_io->get_global_var_names()
306 24 : : _nemesis_io->get_global_var_names());
307 :
308 : // Build nodal/elemental variable lists, limit to variables listed in 'system_variables', if
309 : // provided
310 : // This function could be called more than once, so clear the member variables so we don't keep
311 : // adding to the vectors
312 391 : _nodal_variables.clear();
313 391 : _elemental_variables.clear();
314 391 : _scalar_variables.clear();
315 391 : if (!_system_variables.empty())
316 : {
317 842 : for (const auto & var_name : _system_variables)
318 : {
319 499 : if (std::find(all_nodal.begin(), all_nodal.end(), var_name) != all_nodal.end())
320 373 : _nodal_variables.push_back(var_name);
321 126 : else if (std::find(all_elemental.begin(), all_elemental.end(), var_name) !=
322 252 : all_elemental.end())
323 126 : _elemental_variables.push_back(var_name);
324 0 : else if (std::find(all_scalar.begin(), all_scalar.end(), var_name) != all_scalar.end())
325 : // We checked other variables types first, so if a postprocessor has the same name as a
326 : // field variable, it won't get loaded as a scalar variable
327 0 : _scalar_variables.push_back(var_name);
328 : else
329 0 : paramError("system_variables", "Variable '" + var_name + "' was not found in Exodus file");
330 : }
331 : }
332 : else
333 : {
334 48 : _nodal_variables = all_nodal;
335 48 : _elemental_variables = all_elemental;
336 :
337 84 : for (auto var_name : all_scalar)
338 : // Check if the scalar matches any field variables, and ignore the var if it does. This means
339 : // its a Postprocessor.
340 108 : if (std::find(begin(_nodal_variables), end(_nodal_variables), var_name) ==
341 96 : _nodal_variables.end() &&
342 72 : std::find(begin(_elemental_variables), end(_elemental_variables), var_name) ==
343 60 : _elemental_variables.end())
344 36 : _scalar_variables.push_back(var_name);
345 : }
346 :
347 : // Add the variables to the system
348 836 : for (const auto & var_name : _nodal_variables)
349 445 : _system->add_variable(var_name, Utility::string_to_enum<Order>(_nodal_variable_order));
350 :
351 529 : for (const auto & var_name : _elemental_variables)
352 138 : _system->add_variable(var_name, CONSTANT, MONOMIAL);
353 :
354 415 : for (const auto & var_name : _scalar_variables)
355 24 : _system->add_variable(var_name, FIRST, SCALAR);
356 :
357 : // Helper for copying from the ExodusII/Nemesis file
358 581 : const auto copy_solutions = [this](auto & io, libMesh::System & system, unsigned int index)
359 : {
360 1252 : for (const auto & var_name : _nodal_variables)
361 671 : io.copy_nodal_solution(system, var_name, var_name, index);
362 :
363 767 : for (const auto & var_name : _elemental_variables)
364 186 : io.copy_elemental_solution(system, var_name, var_name, index);
365 :
366 581 : if (_scalar_variables.size() > 0)
367 48 : io.copy_scalar_solution(system, _scalar_variables, _scalar_variables, index);
368 972 : };
369 :
370 : // Initialize the equations systems
371 391 : _es->init();
372 :
373 : // Interpolate times
374 391 : if (_interpolate_times)
375 : {
376 : // Create a second equation system
377 190 : _es2 = std::make_unique<EquationSystems>(*_mesh);
378 190 : _es2->add_system<libMesh::ExplicitSystem>(_system_name);
379 190 : _system2 = &_es2->get_system(_system_name);
380 :
381 : // Add the variables to the system
382 416 : for (const auto & var_name : _nodal_variables)
383 226 : _system2->add_variable(var_name, Utility::string_to_enum<Order>(_nodal_variable_order));
384 :
385 238 : for (const auto & var_name : _elemental_variables)
386 48 : _system2->add_variable(var_name, CONSTANT, MONOMIAL);
387 :
388 214 : for (const auto & var_name : _scalar_variables)
389 24 : _system2->add_variable(var_name, FIRST, SCALAR);
390 :
391 : // Initialize
392 190 : _es2->init();
393 :
394 : // Update the times for interpolation (initially start at 0)
395 190 : updateInterpolationBracketingTimeIndices();
396 :
397 : // Copy the solutions from the first system
398 190 : if (_file_type == "exodusII")
399 : {
400 166 : copy_solutions(*_exodusII_io, *_system, _exodus_index1 + 1);
401 166 : copy_solutions(*_exodusII_io, *_system2, _exodus_index2 + 1);
402 : }
403 : else
404 : {
405 24 : copy_solutions(*_nemesis_io, *_system, _exodus_index1 + 1);
406 24 : copy_solutions(*_nemesis_io, *_system2, _exodus_index2 + 1);
407 : }
408 :
409 : // Update the systems
410 190 : _system->update();
411 190 : _es->update();
412 190 : _system2->update();
413 190 : _es2->update();
414 : }
415 :
416 : // Non-interpolated times
417 : else
418 : {
419 201 : if (_exodus_time_index > num_exo_times)
420 0 : mooseError("In SolutionUserObjectBase, timestep = ",
421 0 : _exodus_time_index,
422 : ", but there are only ",
423 : num_exo_times,
424 : " time steps.");
425 :
426 : // Copy the values from the ExodusII/Nemsis file
427 201 : if (_file_type == "exodusII")
428 201 : copy_solutions(*_exodusII_io, *_system, _exodus_time_index);
429 : else
430 0 : copy_solutions(*_nemesis_io, *_system, _exodus_time_index);
431 :
432 : // Update the equations systems
433 201 : _system->update();
434 201 : _es->update();
435 : }
436 391 : }
437 :
438 : Real
439 40572 : SolutionUserObjectBase::directValue(const Node * node, const std::string & var_name) const
440 : {
441 : // Get the libmesh variable and system numbers
442 40572 : unsigned int var_num = _system->variable_number(var_name);
443 40572 : unsigned int sys_num = _system->number();
444 :
445 : // Get the node id and associated dof
446 40572 : dof_id_type node_id = node->id();
447 40572 : const Node & sys_node = _system->get_mesh().node_ref(node_id);
448 : mooseAssert(sys_node.n_dofs(sys_num, var_num) > 0,
449 : "Variable " << var_name << " has no DoFs on node " << sys_node.id());
450 40572 : dof_id_type dof_id = sys_node.dof_number(sys_num, var_num, 0);
451 :
452 : // Return the desired value for the dof
453 40572 : return directValue(dof_id);
454 : }
455 :
456 : Real
457 0 : SolutionUserObjectBase::directValue(const Elem * elem, const std::string & var_name) const
458 : {
459 : // Get the libmesh variable and system numbers
460 0 : unsigned int var_num = _system->variable_number(var_name);
461 0 : unsigned int sys_num = _system->number();
462 :
463 : // Get the element id and associated dof
464 0 : dof_id_type elem_id = elem->id();
465 0 : const Elem & sys_elem = _system->get_mesh().elem_ref(elem_id);
466 : mooseAssert(sys_elem.n_dofs(sys_num, var_num) > 0,
467 : "Variable " << var_name << " has no DoFs on element " << sys_elem.id());
468 0 : dof_id_type dof_id = sys_elem.dof_number(sys_num, var_num, 0);
469 :
470 : // Return the desired value
471 0 : return directValue(dof_id);
472 : }
473 :
474 : void
475 1688 : SolutionUserObjectBase::initialize()
476 : {
477 1688 : }
478 :
479 : void
480 1688 : SolutionUserObjectBase::finalize()
481 : {
482 1688 : }
483 :
484 : void
485 1743 : SolutionUserObjectBase::timestepSetup()
486 : {
487 : // Update time interpolation for ExodusII or Nemesis solution
488 1743 : if (_interpolate_times)
489 718 : updateTimeInterpolationFromFile();
490 :
491 : // Clear the caches
492 1743 : _cached_p(0) = std::numeric_limits<Real>::max();
493 1743 : _cached_p2(0) = std::numeric_limits<Real>::max();
494 1743 : }
495 :
496 : void
497 1688 : SolutionUserObjectBase::execute()
498 : {
499 1688 : }
500 :
501 : void
502 501 : SolutionUserObjectBase::initialSetup()
503 : {
504 :
505 : // Make sure this only happens once
506 501 : if (_initialized)
507 0 : return;
508 :
509 : // We need to detect exodus files early
510 : const bool has_exodus_extension =
511 1637 : MooseUtils::hasExtension(_mesh_file, "e", /*strip_exodus_ext =*/true) ||
512 769 : MooseUtils::hasExtension(_mesh_file, "exo", true);
513 : const bool has_nemesis_extension =
514 1980 : MooseUtils::hasExtension(_mesh_file, "n", /*strip_exodus_ext =*/true) ||
515 1455 : MooseUtils::hasExtension(_mesh_file, "nem", true);
516 501 : if (has_exodus_extension)
517 734 : _file_type = "exodusII";
518 134 : else if (has_nemesis_extension)
519 48 : _file_type = "nemesis";
520 :
521 : // Create a libmesh::Mesh object for storing the loaded data.
522 507 : if (!_fe_problem.mesh().isDistributedMesh() || (_file_type == "exodusII") ||
523 6 : _force_replicated_source)
524 497 : _mesh = std::make_unique<ReplicatedMesh>(_communicator);
525 : else
526 4 : _mesh = std::make_unique<DistributedMesh>(_communicator);
527 :
528 : // ExodusII or Nemesis mesh file supplied
529 501 : if (has_exodus_extension || has_nemesis_extension)
530 391 : readExodusIIOrNemesis();
531 :
532 : // XDA mesh file supplied
533 220 : else if (MooseUtils::hasExtension(_mesh_file, "xda"))
534 : {
535 110 : _file_type = "xda";
536 110 : readXda();
537 : }
538 :
539 : // XDR mesh file supplied
540 0 : else if (MooseUtils::hasExtension(_mesh_file, "xdr"))
541 : {
542 0 : _file_type = "xdr";
543 0 : readXda();
544 : }
545 :
546 : // Produce an error for an unknown file type
547 : else
548 0 : mooseError(
549 : "In SolutionUserObjectBase, invalid file type: only .xda, .xdr, .exo and .e supported");
550 :
551 : // Initialize the serial solution vector
552 501 : _serialized_solution = NumericVector<Number>::build(_communicator);
553 501 : _serialized_solution->init(_system->n_dofs(), false, libMesh::SERIAL);
554 :
555 : // Pull down a full copy of this vector on every processor so we can get values in parallel
556 501 : _system->solution->localize(*_serialized_solution);
557 :
558 : // Vector of variable numbers to apply the MeshFunction to
559 501 : std::vector<unsigned int> var_nums;
560 :
561 : // If no variables were given, use all of them
562 501 : if (_system_variables.empty())
563 : {
564 48 : _system->get_all_variable_numbers(var_nums);
565 156 : for (const auto & var_num : var_nums)
566 108 : _system_variables.push_back(_system->variable_name(var_num));
567 : }
568 :
569 : // Otherwise, gather the numbers for the variables given
570 : else
571 : {
572 1062 : for (const auto & var_name : _system_variables)
573 609 : var_nums.push_back(_system->variable_number(var_name));
574 : }
575 :
576 : // Create the MeshFunction for working with the solution data
577 1002 : _mesh_function = std::make_unique<libMesh::MeshFunction>(
578 1002 : *_es, *_serialized_solution, _system->get_dof_map(), var_nums);
579 501 : _mesh_function->init();
580 :
581 : // Tell the MeshFunctions that we might be querying them outside the
582 : // mesh, so we can handle any errors at the MOOSE rather than at the
583 : // libMesh level.
584 501 : DenseVector<Number> default_values;
585 501 : _mesh_function->enable_out_of_mesh_mode(default_values);
586 :
587 : // Build second MeshFunction for interpolation
588 501 : if (_interpolate_times)
589 : {
590 : // Need to pull down a full copy of this vector on every processor so we can get values in
591 : // parallel
592 190 : _serialized_solution2 = NumericVector<Number>::build(_communicator);
593 190 : _serialized_solution2->init(_system2->n_dofs(), false, libMesh::SERIAL);
594 190 : _system2->solution->localize(*_serialized_solution2);
595 :
596 : // Create the MeshFunction for the second copy of the data
597 380 : _mesh_function2 = std::make_unique<libMesh::MeshFunction>(
598 380 : *_es2, *_serialized_solution2, _system2->get_dof_map(), var_nums);
599 190 : _mesh_function2->init();
600 190 : _mesh_function2->enable_out_of_mesh_mode(default_values);
601 : }
602 :
603 : // Populate the MeshFunction variable index
604 1218 : for (unsigned int i = 0; i < _system_variables.size(); ++i)
605 : {
606 717 : std::string name = _system_variables[i];
607 717 : _local_variable_index[name] = i;
608 717 : }
609 :
610 : // If the start time is not the same as in the exodus file, we may need this on INITIAL
611 501 : if (_interpolate_times)
612 190 : updateTimeInterpolationFromFile();
613 :
614 : // Set initialization flag
615 501 : _initialized = true;
616 501 : }
617 :
618 : MooseEnum
619 195 : SolutionUserObjectBase::getSolutionFileType() const
620 : {
621 195 : return _file_type;
622 : }
623 :
624 : void
625 908 : SolutionUserObjectBase::updateTimeInterpolationFromFile()
626 : {
627 908 : if (_t != _interpolation_time)
628 : {
629 689 : if (updateInterpolationBracketingTimeIndices())
630 : {
631 : // Helper to copy/update solutions from a file IO to a system
632 458 : const auto update_solutions = [this](auto & io, libMesh::System & system, unsigned int index)
633 : {
634 938 : for (const auto & var_name : _nodal_variables)
635 480 : io.copy_nodal_solution(system, var_name, var_name, index);
636 458 : for (const auto & var_name : _elemental_variables)
637 0 : io.copy_elemental_solution(system, var_name, var_name, index);
638 458 : if (_scalar_variables.size() > 0)
639 44 : io.copy_scalar_solution(system, _scalar_variables, _scalar_variables, index);
640 687 : };
641 :
642 229 : if (_file_type == "exodusII")
643 229 : update_solutions(*_exodusII_io, *_system, _exodus_index1 + 1);
644 : else
645 0 : update_solutions(*_nemesis_io, *_system, _exodus_index1 + 1);
646 :
647 229 : _system->update();
648 229 : _es->update();
649 229 : _system->solution->localize(*_serialized_solution);
650 :
651 229 : if (_file_type == "exodusII")
652 229 : update_solutions(*_exodusII_io, *_system2, _exodus_index2 + 1);
653 : else
654 0 : update_solutions(*_nemesis_io, *_system2, _exodus_index2 + 1);
655 :
656 229 : _system2->update();
657 229 : _es2->update();
658 229 : _system2->solution->localize(*_serialized_solution2);
659 : }
660 689 : _interpolation_time = _t;
661 : }
662 908 : }
663 :
664 : bool
665 879 : SolutionUserObjectBase::updateInterpolationBracketingTimeIndices()
666 : {
667 879 : if (_file_type != "exodusII" && _file_type != "nemesis")
668 0 : mooseError("getTimeInterpolationData only applicable for exodusII or Nemesis file type");
669 :
670 879 : int old_index1 = _exodus_index1;
671 879 : int old_index2 = _exodus_index2;
672 :
673 879 : int num_exo_times = _exodus_times->size();
674 :
675 879 : const auto solution_time = solutionSampleTime();
676 :
677 879 : if (solution_time < (*_exodus_times)[0])
678 : {
679 0 : _exodus_index1 = 0;
680 0 : _exodus_index2 = 0;
681 0 : _interpolation_factor = 0.0;
682 : }
683 : else
684 : {
685 1929 : for (int i = 0; i < num_exo_times - 1; ++i)
686 : {
687 1825 : if (solution_time <= (*_exodus_times)[i + 1])
688 : {
689 775 : _exodus_index1 = i;
690 775 : _exodus_index2 = i + 1;
691 775 : _interpolation_factor =
692 775 : (solution_time - (*_exodus_times)[i]) / ((*_exodus_times)[i + 1] - (*_exodus_times)[i]);
693 775 : break;
694 : }
695 1050 : else if (i == num_exo_times - 2)
696 : {
697 0 : _exodus_index1 = num_exo_times - 1;
698 0 : _exodus_index2 = num_exo_times - 1;
699 0 : _interpolation_factor = 1.0;
700 0 : break;
701 : }
702 : }
703 : }
704 :
705 879 : bool indices_modified(false);
706 :
707 879 : if (_exodus_index1 != old_index1 || _exodus_index2 != old_index2)
708 383 : indices_modified = true;
709 :
710 879 : return indices_modified;
711 : }
712 :
713 : unsigned int
714 581841 : SolutionUserObjectBase::getLocalVarIndex(const std::string & var_name) const
715 : {
716 : // Extract the variable index for the MeshFunction(s)
717 581841 : std::map<std::string, unsigned int>::const_iterator it = _local_variable_index.find(var_name);
718 581841 : if (it == _local_variable_index.end())
719 6 : mooseError("Value requested for nonexistent variable '",
720 : var_name,
721 : "' in the '",
722 3 : name(),
723 : "' SolutionUserObjectBase.\nSystem selected: ",
724 3 : _system_name,
725 : "\nAvailable variables:\n",
726 9 : ConsoleUtils::formatString(Moose::stringify(_system_variables), ""));
727 1163676 : return it->second;
728 : }
729 :
730 : Real
731 0 : SolutionUserObjectBase::pointValueWrapper(Real t,
732 : const Point & p,
733 : const std::string & var_name,
734 : const MooseEnum & weighting_type,
735 : const std::set<subdomain_id_type> * subdomain_ids) const
736 : {
737 : // first check if the FE type is continuous because in that case the value is
738 : // unique and we can take a short cut, the default weighting_type found_first also
739 : // shortcuts out
740 : auto family =
741 0 : _fe_problem
742 0 : .getVariable(
743 0 : _tid, var_name, Moose::VarKindType::VAR_ANY, Moose::VarFieldType::VAR_FIELD_STANDARD)
744 0 : .feType()
745 0 : .family;
746 :
747 0 : if (weighting_type == 1 ||
748 0 : (family != L2_LAGRANGE && family != MONOMIAL && family != L2_HIERARCHIC))
749 0 : return pointValue(t, p, var_name, subdomain_ids);
750 :
751 : // the shape function is discontinuous so we need to compute a suitable unique value
752 0 : std::map<const Elem *, Real> values = discontinuousPointValue(t, p, var_name);
753 0 : switch (weighting_type)
754 : {
755 0 : case 2:
756 : {
757 0 : Real average = 0.0;
758 0 : for (auto & v : values)
759 0 : average += v.second;
760 0 : return average / Real(values.size());
761 : }
762 0 : case 4:
763 : {
764 0 : Real smallest_elem_id_value = std::numeric_limits<Real>::max();
765 0 : dof_id_type smallest_elem_id = _fe_problem.mesh().maxElemId();
766 0 : for (auto & v : values)
767 0 : if (v.first->id() < smallest_elem_id)
768 : {
769 0 : smallest_elem_id = v.first->id();
770 0 : smallest_elem_id_value = v.second;
771 : }
772 0 : return smallest_elem_id_value;
773 : }
774 0 : case 8:
775 : {
776 0 : Real largest_elem_id_value = std::numeric_limits<Real>::lowest();
777 0 : dof_id_type largest_elem_id = 0;
778 0 : for (auto & v : values)
779 0 : if (v.first->id() > largest_elem_id)
780 : {
781 0 : largest_elem_id = v.first->id();
782 0 : largest_elem_id_value = v.second;
783 : }
784 0 : return largest_elem_id_value;
785 : }
786 : }
787 :
788 0 : mooseError("SolutionUserObjectBase::pointValueWrapper reaches line that it should not be able to "
789 : "reach.");
790 : return 0.0;
791 0 : }
792 :
793 : Real
794 581460 : SolutionUserObjectBase::pointValue(Real t,
795 : const Point & p,
796 : const std::string & var_name,
797 : const std::set<subdomain_id_type> * subdomain_ids) const
798 : {
799 581460 : const unsigned int local_var_index = getLocalVarIndex(var_name);
800 581460 : return pointValue(t, p, local_var_index, subdomain_ids);
801 : }
802 :
803 : Real
804 797128 : SolutionUserObjectBase::pointValue(Real libmesh_dbg_var(t),
805 : const Point & p,
806 : const unsigned int local_var_index,
807 : const std::set<subdomain_id_type> * subdomain_ids) const
808 : {
809 : // Create copy of point
810 797128 : Point pt(p);
811 :
812 : // do the transformations
813 2387352 : for (unsigned int trans_num = 0; trans_num < _transformation_order.size(); ++trans_num)
814 : {
815 1590224 : if (_transformation_order[trans_num] == "rotation0")
816 5472 : pt = _r0 * pt;
817 1584752 : else if (_transformation_order[trans_num] == "translation")
818 3165472 : for (const auto i : make_range(Moose::dim))
819 2374104 : pt(i) -= _translation[i];
820 793384 : else if (_transformation_order[trans_num] == "scale")
821 3165472 : for (const auto i : make_range(Moose::dim))
822 2374104 : pt(i) /= _scale[i];
823 2016 : else if (_transformation_order[trans_num] == "scale_multiplier")
824 1152 : for (const auto i : make_range(Moose::dim))
825 864 : pt(i) *= _scale_multiplier[i];
826 1728 : else if (_transformation_order[trans_num] == "rotation1")
827 1728 : pt = _r1 * pt;
828 : }
829 :
830 : // Extract the value at the current point
831 797128 : Real val = evalMeshFunction(pt, local_var_index, 1, subdomain_ids);
832 :
833 : // Interpolate
834 797126 : if (_interpolate_times && (_file_type == "exodusII" || _file_type == "nemesis"))
835 : {
836 : mooseAssert(t == _interpolation_time,
837 : "Time passed into value() must match time at last call to timestepSetup()");
838 180210 : Real val2 = evalMeshFunction(pt, local_var_index, 2, subdomain_ids);
839 180210 : val = val + (val2 - val) * _interpolation_factor;
840 : }
841 :
842 797126 : return val;
843 : }
844 :
845 : std::map<const Elem *, Real>
846 99 : SolutionUserObjectBase::discontinuousPointValue(
847 : Real t,
848 : const Point & p,
849 : const std::string & var_name,
850 : const std::set<subdomain_id_type> * subdomain_ids) const
851 : {
852 99 : const unsigned int local_var_index = getLocalVarIndex(var_name);
853 99 : return discontinuousPointValue(t, p, local_var_index, subdomain_ids);
854 : }
855 :
856 : std::map<const Elem *, Real>
857 99 : SolutionUserObjectBase::discontinuousPointValue(
858 : Real libmesh_dbg_var(t),
859 : Point pt,
860 : const unsigned int local_var_index,
861 : const std::set<subdomain_id_type> * subdomain_ids) const
862 : {
863 : // do the transformations
864 297 : for (unsigned int trans_num = 0; trans_num < _transformation_order.size(); ++trans_num)
865 : {
866 198 : if (_transformation_order[trans_num] == "rotation0")
867 0 : pt = _r0 * pt;
868 198 : else if (_transformation_order[trans_num] == "translation")
869 396 : for (const auto i : make_range(Moose::dim))
870 297 : pt(i) -= _translation[i];
871 99 : else if (_transformation_order[trans_num] == "scale")
872 396 : for (const auto i : make_range(Moose::dim))
873 297 : pt(i) /= _scale[i];
874 0 : else if (_transformation_order[trans_num] == "scale_multiplier")
875 0 : for (const auto i : make_range(Moose::dim))
876 0 : pt(i) *= _scale_multiplier[i];
877 0 : else if (_transformation_order[trans_num] == "rotation1")
878 0 : pt = _r1 * pt;
879 : }
880 :
881 : // Extract the value at the current point
882 : std::map<const Elem *, Real> map =
883 99 : evalMultiValuedMeshFunction(pt, local_var_index, 1, subdomain_ids);
884 :
885 : // Interpolate
886 99 : if (_interpolate_times && (_file_type == "exodusII" || _file_type == "nemesis"))
887 : {
888 : mooseAssert(t == _interpolation_time,
889 : "Time passed into value() must match time at last call to timestepSetup()");
890 99 : std::map<const Elem *, Real> map2 = evalMultiValuedMeshFunction(pt, local_var_index, 2);
891 :
892 99 : if (map.size() != map2.size())
893 0 : mooseError(
894 : "In SolutionUserObjectBase::discontinuousPointValue map and map2 have different size");
895 :
896 : // construct the interpolated map
897 231 : for (auto & k : map)
898 : {
899 132 : if (map2.find(k.first) == map2.end())
900 0 : mooseError(
901 : "In SolutionUserObjectBase::discontinuousPointValue map and map2 have differing keys");
902 132 : Real val = k.second;
903 132 : Real val2 = map2[k.first];
904 132 : map[k.first] = val + (val2 - val) * _interpolation_factor;
905 : }
906 99 : }
907 :
908 99 : return map;
909 0 : }
910 :
911 : RealGradient
912 0 : SolutionUserObjectBase::pointValueGradientWrapper(
913 : Real t,
914 : const Point & p,
915 : const std::string & var_name,
916 : const MooseEnum & weighting_type,
917 : const std::set<subdomain_id_type> * subdomain_ids) const
918 : {
919 : // the default weighting_type found_first shortcuts out
920 0 : if (weighting_type == 1)
921 0 : return pointValueGradient(t, p, var_name, subdomain_ids);
922 :
923 : // the shape function is discontinuous so we need to compute a suitable unique value
924 : std::map<const Elem *, RealGradient> values =
925 0 : discontinuousPointValueGradient(t, p, var_name, subdomain_ids);
926 0 : switch (weighting_type)
927 : {
928 0 : case 2:
929 : {
930 0 : RealGradient average = RealGradient(0.0, 0.0, 0.0);
931 0 : for (auto & v : values)
932 0 : average += v.second;
933 0 : return average / Real(values.size());
934 : }
935 0 : case 4:
936 : {
937 0 : RealGradient smallest_elem_id_value;
938 0 : dof_id_type smallest_elem_id = _fe_problem.mesh().maxElemId();
939 0 : for (auto & v : values)
940 0 : if (v.first->id() < smallest_elem_id)
941 : {
942 0 : smallest_elem_id = v.first->id();
943 0 : smallest_elem_id_value = v.second;
944 : }
945 0 : return smallest_elem_id_value;
946 : }
947 0 : case 8:
948 : {
949 0 : RealGradient largest_elem_id_value;
950 0 : dof_id_type largest_elem_id = 0;
951 0 : for (auto & v : values)
952 0 : if (v.first->id() > largest_elem_id)
953 : {
954 0 : largest_elem_id = v.first->id();
955 0 : largest_elem_id_value = v.second;
956 : }
957 0 : return largest_elem_id_value;
958 : }
959 : }
960 :
961 0 : mooseError("SolutionUserObjectBase::pointValueGradientWrapper reaches line that it should not be "
962 : "able to reach.");
963 : return RealGradient(0.0, 0.0, 0.0);
964 0 : }
965 :
966 : RealGradient
967 0 : SolutionUserObjectBase::pointValueGradient(Real t,
968 : const Point & p,
969 : const std::string & var_name,
970 : const std::set<subdomain_id_type> * subdomain_ids) const
971 : {
972 0 : const unsigned int local_var_index = getLocalVarIndex(var_name);
973 0 : return pointValueGradient(t, p, local_var_index, subdomain_ids);
974 : }
975 :
976 : RealGradient
977 56628 : SolutionUserObjectBase::pointValueGradient(Real libmesh_dbg_var(t),
978 : Point pt,
979 : const unsigned int local_var_index,
980 : const std::set<subdomain_id_type> * subdomain_ids) const
981 : {
982 : // do the transformations
983 169884 : for (unsigned int trans_num = 0; trans_num < _transformation_order.size(); ++trans_num)
984 : {
985 113256 : if (_transformation_order[trans_num] == "rotation0")
986 0 : pt = _r0 * pt;
987 113256 : else if (_transformation_order[trans_num] == "translation")
988 226512 : for (const auto i : make_range(Moose::dim))
989 169884 : pt(i) -= _translation[i];
990 56628 : else if (_transformation_order[trans_num] == "scale")
991 226512 : for (const auto i : make_range(Moose::dim))
992 169884 : pt(i) /= _scale[i];
993 0 : else if (_transformation_order[trans_num] == "scale_multiplier")
994 0 : for (const auto i : make_range(Moose::dim))
995 0 : pt(i) *= _scale_multiplier[i];
996 0 : else if (_transformation_order[trans_num] == "rotation1")
997 0 : pt = _r1 * pt;
998 : }
999 :
1000 : // Extract the value at the current point
1001 56628 : RealGradient val = evalMeshFunctionGradient(pt, local_var_index, 1, subdomain_ids);
1002 :
1003 : // Interpolate
1004 56628 : if (_interpolate_times && (_file_type == "exodusII" || _file_type == "nemesis"))
1005 : {
1006 : mooseAssert(t == _interpolation_time,
1007 : "Time passed into value() must match time at last call to timestepSetup()");
1008 56628 : RealGradient val2 = evalMeshFunctionGradient(pt, local_var_index, 2, subdomain_ids);
1009 56628 : val = val + (val2 - val) * _interpolation_factor;
1010 : }
1011 :
1012 56628 : return val;
1013 : }
1014 :
1015 : std::map<const Elem *, RealGradient>
1016 99 : SolutionUserObjectBase::discontinuousPointValueGradient(
1017 : Real t,
1018 : const Point & p,
1019 : const std::string & var_name,
1020 : const std::set<subdomain_id_type> * subdomain_ids) const
1021 : {
1022 99 : const unsigned int local_var_index = getLocalVarIndex(var_name);
1023 99 : return discontinuousPointValueGradient(t, p, local_var_index, subdomain_ids);
1024 : }
1025 :
1026 : std::map<const Elem *, RealGradient>
1027 99 : SolutionUserObjectBase::discontinuousPointValueGradient(
1028 : Real libmesh_dbg_var(t),
1029 : Point pt,
1030 : const unsigned int local_var_index,
1031 : const std::set<subdomain_id_type> * subdomain_ids) const
1032 : {
1033 : // do the transformations
1034 297 : for (unsigned int trans_num = 0; trans_num < _transformation_order.size(); ++trans_num)
1035 : {
1036 198 : if (_transformation_order[trans_num] == "rotation0")
1037 0 : pt = _r0 * pt;
1038 198 : else if (_transformation_order[trans_num] == "translation")
1039 396 : for (const auto i : make_range(Moose::dim))
1040 297 : pt(i) -= _translation[i];
1041 99 : else if (_transformation_order[trans_num] == "scale")
1042 396 : for (const auto i : make_range(Moose::dim))
1043 297 : pt(i) /= _scale[i];
1044 0 : else if (_transformation_order[trans_num] == "scale_multiplier")
1045 0 : for (const auto i : make_range(Moose::dim))
1046 0 : pt(i) *= _scale_multiplier[i];
1047 0 : else if (_transformation_order[trans_num] == "rotation1")
1048 0 : pt = _r1 * pt;
1049 : }
1050 :
1051 : // Extract the value at the current point
1052 : std::map<const Elem *, RealGradient> map =
1053 99 : evalMultiValuedMeshFunctionGradient(pt, local_var_index, 1, subdomain_ids);
1054 :
1055 : // Interpolate
1056 99 : if (_interpolate_times && (_file_type == "exodusII" || _file_type == "nemesis"))
1057 : {
1058 : mooseAssert(t == _interpolation_time,
1059 : "Time passed into value() must match time at last call to timestepSetup()");
1060 : std::map<const Elem *, RealGradient> map2 =
1061 99 : evalMultiValuedMeshFunctionGradient(pt, local_var_index, 1, subdomain_ids);
1062 :
1063 99 : if (map.size() != map2.size())
1064 0 : mooseError(
1065 : "In SolutionUserObjectBase::discontinuousPointValue map and map2 have different size");
1066 :
1067 : // construct the interpolated map
1068 231 : for (auto & k : map)
1069 : {
1070 132 : if (map2.find(k.first) == map2.end())
1071 0 : mooseError(
1072 : "In SolutionUserObjectBase::discontinuousPointValue map and map2 have differing keys");
1073 132 : RealGradient val = k.second;
1074 132 : RealGradient val2 = map2[k.first];
1075 132 : map[k.first] = val + (val2 - val) * _interpolation_factor;
1076 : }
1077 99 : }
1078 :
1079 99 : return map;
1080 0 : }
1081 :
1082 : Real
1083 40633 : SolutionUserObjectBase::directValue(dof_id_type dof_index) const
1084 : {
1085 40633 : Real val = (*_serialized_solution)(dof_index);
1086 40633 : if (_interpolate_times && (_file_type == "exodusII" || _file_type == "nemesis"))
1087 : {
1088 20311 : Real val2 = (*_serialized_solution2)(dof_index);
1089 20311 : val = val + (val2 - val) * _interpolation_factor;
1090 : }
1091 40633 : return val;
1092 : }
1093 :
1094 : Real
1095 977338 : SolutionUserObjectBase::evalMeshFunction(const Point & p,
1096 : const unsigned int local_var_index,
1097 : unsigned int func_num,
1098 : const std::set<subdomain_id_type> * subdomain_ids) const
1099 : {
1100 : // Storage for mesh function output
1101 977338 : DenseVector<Number> output;
1102 :
1103 : // Extract a value from the _mesh_function
1104 : {
1105 977338 : Threads::spin_mutex::scoped_lock lock(_solution_user_object_mutex);
1106 :
1107 977338 : if (func_num == 1)
1108 : {
1109 : // Check the cache
1110 797128 : if (p == _cached_p && (!subdomain_ids || (*subdomain_ids == _cached_subdomain_ids)))
1111 341448 : return _cached_values(local_var_index);
1112 :
1113 : // else get a new value
1114 455680 : (*_mesh_function)(p, 0.0, output, subdomain_ids);
1115 : // and cache it
1116 455680 : _cached_p = p;
1117 455680 : if (subdomain_ids)
1118 686 : _cached_subdomain_ids = *subdomain_ids;
1119 455680 : _cached_values = output;
1120 : }
1121 :
1122 : // Extract a value from _mesh_function2
1123 180210 : else if (func_num == 2)
1124 : {
1125 : // Check the cache
1126 180210 : if (p == _cached_p2 && (!subdomain_ids || (*subdomain_ids == _cached_subdomain_ids2)))
1127 7521 : return _cached_values2(local_var_index);
1128 :
1129 : // else get a new value
1130 172689 : (*_mesh_function2)(p, 0.0, output, subdomain_ids);
1131 : // and cache it
1132 172689 : _cached_p2 = p;
1133 172689 : if (subdomain_ids)
1134 0 : _cached_subdomain_ids2 = *subdomain_ids;
1135 172689 : _cached_values2 = output;
1136 : }
1137 : else
1138 0 : mooseError("The func_num must be 1 or 2");
1139 977338 : }
1140 :
1141 : // Error if the data is out-of-range, which will be the case if the mesh functions are evaluated
1142 : // outside the domain
1143 628369 : if (output.size() == 0)
1144 : {
1145 2 : std::ostringstream oss;
1146 2 : p.print(oss);
1147 2 : mooseError("Failed to access the data for variable '",
1148 2 : _system_variables[local_var_index],
1149 : "' at point ",
1150 2 : oss.str(),
1151 : " in the '",
1152 2 : name(),
1153 : "' SolutionUserObjectBase");
1154 0 : }
1155 628367 : return output(local_var_index);
1156 977336 : }
1157 :
1158 : std::map<const Elem *, Real>
1159 198 : SolutionUserObjectBase::evalMultiValuedMeshFunction(
1160 : const Point & p,
1161 : const unsigned int local_var_index,
1162 : unsigned int func_num,
1163 : const std::set<subdomain_id_type> * subdomain_ids) const
1164 : {
1165 : // Storage for mesh function output
1166 198 : std::map<const Elem *, DenseVector<Number>> temporary_output;
1167 :
1168 : // Extract a value from the _mesh_function
1169 : {
1170 198 : Threads::spin_mutex::scoped_lock lock(_solution_user_object_mutex);
1171 198 : if (func_num == 1)
1172 99 : _mesh_function->discontinuous_value(p, 0.0, temporary_output, subdomain_ids);
1173 :
1174 : // Extract a value from _mesh_function2
1175 99 : else if (func_num == 2)
1176 99 : _mesh_function2->discontinuous_value(p, 0.0, temporary_output, subdomain_ids);
1177 :
1178 : else
1179 0 : mooseError("The func_num must be 1 or 2");
1180 198 : }
1181 :
1182 : // Error if the data is out-of-range, which will be the case if the mesh functions are evaluated
1183 : // outside the domain
1184 198 : if (temporary_output.size() == 0)
1185 : {
1186 0 : std::ostringstream oss;
1187 0 : p.print(oss);
1188 0 : mooseError("Failed to access the data for variable '",
1189 0 : _system_variables[local_var_index],
1190 : "' at point ",
1191 0 : oss.str(),
1192 : " in the '",
1193 0 : name(),
1194 : "' SolutionUserObjectBase");
1195 0 : }
1196 :
1197 : // Fill the actual map that is returned
1198 198 : std::map<const Elem *, Real> output;
1199 462 : for (auto & k : temporary_output)
1200 : {
1201 : mooseAssert(
1202 : k.second.size() > local_var_index,
1203 : "In SolutionUserObjectBase::evalMultiValuedMeshFunction variable with local_var_index "
1204 : << local_var_index << " does not exist");
1205 264 : output[k.first] = k.second(local_var_index);
1206 : }
1207 :
1208 396 : return output;
1209 198 : }
1210 :
1211 : RealGradient
1212 113256 : SolutionUserObjectBase::evalMeshFunctionGradient(
1213 : const Point & p,
1214 : const unsigned int local_var_index,
1215 : unsigned int func_num,
1216 : const std::set<subdomain_id_type> * subdomain_ids) const
1217 : {
1218 : // Storage for mesh function output
1219 113256 : std::vector<Gradient> output;
1220 :
1221 : // Extract a value from the _mesh_function
1222 : {
1223 113256 : Threads::spin_mutex::scoped_lock lock(_solution_user_object_mutex);
1224 113256 : if (func_num == 1)
1225 56628 : _mesh_function->gradient(p, 0.0, output, subdomain_ids);
1226 :
1227 : // Extract a value from _mesh_function2
1228 56628 : else if (func_num == 2)
1229 56628 : _mesh_function2->gradient(p, 0.0, output, subdomain_ids);
1230 :
1231 : else
1232 0 : mooseError("The func_num must be 1 or 2");
1233 113256 : }
1234 :
1235 : // Error if the data is out-of-range, which will be the case if the mesh functions are evaluated
1236 : // outside the domain
1237 113256 : if (output.size() == 0)
1238 : {
1239 0 : std::ostringstream oss;
1240 0 : p.print(oss);
1241 0 : mooseError("Failed to access the data for variable '",
1242 0 : _system_variables[local_var_index],
1243 : "' at point ",
1244 0 : oss.str(),
1245 : " in the '",
1246 0 : name(),
1247 : "' SolutionUserObjectBase");
1248 0 : }
1249 226512 : return output[local_var_index];
1250 113256 : }
1251 :
1252 : std::map<const Elem *, RealGradient>
1253 198 : SolutionUserObjectBase::evalMultiValuedMeshFunctionGradient(
1254 : const Point & p,
1255 : const unsigned int local_var_index,
1256 : unsigned int func_num,
1257 : const std::set<subdomain_id_type> * subdomain_ids) const
1258 : {
1259 : // Storage for mesh function output
1260 198 : std::map<const Elem *, std::vector<Gradient>> temporary_output;
1261 :
1262 : // Extract a value from the _mesh_function
1263 : {
1264 198 : Threads::spin_mutex::scoped_lock lock(_solution_user_object_mutex);
1265 198 : if (func_num == 1)
1266 198 : _mesh_function->discontinuous_gradient(p, 0.0, temporary_output, subdomain_ids);
1267 :
1268 : // Extract a value from _mesh_function2
1269 0 : else if (func_num == 2)
1270 0 : _mesh_function2->discontinuous_gradient(p, 0.0, temporary_output, subdomain_ids);
1271 :
1272 : else
1273 0 : mooseError("The func_num must be 1 or 2");
1274 198 : }
1275 :
1276 : // Error if the data is out-of-range, which will be the case if the mesh functions are evaluated
1277 : // outside the domain
1278 198 : if (temporary_output.size() == 0)
1279 : {
1280 0 : std::ostringstream oss;
1281 0 : p.print(oss);
1282 0 : mooseError("Failed to access the data for variable '",
1283 0 : _system_variables[local_var_index],
1284 : "' at point ",
1285 0 : oss.str(),
1286 : " in the '",
1287 0 : name(),
1288 : "' SolutionUserObjectBase");
1289 0 : }
1290 :
1291 : // Fill the actual map that is returned
1292 198 : std::map<const Elem *, RealGradient> output;
1293 462 : for (auto & k : temporary_output)
1294 : {
1295 : mooseAssert(
1296 : k.second.size() > local_var_index,
1297 : "In SolutionUserObjectBase::evalMultiValuedMeshFunction variable with local_var_index "
1298 : << local_var_index << " does not exist");
1299 264 : output[k.first] = k.second[local_var_index];
1300 : }
1301 :
1302 396 : return output;
1303 198 : }
1304 :
1305 : const std::vector<std::string> &
1306 300 : SolutionUserObjectBase::variableNames() const
1307 : {
1308 300 : return _system_variables;
1309 : }
1310 :
1311 : bool
1312 0 : SolutionUserObjectBase::isVariableNodal(const std::string & var_name) const
1313 : {
1314 0 : return std::find(_nodal_variables.begin(), _nodal_variables.end(), var_name) !=
1315 0 : _nodal_variables.end();
1316 : }
1317 :
1318 : Real
1319 61 : SolutionUserObjectBase::scalarValue(Real /*t*/, const std::string & var_name) const
1320 : {
1321 61 : unsigned int var_num = _system->variable_number(var_name);
1322 61 : const DofMap & dof_map = _system->get_dof_map();
1323 61 : std::vector<dof_id_type> dofs;
1324 61 : dof_map.SCALAR_dof_indices(dofs, var_num);
1325 : // We can handle only FIRST order scalar variables
1326 122 : return directValue(dofs[0]);
1327 61 : }
1328 :
1329 : void
1330 391 : SolutionUserObjectBase::readBlockIdMapFromExodusII()
1331 : {
1332 : #ifdef LIBMESH_HAVE_EXODUS_API
1333 391 : const auto * exio_helper = (_file_type == "exodusII") ? &_exodusII_io->get_exio_helper()
1334 24 : : &_nemesis_io->get_nemio_helper();
1335 391 : const auto & id_to_block = exio_helper->id_to_block_names;
1336 391 : _block_name_to_id.clear();
1337 391 : _block_id_to_name.clear();
1338 818 : for (const auto & it : id_to_block)
1339 : {
1340 427 : _block_name_to_id[it.second] = it.first;
1341 427 : _block_id_to_name[it.first] = it.second;
1342 : }
1343 : #endif
1344 391 : }
|
