virtual Real computeOutOfPlaneGradDispOld() override;

  /// gets its subblock index for current element
  unsigned int getCurrentSubblockIndex() const
  {
    return _subblock_id_provider ? _subblock_id_provider->getSubblockIndex(*_current_elem) : 0;
  };

  /// A Userobject that carries the subblock ID for all elements

  virtual Real computeOutOfPlaneGradDispOld() override;

  /// gets its subblock index for current element
  unsigned int getCurrentSubblockIndex() const
  {
    return _subblock_id_provider ? _subblock_id_provider->getSubblockIndex(*_current_elem) : 0;
  };

  /// A Userobject that carries the subblock ID for all elements

  virtual ADReal computeOutOfPlaneStrain();

  /// gets its subblock index for current element
  unsigned int getCurrentSubblockIndex() const
  {
    return _subblock_id_provider ? _subblock_id_provider->getSubblockIndex(*_current_elem) : 0;
  };

  /// A Userobject that carries the subblock ID for all elements

  params.addParam<std::vector<TagName>>("absolute_value_vector_tags",
                                        "The tag names for extra vectors that the absolute value "
                                        "of the residual should be accumulated into");
  params.addParam<bool>("use_automatic_differentiation",
                        false,
                        "Use automatic differentiation to assemble the generalized plane strain "
                        "equation and its coupling terms");

  : Action(params),
    _displacements(getParam<std::vector<VariableName>>("displacements")),
    _ndisp(_displacements.size()),
    _out_of_plane_direction(getParam<MooseEnum>("out_of_plane_direction")),
    _use_ad(getParam<bool>("use_automatic_differentiation"))
{
}

unsigned int
GeneralizedPlaneStrainAction::inPlaneDisplacementIndex() const
{
  for (unsigned int i = 0; i < _ndisp; ++i)
    if (i != _out_of_plane_direction)
      return i;

  paramError("displacements", "No in-plane displacement is available to anchor the action");
}

void

  // user object name
  const std::string uo_name = _name + "_GeneralizedPlaneStrainUserObject";

  if (_current_task == "add_variables_physics")
  {
    if (_use_ad)
    {
      std::set<SubdomainID> block_ids;
      if (isParamValid("block"))
        for (const auto & block : getParam<std::vector<SubdomainName>>("block"))
        {
          const auto id = _mesh->getSubdomainID(block);
          if (id == Moose::INVALID_BLOCK_ID)
            paramError("block", "Subdomain '", block, "' was not found in the mesh");
          block_ids.insert(id);
        }

      const auto & subdomains = block_ids.empty() ? _problem->mesh().meshSubdomains() : block_ids;
      if (subdomains.empty())
        mooseError("No subdomains found for the generalized plane strain action");

      const auto coord_system = _problem->getCoordSystem(*subdomains.begin());
      for (const auto subdomain : subdomains)
        if (_problem->getCoordSystem(subdomain) != coord_system)
          paramError("block",
                     "Generalized plane strain requires all selected subdomains to use the same "
                     "coordinate system");

      if (coord_system == Moose::COORD_RZ)
      {
        if (_ndisp != 1)
          paramError("displacements",
                     "One radial displacement is required for 1D axisymmetric generalized plane "
                     "strain");
      }
      else if (coord_system == Moose::COORD_XYZ)
      {
        const unsigned int required_displacements = _out_of_plane_direction == 2 ? 2 : 3;
        if (_ndisp != required_displacements)
          paramError("displacements",
                     required_displacements,
                     " displacement variables are required when the out-of-plane direction is ",
                     getParam<MooseEnum>("out_of_plane_direction"));
      }
      else
        paramError("out_of_plane_direction",
                   "Generalized plane strain supports only Cartesian and axisymmetric coordinate "
                   "systems");
    }

    const auto anchor_displacement = inPlaneDisplacementIndex();
    const auto & anchor_variable = _problem->getVariable(
        0, _displacements[anchor_displacement], Moose::VarKindType::VAR_SOLVER);
    const auto solver_sys_num = anchor_variable.sys().number();
    if (!_problem->isSolverSystemNonlinear(solver_sys_num))
      paramError("displacements", "The in-plane displacements must be nonlinear variables");

    for (unsigned int i = 0; i < _ndisp; ++i)
      if (i != _out_of_plane_direction &&
          _problem->getVariable(0, _displacements[i], Moose::VarKindType::VAR_SOLVER)
                  .sys()
                  .number() != solver_sys_num)
        paramError("displacements",
                   "All in-plane displacements must belong to the same nonlinear system");

    auto & nonlinear_system = _problem->getNonlinearSystemBase(solver_sys_num);
    const auto & scalar_variable = getParam<VariableName>("scalar_out_of_plane_strain");
    if (nonlinear_system.hasScalarVariable(scalar_variable))
      return;
    if (_problem->hasScalarVariable(scalar_variable))
      paramError("scalar_out_of_plane_strain",
                 "Variable '",
                 scalar_variable,
                 "' already exists but is not a nonlinear scalar variable in system '",
                 nonlinear_system.name(),
                 "'");
    if (_problem->hasVariable(scalar_variable))
      paramError("scalar_out_of_plane_strain",
                 "Variable '",
                 scalar_variable,
                 "' already exists as a field variable; a scalar variable is required");

    InputParameters params = _factory.getValidParams("MooseVariableScalar");
    params.set<MooseEnum>("family") = "SCALAR";
    params.set<MooseEnum>("order") = "FIRST";
    params.set<SolverSystemName>("solver_sys") = nonlinear_system.name();
    _problem->addVariable("MooseVariableScalar", scalar_variable, params);
  }

  //
  // Add off diagonal Jacobian kernels
  //
  else if (_current_task == "add_kernel" && _use_ad)
  {
    const std::string k_type = "ADGeneralizedPlaneStrain";
    InputParameters params = _factory.getValidParams(k_type);

    params.applyParameters(parameters(),
                           {"scalar_out_of_plane_strain",
                            "out_of_plane_pressure",
                            "out_of_plane_pressure_function",
                            "factor",
                            "pressure_factor"});
    params.set<std::vector<VariableName>>("scalar_out_of_plane_strain") = {
        getParam<VariableName>("scalar_out_of_plane_strain")};

    if (parameters().isParamSetByUser("out_of_plane_pressure"))
      params.set<FunctionName>("out_of_plane_pressure") =
          getParam<FunctionName>("out_of_plane_pressure");
    if (parameters().isParamSetByUser("out_of_plane_pressure_function"))
      params.set<FunctionName>("out_of_plane_pressure_function") =
          getParam<FunctionName>("out_of_plane_pressure_function");
    if (parameters().isParamSetByUser("factor"))
      params.set<Real>("factor") = getParam<Real>("factor");
    if (parameters().isParamSetByUser("pressure_factor"))
      params.set<Real>("pressure_factor") = getParam<Real>("pressure_factor");

    const auto anchor_displacement = inPlaneDisplacementIndex();
    params.set<NonlinearVariableName>("variable") = _displacements[anchor_displacement];
    _problem->addKernel(k_type, _name + "_ADGeneralizedPlaneStrain", params);
  }
  else if (_use_ad && (_current_task == "add_user_object" || _current_task == "add_scalar_kernel"))
  {
    // ADKernelScalarBase assembles both the elemental resultant and the scalar equation, so the
    // legacy UserObject and ScalarKernel are intentionally unnecessary in AD mode.
  }
  else if (_current_task == "add_kernel")
  {
    std::string k_type = "GeneralizedPlaneStrainOffDiag";
    InputParameters params = _factory.getValidParams(k_type);

registerMooseObject("SolidMechanicsApp", ADGeneralizedPlaneStrain);

InputParameters
ADGeneralizedPlaneStrain::validParams()
{
  InputParameters params = ADKernelScalarBase::validParams();
  params.addClassDescription(
      "Assembles the generalized plane strain scalar equation using automatic differentiation.");
  params.renameCoupledVar("scalar_variable",
                          "scalar_out_of_plane_strain",
                          "Scalar variable for generalized plane strain");
  params.makeParamRequired<std::vector<VariableName>>("scalar_out_of_plane_strain");
  params.addParam<FunctionName>("out_of_plane_pressure_function",
                                "Function used to prescribe pressure (applied toward the body) in "
                                "the out-of-plane direction");
  params.addDeprecatedParam<FunctionName>(
      "out_of_plane_pressure",
      "Function used to prescribe pressure (applied toward the body) in the out-of-plane direction",
      "This has been replaced by 'out_of_plane_pressure_function'");
  params.addParam<MaterialPropertyName>("out_of_plane_pressure_material",
                                        "0",
                                        "Material used to prescribe pressure (applied toward the "
                                        "body) in the out-of-plane direction");
  MooseEnum out_of_plane_direction("x y z", "z");
  params.addParam<MooseEnum>(
      "out_of_plane_direction", out_of_plane_direction, "The direction of the out-of-plane strain");
  params.addDeprecatedParam<Real>(
      "factor",
      "Scale factor applied to prescribed out-of-plane pressure (both material and function)",
      "This has been replaced by 'pressure_factor'");
  params.addParam<Real>(
      "pressure_factor",
      "Scale factor applied to prescribed out-of-plane pressure (both material and function)");
  params.addParam<std::string>("base_name", "Material property base name");
  params.set<bool>("compute_field_residuals") = false;
  params.suppressParameter<bool>("compute_field_residuals");

  return params;
}

ADGeneralizedPlaneStrain::ADGeneralizedPlaneStrain(const InputParameters & parameters)
  : ADKernelScalarBase(parameters),
    _base_name(isParamValid("base_name") ? getParam<std::string>("base_name") + "_" : ""),
    _stress(getADMaterialProperty<RankTwoTensor>(_base_name + "stress")),
    _out_of_plane_pressure_function(parameters.isParamSetByUser("out_of_plane_pressure_function")
                                        ? &getFunction("out_of_plane_pressure_function")
                                    : parameters.isParamSetByUser("out_of_plane_pressure")
                                        ? &getFunction("out_of_plane_pressure")
                                        : nullptr),
    _out_of_plane_pressure_material(getMaterialProperty<Real>("out_of_plane_pressure_material")),
    _pressure_factor(parameters.isParamSetByUser("pressure_factor")
                         ? getParam<Real>("pressure_factor")
                     : parameters.isParamSetByUser("factor") ? getParam<Real>("factor")
                                                             : 1.0),
    _out_of_plane_direction(getParam<MooseEnum>("out_of_plane_direction"))
{
  if (parameters.isParamSetByUser("out_of_plane_pressure_function") &&
      parameters.isParamSetByUser("out_of_plane_pressure"))
    paramError("out_of_plane_pressure_function",
               "Cannot specify both 'out_of_plane_pressure_function' and "
               "'out_of_plane_pressure'");
  if (parameters.isParamSetByUser("pressure_factor") && parameters.isParamSetByUser("factor"))
    paramError("pressure_factor", "Cannot specify both 'pressure_factor' and 'factor'");
}

void
ADGeneralizedPlaneStrain::initialSetup()
{
  if (getBlockCoordSystem() == Moose::COORD_RZ)
    _out_of_plane_direction = 1;
  else if (getBlockCoordSystem() != Moose::COORD_XYZ)
    paramError("out_of_plane_direction",
               "Generalized plane strain supports only Cartesian and axisymmetric coordinate "
               "systems");
}

ADReal
ADGeneralizedPlaneStrain::computeQpResidual()
{
  return 0;
}

ADReal
ADGeneralizedPlaneStrain::computeScalarQpResidual()
{
  const Real out_of_plane_pressure =
      ((_out_of_plane_pressure_function ? _out_of_plane_pressure_function->value(_t, _q_point[_qp])
                                        : 0.0) +
       _out_of_plane_pressure_material[_qp]) *
      _pressure_factor;

  return _stress[_qp](_out_of_plane_direction, _out_of_plane_direction) + out_of_plane_pressure;
}

#include "libmesh/quadrature.h"

InputParameters
ADCompute1DFiniteStrain::validParams()
{
  InputParameters params = ADComputeFiniteStrain::validParams();
  params.addClassDescription("Compute strain increment for finite strain in 1D problem");
  return params;
}

ADCompute1DFiniteStrain::ADCompute1DFiniteStrain(const InputParameters & parameters)
  : ADComputeFiniteStrain(parameters)
{
}

void
ADCompute1DFiniteStrain::computeProperties()
{
  for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
  {
    // Deformation gradient
    auto A = ADRankTwoTensor::initializeFromRows(
        (*_grad_disp[0])[_qp], (*_grad_disp[1])[_qp], (*_grad_disp[2])[_qp]);

    // Old Deformation gradient
    auto Fbar = RankTwoTensor ::initializeFromRows(
        (*_grad_disp_old[0])[_qp], (*_grad_disp_old[1])[_qp], (*_grad_disp_old[2])[_qp]);

    // Compute the displacement gradient dUy/dy and dUz/dz value for 1D problems
    A(1, 1) = computeGradDispYY();
    A(2, 2) = computeGradDispZZ();

    Fbar(1, 1) = computeGradDispYYOld();
    Fbar(2, 2) = computeGradDispZZOld();

    A -= Fbar; // very nearly A = gradU - gradUold, adapted to cylindrical coords

    Fbar.addIa(1.0); // Fbar = ( I + gradUold)

    // Incremental deformation gradient _Fhat = I + A Fbar^-1
    _Fhat[_qp] = A * Fbar.inverse();
    _Fhat[_qp].addIa(1.0);
  }

  for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
    computeQpStrain();
}

#include "libmesh/quadrature.h"

InputParameters
ADCompute1DIncrementalStrain::validParams()
{
  InputParameters params = ADComputeIncrementalStrain::validParams();
  params.addClassDescription("Compute strain increment for small strains in 1D problems.");
  return params;
}

ADCompute1DIncrementalStrain::ADCompute1DIncrementalStrain(const InputParameters & parameters)
  : ADComputeIncrementalStrain(parameters)
{
}

void
ADCompute1DIncrementalStrain::computeTotalStrainIncrement(ADRankTwoTensor & total_strain_increment)
{
  // Deformation gradient calculation for 1D problems
  auto A = ADRankTwoTensor::initializeFromRows(
      (*_grad_disp[0])[_qp], (*_grad_disp[1])[_qp], (*_grad_disp[2])[_qp]);

  // Old Deformation gradient
  auto Fbar = RankTwoTensor ::initializeFromRows(
      (*_grad_disp_old[0])[_qp], (*_grad_disp_old[1])[_qp], (*_grad_disp_old[2])[_qp]);

  // Compute the displacement gradient dUy/dy and dUz/dz value for 1D problems
  A(1, 1) = computeGradDispYY();
  A(2, 2) = computeGradDispZZ();

  Fbar(1, 1) = computeGradDispYYOld();
  Fbar(2, 2) = computeGradDispZZOld();

  A -= Fbar; // very nearly A = gradU - gradUold, adapted to cylindrical coords

  total_strain_increment = 0.5 * (A + A.transpose());
}

#include "libmesh/quadrature.h"

InputParameters
ADCompute1DSmallStrain::validParams()
{
  InputParameters params = ADComputeSmallStrain::validParams();
  params.addClassDescription("Compute a small strain in 1D problem");
  return params;
}

ADCompute1DSmallStrain::ADCompute1DSmallStrain(const InputParameters & parameters)
  : ADComputeSmallStrain(parameters)
{
}

void
ADCompute1DSmallStrain::computeProperties()
{
  for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
  {
    _total_strain[_qp](0, 0) = (*_grad_disp[0])[_qp](0);
    _total_strain[_qp](1, 1) = computeStrainYY();
    _total_strain[_qp](2, 2) = computeStrainZZ();

    _mechanical_strain[_qp] = _total_strain[_qp];

    // Remove the eigenstrain
    for (const auto es : _eigenstrains)
      _mechanical_strain[_qp] -= (*es)[_qp];
  }
}

registerMooseObject("SolidMechanicsApp", ADComputeAxisymmetric1DFiniteStrain);

InputParameters
ADComputeAxisymmetric1DFiniteStrain::validParams()
{
  InputParameters params = ADCompute1DFiniteStrain::validParams();
  params.addClassDescription("Compute a strain increment and rotation increment for finite strains "
                             "in an axisymmetric 1D problem");
  params.addParam<UserObjectName>("subblock_index_provider",
                                  "SubblockIndexProvider user object name");
  params.addCoupledVar("scalar_out_of_plane_strain", "Scalar variable for axisymmetric 1D problem");
  params.addCoupledVar("out_of_plane_strain", "Nonlinear variable for axisymmetric 1D problem");

  return params;
}

ADComputeAxisymmetric1DFiniteStrain::ADComputeAxisymmetric1DFiniteStrain(
    const InputParameters & parameters)
  : ADCompute1DFiniteStrain(parameters),
    _disp_old_0(coupledValueOld("displacements", 0)),
    _subblock_id_provider(isParamValid("subblock_index_provider")
                              ? &getUserObject<SubblockIndexProvider>("subblock_index_provider")
                              : nullptr),
    _has_out_of_plane_strain(isCoupled("out_of_plane_strain")),
    _out_of_plane_strain(_has_out_of_plane_strain ? adCoupledValue("out_of_plane_strain")
                                                  : _ad_zero),
    _out_of_plane_strain_old(_has_out_of_plane_strain ? coupledValueOld("out_of_plane_strain")
                                                      : _zero),
    _has_scalar_out_of_plane_strain(isCoupledScalar("scalar_out_of_plane_strain"))
{
  if (_has_out_of_plane_strain && _has_scalar_out_of_plane_strain)
    mooseError("Must define only one of out_of_plane_strain or scalar_out_of_plane_strain");

  if (_has_scalar_out_of_plane_strain)
  {
    const auto nscalar_strains = coupledScalarComponents("scalar_out_of_plane_strain");
    _scalar_out_of_plane_strain.resize(nscalar_strains);
    _scalar_out_of_plane_strain_old.resize(nscalar_strains);
    for (unsigned int i = 0; i < nscalar_strains; ++i)
    {
      _scalar_out_of_plane_strain[i] = &adCoupledScalarValue("scalar_out_of_plane_strain", i);
      _scalar_out_of_plane_strain_old[i] = &coupledScalarValueOld("scalar_out_of_plane_strain", i);
    }
  }
}

void
ADComputeAxisymmetric1DFiniteStrain::initialSetup()
{
  ADComputeIncrementalStrainBase::initialSetup();

  if (getBlockCoordSystem() != Moose::COORD_RZ)
    mooseError("The coordinate system must be set to RZ for Axisymmetric geometries.");
}

unsigned int
ADComputeAxisymmetric1DFiniteStrain::getCurrentSubblockIndex() const
{
  return _subblock_id_provider ? _subblock_id_provider->getSubblockIndex(*_current_elem) : 0;
}

ADReal
ADComputeAxisymmetric1DFiniteStrain::computeGradDispYY()
{
  using std::exp;
  if (_has_scalar_out_of_plane_strain)
    return exp((*_scalar_out_of_plane_strain[getCurrentSubblockIndex()])[0]) - 1.0;
  else
    return exp(_out_of_plane_strain[_qp]) - 1.0;
}

Real
ADComputeAxisymmetric1DFiniteStrain::computeGradDispYYOld()
{
  using std::exp;
  if (_has_scalar_out_of_plane_strain)
    return exp((*_scalar_out_of_plane_strain_old[getCurrentSubblockIndex()])[0]) - 1.0;
  else
    return exp(_out_of_plane_strain_old[_qp]) - 1.0;
}

ADReal
ADComputeAxisymmetric1DFiniteStrain::computeGradDispZZ()
{
  if (!MooseUtils::absoluteFuzzyEqual(_q_point[_qp](0), 0.0))
    return (*_disp[0])[_qp] / _q_point[_qp](0);
  else
    return 0.0;
}

Real
ADComputeAxisymmetric1DFiniteStrain::computeGradDispZZOld()
{
  if (!MooseUtils::absoluteFuzzyEqual(_q_point[_qp](0), 0.0))
    return _disp_old_0[_qp] / _q_point[_qp](0);
  else
    return 0.0;
}

registerMooseObject("SolidMechanicsApp", ADComputeAxisymmetric1DIncrementalStrain);

InputParameters
ADComputeAxisymmetric1DIncrementalStrain::validParams()
{
  InputParameters params = ADCompute1DIncrementalStrain::validParams();
  params.addClassDescription(
      "Compute strain increment for small strains in an axisymmetric 1D problem");
  params.addParam<UserObjectName>("subblock_index_provider",
                                  "SubblockIndexProvider user object name");
  params.addCoupledVar("scalar_out_of_plane_strain", "Scalar variable for axisymmetric 1D problem");
  params.addCoupledVar("out_of_plane_strain", "Nonlinear variable for axisymmetric 1D problem");

  return params;
}

ADComputeAxisymmetric1DIncrementalStrain::ADComputeAxisymmetric1DIncrementalStrain(
    const InputParameters & parameters)
  : ADCompute1DIncrementalStrain(parameters),
    _disp_old_0(coupledValueOld("displacements", 0)),
    _subblock_id_provider(isParamValid("subblock_index_provider")
                              ? &getUserObject<SubblockIndexProvider>("subblock_index_provider")
                              : nullptr),
    _has_out_of_plane_strain(isCoupled("out_of_plane_strain")),
    _out_of_plane_strain(_has_out_of_plane_strain ? adCoupledValue("out_of_plane_strain")
                                                  : _ad_zero),
    _out_of_plane_strain_old(_has_out_of_plane_strain ? coupledValueOld("out_of_plane_strain")
                                                      : _zero),
    _has_scalar_out_of_plane_strain(isCoupledScalar("scalar_out_of_plane_strain"))
{
  if (_has_out_of_plane_strain && _has_scalar_out_of_plane_strain)
    mooseError("Must define only one of out_of_plane_strain or scalar_out_of_plane_strain");

  if (_has_scalar_out_of_plane_strain)
  {
    const auto nscalar_strains = coupledScalarComponents("scalar_out_of_plane_strain");
    _scalar_out_of_plane_strain.resize(nscalar_strains);
    _scalar_out_of_plane_strain_old.resize(nscalar_strains);
    for (unsigned int i = 0; i < nscalar_strains; ++i)
    {
      _scalar_out_of_plane_strain[i] = &adCoupledScalarValue("scalar_out_of_plane_strain", i);
      _scalar_out_of_plane_strain_old[i] = &coupledScalarValueOld("scalar_out_of_plane_strain", i);
    }
  }
}

void
ADComputeAxisymmetric1DIncrementalStrain::initialSetup()
{
  ADComputeIncrementalStrainBase::initialSetup();

  if (getBlockCoordSystem() != Moose::COORD_RZ)
    mooseError("The coordinate system must be set to RZ for Axisymmetric 1D simulations");
}

unsigned int
ADComputeAxisymmetric1DIncrementalStrain::getCurrentSubblockIndex() const
{
  return _subblock_id_provider ? _subblock_id_provider->getSubblockIndex(*_current_elem) : 0;
}

ADReal
ADComputeAxisymmetric1DIncrementalStrain::computeGradDispYY()
{
  if (_has_scalar_out_of_plane_strain)
    return (*_scalar_out_of_plane_strain[getCurrentSubblockIndex()])[0];
  else
    return _out_of_plane_strain[_qp];
}

Real
ADComputeAxisymmetric1DIncrementalStrain::computeGradDispYYOld()
{
  if (_has_scalar_out_of_plane_strain)
    return (*_scalar_out_of_plane_strain_old[getCurrentSubblockIndex()])[0];
  else
    return _out_of_plane_strain_old[_qp];
}

ADReal
ADComputeAxisymmetric1DIncrementalStrain::computeGradDispZZ()
{
  if (!MooseUtils::absoluteFuzzyEqual(_q_point[_qp](0), 0.0))
    return (*_disp[0])[_qp] / _q_point[_qp](0);
  else
    return 0.0;
}

Real
ADComputeAxisymmetric1DIncrementalStrain::computeGradDispZZOld()
{
  if (!MooseUtils::absoluteFuzzyEqual(_q_point[_qp](0), 0.0))
    return _disp_old_0[_qp] / _q_point[_qp](0);
  else
    return 0.0;
}

registerMooseObject("SolidMechanicsApp", ADComputeAxisymmetric1DSmallStrain);

InputParameters
ADComputeAxisymmetric1DSmallStrain::validParams()
{
  InputParameters params = ADCompute1DSmallStrain::validParams();
  params.addClassDescription("Compute a small strain in an Axisymmetric 1D problem");
  params.addParam<UserObjectName>("subblock_index_provider",
                                  "SubblockIndexProvider user object name");
  params.addCoupledVar("scalar_out_of_plane_strain", "Scalar variable for axisymmetric 1D problem");
  params.addCoupledVar("out_of_plane_strain", "Nonlinear variable for axisymmetric 1D problem");

  return params;
}

ADComputeAxisymmetric1DSmallStrain::ADComputeAxisymmetric1DSmallStrain(
    const InputParameters & parameters)
  : ADCompute1DSmallStrain(parameters),
    _subblock_id_provider(isParamValid("subblock_index_provider")
                              ? &getUserObject<SubblockIndexProvider>("subblock_index_provider")
                              : nullptr),
    _has_out_of_plane_strain(isCoupled("out_of_plane_strain")),
    _out_of_plane_strain(_has_out_of_plane_strain ? adCoupledValue("out_of_plane_strain")
                                                  : _ad_zero),
    _has_scalar_out_of_plane_strain(isCoupledScalar("scalar_out_of_plane_strain"))
{
  if (_has_out_of_plane_strain && _has_scalar_out_of_plane_strain)
    mooseError("Must define only one of out_of_plane_strain or scalar_out_of_plane_strain");

  if (_has_scalar_out_of_plane_strain)
  {
    const auto nscalar_strains = coupledScalarComponents("scalar_out_of_plane_strain");
    _scalar_out_of_plane_strain.resize(nscalar_strains);
    for (unsigned int i = 0; i < nscalar_strains; ++i)
      _scalar_out_of_plane_strain[i] = &adCoupledScalarValue("scalar_out_of_plane_strain", i);
  }
}

void
ADComputeAxisymmetric1DSmallStrain::initialSetup()
{
  ADComputeStrainBase::initialSetup();

  if (getBlockCoordSystem() != Moose::COORD_RZ)
    mooseError("The coordinate system must be set to RZ for Axisymmetric geometries.");
}

unsigned int
ADComputeAxisymmetric1DSmallStrain::getCurrentSubblockIndex() const
{
  return _subblock_id_provider ? _subblock_id_provider->getSubblockIndex(*_current_elem) : 0;
}

ADReal
ADComputeAxisymmetric1DSmallStrain::computeStrainYY()
{
  if (_has_scalar_out_of_plane_strain)
    return (*_scalar_out_of_plane_strain[getCurrentSubblockIndex()])[0];
  else
    return _out_of_plane_strain[_qp];
}

ADReal
ADComputeAxisymmetric1DSmallStrain::computeStrainZZ()
{
  if (!MooseUtils::absoluteFuzzyEqual(_q_point[_qp](0), 0.0))
    return (*_disp[0])[_qp] / _q_point[_qp](0);
  else
    return 0.0;
}

  if (_scalar_out_of_plane_strain_coupled)
  {
    const auto nscalar_strains = coupledScalarComponents("scalar_out_of_plane_strain");
    _scalar_out_of_plane_strain.resize(nscalar_strains);
    _scalar_out_of_plane_strain_old.resize(nscalar_strains);
    for (unsigned int i = 0; i < nscalar_strains; ++i)
    {
      _scalar_out_of_plane_strain[i] = &adCoupledScalarValue("scalar_out_of_plane_strain", i);
      _scalar_out_of_plane_strain_old[i] = &coupledScalarValueOld("scalar_out_of_plane_strain", i);
    }
  }
}

   * D = log(sqrt(Fhat^T * Fhat)) / dt
   */
  if (_scalar_out_of_plane_strain_coupled)
    return exp((*_scalar_out_of_plane_strain[getCurrentSubblockIndex()])[0]) - 1.0;
  else
    return exp(_out_of_plane_strain[_qp]) - 1.0;
}

{
  using std::exp;
  if (_scalar_out_of_plane_strain_coupled)
    return exp((*_scalar_out_of_plane_strain_old[getCurrentSubblockIndex()])[0]) - 1.0;
  else
    return exp(_out_of_plane_strain_old[_qp]) - 1.0;
}

  if (_scalar_out_of_plane_strain_coupled)
  {
    const auto nscalar_strains = coupledScalarComponents("scalar_out_of_plane_strain");
    _scalar_out_of_plane_strain.resize(nscalar_strains);
    _scalar_out_of_plane_strain_old.resize(nscalar_strains);
    for (unsigned int i = 0; i < nscalar_strains; ++i)
    {
      _scalar_out_of_plane_strain[i] = &adCoupledScalarValue("scalar_out_of_plane_strain", i);
      _scalar_out_of_plane_strain_old[i] = &coupledScalarValueOld("scalar_out_of_plane_strain", i);
    }
  }
}

ADComputePlaneIncrementalStrain::computeOutOfPlaneGradDisp()
{
  if (_scalar_out_of_plane_strain_coupled)
    return (*_scalar_out_of_plane_strain[getCurrentSubblockIndex()])[0];
  else
    return _out_of_plane_strain[_qp];
}

ADComputePlaneIncrementalStrain::computeOutOfPlaneGradDispOld()
{
  if (_scalar_out_of_plane_strain_coupled)
    return (*_scalar_out_of_plane_strain_old[getCurrentSubblockIndex()])[0];
  else
    return _out_of_plane_strain_old[_qp];
}

  if (_scalar_out_of_plane_strain_coupled)
  {
    const auto nscalar_strains = coupledScalarComponents("scalar_out_of_plane_strain");
    _scalar_out_of_plane_strain.resize(nscalar_strains);
    for (unsigned int i = 0; i < nscalar_strains; ++i)
      _scalar_out_of_plane_strain[i] = &adCoupledScalarValue("scalar_out_of_plane_strain", i);
  }
}

ADComputePlaneSmallStrain::computeOutOfPlaneStrain()
{
  if (_scalar_out_of_plane_strain_coupled)
    return (*_scalar_out_of_plane_strain[getCurrentSubblockIndex()])[0];
  else
    return _out_of_plane_strain[_qp];
}

                                     "factor",
                                     "pressure_factor"});
      action_params.set<bool>("use_displaced_mesh") = _use_displaced_mesh;
      action_params.set<bool>("use_automatic_differentiation") = _use_ad;

      if (parameters().isParamSetByUser("out_of_plane_pressure"))
        action_params.set<FunctionName>("out_of_plane_pressure") =

           _planar_formulation == PlanarFormulation::PlaneStrain ||
           _planar_formulation == PlanarFormulation::GeneralizedPlaneStrain)
  {
    if (_use_ad && _planar_formulation == PlanarFormulation::PlaneStrain)
      paramError("use_automatic_differentiation", "AD not setup for use with PlaneStrain");

    std::map<std::pair<Moose::CoordinateSystemType, StrainAndIncrement>, std::string> type_map = {
        {{Moose::COORD_XYZ, StrainAndIncrement::SmallTotal}, "ComputePlaneSmallStrain"},