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