www.mooseframework.org
SolutionUserObject.C
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://www.mooseframework.org
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 "SolutionUserObject.h"
11 
12 // MOOSE includes
13 #include "MooseError.h"
14 #include "MooseMesh.h"
15 #include "MooseUtils.h"
16 #include "MooseVariableFE.h"
17 #include "RotationMatrix.h"
18 
19 // libMesh includes
20 #include "libmesh/equation_systems.h"
21 #include "libmesh/mesh_function.h"
22 #include "libmesh/numeric_vector.h"
23 #include "libmesh/nonlinear_implicit_system.h"
24 #include "libmesh/transient_system.h"
25 #include "libmesh/parallel_mesh.h"
26 #include "libmesh/serial_mesh.h"
27 #include "libmesh/exodusII_io.h"
28 #include "libmesh/enum_xdr_mode.h"
29 
31 
32 template <>
35 {
36  // Get the input parameters from the parent class
38 
39  // Add required parameters
40  params.addRequiredParam<MeshFileName>(
41  "mesh", "The name of the mesh file (must be xda or exodusII file).");
42  params.addParam<std::vector<std::string>>(
43  "system_variables",
44  std::vector<std::string>(),
45  "The name of the nodal and elemental variables from the file you want to use for values");
46 
47  // When using XDA files the following must be defined
48  params.addParam<FileName>(
49  "es",
50  "<not supplied>",
51  "The name of the file holding the equation system info in xda format (xda only).");
52  params.addParam<std::string>(
53  "system", "nl0", "The name of the system to pull values out of (xda only).");
54 
55  // When using ExodusII a specific time is extracted
56  params.addParam<std::string>("timestep",
57  "Index of the single timestep used or \"LATEST\" for "
58  "the last timestep (exodusII only). If not supplied, "
59  "time interpolation will occur.");
60 
61  // Add ability to perform coordinate transformation: scale, factor
62  params.addParam<std::vector<Real>>(
63  "scale", std::vector<Real>(LIBMESH_DIM, 1), "Scale factor for points in the simulation");
64  params.addParam<std::vector<Real>>("scale_multiplier",
65  std::vector<Real>(LIBMESH_DIM, 1),
66  "Scale multiplying factor for points in the simulation");
67  params.addParam<std::vector<Real>>("translation",
68  std::vector<Real>(LIBMESH_DIM, 0),
69  "Translation factors for x,y,z coordinates of the simulation");
70  params.addParam<RealVectorValue>("rotation0_vector",
71  RealVectorValue(0, 0, 1),
72  "Vector about which to rotate points of the simulation.");
73  params.addParam<Real>(
74  "rotation0_angle",
75  0.0,
76  "Anticlockwise rotation angle (in degrees) to use for rotation about rotation0_vector.");
77  params.addParam<RealVectorValue>("rotation1_vector",
78  RealVectorValue(0, 0, 1),
79  "Vector about which to rotate points of the simulation.");
80  params.addParam<Real>(
81  "rotation1_angle",
82  0.0,
83  "Anticlockwise rotation angle (in degrees) to use for rotation about rotation1_vector.");
84 
85  // following lines build the default_transformation_order
86  MultiMooseEnum default_transformation_order(
87  "rotation0 translation scale rotation1 scale_multiplier", "translation scale");
88  params.addParam<MultiMooseEnum>(
89  "transformation_order",
90  default_transformation_order,
91  "The order to perform the operations in. Define R0 to be the rotation matrix encoded by "
92  "rotation0_vector and rotation0_angle. Similarly for R1. Denote the scale by s, the "
93  "scale_multiplier by m, and the translation by t. Then, given a point x in the simulation, "
94  "if transformation_order = 'rotation0 scale_multiplier translation scale rotation1' then "
95  "form p = R1*(R0*x*m - t)/s. Then the values provided by the SolutionUserObject at point x "
96  "in the simulation are the variable values at point p in the mesh.");
97  params.addClassDescription("Reads a variable from a mesh in one simulation to another");
98  // Return the parameters
99  return params;
100 }
101 
102 // Static mutex definition
104 
106  : GeneralUserObject(parameters),
107  _file_type(MooseEnum("xda=0 exodusII=1 xdr=2")),
108  _mesh_file(getParam<MeshFileName>("mesh")),
109  _es_file(getParam<FileName>("es")),
110  _system_name(getParam<std::string>("system")),
111  _system_variables(getParam<std::vector<std::string>>("system_variables")),
112  _exodus_time_index(-1),
113  _interpolate_times(false),
114  _system(nullptr),
115  _system2(nullptr),
116  _interpolation_time(0.0),
117  _interpolation_factor(0.0),
118  _exodus_times(nullptr),
119  _exodus_index1(-1),
120  _exodus_index2(-1),
121  _scale(getParam<std::vector<Real>>("scale")),
122  _scale_multiplier(getParam<std::vector<Real>>("scale_multiplier")),
123  _translation(getParam<std::vector<Real>>("translation")),
124  _rotation0_vector(getParam<RealVectorValue>("rotation0_vector")),
125  _rotation0_angle(getParam<Real>("rotation0_angle")),
126  _r0(RealTensorValue()),
127  _rotation1_vector(getParam<RealVectorValue>("rotation1_vector")),
128  _rotation1_angle(getParam<Real>("rotation1_angle")),
129  _r1(RealTensorValue()),
130  _transformation_order(getParam<MultiMooseEnum>("transformation_order")),
131  _initialized(false)
132 {
133  // form rotation matrices with the specified angles
134  Real halfPi = std::acos(0.0);
135  Real a;
136  Real b;
137 
138  a = std::cos(halfPi * -_rotation0_angle / 90);
139  b = std::sin(halfPi * -_rotation0_angle / 90);
140  // the following is an anticlockwise rotation about z
141  RealTensorValue rot0_z(a, -b, 0, b, a, 0, 0, 0, 1);
142  // form the rotation matrix that will take rotation0_vector to the z axis
144  // _r0 is then: rotate points so vec0 lies along z; then rotate about angle0; then rotate points
145  // back
146  _r0 = vec0_to_z.transpose() * (rot0_z * vec0_to_z);
147 
148  a = std::cos(halfPi * -_rotation1_angle / 90);
149  b = std::sin(halfPi * -_rotation1_angle / 90);
150  // the following is an anticlockwise rotation about z
151  RealTensorValue rot1_z(a, -b, 0, b, a, 0, 0, 0, 1);
152  // form the rotation matrix that will take rotation1_vector to the z axis
154  // _r1 is then: rotate points so vec1 lies along z; then rotate about angle1; then rotate points
155  // back
156  _r1 = vec1_to_z.transpose() * (rot1_z * vec1_to_z);
157 
158  if (isParamValid("timestep") && getParam<std::string>("timestep") == "-1")
159  mooseError("A \"timestep\" of -1 is no longer supported for interpolation. Instead simply "
160  "remove this parameter altogether for interpolation");
161 }
162 
164 
165 void
167 {
168  // Check that the required files exist
171 
172  // Read the libmesh::mesh from the xda file
173  _mesh->read(_mesh_file);
174 
175  // Create the libmesh::EquationSystems
176  _es = libmesh_make_unique<EquationSystems>(*_mesh);
177 
178  // Use new read syntax (binary)
179  if (_file_type == "xdr")
180  _es->read(_es_file,
181  DECODE,
182  EquationSystems::READ_HEADER | EquationSystems::READ_DATA |
183  EquationSystems::READ_ADDITIONAL_DATA);
184 
185  // Use new read syntax
186  else if (_file_type == "xda")
187  _es->read(_es_file,
188  READ,
189  EquationSystems::READ_HEADER | EquationSystems::READ_DATA |
190  EquationSystems::READ_ADDITIONAL_DATA);
191 
192  // This should never occur, just in case produce an error
193  else
194  mooseError("Failed to determine proper read method for XDA/XDR equation system file: ",
195  _es_file);
196 
197  // Update and store the EquationSystems name locally
198  _es->update();
199  _system = &_es->get_system(_system_name);
200 }
201 
202 void
204 {
205  // Define a default system name
206  if (_system_name == "")
207  _system_name = "SolutionUserObjectSystem";
208 
209  // Read the Exodus file
210  _exodusII_io = libmesh_make_unique<ExodusII_IO>(*_mesh);
211  _exodusII_io->read(_mesh_file);
212  _exodus_times = &_exodusII_io->get_time_steps();
213 
214  if (isParamValid("timestep"))
215  {
216  std::string s_timestep = getParam<std::string>("timestep");
217  int n_steps = _exodusII_io->get_num_time_steps();
218  if (s_timestep == "LATEST")
219  _exodus_time_index = n_steps;
220  else
221  {
222  std::istringstream ss(s_timestep);
223  if (!((ss >> _exodus_time_index) && ss.eof()) || _exodus_time_index > n_steps)
224  mooseError("Invalid value passed as \"timestep\". Expected \"LATEST\" or a valid integer "
225  "less than ",
226  n_steps,
227  ", received ",
228  s_timestep);
229  }
230  }
231  else
232  // Interpolate between times rather than using values from a set timestep
233  _interpolate_times = true;
234 
235  // Check that the number of time steps is valid
236  int num_exo_times = _exodus_times->size();
237  if (num_exo_times == 0)
238  mooseError("In SolutionUserObject, exodus file contains no timesteps.");
239 
240  // Account for parallel mesh
241  if (dynamic_cast<DistributedMesh *>(_mesh.get()))
242  {
243  _mesh->allow_renumbering(true);
244  _mesh->prepare_for_use(/*false*/);
245  }
246  else
247  {
248  _mesh->allow_renumbering(false);
249  _mesh->prepare_for_use(/*true*/);
250  }
251 
252  // Create EquationSystems object for solution
253  _es = libmesh_make_unique<EquationSystems>(*_mesh);
254  _es->add_system<ExplicitSystem>(_system_name);
255  _system = &_es->get_system(_system_name);
256 
257  // Get the variable name lists as set; these need to be sets to perform set_intersection
258  const std::vector<std::string> & all_nodal(_exodusII_io->get_nodal_var_names());
259  const std::vector<std::string> & all_elemental(_exodusII_io->get_elem_var_names());
260 
261  // Storage for the nodal and elemental variables to consider
262  std::vector<std::string> nodal, elemental;
263 
264  // Build nodal/elemental variable lists, limit to variables listed in 'system_variables', if
265  // provided
266  if (!_system_variables.empty())
267  {
268  for (const auto & var_name : _system_variables)
269  {
270  if (std::find(all_nodal.begin(), all_nodal.end(), var_name) != all_nodal.end())
271  nodal.push_back(var_name);
272  if (std::find(all_elemental.begin(), all_elemental.end(), var_name) != all_elemental.end())
273  elemental.push_back(var_name);
274  }
275  }
276  else
277  {
278  nodal = all_nodal;
279  elemental = all_elemental;
280  }
281 
282  // Add the variables to the system
283  for (const auto & var_name : nodal)
284  _system->add_variable(var_name, FIRST);
285 
286  for (const auto & var_name : elemental)
287  _system->add_variable(var_name, CONSTANT, MONOMIAL);
288 
289  // Initialize the equations systems
290  _es->init();
291 
292  // Interpolate times
293  if (_interpolate_times)
294  {
295  // Create a second equation system
296  _es2 = libmesh_make_unique<EquationSystems>(*_mesh);
297  _es2->add_system<ExplicitSystem>(_system_name);
298  _system2 = &_es2->get_system(_system_name);
299 
300  // Add the variables to the system
301  for (const auto & var_name : nodal)
302  _system2->add_variable(var_name, FIRST);
303 
304  for (const auto & var_name : elemental)
305  _system2->add_variable(var_name, CONSTANT, MONOMIAL);
306 
307  // Initialize
308  _es2->init();
309 
310  // Update the times for interpolation (initially start at 0)
312 
313  // Copy the solutions from the first system
314  for (const auto & var_name : nodal)
315  {
316  _exodusII_io->copy_nodal_solution(*_system, var_name, var_name, _exodus_index1 + 1);
317  _exodusII_io->copy_nodal_solution(*_system2, var_name, var_name, _exodus_index2 + 1);
318  }
319 
320  for (const auto & var_name : elemental)
321  {
322  _exodusII_io->copy_elemental_solution(*_system, var_name, var_name, _exodus_index1 + 1);
323  _exodusII_io->copy_elemental_solution(*_system2, var_name, var_name, _exodus_index2 + 1);
324  }
325 
326  // Update the systems
327  _system->update();
328  _es->update();
329  _system2->update();
330  _es2->update();
331  }
332 
333  // Non-interpolated times
334  else
335  {
336  if (_exodus_time_index > num_exo_times)
337  mooseError("In SolutionUserObject, timestep = ",
339  ", but there are only ",
340  num_exo_times,
341  " time steps.");
342 
343  // Copy the values from the ExodusII file
344  for (const auto & var_name : nodal)
345  _exodusII_io->copy_nodal_solution(*_system, var_name, var_name, _exodus_time_index);
346 
347  for (const auto & var_name : elemental)
348  _exodusII_io->copy_elemental_solution(*_system, var_name, var_name, _exodus_time_index);
349 
350  // Update the equations systems
351  _system->update();
352  _es->update();
353  }
354 }
355 
356 Real
357 SolutionUserObject::directValue(const Node * node, const std::string & var_name) const
358 {
359  // Get the libmesh variable and system numbers
360  unsigned int var_num = _system->variable_number(var_name);
361  unsigned int sys_num = _system->number();
362 
363  // Get the node id and associated dof
364  dof_id_type node_id = node->id();
365  dof_id_type dof_id = _system->get_mesh().node_ref(node_id).dof_number(sys_num, var_num, 0);
366 
367  // Return the desired value for the dof
368  return directValue(dof_id);
369 }
370 
371 Real
372 SolutionUserObject::directValue(const Elem * elem, const std::string & var_name) const
373 {
374  // Get the libmesh variable and system numbers
375  unsigned int var_num = _system->variable_number(var_name);
376  unsigned int sys_num = _system->number();
377 
378  // Get the element id and associated dof
379  dof_id_type elem_id = elem->id();
380  dof_id_type dof_id = _system->get_mesh().elem_ptr(elem_id)->dof_number(sys_num, var_num, 0);
381 
382  // Return the desired value
383  return directValue(dof_id);
384 }
385 
386 void
388 {
389 }
390 
391 void
393 {
394 }
395 
396 void
398 {
399  // Update time interpolation for ExodusII solution
400  if (_file_type == 1 && _interpolate_times)
402 }
403 
404 void
406 {
407 }
408 
409 void
411 {
412 
413  // Make sure this only happens once
414  if (_initialized)
415  return;
416 
417  // Create a libmesh::Mesh object for storing the loaded data.
418  // Several aspects of SolutionUserObject won't work with a DistributedMesh:
419  // .) ExodusII_IO::copy_nodal_solution() doesn't work in parallel.
420  // .) We don't know if directValue will be used, which may request
421  // a value on a Node we don't have.
422  // We force the Mesh used here to be a ReplicatedMesh.
423  _mesh = libmesh_make_unique<ReplicatedMesh>(_communicator);
424 
425  // ExodusII mesh file supplied
426  if (MooseUtils::hasExtension(_mesh_file, "e", /*strip_exodus_ext =*/true))
427  {
428  _file_type = "exodusII";
429  readExodusII();
430  }
431 
432  // XDA mesh file supplied
433  else if (MooseUtils::hasExtension(_mesh_file, "xda"))
434  {
435  _file_type = "xda";
436  readXda();
437  }
438 
439  else if (MooseUtils::hasExtension(_mesh_file, "xdr"))
440  {
441  _file_type = "xdr";
442  readXda();
443  }
444 
445  // Produce an error for an unknown file type
446  else
447  mooseError("In SolutionUserObject, invalid file type (only .xda, .xdr, and .e supported)");
448 
449  // Intilize the serial solution vector
450  _serialized_solution = NumericVector<Number>::build(_communicator);
451  _serialized_solution->init(_system->n_dofs(), false, SERIAL);
452 
453  // Pull down a full copy of this vector on every processor so we can get values in parallel
454  _system->solution->localize(*_serialized_solution);
455 
456  // Vector of variable numbers to apply the MeshFunction to
457  std::vector<unsigned int> var_nums;
458 
459  // If no variables were given, use all of them
460  if (_system_variables.empty())
461  {
462  _system->get_all_variable_numbers(var_nums);
463  for (const auto & var_num : var_nums)
464  _system_variables.push_back(_system->variable_name(var_num));
465  }
466 
467  // Otherwise, gather the numbers for the variables given
468  else
469  {
470  for (const auto & var_name : _system_variables)
471  var_nums.push_back(_system->variable_number(var_name));
472  }
473 
474  // Create the MeshFunction for working with the solution data
475  _mesh_function = libmesh_make_unique<MeshFunction>(
476  *_es, *_serialized_solution, _system->get_dof_map(), var_nums);
477  _mesh_function->init();
478 
479  // Tell the MeshFunctions that we might be querying them outside the
480  // mesh, so we can handle any errors at the MOOSE rather than at the
481  // libMesh level.
482  DenseVector<Number> default_values;
483  _mesh_function->enable_out_of_mesh_mode(default_values);
484 
485  // Build second MeshFunction for interpolation
486  if (_interpolate_times)
487  {
488  // Need to pull down a full copy of this vector on every processor so we can get values in
489  // parallel
490  _serialized_solution2 = NumericVector<Number>::build(_communicator);
491  _serialized_solution2->init(_system2->n_dofs(), false, SERIAL);
492  _system2->solution->localize(*_serialized_solution2);
493 
494  // Create the MeshFunction for the second copy of the data
495  _mesh_function2 = libmesh_make_unique<MeshFunction>(
496  *_es2, *_serialized_solution2, _system2->get_dof_map(), var_nums);
497  _mesh_function2->init();
498  _mesh_function2->enable_out_of_mesh_mode(default_values);
499  }
500 
501  // Populate the data maps that indicate if the variable is nodal and the MeshFunction variable
502  // index
503  for (unsigned int i = 0; i < _system_variables.size(); ++i)
504  {
505  std::string name = _system_variables[i];
506  FEType type = _system->variable_type(name);
507  if (type.order == CONSTANT)
508  _local_variable_nodal[name] = false;
509  else
510  _local_variable_nodal[name] = true;
511 
513  }
514 
515  // Set initialization flag
516  _initialized = true;
517 }
518 
519 MooseEnum
521 {
522  return _file_type;
523 }
524 
525 void
527 {
528  if (time != _interpolation_time)
529  {
531  {
532 
533  for (const auto & var_name : _system_variables)
534  {
535  if (_local_variable_nodal[var_name])
536  _exodusII_io->copy_nodal_solution(*_system, var_name, var_name, _exodus_index1 + 1);
537  else
538  _exodusII_io->copy_elemental_solution(*_system, var_name, var_name, _exodus_index1 + 1);
539  }
540 
541  _system->update();
542  _es->update();
543  _system->solution->localize(*_serialized_solution);
544 
545  for (const auto & var_name : _system_variables)
546  {
547  if (_local_variable_nodal[var_name])
548  _exodusII_io->copy_nodal_solution(*_system2, var_name, var_name, _exodus_index2 + 1);
549  else
550  _exodusII_io->copy_elemental_solution(*_system2, var_name, var_name, _exodus_index2 + 1);
551  }
552 
553  _system2->update();
554  _es2->update();
555  _system2->solution->localize(*_serialized_solution2);
556  }
557  _interpolation_time = time;
558  }
559 }
560 
561 bool
563 {
564  if (_file_type != 1)
565  mooseError(
566  "In SolutionUserObject, getTimeInterpolationData only applicable for exodusII file type");
567 
568  int old_index1 = _exodus_index1;
569  int old_index2 = _exodus_index2;
570 
571  int num_exo_times = _exodus_times->size();
572 
573  if (time < (*_exodus_times)[0])
574  {
575  _exodus_index1 = 0;
576  _exodus_index2 = 0;
577  _interpolation_factor = 0.0;
578  }
579  else
580  {
581  for (int i = 0; i < num_exo_times - 1; ++i)
582  {
583  if (time <= (*_exodus_times)[i + 1])
584  {
585  _exodus_index1 = i;
586  _exodus_index2 = i + 1;
588  (time - (*_exodus_times)[i]) / ((*_exodus_times)[i + 1] - (*_exodus_times)[i]);
589  break;
590  }
591  else if (i == num_exo_times - 2)
592  {
593  _exodus_index1 = num_exo_times - 1;
594  _exodus_index2 = num_exo_times - 1;
595  _interpolation_factor = 1.0;
596  break;
597  }
598  }
599  }
600 
601  bool indices_modified(false);
602 
603  if (_exodus_index1 != old_index1 || _exodus_index2 != old_index2)
604  indices_modified = true;
605 
606  return indices_modified;
607 }
608 
609 unsigned int
610 SolutionUserObject::getLocalVarIndex(const std::string & var_name) const
611 {
612  // Extract the variable index for the MeshFunction(s)
613  std::map<std::string, unsigned int>::const_iterator it = _local_variable_index.find(var_name);
614  if (it == _local_variable_index.end())
615  mooseError("Value requested for nonexistent variable '",
616  var_name,
617  "' in the '",
618  name(),
619  "' SolutionUserObject");
620  return it->second;
621 }
622 
623 Real
625  const Point & p,
626  const std::string & var_name,
627  const MooseEnum & weighting_type) const
628 {
629  // first check if the FE type is continuous because in that case the value is
630  // unique and we can take a short cut, the default weighting_type found_first also
631  // shortcuts out
632  auto family =
634  .getVariable(
636  .feType()
637  .family;
638 
639  if (weighting_type == 1 ||
640  (family != L2_LAGRANGE && family != MONOMIAL && family != L2_HIERARCHIC))
641  return pointValue(t, p, var_name);
642 
643  // the shape function is discontinuous so we need to compute a suitable unique value
644  std::map<const Elem *, Real> values = discontinuousPointValue(t, p, var_name);
645  switch (weighting_type)
646  {
647  case 2:
648  {
649  Real average = 0.0;
650  for (auto & v : values)
651  average += v.second;
652  return average / Real(values.size());
653  }
654  case 4:
655  {
656  Real smallest_elem_id_value = std::numeric_limits<Real>::max();
657  dof_id_type smallest_elem_id = _fe_problem.mesh().maxElemId();
658  for (auto & v : values)
659  if (v.first->id() < smallest_elem_id)
660  {
661  smallest_elem_id = v.first->id();
662  smallest_elem_id_value = v.second;
663  }
664  return smallest_elem_id_value;
665  }
666  case 8:
667  {
668  Real largest_elem_id_value = std::numeric_limits<Real>::lowest();
669  dof_id_type largest_elem_id = 0;
670  for (auto & v : values)
671  if (v.first->id() > largest_elem_id)
672  {
673  largest_elem_id = v.first->id();
674  largest_elem_id_value = v.second;
675  }
676  return largest_elem_id_value;
677  }
678  }
679 
680  mooseError(
681  "SolutionUserObject::pointValueWrapper reaches line that it should not be able to reach.");
682  return 0.0;
683 }
684 
685 Real
686 SolutionUserObject::pointValue(Real t, const Point & p, const std::string & var_name) const
687 {
688  const unsigned int local_var_index = getLocalVarIndex(var_name);
689  return pointValue(t, p, local_var_index);
690 }
691 
692 Real
693 SolutionUserObject::pointValue(Real libmesh_dbg_var(t),
694  const Point & p,
695  const unsigned int local_var_index) const
696 {
697  // Create copy of point
698  Point pt(p);
699 
700  // do the transformations
701  for (unsigned int trans_num = 0; trans_num < _transformation_order.size(); ++trans_num)
702  {
703  if (_transformation_order[trans_num] == "rotation0")
704  pt = _r0 * pt;
705  else if (_transformation_order[trans_num] == "translation")
706  for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
707  pt(i) -= _translation[i];
708  else if (_transformation_order[trans_num] == "scale")
709  for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
710  pt(i) /= _scale[i];
711  else if (_transformation_order[trans_num] == "scale_multiplier")
712  for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
713  pt(i) *= _scale_multiplier[i];
714  else if (_transformation_order[trans_num] == "rotation1")
715  pt = _r1 * pt;
716  }
717 
718  // Extract the value at the current point
719  Real val = evalMeshFunction(pt, local_var_index, 1);
720 
721  // Interpolate
722  if (_file_type == 1 && _interpolate_times)
723  {
724  mooseAssert(t == _interpolation_time,
725  "Time passed into value() must match time at last call to timestepSetup()");
726  Real val2 = evalMeshFunction(pt, local_var_index, 2);
727  val = val + (val2 - val) * _interpolation_factor;
728  }
729 
730  return val;
731 }
732 
733 std::map<const Elem *, Real>
735  const Point & p,
736  const std::string & var_name) const
737 {
738  const unsigned int local_var_index = getLocalVarIndex(var_name);
739  return discontinuousPointValue(t, p, local_var_index);
740 }
741 
742 std::map<const Elem *, Real>
743 SolutionUserObject::discontinuousPointValue(Real libmesh_dbg_var(t),
744  Point pt,
745  const unsigned int local_var_index) const
746 {
747  // do the transformations
748  for (unsigned int trans_num = 0; trans_num < _transformation_order.size(); ++trans_num)
749  {
750  if (_transformation_order[trans_num] == "rotation0")
751  pt = _r0 * pt;
752  else if (_transformation_order[trans_num] == "translation")
753  for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
754  pt(i) -= _translation[i];
755  else if (_transformation_order[trans_num] == "scale")
756  for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
757  pt(i) /= _scale[i];
758  else if (_transformation_order[trans_num] == "scale_multiplier")
759  for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
760  pt(i) *= _scale_multiplier[i];
761  else if (_transformation_order[trans_num] == "rotation1")
762  pt = _r1 * pt;
763  }
764 
765  // Extract the value at the current point
766  std::map<const Elem *, Real> map = evalMultiValuedMeshFunction(pt, local_var_index, 1);
767 
768  // Interpolate
769  if (_file_type == 1 && _interpolate_times)
770  {
771  mooseAssert(t == _interpolation_time,
772  "Time passed into value() must match time at last call to timestepSetup()");
773  std::map<const Elem *, Real> map2 = evalMultiValuedMeshFunction(pt, local_var_index, 2);
774 
775  if (map.size() != map2.size())
776  mooseError("In SolutionUserObject::discontinuousPointValue map and map2 have different size");
777 
778  // construct the interpolated map
779  for (auto & k : map)
780  {
781  if (map2.find(k.first) == map2.end())
782  mooseError(
783  "In SolutionUserObject::discontinuousPointValue map and map2 have differing keys");
784  Real val = k.second;
785  Real val2 = map2[k.first];
786  map[k.first] = val + (val2 - val) * _interpolation_factor;
787  }
788  }
789 
790  return map;
791 }
792 
793 RealGradient
795  const Point & p,
796  const std::string & var_name,
797  const MooseEnum & weighting_type) const
798 {
799  // the default weighting_type found_first shortcuts out
800  if (weighting_type == 1)
801  return pointValueGradient(t, p, var_name);
802 
803  // the shape function is discontinuous so we need to compute a suitable unique value
804  std::map<const Elem *, RealGradient> values = discontinuousPointValueGradient(t, p, var_name);
805  switch (weighting_type)
806  {
807  case 2:
808  {
809  RealGradient average = RealGradient(0.0, 0.0, 0.0);
810  for (auto & v : values)
811  average += v.second;
812  return average / Real(values.size());
813  }
814  case 4:
815  {
816  RealGradient smallest_elem_id_value;
817  dof_id_type smallest_elem_id = _fe_problem.mesh().maxElemId();
818  for (auto & v : values)
819  if (v.first->id() < smallest_elem_id)
820  {
821  smallest_elem_id = v.first->id();
822  smallest_elem_id_value = v.second;
823  }
824  return smallest_elem_id_value;
825  }
826  case 8:
827  {
828  RealGradient largest_elem_id_value;
829  dof_id_type largest_elem_id = 0;
830  for (auto & v : values)
831  if (v.first->id() > largest_elem_id)
832  {
833  largest_elem_id = v.first->id();
834  largest_elem_id_value = v.second;
835  }
836  return largest_elem_id_value;
837  }
838  }
839 
840  mooseError("SolutionUserObject::pointValueGradientWrapper reaches line that it should not be "
841  "able to reach.");
842  return RealGradient(0.0, 0.0, 0.0);
843 }
844 
845 RealGradient
846 SolutionUserObject::pointValueGradient(Real t, const Point & p, const std::string & var_name) const
847 {
848  const unsigned int local_var_index = getLocalVarIndex(var_name);
849  return pointValueGradient(t, p, local_var_index);
850 }
851 
852 RealGradient
853 SolutionUserObject::pointValueGradient(Real libmesh_dbg_var(t),
854  Point pt,
855  const unsigned int local_var_index) const
856 {
857  // do the transformations
858  for (unsigned int trans_num = 0; trans_num < _transformation_order.size(); ++trans_num)
859  {
860  if (_transformation_order[trans_num] == "rotation0")
861  pt = _r0 * pt;
862  else if (_transformation_order[trans_num] == "translation")
863  for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
864  pt(i) -= _translation[i];
865  else if (_transformation_order[trans_num] == "scale")
866  for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
867  pt(i) /= _scale[i];
868  else if (_transformation_order[trans_num] == "scale_multiplier")
869  for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
870  pt(i) *= _scale_multiplier[i];
871  else if (_transformation_order[trans_num] == "rotation1")
872  pt = _r1 * pt;
873  }
874 
875  // Extract the value at the current point
876  RealGradient val = evalMeshFunctionGradient(pt, local_var_index, 1);
877 
878  // Interpolate
879  if (_file_type == 1 && _interpolate_times)
880  {
881  mooseAssert(t == _interpolation_time,
882  "Time passed into value() must match time at last call to timestepSetup()");
883  RealGradient val2 = evalMeshFunctionGradient(pt, local_var_index, 2);
884  val = val + (val2 - val) * _interpolation_factor;
885  }
886 
887  return val;
888 }
889 
890 std::map<const Elem *, RealGradient>
892  const Point & p,
893  const std::string & var_name) const
894 {
895  const unsigned int local_var_index = getLocalVarIndex(var_name);
896  return discontinuousPointValueGradient(t, p, local_var_index);
897 }
898 
899 std::map<const Elem *, RealGradient>
901  Point pt,
902  const unsigned int local_var_index) const
903 {
904  // do the transformations
905  for (unsigned int trans_num = 0; trans_num < _transformation_order.size(); ++trans_num)
906  {
907  if (_transformation_order[trans_num] == "rotation0")
908  pt = _r0 * pt;
909  else if (_transformation_order[trans_num] == "translation")
910  for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
911  pt(i) -= _translation[i];
912  else if (_transformation_order[trans_num] == "scale")
913  for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
914  pt(i) /= _scale[i];
915  else if (_transformation_order[trans_num] == "scale_multiplier")
916  for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
917  pt(i) *= _scale_multiplier[i];
918  else if (_transformation_order[trans_num] == "rotation1")
919  pt = _r1 * pt;
920  }
921 
922  // Extract the value at the current point
923  std::map<const Elem *, RealGradient> map =
924  evalMultiValuedMeshFunctionGradient(pt, local_var_index, 1);
925 
926  // Interpolate
927  if (_file_type == 1 && _interpolate_times)
928  {
929  mooseAssert(t == _interpolation_time,
930  "Time passed into value() must match time at last call to timestepSetup()");
931  std::map<const Elem *, RealGradient> map2 =
932  evalMultiValuedMeshFunctionGradient(pt, local_var_index, 1);
933 
934  if (map.size() != map2.size())
935  mooseError("In SolutionUserObject::discontinuousPointValue map and map2 have different size");
936 
937  // construct the interpolated map
938  for (auto & k : map)
939  {
940  if (map2.find(k.first) == map2.end())
941  mooseError(
942  "In SolutionUserObject::discontinuousPointValue map and map2 have differing keys");
943  RealGradient val = k.second;
944  RealGradient val2 = map2[k.first];
945  map[k.first] = val + (val2 - val) * _interpolation_factor;
946  }
947  }
948 
949  return map;
950 }
951 
952 Real
953 SolutionUserObject::directValue(dof_id_type dof_index) const
954 {
955  Real val = (*_serialized_solution)(dof_index);
956  if (_file_type == 1 && _interpolate_times)
957  {
958  Real val2 = (*_serialized_solution2)(dof_index);
959  val = val + (val2 - val) * _interpolation_factor;
960  }
961  return val;
962 }
963 
964 Real
966  const unsigned int local_var_index,
967  unsigned int func_num) const
968 {
969  // Storage for mesh function output
970  DenseVector<Number> output;
971 
972  // Extract a value from the _mesh_function
973  {
974  Threads::spin_mutex::scoped_lock lock(_solution_user_object_mutex);
975  if (func_num == 1)
976  (*_mesh_function)(p, 0.0, output);
977 
978  // Extract a value from _mesh_function2
979  else if (func_num == 2)
980  (*_mesh_function2)(p, 0.0, output);
981 
982  else
983  mooseError("The func_num must be 1 or 2");
984  }
985 
986  // Error if the data is out-of-range, which will be the case if the mesh functions are evaluated
987  // outside the domain
988  if (output.size() == 0)
989  {
990  std::ostringstream oss;
991  p.print(oss);
992  mooseError("Failed to access the data for variable '",
993  _system_variables[local_var_index],
994  "' at point ",
995  oss.str(),
996  " in the '",
997  name(),
998  "' SolutionUserObject");
999  }
1000  return output(local_var_index);
1001 }
1002 
1003 std::map<const Elem *, Real>
1005  const unsigned int local_var_index,
1006  unsigned int func_num) const
1007 {
1008  // Storage for mesh function output
1009  std::map<const Elem *, DenseVector<Number>> temporary_output;
1010 
1011  // Extract a value from the _mesh_function
1012  {
1013  Threads::spin_mutex::scoped_lock lock(_solution_user_object_mutex);
1014  if (func_num == 1)
1015  _mesh_function->discontinuous_value(p, 0.0, temporary_output);
1016 
1017  // Extract a value from _mesh_function2
1018  else if (func_num == 2)
1019  _mesh_function2->discontinuous_value(p, 0.0, temporary_output);
1020 
1021  else
1022  mooseError("The func_num must be 1 or 2");
1023  }
1024 
1025  // Error if the data is out-of-range, which will be the case if the mesh functions are evaluated
1026  // outside the domain
1027  if (temporary_output.size() == 0)
1028  {
1029  std::ostringstream oss;
1030  p.print(oss);
1031  mooseError("Failed to access the data for variable '",
1032  _system_variables[local_var_index],
1033  "' at point ",
1034  oss.str(),
1035  " in the '",
1036  name(),
1037  "' SolutionUserObject");
1038  }
1039 
1040  // Fill the actual map that is returned
1041  std::map<const Elem *, Real> output;
1042  for (auto & k : temporary_output)
1043  {
1044  mooseAssert(k.second.size() > local_var_index,
1045  "In SolutionUserObject::evalMultiValuedMeshFunction variable with local_var_index "
1046  << local_var_index
1047  << " does not exist");
1048  output[k.first] = k.second(local_var_index);
1049  }
1050 
1051  return output;
1052 }
1053 
1054 RealGradient
1056  const unsigned int local_var_index,
1057  unsigned int func_num) const
1058 {
1059  // Storage for mesh function output
1060  std::vector<Gradient> output;
1061 
1062  // Extract a value from the _mesh_function
1063  {
1064  Threads::spin_mutex::scoped_lock lock(_solution_user_object_mutex);
1065  if (func_num == 1)
1066  _mesh_function->gradient(p, 0.0, output, libmesh_nullptr);
1067 
1068  // Extract a value from _mesh_function2
1069  else if (func_num == 2)
1070  _mesh_function2->gradient(p, 0.0, output, libmesh_nullptr);
1071 
1072  else
1073  mooseError("The func_num must be 1 or 2");
1074  }
1075 
1076  // Error if the data is out-of-range, which will be the case if the mesh functions are evaluated
1077  // outside the domain
1078  if (output.size() == 0)
1079  {
1080  std::ostringstream oss;
1081  p.print(oss);
1082  mooseError("Failed to access the data for variable '",
1083  _system_variables[local_var_index],
1084  "' at point ",
1085  oss.str(),
1086  " in the '",
1087  name(),
1088  "' SolutionUserObject");
1089  }
1090  return output[local_var_index];
1091 }
1092 
1093 std::map<const Elem *, RealGradient>
1095  const unsigned int local_var_index,
1096  unsigned int func_num) const
1097 {
1098  // Storage for mesh function output
1099  std::map<const Elem *, std::vector<Gradient>> temporary_output;
1100 
1101  // Extract a value from the _mesh_function
1102  {
1103  Threads::spin_mutex::scoped_lock lock(_solution_user_object_mutex);
1104  if (func_num == 1)
1105  _mesh_function->discontinuous_gradient(p, 0.0, temporary_output);
1106 
1107  // Extract a value from _mesh_function2
1108  else if (func_num == 2)
1109  _mesh_function2->discontinuous_gradient(p, 0.0, temporary_output);
1110 
1111  else
1112  mooseError("The func_num must be 1 or 2");
1113  }
1114 
1115  // Error if the data is out-of-range, which will be the case if the mesh functions are evaluated
1116  // outside the domain
1117  if (temporary_output.size() == 0)
1118  {
1119  std::ostringstream oss;
1120  p.print(oss);
1121  mooseError("Failed to access the data for variable '",
1122  _system_variables[local_var_index],
1123  "' at point ",
1124  oss.str(),
1125  " in the '",
1126  name(),
1127  "' SolutionUserObject");
1128  }
1129 
1130  // Fill the actual map that is returned
1131  std::map<const Elem *, RealGradient> output;
1132  for (auto & k : temporary_output)
1133  {
1134  mooseAssert(k.second.size() > local_var_index,
1135  "In SolutionUserObject::evalMultiValuedMeshFunction variable with local_var_index "
1136  << local_var_index
1137  << " does not exist");
1138  output[k.first] = k.second[local_var_index];
1139  }
1140 
1141  return output;
1142 }
1143 
1144 const std::vector<std::string> &
1146 {
1147  return _system_variables;
1148 }
1149 
1150 bool
1151 SolutionUserObject::isVariableNodal(const std::string & var_name) const
1152 {
1153  // Use iterator method, [] is not marked const
1154  std::map<std::string, bool>::const_iterator it = _local_variable_nodal.find(var_name);
1155  return it->second;
1156 }
void readXda()
Method for reading XDA mesh and equation systems file(s) This method is called by the constructor whe...
std::unique_ptr< ExodusII_IO > _exodusII_io
Pointer to the libMesh::ExodusII used to read the files.
const std::vector< std::string > & variableNames() const
Real pointValue(Real t, const Point &p, const unsigned int local_var_index) const
Returns a value at a specific location and variable (see SolutionFunction)
std::vector< Real > _scale_multiplier
scale_multiplier parameter
std::unique_ptr< NumericVector< Number > > _serialized_solution
Pointer to the serial solution vector.
virtual MooseVariableFEBase & getVariable(THREAD_ID tid, const std::string &var_name, Moose::VarKindType expected_var_type=Moose::VarKindType::VAR_ANY, Moose::VarFieldType expected_var_field_type=Moose::VarFieldType::VAR_FIELD_ANY) override
Returns the variable reference for requested variable which must be of the expected_var_type (Nonline...
std::vector< std::string > _system_variables
A list of variables to extract from the read system.
VectorValue< Real > RealVectorValue
Definition: Assembly.h:31
std::string _es_file
The XDA file that contians the EquationSystems data (xda only)
RealTensorValue rotVecToZ(RealVectorValue vec)
provides a rotation matrix that will rotate the vector vec to the z axis (the "2" direction) ...
registerMooseObject("MooseApp", SolutionUserObject)
unsigned int getLocalVarIndex(const std::string &var_name) const
Returns the local index for a given variable name.
std::unique_ptr< EquationSystems > _es
Pointer to the libmesh::EquationSystems object.
void updateExodusTimeInterpolation(Real time)
Updates the times for interpolating ExodusII data.
unsigned int size() const
Return the number of active items in the MultiMooseEnum.
The main MOOSE class responsible for handling user-defined parameters in almost every MOOSE system...
virtual dof_id_type maxElemId() const
Definition: MooseMesh.C:2247
RealVectorValue _rotation1_vector
vector about which to rotate
RealGradient pointValueGradientWrapper(Real t, const Point &p, const std::string &var_name, const MooseEnum &weighting_type=weightingType()) const
Returns the gradient at a specific location and variable checking for multiple values and weighting t...
THREAD_ID _tid
Thread ID of this postprocessor.
Definition: UserObject.h:144
void mooseError(Args &&... args) const
Definition: MooseObject.h:147
virtual void finalize() override
Finalize.
InputParameters validParams< GeneralUserObject >()
const std::vector< Real > * _exodus_times
The times available in the ExodusII file.
RealVectorValue _rotation0_vector
vector about which to rotate
const FEType & feType() const
Get the type of finite element object.
void addRequiredParam(const std::string &name, const std::string &doc_string)
This method adds a parameter and documentation string to the InputParameters object that will be extr...
int _exodus_index1
Time index 1, used for interpolation.
Real directValue(const Node *node, const std::string &var_name) const
Return a value directly from a Node.
const std::string & type() const
Get the type of this object.
Definition: MooseObject.h:53
virtual void timestepSetup() override
When reading ExodusII files, this will update the interpolation times.
bool checkFileReadable(const std::string &filename, bool check_line_endings=false, bool throw_on_unreadable=true)
Checks to see if a file is readable (exists and permissions)
Definition: MooseUtils.C:146
int _exodus_time_index
Current ExodusII time index.
Real _rotation0_angle
angle (in degrees) which to rotate through about vector _rotation0_vector
RealTensorValue _r0
Rotation matrix that performs the "_rotation0_angle about rotation0_vector".
Real pointValueWrapper(Real t, const Point &p, const std::string &var_name, const MooseEnum &weighting_type=weightingType()) const
Returns a value at a specific location and variable checking for multiple values and weighting these ...
void readExodusII()
Method for reading an ExodusII file, which is called when &#39;file_type = exodusII is set in the input f...
std::map< const Elem *, RealGradient > evalMultiValuedMeshFunctionGradient(const Point &p, const unsigned int local_var_index, unsigned int func_num) const
A wrapper method interfacing with the libMesh mesh function that calls the gradient functionality for...
std::string _system_name
The system name to extract from the XDA file (xda only)
Real evalMeshFunction(const Point &p, const unsigned int local_var_index, unsigned int func_num) const
A wrapper method for calling the various MeshFunctions used for reading the data. ...
std::unique_ptr< NumericVector< Number > > _serialized_solution2
Pointer to second serial solution, used for interpolation.
RealTensorValue _r1
Rotation matrix that performs the "_rotation1_angle about rotation1_vector".
This is a "smart" enum class intended to replace many of the shortcomings in the C++ enum type It sho...
Definition: MooseEnum.h:31
bool hasExtension(const std::string &filename, std::string ext, bool strip_exodus_ext=false)
Function tests if the supplied filename as the desired extension.
Definition: MooseUtils.C:256
std::map< const Elem *, RealGradient > discontinuousPointValueGradient(Real t, const Point &p, const std::string &var_name) const
Returns the gradient at a specific location and variable for cases where the gradient is multivalued ...
std::unique_ptr< MeshFunction > _mesh_function2
Pointer to second libMesh::MeshFuntion, used for interpolation.
Real _interpolation_time
Interpolation time.
static Threads::spin_mutex _solution_user_object_mutex
System * _system2
Pointer to a second libMesh::System object, used for interpolation.
SolutionUserObject(const InputParameters &parameters)
MooseEnum getSolutionFileType()
Get the type of file that was read.
virtual void initialSetup() override
Initialize the System and Mesh objects for the solution being read.
std::map< std::string, unsigned int > _local_variable_index
Stores the local index need by MeshFunction.
RealGradient pointValueGradient(Real t, const Point &p, const std::string &var_name) const
Returns the gradient at a specific location and variable (see SolutionFunction)
std::map< const Elem *, Real > evalMultiValuedMeshFunction(const Point &p, const unsigned int local_var_index, unsigned int func_num) const
A wrapper method for calling the various MeshFunctions that calls the mesh function functionality for...
std::map< const Elem *, Real > discontinuousPointValue(Real t, Point pt, const unsigned int local_var_index) const
Returns a value at a specific location and variable for cases where the solution is multivalued at el...
bool _initialized
True if initial_setup has executed.
FEProblemBase & _fe_problem
Reference to the FEProblemBase for this user object.
Definition: UserObject.h:141
MultiMooseEnum _transformation_order
transformations (rotations, translation, scales) are performed in this order
std::unique_ptr< MeshFunction > _mesh_function
Pointer the libMesh::MeshFunction object that the read data is stored.
std::vector< Real > _scale
Scale parameter.
bool _interpolate_times
Flag for triggering interpolation of ExodusII data.
virtual MooseMesh & mesh() override
const std::string & name() const
Get the name of the object.
Definition: MooseObject.h:59
void addClassDescription(const std::string &doc_string)
This method adds a description of the class that will be displayed in the input file syntax dump...
void addParam(const std::string &name, const S &value, const std::string &doc_string)
These methods add an option parameter and a documentation string to the InputParameters object...
Real _interpolation_factor
Interpolation weight factor.
std::unique_ptr< EquationSystems > _es2
Pointer to second libMesh::EquationSystems object, used for interpolation.
virtual void execute() override
Execute method.
This is a "smart" enum class intended to replace many of the shortcomings in the C++ enum type It sho...
bool isParamValid(const std::string &name) const
Test if the supplied parameter is valid.
Definition: MooseObject.h:89
User object that reads an existing solution from an input file and uses it in the current simulation...
bool updateExodusBracketingTimeIndices(Real time)
Updates the time indices to interpolate between for ExodusII data.
TensorValue< Real > RealTensorValue
Definition: MooseTypes.h:132
MooseEnum _file_type
File type to read (0 = xda; 1 = ExodusII)
bool isVariableNodal(const std::string &var_name) const
int _exodus_index2
Time index 2, used for interpolation.
std::map< std::string, bool > _local_variable_nodal
Stores flag indicating if the variable is nodal.
InputParameters validParams< SolutionUserObject >()
std::string _mesh_file
The XDA or ExodusII file that is being read.
System * _system
Pointer libMesh::System class storing the read solution.
std::unique_ptr< MeshBase > _mesh
Pointer the libmesh::mesh object.
Real _rotation1_angle
angle (in degrees) which to rotate through about vector _rotation1_vector
virtual void initialize() override
Called before execute() is ever called so that data can be cleared.
std::vector< Real > _translation
Translation.
RealGradient evalMeshFunctionGradient(const Point &p, const unsigned int local_var_index, unsigned int func_num) const
A wrapper method interfacing with the libMesh mesh function for evaluating the gradient.