www.mooseframework.org
SolutionUserObject.C
Go to the documentation of this file.
1 /****************************************************************/
2 /* DO NOT MODIFY THIS HEADER */
3 /* MOOSE - Multiphysics Object Oriented Simulation Environment */
4 /* */
5 /* (c) 2010 Battelle Energy Alliance, LLC */
6 /* ALL RIGHTS RESERVED */
7 /* */
8 /* Prepared by Battelle Energy Alliance, LLC */
9 /* Under Contract No. DE-AC07-05ID14517 */
10 /* With the U. S. Department of Energy */
11 /* */
12 /* See COPYRIGHT for full restrictions */
13 /****************************************************************/
14 
15 #include "SolutionUserObject.h"
16 
17 // MOOSE includes
18 #include "MooseError.h"
19 #include "MooseMesh.h"
20 #include "MooseUtils.h"
21 #include "MooseVariable.h"
22 #include "RotationMatrix.h"
23 
24 #include "libmesh/equation_systems.h"
25 #include "libmesh/mesh_function.h"
26 #include "libmesh/numeric_vector.h"
27 #include "libmesh/nonlinear_implicit_system.h"
28 #include "libmesh/transient_system.h"
29 #include "libmesh/parallel_mesh.h"
30 #include "libmesh/serial_mesh.h"
31 #include "libmesh/exodusII_io.h"
32 
33 template <>
36 {
37  // Get the input parameters from the parent class
39 
40  // Add required parameters
41  params.addRequiredParam<MeshFileName>(
42  "mesh", "The name of the mesh file (must be xda or exodusII file).");
43  params.addParam<std::vector<std::string>>(
44  "system_variables",
45  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 files the following must be defined
49  params.addParam<FileName>(
50  "es",
51  "<not supplied>",
52  "The name of the file holding the equation system info in xda format (xda only).");
53  params.addParam<std::string>(
54  "system", "nl0", "The name of the system to pull values out of (xda only).");
55 
56  // When using ExodusII a specific time is extracted
57  params.addParam<std::string>("timestep",
58  "Index of the single timestep used or \"LATEST\" for "
59  "the last timestep (exodusII only). If not supplied, "
60  "time interpolation will occur.");
61 
62  // Add ability to perform coordinate transformation: scale, factor
63  params.addParam<std::vector<Real>>(
64  "scale", std::vector<Real>(LIBMESH_DIM, 1), "Scale factor for points in the simulation");
65  params.addParam<std::vector<Real>>("scale_multiplier",
66  std::vector<Real>(LIBMESH_DIM, 1),
67  "Scale multiplying factor for points in the simulation");
68  params.addParam<std::vector<Real>>("translation",
69  std::vector<Real>(LIBMESH_DIM, 0),
70  "Translation factors for x,y,z coordinates of the simulation");
71  params.addParam<RealVectorValue>("rotation0_vector",
72  RealVectorValue(0, 0, 1),
73  "Vector about which to rotate points of the simulation.");
74  params.addParam<Real>(
75  "rotation0_angle",
76  0.0,
77  "Anticlockwise rotation angle (in degrees) to use for rotation about rotation0_vector.");
78  params.addParam<RealVectorValue>("rotation1_vector",
79  RealVectorValue(0, 0, 1),
80  "Vector about which to rotate points of the simulation.");
81  params.addParam<Real>(
82  "rotation1_angle",
83  0.0,
84  "Anticlockwise rotation angle (in degrees) to use for rotation about rotation1_vector.");
85 
86  // following lines build the default_transformation_order
87  MultiMooseEnum default_transformation_order(
88  "rotation0 translation scale rotation1 scale_multiplier", "translation scale");
89  params.addParam<MultiMooseEnum>(
90  "transformation_order",
91  default_transformation_order,
92  "The order to perform the operations in. Define R0 to be the rotation matrix encoded by "
93  "rotation0_vector and rotation0_angle. Similarly for R1. Denote the scale by s, the "
94  "scale_multiplier by m, and the translation by t. Then, given a point x in the simulation, "
95  "if transformation_order = 'rotation0 scale_multiplier translation scale rotation1' then "
96  "form p = R1*(R0*x*m - t)/s. Then the values provided by the SolutionUserObject at point x "
97  "in the simulation are the variable values at point p in the mesh.");
98  params.addClassDescription("Reads a variable from a mesh in one simulation to another");
99  // Return the parameters
100  return params;
101 }
102 
103 // Static mutex definition
105 
107  : GeneralUserObject(parameters),
108  _file_type(MooseEnum("xda=0 exodusII=1 xdr=2")),
109  _mesh_file(getParam<MeshFileName>("mesh")),
110  _es_file(getParam<FileName>("es")),
111  _system_name(getParam<std::string>("system")),
112  _system_variables(getParam<std::vector<std::string>>("system_variables")),
113  _exodus_time_index(-1),
114  _interpolate_times(false),
115  _system(nullptr),
116  _system2(nullptr),
117  _interpolation_time(0.0),
118  _interpolation_factor(0.0),
119  _exodus_times(nullptr),
120  _exodus_index1(-1),
121  _exodus_index2(-1),
122  _scale(getParam<std::vector<Real>>("scale")),
123  _scale_multiplier(getParam<std::vector<Real>>("scale_multiplier")),
124  _translation(getParam<std::vector<Real>>("translation")),
125  _rotation0_vector(getParam<RealVectorValue>("rotation0_vector")),
126  _rotation0_angle(getParam<Real>("rotation0_angle")),
127  _r0(RealTensorValue()),
128  _rotation1_vector(getParam<RealVectorValue>("rotation1_vector")),
129  _rotation1_angle(getParam<Real>("rotation1_angle")),
130  _r1(RealTensorValue()),
131  _transformation_order(getParam<MultiMooseEnum>("transformation_order")),
132  _initialized(false)
133 {
134  // form rotation matrices with the specified angles
135  Real halfPi = std::acos(0.0);
136  Real a;
137  Real b;
138 
139  a = std::cos(halfPi * -_rotation0_angle / 90);
140  b = std::sin(halfPi * -_rotation0_angle / 90);
141  // the following is an anticlockwise rotation about z
142  RealTensorValue rot0_z(a, -b, 0, b, a, 0, 0, 0, 1);
143  // form the rotation matrix that will take rotation0_vector to the z axis
145  // _r0 is then: rotate points so vec0 lies along z; then rotate about angle0; then rotate points
146  // back
147  _r0 = vec0_to_z.transpose() * (rot0_z * vec0_to_z);
148 
149  a = std::cos(halfPi * -_rotation1_angle / 90);
150  b = std::sin(halfPi * -_rotation1_angle / 90);
151  // the following is an anticlockwise rotation about z
152  RealTensorValue rot1_z(a, -b, 0, b, a, 0, 0, 0, 1);
153  // form the rotation matrix that will take rotation1_vector to the z axis
155  // _r1 is then: rotate points so vec1 lies along z; then rotate about angle1; then rotate points
156  // back
157  _r1 = vec1_to_z.transpose() * (rot1_z * vec1_to_z);
158 
159  if (isParamValid("timestep") && getParam<std::string>("timestep") == "-1")
160  mooseError("A \"timestep\" of -1 is no longer supported for interpolation. Instead simply "
161  "remove this parameter altogether for interpolation");
162 }
163 
165 
166 void
168 {
169  // Check that the required files exist
172 
173  // Read the libmesh::mesh from the xda file
174  _mesh->read(_mesh_file);
175 
176  // Create the libmesh::EquationSystems
177  _es = libmesh_make_unique<EquationSystems>(*_mesh);
178 
179  // Use new read syntax (binary)
180  if (_file_type == "xdr")
181  _es->read(_es_file,
182  DECODE,
183  EquationSystems::READ_HEADER | EquationSystems::READ_DATA |
184  EquationSystems::READ_ADDITIONAL_DATA);
185 
186  // Use new read syntax
187  else if (_file_type == "xda")
188  _es->read(_es_file,
189  READ,
190  EquationSystems::READ_HEADER | EquationSystems::READ_DATA |
191  EquationSystems::READ_ADDITIONAL_DATA);
192 
193  // This should never occur, just in case produce an error
194  else
195  mooseError("Failed to determine proper read method for XDA/XDR equation system file: ",
196  _es_file);
197 
198  // Update and store the EquationSystems name locally
199  _es->update();
200  _system = &_es->get_system(_system_name);
201 }
202 
203 void
205 {
206  // Define a default system name
207  if (_system_name == "")
208  _system_name = "SolutionUserObjectSystem";
209 
210  // Read the Exodus file
211  _exodusII_io = libmesh_make_unique<ExodusII_IO>(*_mesh);
212  _exodusII_io->read(_mesh_file);
213  _exodus_times = &_exodusII_io->get_time_steps();
214 
215  if (isParamValid("timestep"))
216  {
217  std::string s_timestep = getParam<std::string>("timestep");
218  int n_steps = _exodusII_io->get_num_time_steps();
219  if (s_timestep == "LATEST")
220  _exodus_time_index = n_steps;
221  else
222  {
223  std::istringstream ss(s_timestep);
224  if (!((ss >> _exodus_time_index) && ss.eof()) || _exodus_time_index > n_steps)
225  mooseError("Invalid value passed as \"timestep\". Expected \"LATEST\" or a valid integer "
226  "less than ",
227  n_steps,
228  ", received ",
229  s_timestep);
230  }
231  }
232  else
233  // Interpolate between times rather than using values from a set timestep
234  _interpolate_times = true;
235 
236  // Check that the number of time steps is valid
237  int num_exo_times = _exodus_times->size();
238  if (num_exo_times == 0)
239  mooseError("In SolutionUserObject, exodus file contains no timesteps.");
240 
241  // Account for parallel mesh
242  if (dynamic_cast<DistributedMesh *>(_mesh.get()))
243  {
244  _mesh->allow_renumbering(true);
245  _mesh->prepare_for_use(/*false*/);
246  }
247  else
248  {
249  _mesh->allow_renumbering(false);
250  _mesh->prepare_for_use(/*true*/);
251  }
252 
253  // Create EquationSystems object for solution
254  _es = libmesh_make_unique<EquationSystems>(*_mesh);
255  _es->add_system<ExplicitSystem>(_system_name);
256  _system = &_es->get_system(_system_name);
257 
258  // Get the variable name lists as set; these need to be sets to perform set_intersection
259  const std::vector<std::string> & all_nodal(_exodusII_io->get_nodal_var_names());
260  const std::vector<std::string> & all_elemental(_exodusII_io->get_elem_var_names());
261 
262  // Storage for the nodal and elemental variables to consider
263  std::vector<std::string> nodal, elemental;
264 
265  // Build nodal/elemental variable lists, limit to variables listed in 'system_variables', if
266  // provided
267  if (!_system_variables.empty())
268  {
269  for (const auto & var_name : _system_variables)
270  {
271  if (std::find(all_nodal.begin(), all_nodal.end(), var_name) != all_nodal.end())
272  nodal.push_back(var_name);
273  if (std::find(all_elemental.begin(), all_elemental.end(), var_name) != all_elemental.end())
274  elemental.push_back(var_name);
275  }
276  }
277  else
278  {
279  nodal = all_nodal;
280  elemental = all_elemental;
281  }
282 
283  // Add the variables to the system
284  for (const auto & var_name : nodal)
285  _system->add_variable(var_name, FIRST);
286 
287  for (const auto & var_name : elemental)
288  _system->add_variable(var_name, CONSTANT, MONOMIAL);
289 
290  // Initialize the equations systems
291  _es->init();
292 
293  // Interpolate times
294  if (_interpolate_times)
295  {
296  // Create a second equation system
297  _es2 = libmesh_make_unique<EquationSystems>(*_mesh);
298  _es2->add_system<ExplicitSystem>(_system_name);
299  _system2 = &_es2->get_system(_system_name);
300 
301  // Add the variables to the system
302  for (const auto & var_name : nodal)
303  _system2->add_variable(var_name, FIRST);
304 
305  for (const auto & var_name : elemental)
306  _system2->add_variable(var_name, CONSTANT, MONOMIAL);
307 
308  // Initialize
309  _es2->init();
310 
311  // Update the times for interpolation (initially start at 0)
313 
314  // Copy the solutions from the first system
315  for (const auto & var_name : nodal)
316  {
317  _exodusII_io->copy_nodal_solution(*_system, var_name, var_name, _exodus_index1 + 1);
318  _exodusII_io->copy_nodal_solution(*_system2, var_name, var_name, _exodus_index2 + 1);
319  }
320 
321  for (const auto & var_name : elemental)
322  {
323  _exodusII_io->copy_elemental_solution(*_system, var_name, var_name, _exodus_index1 + 1);
324  _exodusII_io->copy_elemental_solution(*_system2, var_name, var_name, _exodus_index2 + 1);
325  }
326 
327  // Update the systems
328  _system->update();
329  _es->update();
330  _system2->update();
331  _es2->update();
332  }
333 
334  // Non-interpolated times
335  else
336  {
337  if (_exodus_time_index > num_exo_times)
338  mooseError("In SolutionUserObject, timestep = ",
340  ", but there are only ",
341  num_exo_times,
342  " time steps.");
343 
344  // Copy the values from the ExodusII file
345  for (const auto & var_name : nodal)
346  _exodusII_io->copy_nodal_solution(*_system, var_name, var_name, _exodus_time_index);
347 
348  for (const auto & var_name : elemental)
349  _exodusII_io->copy_elemental_solution(*_system, var_name, var_name, _exodus_time_index);
350 
351  // Update the equations systems
352  _system->update();
353  _es->update();
354  }
355 }
356 
357 Real
358 SolutionUserObject::directValue(const Node * node, const std::string & var_name) const
359 {
360  // Get the libmesh variable and system numbers
361  unsigned int var_num = _system->variable_number(var_name);
362  unsigned int sys_num = _system->number();
363 
364  // Get the node id and associated dof
365  dof_id_type node_id = node->id();
366  dof_id_type dof_id = _system->get_mesh().node_ref(node_id).dof_number(sys_num, var_num, 0);
367 
368  // Return the desired value for the dof
369  return directValue(dof_id);
370 }
371 
372 Real
373 SolutionUserObject::directValue(const Elem * elem, const std::string & var_name) const
374 {
375  // Get the libmesh variable and system numbers
376  unsigned int var_num = _system->variable_number(var_name);
377  unsigned int sys_num = _system->number();
378 
379  // Get the element id and associated dof
380  dof_id_type elem_id = elem->id();
381  dof_id_type dof_id = _system->get_mesh().elem(elem_id)->dof_number(sys_num, var_num, 0);
382 
383  // Return the desired value
384  return directValue(dof_id);
385 }
386 
387 void
389 {
390 }
391 
392 void
394 {
395 }
396 
397 void
399 {
400  // Update time interpolation for ExodusII solution
401  if (_file_type == 1 && _interpolate_times)
403 }
404 
405 void
407 {
408 }
409 
410 void
412 {
413 
414  // Make sure this only happens once
415  if (_initialized)
416  return;
417 
418  // Several aspects of SolutionUserObject won't work if the FEProblemBase's MooseMesh is
419  // a DistributedMesh:
420  // .) ExodusII_IO::copy_nodal_solution() doesn't work in parallel.
421  // .) We don't know if directValue will be used, which may request
422  // a value on a Node we don't have.
423  _fe_problem.mesh().errorIfDistributedMesh("SolutionUserObject");
424 
425  // Create a libmesh::Mesh object for storing the loaded data. Since
426  // SolutionUserObject is restricted to only work with ReplicatedMesh
427  // (see above) we can force the Mesh used here to be a ReplicatedMesh.
428  _mesh = libmesh_make_unique<ReplicatedMesh>(_communicator);
429 
430  // ExodusII mesh file supplied
431  if (MooseUtils::hasExtension(_mesh_file, "e", /*strip_exodus_ext =*/true))
432  {
433  _file_type = "exodusII";
434  readExodusII();
435  }
436 
437  // XDA mesh file supplied
438  else if (MooseUtils::hasExtension(_mesh_file, "xda"))
439  {
440  _file_type = "xda";
441  readXda();
442  }
443 
444  else if (MooseUtils::hasExtension(_mesh_file, "xdr"))
445  {
446  _file_type = "xdr";
447  readXda();
448  }
449 
450  // Produce an error for an unknown file type
451  else
452  mooseError("In SolutionUserObject, invalid file type (only .xda, .xdr, and .e supported)");
453 
454  // Intilize the serial solution vector
455  _serialized_solution = NumericVector<Number>::build(_communicator);
456  _serialized_solution->init(_system->n_dofs(), false, SERIAL);
457 
458  // Pull down a full copy of this vector on every processor so we can get values in parallel
459  _system->solution->localize(*_serialized_solution);
460 
461  // Vector of variable numbers to apply the MeshFunction to
462  std::vector<unsigned int> var_nums;
463 
464  // If no variables were given, use all of them
465  if (_system_variables.empty())
466  {
467  _system->get_all_variable_numbers(var_nums);
468  for (const auto & var_num : var_nums)
469  _system_variables.push_back(_system->variable_name(var_num));
470  }
471 
472  // Otherwise, gather the numbers for the variables given
473  else
474  {
475  for (const auto & var_name : _system_variables)
476  var_nums.push_back(_system->variable_number(var_name));
477  }
478 
479  // Create the MeshFunction for working with the solution data
480  _mesh_function = libmesh_make_unique<MeshFunction>(
481  *_es, *_serialized_solution, _system->get_dof_map(), var_nums);
482  _mesh_function->init();
483 
484  // Tell the MeshFunctions that we might be querying them outside the
485  // mesh, so we can handle any errors at the MOOSE rather than at the
486  // libMesh level.
487  DenseVector<Number> default_values;
488  _mesh_function->enable_out_of_mesh_mode(default_values);
489 
490  // Build second MeshFunction for interpolation
491  if (_interpolate_times)
492  {
493  // Need to pull down a full copy of this vector on every processor so we can get values in
494  // parallel
495  _serialized_solution2 = NumericVector<Number>::build(_communicator);
496  _serialized_solution2->init(_system2->n_dofs(), false, SERIAL);
497  _system2->solution->localize(*_serialized_solution2);
498 
499  // Create the MeshFunction for the second copy of the data
500  _mesh_function2 = libmesh_make_unique<MeshFunction>(
501  *_es2, *_serialized_solution2, _system2->get_dof_map(), var_nums);
502  _mesh_function2->init();
503  _mesh_function2->enable_out_of_mesh_mode(default_values);
504  }
505 
506  // Populate the data maps that indicate if the variable is nodal and the MeshFunction variable
507  // index
508  for (unsigned int i = 0; i < _system_variables.size(); ++i)
509  {
510  std::string name = _system_variables[i];
511  FEType type = _system->variable_type(name);
512  if (type.order == CONSTANT)
513  _local_variable_nodal[name] = false;
514  else
515  _local_variable_nodal[name] = true;
516 
518  }
519 
520  // Set initialization flag
521  _initialized = true;
522 }
523 
524 MooseEnum
526 {
527  return _file_type;
528 }
529 
530 void
532 {
533  if (time != _interpolation_time)
534  {
536  {
537 
538  for (const auto & var_name : _system_variables)
539  {
540  if (_local_variable_nodal[var_name])
541  _exodusII_io->copy_nodal_solution(*_system, var_name, _exodus_index1 + 1);
542  else
543  _exodusII_io->copy_elemental_solution(*_system, var_name, var_name, _exodus_index1 + 1);
544  }
545 
546  _system->update();
547  _es->update();
548  _system->solution->localize(*_serialized_solution);
549 
550  for (const auto & var_name : _system_variables)
551  {
552  if (_local_variable_nodal[var_name])
553  _exodusII_io->copy_nodal_solution(*_system2, var_name, _exodus_index2 + 1);
554  else
555  _exodusII_io->copy_elemental_solution(*_system2, var_name, var_name, _exodus_index2 + 1);
556  }
557 
558  _system2->update();
559  _es2->update();
560  _system2->solution->localize(*_serialized_solution2);
561  }
562  _interpolation_time = time;
563  }
564 }
565 
566 bool
568 {
569  if (_file_type != 1)
570  mooseError(
571  "In SolutionUserObject, getTimeInterpolationData only applicable for exodusII file type");
572 
573  int old_index1 = _exodus_index1;
574  int old_index2 = _exodus_index2;
575 
576  int num_exo_times = _exodus_times->size();
577 
578  if (time < (*_exodus_times)[0])
579  {
580  _exodus_index1 = 0;
581  _exodus_index2 = 0;
582  _interpolation_factor = 0.0;
583  }
584  else
585  {
586  for (int i = 0; i < num_exo_times - 1; ++i)
587  {
588  if (time <= (*_exodus_times)[i + 1])
589  {
590  _exodus_index1 = i;
591  _exodus_index2 = i + 1;
593  (time - (*_exodus_times)[i]) / ((*_exodus_times)[i + 1] - (*_exodus_times)[i]);
594  break;
595  }
596  else if (i == num_exo_times - 2)
597  {
598  _exodus_index1 = num_exo_times - 1;
599  _exodus_index2 = num_exo_times - 1;
600  _interpolation_factor = 1.0;
601  break;
602  }
603  }
604  }
605 
606  bool indices_modified(false);
607 
608  if (_exodus_index1 != old_index1 || _exodus_index2 != old_index2)
609  indices_modified = true;
610 
611  return indices_modified;
612 }
613 
614 unsigned int
615 SolutionUserObject::getLocalVarIndex(const std::string & var_name) const
616 {
617  // Extract the variable index for the MeshFunction(s)
618  std::map<std::string, unsigned int>::const_iterator it = _local_variable_index.find(var_name);
619  if (it == _local_variable_index.end())
620  mooseError("Value requested for nonexistent variable '",
621  var_name,
622  "' in the '",
623  name(),
624  "' SolutionUserObject");
625  return it->second;
626 }
627 
628 Real
630  const Point & p,
631  const std::string & var_name,
632  const MooseEnum & weighting_type) const
633 {
634  // first check if the FE type is continuous because in that case the value is
635  // unique and we can take a short cut, the default weighting_type found_first also
636  // shortcuts out
637  if (weighting_type == 1 ||
638  (_fe_problem.getVariable(_tid, var_name).feType().family != L2_LAGRANGE &&
639  _fe_problem.getVariable(_tid, var_name).feType().family != MONOMIAL &&
640  _fe_problem.getVariable(_tid, var_name).feType().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 
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 
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 
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 
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 }
const std::string & name() const
Get the name of the object.
Definition: MooseObject.h:47
RealVectorValue RealGradient
Definition: Assembly.h:43
void readXda()
Method for reading XDA mesh and equation systems file(s) This method is called by the constructor whe...
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::unique_ptr< ExodusII_IO > _exodusII_io
Pointer to the libMesh::ExodusII used to read the files.
VectorValue< Real > RealVectorValue
Definition: Assembly.h:40
const FEType & feType() const
Get the type of finite element object.
std::vector< Real > _scale_multiplier
scale_multiplier parameter
std::unique_ptr< NumericVector< Number > > _serialized_solution
Pointer to the serial solution vector.
std::vector< std::string > _system_variables
A list of variables to extract from the read system.
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) ...
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::unique_ptr< EquationSystems > _es
Pointer to the libmesh::EquationSystems object.
void updateExodusTimeInterpolation(Real time)
Updates the times for interpolating ExodusII data.
The main MOOSE class responsible for handling user-defined parameters in almost every MOOSE system...
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 ...
RealVectorValue _rotation1_vector
vector about which to rotate
THREAD_ID _tid
Thread ID of this postprocessor.
Definition: UserObject.h:152
void errorIfDistributedMesh(std::string name) const
Generate a unified error message if the underlying libMesh mesh is a DistributedMesh.
Definition: MooseMesh.C:2473
virtual void finalize() override
Finalize.
InputParameters validParams< GeneralUserObject >()
bool isVariableNodal(const std::string &var_name) const
const std::vector< Real > * _exodus_times
The times available in the ExodusII file.
RealVectorValue _rotation0_vector
vector about which to rotate
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...
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 pointValue(Real t, const Point &p, const unsigned int local_var_index) const
Returns a value at a specific location and variable (see SolutionFunction)
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:121
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".
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 ...
void readExodusII()
Method for reading an ExodusII file, which is called when &#39;file_type = exodusII is set in the input f...
Real directValue(const Node *node, const std::string &var_name) const
Return a value directly from a Node.
std::string _system_name
The system name to extract from the XDA file (xda only)
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:37
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:227
std::unique_ptr< MeshFunction > _mesh_function2
Pointer to second libMesh::MeshFuntion, used for interpolation.
bool isParamValid(const std::string &name) const
Test if the supplied parameter is valid.
Definition: MooseObject.h:67
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.
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. ...
MatType type
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...
bool _initialized
True if initial_setup has executed.
FEProblemBase & _fe_problem
Reference to the FEProblemBase for this user object.
Definition: UserObject.h:149
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
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...
unsigned int size() const
Return the number of items in the MultiMooseEnum.
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...
void mooseError(Args &&...args) const
Definition: MooseObject.h:80
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...
virtual dof_id_type maxElemId() const
Definition: MooseMesh.C:2042
unsigned int getLocalVarIndex(const std::string &var_name) const
Returns the local index for a given variable name.
const std::vector< std::string > & variableNames() const
bool updateExodusBracketingTimeIndices(Real time)
Updates the time indices to interpolate between for ExodusII data.
TensorValue< Real > RealTensorValue
Definition: Assembly.h:45
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...
MooseEnum _file_type
File type to read (0 = xda; 1 = ExodusII)
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.
int _exodus_index2
Time index 2, used for interpolation.
virtual MooseVariable & getVariable(THREAD_ID tid, const std::string &var_name) override
Returns the variable reference for requested variable which may be in any system. ...
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.