SolutionUserObjectBase.C

Go to the documentation of this file.

//* This file is part of the MOOSE framework
//* https://mooseframework.inl.gov
//*
//* All rights reserved, see COPYRIGHT for full restrictions
//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
//*
//* Licensed under LGPL 2.1, please see LICENSE for details
//* https://www.gnu.org/licenses/lgpl-2.1.html
 
#include "SolutionUserObjectBase.h"
 
#include <iterator>
 
// MOOSE includes
#include "ConsoleUtils.h"
#include "MooseError.h"
#include "MooseMesh.h"
#include "MooseUtils.h"
#include "MooseVariableFE.h"
#include "RotationMatrix.h"
#include "Function.h"
 
// libMesh includes
#include "libmesh/equation_systems.h"
#include "libmesh/mesh_function.h"
#include "libmesh/numeric_vector.h"
#include "libmesh/nonlinear_implicit_system.h"
#include "libmesh/transient_system.h"
#include "libmesh/parallel_mesh.h"
#include "libmesh/serial_mesh.h"
#include "libmesh/exodusII_io.h"
#include "libmesh/exodusII_io_helper.h"
#include "libmesh/nemesis_io.h"
#include "libmesh/nemesis_io_helper.h"
#include "libmesh/enum_xdr_mode.h"
#include "libmesh/string_to_enum.h"
#include "libmesh/fe_interface.h"
 
InputParameters
SolutionUserObjectBase::validParams()
{
  // Get the input parameters from the parent class
  InputParameters params = GeneralUserObject::validParams();
 
  // Add required parameters
  params.addRequiredParam<MeshFileName>(
      "mesh", "The name of the mesh file (must be xda/xdr, exodusII or nemesis file).");
  params.addParam<std::vector<std::string>>(
      "system_variables",
      std::vector<std::string>(),
      "The name of the nodal and elemental variables from the file you want to use for values");
  params.addParam<bool>(
      "force_replicated_source_mesh",
      false,
      "Whether to force the serialization of the source mesh. This can be useful for discontinuous "
      "variables or if the partitioning of the source mesh and the current mesh are different");
 
  // When using XDA/XDR files the following must be defined
  params.addParam<FileName>(
      "es",
      "<not supplied>",
      "The name of the file holding the equation system info in xda/xdr format (xda/xdr only).");
  params.addParam<std::string>(
      "system",
      "nl0",
      "The name of the system to pull values out of (xda/xdr only). The default name for the "
      "nonlinear system is 'nl0', auxiliary system is 'aux0'");
 
  // When using ExodusII a specific time is extracted
  params.addParam<std::string>("timestep",
                               "Index of the single timestep used or \"LATEST\" for "
                               "the last timestep (exodusII only).  If not supplied, "
                               "time interpolation will occur.");
 
  // Choose the interpolation order for incoming nodal variables
  params.addParam<MooseEnum>("nodal_variable_order",
                             MooseEnum("FIRST SECOND", "FIRST"),
                             "Specifies the order of the nodal solution data.");
 
  // Add ability to perform coordinate transformation: scale, factor
  params.addParam<std::vector<Real>>(
      "scale", std::vector<Real>(LIBMESH_DIM, 1), "Scale factor for points in the simulation");
  params.addParam<std::vector<Real>>("scale_multiplier",
                                     std::vector<Real>(LIBMESH_DIM, 1),
                                     "Scale multiplying factor for points in the simulation");
  params.addParam<std::vector<Real>>("translation",
                                     std::vector<Real>(LIBMESH_DIM, 0),
                                     "Translation factors for x,y,z coordinates of the simulation");
  params.addParam<RealVectorValue>("rotation0_vector",
                                   RealVectorValue(0, 0, 1),
                                   "Vector about which to rotate points of the simulation.");
  params.addParam<Real>(
      "rotation0_angle",
      0.0,
      "Anticlockwise rotation angle (in degrees) to use for rotation about rotation0_vector.");
  params.addParam<RealVectorValue>("rotation1_vector",
                                   RealVectorValue(0, 0, 1),
                                   "Vector about which to rotate points of the simulation.");
  params.addParam<Real>(
      "rotation1_angle",
      0.0,
      "Anticlockwise rotation angle (in degrees) to use for rotation about rotation1_vector.");
 
  // following lines build the default_transformation_order
  MultiMooseEnum default_transformation_order(
      "rotation0 translation scale rotation1 scale_multiplier", "translation scale");
  params.addParam<MultiMooseEnum>(
      "transformation_order",
      default_transformation_order,
      "The order to perform the operations in.  Define R0 to be the rotation matrix encoded by "
      "rotation0_vector and rotation0_angle.  Similarly for R1.  Denote the scale by s, the "
      "scale_multiplier by m, and the translation by t.  Then, given a point x in the simulation, "
      "if transformation_order = 'rotation0 scale_multiplier translation scale rotation1' then "
      "form p = R1*(R0*x*m - t)/s.  Then the values provided by the SolutionUserObjectBase at "
      "point x "
      "in the simulation are the variable values at point p in the mesh.");
  params.addParamNamesToGroup("scale scale_multiplier translation rotation0_vector rotation0_angle "
                              "rotation1_angle transformation_order",
                              "Coordinate system transformation");
  params.addClassDescription("Reads a variable from a mesh in one simulation to another");
  return params;
}
 
// Static mutex definition
Threads::spin_mutex SolutionUserObjectBase::_solution_user_object_mutex;
 
SolutionUserObjectBase::SolutionUserObjectBase(const InputParameters & parameters)
  : GeneralUserObject(parameters),
    _file_type(MooseEnum("xda=0 exodusII=1 xdr=2 nemesis=3")),
    _mesh_file(getParam<MeshFileName>("mesh")),
    _es_file(getParam<FileName>("es")),
    _system_name(getParam<std::string>("system")),
    _system_variables(getParam<std::vector<std::string>>("system_variables")),
    _exodus_time_index(-1),
    _interpolate_times(false),
    _system(nullptr),
    _system2(nullptr),
    _interpolation_time(0.0),
    _interpolation_factor(0.0),
    _exodus_times(nullptr),
    _exodus_index1(0),
    _exodus_index2(0),
    _nodal_variable_order(getParam<MooseEnum>("nodal_variable_order")),
    _scale(getParam<std::vector<Real>>("scale")),
    _scale_multiplier(getParam<std::vector<Real>>("scale_multiplier")),
    _translation(getParam<std::vector<Real>>("translation")),
    _rotation0_vector(getParam<RealVectorValue>("rotation0_vector")),
    _rotation0_angle(getParam<Real>("rotation0_angle")),
    _r0(RealTensorValue()),
    _rotation1_vector(getParam<RealVectorValue>("rotation1_vector")),
    _rotation1_angle(getParam<Real>("rotation1_angle")),
    _r1(RealTensorValue()),
    _transformation_order(getParam<MultiMooseEnum>("transformation_order")),
    _force_replicated_source(getParam<bool>("force_replicated_source_mesh")),
    _initialized(false)
{
  // form rotation matrices with the specified angles
  Real halfPi = libMesh::pi / 2.0;
  Real a;
  Real b;
 
  a = std::cos(halfPi * -_rotation0_angle / 90.0);
  b = std::sin(halfPi * -_rotation0_angle / 90.0);
  // the following is an anticlockwise rotation about z
  RealTensorValue rot0_z(a, -b, 0, b, a, 0, 0, 0, 1);
  // form the rotation matrix that will take rotation0_vector to the z axis
  RealTensorValue vec0_to_z = RotationMatrix::rotVecToZ(_rotation0_vector);
  // _r0 is then: rotate points so vec0 lies along z; then rotate about angle0; then rotate points
  // back
  _r0 = vec0_to_z.transpose() * (rot0_z * vec0_to_z);
 
  a = std::cos(halfPi * -_rotation1_angle / 90.0);
  b = std::sin(halfPi * -_rotation1_angle / 90.0);
  // the following is an anticlockwise rotation about z
  RealTensorValue rot1_z(a, -b, 0, b, a, 0, 0, 0, 1);
  // form the rotation matrix that will take rotation1_vector to the z axis
  RealTensorValue vec1_to_z = RotationMatrix::rotVecToZ(_rotation1_vector);
  // _r1 is then: rotate points so vec1 lies along z; then rotate about angle1; then rotate points
  // back
  _r1 = vec1_to_z.transpose() * (rot1_z * vec1_to_z);
 
  if (isParamValid("timestep") && getParam<std::string>("timestep") == "-1")
    mooseError("A \"timestep\" of -1 is no longer supported for interpolation. Instead simply "
               "remove this parameter altogether for interpolation");
 
  // check for possible inconsistencies between mesh and input data
  if (_fe_problem.mesh().hasSecondOrderElements())
    mooseInfo(
        "The mesh has second-order elements, be sure to set 'nodal_variable_order' if needed.");
}
 
void
SolutionUserObjectBase::readXda()
{
  if (!isParamSetByUser("es"))
    paramError("es", "Equation system file (.xda or .xdr) should have been specified");
 
  // Check that the required files exist
  MooseUtils::checkFileReadable(_es_file);
  MooseUtils::checkFileReadable(_mesh_file);
 
  // Read the libmesh::mesh from the xda file
  _mesh->read(_mesh_file);
 
  // Create the libmesh::EquationSystems
  _es = std::make_unique<EquationSystems>(*_mesh);
 
  // Use new read syntax (binary)
  if (_file_type == "xdr")
    _es->read(_es_file,
              libMesh::DECODE,
              EquationSystems::READ_HEADER | EquationSystems::READ_DATA |
                  EquationSystems::READ_ADDITIONAL_DATA);
 
  // Use new read syntax
  else if (_file_type == "xda")
    _es->read(_es_file,
              libMesh::READ,
              EquationSystems::READ_HEADER | EquationSystems::READ_DATA |
                  EquationSystems::READ_ADDITIONAL_DATA);
 
  // This should never occur, just in case produce an error
  else
    mooseError("Failed to determine proper read method for XDA/XDR equation system file: ",
               _es_file);
 
  // Update and store the EquationSystems name locally
  _es->update();
  _system = &_es->get_system(_system_name);
}
 
void
SolutionUserObjectBase::readExodusIIOrNemesis()
{
  // Define a default system name
  if (_system_name == "")
    _system_name = "SolutionUserObjectSystem";
 
  // Read the Exodus or Nemesis file
  if (_file_type == "exodusII")
  {
    _exodusII_io = std::make_unique<libMesh::ExodusII_IO>(*_mesh);
    _exodusII_io->read(_mesh_file);
    _exodus_times = &_exodusII_io->get_time_steps();
  }
  else
  {
    _nemesis_io = std::make_unique<libMesh::Nemesis_IO>(*_mesh);
    _nemesis_io->read(_mesh_file);
    _exodus_times = &_nemesis_io->get_time_steps();
  }
  // Same routine for Nemesis
  readBlockIdMapFromExodusII();
 
  if (isParamValid("timestep"))
  {
    std::string s_timestep = getParam<std::string>("timestep");
    int n_steps = (_file_type == "exodusII") ? _exodusII_io->get_num_time_steps()
                                             : _nemesis_io->get_num_time_steps();
    if (s_timestep == "LATEST")
      _exodus_time_index = n_steps;
    else
    {
      std::istringstream ss(s_timestep);
      if (!((ss >> _exodus_time_index) && ss.eof()) || _exodus_time_index > n_steps)
        mooseError("Invalid value passed as \"timestep\". Expected \"LATEST\" or a valid integer "
                   "less than ",
                   n_steps,
                   ", received ",
                   s_timestep);
    }
  }
  else
    // Interpolate between times rather than using values from a set timestep
    _interpolate_times = true;
 
  // Check that the number of time steps is valid
  int num_exo_times = _exodus_times->size();
  if (num_exo_times == 0)
    mooseError("In SolutionUserObjectBase, exodus file contains no timesteps.");
 
  // Account for parallel mesh
  if (_fe_problem.mesh().isDistributedMesh() && (_file_type != "exodusII") &&
      !_force_replicated_source)
  {
    _mesh->allow_renumbering(true);
    _mesh->prepare_for_use();
  }
  else
  {
    _mesh->allow_renumbering(false);
    _mesh->prepare_for_use();
  }
 
  // Create EquationSystems object for solution
  _es = std::make_unique<EquationSystems>(*_mesh);
  _es->add_system<libMesh::ExplicitSystem>(_system_name);
  _system = &_es->get_system(_system_name);
 
  // Get the variable name lists as set; these need to be sets to perform set_intersection
  const std::vector<std::string> & all_nodal((_file_type == "exodusII")
                                                 ? _exodusII_io->get_nodal_var_names()
                                                 : _nemesis_io->get_nodal_var_names());
  const std::vector<std::string> & all_elemental((_file_type == "exodusII")
                                                     ? _exodusII_io->get_elem_var_names()
                                                     : _nemesis_io->get_elem_var_names());
  const std::vector<std::string> & all_scalar((_file_type == "exodusII")
                                                  ? _exodusII_io->get_global_var_names()
                                                  : _nemesis_io->get_global_var_names());
 
  // Build nodal/elemental variable lists, limit to variables listed in 'system_variables', if
  // provided
  // This function could be called more than once, so clear the member variables so we don't keep
  // adding to the vectors
  _nodal_variables.clear();
  _elemental_variables.clear();
  _scalar_variables.clear();
  if (!_system_variables.empty())
  {
    for (const auto & var_name : _system_variables)
    {
      if (std::find(all_nodal.begin(), all_nodal.end(), var_name) != all_nodal.end())
        _nodal_variables.push_back(var_name);
      else if (std::find(all_elemental.begin(), all_elemental.end(), var_name) !=
               all_elemental.end())
        _elemental_variables.push_back(var_name);
      else if (std::find(all_scalar.begin(), all_scalar.end(), var_name) != all_scalar.end())
        // We checked other variables types first, so if a postprocessor has the same name as a
        // field variable, it won't get loaded as a scalar variable
        _scalar_variables.push_back(var_name);
      else
        paramError("system_variables", "Variable '" + var_name + "' was not found in Exodus file");
    }
  }
  else
  {
    _nodal_variables = all_nodal;
    _elemental_variables = all_elemental;
 
    for (auto var_name : all_scalar)
      // Check if the scalar matches any field variables, and ignore the var if it does. This means
      // its a Postprocessor.
      if (std::find(begin(_nodal_variables), end(_nodal_variables), var_name) ==
              _nodal_variables.end() &&
          std::find(begin(_elemental_variables), end(_elemental_variables), var_name) ==
              _elemental_variables.end())
        _scalar_variables.push_back(var_name);
  }
 
  // Add the variables to the system
  for (const auto & var_name : _nodal_variables)
    _system->add_variable(var_name, Utility::string_to_enum<Order>(_nodal_variable_order));
 
  for (const auto & var_name : _elemental_variables)
    _system->add_variable(var_name, CONSTANT, MONOMIAL);
 
  for (const auto & var_name : _scalar_variables)
    _system->add_variable(var_name, FIRST, SCALAR);
 
  // Helper for copying from the ExodusII/Nemesis file
  const auto copy_solutions = [this](auto & io, libMesh::System & system, unsigned int index)
  {
    for (const auto & var_name : _nodal_variables)
      io.copy_nodal_solution(system, var_name, var_name, index);
 
    for (const auto & var_name : _elemental_variables)
      io.copy_elemental_solution(system, var_name, var_name, index);
 
    if (_scalar_variables.size() > 0)
      io.copy_scalar_solution(system, _scalar_variables, _scalar_variables, index);
  };
 
  // Initialize the equations systems
  _es->init();
 
  // Interpolate times
  if (_interpolate_times)
  {
    // Create a second equation system
    _es2 = std::make_unique<EquationSystems>(*_mesh);
    _es2->add_system<libMesh::ExplicitSystem>(_system_name);
    _system2 = &_es2->get_system(_system_name);
 
    // Add the variables to the system
    for (const auto & var_name : _nodal_variables)
      _system2->add_variable(var_name, Utility::string_to_enum<Order>(_nodal_variable_order));
 
    for (const auto & var_name : _elemental_variables)
      _system2->add_variable(var_name, CONSTANT, MONOMIAL);
 
    for (const auto & var_name : _scalar_variables)
      _system2->add_variable(var_name, FIRST, SCALAR);
 
    // Initialize
    _es2->init();
 
    // Update the times for interpolation (initially start at 0)
    updateInterpolationBracketingTimeIndices();
 
    // Copy the solutions from the first system
    if (_file_type == "exodusII")
    {
      copy_solutions(*_exodusII_io, *_system, _exodus_index1 + 1);
      copy_solutions(*_exodusII_io, *_system2, _exodus_index2 + 1);
    }
    else
    {
      copy_solutions(*_nemesis_io, *_system, _exodus_index1 + 1);
      copy_solutions(*_nemesis_io, *_system2, _exodus_index2 + 1);
    }
 
    // Update the systems
    _system->update();
    _es->update();
    _system2->update();
    _es2->update();
  }
 
  // Non-interpolated times
  else
  {
    if (_exodus_time_index > num_exo_times)
      mooseError("In SolutionUserObjectBase, timestep = ",
                 _exodus_time_index,
                 ", but there are only ",
                 num_exo_times,
                 " time steps.");
 
    // Copy the values from the ExodusII/Nemsis file
    if (_file_type == "exodusII")
      copy_solutions(*_exodusII_io, *_system, _exodus_time_index);
    else
      copy_solutions(*_nemesis_io, *_system, _exodus_time_index);
 
    // Update the equations systems
    _system->update();
    _es->update();
  }
}
 
Real
SolutionUserObjectBase::directValue(const Node * node, const std::string & var_name) const
{
  // Get the libmesh variable and system numbers
  unsigned int var_num = _system->variable_number(var_name);
  unsigned int sys_num = _system->number();
 
  // Get the node id and associated dof
  dof_id_type node_id = node->id();
  const Node & sys_node = _system->get_mesh().node_ref(node_id);
  mooseAssert(sys_node.n_dofs(sys_num, var_num) > 0,
              "Variable " << var_name << " has no DoFs on node " << sys_node.id());
  dof_id_type dof_id = sys_node.dof_number(sys_num, var_num, 0);
 
  // Return the desired value for the dof
  return directValue(dof_id);
}
 
Real
SolutionUserObjectBase::directValue(const Elem * elem, const std::string & var_name) const
{
  // Get the libmesh variable and system numbers
  unsigned int var_num = _system->variable_number(var_name);
  unsigned int sys_num = _system->number();
 
  // Get the element id and associated dof
  dof_id_type elem_id = elem->id();
  const Elem & sys_elem = _system->get_mesh().elem_ref(elem_id);
  mooseAssert(sys_elem.n_dofs(sys_num, var_num) > 0,
              "Variable " << var_name << " has no DoFs on element " << sys_elem.id());
  dof_id_type dof_id = sys_elem.dof_number(sys_num, var_num, 0);
 
  // Return the desired value
  return directValue(dof_id);
}
 
void
SolutionUserObjectBase::initialize()
{
}
 
void
SolutionUserObjectBase::finalize()
{
}
 
void
SolutionUserObjectBase::timestepSetup()
{
  // Update time interpolation for ExodusII or Nemesis solution
  if (_interpolate_times)
    updateTimeInterpolationFromFile();
 
  // Clear the caches
  _cached_p(0) = std::numeric_limits<Real>::max();
  _cached_p2(0) = std::numeric_limits<Real>::max();
}
 
void
SolutionUserObjectBase::execute()
{
}
 
void
SolutionUserObjectBase::initialSetup()
{
 
  // Make sure this only happens once
  if (_initialized)
    return;
 
  // We need to detect exodus files early
  const bool has_exodus_extension =
      MooseUtils::hasExtension(_mesh_file, "e", /*strip_exodus_ext =*/true) ||
      MooseUtils::hasExtension(_mesh_file, "exo", true);
  const bool has_nemesis_extension =
      MooseUtils::hasExtension(_mesh_file, "n", /*strip_exodus_ext =*/true) ||
      MooseUtils::hasExtension(_mesh_file, "nem", true);
  if (has_exodus_extension)
    _file_type = "exodusII";
  else if (has_nemesis_extension)
    _file_type = "nemesis";
 
  // Create a libmesh::Mesh object for storing the loaded data.
  if (!_fe_problem.mesh().isDistributedMesh() || (_file_type == "exodusII") ||
      _force_replicated_source)
    _mesh = std::make_unique<ReplicatedMesh>(_communicator);
  else
    _mesh = std::make_unique<DistributedMesh>(_communicator);
 
  // ExodusII or Nemesis mesh file supplied
  if (has_exodus_extension || has_nemesis_extension)
    readExodusIIOrNemesis();
 
  // XDA mesh file supplied
  else if (MooseUtils::hasExtension(_mesh_file, "xda"))
  {
    _file_type = "xda";
    readXda();
  }
 
  // XDR mesh file supplied
  else if (MooseUtils::hasExtension(_mesh_file, "xdr"))
  {
    _file_type = "xdr";
    readXda();
  }
 
  // Produce an error for an unknown file type
  else
    mooseError(
        "In SolutionUserObjectBase, invalid file type: only .xda, .xdr, .exo and .e supported");
 
  // Initialize the serial solution vector
  _serialized_solution = NumericVector<Number>::build(_communicator);
  _serialized_solution->init(_system->n_dofs(), false, libMesh::SERIAL);
 
  // Pull down a full copy of this vector on every processor so we can get values in parallel
  _system->solution->localize(*_serialized_solution);
 
  // Vector of variable numbers to apply the MeshFunction to
  std::vector<unsigned int> var_nums;
 
  // If no variables were given, use all of them
  if (_system_variables.empty())
  {
    _system->get_all_variable_numbers(var_nums);
    for (const auto & var_num : var_nums)
      _system_variables.push_back(_system->variable_name(var_num));
  }
 
  // Otherwise, gather the numbers for the variables given
  else
  {
    for (const auto & var_name : _system_variables)
      var_nums.push_back(_system->variable_number(var_name));
  }
 
  // Create the MeshFunction for working with the solution data
  _mesh_function = std::make_unique<libMesh::MeshFunction>(
      *_es, *_serialized_solution, _system->get_dof_map(), var_nums);
  _mesh_function->init();
 
  // Tell the MeshFunctions that we might be querying them outside the
  // mesh, so we can handle any errors at the MOOSE rather than at the
  // libMesh level.
  DenseVector<Number> default_values;
  _mesh_function->enable_out_of_mesh_mode(default_values);
 
  // Build second MeshFunction for interpolation
  if (_interpolate_times)
  {
    // Need to pull down a full copy of this vector on every processor so we can get values in
    // parallel
    _serialized_solution2 = NumericVector<Number>::build(_communicator);
    _serialized_solution2->init(_system2->n_dofs(), false, libMesh::SERIAL);
    _system2->solution->localize(*_serialized_solution2);
 
    // Create the MeshFunction for the second copy of the data
    _mesh_function2 = std::make_unique<libMesh::MeshFunction>(
        *_es2, *_serialized_solution2, _system2->get_dof_map(), var_nums);
    _mesh_function2->init();
    _mesh_function2->enable_out_of_mesh_mode(default_values);
  }
 
  // Populate the MeshFunction variable index
  for (unsigned int i = 0; i < _system_variables.size(); ++i)
  {
    std::string name = _system_variables[i];
    _local_variable_index[name] = i;
  }
 
  // If the start time is not the same as in the exodus file, we may need this on INITIAL
  if (_interpolate_times)
    updateTimeInterpolationFromFile();
 
  // Set initialization flag
  _initialized = true;
}
 
bool
SolutionUserObjectBase::isVariableADiscontinuousScalarField(const std::string & var_name) const
{
  const auto & fe_type = _system->variable_type(var_name);
 
  return FEInterface::field_type(fe_type) == libMesh::TYPE_SCALAR &&
         fe_type.family != libMesh::SCALAR &&
         FEInterface::get_continuity(fe_type) == libMesh::DISCONTINUOUS;
}
 
bool
SolutionUserObjectBase::isVariableScalarValued(const std::string & var_name) const
{
  return isVariableScalarValued(getLocalVarIndex(var_name));
}
 
bool
SolutionUserObjectBase::isVariableScalarValued(const unsigned int local_var_index) const
{
  mooseAssert(local_var_index < _system_variables.size(),
              "The local variable index is outside the range of imported solution variables.");
 
  const auto & fe_type = _system->variable_type(_system_variables[local_var_index]);
  return FEInterface::field_type(fe_type) == libMesh::TYPE_SCALAR;
}
 
MooseEnum
SolutionUserObjectBase::getSolutionFileType() const
{
  return _file_type;
}
 
void
SolutionUserObjectBase::updateTimeInterpolationFromFile()
{
  if (_t != _interpolation_time)
  {
    if (updateInterpolationBracketingTimeIndices())
    {
      // Helper to copy/update solutions from a file IO to a system
      const auto update_solutions = [this](auto & io, libMesh::System & system, unsigned int index)
      {
        for (const auto & var_name : _nodal_variables)
          io.copy_nodal_solution(system, var_name, var_name, index);
        for (const auto & var_name : _elemental_variables)
          io.copy_elemental_solution(system, var_name, var_name, index);
        if (_scalar_variables.size() > 0)
          io.copy_scalar_solution(system, _scalar_variables, _scalar_variables, index);
      };
 
      if (_file_type == "exodusII")
        update_solutions(*_exodusII_io, *_system, _exodus_index1 + 1);
      else
        update_solutions(*_nemesis_io, *_system, _exodus_index1 + 1);
 
      _system->update();
      _es->update();
      _system->solution->localize(*_serialized_solution);
 
      if (_file_type == "exodusII")
        update_solutions(*_exodusII_io, *_system2, _exodus_index2 + 1);
      else
        update_solutions(*_nemesis_io, *_system2, _exodus_index2 + 1);
 
      _system2->update();
      _es2->update();
      _system2->solution->localize(*_serialized_solution2);
    }
    _interpolation_time = _t;
  }
}
 
bool
SolutionUserObjectBase::updateInterpolationBracketingTimeIndices()
{
  if (_file_type != "exodusII" && _file_type != "nemesis")
    mooseError("getTimeInterpolationData only applicable for exodusII or Nemesis file type");
 
  int old_index1 = _exodus_index1;
  int old_index2 = _exodus_index2;
 
  int num_exo_times = _exodus_times->size();
 
  const auto solution_time = solutionSampleTime();
 
  if (solution_time < (*_exodus_times)[0])
  {
    _exodus_index1 = 0;
    _exodus_index2 = 0;
    _interpolation_factor = 0.0;
  }
  else
  {
    for (int i = 0; i < num_exo_times - 1; ++i)
    {
      if (solution_time <= (*_exodus_times)[i + 1])
      {
        _exodus_index1 = i;
        _exodus_index2 = i + 1;
        _interpolation_factor =
            (solution_time - (*_exodus_times)[i]) / ((*_exodus_times)[i + 1] - (*_exodus_times)[i]);
        break;
      }
      else if (i == num_exo_times - 2)
      {
        _exodus_index1 = num_exo_times - 1;
        _exodus_index2 = num_exo_times - 1;
        _interpolation_factor = 1.0;
        break;
      }
    }
  }
 
  bool indices_modified(false);
 
  if (_exodus_index1 != old_index1 || _exodus_index2 != old_index2)
    indices_modified = true;
 
  return indices_modified;
}
 
unsigned int
SolutionUserObjectBase::getLocalVarIndex(const std::string & var_name) const
{
  // Extract the variable index for the MeshFunction(s)
  std::map<std::string, unsigned int>::const_iterator it = _local_variable_index.find(var_name);
  if (it == _local_variable_index.end())
    mooseError("Value requested for nonexistent variable '",
               var_name,
               "' in the '",
               name(),
               "' SolutionUserObjectBase.\nSystem selected: ",
               _system_name,
               "\nAvailable variables:\n",
               ConsoleUtils::formatString(Moose::stringify(_system_variables), ""));
  return it->second;
}
 
Real
SolutionUserObjectBase::pointValue(Real t,
                                   const Point & p,
                                   const std::string & var_name,
                                   WeightingType weighting_type,
                                   const std::set<subdomain_id_type> * subdomain_ids) const
{
  // Use multivalued evaluation only for spatial discontinuous scalar fields.
  const auto & fe_type = _system->variable_type(var_name);
  const auto continuity = FEInterface::get_continuity(fe_type);
  const auto field_type = FEInterface::field_type(fe_type);
 
  if (field_type != libMesh::TYPE_SCALAR || continuity != libMesh::DISCONTINUOUS ||
      fe_type.family == libMesh::SCALAR)
    return pointValue(t, p, var_name, subdomain_ids);
 
  // the shape function is discontinuous so we need to compute a suitable unique value
  std::map<const Elem *, Real> values = discontinuousPointValue(t, p, var_name, subdomain_ids);
 
  mooseAssert(!values.empty(),
              "discontinuousPointValue() should error rather than return an empty map.");
 
  switch (weighting_type)
  {
    case WeightingType::AVERAGE:
    {
      Real average = 0.0;
      for (auto & v : values)
        average += v.second;
      return average / Real(values.size());
    }
  }
 
  mooseError("SolutionUserObjectBase::pointValue reaches line that it should not be able to "
             "reach.");
  return 0.0;
}
 
Real
SolutionUserObjectBase::pointValue(Real t,
                                   const Point & p,
                                   const std::string & var_name,
                                   const std::set<subdomain_id_type> * subdomain_ids) const
{
  const unsigned int local_var_index = getLocalVarIndex(var_name);
  return pointValue(t, p, local_var_index, subdomain_ids);
}
 
Real
SolutionUserObjectBase::pointValue(Real libmesh_dbg_var(t),
                                   const Point & p,
                                   const unsigned int local_var_index,
                                   const std::set<subdomain_id_type> * subdomain_ids) const
{
  mooseAssert(isVariableScalarValued(local_var_index),
              "Point value evaluation requires a scalar-valued imported variable.");
 
  // Create copy of point
  Point pt(p);
 
  // do the transformations
  for (unsigned int trans_num = 0; trans_num < _transformation_order.size(); ++trans_num)
  {
    if (_transformation_order[trans_num] == "rotation0")
      pt = _r0 * pt;
    else if (_transformation_order[trans_num] == "translation")
      for (const auto i : make_range(Moose::dim))
        pt(i) -= _translation[i];
    else if (_transformation_order[trans_num] == "scale")
      for (const auto i : make_range(Moose::dim))
        pt(i) /= _scale[i];
    else if (_transformation_order[trans_num] == "scale_multiplier")
      for (const auto i : make_range(Moose::dim))
        pt(i) *= _scale_multiplier[i];
    else if (_transformation_order[trans_num] == "rotation1")
      pt = _r1 * pt;
  }
 
  // Extract the value at the current point
  Real val = evalMeshFunction(pt, local_var_index, 1, subdomain_ids);
 
  // Interpolate
  if (_interpolate_times && (_file_type == "exodusII" || _file_type == "nemesis"))
  {
    mooseAssert(t == _interpolation_time,
                "Time passed into value() must match time at last call to timestepSetup()");
    Real val2 = evalMeshFunction(pt, local_var_index, 2, subdomain_ids);
    val = val + (val2 - val) * _interpolation_factor;
  }
 
  return val;
}
 
std::map<const Elem *, Real>
SolutionUserObjectBase::discontinuousPointValue(
    Real t,
    const Point & p,
    const std::string & var_name,
    const std::set<subdomain_id_type> * subdomain_ids) const
{
  const unsigned int local_var_index = getLocalVarIndex(var_name);
  return discontinuousPointValue(t, p, local_var_index, subdomain_ids);
}
 
std::map<const Elem *, Real>
SolutionUserObjectBase::discontinuousPointValue(
    Real libmesh_dbg_var(t),
    Point pt,
    const unsigned int local_var_index,
    const std::set<subdomain_id_type> * subdomain_ids) const
{
  mooseAssert(isVariableScalarValued(local_var_index),
              "Point value evaluation requires a scalar-valued imported variable.");
 
  // do the transformations
  for (unsigned int trans_num = 0; trans_num < _transformation_order.size(); ++trans_num)
  {
    if (_transformation_order[trans_num] == "rotation0")
      pt = _r0 * pt;
    else if (_transformation_order[trans_num] == "translation")
      for (const auto i : make_range(Moose::dim))
        pt(i) -= _translation[i];
    else if (_transformation_order[trans_num] == "scale")
      for (const auto i : make_range(Moose::dim))
        pt(i) /= _scale[i];
    else if (_transformation_order[trans_num] == "scale_multiplier")
      for (const auto i : make_range(Moose::dim))
        pt(i) *= _scale_multiplier[i];
    else if (_transformation_order[trans_num] == "rotation1")
      pt = _r1 * pt;
  }
 
  // Extract the value at the current point
  std::map<const Elem *, Real> map =
      evalMultiValuedMeshFunction(pt, local_var_index, 1, subdomain_ids);
 
  // Interpolate
  if (_interpolate_times && (_file_type == "exodusII" || _file_type == "nemesis"))
  {
    mooseAssert(t == _interpolation_time,
                "Time passed into value() must match time at last call to timestepSetup()");
    std::map<const Elem *, Real> map2 =
        evalMultiValuedMeshFunction(pt, local_var_index, 2, subdomain_ids);
 
    if (map.size() != map2.size())
      mooseError(
          "In SolutionUserObjectBase::discontinuousPointValue map and map2 have different size");
 
    // construct the interpolated map
    for (auto & k : map)
    {
      if (map2.find(k.first) == map2.end())
        mooseError(
            "In SolutionUserObjectBase::discontinuousPointValue map and map2 have differing keys");
      Real val = k.second;
      Real val2 = map2[k.first];
      map[k.first] = val + (val2 - val) * _interpolation_factor;
    }
  }
 
  return map;
}
 
RealGradient
SolutionUserObjectBase::pointValueGradient(Real t,
                                           const Point & p,
                                           const std::string & var_name,
                                           WeightingType weighting_type,
                                           const std::set<subdomain_id_type> * subdomain_ids) const
{
 
  // the shape function is discontinuous so we need to compute a suitable unique value
  std::map<const Elem *, RealGradient> values =
      discontinuousPointValueGradient(t, p, var_name, subdomain_ids);
 
  mooseAssert(!values.empty(),
              "discontinuousPointValueGradient() should error rather than return an empty map.");
 
  switch (weighting_type)
  {
    case WeightingType::AVERAGE:
    {
      RealGradient average = RealGradient(0.0, 0.0, 0.0);
      for (auto & v : values)
        average += v.second;
      return average / Real(values.size());
    }
  }
 
  mooseError("SolutionUserObjectBase::pointValueGradient reaches line that it should not be "
             "able to reach.");
  return RealGradient(0.0, 0.0, 0.0);
}
 
RealGradient
SolutionUserObjectBase::pointValueGradient(Real t,
                                           const Point & p,
                                           const std::string & var_name,
                                           const std::set<subdomain_id_type> * subdomain_ids) const
{
  const unsigned int local_var_index = getLocalVarIndex(var_name);
  return pointValueGradient(t, p, local_var_index, subdomain_ids);
}
 
RealGradient
SolutionUserObjectBase::pointValueGradient(Real libmesh_dbg_var(t),
                                           Point pt,
                                           const unsigned int local_var_index,
                                           const std::set<subdomain_id_type> * subdomain_ids) const
{
  mooseAssert(isVariableScalarValued(local_var_index),
              "Point gradient evaluation requires a scalar-valued imported variable.");
 
  // do the transformations
  for (unsigned int trans_num = 0; trans_num < _transformation_order.size(); ++trans_num)
  {
    if (_transformation_order[trans_num] == "rotation0")
      pt = _r0 * pt;
    else if (_transformation_order[trans_num] == "translation")
      for (const auto i : make_range(Moose::dim))
        pt(i) -= _translation[i];
    else if (_transformation_order[trans_num] == "scale")
      for (const auto i : make_range(Moose::dim))
        pt(i) /= _scale[i];
    else if (_transformation_order[trans_num] == "scale_multiplier")
      for (const auto i : make_range(Moose::dim))
        pt(i) *= _scale_multiplier[i];
    else if (_transformation_order[trans_num] == "rotation1")
      pt = _r1 * pt;
  }
 
  // Extract the value at the current point
  RealGradient val = evalMeshFunctionGradient(pt, local_var_index, 1, subdomain_ids);
 
  // Interpolate
  if (_interpolate_times && (_file_type == "exodusII" || _file_type == "nemesis"))
  {
    mooseAssert(t == _interpolation_time,
                "Time passed into value() must match time at last call to timestepSetup()");
    RealGradient val2 = evalMeshFunctionGradient(pt, local_var_index, 2, subdomain_ids);
    val = val + (val2 - val) * _interpolation_factor;
  }
 
  return val;
}
 
std::map<const Elem *, RealGradient>
SolutionUserObjectBase::discontinuousPointValueGradient(
    Real t,
    const Point & p,
    const std::string & var_name,
    const std::set<subdomain_id_type> * subdomain_ids) const
{
  const unsigned int local_var_index = getLocalVarIndex(var_name);
  return discontinuousPointValueGradient(t, p, local_var_index, subdomain_ids);
}
 
std::map<const Elem *, RealGradient>
SolutionUserObjectBase::discontinuousPointValueGradient(
    Real libmesh_dbg_var(t),
    Point pt,
    const unsigned int local_var_index,
    const std::set<subdomain_id_type> * subdomain_ids) const
{
  mooseAssert(isVariableScalarValued(local_var_index),
              "Point gradient evaluation requires a scalar-valued imported variable.");
 
  // do the transformations
  for (unsigned int trans_num = 0; trans_num < _transformation_order.size(); ++trans_num)
  {
    if (_transformation_order[trans_num] == "rotation0")
      pt = _r0 * pt;
    else if (_transformation_order[trans_num] == "translation")
      for (const auto i : make_range(Moose::dim))
        pt(i) -= _translation[i];
    else if (_transformation_order[trans_num] == "scale")
      for (const auto i : make_range(Moose::dim))
        pt(i) /= _scale[i];
    else if (_transformation_order[trans_num] == "scale_multiplier")
      for (const auto i : make_range(Moose::dim))
        pt(i) *= _scale_multiplier[i];
    else if (_transformation_order[trans_num] == "rotation1")
      pt = _r1 * pt;
  }
 
  // Extract the value at the current point
  std::map<const Elem *, RealGradient> map =
      evalMultiValuedMeshFunctionGradient(pt, local_var_index, 1, subdomain_ids);
 
  // Interpolate
  if (_interpolate_times && (_file_type == "exodusII" || _file_type == "nemesis"))
  {
    mooseAssert(t == _interpolation_time,
                "Time passed into value() must match time at last call to timestepSetup()");
    std::map<const Elem *, RealGradient> map2 =
        evalMultiValuedMeshFunctionGradient(pt, local_var_index, 2, subdomain_ids);
 
    if (map.size() != map2.size())
      mooseError("In SolutionUserObjectBase::discontinuousPointValueGradient map and map2 have "
                 "different size");
 
    // construct the interpolated map
    for (auto & k : map)
    {
      if (map2.find(k.first) == map2.end())
        mooseError("In SolutionUserObjectBase::discontinuousPointValueGradient map and map2 have "
                   "differing keys");
      RealGradient val = k.second;
      RealGradient val2 = map2[k.first];
      map[k.first] = val + (val2 - val) * _interpolation_factor;
    }
  }
 
  return map;
}
 
Real
SolutionUserObjectBase::directValue(dof_id_type dof_index) const
{
  Real val = (*_serialized_solution)(dof_index);
  if (_interpolate_times && (_file_type == "exodusII" || _file_type == "nemesis"))
  {
    Real val2 = (*_serialized_solution2)(dof_index);
    val = val + (val2 - val) * _interpolation_factor;
  }
  return val;
}
 
Real
SolutionUserObjectBase::evalMeshFunction(const Point & p,
                                         const unsigned int local_var_index,
                                         unsigned int func_num,
                                         const std::set<subdomain_id_type> * subdomain_ids) const
{
  // Storage for mesh function output
  DenseVector<Number> output;
 
  // Extract a value from the _mesh_function
  {
    Threads::spin_mutex::scoped_lock lock(_solution_user_object_mutex);
 
    if (func_num == 1)
    {
      // Check the cache
      if (p == _cached_p && (!subdomain_ids || (*subdomain_ids == _cached_subdomain_ids)))
        return _cached_values(local_var_index);
 
      // else get a new value
      (*_mesh_function)(p, 0.0, output, subdomain_ids);
      // and cache it
      _cached_p = p;
      if (subdomain_ids)
        _cached_subdomain_ids = *subdomain_ids;
      else
        _cached_subdomain_ids.clear();
      _cached_values = output;
    }
 
    // Extract a value from _mesh_function2
    else if (func_num == 2)
    {
      // Check the cache
      if (p == _cached_p2 && (!subdomain_ids || (*subdomain_ids == _cached_subdomain_ids2)))
        return _cached_values2(local_var_index);
 
      // else get a new value
      (*_mesh_function2)(p, 0.0, output, subdomain_ids);
      // and cache it
      _cached_p2 = p;
      if (subdomain_ids)
        _cached_subdomain_ids2 = *subdomain_ids;
      else
        _cached_subdomain_ids2.clear();
      _cached_values2 = output;
    }
    else
      mooseError("The func_num must be 1 or 2");
  }
 
  // Error if the data is out-of-range, which will be the case if the mesh functions are evaluated
  // outside the domain
  if (output.size() == 0)
  {
    std::ostringstream oss;
    p.print(oss);
    mooseError("Failed to access the data for variable '",
               _system_variables[local_var_index],
               "' at point ",
               oss.str(),
               " in the '",
               name(),
               "' SolutionUserObjectBase");
  }
  return output(local_var_index);
}
 
std::map<const Elem *, Real>
SolutionUserObjectBase::evalMultiValuedMeshFunction(
    const Point & p,
    const unsigned int local_var_index,
    unsigned int func_num,
    const std::set<subdomain_id_type> * subdomain_ids) const
{
  // Storage for mesh function output
  std::map<const Elem *, DenseVector<Number>> temporary_output;
 
  // Extract a value from the _mesh_function
  {
    Threads::spin_mutex::scoped_lock lock(_solution_user_object_mutex);
    if (func_num == 1)
      _mesh_function->discontinuous_value(p, 0.0, temporary_output, subdomain_ids);
 
    // Extract a value from _mesh_function2
    else if (func_num == 2)
      _mesh_function2->discontinuous_value(p, 0.0, temporary_output, subdomain_ids);
 
    else
      mooseError("The func_num must be 1 or 2");
  }
 
  // Error if the data is out-of-range, which will be the case if the mesh functions are evaluated
  // outside the domain
  if (temporary_output.size() == 0)
  {
    std::ostringstream oss;
    p.print(oss);
    mooseError("Failed to access the data for variable '",
               _system_variables[local_var_index],
               "' at point ",
               oss.str(),
               " in the '",
               name(),
               "' SolutionUserObjectBase");
  }
 
  // Fill the actual map that is returned
  std::map<const Elem *, Real> output;
  for (auto & k : temporary_output)
  {
    mooseAssert(
        k.second.size() > local_var_index,
        "In SolutionUserObjectBase::evalMultiValuedMeshFunction variable with local_var_index "
            << local_var_index << " does not exist");
    output[k.first] = k.second(local_var_index);
  }
 
  return output;
}
 
RealGradient
SolutionUserObjectBase::evalMeshFunctionGradient(
    const Point & p,
    const unsigned int local_var_index,
    unsigned int func_num,
    const std::set<subdomain_id_type> * subdomain_ids) const
{
  // Storage for mesh function output
  std::vector<Gradient> output;
 
  // Extract a value from the _mesh_function
  {
    Threads::spin_mutex::scoped_lock lock(_solution_user_object_mutex);
    if (func_num == 1)
      _mesh_function->gradient(p, 0.0, output, subdomain_ids);
 
    // Extract a value from _mesh_function2
    else if (func_num == 2)
      _mesh_function2->gradient(p, 0.0, output, subdomain_ids);
 
    else
      mooseError("The func_num must be 1 or 2");
  }
 
  // Error if the data is out-of-range, which will be the case if the mesh functions are evaluated
  // outside the domain
  if (output.size() == 0)
  {
    std::ostringstream oss;
    p.print(oss);
    mooseError("Failed to access the data for variable '",
               _system_variables[local_var_index],
               "' at point ",
               oss.str(),
               " in the '",
               name(),
               "' SolutionUserObjectBase");
  }
  return output[local_var_index];
}
 
std::map<const Elem *, RealGradient>
SolutionUserObjectBase::evalMultiValuedMeshFunctionGradient(
    const Point & p,
    const unsigned int local_var_index,
    unsigned int func_num,
    const std::set<subdomain_id_type> * subdomain_ids) const
{
  // Storage for mesh function output
  std::map<const Elem *, std::vector<Gradient>> temporary_output;
 
  // Extract a value from the _mesh_function
  {
    Threads::spin_mutex::scoped_lock lock(_solution_user_object_mutex);
    if (func_num == 1)
      _mesh_function->discontinuous_gradient(p, 0.0, temporary_output, subdomain_ids);
 
    // Extract a value from _mesh_function2
    else if (func_num == 2)
      _mesh_function2->discontinuous_gradient(p, 0.0, temporary_output, subdomain_ids);
 
    else
      mooseError("The func_num must be 1 or 2");
  }
 
  // Error if the data is out-of-range, which will be the case if the mesh functions are evaluated
  // outside the domain
  if (temporary_output.size() == 0)
  {
    std::ostringstream oss;
    p.print(oss);
    mooseError("Failed to access the data for variable '",
               _system_variables[local_var_index],
               "' at point ",
               oss.str(),
               " in the '",
               name(),
               "' SolutionUserObjectBase");
  }
 
  // Fill the actual map that is returned
  std::map<const Elem *, RealGradient> output;
  for (auto & k : temporary_output)
  {
    mooseAssert(
        k.second.size() > local_var_index,
        "In SolutionUserObjectBase::evalMultiValuedMeshFunction variable with local_var_index "
            << local_var_index << " does not exist");
    output[k.first] = k.second[local_var_index];
  }
 
  return output;
}
 
const std::vector<std::string> &
SolutionUserObjectBase::variableNames() const
{
  return _system_variables;
}
 
bool
SolutionUserObjectBase::isVariableNodal(const std::string & var_name) const
{
  return std::find(_nodal_variables.begin(), _nodal_variables.end(), var_name) !=
         _nodal_variables.end();
}
 
Real
SolutionUserObjectBase::scalarValue(Real /*t*/, const std::string & var_name) const
{
  unsigned int var_num = _system->variable_number(var_name);
  const DofMap & dof_map = _system->get_dof_map();
  std::vector<dof_id_type> dofs;
  dof_map.SCALAR_dof_indices(dofs, var_num);
  // We can handle only FIRST order scalar variables
  return directValue(dofs[0]);
}
 
void
SolutionUserObjectBase::readBlockIdMapFromExodusII()
{
#ifdef LIBMESH_HAVE_EXODUS_API
  const auto * exio_helper = (_file_type == "exodusII") ? &_exodusII_io->get_exio_helper()
                                                        : &_nemesis_io->get_nemio_helper();
  const auto & id_to_block = exio_helper->id_to_block_names;
  _block_name_to_id.clear();
  _block_id_to_name.clear();
  for (const auto & it : id_to_block)
  {
    _block_name_to_id[it.second] = it.first;
    _block_id_to_name[it.first] = it.second;
  }
#endif
}