NodalExtremeValue

Reports the minimum or maximum value of a variable.

You can optionally provide a "proxy_variable", which will change the behavior of this postprocessor to find the node at which the proxy variable reaches the max/min value, and then return the value of the specified variable at that node.

The corresponding postprocessor that find extreme values of variables evaluated inside elements (at quadrature points) is ElementExtremeValue

Example Input File Syntax

Input Parameters

variableThe name of the variable that this postprocessor operates on
C++ Type:std::vector<VariableName>
Controllable:No
Description:The name of the variable that this postprocessor operates on

Required Parameters

blockThe list of blocks (ids or names) that this object will be applied
C++ Type:std::vector<SubdomainName>
Controllable:No
Description:The list of blocks (ids or names) that this object will be applied
boundaryThe list of boundaries (ids or names) from the mesh where this boundary condition applies
C++ Type:std::vector<BoundaryName>
Controllable:No
Description:The list of boundaries (ids or names) from the mesh where this boundary condition applies
execute_onTIMESTEP_ENDThe list of flag(s) indicating when this object should be executed, the available options include NONE, INITIAL, LINEAR, NONLINEAR, TIMESTEP_END, TIMESTEP_BEGIN, FINAL, CUSTOM, ALWAYS.
Default:TIMESTEP_END
C++ Type:ExecFlagEnum
Options:NONE, INITIAL, LINEAR, NONLINEAR, TIMESTEP_END, TIMESTEP_BEGIN, FINAL, CUSTOM, TRANSFER, ALWAYS
Controllable:No
Description:The list of flag(s) indicating when this object should be executed, the available options include NONE, INITIAL, LINEAR, NONLINEAR, TIMESTEP_END, TIMESTEP_BEGIN, FINAL, CUSTOM, ALWAYS.
proxy_variableThe name of the variable to use to identify the location at which the variable value should be taken; if not provided, this defaults to the 'variable'.
C++ Type:std::vector<VariableName>
Controllable:No
Description:The name of the variable to use to identify the location at which the variable value should be taken; if not provided, this defaults to the 'variable'.
unique_node_executeFalseWhen false (default), block restricted objects will have the execute method called multiple times on a single node if the node lies on a interface between two subdomains.
Default:False
C++ Type:bool
Controllable:No
Description:When false (default), block restricted objects will have the execute method called multiple times on a single node if the node lies on a interface between two subdomains.
value_typemaxType of extreme value to return. 'max' returns the maximum value. 'min' returns the minimum value.
Default:max
C++ Type:MooseEnum
Options:max, min
Controllable:No
Description:Type of extreme value to return. 'max' returns the maximum value. 'min' returns the minimum value.

Optional Parameters

allow_duplicate_execution_on_initialFalseIn the case where this UserObject is depended upon by an initial condition, allow it to be executed twice during the initial setup (once before the IC and again after mesh adaptivity (if applicable).
Default:False
C++ Type:bool
Controllable:No
Description:In the case where this UserObject is depended upon by an initial condition, allow it to be executed twice during the initial setup (once before the IC and again after mesh adaptivity (if applicable).
control_tagsAdds user-defined labels for accessing object parameters via control logic.
C++ Type:std::vector<std::string>
Controllable:No
Description:Adds user-defined labels for accessing object parameters via control logic.
enableTrueSet the enabled status of the MooseObject.
Default:True
C++ Type:bool
Controllable:Yes
Description:Set the enabled status of the MooseObject.
force_postauxFalseForces the UserObject to be executed in POSTAUX
Default:False
C++ Type:bool
Controllable:No
Description:Forces the UserObject to be executed in POSTAUX
force_preauxFalseForces the UserObject to be executed in PREAUX
Default:False
C++ Type:bool
Controllable:No
Description:Forces the UserObject to be executed in PREAUX
force_preicFalseForces the UserObject to be executed in PREIC during initial setup
Default:False
C++ Type:bool
Controllable:No
Description:Forces the UserObject to be executed in PREIC during initial setup
outputsVector of output names were you would like to restrict the output of variables(s) associated with this object
C++ Type:std::vector<OutputName>
Controllable:No
Description:Vector of output names were you would like to restrict the output of variables(s) associated with this object
seed0The seed for the master random number generator
Default:0
C++ Type:unsigned int
Controllable:No
Description:The seed for the master random number generator
use_displaced_meshFalseWhether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used.
Default:False
C++ Type:bool
Controllable:No
Description:Whether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used.

Advanced Parameters

Input Files

(modules/tensor_mechanics/test/tests/thermal_expansion_function/small_linear.i)
(modules/stochastic_tools/examples/surrogates/combined/trans_diff_2d/trans_diff_sub.i)
(modules/tensor_mechanics/test/tests/dynamics/time_integration/hht_test.i)
(modules/tensor_mechanics/test/tests/dynamics/time_integration/newmark.i)
(modules/tensor_mechanics/test/tests/plane_stress/weak_plane_stress_incremental.i)
(modules/tensor_mechanics/test/tests/thermal_expansion_function/finite_const.i)
(modules/combined/test/tests/restart-transient-from-ss-with-stateful/sub_tr.i)
(modules/tensor_mechanics/test/tests/ad_thermal_expansion_function/small_const.i)
(modules/tensor_mechanics/test/tests/dynamics/time_integration/hht_test_action.i)
(modules/tensor_mechanics/examples/coal_mining/fine.i)
(test/tests/userobjects/Terminator/terminator.i)
(modules/porous_flow/examples/coal_mining/fine_with_fluid.i)
(modules/porous_flow/examples/coal_mining/coarse_with_fluid.i)
(modules/tensor_mechanics/examples/coal_mining/coarse.i)
(modules/stochastic_tools/test/tests/surrogates/load_store/sub.i)
(modules/tensor_mechanics/test/tests/ad_thermal_expansion_function/finite_linear.i)
(modules/tensor_mechanics/test/tests/dynamics/rayleigh_damping/rayleigh_hht_ti.i)
(modules/tensor_mechanics/test/tests/volumetric_eigenstrain/volumetric_mechanical.i)
(modules/tensor_mechanics/test/tests/plane_stress/weak_plane_stress_finite.i)
(modules/tensor_mechanics/test/tests/dynamics/rayleigh_damping/rayleigh_newmark.i)
(modules/tensor_mechanics/test/tests/ad_thermal_expansion_function/finite_const.i)
(modules/tensor_mechanics/test/tests/plane_stress/ad_weak_plane_stress_finite.i)
(modules/tensor_mechanics/test/tests/dynamics/rayleigh_damping/rayleigh_newmark_material_dependent.i)
(modules/tensor_mechanics/test/tests/dynamics/time_integration/newmark_action.i)
(modules/stochastic_tools/examples/surrogates/sub.i)
(modules/tensor_mechanics/test/tests/volumetric_eigenstrain/volumetric_eigenstrain.i)
(modules/tensor_mechanics/test/tests/dynamics/rayleigh_damping/rayleigh_hht.i)
(test/tests/misc/block_user_object_check/block_check.i)
(modules/tensor_mechanics/test/tests/thermal_expansion_function/finite_linear.i)
(modules/tensor_mechanics/test/tests/ad_thermal_expansion_function/small_linear.i)
(test/tests/postprocessors/nodal_extreme_value/nodal_extreme_pps_test.i)
(modules/combined/test/tests/gap_heat_transfer_mortar/finite-2d/varied_pressure_thermomechanical_mortar.i)
(test/tests/misc/check_error/double_restrict_uo.i)
(modules/stochastic_tools/examples/surrogates/polynomial_regression/sub.i)
(modules/stochastic_tools/test/tests/vectorpostprocessors/multiple_stochastic_results/sub.i)
(modules/tensor_mechanics/test/tests/plane_stress/ad_weak_plane_stress_incremental.i)
(modules/tensor_mechanics/test/tests/plane_stress/weak_plane_stress_small.i)
(modules/tensor_mechanics/test/tests/thermal_expansion_function/small_const.i)
(modules/tensor_mechanics/test/tests/action/action_multi_eigenstrain_same_conditions.i)
(modules/tensor_mechanics/test/tests/action/ad_converter_action_multi_eigenstrain.i)
(modules/tensor_mechanics/test/tests/action/action_multi_eigenstrain.i)
(modules/contact/examples/3d_berkovich/indenter_berkovich_friction.i)
(modules/stochastic_tools/test/tests/surrogates/poly_chaos/sub.i)
(modules/tensor_mechanics/test/tests/plane_stress/ad_weak_plane_stress_small.i)
(modules/stochastic_tools/examples/surrogates/gaussian_process/sub.i)
(modules/combined/test/tests/restart-transient-from-ss-with-stateful/sub_ss.i)
(modules/tensor_mechanics/tutorials/introduction/mech_step04a.i)
(modules/tensor_mechanics/test/tests/dynamics/time_integration/hht_test_ti.i)

(modules/tensor_mechanics/test/tests/thermal_expansion_function/small_linear.i)

# This tests the thermal expansion coefficient function using both
# options to specify that function: mean and instantaneous.  There
# two blocks, each containing a single element, and these use the
# two variants of the function.

# In this test, the instantaneous CTE function is a linear function
# while the mean CTE function is an analytic function designed to
# give the same response.  If \bar{alpha}(T) is the mean CTE function,
# and \alpha(T) is the instantaneous CTE function,

# \bar{\alpha}(T) = 1/(T-Tref) \intA^{T}_{Tsf} \alpha(T) dT

# where Tref is the reference temperature used to define the mean CTE
# function, and Tsf is the stress-free temperature.

# This version of the test uses small deformation theory.  The results
# from the two models are identical.

[Mesh]
  file = 'blocks.e'
[]

[GlobalParams]
  displacements = 'disp_x disp_y disp_z'
[]

[AuxVariables]
  [./temp]
    order = FIRST
    family = LAGRANGE
  [../]
[]

[Modules/TensorMechanics/Master]
  [./all]
    strain = SMALL
    add_variables = true
    eigenstrain_names = eigenstrain
    generate_output = 'strain_xx strain_yy strain_zz'
  [../]
[]

[BCs]
  [./left]
    type = DirichletBC
    variable = disp_x
    boundary = 3
    value = 0.0
  [../]

  [./bottom]
    type = DirichletBC
    variable = disp_y
    boundary = 2
    value = 0.0
  [../]

  [./back]
    type = DirichletBC
    variable = disp_z
    boundary = 1
    value = 0.0
  [../]
[]

[AuxKernels]
  [./temp]
    type = FunctionAux
    variable = temp
    block = '1 2'
    function = temp_func
  [../]
[]

[Materials]
  [./elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 1e6
    poissons_ratio = 0.3
  [../]
  [./small_stress]
    type = ComputeLinearElasticStress
  [../]
  [./thermal_expansion_strain1]
    type = ComputeMeanThermalExpansionFunctionEigenstrain
    block = 1
    thermal_expansion_function = cte_func_mean
    thermal_expansion_function_reference_temperature = 0.5
    stress_free_temperature = 0.0
    temperature = temp
    eigenstrain_name = eigenstrain
  [../]
  [./thermal_expansion_strain2]
    type = ComputeInstantaneousThermalExpansionFunctionEigenstrain
    block = 2
    thermal_expansion_function = cte_func_inst
    stress_free_temperature = 0.0
    temperature = temp
    eigenstrain_name = eigenstrain
  [../]
[]

[Functions]
  [./cte_func_mean]
    type = ParsedFunction
    vars = 'tsf tref scale' #stress free temp, reference temp, scale factor
    vals = '0.0 0.5  1e-4'
    value = 'scale * (0.5 * t^2 - 0.5 * tsf^2) / (t - tref)'
  [../]
  [./cte_func_inst]
    type = PiecewiseLinear
    xy_data = '0 0.0
               2 2.0'
    scale_factor = 1e-4
  [../]

  [./temp_func]
    type = PiecewiseLinear
    xy_data = '0 1
               1 2'
  [../]
[]

[Postprocessors]
  [./disp_1]
    type = NodalExtremeValue
    variable = disp_x
    boundary = 101
  [../]

  [./disp_2]
    type = NodalExtremeValue
    variable = disp_x
    boundary = 102
  [../]
[]

[Executioner]
  type = Transient

  solve_type = PJFNK
  l_max_its = 100
  l_tol = 1e-4
  nl_abs_tol = 1e-8
  nl_rel_tol = 1e-12

  start_time = 0.0
  end_time = 1.0
  dt = 0.1
[]

[Outputs]
  csv = true
[]

(modules/stochastic_tools/examples/surrogates/combined/trans_diff_2d/trans_diff_sub.i)

[Functions]
  [src_func]
    type = ParsedFunction
    value = "1000*sin(f*t)"
    vars = 'f'
    vals = '20'
  []
[]

[Mesh]
  [msh]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 100
    xmin = -0.5
    xmax = 0.5
    ny = 100
    ymin = -0.5
    ymax = 0.5
  []
  [source_domain]
    type = ParsedSubdomainMeshGenerator
    input = msh
    combinatorial_geometry = '(x<0.1 & x>-0.1) & (y<0.1 & y>-0.1)'
    block_id=1
  []
[]

[Variables]
  [T]
    initial_condition = 300
  []
[]

[Kernels]
  [diffusion]
    type = MatDiffusion
    variable = T
    diffusivity = diff_coeff
  []
  [source]
    type = BodyForce
    variable = T
    function = src_func
    block = 1
  []
  [time_deriv]
    type = TimeDerivative
    variable = T
  []
[]

[Materials]
  [diff_coeff]
    type = ParsedMaterial
    f_name = diff_coeff
    args = 'T'
    constant_names = 'C'
    constant_expressions = 0.02
    function = 'C * pow(300/T, 2)'
  []
[]

[BCs]
  [neumann_all]
    type = NeumannBC
    variable = T
    boundary = 'left right top bottom'
    value = 0
  []
[]

[Executioner]
  type = Transient
  num_steps = 100
  dt = 0.01
  solve_type = PJFNK
  petsc_options_iname = '-pc_type -pc_hypre_type'
  petsc_options_value = 'hypre boomeramg'
  nl_rel_tol = 1e-6
  l_abs_tol = 1e-6
  timestep_tolerance = 1e-6
[]

[Postprocessors]
  [max]
    type = NodalExtremeValue
    variable = T
  []
  [min]
    type = NodalExtremeValue
    variable = T
    value_type = min
  []
  [time_max]
    type = TimeExtremeValue
    postprocessor = max
  []
  [time_min]
    type = TimeExtremeValue
    postprocessor = min
    value_type = min
  []
[]

(modules/tensor_mechanics/test/tests/dynamics/time_integration/hht_test.i)

# Test for  HHT time integration

# The test is for an 1D bar element of  unit length fixed on one end
# with a ramped pressure boundary condition applied to the other end.
# alpha, beta and gamma are HHT time integration parameters
# The equation of motion in terms of matrices is:
#
# M*accel + alpha*(K*disp - K*disp_old) + K*disp = P(t+alpha dt)*Area
#
# Here M is the mass matrix, K is the stiffness matrix, P is the applied pressure
#
# This equation is equivalent to:
#
# density*accel + alpha*(Div stress - Div stress_old) +Div Stress= P(t+alpha dt)
#
# The first term on the left is evaluated using the Inertial force kernel
# The next two terms on the left involving alpha are evaluated using the
# DynamicStressDivergenceTensors Kernel
# The residual due to Pressure is evaluated using Pressure boundary condition
#
# The system will come to steady state slowly after the pressure becomes constant.
# Alpha equal to zero will result in Newmark integration.

[Mesh]
  type = GeneratedMesh
  dim = 3
  nx = 1
  ny = 1
  nz = 1
  xmin = 0.0
  xmax = 0.1
  ymin = 0.0
  ymax = 1.0
  zmin = 0.0
  zmax = 0.1
[]


[Variables]
  [./disp_x]
  [../]
  [./disp_y]
  [../]
  [./disp_z]
  [../]
[]

[AuxVariables]
  [./vel_x]
  [../]
  [./accel_x]
  [../]
  [./vel_y]
  [../]
  [./accel_y]
  [../]
  [./vel_z]
  [../]
  [./accel_z]
  [../]
  [./stress_yy]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./strain_yy]
    order = CONSTANT
    family = MONOMIAL
  [../]

[]

[Kernels]
  [./DynamicTensorMechanics]
    displacements = 'disp_x disp_y disp_z'
    hht_alpha = 0.11
  [../]
  [./inertia_x]
    type = InertialForce
    variable = disp_x
    velocity = vel_x
    acceleration = accel_x
    beta = 0.25
    gamma = 0.5
  [../]
  [./inertia_y]
    type = InertialForce
    variable = disp_y
    velocity = vel_y
    acceleration = accel_y
    beta = 0.25
    gamma = 0.5
  [../]
  [./inertia_z]
    type = InertialForce
    variable = disp_z
    velocity = vel_z
    acceleration = accel_z
    beta = 0.25
    gamma = 0.5
  [../]

[]

[AuxKernels]
  [./accel_x]
    type = NewmarkAccelAux
    variable = accel_x
    displacement = disp_x
    velocity = vel_x
    beta = 0.25
    execute_on = timestep_end
  [../]
  [./vel_x]
    type = NewmarkVelAux
    variable = vel_x
    acceleration = accel_x
    gamma = 0.5
    execute_on = timestep_end
  [../]
  [./accel_y]
    type = NewmarkAccelAux
    variable = accel_y
    displacement = disp_y
    velocity = vel_y
    beta = 0.25
    execute_on = timestep_end
  [../]
  [./vel_y]
    type = NewmarkVelAux
    variable = vel_y
    acceleration = accel_y
    gamma = 0.5
    execute_on = timestep_end
  [../]
  [./accel_z]
    type = NewmarkAccelAux
    variable = accel_z
    displacement = disp_z
    velocity = vel_z
    beta = 0.25
    execute_on = timestep_end
  [../]
  [./vel_z]
    type = NewmarkVelAux
    variable = vel_z
    acceleration = accel_z
    gamma = 0.5
    execute_on = timestep_end
  [../]
  [./stress_yy]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_yy
    index_i = 0
    index_j = 1
  [../]
  [./strain_yy]
    type = RankTwoAux
    rank_two_tensor = total_strain
    variable = strain_yy
    index_i = 0
    index_j = 1
  [../]

[]


[BCs]
  [./top_y]
    type = DirichletBC
    variable = disp_y
    boundary = top
    value=0.0
  [../]
  [./top_x]
    type = DirichletBC
    variable = disp_x
    boundary = top
    value=0.0
  [../]
  [./top_z]
    type = DirichletBC
    variable = disp_z
    boundary = top
    value=0.0
  [../]
  [./bottom_x]
    type = DirichletBC
    variable = disp_x
    boundary = bottom
    value=0.0
  [../]
  [./bottom_z]
    type = DirichletBC
    variable = disp_z
    boundary = bottom
    value=0.0
  [../]
  [./Pressure]
    [./Side1]
      boundary = bottom
      function = pressure
      factor = 1
      hht_alpha = 0.11
      displacements = 'disp_x disp_y disp_z'
    [../]
  [../]
[]

[Materials]
  [./Elasticity_tensor]
    type = ComputeElasticityTensor
    block = 0
    fill_method = symmetric_isotropic
    C_ijkl = '210e9 0'
  [../]

  [./strain]
    type = ComputeSmallStrain
    block = 0
    displacements = 'disp_x disp_y disp_z'
  [../]

  [./stress]
    type = ComputeLinearElasticStress
    block = 0
  [../]
  [./density]
    type = GenericConstantMaterial
    block = 0
    prop_names = 'density'
    prop_values = '7750'
  [../]

[]

[Executioner]
  type = Transient
  start_time = 0
  end_time = 2
  dt = 0.1
[]


[Functions]
  [./pressure]
    type = PiecewiseLinear
    x = '0.0 0.1 0.2 1.0 2.0 5.0'
    y = '0.0 0.1 0.2 1.0 1.0 1.0'
    scale_factor = 1e9
  [../]
[]

[Postprocessors]
  [./_dt]
    type = TimestepSize
  [../]
  [./disp]
    type = NodalExtremeValue
    variable = disp_y
    boundary = bottom
  [../]
  [./vel]
    type = NodalExtremeValue
    variable = vel_y
    boundary = bottom
  [../]
  [./accel]
    type = NodalExtremeValue
    variable = accel_y
    boundary = bottom
  [../]
  [./stress_yy]
    type = ElementAverageValue
    variable = stress_yy
  [../]
  [./strain_yy]
    type = ElementAverageValue
    variable = strain_yy
  [../]
[]

[Outputs]
  exodus = true
  perf_graph = true
[]

(modules/tensor_mechanics/test/tests/dynamics/time_integration/newmark.i)

# Test for  Newmark time integration

# The test is for an 1D bar element of  unit length fixed on one end
# with a ramped pressure boundary condition applied to the other end.
# beta and gamma are Newmark time integration parameters
# The equation of motion in terms of matrices is:
#
# M*accel + K*disp = P*Area
#
# Here M is the mass matrix, K is the stiffness matrix, P is the applied pressure
#
# This equation is equivalent to:
#
# density*accel + Div Stress = P
#
# The first term on the left is evaluated using the Inertial force kernel
# The last term on the left is evaluated using StressDivergenceTensors
# The residual due to Pressure is evaluated using Pressure boundary condition

[Mesh]
  type = GeneratedMesh
  dim = 3
  xmax = 0.1
  ymax = 1.0
  zmax = 0.1
[]

[Variables]
  [disp_x]
  []
  [disp_y]
  []
  [disp_z]
  []
[]

[AuxVariables]
  [vel_x]
  []
  [accel_x]
  []
  [vel_y]
  []
  [accel_y]
  []
  [vel_z]
  []
  [accel_z]
  []
  [stress_yy]
    order = CONSTANT
    family = MONOMIAL
  []
  [strain_yy]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[Kernels]
  [TensorMechanics]
    displacements = 'disp_x disp_y disp_z'
  []
  [inertia_x]
    type = InertialForce
    variable = disp_x
    velocity = vel_x
    acceleration = accel_x
    beta = 0.25
    gamma = 0.5
  []
  [inertia_y]
    type = InertialForce
    variable = disp_y
    velocity = vel_y
    acceleration = accel_y
    beta = 0.25
    gamma = 0.5
  []
  [inertia_z]
    type = InertialForce
    variable = disp_z
    velocity = vel_z
    acceleration = accel_z
    beta = 0.25
    gamma = 0.5
  []

[]

[AuxKernels]
  [accel_x]
    type = NewmarkAccelAux
    variable = accel_x
    displacement = disp_x
    velocity = vel_x
    beta = 0.25
    execute_on = timestep_end
  []
  [vel_x]
    type = NewmarkVelAux
    variable = vel_x
    acceleration = accel_x
    gamma = 0.5
    execute_on = timestep_end
  []
  [accel_y]
    type = NewmarkAccelAux
    variable = accel_y
    displacement = disp_y
    velocity = vel_y
    beta = 0.25
    execute_on = timestep_end
  []
  [vel_y]
    type = NewmarkVelAux
    variable = vel_y
    acceleration = accel_y
    gamma = 0.5
    execute_on = timestep_end
  []
  [accel_z]
    type = NewmarkAccelAux
    variable = accel_z
    displacement = disp_z
    velocity = vel_z
    beta = 0.25
    execute_on = timestep_end
  []
  [vel_z]
    type = NewmarkVelAux
    variable = vel_z
    acceleration = accel_z
    gamma = 0.5
    execute_on = timestep_end
  []
  [stress_yy]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_yy
    index_i = 1
    index_j = 1
  []
  [strain_yy]
    type = RankTwoAux
    rank_two_tensor = total_strain
    variable = strain_yy
    index_i = 1
    index_j = 1
  []
[]

[BCs]
  [top_x]
    type = DirichletBC
    variable = disp_x
    boundary = top
    value = 0.0
  []
  [top_y]
    type = DirichletBC
    variable = disp_y
    boundary = top
    value = 0.0
  []
  [top_z]
    type = DirichletBC
    variable = disp_z
    boundary = top
    value = 0.0
  []
  [Pressure]
    [Side1]
      boundary = bottom
      function = pressure
      factor = 1
      displacements = 'disp_x disp_y disp_z'
    []
  []
[]

[Materials]
  [Elasticity_tensor]
    type = ComputeElasticityTensor
    fill_method = symmetric_isotropic
    C_ijkl = '210 0'
  []

  [strain]
    type = ComputeSmallStrain
    displacements = 'disp_x disp_y disp_z'
  []

  [stress]
    type = ComputeLinearElasticStress
  []
  [density]
    type = GenericConstantMaterial
    prop_names = 'density'
    prop_values = '7750'
  []
[]

[Executioner]
  type = Transient
  start_time = 0
  end_time = 2
  dt = 0.1
[]

[Functions]
  [pressure]
    type = PiecewiseLinear
    x = '0.0 0.2 1.0 5.0'
    y = '0.0 0.2 1.0 1.0'
    scale_factor = 1e3
  []
[]

[Postprocessors]
  [dt]
    type = TimestepSize
  []
  [disp]
    type = NodalExtremeValue
    variable = disp_y
    boundary = bottom
  []
  [vel]
    type = NodalExtremeValue
    variable = vel_y
    boundary = bottom
  []
  [accel]
    type = NodalExtremeValue
    variable = accel_y
    boundary = bottom
  []
  [stress_yy]
    type = ElementAverageValue
    variable = stress_yy
  []
  [strain_yy]
    type = ElementAverageValue
    variable = strain_yy
  []

[]

[Outputs]
  exodus = true
  perf_graph = true
[]

(modules/tensor_mechanics/test/tests/plane_stress/weak_plane_stress_incremental.i)

[GlobalParams]
  order = FIRST
  family = LAGRANGE
  displacements = 'disp_x disp_y'
  temperature = temp
  out_of_plane_strain = strain_zz
[]

[Mesh]
  [./square]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 2
    ny = 2
  [../]
[]

[Variables]
  [./disp_x]
  [../]
  [./disp_y]
  [../]
  [./strain_zz]
  [../]
[]

[AuxVariables]
  [./temp]
  [../]
  [./nl_strain_zz]
    order = CONSTANT
    family = MONOMIAL
  [../]
[]

[Postprocessors]
  [./react_z]
    type = MaterialTensorIntegral
    rank_two_tensor = stress
    index_i = 2
    index_j = 2
  [../]
  [./min_strain_zz]
    type = NodalExtremeValue
    variable = strain_zz
    value_type = min
  [../]
  [./max_strain_zz]
    type = NodalExtremeValue
    variable = strain_zz
    value_type = max
  [../]
[]

[Modules/TensorMechanics/Master]
  [plane_stress]
    planar_formulation = WEAK_PLANE_STRESS
    strain = SMALL
    incremental = true
    generate_output = 'stress_xx stress_xy stress_yy stress_zz strain_xx strain_xy strain_yy'
    eigenstrain_names = eigenstrain
  []
[]

[AuxKernels]
  [./tempfuncaux]
    type = FunctionAux
    variable = temp
    function = tempfunc
    use_displaced_mesh = false
  [../]
  [./strain_zz]
    type = RankTwoAux
    rank_two_tensor = total_strain
    variable = nl_strain_zz
    index_i = 2
    index_j = 2
  [../]
[]

[Functions]
  [./pull]
    type = PiecewiseLinear
    x='0     1   100'
    y='0  0.00  0.00'
  [../]
  [./tempfunc]
    type = ParsedFunction
    value = '(1 - x) * t'
  [../]
[]

[BCs]
  [./bottomx]
    type = DirichletBC
    boundary = 0
    variable = disp_x
    value = 0.0
  [../]
  [./bottomy]
    type = DirichletBC
    boundary = 0
    variable = disp_y
    value = 0.0
  [../]
[]

[Materials]
  [./elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    poissons_ratio = 0.3
    youngs_modulus = 1e6
  [../]
  [./thermal_strain]
    type = ComputeThermalExpansionEigenstrain
    thermal_expansion_coeff = 0.02
    stress_free_temperature = 0.5
    eigenstrain_name = eigenstrain
  [../]
  [./stress]
    type = ComputeStrainIncrementBasedStress
  [../]
[]

[Executioner]
  type = Transient

  solve_type = PJFNK
  line_search = none

# controls for linear iterations
  l_max_its = 100
  l_tol = 1e-06

# controls for nonlinear iterations
  nl_max_its = 15
  nl_rel_tol = 1e-14
  nl_abs_tol = 1e-12

# time control
  start_time = 0.0
  dt = 1.0
  dtmin = 1.0
  end_time = 2.0
[]

[Outputs]
  exodus = true
[]

(modules/tensor_mechanics/test/tests/thermal_expansion_function/finite_const.i)

# This tests the thermal expansion coefficient function using both
# options to specify that function: mean and instantaneous.  There
# two blocks, each containing a single element, and these use the
# two variants of the function.

# In this test, the instantaneous CTE function has a constant value,
# while the mean CTE function is an analytic function designed to
# give the same response.  If \bar{alpha}(T) is the mean CTE function,
# and \alpha(T) is the instantaneous CTE function,

# \bar{\alpha}(T) = 1/(T-Tref) \intA^{T}_{Tsf} \alpha(T) dT

# where Tref is the reference temperature used to define the mean CTE
# function, and Tsf is the stress-free temperature.

# This version of the test uses finite deformation theory.
# The two models produce very similar results.  There are slight
# differences due to the large deformation treatment.

[Mesh]
  file = 'blocks.e'
[]

[GlobalParams]
  displacements = 'disp_x disp_y disp_z'
[]

[AuxVariables]
  [./temp]
    order = FIRST
    family = LAGRANGE
  [../]
[]

[Modules/TensorMechanics/Master]
  [./all]
    strain = FINITE
    add_variables = true
    eigenstrain_names = eigenstrain
    generate_output = 'strain_xx strain_yy strain_zz'
  [../]
[]

[BCs]
  [./left]
    type = DirichletBC
    variable = disp_x
    boundary = 3
    value = 0.0
  [../]

  [./bottom]
    type = DirichletBC
    variable = disp_y
    boundary = 2
    value = 0.0
  [../]

  [./back]
    type = DirichletBC
    variable = disp_z
    boundary = 1
    value = 0.0
  [../]
[]

[AuxKernels]
  [./temp]
    type = FunctionAux
    variable = temp
    block = '1 2'
    function = temp_func
  [../]
[]

[Materials]
  [./elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 1e6
    poissons_ratio = 0.3
  [../]
  [./small_stress]
    type = ComputeFiniteStrainElasticStress
  [../]
  [./thermal_expansion_strain1]
    type = ComputeMeanThermalExpansionFunctionEigenstrain
    block = 1
    thermal_expansion_function = cte_func_mean
    thermal_expansion_function_reference_temperature = 0.5
    stress_free_temperature = 0.0
    temperature = temp
    eigenstrain_name = eigenstrain
  [../]
  [./thermal_expansion_strain2]
    type = ComputeInstantaneousThermalExpansionFunctionEigenstrain
    block = 2
    thermal_expansion_function = cte_func_inst
    stress_free_temperature = 0.0
    temperature = temp
    eigenstrain_name = eigenstrain
  [../]
[]

[Functions]
  [./cte_func_mean]
    type = ParsedFunction
    vars = 'tsf tref scale' #stress free temp, reference temp, scale factor
    vals = '0.0 0.5  1e-4'
    value = 'scale * (t - tsf) / (t - tref)'
  [../]
  [./cte_func_inst]
    type = PiecewiseLinear
    xy_data = '0 1.0
               2 1.0'
    scale_factor = 1e-4
  [../]

  [./temp_func]
    type = PiecewiseLinear
    xy_data = '0 1
               1 2'
  [../]
[]

[Postprocessors]
  [./disp_1]
    type = NodalExtremeValue
    variable = disp_x
    boundary = 101
  [../]

  [./disp_2]
    type = NodalExtremeValue
    variable = disp_x
    boundary = 102
  [../]
[]

[Executioner]
  type = Transient

  solve_type = PJFNK
  l_max_its = 100
  l_tol = 1e-4
  nl_abs_tol = 1e-8
  nl_rel_tol = 1e-12

  start_time = 0.0
  end_time = 1.0
  dt = 0.1
[]

[Outputs]
  csv = true
[]

(modules/combined/test/tests/restart-transient-from-ss-with-stateful/sub_tr.i)

[GlobalParams]
  displacements = 'disp_x disp_y'
[]

[Mesh]
  [gen]
    type = GeneratedMeshGenerator
    nx = 8
    ny = 8
    xmin = -82.627
    xmax = 82.627
    ymin = -82.627
    ymax = 82.627
    dim = 2
  []
  [./extra_nodes_x]
    type = ExtraNodesetGenerator
    input = 'gen'
    new_boundary = 'no_x'
    coord = '0 82.627 0'
  [../]
  [./extra_nodes_y]
    type = ExtraNodesetGenerator
    input = 'extra_nodes_x'
    new_boundary = 'no_y'
    coord = '-82.627 0 0'
  [../]
[]

[Problem]
  type = ReferenceResidualProblem
  reference_vector = 'ref'
  extra_tag_vectors = 'ref'
[]

[AuxVariables]
  [./temp]
  [../]
[]

[Modules/TensorMechanics/Master]
  # FINITE strain when strain is large, i.e., visible movement.
  # SMALL strain when things are stressed, but may not move.
  [./fuel]
    add_variables = true
    strain = FINITE
    temperature = temp
    eigenstrain_names = 'thermal_eigenstrain'
    generate_output = 'vonmises_stress stress_xx stress_yy hydrostatic_stress max_principal_stress strain_xy elastic_strain_xx stress_xy'
    extra_vector_tags = 'ref'
    use_finite_deform_jacobian = true
    incremental = true
  [../]
[]

[BCs]
  [./no_x]
    type = DirichletBC
    variable = disp_x
    boundary = 'no_x'
    value = 0.0
    preset = true
  [../]
  [./no_y]
    type = DirichletBC
    preset = true
    variable = disp_y
    boundary = 'no_y'
    value = 0.0
  [../]
[]

[Materials]
  [./elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 3e10   # Pa
    poissons_ratio = 0.33    # unitless
  [../]
  [./thermal_strains]
    type = ComputeThermalExpansionEigenstrain
    temperature = temp
    thermal_expansion_coeff = 2e-6 # 1/K
    stress_free_temperature = 500 # K
    eigenstrain_name = 'thermal_eigenstrain'
  [../]
  [./stress_finite] # goes with FINITE strain formulation
    type = ComputeFiniteStrainElasticStress
  [../]
[]

[Postprocessors]
  [./avg_temp]
    type = ElementAverageValue
    variable = temp
    execute_on = 'initial timestep_end'
  [../]
  [./disp_x_max_element]
    type = ElementExtremeValue
    value_type = max
    variable = disp_x
    execute_on = 'initial timestep_end'
  [../]
  [./disp_y_max_element]
    type = ElementExtremeValue
    value_type = max
    variable = disp_y
    execute_on = 'initial timestep_end'
  [../]
  [./disp_x_max_nodal]
    type = NodalExtremeValue
    value_type = max
    variable = disp_x
    execute_on = 'initial timestep_end'
  [../]
  [./disp_y_max_nodal]
    type = NodalExtremeValue
    value_type = max
    variable = disp_y
    execute_on = 'initial timestep_end'
  [../]
[]

[Executioner]
  type = Transient
  petsc_options_iname = '-pc_type -pc_hypre_type -ksp_gmres_restart'
  petsc_options_value = 'hypre boomeramg 300'
  line_search = 'none'

  l_tol = 1e-02
  nl_rel_tol = 5e-04
  nl_abs_tol = 1e-2

  l_max_its = 50
  nl_max_its = 25

  start_time = 0
  end_time = 40
  dt = 10
[]

[Outputs]
  print_linear_residuals = false
  perf_graph = true
  exodus = true
[]

(modules/tensor_mechanics/test/tests/ad_thermal_expansion_function/small_const.i)

# This tests the thermal expansion coefficient function using both
# options to specify that function: mean and instantaneous.  There
# two blocks, each containing a single element, and these use the
# two variants of the function.

# In this test, the instantaneous CTE function has a constant value,
# while the mean CTE function is an analytic function designed to
# give the same response.  If \bar{alpha}(T) is the mean CTE function,
# and \alpha(T) is the instantaneous CTE function,

# \bar{\alpha}(T) = 1/(T-Tref) \intA^{T}_{Tsf} \alpha(T) dT

# where Tref is the reference temperature used to define the mean CTE
# function, and Tsf is the stress-free temperature.

# This version of the test uses small deformation theory.  The results
# from the two models are identical.

[Mesh]
  file = 'blocks.e'
[]

[GlobalParams]
  displacements = 'disp_x disp_y disp_z'
[]

[AuxVariables]
  [./temp]
    order = FIRST
    family = LAGRANGE
  [../]
[]

[Modules/TensorMechanics/Master]
  [./all]
    strain = SMALL
    add_variables = true
    eigenstrain_names = eigenstrain
    generate_output = 'strain_xx strain_yy strain_zz'
    use_automatic_differentiation = true
  [../]
[]

[BCs]
  [./left]
    type = DirichletBC
    variable = disp_x
    boundary = 3
    value = 0.0
  [../]

  [./bottom]
    type = DirichletBC
    variable = disp_y
    boundary = 2
    value = 0.0
  [../]

  [./back]
    type = DirichletBC
    variable = disp_z
    boundary = 1
    value = 0.0
  [../]
[]

[AuxKernels]
  [./temp]
    type = FunctionAux
    variable = temp
    block = '1 2'
    function = temp_func
  [../]
[]

[Materials]
  [./elasticity_tensor]
    type = ADComputeIsotropicElasticityTensor
    youngs_modulus = 1e6
    poissons_ratio = 0.3
  [../]
  [./small_stress]
    type = ADComputeLinearElasticStress
  [../]
  [./thermal_expansion_strain1]
    type = ADComputeMeanThermalExpansionFunctionEigenstrain
    block = 1
    thermal_expansion_function = cte_func_mean
    thermal_expansion_function_reference_temperature = 0.5
    stress_free_temperature = 0.0
    temperature = temp
    eigenstrain_name = eigenstrain
  [../]
  [./thermal_expansion_strain2]
    type = ADComputeInstantaneousThermalExpansionFunctionEigenstrain
    block = 2
    thermal_expansion_function = cte_func_inst
    stress_free_temperature = 0.0
    temperature = temp
    eigenstrain_name = eigenstrain
  [../]
[]

[Functions]
  [./cte_func_mean]
    type = ParsedFunction
    vars = 'tsf tref scale' #stress free temp, reference temp, scale factor
    vals = '0.0 0.5  1e-4'
    value = 'scale * (t - tsf) / (t - tref)'
  [../]
  [./cte_func_inst]
    type = PiecewiseLinear
    xy_data = '0 1.0
               2 1.0'
    scale_factor = 1e-4
  [../]

  [./temp_func]
    type = PiecewiseLinear
    xy_data = '0 1
               1 2'
  [../]
[]

[Postprocessors]
  [./disp_1]
    type = NodalExtremeValue
    variable = disp_x
    boundary = 101
  [../]

  [./disp_2]
    type = NodalExtremeValue
    variable = disp_x
    boundary = 102
  [../]
[]

[Executioner]
  type = Transient

  solve_type = NEWTON
  l_max_its = 100
  l_tol = 1e-4
  nl_abs_tol = 1e-8
  nl_rel_tol = 1e-12

  start_time = 0.0
  end_time = 1.0
  dt = 0.1
[]

[Outputs]
  csv = true
[]

(modules/tensor_mechanics/test/tests/dynamics/time_integration/hht_test_action.i)

# Test for  HHT time integration

# The test is for an 1D bar element of  unit length fixed on one end
# with a ramped pressure boundary condition applied to the other end.
# alpha, beta and gamma are HHT time integration parameters
# The equation of motion in terms of matrices is:
#
# M*accel + alpha*(K*disp - K*disp_old) + K*disp = P(t+alpha dt)*Area
#
# Here M is the mass matrix, K is the stiffness matrix, P is the applied pressure
#
# This equation is equivalent to:
#
# density*accel + alpha*(Div stress - Div stress_old) +Div Stress= P(t+alpha dt)
#
# The first term on the left is evaluated using the Inertial force kernel
# The next two terms on the left involving alpha are evaluated using the
# DynamicStressDivergenceTensors Kernel
# The residual due to Pressure is evaluated using Pressure boundary condition
#
# The system will come to steady state slowly after the pressure becomes constant.
# Alpha equal to zero will result in Newmark integration.

[Mesh]
  type = GeneratedMesh
  dim = 3
  nx = 1
  ny = 1
  nz = 1
  xmin = 0.0
  xmax = 0.1
  ymin = 0.0
  ymax = 1.0
  zmin = 0.0
  zmax = 0.1
[]

[GlobalParams]
  displacements = 'disp_x disp_y disp_z'
[]

[AuxVariables]
  [stress_yy]
    order = CONSTANT
    family = MONOMIAL
  []
  [strain_yy]
    order = CONSTANT
    family = MONOMIAL
  []

[]

[Modules/TensorMechanics/DynamicMaster]
  [all]
    add_variables = true
    hht_alpha = 0.11
    newmark_beta = 0.25
    newmark_gamma = 0.5
    density = 7750
  []
[]

[AuxKernels]
  [stress_yy]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_yy
    index_i = 0
    index_j = 1
  []
  [strain_yy]
    type = RankTwoAux
    rank_two_tensor = total_strain
    variable = strain_yy
    index_i = 0
    index_j = 1
  []
[]

[BCs]
  [top_y]
    type = DirichletBC
    variable = disp_y
    boundary = top
    value = 0.0
  []
  [top_x]
    type = DirichletBC
    variable = disp_x
    boundary = top
    value = 0.0
  []
  [top_z]
    type = DirichletBC
    variable = disp_z
    boundary = top
    value = 0.0
  []
  [bottom_x]
    type = DirichletBC
    variable = disp_x
    boundary = bottom
    value = 0.0
  []
  [bottom_z]
    type = DirichletBC
    variable = disp_z
    boundary = bottom
    value = 0.0
  []
  [Pressure]
    [Side1]
      boundary = bottom
      function = pressure
      factor = 1
      hht_alpha = 0.11
      displacements = 'disp_x disp_y disp_z'
    []
  []
[]

[Materials]
  [Elasticity_tensor]
    type = ComputeElasticityTensor
    block = 0
    fill_method = symmetric_isotropic
    C_ijkl = '210e9 0'
  []
  [stress]
    type = ComputeLinearElasticStress
    block = 0
  []
[]

[Executioner]
  type = Transient
  start_time = 0
  end_time = 2
  dt = 0.1
[]

[Functions]
  [pressure]
    type = PiecewiseLinear
    x = '0.0 0.1 0.2 1.0 2.0 5.0'
    y = '0.0 0.1 0.2 1.0 1.0 1.0'
    scale_factor = 1e9
  []
[]

[Postprocessors]
  [_dt]
    type = TimestepSize
  []
  [disp]
    type = NodalExtremeValue
    variable = disp_y
    boundary = bottom
  []
  [vel]
    type = NodalExtremeValue
    variable = vel_y
    boundary = bottom
  []
  [accel]
    type = NodalExtremeValue
    variable = accel_y
    boundary = bottom
  []
  [stress_yy]
    type = ElementAverageValue
    variable = stress_yy
  []
  [strain_yy]
    type = ElementAverageValue
    variable = strain_yy
  []
[]

[Outputs]
  exodus = true
  perf_graph = true
[]

(modules/tensor_mechanics/examples/coal_mining/fine.i)

# Strata deformation and fracturing around a coal mine - 3D model
#
# A "half model" is used.  The mine is 400m deep and
# just the roof is studied (-400<=z<=0).  The mining panel
# sits between 0<=x<=150, and 0<=y<=1000, so this simulates
# a coal panel that is 300m wide and 1000m long.  The outer boundaries
# are 1km from the excavation boundaries.
#
# Time is meaningless in this example
# as quasi-static solutions are sought at each timestep, but
# the number of timesteps controls the resolution of the
# process.
#
# The boundary conditions for this simulation are:
#  - disp_x = 0 at x=0 and x=1150
#  - disp_y = 0 at y=-1000 and y=1000
#  - disp_z = 0 at z=-400, but there is a time-dependent
#               Young's modulus that simulates excavation
#  - wc_x = 0 at y=-1000 and y=1000
#  - wc_y = 0 at x=0 and x=1150
# That is, rollers on the sides, free at top,
# and prescribed at bottom in the unexcavated portion.
#
# The small strain formulation is used.
#
# All stresses are measured in MPa.  The initial stress is consistent with
# the weight force from density 2500 kg/m^3, ie, stress_zz = 0.025*z MPa
# where gravity = 10 m.s^-2 = 1E-5 MPa m^2/kg.  The maximum and minimum
# principal horizontal stresses are assumed to be equal to 0.8*stress_zz.
#
# Material properties:
# Young's modulus = 8 GPa
# Poisson's ratio = 0.25
# Cosserat layer thickness = 1 m
# Cosserat-joint normal stiffness = large
# Cosserat-joint shear stiffness = 1 GPa
# MC cohesion = 3 MPa
# MC friction angle = 37 deg
# MC dilation angle = 8 deg
# MC tensile strength = 1 MPa
# MC compressive strength = 100 MPa
# WeakPlane cohesion = 0.1 MPa
# WeakPlane friction angle = 30 deg
# WeakPlane dilation angle = 10 deg
# WeakPlane tensile strength = 0.1 MPa
# WeakPlane compressive strength = 100 MPa softening to 1 MPa at strain = 1
#
[Mesh]
  [file]
    type = FileMeshGenerator
    file = mesh/fine.e
  []
  [./xmin]
    input = file
    type = SideSetsAroundSubdomainGenerator
    block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30'
    new_boundary = xmin
    normal = '-1 0 0'
  [../]
  [./xmax]
    input = xmin
    type = SideSetsAroundSubdomainGenerator
    block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30'
    new_boundary = xmax
    normal = '1 0 0'
  [../]
  [./ymin]
    input = xmax
    type = SideSetsAroundSubdomainGenerator
    block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30'
    new_boundary = ymin
    normal = '0 -1 0'
  [../]
  [./ymax]
    input = ymin
    type = SideSetsAroundSubdomainGenerator
    block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30'
    new_boundary = ymax
    normal = '0 1 0'
  [../]
  [./zmax]
    input = ymax
    type = SideSetsAroundSubdomainGenerator
    block = 30
    new_boundary = zmax
    normal = '0 0 1'
  [../]
  [./zmin]
    input = zmax
    type = SideSetsAroundSubdomainGenerator
    block = 2
    new_boundary = zmin
    normal = '0 0 -1'
  [../]
  [./excav]
    type = SubdomainBoundingBoxGenerator
    input = zmin
    block_id = 1
    bottom_left = '0 0 -400'
    top_right = '150 1000 -397'
  [../]
  [./roof]
    type = SideSetsAroundSubdomainGenerator
    block = 1
    input = excav
    new_boundary = roof
    normal = '0 0 1'
  [../]
[]

[GlobalParams]
  perform_finite_strain_rotations = false
  displacements = 'disp_x disp_y disp_z'
  Cosserat_rotations = 'wc_x wc_y wc_z'
[]

[Variables]
  [./disp_x]
  [../]
  [./disp_y]
  [../]
  [./disp_z]
  [../]
  [./wc_x]
  [../]
  [./wc_y]
  [../]
[]

[Kernels]
  [./cx_elastic]
    type = CosseratStressDivergenceTensors
    use_displaced_mesh = false
    variable = disp_x
    component = 0
  [../]
  [./cy_elastic]
    type = CosseratStressDivergenceTensors
    use_displaced_mesh = false
    variable = disp_y
    component = 1
  [../]
  [./cz_elastic]
    type = CosseratStressDivergenceTensors
    use_displaced_mesh = false
    variable = disp_z
    component = 2
  [../]
  [./x_couple]
    type = StressDivergenceTensors
    use_displaced_mesh = false
    variable = wc_x
    displacements = 'wc_x wc_y wc_z'
    component = 0
    base_name = couple
  [../]
  [./y_couple]
    type = StressDivergenceTensors
    use_displaced_mesh = false
    variable = wc_y
    displacements = 'wc_x wc_y wc_z'
    component = 1
    base_name = couple
  [../]
  [./x_moment]
    type = MomentBalancing
    use_displaced_mesh = false
    variable = wc_x
    component = 0
  [../]
  [./y_moment]
    type = MomentBalancing
    use_displaced_mesh = false
    variable = wc_y
    component = 1
  [../]
  [./gravity]
    type = Gravity
    use_displaced_mesh = false
    variable = disp_z
    value = -10E-6 # remember this is in MPa
  [../]
[]


[AuxVariables]
  [./wc_z]
  [../]
  [./stress_xx]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./stress_xy]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./stress_xz]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./stress_yx]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./stress_yy]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./stress_yz]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./stress_zx]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./stress_zy]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./stress_zz]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./mc_shear]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./mc_tensile]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./wp_shear]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./wp_tensile]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./wp_shear_f]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./wp_tensile_f]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./mc_shear_f]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./mc_tensile_f]
    order = CONSTANT
    family = MONOMIAL
  [../]
[]

[AuxKernels]
  [./stress_xx]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_xx
    index_i = 0
    index_j = 0
  [../]
  [./stress_xy]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_xy
    index_i = 0
    index_j = 1
  [../]
  [./stress_xz]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_xz
    index_i = 0
    index_j = 2
  [../]
  [./stress_yx]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_yx
    index_i = 1
    index_j = 0
  [../]
  [./stress_yy]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_yy
    index_i = 1
    index_j = 1
  [../]
  [./stress_yz]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_yz
    index_i = 1
    index_j = 2
  [../]
  [./stress_zx]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_zx
    index_i = 2
    index_j = 0
  [../]
  [./stress_zy]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_zy
    index_i = 2
    index_j = 1
  [../]
  [./stress_zz]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_zz
    index_i = 2
    index_j = 2
  [../]
  [./mc_shear]
    type = MaterialStdVectorAux
    index = 0
    property = mc_plastic_internal_parameter
    variable = mc_shear
  [../]
  [./mc_tensile]
    type = MaterialStdVectorAux
    index = 1
    property = mc_plastic_internal_parameter
    variable = mc_tensile
  [../]
  [./wp_shear]
    type = MaterialStdVectorAux
    index = 0
    property = wp_plastic_internal_parameter
    variable = wp_shear
  [../]
  [./wp_tensile]
    type = MaterialStdVectorAux
    index = 1
    property = wp_plastic_internal_parameter
    variable = wp_tensile
  [../]
  [./mc_shear_f]
    type = MaterialStdVectorAux
    index = 6
    property = mc_plastic_yield_function
    variable = mc_shear_f
  [../]
  [./mc_tensile_f]
    type = MaterialStdVectorAux
    index = 0
    property = mc_plastic_yield_function
    variable = mc_tensile_f
  [../]
  [./wp_shear_f]
    type = MaterialStdVectorAux
    index = 0
    property = wp_plastic_yield_function
    variable = wp_shear_f
  [../]
  [./wp_tensile_f]
    type = MaterialStdVectorAux
    index = 1
    property = wp_plastic_yield_function
    variable = wp_tensile_f
  [../]
[]



[BCs]
  [./no_x]
    type = DirichletBC
    variable = disp_x
    boundary = 'xmin xmax'
    value = 0.0
  [../]
  [./no_y]
    type = DirichletBC
    variable = disp_y
    boundary = 'ymin ymax'
    value = 0.0
  [../]
  [./no_z]
    type = DirichletBC
    variable = disp_z
    boundary = zmin
    value = 0.0
  [../]
  [./no_wc_x]
    type = DirichletBC
    variable = wc_x
    boundary = 'ymin ymax'
    value = 0.0
  [../]
  [./no_wc_y]
    type = DirichletBC
    variable = wc_y
    boundary = 'xmin xmax'
    value = 0.0
  [../]
  [./roof]
    type = StickyBC
    variable = disp_z
    min_value = -3.0
    boundary = roof
  [../]
[]

[Functions]
  [./ini_xx]
    type = ParsedFunction
    value = '0.8*2500*10E-6*z'
  [../]
  [./ini_zz]
    type = ParsedFunction
    value = '2500*10E-6*z'
  [../]
  [./excav_sideways]
    type = ParsedFunction
    vars = 'end_t ymin ymax  minval maxval slope'
    vals = '100.0   0    1000.0 1E-9 1 10'
    # excavation face at ymin+(ymax-ymin)*min(t/end_t,1)
    # slope is the distance over which the modulus reduces from maxval to minval
    value = 'if(y<ymin+(ymax-ymin)*min(t/end_t,1),minval,if(y<ymin+(ymax-ymin)*min(t/end_t,1)+slope,minval+(maxval-minval)*(y-(ymin+(ymax-ymin)*min(t/end_t,1)))/slope,maxval))'
  [../]
  [./density_sideways]
    type = ParsedFunction
    vars = 'end_t ymin ymax  minval maxval'
    vals = '100.0   0    1000.0 0 2500'
    value = 'if(y<ymin+(ymax-ymin)*min(t/end_t,1),minval,maxval)'
  [../]
[]

[UserObjects]
  [./mc_coh_strong_harden]
    type = TensorMechanicsHardeningExponential
    value_0 = 2.99 # MPa
    value_residual = 3.01 # MPa
    rate = 1.0
  [../]
  [./mc_fric]
    type = TensorMechanicsHardeningConstant
    value = 0.65 # 37deg
  [../]
  [./mc_dil]
    type = TensorMechanicsHardeningConstant
    value = 0.15 # 8deg
  [../]

  [./mc_tensile_str_strong_harden]
    type = TensorMechanicsHardeningExponential
    value_0 = 1.0 # MPa
    value_residual = 1.0 # MPa
    rate = 1.0
  [../]
  [./mc_compressive_str]
    type = TensorMechanicsHardeningCubic
    value_0 = 100 # Large!
    value_residual = 100
    internal_limit = 0.1
  [../]

  [./wp_coh_harden]
    type = TensorMechanicsHardeningCubic
    value_0 = 0.1
    value_residual = 0.1
    internal_limit = 10
  [../]
  [./wp_tan_fric]
    type = TensorMechanicsHardeningConstant
    value = 0.36 # 20deg
  [../]
  [./wp_tan_dil]
    type = TensorMechanicsHardeningConstant
    value = 0.18 # 10deg
  [../]

  [./wp_tensile_str_harden]
    type = TensorMechanicsHardeningCubic
    value_0 = 0.1
    value_residual = 0.1
    internal_limit = 10
  [../]
  [./wp_compressive_str_soften]
    type = TensorMechanicsHardeningCubic
    value_0 = 100
    value_residual = 1
    internal_limit = 1.0
  [../]
[]

[Materials]
  [./elasticity_tensor_0]
    type = ComputeLayeredCosseratElasticityTensor
    block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30'
    young = 8E3 # MPa
    poisson = 0.25
    layer_thickness = 1.0
    joint_normal_stiffness = 1E9 # huge
    joint_shear_stiffness = 1E3 # MPa
  [../]
  [./elasticity_tensor_1]
    type = ComputeLayeredCosseratElasticityTensor
    block = 1
    young = 8E3 # MPa
    poisson = 0.25
    layer_thickness = 1.0
    joint_normal_stiffness = 1E9 # huge
    joint_shear_stiffness = 1E3 # MPa
    elasticity_tensor_prefactor = excav_sideways
  [../]
  [./strain]
    type = ComputeCosseratIncrementalSmallStrain
    eigenstrain_names = ini_stress
  [../]
  [./ini_stress]
    type = ComputeEigenstrainFromInitialStress
    eigenstrain_name = ini_stress
    initial_stress = 'ini_xx 0 0  0 ini_xx 0  0 0 ini_zz'
  [../]

  [./stress_0]
    type = ComputeMultipleInelasticCosseratStress
    block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30'
    inelastic_models = 'mc wp'
    cycle_models = true
    relative_tolerance = 2.0
    absolute_tolerance = 1E6
    max_iterations = 1
    tangent_operator = nonlinear
    perform_finite_strain_rotations = false
  [../]
  [./stress_1]
    type = ComputeMultipleInelasticCosseratStress
    block = 1
    inelastic_models = ''
    relative_tolerance = 2.0
    absolute_tolerance = 1E6
    max_iterations = 1
    tangent_operator = nonlinear
    perform_finite_strain_rotations = false
  [../]
  [./mc]
    type = CappedMohrCoulombCosseratStressUpdate
    warn_about_precision_loss = false
    host_youngs_modulus = 8E3
    host_poissons_ratio = 0.25
    base_name = mc
    tensile_strength = mc_tensile_str_strong_harden
    compressive_strength = mc_compressive_str
    cohesion = mc_coh_strong_harden
    friction_angle = mc_fric
    dilation_angle = mc_dil
    max_NR_iterations = 100000
    smoothing_tol = 0.1 # MPa  # Must be linked to cohesion
    yield_function_tol = 1E-9 # MPa.  this is essentially the lowest possible without lots of precision loss
    perfect_guess = true
    min_step_size = 1.0
  [../]
  [./wp]
    type = CappedWeakPlaneCosseratStressUpdate
    warn_about_precision_loss = false
    base_name = wp
    cohesion = wp_coh_harden
    tan_friction_angle = wp_tan_fric
    tan_dilation_angle = wp_tan_dil
    tensile_strength = wp_tensile_str_harden
    compressive_strength = wp_compressive_str_soften
    max_NR_iterations = 10000
    tip_smoother = 0.1
    smoothing_tol = 0.1 # MPa  # Note, this must be tied to cohesion, otherwise get no possible return at cone apex
    yield_function_tol = 1E-11 # MPa.  this is essentially the lowest possible without lots of precision loss
    perfect_guess = true
    min_step_size = 1.0E-3
  [../]


  [./density_0]
    type = GenericConstantMaterial
    block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30'
    prop_names = density
    prop_values = 2500
  [../]
  [./density_1]
    type = GenericFunctionMaterial
    block = 1
    prop_names = density
    prop_values = density_sideways
  [../]
[]

[Preconditioning]
  [./SMP]
    type = SMP
    full = true
  []
[]

[Postprocessors]
  [./min_roof_disp]
    type = NodalExtremeValue
    boundary = roof
    value_type = min
    variable = disp_z
  [../]
  [./min_surface_disp]
    type = NodalExtremeValue
    boundary = zmax
    value_type = min
    variable = disp_z
  [../]
[]

[Executioner]
  type = Transient

  solve_type = 'NEWTON'

  petsc_options = '-snes_converged_reason'
  petsc_options_iname = '-pc_type -ksp_type -ksp_gmres_restart'
  petsc_options_value = ' bjacobi  gmres     200'

  line_search = bt

  nl_abs_tol = 1e-3
  nl_rel_tol = 1e-5

  l_max_its = 30
  nl_max_its = 1000

  start_time = 0.0
  dt = 0.5
  end_time = 100.0

[]

[Outputs]
  interval = 1
  print_linear_residuals = false
  exodus = true
  csv = true
  console = true
[]

(test/tests/userobjects/Terminator/terminator.i)

###########################################################
# This is a test of the UserObject System. The
# Terminator UserObject executes independently after
# each solve and can terminate the solve early due to
# user-defined criteria. (Type: GeneralUserObject)
#
# @Requirement F6.40
###########################################################


[Mesh]
  type = GeneratedMesh
  dim = 2
  nx = 30
  ny = 6
  xmin = -15.0
  xmax = 15.0
  ymin = -3.0
  ymax = 3.0
  elem_type = QUAD4
[]

[Variables]
  [./c]
    order = FIRST
    family = LAGRANGE
    initial_condition = 1
  [../]
[]

[Postprocessors]
  [./max_c]
    type = NodalExtremeValue
    variable = c
    execute_on = 'initial timestep_end'
  [../]
[]

[UserObjects]
  [./arnold]
    type = Terminator
    expression = 'max_c < 0.5'
  [../]
[]

[Kernels]
  [./cres]
    type = Diffusion
    variable = c
  [../]

  [./time]
    type = TimeDerivative
    variable = c
  [../]
[]

[BCs]
  [./c]
    type = DirichletBC
    variable = c
    boundary = left
    value = 0
  [../]
[]

[Executioner]
  type = Transient
  dt = 100
  num_steps = 6
[]

[Outputs]
  exodus = true
[]

(modules/porous_flow/examples/coal_mining/fine_with_fluid.i)

#################################################################
#
#  NOTE:
#  The mesh for this model is too large for the MOOSE repository
#  so is kept in the the large_media submodule
#
#################################################################
#
# Strata deformation and fluid flow aaround a coal mine - 3D model
#
# A "half model" is used.  The mine is 400m deep and
# just the roof is studied (-400<=z<=0).  The mining panel
# sits between 0<=x<=150, and 0<=y<=1000, so this simulates
# a coal panel that is 300m wide and 1000m long.  The outer boundaries
# are 1km from the excavation boundaries.
#
# The excavation takes 0.5 years.
#
# The boundary conditions for this simulation are:
#  - disp_x = 0 at x=0 and x=1150
#  - disp_y = 0 at y=-1000 and y=1000
#  - disp_z = 0 at z=-400, but there is a time-dependent
#               Young modulus that simulates excavation
#  - wc_x = 0 at y=-1000 and y=1000
#  - wc_y = 0 at x=0 and x=1150
#  - no flow at x=0, z=-400 and z=0
#  - fixed porepressure at y=-1000, y=1000 and x=1150
# That is, rollers on the sides, free at top,
# and prescribed at bottom in the unexcavated portion.
#
# A single-phase unsaturated fluid is used.
#
# The small strain formulation is used.
#
# All stresses are measured in MPa, and time units are measured in years.
#
# The initial porepressure is hydrostatic with P=0 at z=0, so
# Porepressure ~ - 0.01*z MPa, where the fluid has density 1E3 kg/m^3 and
# gravity = = 10 m.s^-2 = 1E-5 MPa m^2/kg.
# To be more accurate, i use
# Porepressure = -bulk * log(1 + g*rho0*z/bulk)
# where bulk=2E3 MPa and rho0=1Ee kg/m^3.
# The initial stress is consistent with the weight force from undrained
# density 2500 kg/m^3, and fluid porepressure, and a Biot coefficient of 0.7, ie,
# stress_zz^effective = 0.025*z + 0.7 * initial_porepressure
# The maximum and minimum principal horizontal effective stresses are
# assumed to be equal to 0.8*stress_zz.
#
# Material properties:
# Young's modulus = 8 GPa
# Poisson's ratio = 0.25
# Cosserat layer thickness = 1 m
# Cosserat-joint normal stiffness = large
# Cosserat-joint shear stiffness = 1 GPa
# MC cohesion = 2 MPa
# MC friction angle = 35 deg
# MC dilation angle = 8 deg
# MC tensile strength = 1 MPa
# MC compressive strength = 100 MPa
# WeakPlane cohesion = 0.1 MPa
# WeakPlane friction angle = 30 deg
# WeakPlane dilation angle = 10 deg
# WeakPlane tensile strength = 0.1 MPa
# WeakPlane compressive strength = 100 MPa softening to 1 MPa at strain = 1
# Fluid density at zero porepressure = 1E3 kg/m^3
# Fluid bulk modulus = 2E3 MPa
# Fluid viscosity = 1.1E-3 Pa.s = 1.1E-9 MPa.s = 3.5E-17 MPa.year
#
[GlobalParams]
  perform_finite_strain_rotations = false
  displacements = 'disp_x disp_y disp_z'
  Cosserat_rotations = 'wc_x wc_y wc_z'
  PorousFlowDictator = dictator
  biot_coefficient = 0.7
[]

[Mesh]
  [file]
    type = FileMeshGenerator
    file = fine.e
  []
  [xmin]
    type = SideSetsAroundSubdomainGenerator
     block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30'
    new_boundary = xmin
    normal = '-1 0 0'
    input = file
  []
  [xmax]
    type = SideSetsAroundSubdomainGenerator
     block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30'
    new_boundary = xmax
    normal = '1 0 0'
    input = xmin
  []
  [ymin]
    type = SideSetsAroundSubdomainGenerator
     block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30'
    new_boundary = ymin
    normal = '0 -1 0'
    input = xmax
  []
  [ymax]
    type = SideSetsAroundSubdomainGenerator
     block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30'
    new_boundary = ymax
    normal = '0 1 0'
    input = ymin
  []
  [zmax]
    type = SideSetsAroundSubdomainGenerator
    block = 30
    new_boundary = zmax
    normal = '0 0 1'
    input = ymax
  []
  [zmin]
    type = SideSetsAroundSubdomainGenerator
    block = 2
    new_boundary = zmin
    normal = '0 0 -1'
    input = zmax
  []
  [excav]
    type = SubdomainBoundingBoxGenerator
    input = zmin
    block_id = 1
    bottom_left = '0 0 -400'
    top_right = '150 1000 -397'
  []
  [roof]
    type = SideSetsBetweenSubdomainsGenerator
    primary_block = 3
    paired_block = 1
    input = excav
    new_boundary = roof
  []
[]

[Variables]
  [disp_x]
  []
  [disp_y]
  []
  [disp_z]
  []
  [wc_x]
  []
  [wc_y]
  []
  [porepressure]
    scaling = 1E-5
  []
[]

[ICs]
  [porepressure]
    type = FunctionIC
    variable = porepressure
    function = ini_pp
  []
[]

[Kernels]
  [cx_elastic]
    type = CosseratStressDivergenceTensors
    use_displaced_mesh = false
    variable = disp_x
    component = 0
  []
  [cy_elastic]
    type = CosseratStressDivergenceTensors
    use_displaced_mesh = false
    variable = disp_y
    component = 1
  []
  [cz_elastic]
    type = CosseratStressDivergenceTensors
    use_displaced_mesh = false
    variable = disp_z
    component = 2
  []
  [x_couple]
    type = StressDivergenceTensors
    use_displaced_mesh = false
    variable = wc_x
    displacements = 'wc_x wc_y wc_z'
    component = 0
    base_name = couple
  []
  [y_couple]
    type = StressDivergenceTensors
    use_displaced_mesh = false
    variable = wc_y
    displacements = 'wc_x wc_y wc_z'
    component = 1
    base_name = couple
  []
  [x_moment]
    type = MomentBalancing
    use_displaced_mesh = false
    variable = wc_x
    component = 0
  []
  [y_moment]
    type = MomentBalancing
    use_displaced_mesh = false
    variable = wc_y
    component = 1
  []
  [gravity]
    type = Gravity
    use_displaced_mesh = false
    variable = disp_z
    value = -10E-6 # remember this is in MPa
  []
  [poro_x]
    type = PorousFlowEffectiveStressCoupling
    use_displaced_mesh = false
    variable = disp_x
    component = 0
  []
  [poro_y]
    type = PorousFlowEffectiveStressCoupling
    use_displaced_mesh = false
    variable = disp_y
    component = 1
  []
  [poro_z]
    type = PorousFlowEffectiveStressCoupling
    use_displaced_mesh = false
    component = 2
    variable = disp_z
  []
  [poro_vol_exp]
    type = PorousFlowMassVolumetricExpansion
    use_displaced_mesh = false
     block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30'
    variable = porepressure
    fluid_component = 0
  []
  [mass0]
    type = PorousFlowMassTimeDerivative
    use_displaced_mesh = false
    fluid_component = 0
    variable = porepressure
  []
  [flux]
    type = PorousFlowAdvectiveFlux
    use_displaced_mesh = false
    variable = porepressure
    gravity = '0 0 -10E-6'
    fluid_component = 0
  []
[]


[AuxVariables]
  [saturation]
    order = CONSTANT
    family = MONOMIAL
  []
  [darcy_x]
    order = CONSTANT
    family = MONOMIAL
  []
  [darcy_y]
    order = CONSTANT
    family = MONOMIAL
  []
  [darcy_z]
    order = CONSTANT
    family = MONOMIAL
  []
  [porosity]
    order = CONSTANT
    family = MONOMIAL
  []
  [wc_z]
  []
  [stress_xx]
    order = CONSTANT
    family = MONOMIAL
  []
  [stress_xy]
    order = CONSTANT
    family = MONOMIAL
  []
  [stress_xz]
    order = CONSTANT
    family = MONOMIAL
  []
  [stress_yx]
    order = CONSTANT
    family = MONOMIAL
  []
  [stress_yy]
    order = CONSTANT
    family = MONOMIAL
  []
  [stress_yz]
    order = CONSTANT
    family = MONOMIAL
  []
  [stress_zx]
    order = CONSTANT
    family = MONOMIAL
  []
  [stress_zy]
    order = CONSTANT
    family = MONOMIAL
  []
  [stress_zz]
    order = CONSTANT
    family = MONOMIAL
  []
  [total_strain_xx]
    order = CONSTANT
    family = MONOMIAL
  []
  [total_strain_xy]
    order = CONSTANT
    family = MONOMIAL
  []
  [total_strain_xz]
    order = CONSTANT
    family = MONOMIAL
  []
  [total_strain_yx]
    order = CONSTANT
    family = MONOMIAL
  []
  [total_strain_yy]
    order = CONSTANT
    family = MONOMIAL
  []
  [total_strain_yz]
    order = CONSTANT
    family = MONOMIAL
  []
  [total_strain_zx]
    order = CONSTANT
    family = MONOMIAL
  []
  [total_strain_zy]
    order = CONSTANT
    family = MONOMIAL
  []
  [total_strain_zz]
    order = CONSTANT
    family = MONOMIAL
  []
  [perm_xx]
    order = CONSTANT
    family = MONOMIAL
  []
  [perm_yy]
    order = CONSTANT
    family = MONOMIAL
  []
  [perm_zz]
    order = CONSTANT
    family = MONOMIAL
  []
  [mc_shear]
    order = CONSTANT
    family = MONOMIAL
  []
  [mc_tensile]
    order = CONSTANT
    family = MONOMIAL
  []
  [wp_shear]
    order = CONSTANT
    family = MONOMIAL
  []
  [wp_tensile]
    order = CONSTANT
    family = MONOMIAL
  []
  [wp_shear_f]
    order = CONSTANT
    family = MONOMIAL
  []
  [wp_tensile_f]
    order = CONSTANT
    family = MONOMIAL
  []
  [mc_shear_f]
    order = CONSTANT
    family = MONOMIAL
  []
  [mc_tensile_f]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[AuxKernels]
  [saturation_water]
    type = PorousFlowPropertyAux
    variable = saturation
    property = saturation
    phase = 0
    execute_on = timestep_end
  []
  [darcy_x]
    type = PorousFlowDarcyVelocityComponent
    variable = darcy_x
    gravity = '0 0 -10E-6'
    component = x
  []
  [darcy_y]
    type = PorousFlowDarcyVelocityComponent
    variable = darcy_y
    gravity = '0 0 -10E-6'
    component = y
  []
  [darcy_z]
    type = PorousFlowDarcyVelocityComponent
    variable = darcy_z
    gravity = '0 0 -10E-6'
    component = z
  []
  [porosity]
    type = PorousFlowPropertyAux
    property = porosity
    variable = porosity
    execute_on = timestep_end
  []
  [stress_xx]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_xx
    index_i = 0
    index_j = 0
    execute_on = timestep_end
  []
  [stress_xy]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_xy
    index_i = 0
    index_j = 1
    execute_on = timestep_end
  []
  [stress_xz]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_xz
    index_i = 0
    index_j = 2
    execute_on = timestep_end
  []
  [stress_yx]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_yx
    index_i = 1
    index_j = 0
    execute_on = timestep_end
  []
  [stress_yy]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_yy
    index_i = 1
    index_j = 1
    execute_on = timestep_end
  []
  [stress_yz]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_yz
    index_i = 1
    index_j = 2
    execute_on = timestep_end
  []
  [stress_zx]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_zx
    index_i = 2
    index_j = 0
    execute_on = timestep_end
  []
  [stress_zy]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_zy
    index_i = 2
    index_j = 1
    execute_on = timestep_end
  []
  [stress_zz]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_zz
    index_i = 2
    index_j = 2
    execute_on = timestep_end
  []
  [total_strain_xx]
    type = RankTwoAux
    rank_two_tensor = total_strain
    variable = total_strain_xx
    index_i = 0
    index_j = 0
    execute_on = timestep_end
  []
  [total_strain_xy]
    type = RankTwoAux
    rank_two_tensor = total_strain
    variable = total_strain_xy
    index_i = 0
    index_j = 1
    execute_on = timestep_end
  []
  [total_strain_xz]
    type = RankTwoAux
    rank_two_tensor = total_strain
    variable = total_strain_xz
    index_i = 0
    index_j = 2
    execute_on = timestep_end
  []
  [total_strain_yx]
    type = RankTwoAux
    rank_two_tensor = total_strain
    variable = total_strain_yx
    index_i = 1
    index_j = 0
    execute_on = timestep_end
  []
  [total_strain_yy]
    type = RankTwoAux
    rank_two_tensor = total_strain
    variable = total_strain_yy
    index_i = 1
    index_j = 1
    execute_on = timestep_end
  []
  [total_strain_yz]
    type = RankTwoAux
    rank_two_tensor = total_strain
    variable = total_strain_yz
    index_i = 1
    index_j = 2
    execute_on = timestep_end
  []
  [total_strain_zx]
    type = RankTwoAux
    rank_two_tensor = total_strain
    variable = total_strain_zx
    index_i = 2
    index_j = 0
    execute_on = timestep_end
  []
  [total_strain_zy]
    type = RankTwoAux
    rank_two_tensor = total_strain
    variable = total_strain_zy
    index_i = 2
    index_j = 1
    execute_on = timestep_end
  []
  [total_strain_zz]
    type = RankTwoAux
    rank_two_tensor = total_strain
    variable = total_strain_zz
    index_i = 2
    index_j = 2
    execute_on = timestep_end
  []
  [perm_xx]
    type = PorousFlowPropertyAux
    property = permeability
    variable = perm_xx
    row = 0
    column = 0
    execute_on = timestep_end
  []
  [perm_yy]
    type = PorousFlowPropertyAux
    property = permeability
    variable = perm_yy
    row = 1
    column = 1
    execute_on = timestep_end
  []
  [perm_zz]
    type = PorousFlowPropertyAux
    property = permeability
    variable = perm_zz
    row = 2
    column = 2
    execute_on = timestep_end
  []
  [mc_shear]
    type = MaterialStdVectorAux
    index = 0
    property = mc_plastic_internal_parameter
    variable = mc_shear
    execute_on = timestep_end
  []
  [mc_tensile]
    type = MaterialStdVectorAux
    index = 1
    property = mc_plastic_internal_parameter
    variable = mc_tensile
    execute_on = timestep_end
  []
  [wp_shear]
    type = MaterialStdVectorAux
    index = 0
    property = wp_plastic_internal_parameter
    variable = wp_shear
    execute_on = timestep_end
  []
  [wp_tensile]
    type = MaterialStdVectorAux
    index = 1
    property = wp_plastic_internal_parameter
    variable = wp_tensile
    execute_on = timestep_end
  []
  [mc_shear_f]
    type = MaterialStdVectorAux
    index = 6
    property = mc_plastic_yield_function
    variable = mc_shear_f
    execute_on = timestep_end
  []
  [mc_tensile_f]
    type = MaterialStdVectorAux
    index = 0
    property = mc_plastic_yield_function
    variable = mc_tensile_f
    execute_on = timestep_end
  []
  [wp_shear_f]
    type = MaterialStdVectorAux
    index = 0
    property = wp_plastic_yield_function
    variable = wp_shear_f
    execute_on = timestep_end
  []
  [wp_tensile_f]
    type = MaterialStdVectorAux
    index = 1
    property = wp_plastic_yield_function
    variable = wp_tensile_f
    execute_on = timestep_end
  []
[]



[BCs]
  [no_x]
    type = DirichletBC
    variable = disp_x
    boundary = 'xmin xmax'
    value = 0.0
  []
  [no_y]
    type = DirichletBC
    variable = disp_y
    boundary = 'ymin ymax'
    value = 0.0
  []
  [no_z]
    type = DirichletBC
    variable = disp_z
    boundary = zmin
    value = 0.0
  []
  [no_wc_x]
    type = DirichletBC
    variable = wc_x
    boundary = 'ymin ymax'
    value = 0.0
  []
  [no_wc_y]
    type = DirichletBC
    variable = wc_y
    boundary = 'xmin xmax'
    value = 0.0
  []
  [fix_porepressure]
    type = FunctionDirichletBC
    variable = porepressure
    boundary = 'ymin ymax xmax'
    function = ini_pp
  []
  [roof_porepressure]
    type = PorousFlowPiecewiseLinearSink
    variable = porepressure
    pt_vals = '-1E3 1E3'
    multipliers = '-1 1'
    fluid_phase = 0
    flux_function = roof_conductance
    boundary = roof
  []
  [roof]
    type = StickyBC
    variable = disp_z
    min_value = -3.0
    boundary = roof
  []
[]

[Functions]
  [ini_pp]
    type = ParsedFunction
    vars = 'bulk p0 g    rho0'
    vals = '2E3 0.0 1E-5 1E3'
    value = '-bulk*log(exp(-p0/bulk)+g*rho0*z/bulk)'
  []
  [ini_xx]
    type = ParsedFunction
    vars = 'bulk p0 g    rho0 biot'
    vals = '2E3 0.0 1E-5 1E3  0.7'
    value = '0.8*(2500*10E-6*z+biot*(-bulk*log(exp(-p0/bulk)+g*rho0*z/bulk)))'
  []
  [ini_zz]
    type = ParsedFunction
    vars = 'bulk p0 g    rho0 biot'
    vals = '2E3 0.0 1E-5 1E3  0.7'
    value = '2500*10E-6*z+biot*(-bulk*log(exp(-p0/bulk)+g*rho0*z/bulk))'
  []
  [excav_sideways]
    type = ParsedFunction
    vars = 'end_t ymin ymax  minval maxval slope'
    vals = '0.5   0    1000.0 1E-9 1 10'
    # excavation face at ymin+(ymax-ymin)*min(t/end_t,1)
    # slope is the distance over which the modulus reduces from maxval to minval
    value = 'if(y<ymin+(ymax-ymin)*min(t/end_t,1),minval,if(y<ymin+(ymax-ymin)*min(t/end_t,1)+slope,minval+(maxval-minval)*(y-(ymin+(ymax-ymin)*min(t/end_t,1)))/slope,maxval))'
  []
  [density_sideways]
    type = ParsedFunction
    vars = 'end_t ymin ymax  minval maxval'
    vals = '0.5   0    1000.0 0 2500'
    value = 'if(y<ymin+(ymax-ymin)*min(t/end_t,1),minval,maxval)'
  []
  [roof_conductance]
    type = ParsedFunction
    vars = 'end_t ymin ymax   maxval minval'
    vals = '0.5   0    1000.0 1E7      0'
    value = 'if(y<ymin+(ymax-ymin)*min(t/end_t,1),maxval,minval)'
  []
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'porepressure disp_x disp_y disp_z'
    number_fluid_phases = 1
    number_fluid_components = 1
  []
  [pc]
    type = PorousFlowCapillaryPressureVG
    m = 0.5
    alpha = 1 # MPa^-1
  []

  [mc_coh_strong_harden]
    type = TensorMechanicsHardeningExponential
    value_0 = 1.99 # MPa
    value_residual = 2.01 # MPa
    rate = 1.0
  []
  [mc_fric]
    type = TensorMechanicsHardeningConstant
    value = 0.61 # 35deg
  []
  [mc_dil]
    type = TensorMechanicsHardeningConstant
    value = 0.15 # 8deg
  []

  [mc_tensile_str_strong_harden]
    type = TensorMechanicsHardeningExponential
    value_0 = 1.0 # MPa
    value_residual = 1.0 # MPa
    rate = 1.0
  []
  [mc_compressive_str]
    type = TensorMechanicsHardeningCubic
    value_0 = 100 # Large!
    value_residual = 100
    internal_limit = 0.1
  []

  [wp_coh_harden]
    type = TensorMechanicsHardeningCubic
    value_0 = 0.05
    value_residual = 0.05
    internal_limit = 10
  []
  [wp_tan_fric]
    type = TensorMechanicsHardeningConstant
    value = 0.26 # 15deg
  []
  [wp_tan_dil]
    type = TensorMechanicsHardeningConstant
    value = 0.18 # 10deg
  []

  [wp_tensile_str_harden]
    type = TensorMechanicsHardeningCubic
    value_0 = 0.05
    value_residual = 0.05
    internal_limit = 10
  []
  [wp_compressive_str_soften]
    type = TensorMechanicsHardeningCubic
    value_0 = 100
    value_residual = 1
    internal_limit = 1.0
  []
[]

[Modules]
  [FluidProperties]
    [simple_fluid]
      type = SimpleFluidProperties
      bulk_modulus = 2E3
      density0 = 1000
      thermal_expansion = 0
      viscosity = 3.5E-17
    []
  []
[]

[Materials]
  [temperature]
    type = PorousFlowTemperature
  []
  [eff_fluid_pressure]
    type = PorousFlowEffectiveFluidPressure
  []
  [vol_strain]
    type = PorousFlowVolumetricStrain
  []
  [ppss]
    type = PorousFlow1PhaseP
    porepressure = porepressure
    capillary_pressure = pc
  []
  [massfrac]
    type = PorousFlowMassFraction
  []
  [simple_fluid]
    type = PorousFlowSingleComponentFluid
    fp = simple_fluid
    phase = 0
  []
  [porosity_for_aux]
    type = PorousFlowPorosity
    at_nodes = false
    fluid = true
    mechanical = true
    ensure_positive = true
    porosity_zero = 0.02
    solid_bulk = 5.3333E3
  []
  [porosity_bulk]
    type = PorousFlowPorosity
    fluid = true
    mechanical = true
    block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30'
    ensure_positive = true
    porosity_zero = 0.02
    solid_bulk = 5.3333E3
  []
  [porosity_excav]
    type = PorousFlowPorosityConst
    block = 1
    porosity = 1.0
  []
  [permeability_bulk]
    type = PorousFlowPermeabilityKozenyCarman
    block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30'
    poroperm_function = kozeny_carman_phi0
    k0 = 1E-15
    phi0 = 0.02
    n = 2
    m = 2
  []
  [permeability_excav]
    type = PorousFlowPermeabilityConst
    block = 1
    permeability = '0 0 0   0 0 0   0 0 0'
  []
  [relperm]
    type = PorousFlowRelativePermeabilityCorey
    n = 4
    s_res = 0.4
    sum_s_res = 0.4
    phase = 0
  []

  [elasticity_tensor_0]
    type = ComputeLayeredCosseratElasticityTensor
     block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30'
    young = 8E3 # MPa
    poisson = 0.25
    layer_thickness = 1.0
    joint_normal_stiffness = 1E9 # huge
    joint_shear_stiffness = 1E3 # MPa
  []
  [elasticity_tensor_1]
    type = ComputeLayeredCosseratElasticityTensor
    block = 1
    young = 8E3 # MPa
    poisson = 0.25
    layer_thickness = 1.0
    joint_normal_stiffness = 1E9 # huge
    joint_shear_stiffness = 1E3 # MPa
    elasticity_tensor_prefactor = excav_sideways
  []
  [strain]
    type = ComputeCosseratIncrementalSmallStrain
    eigenstrain_names = ini_stress
  []
  [ini_stress]
    type = ComputeEigenstrainFromInitialStress
    eigenstrain_name = ini_stress
    initial_stress = 'ini_xx 0 0  0 ini_xx 0  0 0 ini_zz'
  []

  [stress_0]
    type = ComputeMultipleInelasticCosseratStress
     block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30'
    inelastic_models = 'mc wp'
    cycle_models = true
    relative_tolerance = 2.0
    absolute_tolerance = 1E6
    max_iterations = 1
    tangent_operator = nonlinear
    perform_finite_strain_rotations = false
  []
  [stress_1]
    type = ComputeMultipleInelasticCosseratStress
    block = 1
    inelastic_models = ''
    relative_tolerance = 2.0
    absolute_tolerance = 1E6
    max_iterations = 1
    tangent_operator = nonlinear
    perform_finite_strain_rotations = false
  []
  [mc]
    type = CappedMohrCoulombCosseratStressUpdate
    warn_about_precision_loss = false
    host_youngs_modulus = 8E3
    host_poissons_ratio = 0.25
    base_name = mc
    tensile_strength = mc_tensile_str_strong_harden
    compressive_strength = mc_compressive_str
    cohesion = mc_coh_strong_harden
    friction_angle = mc_fric
    dilation_angle = mc_dil
    max_NR_iterations = 100000
    smoothing_tol = 0.1 # MPa  # Must be linked to cohesion
    yield_function_tol = 1E-9 # MPa.  this is essentially the lowest possible without lots of precision loss
    perfect_guess = true
    min_step_size = 1.0
  []
  [wp]
    type = CappedWeakPlaneCosseratStressUpdate
    warn_about_precision_loss = false
    base_name = wp
    cohesion = wp_coh_harden
    tan_friction_angle = wp_tan_fric
    tan_dilation_angle = wp_tan_dil
    tensile_strength = wp_tensile_str_harden
    compressive_strength = wp_compressive_str_soften
    max_NR_iterations = 10000
    tip_smoother = 0.05
    smoothing_tol = 0.05 # MPa  # Note, this must be tied to cohesion, otherwise get no possible return at cone apex
    yield_function_tol = 1E-11 # MPa.  this is essentially the lowest possible without lots of precision loss
    perfect_guess = true
    min_step_size = 1.0E-3
  []


  [undrained_density_0]
    type = GenericConstantMaterial
     block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30'
    prop_names = density
    prop_values = 2500
  []
  [undrained_density_1]
    type = GenericFunctionMaterial
    block = 1
    prop_names = density
    prop_values = density_sideways
  []
[]

[Preconditioning]
  [SMP]
    type = SMP
    full = true
  []
[]

[Postprocessors]
  [min_roof_disp]
    type = NodalExtremeValue
    boundary = roof
    value_type = min
    variable = disp_z
  []
  [min_roof_pp]
    type = NodalExtremeValue
    boundary = roof
    value_type = min
    variable = porepressure
  []
  [min_surface_disp]
    type = NodalExtremeValue
    boundary = zmax
    value_type = min
    variable = disp_z
  []
  [min_surface_pp]
    type = NodalExtremeValue
    boundary = zmax
    value_type = min
    variable = porepressure
  []
  [max_perm_zz]
    type = ElementExtremeValue
     block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30'
    variable = perm_zz
  []
[]

[Executioner]
  type = Transient

  solve_type = 'NEWTON'

  petsc_options = '-snes_converged_reason'

  # best overall
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
  petsc_options_value = ' lu       mumps'

  # best if you don't have mumps:
  #petsc_options_iname = '-pc_type -pc_asm_overlap -sub_pc_type -ksp_type -ksp_gmres_restart'
  #petsc_options_value = ' asm      2              lu            gmres     200'

  # very basic:
  #petsc_options_iname = '-pc_type -ksp_type -ksp_gmres_restart'
  #petsc_options_value = ' bjacobi  gmres     200'

  line_search = bt

  nl_abs_tol = 1e-3
  nl_rel_tol = 1e-5

  l_max_its = 200
  nl_max_its = 30

  start_time = 0.0
  dt = 0.0025
  end_time = 0.5
[]

[Outputs]
  interval = 1
  print_linear_residuals = true
  exodus = true
  csv = true
  console = true
[]

(modules/porous_flow/examples/coal_mining/coarse_with_fluid.i)

# Strata deformation and fluid flow aaround a coal mine - 3D model
#
# A "half model" is used.  The mine is 400m deep and
# just the roof is studied (-400<=z<=0).  The mining panel
# sits between 0<=x<=150, and 0<=y<=1000, so this simulates
# a coal panel that is 300m wide and 1000m long.  The outer boundaries
# are 1km from the excavation boundaries.
#
# The excavation takes 0.5 years.
#
# The boundary conditions for this simulation are:
#  - disp_x = 0 at x=0 and x=1150
#  - disp_y = 0 at y=-1000 and y=1000
#  - disp_z = 0 at z=-400, but there is a time-dependent
#               Young modulus that simulates excavation
#  - wc_x = 0 at y=-1000 and y=1000
#  - wc_y = 0 at x=0 and x=1150
#  - no flow at x=0, z=-400 and z=0
#  - fixed porepressure at y=-1000, y=1000 and x=1150
# That is, rollers on the sides, free at top,
# and prescribed at bottom in the unexcavated portion.
#
# A single-phase unsaturated fluid is used.
#
# The small strain formulation is used.
#
# All stresses are measured in MPa, and time units are measured in years.
#
# The initial porepressure is hydrostatic with P=0 at z=0, so
# Porepressure ~ - 0.01*z MPa, where the fluid has density 1E3 kg/m^3 and
# gravity = = 10 m.s^-2 = 1E-5 MPa m^2/kg.
# To be more accurate, i use
# Porepressure = -bulk * log(1 + g*rho0*z/bulk)
# where bulk=2E3 MPa and rho0=1Ee kg/m^3.
# The initial stress is consistent with the weight force from undrained
# density 2500 kg/m^3, and fluid porepressure, and a Biot coefficient of 0.7, ie,
# stress_zz^effective = 0.025*z + 0.7 * initial_porepressure
# The maximum and minimum principal horizontal effective stresses are
# assumed to be equal to 0.8*stress_zz.
#
# Material properties:
# Young's modulus = 8 GPa
# Poisson's ratio = 0.25
# Cosserat layer thickness = 1 m
# Cosserat-joint normal stiffness = large
# Cosserat-joint shear stiffness = 1 GPa
# MC cohesion = 2 MPa
# MC friction angle = 35 deg
# MC dilation angle = 8 deg
# MC tensile strength = 1 MPa
# MC compressive strength = 100 MPa
# WeakPlane cohesion = 0.1 MPa
# WeakPlane friction angle = 30 deg
# WeakPlane dilation angle = 10 deg
# WeakPlane tensile strength = 0.1 MPa
# WeakPlane compressive strength = 100 MPa softening to 1 MPa at strain = 1
# Fluid density at zero porepressure = 1E3 kg/m^3
# Fluid bulk modulus = 2E3 MPa
# Fluid viscosity = 1.1E-3 Pa.s = 1.1E-9 MPa.s = 3.5E-17 MPa.year
#
[GlobalParams]
  perform_finite_strain_rotations = false
  displacements = 'disp_x disp_y disp_z'
  Cosserat_rotations = 'wc_x wc_y wc_z'
  PorousFlowDictator = dictator
  biot_coefficient = 0.7
[]

[Mesh]
  [file]
    type = FileMeshGenerator
    file = mesh/coarse.e
  []
  [xmin]
    type = SideSetsAroundSubdomainGenerator
     block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16'
    new_boundary = xmin
    normal = '-1 0 0'
    input = file
  []
  [xmax]
    type = SideSetsAroundSubdomainGenerator
     block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16'
    new_boundary = xmax
    normal = '1 0 0'
    input = xmin
  []
  [ymin]
    type = SideSetsAroundSubdomainGenerator
     block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16'
    new_boundary = ymin
    normal = '0 -1 0'
    input = xmax
  []
  [ymax]
    type = SideSetsAroundSubdomainGenerator
     block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16'
    new_boundary = ymax
    normal = '0 1 0'
    input = ymin
  []
  [zmax]
    type = SideSetsAroundSubdomainGenerator
    block = 16
    new_boundary = zmax
    normal = '0 0 1'
    input = ymax
  []
  [zmin]
    type = SideSetsAroundSubdomainGenerator
    block = 2
    new_boundary = zmin
    normal = '0 0 -1'
    input = zmax
  []
  [excav]
    type = SubdomainBoundingBoxGenerator
    input = zmin
    block_id = 1
    bottom_left = '0 0 -400'
    top_right = '150 1000 -397'
  []
  [roof]
    type = SideSetsBetweenSubdomainsGenerator
    primary_block = 3
    paired_block = 1
    input = excav
    new_boundary = roof
  []
[]

[Variables]
  [disp_x]
  []
  [disp_y]
  []
  [disp_z]
  []
  [wc_x]
  []
  [wc_y]
  []
  [porepressure]
    scaling = 1E-5
  []
[]

[ICs]
  [porepressure]
    type = FunctionIC
    variable = porepressure
    function = ini_pp
  []
[]

[Kernels]
  [cx_elastic]
    type = CosseratStressDivergenceTensors
    use_displaced_mesh = false
    variable = disp_x
    component = 0
  []
  [cy_elastic]
    type = CosseratStressDivergenceTensors
    use_displaced_mesh = false
    variable = disp_y
    component = 1
  []
  [cz_elastic]
    type = CosseratStressDivergenceTensors
    use_displaced_mesh = false
    variable = disp_z
    component = 2
  []
  [x_couple]
    type = StressDivergenceTensors
    use_displaced_mesh = false
    variable = wc_x
    displacements = 'wc_x wc_y wc_z'
    component = 0
    base_name = couple
  []
  [y_couple]
    type = StressDivergenceTensors
    use_displaced_mesh = false
    variable = wc_y
    displacements = 'wc_x wc_y wc_z'
    component = 1
    base_name = couple
  []
  [x_moment]
    type = MomentBalancing
    use_displaced_mesh = false
    variable = wc_x
    component = 0
  []
  [y_moment]
    type = MomentBalancing
    use_displaced_mesh = false
    variable = wc_y
    component = 1
  []
  [gravity]
    type = Gravity
    use_displaced_mesh = false
    variable = disp_z
    value = -10E-6 # remember this is in MPa
  []
  [poro_x]
    type = PorousFlowEffectiveStressCoupling
    use_displaced_mesh = false
    variable = disp_x
    component = 0
  []
  [poro_y]
    type = PorousFlowEffectiveStressCoupling
    use_displaced_mesh = false
    variable = disp_y
    component = 1
  []
  [poro_z]
    type = PorousFlowEffectiveStressCoupling
    use_displaced_mesh = false
    component = 2
    variable = disp_z
  []
  [mass0]
    type = PorousFlowMassTimeDerivative
    fluid_component = 0
    variable = porepressure
  []
  [flux]
    type = PorousFlowAdvectiveFlux
    use_displaced_mesh = false
    variable = porepressure
    gravity = '0 0 -10E-6'
    fluid_component = 0
  []
  [poro_vol_exp]
    type = PorousFlowMassVolumetricExpansion
    block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16'
    variable = porepressure
    fluid_component = 0
  []
[]


[AuxVariables]
  [saturation]
    order = CONSTANT
    family = MONOMIAL
  []
  [darcy_x]
    order = CONSTANT
    family = MONOMIAL
  []
  [darcy_y]
    order = CONSTANT
    family = MONOMIAL
  []
  [darcy_z]
    order = CONSTANT
    family = MONOMIAL
  []
  [porosity]
    order = CONSTANT
    family = MONOMIAL
  []
  [wc_z]
  []
  [stress_xx]
    order = CONSTANT
    family = MONOMIAL
  []
  [stress_xy]
    order = CONSTANT
    family = MONOMIAL
  []
  [stress_xz]
    order = CONSTANT
    family = MONOMIAL
  []
  [stress_yx]
    order = CONSTANT
    family = MONOMIAL
  []
  [stress_yy]
    order = CONSTANT
    family = MONOMIAL
  []
  [stress_yz]
    order = CONSTANT
    family = MONOMIAL
  []
  [stress_zx]
    order = CONSTANT
    family = MONOMIAL
  []
  [stress_zy]
    order = CONSTANT
    family = MONOMIAL
  []
  [stress_zz]
    order = CONSTANT
    family = MONOMIAL
  []
  [total_strain_xx]
    order = CONSTANT
    family = MONOMIAL
  []
  [total_strain_xy]
    order = CONSTANT
    family = MONOMIAL
  []
  [total_strain_xz]
    order = CONSTANT
    family = MONOMIAL
  []
  [total_strain_yx]
    order = CONSTANT
    family = MONOMIAL
  []
  [total_strain_yy]
    order = CONSTANT
    family = MONOMIAL
  []
  [total_strain_yz]
    order = CONSTANT
    family = MONOMIAL
  []
  [total_strain_zx]
    order = CONSTANT
    family = MONOMIAL
  []
  [total_strain_zy]
    order = CONSTANT
    family = MONOMIAL
  []
  [total_strain_zz]
    order = CONSTANT
    family = MONOMIAL
  []
  [perm_xx]
    order = CONSTANT
    family = MONOMIAL
  []
  [perm_yy]
    order = CONSTANT
    family = MONOMIAL
  []
  [perm_zz]
    order = CONSTANT
    family = MONOMIAL
  []
  [mc_shear]
    order = CONSTANT
    family = MONOMIAL
  []
  [mc_tensile]
    order = CONSTANT
    family = MONOMIAL
  []
  [wp_shear]
    order = CONSTANT
    family = MONOMIAL
  []
  [wp_tensile]
    order = CONSTANT
    family = MONOMIAL
  []
  [wp_shear_f]
    order = CONSTANT
    family = MONOMIAL
  []
  [wp_tensile_f]
    order = CONSTANT
    family = MONOMIAL
  []
  [mc_shear_f]
    order = CONSTANT
    family = MONOMIAL
  []
  [mc_tensile_f]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[AuxKernels]
  [saturation_water]
    type = PorousFlowPropertyAux
    variable = saturation
    property = saturation
    phase = 0
    execute_on = timestep_end
  []
  [darcy_x]
    type = PorousFlowDarcyVelocityComponent
    variable = darcy_x
    gravity = '0 0 -10E-6'
    component = x
  []
  [darcy_y]
    type = PorousFlowDarcyVelocityComponent
    variable = darcy_y
    gravity = '0 0 -10E-6'
    component = y
  []
  [darcy_z]
    type = PorousFlowDarcyVelocityComponent
    variable = darcy_z
    gravity = '0 0 -10E-6'
    component = z
  []
  [porosity]
    type = PorousFlowPropertyAux
    property = porosity
    variable = porosity
    execute_on = timestep_end
  []
  [stress_xx]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_xx
    index_i = 0
    index_j = 0
    execute_on = timestep_end
  []
  [stress_xy]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_xy
    index_i = 0
    index_j = 1
    execute_on = timestep_end
  []
  [stress_xz]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_xz
    index_i = 0
    index_j = 2
    execute_on = timestep_end
  []
  [stress_yx]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_yx
    index_i = 1
    index_j = 0
    execute_on = timestep_end
  []
  [stress_yy]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_yy
    index_i = 1
    index_j = 1
    execute_on = timestep_end
  []
  [stress_yz]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_yz
    index_i = 1
    index_j = 2
    execute_on = timestep_end
  []
  [stress_zx]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_zx
    index_i = 2
    index_j = 0
    execute_on = timestep_end
  []
  [stress_zy]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_zy
    index_i = 2
    index_j = 1
    execute_on = timestep_end
  []
  [stress_zz]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_zz
    index_i = 2
    index_j = 2
    execute_on = timestep_end
  []
  [total_strain_xx]
    type = RankTwoAux
    rank_two_tensor = total_strain
    variable = total_strain_xx
    index_i = 0
    index_j = 0
    execute_on = timestep_end
  []
  [total_strain_xy]
    type = RankTwoAux
    rank_two_tensor = total_strain
    variable = total_strain_xy
    index_i = 0
    index_j = 1
    execute_on = timestep_end
  []
  [total_strain_xz]
    type = RankTwoAux
    rank_two_tensor = total_strain
    variable = total_strain_xz
    index_i = 0
    index_j = 2
    execute_on = timestep_end
  []
  [total_strain_yx]
    type = RankTwoAux
    rank_two_tensor = total_strain
    variable = total_strain_yx
    index_i = 1
    index_j = 0
    execute_on = timestep_end
  []
  [total_strain_yy]
    type = RankTwoAux
    rank_two_tensor = total_strain
    variable = total_strain_yy
    index_i = 1
    index_j = 1
    execute_on = timestep_end
  []
  [total_strain_yz]
    type = RankTwoAux
    rank_two_tensor = total_strain
    variable = total_strain_yz
    index_i = 1
    index_j = 2
    execute_on = timestep_end
  []
  [total_strain_zx]
    type = RankTwoAux
    rank_two_tensor = total_strain
    variable = total_strain_zx
    index_i = 2
    index_j = 0
    execute_on = timestep_end
  []
  [total_strain_zy]
    type = RankTwoAux
    rank_two_tensor = total_strain
    variable = total_strain_zy
    index_i = 2
    index_j = 1
    execute_on = timestep_end
  []
  [total_strain_zz]
    type = RankTwoAux
    rank_two_tensor = total_strain
    variable = total_strain_zz
    index_i = 2
    index_j = 2
    execute_on = timestep_end
  []
  [perm_xx]
    type = PorousFlowPropertyAux
    property = permeability
    variable = perm_xx
    row = 0
    column = 0
    execute_on = timestep_end
  []
  [perm_yy]
    type = PorousFlowPropertyAux
    property = permeability
    variable = perm_yy
    row = 1
    column = 1
    execute_on = timestep_end
  []
  [perm_zz]
    type = PorousFlowPropertyAux
    property = permeability
    variable = perm_zz
    row = 2
    column = 2
    execute_on = timestep_end
  []
  [mc_shear]
    type = MaterialStdVectorAux
    index = 0
    property = mc_plastic_internal_parameter
    variable = mc_shear
    execute_on = timestep_end
  []
  [mc_tensile]
    type = MaterialStdVectorAux
    index = 1
    property = mc_plastic_internal_parameter
    variable = mc_tensile
    execute_on = timestep_end
  []
  [wp_shear]
    type = MaterialStdVectorAux
    index = 0
    property = wp_plastic_internal_parameter
    variable = wp_shear
    execute_on = timestep_end
  []
  [wp_tensile]
    type = MaterialStdVectorAux
    index = 1
    property = wp_plastic_internal_parameter
    variable = wp_tensile
    execute_on = timestep_end
  []
  [mc_shear_f]
    type = MaterialStdVectorAux
    index = 6
    property = mc_plastic_yield_function
    variable = mc_shear_f
    execute_on = timestep_end
  []
  [mc_tensile_f]
    type = MaterialStdVectorAux
    index = 0
    property = mc_plastic_yield_function
    variable = mc_tensile_f
    execute_on = timestep_end
  []
  [wp_shear_f]
    type = MaterialStdVectorAux
    index = 0
    property = wp_plastic_yield_function
    variable = wp_shear_f
    execute_on = timestep_end
  []
  [wp_tensile_f]
    type = MaterialStdVectorAux
    index = 1
    property = wp_plastic_yield_function
    variable = wp_tensile_f
    execute_on = timestep_end
  []
[]



[BCs]
  [no_x]
    type = DirichletBC
    variable = disp_x
    boundary = 'xmin xmax'
    value = 0.0
  []
  [no_y]
    type = DirichletBC
    variable = disp_y
    boundary = 'ymin ymax'
    value = 0.0
  []
  [no_z]
    type = DirichletBC
    variable = disp_z
    boundary = zmin
    value = 0.0
  []
  [no_wc_x]
    type = DirichletBC
    variable = wc_x
    boundary = 'ymin ymax'
    value = 0.0
  []
  [no_wc_y]
    type = DirichletBC
    variable = wc_y
    boundary = 'xmin xmax'
    value = 0.0
  []
  [fix_porepressure]
    type = FunctionDirichletBC
    variable = porepressure
    boundary = 'ymin ymax xmax'
    function = ini_pp
  []
  [roof_porepressure]
    type = PorousFlowPiecewiseLinearSink
    variable = porepressure
    pt_vals = '-1E3 1E3'
    multipliers = '-1 1'
    fluid_phase = 0
    flux_function = roof_conductance
    boundary = roof
  []
  [roof_bcs]
    type = StickyBC
    variable = disp_z
    min_value = -3.0
    boundary = roof
  []
[]

[Functions]
  [ini_pp]
    type = ParsedFunction
    vars = 'bulk p0 g    rho0'
    vals = '2E3 0.0 1E-5 1E3'
    value = '-bulk*log(exp(-p0/bulk)+g*rho0*z/bulk)'
  []
  [ini_xx]
    type = ParsedFunction
    vars = 'bulk p0 g    rho0 biot'
    vals = '2E3 0.0 1E-5 1E3  0.7'
    value = '0.8*(2500*10E-6*z+biot*(-bulk*log(exp(-p0/bulk)+g*rho0*z/bulk)))'
  []
  [ini_zz]
    type = ParsedFunction
    vars = 'bulk p0 g    rho0 biot'
    vals = '2E3 0.0 1E-5 1E3  0.7'
    value = '2500*10E-6*z+biot*(-bulk*log(exp(-p0/bulk)+g*rho0*z/bulk))'
  []
  [excav_sideways]
    type = ParsedFunction
    vars = 'end_t ymin ymax  minval maxval slope'
    vals = '0.5   0    1000.0 1E-9 1 60'
    # excavation face at ymin+(ymax-ymin)*min(t/end_t,1)
    # slope is the distance over which the modulus reduces from maxval to minval
    value = 'if(y<ymin+(ymax-ymin)*min(t/end_t,1),minval,if(y<ymin+(ymax-ymin)*min(t/end_t,1)+slope,minval+(maxval-minval)*(y-(ymin+(ymax-ymin)*min(t/end_t,1)))/slope,maxval))'
  []
  [density_sideways]
    type = ParsedFunction
    vars = 'end_t ymin ymax  minval maxval'
    vals = '0.5   0    1000.0 0 2500'
    value = 'if(y<ymin+(ymax-ymin)*min(t/end_t,1),minval,maxval)'
  []
  [roof_conductance]
    type = ParsedFunction
    vars = 'end_t ymin ymax   maxval minval'
    vals = '0.5   0    1000.0 1E7      0'
    value = 'if(y<ymin+(ymax-ymin)*min(t/end_t,1),maxval,minval)'
  []
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'porepressure disp_x disp_y disp_z'
    number_fluid_phases = 1
    number_fluid_components = 1
  []
  [pc]
    type = PorousFlowCapillaryPressureVG
    m = 0.5
    alpha = 1 # MPa^-1
  []

  [mc_coh_strong_harden]
    type = TensorMechanicsHardeningExponential
    value_0 = 1.99 # MPa
    value_residual = 2.01 # MPa
    rate = 1.0
  []
  [mc_fric]
    type = TensorMechanicsHardeningConstant
    value = 0.61 # 35deg
  []
  [mc_dil]
    type = TensorMechanicsHardeningConstant
    value = 0.15 # 8deg
  []

  [mc_tensile_str_strong_harden]
    type = TensorMechanicsHardeningExponential
    value_0 = 1.0 # MPa
    value_residual = 1.0 # MPa
    rate = 1.0
  []
  [mc_compressive_str]
    type = TensorMechanicsHardeningCubic
    value_0 = 100 # Large!
    value_residual = 100
    internal_limit = 0.1
  []

  [wp_coh_harden]
    type = TensorMechanicsHardeningCubic
    value_0 = 0.05
    value_residual = 0.05
    internal_limit = 10
  []
  [wp_tan_fric]
    type = TensorMechanicsHardeningConstant
    value = 0.26 # 15deg
  []
  [wp_tan_dil]
    type = TensorMechanicsHardeningConstant
    value = 0.18 # 10deg
  []

  [wp_tensile_str_harden]
    type = TensorMechanicsHardeningCubic
    value_0 = 0.05
    value_residual = 0.05
    internal_limit = 10
  []
  [wp_compressive_str_soften]
    type = TensorMechanicsHardeningCubic
    value_0 = 100
    value_residual = 1
    internal_limit = 1.0
  []
[]

[Modules]
  [FluidProperties]
    [simple_fluid]
      type = SimpleFluidProperties
      bulk_modulus = 2E3
      density0 = 1000
      thermal_expansion = 0
      viscosity = 3.5E-17
    []
  []
[]

[Materials]
  [temperature]
    type = PorousFlowTemperature
  []
  [eff_fluid_pressure]
    type = PorousFlowEffectiveFluidPressure
  []
  [vol_strain]
    type = PorousFlowVolumetricStrain
  []
  [ppss]
    type = PorousFlow1PhaseP
    porepressure = porepressure
    capillary_pressure = pc
  []
  [massfrac]
    type = PorousFlowMassFraction
  []
  [simple_fluid]
    type = PorousFlowSingleComponentFluid
    fp = simple_fluid
    phase = 0
  []
  [porosity_bulk]
    type = PorousFlowPorosity
    fluid = true
    mechanical = true
    block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16'
    ensure_positive = true
    porosity_zero = 0.02
    solid_bulk = 5.3333E3
  []
  [porosity_excav]
    type = PorousFlowPorosityConst
    block = 1
    porosity = 1.0
  []
  [permeability_bulk]
    type = PorousFlowPermeabilityKozenyCarman
    block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16'
    poroperm_function = kozeny_carman_phi0
    k0 = 1E-15
    phi0 = 0.02
    n = 2
    m = 2
  []
  [permeability_excav]
    type = PorousFlowPermeabilityConst
    block = 1
    permeability = '0 0 0   0 0 0   0 0 0'
  []
  [relperm]
    type = PorousFlowRelativePermeabilityCorey
    n = 4
    s_res = 0.4
    sum_s_res = 0.4
    phase = 0
  []

  [elasticity_tensor_0]
    type = ComputeLayeredCosseratElasticityTensor
     block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16'
    young = 8E3 # MPa
    poisson = 0.25
    layer_thickness = 1.0
    joint_normal_stiffness = 1E9 # huge
    joint_shear_stiffness = 1E3 # MPa
  []
  [elasticity_tensor_1]
    type = ComputeLayeredCosseratElasticityTensor
    block = 1
    young = 8E3 # MPa
    poisson = 0.25
    layer_thickness = 1.0
    joint_normal_stiffness = 1E9 # huge
    joint_shear_stiffness = 1E3 # MPa
    elasticity_tensor_prefactor = excav_sideways
  []
  [strain]
    type = ComputeCosseratIncrementalSmallStrain
    eigenstrain_names = ini_stress
  []
  [ini_stress]
    type = ComputeEigenstrainFromInitialStress
    eigenstrain_name = ini_stress
    initial_stress = 'ini_xx 0 0  0 ini_xx 0  0 0 ini_zz'
  []

  [stress_0]
    type = ComputeMultipleInelasticCosseratStress
     block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16'
    inelastic_models = 'mc wp'
    cycle_models = true
    relative_tolerance = 2.0
    absolute_tolerance = 1E6
    max_iterations = 1
    tangent_operator = nonlinear
    perform_finite_strain_rotations = false
  []
  [stress_1]
    type = ComputeMultipleInelasticCosseratStress
    block = 1
    inelastic_models = ''
    relative_tolerance = 2.0
    absolute_tolerance = 1E6
    max_iterations = 1
    tangent_operator = nonlinear
    perform_finite_strain_rotations = false
  []
  [mc]
    type = CappedMohrCoulombCosseratStressUpdate
    warn_about_precision_loss = false
    host_youngs_modulus = 8E3
    host_poissons_ratio = 0.25
    base_name = mc
    tensile_strength = mc_tensile_str_strong_harden
    compressive_strength = mc_compressive_str
    cohesion = mc_coh_strong_harden
    friction_angle = mc_fric
    dilation_angle = mc_dil
    max_NR_iterations = 100000
    smoothing_tol = 0.1 # MPa  # Must be linked to cohesion
    yield_function_tol = 1E-9 # MPa.  this is essentially the lowest possible without lots of precision loss
    perfect_guess = true
    min_step_size = 1.0
  []
  [wp]
    type = CappedWeakPlaneCosseratStressUpdate
    warn_about_precision_loss = false
    base_name = wp
    cohesion = wp_coh_harden
    tan_friction_angle = wp_tan_fric
    tan_dilation_angle = wp_tan_dil
    tensile_strength = wp_tensile_str_harden
    compressive_strength = wp_compressive_str_soften
    max_NR_iterations = 10000
    tip_smoother = 0.05
    smoothing_tol = 0.05 # MPa  # Note, this must be tied to cohesion, otherwise get no possible return at cone apex
    yield_function_tol = 1E-11 # MPa.  this is essentially the lowest possible without lots of precision loss
    perfect_guess = true
    min_step_size = 1.0E-3
  []


  [undrained_density_0]
    type = GenericConstantMaterial
     block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16'
    prop_names = density
    prop_values = 2500
  []
  [undrained_density_1]
    type = GenericFunctionMaterial
    block = 1
    prop_names = density
    prop_values = density_sideways
  []
[]

[Preconditioning]
  [SMP]
    type = SMP
    full = true
  []
[]

[Postprocessors]
  [min_roof_disp]
    type = NodalExtremeValue
    boundary = roof
    value_type = min
    variable = disp_z
  []
  [min_roof_pp]
    type = NodalExtremeValue
    boundary = roof
    value_type = min
    variable = porepressure
  []
  [min_surface_disp]
    type = NodalExtremeValue
    boundary = zmax
    value_type = min
    variable = disp_z
  []
  [min_surface_pp]
    type = NodalExtremeValue
    boundary = zmax
    value_type = min
    variable = porepressure
  []
  [max_perm_zz]
    type = ElementExtremeValue
     block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16'
    variable = perm_zz
  []
[]

[Executioner]
  type = Transient

  solve_type = 'NEWTON'

  petsc_options = '-snes_converged_reason'

  # best overall
  # petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
  # petsc_options_value = ' lu       mumps'

  # best if you do not have mumps:
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
  petsc_options_value = ' lu       superlu_dist'

  # best if you do not have mumps or superlu_dist:
  #petsc_options_iname = '-pc_type -pc_asm_overlap -sub_pc_type -ksp_type -ksp_gmres_restart'
  #petsc_options_value = ' asm      2              lu            gmres     200'

  # very basic:
  #petsc_options_iname = '-pc_type -ksp_type -ksp_gmres_restart'
  #petsc_options_value = ' bjacobi  gmres     200'

  line_search = bt

  nl_abs_tol = 1e-3
  nl_rel_tol = 1e-5

  l_max_its = 200
  nl_max_its = 30

  start_time = 0.0
  dt = 0.014706
  end_time = 0.014706 #0.5
[]

[Outputs]
  interval = 1
  print_linear_residuals = true
  exodus = true
  csv = true
  console = true
[]

(modules/tensor_mechanics/examples/coal_mining/coarse.i)

# Strata deformation and fracturing around a coal mine - 3D model
#
# A "half model" is used.  The mine is 400m deep and
# just the roof is studied (-400<=z<=0).  The mining panel
# sits between 0<=x<=150, and 0<=y<=1000, so this simulates
# a coal panel that is 300m wide and 1000m long.  The outer boundaries
# are 1km from the excavation boundaries.
#
# Time is meaningless in this example
# as quasi-static solutions are sought at each timestep, but
# the number of timesteps controls the resolution of the
# process.
#
# The boundary conditions for this simulation are:
#  - disp_x = 0 at x=0 and x=1150
#  - disp_y = 0 at y=-1000 and y=1000
#  - disp_z = 0 at z=-400, but there is a time-dependent
#               Young's modulus that simulates excavation
#  - wc_x = 0 at y=-1000 and y=1000
#  - wc_y = 0 at x=0 and x=1150
# That is, rollers on the sides, free at top,
# and prescribed at bottom in the unexcavated portion.
#
# The small strain formulation is used.
#
# All stresses are measured in MPa.  The initial stress is consistent with
# the weight force from density 2500 kg/m^3, ie, stress_zz = 0.025*z MPa
# where gravity = 10 m.s^-2 = 1E-5 MPa m^2/kg.  The maximum and minimum
# principal horizontal stresses are assumed to be equal to 0.8*stress_zz.
#
# Material properties:
# Young's modulus = 8 GPa
# Poisson's ratio = 0.25
# Cosserat layer thickness = 1 m
# Cosserat-joint normal stiffness = large
# Cosserat-joint shear stiffness = 1 GPa
# MC cohesion = 3 MPa
# MC friction angle = 37 deg
# MC dilation angle = 8 deg
# MC tensile strength = 1 MPa
# MC compressive strength = 100 MPa
# WeakPlane cohesion = 0.1 MPa
# WeakPlane friction angle = 30 deg
# WeakPlane dilation angle = 10 deg
# WeakPlane tensile strength = 0.1 MPa
# WeakPlane compressive strength = 100 MPa softening to 1 MPa at strain = 1
#
[Mesh]
  [file]
    type = FileMeshGenerator
    file = mesh/coarse.e
  []
  [./xmin]
    input = file
    type = SideSetsAroundSubdomainGenerator
    block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16'
    new_boundary = xmin
    normal = '-1 0 0'
  [../]
  [./xmax]
    input = xmin
    type = SideSetsAroundSubdomainGenerator
    block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16'
    new_boundary = xmax
    normal = '1 0 0'
  [../]
  [./ymin]
    input = xmax
    type = SideSetsAroundSubdomainGenerator
    block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16'
    new_boundary = ymin
    normal = '0 -1 0'
  [../]
  [./ymax]
    input = ymin
    type = SideSetsAroundSubdomainGenerator
    block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16'
    new_boundary = ymax
    normal = '0 1 0'
  [../]
  [./zmax]
    input = ymax
    type = SideSetsAroundSubdomainGenerator
    block = 16
    new_boundary = zmax
    normal = '0 0 1'
  [../]
  [./zmin]
    input = zmax
    type = SideSetsAroundSubdomainGenerator
    block = 2
    new_boundary = zmin
    normal = '0 0 -1'
  [../]
  [./excav]
    type = SubdomainBoundingBoxGenerator
    input = zmin
    block_id = 1
    bottom_left = '0 0 -400'
    top_right = '150 1000 -397'
  [../]
  [./roof]
    type = SideSetsAroundSubdomainGenerator
    block = 1
    input = excav
    new_boundary = roof
    normal = '0 0 1'
  [../]
[]

[GlobalParams]
  perform_finite_strain_rotations = false
  displacements = 'disp_x disp_y disp_z'
  Cosserat_rotations = 'wc_x wc_y wc_z'
[]

[Variables]
  [./disp_x]
  [../]
  [./disp_y]
  [../]
  [./disp_z]
  [../]
  [./wc_x]
  [../]
  [./wc_y]
  [../]
[]

[Kernels]
  [./cx_elastic]
    type = CosseratStressDivergenceTensors
    use_displaced_mesh = false
    variable = disp_x
    component = 0
  [../]
  [./cy_elastic]
    type = CosseratStressDivergenceTensors
    use_displaced_mesh = false
    variable = disp_y
    component = 1
  [../]
  [./cz_elastic]
    type = CosseratStressDivergenceTensors
    use_displaced_mesh = false
    variable = disp_z
    component = 2
  [../]
  [./x_couple]
    type = StressDivergenceTensors
    use_displaced_mesh = false
    variable = wc_x
    displacements = 'wc_x wc_y wc_z'
    component = 0
    base_name = couple
  [../]
  [./y_couple]
    type = StressDivergenceTensors
    use_displaced_mesh = false
    variable = wc_y
    displacements = 'wc_x wc_y wc_z'
    component = 1
    base_name = couple
  [../]
  [./x_moment]
    type = MomentBalancing
    use_displaced_mesh = false
    variable = wc_x
    component = 0
  [../]
  [./y_moment]
    type = MomentBalancing
    use_displaced_mesh = false
    variable = wc_y
    component = 1
  [../]
  [./gravity]
    type = Gravity
    use_displaced_mesh = false
    variable = disp_z
    value = -10E-6 # remember this is in MPa
  [../]
[]


[AuxVariables]
  [./wc_z]
  [../]
  [./stress_xx]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./stress_xy]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./stress_xz]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./stress_yx]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./stress_yy]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./stress_yz]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./stress_zx]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./stress_zy]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./stress_zz]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./mc_shear]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./mc_tensile]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./wp_shear]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./wp_tensile]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./wp_shear_f]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./wp_tensile_f]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./mc_shear_f]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./mc_tensile_f]
    order = CONSTANT
    family = MONOMIAL
  [../]
[]

[AuxKernels]
  [./stress_xx]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_xx
    index_i = 0
    index_j = 0
  [../]
  [./stress_xy]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_xy
    index_i = 0
    index_j = 1
  [../]
  [./stress_xz]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_xz
    index_i = 0
    index_j = 2
  [../]
  [./stress_yx]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_yx
    index_i = 1
    index_j = 0
  [../]
  [./stress_yy]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_yy
    index_i = 1
    index_j = 1
  [../]
  [./stress_yz]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_yz
    index_i = 1
    index_j = 2
  [../]
  [./stress_zx]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_zx
    index_i = 2
    index_j = 0
  [../]
  [./stress_zy]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_zy
    index_i = 2
    index_j = 1
  [../]
  [./stress_zz]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_zz
    index_i = 2
    index_j = 2
  [../]
  [./mc_shear]
    type = MaterialStdVectorAux
    index = 0
    property = mc_plastic_internal_parameter
    variable = mc_shear
  [../]
  [./mc_tensile]
    type = MaterialStdVectorAux
    index = 1
    property = mc_plastic_internal_parameter
    variable = mc_tensile
  [../]
  [./wp_shear]
    type = MaterialStdVectorAux
    index = 0
    property = wp_plastic_internal_parameter
    variable = wp_shear
  [../]
  [./wp_tensile]
    type = MaterialStdVectorAux
    index = 1
    property = wp_plastic_internal_parameter
    variable = wp_tensile
  [../]
  [./mc_shear_f]
    type = MaterialStdVectorAux
    index = 6
    property = mc_plastic_yield_function
    variable = mc_shear_f
  [../]
  [./mc_tensile_f]
    type = MaterialStdVectorAux
    index = 0
    property = mc_plastic_yield_function
    variable = mc_tensile_f
  [../]
  [./wp_shear_f]
    type = MaterialStdVectorAux
    index = 0
    property = wp_plastic_yield_function
    variable = wp_shear_f
  [../]
  [./wp_tensile_f]
    type = MaterialStdVectorAux
    index = 1
    property = wp_plastic_yield_function
    variable = wp_tensile_f
  [../]
[]



[BCs]
  [./no_x]
    type = DirichletBC
    variable = disp_x
    boundary = 'xmin xmax'
    value = 0.0
  [../]
  [./no_y]
    type = DirichletBC
    variable = disp_y
    boundary = 'ymin ymax'
    value = 0.0
  [../]
  [./no_z]
    type = DirichletBC
    variable = disp_z
    boundary = zmin
    value = 0.0
  [../]
  [./no_wc_x]
    type = DirichletBC
    variable = wc_x
    boundary = 'ymin ymax'
    value = 0.0
  [../]
  [./no_wc_y]
    type = DirichletBC
    variable = wc_y
    boundary = 'xmin xmax'
    value = 0.0
  [../]
  [./roof]
    type = StickyBC
    variable = disp_z
    min_value = -3.0
    boundary = roof
  [../]
[]

[Functions]
  [./ini_xx]
    type = ParsedFunction
    value = '0.8*2500*10E-6*z'
  [../]
  [./ini_zz]
    type = ParsedFunction
    value = '2500*10E-6*z'
  [../]
  [./excav_sideways]
    type = ParsedFunction
    vars = 'end_t ymin ymax  minval maxval slope'
    vals = '17.0   0    1000.0 1E-9 1 60'
    # excavation face at ymin+(ymax-ymin)*min(t/end_t,1)
    # slope is the distance over which the modulus reduces from maxval to minval
    value = 'if(y<ymin+(ymax-ymin)*min(t/end_t,1),minval,if(y<ymin+(ymax-ymin)*min(t/end_t,1)+slope,minval+(maxval-minval)*(y-(ymin+(ymax-ymin)*min(t/end_t,1)))/slope,maxval))'
  [../]
  [./density_sideways]
    type = ParsedFunction
    vars = 'end_t ymin ymax  minval maxval'
    vals = '17.0   0    1000.0 0 2500'
    value = 'if(y<ymin+(ymax-ymin)*min(t/end_t,1),minval,maxval)'
  [../]
[]

[UserObjects]
  [./mc_coh_strong_harden]
    type = TensorMechanicsHardeningExponential
    value_0 = 2.99 # MPa
    value_residual = 3.01 # MPa
    rate = 1.0
  [../]
  [./mc_fric]
    type = TensorMechanicsHardeningConstant
    value = 0.65 # 37deg
  [../]
  [./mc_dil]
    type = TensorMechanicsHardeningConstant
    value = 0.15 # 8deg
  [../]

  [./mc_tensile_str_strong_harden]
    type = TensorMechanicsHardeningExponential
    value_0 = 1.0 # MPa
    value_residual = 1.0 # MPa
    rate = 1.0
  [../]
  [./mc_compressive_str]
    type = TensorMechanicsHardeningCubic
    value_0 = 100 # Large!
    value_residual = 100
    internal_limit = 0.1
  [../]

  [./wp_coh_harden]
    type = TensorMechanicsHardeningCubic
    value_0 = 0.1
    value_residual = 0.1
    internal_limit = 10
  [../]
  [./wp_tan_fric]
    type = TensorMechanicsHardeningConstant
    value = 0.36 # 20deg
  [../]
  [./wp_tan_dil]
    type = TensorMechanicsHardeningConstant
    value = 0.18 # 10deg
  [../]

  [./wp_tensile_str_harden]
    type = TensorMechanicsHardeningCubic
    value_0 = 0.1
    value_residual = 0.1
    internal_limit = 10
  [../]
  [./wp_compressive_str_soften]
    type = TensorMechanicsHardeningCubic
    value_0 = 100
    value_residual = 1
    internal_limit = 1.0
  [../]
[]

[Materials]
  [./elasticity_tensor_0]
    type = ComputeLayeredCosseratElasticityTensor
    block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16'
    young = 8E3 # MPa
    poisson = 0.25
    layer_thickness = 1.0
    joint_normal_stiffness = 1E9 # huge
    joint_shear_stiffness = 1E3 # MPa
  [../]
  [./elasticity_tensor_1]
    type = ComputeLayeredCosseratElasticityTensor
    block = 1
    young = 8E3 # MPa
    poisson = 0.25
    layer_thickness = 1.0
    joint_normal_stiffness = 1E9 # huge
    joint_shear_stiffness = 1E3 # MPa
    elasticity_tensor_prefactor = excav_sideways
  [../]
  [./strain]
    type = ComputeCosseratIncrementalSmallStrain
    eigenstrain_names = ini_stress
  [../]
  [./ini_stress]
    type = ComputeEigenstrainFromInitialStress
    eigenstrain_name = ini_stress
    initial_stress = 'ini_xx 0 0  0 ini_xx 0  0 0 ini_zz'
  [../]

  [./stress_0]
    type = ComputeMultipleInelasticCosseratStress
    block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16'
    inelastic_models = 'mc wp'
    cycle_models = true
    relative_tolerance = 2.0
    absolute_tolerance = 1E6
    max_iterations = 1
    tangent_operator = nonlinear
    perform_finite_strain_rotations = false
  [../]
  [./stress_1]
    type = ComputeMultipleInelasticCosseratStress
    block = 1
    inelastic_models = ''
    relative_tolerance = 2.0
    absolute_tolerance = 1E6
    max_iterations = 1
    tangent_operator = nonlinear
    perform_finite_strain_rotations = false
  [../]
  [./mc]
    type = CappedMohrCoulombCosseratStressUpdate
    warn_about_precision_loss = false
    host_youngs_modulus = 8E3
    host_poissons_ratio = 0.25
    base_name = mc
    tensile_strength = mc_tensile_str_strong_harden
    compressive_strength = mc_compressive_str
    cohesion = mc_coh_strong_harden
    friction_angle = mc_fric
    dilation_angle = mc_dil
    max_NR_iterations = 100000
    smoothing_tol = 0.1 # MPa  # Must be linked to cohesion
    yield_function_tol = 1E-9 # MPa.  this is essentially the lowest possible without lots of precision loss
    perfect_guess = true
    min_step_size = 1.0
  [../]
  [./wp]
    type = CappedWeakPlaneCosseratStressUpdate
    warn_about_precision_loss = false
    base_name = wp
    cohesion = wp_coh_harden
    tan_friction_angle = wp_tan_fric
    tan_dilation_angle = wp_tan_dil
    tensile_strength = wp_tensile_str_harden
    compressive_strength = wp_compressive_str_soften
    max_NR_iterations = 10000
    tip_smoother = 0.1
    smoothing_tol = 0.1 # MPa  # Note, this must be tied to cohesion, otherwise get no possible return at cone apex
    yield_function_tol = 1E-11 # MPa.  this is essentially the lowest possible without lots of precision loss
    perfect_guess = true
    min_step_size = 1.0E-3
  [../]


  [./density_0]
    type = GenericConstantMaterial
    block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16'
    prop_names = density
    prop_values = 2500
  [../]
  [./density_1]
    type = GenericFunctionMaterial
    block = 1
    prop_names = density
    prop_values = density_sideways
  [../]
[]

[Preconditioning]
  [./SMP]
    type = SMP
    full = true
  []
[]

[Postprocessors]
  [./min_roof_disp]
    type = NodalExtremeValue
    boundary = roof
    value_type = min
    variable = disp_z
  [../]
  [./min_surface_disp]
    type = NodalExtremeValue
    boundary = zmax
    value_type = min
    variable = disp_z
  [../]
[]

[Executioner]
  type = Transient

  solve_type = 'NEWTON'

  petsc_options = '-snes_converged_reason'
  petsc_options_iname = '-pc_type -ksp_type -ksp_gmres_restart'
  petsc_options_value = ' bjacobi  gmres     200'

  line_search = bt

  nl_abs_tol = 1e-3
  nl_rel_tol = 1e-5

  l_max_its = 30
  nl_max_its = 1000

  start_time = 0.0
  dt = 0.5 # this gives min(disp_z)=-4.3, use dt=0.0625 if you want to restrict disp_z>=-3.2
  end_time = 17.0

[]

[Outputs]
  interval = 1
  print_linear_residuals = false
  exodus = true
  csv = true
  console = true
[]

(modules/stochastic_tools/test/tests/surrogates/load_store/sub.i)

[Mesh]
  type = GeneratedMesh
  dim = 1
  nx = 100
  xmax = 10
[]

[Variables]
  [u]
  []
[]

[Kernels]
  [diffusion]
    type = MatDiffusion
    variable = u
    diffusivity = D
  []
  [absorption]
    type = MaterialReaction
    variable = u
    coefficient = sig
  []
  [source]
    type = BodyForce
    variable = u
    value = 1.0
  []
[]

[Materials]
  [diffusivity]
    type = GenericConstantMaterial
    prop_names = D
    prop_values = 2.0
  []
  [xs]
    type = GenericConstantMaterial
    prop_names = sig
    prop_values = 2.0
  []
[]

[BCs]
  [left]
    type = DirichletBC
    variable = u
    boundary = left
    value = 0
  []
  [right]
    type = DirichletBC
    variable = u
    boundary = right
    value = 0
  []
[]

[Executioner]
  type = Steady
  solve_type = PJFNK
  petsc_options_iname = '-pc_type -pc_hypre_type'
  petsc_options_value = 'hypre boomeramg'
[]

[Controls]
  [stochastic]
    type = SamplerReceiver
  []
[]

[Postprocessors]
  [avg]
    type = AverageNodalVariableValue
    variable = u
  []
  [max]
    type = NodalExtremeValue
    variable = u
    value_type = max
  []
[]

[Outputs]
[]

(modules/tensor_mechanics/test/tests/ad_thermal_expansion_function/finite_linear.i)

# This tests the thermal expansion coefficient function using both
# options to specify that function: mean and instantaneous.  There
# two blocks, each containing a single element, and these use the
# two variants of the function.

# In this test, the instantaneous CTE function is a linear function
# while the mean CTE function is an analytic function designed to
# give the same response.  If \bar{alpha}(T) is the mean CTE function,
# and \alpha(T) is the instantaneous CTE function,

# \bar{\alpha}(T) = 1/(T-Tref) \intA^{T}_{Tsf} \alpha(T) dT

# where Tref is the reference temperature used to define the mean CTE
# function, and Tsf is the stress-free temperature.

# This version of the test uses finite deformation theory.
# The two models produce very similar results.  There are slight
# differences due to the large deformation treatment.

[Mesh]
  file = 'blocks.e'
[]

[GlobalParams]
  displacements = 'disp_x disp_y disp_z'
[]

[AuxVariables]
  [./temp]
    order = FIRST
    family = LAGRANGE
  [../]
[]

[Modules/TensorMechanics/Master]
  [./all]
    strain = FINITE
    add_variables = true
    eigenstrain_names = eigenstrain
    generate_output = 'strain_xx strain_yy strain_zz'
    use_automatic_differentiation = true
  [../]
[]

[BCs]
  [./left]
    type = DirichletBC
    variable = disp_x
    boundary = 3
    value = 0.0
  [../]

  [./bottom]
    type = DirichletBC
    variable = disp_y
    boundary = 2
    value = 0.0
  [../]

  [./back]
    type = DirichletBC
    variable = disp_z
    boundary = 1
    value = 0.0
  [../]
[]

[AuxKernels]
  [./temp]
    type = FunctionAux
    variable = temp
    block = '1 2'
    function = temp_func
  [../]
[]

[Materials]
  [./elasticity_tensor]
    type = ADComputeIsotropicElasticityTensor
    youngs_modulus = 1e6
    poissons_ratio = 0.3
  [../]
  [./small_stress]
    type = ADComputeFiniteStrainElasticStress
  [../]
  [./thermal_expansion_strain1]
    type = ADComputeMeanThermalExpansionFunctionEigenstrain
    block = 1
    thermal_expansion_function = cte_func_mean
    thermal_expansion_function_reference_temperature = 0.5
    stress_free_temperature = 0.0
    temperature = temp
    eigenstrain_name = eigenstrain
  [../]
  [./thermal_expansion_strain2]
    type = ADComputeInstantaneousThermalExpansionFunctionEigenstrain
    block = 2
    thermal_expansion_function = cte_func_inst
    stress_free_temperature = 0.0
    temperature = temp
    eigenstrain_name = eigenstrain
  [../]
[]

[Functions]
  [./cte_func_mean]
    type = ParsedFunction
    vars = 'tsf tref scale' #stress free temp, reference temp, scale factor
    vals = '0.0 0.5  1e-4'
    value = 'scale * (0.5 * t^2 - 0.5 * tsf^2) / (t - tref)'
  [../]
  [./cte_func_inst]
    type = PiecewiseLinear
    xy_data = '0 0.0
               2 2.0'
    scale_factor = 1e-4
  [../]

  [./temp_func]
    type = PiecewiseLinear
    xy_data = '0 1
               1 2'
  [../]
[]

[Postprocessors]
  [./disp_1]
    type = NodalExtremeValue
    variable = disp_x
    boundary = 101
  [../]

  [./disp_2]
    type = NodalExtremeValue
    variable = disp_x
    boundary = 102
  [../]
[]

[Executioner]
  type = Transient

  solve_type = NEWTON
  l_max_its = 100
  l_tol = 1e-4
  nl_abs_tol = 1e-8
  nl_rel_tol = 1e-12

  start_time = 0.0
  end_time = 1.0
  dt = 0.1
[]

[Outputs]
  csv = true
[]

(modules/tensor_mechanics/test/tests/dynamics/rayleigh_damping/rayleigh_hht_ti.i)

# Test for rayleigh damping implemented using HHT time integration
#
# The test is for an 1D bar element of  unit length fixed on one end
# with a ramped pressure boundary condition applied to the other end.
# zeta and eta correspond to the stiffness and mass proportional rayleigh damping
# alpha, beta and gamma are HHT time integration parameters
# The equation of motion in terms of matrices is:
#
# M*accel + (eta*M+zeta*K)*[(1+alpha)vel-alpha vel_old]
# + alpha*(K*disp - K*disp_old) + K*disp = P(t+alpha dt)*Area
#
# Here M is the mass matrix, K is the stiffness matrix, P is the applied pressure
#
# This equation is equivalent to:
#
# density*accel + eta*density*[(1+alpha)vel-alpha vel_old]
# + zeta*[(1+alpha)*d/dt(Div stress)- alpha*d/dt(Div stress_old)]
# + alpha *(Div stress - Div stress_old) +Div Stress= P(t+alpha dt)
#
# The first two terms on the left are evaluated using the Inertial force kernel
# The next three terms on the left involving zeta and alpha are evaluated using
# the DynamicStressDivergenceTensors Kernel
# The residual due to Pressure is evaluated using Pressure boundary condition
#
# The system will come to steady state slowly after the pressure becomes constant.
# Alpha equal to zero will result in Newmark integration.

[Mesh]
  type = GeneratedMesh
  dim = 3
  nx = 1
  ny = 1
  nz = 1
  xmin = 0.0
  xmax = 0.1
  ymin = 0.0
  ymax = 1.0
  zmin = 0.0
  zmax = 0.1
[]

[Variables]
  [disp_x]
  []
  [disp_y]
  []
  [disp_z]
  []
[]

[AuxVariables]
  [vel_x]
  []
  [accel_x]
  []
  [vel_y]
  []
  [accel_y]
  []
  [vel_z]
  []
  [accel_z]
  []
  [stress_yy]
    order = CONSTANT
    family = MONOMIAL
  []
  [strain_yy]
    order = CONSTANT
    family = MONOMIAL
  []

[]

[Kernels]
  [DynamicTensorMechanics]
    displacements = 'disp_x disp_y disp_z'
    stiffness_damping_coefficient = 0.1
    hht_alpha = 0.11
  []
  [inertia_x]
    type = InertialForce
    variable = disp_x
    eta = 0.1
    alpha = 0.11
  []
  [inertia_y]
    type = InertialForce
    variable = disp_y
    eta = 0.1
    alpha = 0.11
  []
  [inertia_z]
    type = InertialForce
    variable = disp_z
    eta = 0.1
    alpha = 0.11
  []

[]

[AuxKernels]
  [accel_x] # These auxkernels are only to check output
    type = TestNewmarkTI
    displacement = disp_x
    variable = accel_x
    first = false
  []
  [accel_y]
    type = TestNewmarkTI
    displacement = disp_y
    variable = accel_y
    first = false
  []
  [accel_z]
    type = TestNewmarkTI
    displacement = disp_z
    variable = accel_z
    first = false
  []
  [vel_x]
    type = TestNewmarkTI
    displacement = disp_x
    variable = vel_x
  []
  [vel_y]
    type = TestNewmarkTI
    displacement = disp_y
    variable = vel_y
  []
  [vel_z]
    type = TestNewmarkTI
    displacement = disp_z
    variable = vel_z
  []
  [stress_yy]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_yy
    index_i = 1
    index_j = 1
  []
  [strain_yy]
    type = RankTwoAux
    rank_two_tensor = total_strain
    variable = strain_yy
    index_i = 1
    index_j = 1
  []

[]

[BCs]
  [top_y]
    type = DirichletBC
    variable = disp_y
    boundary = top
    value = 0.0
  []
  [top_x]
    type = DirichletBC
    variable = disp_x
    boundary = top
    value = 0.0
  []
  [top_z]
    type = DirichletBC
    variable = disp_z
    boundary = top
    value = 0.0
  []
  [bottom_x]
    type = DirichletBC
    variable = disp_x
    boundary = bottom
    value = 0.0
  []
  [bottom_z]
    type = DirichletBC
    variable = disp_z
    boundary = bottom
    value = 0.0
  []
  [Pressure]
    [Side1]
      boundary = bottom
      function = pressure
      displacements = 'disp_x disp_y disp_z'
      factor = 1
      hht_alpha = 0.11
    []
  []
[]

[Materials]
  [Elasticity_tensor]
    type = ComputeElasticityTensor
    block = 0
    fill_method = symmetric_isotropic
    C_ijkl = '210e9 0'
  []

  [strain]
    type = ComputeSmallStrain
    block = 0
    displacements = 'disp_x disp_y disp_z'
  []

  [stress]
    type = ComputeLinearElasticStress
    block = 0
  []

  [density]
    type = GenericConstantMaterial
    block = 0
    prop_names = 'density'
    prop_values = '7750'
  []

[]

[Executioner]
  type = Transient
  start_time = 0
  end_time = 2
  dt = 0.1

  # Time integrator scheme
  scheme = "newmark-beta"
[]

[Functions]
  [pressure]
    type = PiecewiseLinear
    x = '0.0 0.1 0.2 1.0 2.0 5.0'
    y = '0.0 0.1 0.2 1.0 1.0 1.0'
    scale_factor = 1e9
  []
[]

[Postprocessors]
  [_dt]
    type = TimestepSize
  []
  [disp]
    type = NodalExtremeValue
    variable = disp_y
    boundary = bottom
  []
  [vel]
    type = NodalExtremeValue
    variable = vel_y
    boundary = bottom
  []
  [accel]
    type = NodalExtremeValue
    variable = accel_y
    boundary = bottom
  []
  [stress_yy]
    type = ElementAverageValue
    variable = stress_yy
  []
  [strain_yy]
    type = ElementAverageValue
    variable = strain_yy
  []
[]

[Outputs]
  file_base = 'rayleigh_hht_out'
  exodus = true
  perf_graph = true
[]

(modules/tensor_mechanics/test/tests/volumetric_eigenstrain/volumetric_mechanical.i)

# This test ensures that the reported volumetric strain for a cube with
# mechanically imposed displacements (through Dirichlet BCs) exactly
# matches that from a version of this test that experiences the same
# defomation, but due to imposed eigenstrains.

[Mesh]
  type = GeneratedMesh
  dim = 3
  nx = 1
  ny = 1
  nz = 1
[]

[Variables]
  [./disp_x]
  [../]
  [./disp_y]
  [../]
  [./disp_z]
  [../]
[]

[AuxVariables]
  [./volumetric_strain]
    order = CONSTANT
    family = MONOMIAL
  [../]
[]

[GlobalParams]
  displacements = 'disp_x disp_y disp_z'
[]

[Modules/TensorMechanics/Master]
  [./master]
    strain = FINITE
    decomposition_method = EigenSolution #Necessary for exact solution
  [../]
[]

[AuxKernels]
  [./volumetric_strain]
    type = RankTwoScalarAux
    scalar_type = VolumetricStrain
    rank_two_tensor = total_strain
    variable = volumetric_strain
  [../]
[]

[Functions]
  [pres_disp]
    type = PiecewiseLinear
    # These values are taken from the displacements in the eigenstrain
    # version of this test. The volume of the cube (which starts out as
    # a 1x1x1 cube) is (1 + disp)^3. At time 2, this is
    # (1.44224957030741)^3, which is 3.0.
    xy_data = '0 0
               1 0.25992104989487
               2 0.44224957030741'
  []
[]

[BCs]
  [./left]
    type = DirichletBC
    variable = disp_x
    boundary = left
    value = 0.0
  [../]
  [./bottom]
    type = DirichletBC
    variable = disp_y
    boundary = bottom
    value = 0.0
  [../]
  [./back]
    type = DirichletBC
    variable = disp_z
    boundary = back
    value = 0.0
  [../]
  [./right]
    type = FunctionDirichletBC
    variable = disp_x
    boundary = right
    function = pres_disp
  [../]
  [./top]
    type = FunctionDirichletBC
    variable = disp_y
    boundary = top
    function = pres_disp
  [../]
  [./front]
    type = FunctionDirichletBC
    variable = disp_z
    boundary = front
    function = pres_disp
  [../]
[]

[Materials]
  [./elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 1e6
    poissons_ratio = 0.3
  [../]
  [./finite_strain_stress]
    type = ComputeFiniteStrainElasticStress
  [../]
  [./volumetric_change]
    type = GenericFunctionMaterial
    prop_names = volumetric_change
    prop_values = t
  [../]
[]

[Postprocessors]
  [./vol]
    type = VolumePostprocessor
    use_displaced_mesh = true
    execute_on = 'initial timestep_end'
  [../]
  [./volumetric_strain]
    type = ElementalVariableValue
    variable = volumetric_strain
    elementid = 0
  [../]
  [./disp_right]
    type = NodalExtremeValue
    variable = disp_x
    boundary = right
  [../]
[]

[Executioner]
  type = Transient
  solve_type = PJFNK
  l_max_its = 100
  l_tol = 1e-4
  nl_abs_tol = 1e-8
  nl_rel_tol = 1e-12

  start_time = 0.0
  end_time = 2.0
  dt = 1.0
[]

[Outputs]
  exodus = true
  csv = true
[]

(modules/tensor_mechanics/test/tests/plane_stress/weak_plane_stress_finite.i)

[GlobalParams]
  order = FIRST
  family = LAGRANGE
  displacements = 'disp_x disp_y'
  temperature = temp
  out_of_plane_strain = strain_zz
[]

[Mesh]
  [./square]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 2
    ny = 2
  [../]
[]

[Variables]
  [./disp_x]
  [../]
  [./disp_y]
  [../]
  [./strain_zz]
  [../]
[]

[AuxVariables]
  [./temp]
  [../]
  [./nl_strain_zz]
    order = CONSTANT
    family = MONOMIAL
  [../]
[]

[Postprocessors]
  [./react_z]
    type = MaterialTensorIntegral
    rank_two_tensor = stress
    index_i = 2
    index_j = 2
  [../]
  [./min_strain_zz]
    type = NodalExtremeValue
    variable = strain_zz
    value_type = min
  [../]
  [./max_strain_zz]
    type = NodalExtremeValue
    variable = strain_zz
    value_type = max
  [../]
[]

[Modules/TensorMechanics/Master]
  [plane_stress]
    planar_formulation = WEAK_PLANE_STRESS
    strain = FINITE
    generate_output = 'stress_xx stress_xy stress_yy stress_zz strain_xx strain_xy strain_yy'
    eigenstrain_names = eigenstrain
  []
[]

[AuxKernels]
  [./tempfuncaux]
    type = FunctionAux
    variable = temp
    function = tempfunc
    use_displaced_mesh = false
  [../]
  [./strain_zz]
    type = RankTwoAux
    rank_two_tensor = total_strain
    variable = nl_strain_zz
    index_i = 2
    index_j = 2
  [../]
[]

[Functions]
  [./pull]
    type = PiecewiseLinear
    x='0     1   100'
    y='0  0.00  0.00'
  [../]
  [./tempfunc]
    type = ParsedFunction
    value = '(1 - x) * t'
  [../]
[]

[BCs]
  [./bottomx]
    type = DirichletBC
    boundary = 0
    variable = disp_x
    value = 0.0
  [../]
  [./bottomy]
    type = DirichletBC
    boundary = 0
    variable = disp_y
    value = 0.0
  [../]
[]

[Materials]
  [./elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    poissons_ratio = 0.3
    youngs_modulus = 1e6
  [../]
  [./thermal_strain]
    type = ComputeThermalExpansionEigenstrain
    thermal_expansion_coeff = 0.02
    stress_free_temperature = 0.5
    eigenstrain_name = eigenstrain
  [../]
  [./stress]
    type = ComputeFiniteStrainElasticStress
  [../]
[]

[Executioner]
  type = Transient

  solve_type = PJFNK
  line_search = none

# controls for linear iterations
  l_max_its = 100
  l_tol = 1e-06

# controls for nonlinear iterations
  nl_max_its = 15
  nl_rel_tol = 1e-12
  nl_abs_tol = 1e-10

# time control
  start_time = 0.0
  dt = 1.0
  dtmin = 1.0
  end_time = 2.0
[]

[Outputs]
  exodus = true
[]

(modules/tensor_mechanics/test/tests/dynamics/rayleigh_damping/rayleigh_newmark.i)

# Test for rayleigh damping implemented using Newmark time integration

# The test is for an 1D bar element of  unit length fixed on one end
# with a ramped pressure boundary condition applied to the other end.
# zeta and eta correspond to the stiffness and mass proportional rayleigh damping
# beta and gamma are Newmark time integration parameters
# The equation of motion in terms of matrices is:
#
# M*accel + eta*M*vel + zeta*K*vel + K*disp = P*Area
#
# Here M is the mass matrix, K is the stiffness matrix, P is the applied pressure
#
# This equation is equivalent to:
#
# density*accel + eta*density*vel + zeta*d/dt(Div stress) + Div stress = P
#
# The first two terms on the left are evaluated using the Inertial force kernel
# The next two terms on the left involving zeta are evaluated using the
# DynamicStressDivergenceTensors Kernel
# The residual due to Pressure is evaluated using Pressure boundary condition
#
# The system will come to steady state slowly after the pressure becomes constant.

[Mesh]
  type = GeneratedMesh
  dim = 3
  nx = 1
  ny = 1
  nz = 1
  xmin = 0.0
  xmax = 0.1
  ymin = 0.0
  ymax = 1.0
  zmin = 0.0
  zmax = 0.1
[]

[Variables]
  [disp_x]
  []
  [disp_y]
  []
  [disp_z]
  []
[]

[AuxVariables]
  [vel_x]
  []
  [accel_x]
  []
  [vel_y]
  []
  [accel_y]
  []
  [vel_z]
  []
  [accel_z]
  []
  [stress_yy]
    order = CONSTANT
    family = MONOMIAL
  []
  [strain_yy]
    order = CONSTANT
    family = MONOMIAL
  []

[]

[Kernels]
  [DynamicTensorMechanics]
    displacements = 'disp_x disp_y disp_z'
    stiffness_damping_coefficient = 0.1
  []
  [inertia_x]
    type = InertialForce
    variable = disp_x
    velocity = vel_x
    acceleration = accel_x
    beta = 0.25
    gamma = 0.5
    eta = 0.1
  []
  [inertia_y]
    type = InertialForce
    variable = disp_y
    velocity = vel_y
    acceleration = accel_y
    beta = 0.25
    gamma = 0.5
    eta = 0.1
  []
  [inertia_z]
    type = InertialForce
    variable = disp_z
    velocity = vel_z
    acceleration = accel_z
    beta = 0.25
    gamma = 0.5
    eta = 0.1
  []

[]

[AuxKernels]
  [accel_x]
    type = NewmarkAccelAux
    variable = accel_x
    displacement = disp_x
    velocity = vel_x
    beta = 0.25
    execute_on = timestep_end
  []
  [vel_x]
    type = NewmarkVelAux
    variable = vel_x
    acceleration = accel_x
    gamma = 0.5
    execute_on = timestep_end
  []
  [accel_y]
    type = NewmarkAccelAux
    variable = accel_y
    displacement = disp_y
    velocity = vel_y
    beta = 0.25
    execute_on = timestep_end
  []
  [vel_y]
    type = NewmarkVelAux
    variable = vel_y
    acceleration = accel_y
    gamma = 0.5
    execute_on = timestep_end
  []
  [accel_z]
    type = NewmarkAccelAux
    variable = accel_z
    displacement = disp_z
    velocity = vel_z
    beta = 0.25
    execute_on = timestep_end
  []
  [vel_z]
    type = NewmarkVelAux
    variable = vel_z
    acceleration = accel_z
    gamma = 0.5
    execute_on = timestep_end
  []
  [stress_yy]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_yy
    index_i = 1
    index_j = 1
  []
  [strain_yy]
    type = RankTwoAux
    rank_two_tensor = total_strain
    variable = strain_yy
    index_i = 1
    index_j = 1
  []

[]

[BCs]
  [top_y]
    type = DirichletBC
    variable = disp_y
    boundary = top
    value = 0.0
  []
  [top_x]
    type = DirichletBC
    variable = disp_x
    boundary = top
    value = 0.0
  []
  [top_z]
    type = DirichletBC
    variable = disp_z
    boundary = top
    value = 0.0
  []
  [bottom_x]
    type = DirichletBC
    variable = disp_x
    boundary = bottom
    value = 0.0
  []
  [bottom_z]
    type = DirichletBC
    variable = disp_z
    boundary = bottom
    value = 0.0
  []
  [Pressure]
    [Side1]
      boundary = bottom
      function = pressure
      factor = 1
      displacements = 'disp_x disp_y disp_z'
    []
  []
[]

[Materials]
  [Elasticity_tensor]
    type = ComputeElasticityTensor
    block = 0
    fill_method = symmetric_isotropic
    C_ijkl = '210e9 0'
  []

  [strain]
    type = ComputeSmallStrain
    block = 0
    displacements = 'disp_x disp_y disp_z'
  []

  [stress]
    type = ComputeLinearElasticStress
    block = 0
  []

  [density]
    type = GenericConstantMaterial
    block = 0
    prop_names = 'density'
    prop_values = '7750'
  []

[]

[Executioner]
  type = Transient
  start_time = 0
  end_time = 2
  dt = 0.1
[]

[Functions]
  [pressure]
    type = PiecewiseLinear
    x = '0.0 0.1 0.2 1.0 2.0 5.0'
    y = '0.0 0.1 0.2 1.0 1.0 1.0'
    scale_factor = 1e9
  []
[]

[Postprocessors]
  [_dt]
    type = TimestepSize
  []
  [disp]
    type = NodalExtremeValue
    variable = disp_y
    boundary = bottom
  []
  [vel]
    type = NodalExtremeValue
    variable = vel_y
    boundary = bottom
  []
  [accel]
    type = NodalExtremeValue
    variable = accel_y
    boundary = bottom
  []
  [stress_yy]
    type = ElementAverageValue
    variable = stress_yy
  []
  [strain_yy]
    type = ElementAverageValue
    variable = strain_yy
  []

[]

[Outputs]
  exodus = true
  perf_graph = true
[]

(modules/tensor_mechanics/test/tests/ad_thermal_expansion_function/finite_const.i)

# This tests the thermal expansion coefficient function using both
# options to specify that function: mean and instantaneous.  There
# two blocks, each containing a single element, and these use the
# two variants of the function.

# In this test, the instantaneous CTE function has a constant value,
# while the mean CTE function is an analytic function designed to
# give the same response.  If \bar{alpha}(T) is the mean CTE function,
# and \alpha(T) is the instantaneous CTE function,

# \bar{\alpha}(T) = 1/(T-Tref) \intA^{T}_{Tsf} \alpha(T) dT

# where Tref is the reference temperature used to define the mean CTE
# function, and Tsf is the stress-free temperature.

# This version of the test uses finite deformation theory.
# The two models produce very similar results.  There are slight
# differences due to the large deformation treatment.

[Mesh]
  file = 'blocks.e'
[]

[GlobalParams]
  displacements = 'disp_x disp_y disp_z'
[]

[AuxVariables]
  [./temp]
    order = FIRST
    family = LAGRANGE
  [../]
[]

[Modules/TensorMechanics/Master]
  [./all]
    strain = FINITE
    add_variables = true
    eigenstrain_names = eigenstrain
    generate_output = 'strain_xx strain_yy strain_zz'
    use_automatic_differentiation = true
  [../]
[]

[BCs]
  [./left]
    type = DirichletBC
    variable = disp_x
    boundary = 3
    value = 0.0
  [../]

  [./bottom]
    type = DirichletBC
    variable = disp_y
    boundary = 2
    value = 0.0
  [../]

  [./back]
    type = DirichletBC
    variable = disp_z
    boundary = 1
    value = 0.0
  [../]
[]

[AuxKernels]
  [./temp]
    type = FunctionAux
    variable = temp
    block = '1 2'
    function = temp_func
  [../]
[]

[Materials]
  [./elasticity_tensor]
    type = ADComputeIsotropicElasticityTensor
    youngs_modulus = 1e6
    poissons_ratio = 0.3
  [../]
  [./small_stress]
    type = ADComputeFiniteStrainElasticStress
  [../]
  [./thermal_expansion_strain1]
    type = ADComputeMeanThermalExpansionFunctionEigenstrain
    block = 1
    thermal_expansion_function = cte_func_mean
    thermal_expansion_function_reference_temperature = 0.5
    stress_free_temperature = 0.0
    temperature = temp
    eigenstrain_name = eigenstrain
  [../]
  [./thermal_expansion_strain2]
    type = ADComputeInstantaneousThermalExpansionFunctionEigenstrain
    block = 2
    thermal_expansion_function = cte_func_inst
    stress_free_temperature = 0.0
    temperature = temp
    eigenstrain_name = eigenstrain
  [../]
[]

[Functions]
  [./cte_func_mean]
    type = ParsedFunction
    vars = 'tsf tref scale' #stress free temp, reference temp, scale factor
    vals = '0.0 0.5  1e-4'
    value = 'scale * (t - tsf) / (t - tref)'
  [../]
  [./cte_func_inst]
    type = PiecewiseLinear
    xy_data = '0 1.0
               2 1.0'
    scale_factor = 1e-4
  [../]

  [./temp_func]
    type = PiecewiseLinear
    xy_data = '0 1
               1 2'
  [../]
[]

[Postprocessors]
  [./disp_1]
    type = NodalExtremeValue
    variable = disp_x
    boundary = 101
  [../]

  [./disp_2]
    type = NodalExtremeValue
    variable = disp_x
    boundary = 102
  [../]
[]

[Executioner]
  type = Transient

  solve_type = NEWTON
  l_max_its = 100
  l_tol = 1e-4
  nl_abs_tol = 1e-8
  nl_rel_tol = 1e-12

  start_time = 0.0
  end_time = 1.0
  dt = 0.1
[]

[Outputs]
  csv = true
[]

(modules/tensor_mechanics/test/tests/plane_stress/ad_weak_plane_stress_finite.i)

[GlobalParams]
  displacements = 'disp_x disp_y'
  temperature = temp
  out_of_plane_strain = strain_zz
[]

[Mesh]
  [./square]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 2
    ny = 2
  [../]
[]

[Variables]
  [./disp_x]
  [../]
  [./disp_y]
  [../]
  [./strain_zz]
  [../]
[]

[AuxVariables]
  [./temp]
  [../]
  [./nl_strain_zz]
    order = CONSTANT
    family = MONOMIAL
  [../]
[]

[Postprocessors]
  [./react_z]
    type = ADMaterialTensorIntegral
    rank_two_tensor = stress
    index_i = 2
    index_j = 2
  [../]
  [./min_strain_zz]
    type = NodalExtremeValue
    variable = strain_zz
    value_type = min
  [../]
  [./max_strain_zz]
    type = NodalExtremeValue
    variable = strain_zz
    value_type = max
  [../]
[]

[Modules/TensorMechanics/Master]
  [./plane_stress]
    planar_formulation = WEAK_PLANE_STRESS
    strain = FINITE
    generate_output = 'stress_xx stress_xy stress_yy stress_zz strain_xx strain_xy strain_yy'
    eigenstrain_names = eigenstrain
    use_automatic_differentiation = true
  [../]
[]

[AuxKernels]
  [./tempfuncaux]
    type = FunctionAux
    variable = temp
    function = tempfunc
    use_displaced_mesh = false
  [../]
  [./strain_zz]
    type = ADRankTwoAux
    rank_two_tensor = total_strain
    variable = nl_strain_zz
    index_i = 2
    index_j = 2
  [../]
[]

[Functions]
  [./pull]
    type = PiecewiseLinear
    x='0     1   100'
    y='0  0.00  0.00'
  [../]
  [./tempfunc]
    type = ParsedFunction
    value = '(1 - x) * t'
  [../]
[]

[BCs]
  [./bottomx]
    type = DirichletBC
    boundary = 0
    variable = disp_x
    value = 0.0
  [../]
  [./bottomy]
    type = DirichletBC
    boundary = 0
    variable = disp_y
    value = 0.0
  [../]
[]

[Materials]
  [./elasticity_tensor]
    type = ADComputeIsotropicElasticityTensor
    poissons_ratio = 0.3
    youngs_modulus = 1e6
  [../]
  [./thermal_strain]
    type = ADComputeThermalExpansionEigenstrain
    thermal_expansion_coeff = 0.02
    stress_free_temperature = 0.5
    eigenstrain_name = eigenstrain
  [../]
  [./stress]
    type = ADComputeFiniteStrainElasticStress
  [../]
[]

[Executioner]
  type = Transient

  solve_type = 'NEWTON'

  # controls for nonlinear iterations
  nl_max_its = 15
  nl_rel_tol = 1e-14
  nl_abs_tol = 1e-12

  # time control
  start_time = 0.0
  dt = 1.0
  dtmin = 1.0
  end_time = 2.0
[]

[Outputs]
  file_base = 'weak_plane_stress_finite_out'
  exodus = true
[]

(modules/tensor_mechanics/test/tests/dynamics/rayleigh_damping/rayleigh_newmark_material_dependent.i)

# Test for rayleigh damping implemented using Newmark time integration

# The test is for an 1D bar element of  unit length fixed on one end
# with a ramped pressure boundary condition applied to the other end.
# zeta and eta correspond to the stiffness and mass proportional rayleigh damping
# beta and gamma are Newmark time integration parameters
# The equation of motion in terms of matrices is:
#
# M*accel + eta*M*vel + zeta*K*vel + K*disp = P*Area
#
# Here M is the mass matrix, K is the stiffness matrix, P is the applied pressure
#
# This equation is equivalent to:
#
# density*accel + eta*density*vel + zeta*d/dt(Div stress) + Div stress = P
#
# The first two terms on the left are evaluated using the Inertial force kernel
# The next two terms on the left involving zeta are evaluated using the
# DynamicStressDivergenceTensors Kernel
# The residual due to Pressure is evaluated using Pressure boundary condition
#
# The system will come to steady state slowly after the pressure becomes constant.

[Mesh]
  type = GeneratedMesh
  dim = 3
  nx = 1
  ny = 1
  nz = 1
  xmin = 0.0
  xmax = 0.1
  ymin = 0.0
  ymax = 1.0
  zmin = 0.0
  zmax = 0.1
[]

[Variables]
  [disp_x]
  []
  [disp_y]
  []
  [disp_z]
  []
[]

[AuxVariables]
  [vel_x]
  []
  [accel_x]
  []
  [vel_y]
  []
  [accel_y]
  []
  [vel_z]
  []
  [accel_z]
  []
  [stress_yy]
    order = CONSTANT
    family = MONOMIAL
  []
  [strain_yy]
    order = CONSTANT
    family = MONOMIAL
  []

[]

[Kernels]
  [DynamicTensorMechanics]
    displacements = 'disp_x disp_y disp_z'
    stiffness_damping_coefficient = 'zeta_rayleigh'
  []
  [inertia_x]
    type = InertialForce
    variable = disp_x
    velocity = vel_x
    acceleration = accel_x
    beta = 0.25
    gamma = 0.5
    eta = 'eta_rayleigh'
  []
  [inertia_y]
    type = InertialForce
    variable = disp_y
    velocity = vel_y
    acceleration = accel_y
    beta = 0.25
    gamma = 0.5
    eta = 'eta_rayleigh'
  []
  [inertia_z]
    type = InertialForce
    variable = disp_z
    velocity = vel_z
    acceleration = accel_z
    beta = 0.25
    gamma = 0.5
    eta = 'eta_rayleigh'
  []
[]

[AuxKernels]
  [accel_x]
    type = NewmarkAccelAux
    variable = accel_x
    displacement = disp_x
    velocity = vel_x
    beta = 0.25
    execute_on = timestep_end
  []
  [vel_x]
    type = NewmarkVelAux
    variable = vel_x
    acceleration = accel_x
    gamma = 0.5
    execute_on = timestep_end
  []
  [accel_y]
    type = NewmarkAccelAux
    variable = accel_y
    displacement = disp_y
    velocity = vel_y
    beta = 0.25
    execute_on = timestep_end
  []
  [vel_y]
    type = NewmarkVelAux
    variable = vel_y
    acceleration = accel_y
    gamma = 0.5
    execute_on = timestep_end
  []
  [accel_z]
    type = NewmarkAccelAux
    variable = accel_z
    displacement = disp_z
    velocity = vel_z
    beta = 0.25
    execute_on = timestep_end
  []
  [vel_z]
    type = NewmarkVelAux
    variable = vel_z
    acceleration = accel_z
    gamma = 0.5
    execute_on = timestep_end
  []
  [stress_yy]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_yy
    index_i = 1
    index_j = 1
  []
  [strain_yy]
    type = RankTwoAux
    rank_two_tensor = total_strain
    variable = strain_yy
    index_i = 1
    index_j = 1
  []
[]

[BCs]
  [top_y]
    type = DirichletBC
    variable = disp_y
    boundary = top
    value = 0.0
  []
  [top_x]
    type = DirichletBC
    variable = disp_x
    boundary = top
    value = 0.0
  []
  [top_z]
    type = DirichletBC
    variable = disp_z
    boundary = top
    value = 0.0
  []
  [bottom_x]
    type = DirichletBC
    variable = disp_x
    boundary = bottom
    value = 0.0
  []
  [bottom_z]
    type = DirichletBC
    variable = disp_z
    boundary = bottom
    value = 0.0
  []
  [Pressure]
    [Side1]
      boundary = bottom
      function = pressure
      displacements = 'disp_x disp_y disp_z'
      factor = 1
    []
  []
[]

[Materials]
  [Elasticity_tensor]
    type = ComputeElasticityTensor
    block = 0
    fill_method = symmetric_isotropic
    C_ijkl = '210e9 0'
  []

  [strain]
    type = ComputeSmallStrain
    block = 0
    displacements = 'disp_x disp_y disp_z'
  []

  [stress]
    type = ComputeLinearElasticStress
    block = 0
  []

  [density]
    type = GenericConstantMaterial
    block = 0
    prop_names = 'density'
    prop_values = '7750'
  []
  [material_zeta]
    type = GenericConstantMaterial
    block = 0
    prop_names = 'zeta_rayleigh'
    prop_values = '0.1'
  []
  [material_eta]
    type = GenericConstantMaterial
    block = 0
    prop_names = 'eta_rayleigh'
    prop_values = '0.1'
  []
[]

[Executioner]
  type = Transient
  start_time = 0
  end_time = 2
  dt = 0.1
[]

[Functions]
  [pressure]
    type = PiecewiseLinear
    x = '0.0 0.1 0.2 1.0 2.0 5.0'
    y = '0.0 0.1 0.2 1.0 1.0 1.0'
    scale_factor = 1e9
  []
[]

[Postprocessors]
  [_dt]
    type = TimestepSize
  []
  [disp]
    type = NodalExtremeValue
    variable = disp_y
    boundary = bottom
  []
  [vel]
    type = NodalExtremeValue
    variable = vel_y
    boundary = bottom
  []
  [accel]
    type = NodalExtremeValue
    variable = accel_y
    boundary = bottom
  []
  [stress_yy]
    type = ElementAverageValue
    variable = stress_yy
  []
  [strain_yy]
    type = ElementAverageValue
    variable = strain_yy
  []

[]

[Outputs]
  file_base = 'rayleigh_newmark_out'
  exodus = true
  perf_graph = true
[]

(modules/tensor_mechanics/test/tests/dynamics/time_integration/newmark_action.i)

# Test for  Newmark time integration

# The test is for an 1D bar element of  unit length fixed on one end
# with a ramped pressure boundary condition applied to the other end.
# beta and gamma are Newmark time integration parameters
# The equation of motion in terms of matrices is:
#
# M*accel + K*disp = P*Area
#
# Here M is the mass matrix, K is the stiffness matrix, P is the applied pressure
#
# This equation is equivalent to:
#
# density*accel + Div Stress = P
#
# The first term on the left is evaluated using the Inertial force kernel
# The last term on the left is evaluated using StressDivergenceTensors
# The residual due to Pressure is evaluated using Pressure boundary condition

[Mesh]
  type = GeneratedMesh
  dim = 3
  xmax = 0.1
  ymax = 1.0
  zmax = 0.1
  use_displaced_mesh = false
[]

[GlobalParams]
  displacements = 'disp_x disp_y disp_z'
[]

[Modules/TensorMechanics/DynamicMaster]
  [all]
    add_variables = true
    newmark_beta = 0.25
    newmark_gamma = 0.5
    strain = SMALL
    density = 7750
    generate_output = 'stress_yy strain_yy'
  []
[]

[BCs]
  [top_x]
    type = DirichletBC
    variable = disp_x
    boundary = top
    value = 0.0
  []
  [top_y]
    type = DirichletBC
    variable = disp_y
    boundary = top
    value = 0.0
  []
  [top_z]
    type = DirichletBC
    variable = disp_z
    boundary = top
    value = 0.0
  []
  [Pressure]
    [Side1]
      boundary = bottom
      function = pressure
      factor = 1
    []
  []
[]

[Materials]
  [Elasticity_tensor]
    type = ComputeElasticityTensor
    fill_method = symmetric_isotropic
    C_ijkl = '210 0'
  []
  [stress]
    type = ComputeLinearElasticStress
  []
[]

[Executioner]
  type = Transient
  start_time = 0
  end_time = 2
  dt = 0.1
[]

[Functions]
  [pressure]
    type = PiecewiseLinear
    x = '0.0 0.2 1.0 5.0'
    y = '0.0 0.2 1.0 1.0'
    scale_factor = 1e3
  []
[]

[Postprocessors]
  [dt]
    type = TimestepSize
  []
  [disp]
    type = NodalExtremeValue
    variable = disp_y
    boundary = bottom
  []
  [vel]
    type = NodalExtremeValue
    variable = vel_y
    boundary = bottom
  []
  [accel]
    type = NodalExtremeValue
    variable = accel_y
    boundary = bottom
  []
  [stress_yy]
    type = ElementAverageValue
    variable = stress_yy
  []
  [strain_yy]
    type = ElementAverageValue
    variable = strain_yy
  []

[]

[Outputs]
  exodus = true
  perf_graph = true
[]

(modules/stochastic_tools/examples/surrogates/sub.i)

[Mesh]
  type = GeneratedMesh
  dim = 1
  nx = 100
  xmax = 1
  elem_type = EDGE3
[]

[Variables]
  [T]
    order = SECOND
    family = LAGRANGE
  []
[]

[Kernels]
  [diffusion]
    type = MatDiffusion
    variable = T
    diffusivity = k
  []
  [source]
    type = BodyForce
    variable = T
    value = 1.0
  []
[]

[Materials]
  [conductivity]
    type = GenericConstantMaterial
    prop_names = k
    prop_values = 2.0
  []
[]

[BCs]
  [right]
    type = DirichletBC
    variable = T
    boundary = right
    value = 300
  []
[]

[Executioner]
  type = Steady
  solve_type = PJFNK
  petsc_options_iname = '-pc_type -pc_hypre_type'
  petsc_options_value = 'hypre boomeramg'
[]

[Postprocessors]
  [avg]
    type = AverageNodalVariableValue
    variable = T
  []
  [max]
    type = NodalExtremeValue
    variable = T
    value_type = max
  []
[]

[Outputs]
[]

(modules/tensor_mechanics/test/tests/volumetric_eigenstrain/volumetric_eigenstrain.i)

# This tests the ability of the ComputeVolumetricEigenstrain material
# to compute an eigenstrain tensor that results in a solution that exactly
# recovers the specified volumetric expansion.
# This model applies volumetric strain that ramps from 0 to 2 to a unit cube
# and computes the final volume, which should be exactly 3.  Note that the default
# TaylorExpansion option for decomposition_method gives a small (~4%) error
# with this very large incremental strain, but decomposition_method=EigenSolution
# gives the exact solution.

[Mesh]
  type = GeneratedMesh
  dim = 3
  nx = 1
  ny = 1
  nz = 1
[]

[Variables]
  [./disp_x]
  [../]
  [./disp_y]
  [../]
  [./disp_z]
  [../]
[]

[AuxVariables]
  [./volumetric_strain]
    order = CONSTANT
    family = MONOMIAL
  [../]
[]

[GlobalParams]
  displacements = 'disp_x disp_y disp_z'
[]

[Modules/TensorMechanics/Master]
  [./master]
    strain = FINITE
    eigenstrain_names = eigenstrain
    decomposition_method = EigenSolution #Necessary for exact solution
  [../]
[]

[AuxKernels]
  [./volumetric_strain]
    type = RankTwoScalarAux
    scalar_type = VolumetricStrain
    rank_two_tensor = total_strain
    variable = volumetric_strain
  [../]
[]

[BCs]
  [./left]
    type = DirichletBC
    variable = disp_x
    boundary = left
    value = 0.0
  [../]
  [./bottom]
    type = DirichletBC
    variable = disp_y
    boundary = bottom
    value = 0.0
  [../]
  [./back]
    type = DirichletBC
    variable = disp_z
    boundary = back
    value = 0.0
  [../]
[]

[Materials]
  [./elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 1e6
    poissons_ratio = 0.3
  [../]
  [./finite_strain_stress]
    type = ComputeFiniteStrainElasticStress
  [../]
  [./volumetric_eigenstrain]
    type = ComputeVolumetricEigenstrain
    volumetric_materials = volumetric_change
    eigenstrain_name = eigenstrain
    args = ''
  [../]
  [./volumetric_change]
    type = GenericFunctionMaterial
    prop_names = volumetric_change
    prop_values = t
  [../]
[]

[Postprocessors]
  [./vol]
    type = VolumePostprocessor
    use_displaced_mesh = true
    execute_on = 'initial timestep_end'
  [../]
  [./volumetric_strain]
    type = ElementalVariableValue
    variable = volumetric_strain
    elementid = 0
  [../]
  [./disp_right]
    type = NodalExtremeValue
    variable = disp_x
    boundary = right
  [../]
[]

[Executioner]
  type = Transient
  solve_type = PJFNK
  l_max_its = 100
  l_tol = 1e-4
  nl_abs_tol = 1e-8
  nl_rel_tol = 1e-12

  start_time = 0.0
  end_time = 2.0
  dt = 1.0
[]

[Outputs]
  exodus = true
  csv = true
[]

(modules/tensor_mechanics/test/tests/dynamics/rayleigh_damping/rayleigh_hht.i)

# Test for rayleigh damping implemented using HHT time integration
#
# The test is for an 1D bar element of  unit length fixed on one end
# with a ramped pressure boundary condition applied to the other end.
# zeta and eta correspond to the stiffness and mass proportional rayleigh damping
# alpha, beta and gamma are HHT time integration parameters
# The equation of motion in terms of matrices is:
#
# M*accel + (eta*M+zeta*K)*[(1+alpha)vel-alpha vel_old]
# + alpha*(K*disp - K*disp_old) + K*disp = P(t+alpha dt)*Area
#
# Here M is the mass matrix, K is the stiffness matrix, P is the applied pressure
#
# This equation is equivalent to:
#
# density*accel + eta*density*[(1+alpha)vel-alpha vel_old]
# + zeta*[(1+alpha)*d/dt(Div stress)- alpha*d/dt(Div stress_old)]
# + alpha *(Div stress - Div stress_old) +Div Stress= P(t+alpha dt)
#
# The first two terms on the left are evaluated using the Inertial force kernel
# The next three terms on the left involving zeta and alpha are evaluated using
# the DynamicStressDivergenceTensors Kernel
# The residual due to Pressure is evaluated using Pressure boundary condition
#
# The system will come to steady state slowly after the pressure becomes constant.
# Alpha equal to zero will result in Newmark integration.

[Mesh]
  type = GeneratedMesh
  dim = 3
  nx = 1
  ny = 1
  nz = 1
  xmin = 0.0
  xmax = 0.1
  ymin = 0.0
  ymax = 1.0
  zmin = 0.0
  zmax = 0.1
  use_displaced_mesh = false
[]

[GlobalParams]
  displacements = 'disp_x disp_y disp_z'
[]

[AuxVariables]
  [stress_yy]
    order = CONSTANT
    family = MONOMIAL
  []
  [strain_yy]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[AuxKernels]
  [stress_yy]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_yy
    index_i = 1
    index_j = 1
  []
  [strain_yy]
    type = RankTwoAux
    rank_two_tensor = total_strain
    variable = strain_yy
    index_i = 1
    index_j = 1
  []
[]

[Modules/TensorMechanics/DynamicMaster]
  [all]
    add_variables = true
    hht_alpha = 0.11
    newmark_beta = 0.25
    newmark_gamma = 0.5
    mass_damping_coefficient = 0.1
    stiffness_damping_coefficient = 0.1
    density = 7750
  []
[]

[BCs]
  [top_y]
    type = DirichletBC
    variable = disp_y
    boundary = top
    value = 0.0
  []
  [top_x]
    type = DirichletBC
    variable = disp_x
    boundary = top
    value = 0.0
  []
  [top_z]
    type = DirichletBC
    variable = disp_z
    boundary = top
    value = 0.0
  []
  [bottom_x]
    type = DirichletBC
    variable = disp_x
    boundary = bottom
    value = 0.0
  []
  [bottom_z]
    type = DirichletBC
    variable = disp_z
    boundary = bottom
    value = 0.0
  []
  [Pressure]
    [Side1]
      boundary = bottom
      function = pressure
      factor = 1
      hht_alpha = 0.11
    []
  []
[]

[Materials]
  [Elasticity_tensor]
    type = ComputeElasticityTensor
    fill_method = symmetric_isotropic
    C_ijkl = '210e9 0'
  []
  [stress]
    type = ComputeLinearElasticStress
  []
[]

[Executioner]
  type = Transient
  start_time = 0
  end_time = 2
  dt = 0.1
[]

[Functions]
  [pressure]
    type = PiecewiseLinear
    x = '0.0 0.1 0.2 1.0 2.0 5.0'
    y = '0.0 0.1 0.2 1.0 1.0 1.0'
    scale_factor = 1e9
  []
[]

[Postprocessors]
  [_dt]
    type = TimestepSize
  []
  [disp]
    type = NodalExtremeValue
    variable = disp_y
    boundary = bottom
  []
  [vel]
    type = NodalExtremeValue
    variable = vel_y
    boundary = bottom
  []
  [accel]
    type = NodalExtremeValue
    variable = accel_y
    boundary = bottom
  []
  [stress_yy]
    type = ElementAverageValue
    variable = stress_yy
  []
  [strain_yy]
    type = ElementAverageValue
    variable = strain_yy
  []
[]

[Outputs]
  exodus = true
  perf_graph = true
[]

(test/tests/misc/block_user_object_check/block_check.i)

[Mesh]
  [./generator]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 10
    ny = 5
  [../]
  [./left_block]
    type = SubdomainBoundingBoxGenerator
    input = generator
    block_id = 1
    bottom_left = '0 0 0'
    top_right = '0.5 1 0'
  [../]
  [./right_block]
    type = SubdomainBoundingBoxGenerator
    input = left_block
    block_id = 2
    bottom_left = '0.5 0 0'
    top_right = '1 1 0'
  [../]
[]

[Variables]
  [./var_1]
    block = 1
    initial_condition = 100
  [../]
  [./var_2]
    block = 2
    initial_condition = 200
  [../]
[]

[Problem]
  type = FEProblem
  kernel_coverage_check = true
  solve = false
[]

[Executioner]
  type = Steady
[]

[Postprocessors]
  [./obj]
    type = NodalExtremeValue
    variable = var_1
    #block = 1 # this is what being tested, see the test spec
    execute_on = 'initial'
  [../]
[]

(modules/tensor_mechanics/test/tests/thermal_expansion_function/finite_linear.i)

# This tests the thermal expansion coefficient function using both
# options to specify that function: mean and instantaneous.  There
# two blocks, each containing a single element, and these use the
# two variants of the function.

# In this test, the instantaneous CTE function is a linear function
# while the mean CTE function is an analytic function designed to
# give the same response.  If \bar{alpha}(T) is the mean CTE function,
# and \alpha(T) is the instantaneous CTE function,

# \bar{\alpha}(T) = 1/(T-Tref) \intA^{T}_{Tsf} \alpha(T) dT

# where Tref is the reference temperature used to define the mean CTE
# function, and Tsf is the stress-free temperature.

# This version of the test uses finite deformation theory.
# The two models produce very similar results.  There are slight
# differences due to the large deformation treatment.

[Mesh]
  file = 'blocks.e'
[]

[GlobalParams]
  displacements = 'disp_x disp_y disp_z'
[]

[AuxVariables]
  [./temp]
    order = FIRST
    family = LAGRANGE
  [../]
[]

[Modules/TensorMechanics/Master]
  [./all]
    strain = FINITE
    add_variables = true
    eigenstrain_names = eigenstrain
    generate_output = 'strain_xx strain_yy strain_zz'
  [../]
[]

[BCs]
  [./left]
    type = DirichletBC
    variable = disp_x
    boundary = 3
    value = 0.0
  [../]

  [./bottom]
    type = DirichletBC
    variable = disp_y
    boundary = 2
    value = 0.0
  [../]

  [./back]
    type = DirichletBC
    variable = disp_z
    boundary = 1
    value = 0.0
  [../]
[]

[AuxKernels]
  [./temp]
    type = FunctionAux
    variable = temp
    block = '1 2'
    function = temp_func
  [../]
[]

[Materials]
  [./elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 1e6
    poissons_ratio = 0.3
  [../]
  [./small_stress]
    type = ComputeFiniteStrainElasticStress
  [../]
  [./thermal_expansion_strain1]
    type = ComputeMeanThermalExpansionFunctionEigenstrain
    block = 1
    thermal_expansion_function = cte_func_mean
    thermal_expansion_function_reference_temperature = 0.5
    stress_free_temperature = 0.0
    temperature = temp
    eigenstrain_name = eigenstrain
  [../]
  [./thermal_expansion_strain2]
    type = ComputeInstantaneousThermalExpansionFunctionEigenstrain
    block = 2
    thermal_expansion_function = cte_func_inst
    stress_free_temperature = 0.0
    temperature = temp
    eigenstrain_name = eigenstrain
  [../]
[]

[Functions]
  [./cte_func_mean]
    type = ParsedFunction
    vars = 'tsf tref scale' #stress free temp, reference temp, scale factor
    vals = '0.0 0.5  1e-4'
    value = 'scale * (0.5 * t^2 - 0.5 * tsf^2) / (t - tref)'
  [../]
  [./cte_func_inst]
    type = PiecewiseLinear
    xy_data = '0 0.0
               2 2.0'
    scale_factor = 1e-4
  [../]

  [./temp_func]
    type = PiecewiseLinear
    xy_data = '0 1
               1 2'
  [../]
[]

[Postprocessors]
  [./disp_1]
    type = NodalExtremeValue
    variable = disp_x
    boundary = 101
  [../]

  [./disp_2]
    type = NodalExtremeValue
    variable = disp_x
    boundary = 102
  [../]
[]

[Executioner]
  type = Transient

  solve_type = PJFNK
  l_max_its = 100
  l_tol = 1e-4
  nl_abs_tol = 1e-8
  nl_rel_tol = 1e-12

  start_time = 0.0
  end_time = 1.0
  dt = 0.1
[]

[Outputs]
  csv = true
[]

(modules/tensor_mechanics/test/tests/ad_thermal_expansion_function/small_linear.i)

# This tests the thermal expansion coefficient function using both
# options to specify that function: mean and instantaneous.  There
# two blocks, each containing a single element, and these use the
# two variants of the function.

# In this test, the instantaneous CTE function is a linear function
# while the mean CTE function is an analytic function designed to
# give the same response.  If \bar{alpha}(T) is the mean CTE function,
# and \alpha(T) is the instantaneous CTE function,

# \bar{\alpha}(T) = 1/(T-Tref) \intA^{T}_{Tsf} \alpha(T) dT

# where Tref is the reference temperature used to define the mean CTE
# function, and Tsf is the stress-free temperature.

# This version of the test uses small deformation theory.  The results
# from the two models are identical.

[Mesh]
  file = 'blocks.e'
[]

[GlobalParams]
  displacements = 'disp_x disp_y disp_z'
[]

[AuxVariables]
  [./temp]
    order = FIRST
    family = LAGRANGE
  [../]
[]

[Modules/TensorMechanics/Master]
  [./all]
    strain = SMALL
    add_variables = true
    eigenstrain_names = eigenstrain
    generate_output = 'strain_xx strain_yy strain_zz'
    use_automatic_differentiation = true
  [../]
[]

[BCs]
  [./left]
    type = DirichletBC
    variable = disp_x
    boundary = 3
    value = 0.0
  [../]

  [./bottom]
    type = DirichletBC
    variable = disp_y
    boundary = 2
    value = 0.0
  [../]

  [./back]
    type = DirichletBC
    variable = disp_z
    boundary = 1
    value = 0.0
  [../]
[]

[AuxKernels]
  [./temp]
    type = FunctionAux
    variable = temp
    block = '1 2'
    function = temp_func
  [../]
[]

[Materials]
  [./elasticity_tensor]
    type = ADComputeIsotropicElasticityTensor
    youngs_modulus = 1e6
    poissons_ratio = 0.3
  [../]
  [./small_stress]
    type = ADComputeLinearElasticStress
  [../]
  [./thermal_expansion_strain1]
    type = ADComputeMeanThermalExpansionFunctionEigenstrain
    block = 1
    thermal_expansion_function = cte_func_mean
    thermal_expansion_function_reference_temperature = 0.5
    stress_free_temperature = 0.0
    temperature = temp
    eigenstrain_name = eigenstrain
  [../]
  [./thermal_expansion_strain2]
    type = ADComputeInstantaneousThermalExpansionFunctionEigenstrain
    block = 2
    thermal_expansion_function = cte_func_inst
    stress_free_temperature = 0.0
    temperature = temp
    eigenstrain_name = eigenstrain
  [../]
[]

[Functions]
  [./cte_func_mean]
    type = ParsedFunction
    vars = 'tsf tref scale' #stress free temp, reference temp, scale factor
    vals = '0.0 0.5  1e-4'
    value = 'scale * (0.5 * t^2 - 0.5 * tsf^2) / (t - tref)'
  [../]
  [./cte_func_inst]
    type = PiecewiseLinear
    xy_data = '0 0.0
               2 2.0'
    scale_factor = 1e-4
  [../]

  [./temp_func]
    type = PiecewiseLinear
    xy_data = '0 1
               1 2'
  [../]
[]

[Postprocessors]
  [./disp_1]
    type = NodalExtremeValue
    variable = disp_x
    boundary = 101
  [../]

  [./disp_2]
    type = NodalExtremeValue
    variable = disp_x
    boundary = 102
  [../]
[]

[Executioner]
  type = Transient

  solve_type = PJFNK
  l_max_its = 100
  l_tol = 1e-4
  nl_abs_tol = 1e-8
  nl_rel_tol = 1e-12

  start_time = 0.0
  end_time = 1.0
  dt = 0.1
[]

[Outputs]
  csv = true
[]

(test/tests/postprocessors/nodal_extreme_value/nodal_extreme_pps_test.i)

[Mesh]
  type = FileMesh
  file = trapezoid.e
  uniform_refine = 1
  # This test will not work in parallel with DistributedMesh enabled
  # due to a bug in PeriodicBCs.
  parallel_type = replicated
[]

[Functions]
  [./tr_x]
    type = ParsedFunction
    value = -x*cos(pi/3)
  [../]
  [./tr_y]
    type = ParsedFunction
    value = x*sin(pi/3)
  [../]
  [./itr_x]
    type = ParsedFunction
    value = -x/cos(pi/3)
  [../]
  [./itr_y]
    type = ParsedFunction
    value = 0
  [../]
[]

[Variables]
  [./u]
    order = FIRST
    family = LAGRANGE
  [../]
[]

[Kernels]
  [./diff]
    type = Diffusion
    variable = u
  [../]
  [./forcing]
    type = GaussContForcing
    variable = u
    x_center = 2
    y_center = -1
    x_spread = 0.25
    y_spread = 0.5
  [../]
  [./dot]
    type = TimeDerivative
    variable = u
  [../]
[]

[BCs]
  # active = ' '
  [./Periodic]
    [./x]
      primary = 1
      secondary = 4
      transform_func = 'tr_x tr_y'
      inv_transform_func = 'itr_x itr_y'
    [../]
  [../]
[]

[Postprocessors]
  [./max_nodal_pps]
    type = NodalExtremeValue
    variable = u
  [../]
  [./max_node_id]
    type = NodalMaxValueId
    variable = u
  [../]
  [./min_nodal_pps]
    type = NodalExtremeValue
    variable = u
    value_type = min
  [../]
[]

[Executioner]
  type = Transient
  dt = 0.5
  num_steps = 6
[]

[Outputs]
  execute_on = 'timestep_end'
  exodus = true
[]

(modules/combined/test/tests/gap_heat_transfer_mortar/finite-2d/varied_pressure_thermomechanical_mortar.i)

## Units in the input file: m-Pa-s-K

[GlobalParams]
  displacements = 'disp_x disp_y'
[]

[Mesh]
  [left_rectangle]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 20
    ny = 10
    xmax = 0.25
    ymin = 0
    ymax = 0.5
    boundary_name_prefix = moving_block
  []
  [left_block]
    type = SubdomainIDGenerator
    input = left_rectangle
    subdomain_id = 1
  []
  [right_rectangle]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 20
    ny = 10
    xmin = 0.25
    xmax = 0.5
    ymin = 0
    ymax = 0.5
    boundary_name_prefix = fixed_block
    boundary_id_offset = 4
  []
  [right_block]
    type = SubdomainIDGenerator
    input = right_rectangle
    subdomain_id = 2
  []
  [two_blocks]
    type = MeshCollectionGenerator
    inputs = 'left_block right_block'
  []
  [block_rename]
    type = RenameBlockGenerator
    input = two_blocks
    old_block = '1 2'
    new_block = 'left_block right_block'
  []
  patch_update_strategy = iteration
[]

[Variables]
  [disp_x]
    block = 'left_block right_block'
  []
  [disp_y]
    block = 'left_block right_block'
  []
  [temperature]
    initial_condition = 300.0
  []
  [temperature_interface_lm]
    block = 'interface_secondary_subdomain'
  []
[]

[Modules]
  [TensorMechanics/Master]
    [steel]
      strain = FINITE
      add_variables = false
      use_automatic_differentiation = true
      generate_output = 'strain_xx strain_xy strain_yy stress_xx stress_xy stress_yy'
      additional_generate_output = 'vonmises_stress'
      additional_material_output_family = 'MONOMIAL'
      additional_material_output_order = 'FIRST'
      block = 'left_block'
    []
    [aluminum]
      strain = FINITE
      add_variables = false
      use_automatic_differentiation = true
      generate_output = 'strain_xx strain_xy strain_yy stress_xx stress_xy stress_yy'
      additional_generate_output = 'vonmises_stress'
      additional_material_output_family = 'MONOMIAL'
      additional_material_output_order = 'FIRST'
      block = 'right_block'
    []
  []
[]

[Kernels]
  [HeatDiff_steel]
    type = ADHeatConduction
    variable = temperature
    thermal_conductivity = steel_thermal_conductivity
    block = 'left_block'
  []
  [HeatTdot_steel]
    type = ADHeatConductionTimeDerivative
    variable = temperature
    specific_heat = steel_heat_capacity
    density_name = steel_density
    block = 'left_block'
  []
  [HeatDiff_aluminum]
    type = ADHeatConduction
    variable = temperature
    thermal_conductivity = aluminum_thermal_conductivity
    block = 'right_block'
  []
  [HeatTdot_aluminum]
    type = ADHeatConductionTimeDerivative
    variable = temperature
    specific_heat = aluminum_heat_capacity
    density_name = aluminum_density
    block = 'right_block'
  []
[]

[BCs]
  [fixed_bottom_edge]
    type = ADDirichletBC
    variable = disp_y
    value = 0
    boundary = 'moving_block_bottom fixed_block_bottom'
  []
  [fixed_outer_edge]
    type = ADDirichletBC
    variable = disp_x
    value = 0
    boundary = 'fixed_block_right'
  []
  [pressure_left_block]
    type = ADPressure
    variable = disp_x
    component = 0
    boundary = 'moving_block_left'
    function = '1e4*t*y'
  []
  [temperature_left]
    type = ADDirichletBC
    variable = temperature
    value = 300
    boundary = 'moving_block_left'
  []
  [temperature_right]
    type = ADDirichletBC
    variable = temperature
    value = 800
    boundary = 'fixed_block_right'
  []
[]

[Contact]
  [interface]
    primary = moving_block_right
    secondary = fixed_block_left
    model = frictionless
    formulation = mortar
    correct_edge_dropping = true
  []
[]

[Constraints]
  [thermal_contact]
    type = ModularGapConductanceConstraint
    variable = temperature_interface_lm
    secondary_variable = temperature
    primary_boundary = moving_block_right
    primary_subdomain = interface_primary_subdomain
    secondary_boundary = fixed_block_left
    secondary_subdomain = interface_secondary_subdomain
    gap_flux_models = 'closed'
    use_displaced_mesh = true
  []
[]

[Materials]
  [steel_elasticity_tensor]
    type = ADComputeIsotropicElasticityTensor
    youngs_modulus = 1.93e11 #in Pa, 193 GPa, stainless steel 304
    poissons_ratio = 0.29
    block = 'left_block'
  []
  [steel_stress]
    type = ADComputeFiniteStrainElasticStress
    block = 'left_block'
  []
  [steel_thermal_properties]
    type = ADGenericConstantMaterial
    prop_names = 'steel_density steel_thermal_conductivity steel_heat_capacity steel_hardness'
    prop_values = ' 8e3            16.2                     0.5                 129' ## for stainless steel 304
    block = 'left_block'
  []

  [aluminum_elasticity_tensor]
    type = ADComputeIsotropicElasticityTensor
    youngs_modulus = 6.8e10 #in Pa, 68 GPa, aluminum
    poissons_ratio = 0.36
    block = 'right_block'
  []
  [aluminum_stress]
    type = ADComputeFiniteStrainElasticStress
    block = 'right_block'
  []
  [aluminum_thermal_properties]
    type = ADGenericConstantMaterial
    prop_names = 'aluminum_density aluminum_thermal_conductivity aluminum_heat_capacity aluminum_hardness'
    prop_values = ' 2.7e3            210                           0.9                   15' #for 99% pure Al
    block = 'right_block'
  []
[]

[UserObjects]
  [closed]
    type = GapFluxModelPressureDependentConduction
    primary_conductivity = steel_thermal_conductivity
    secondary_conductivity = aluminum_thermal_conductivity
    temperature = temperature
    contact_pressure = interface_normal_lm
    primary_hardness = steel_hardness
    secondary_hardness = aluminum_hardness
    boundary = moving_block_right
  []
[]

[Postprocessors]
  [contact_pressure_max]
    type = NodalExtremeValue
    variable = interface_normal_lm
    block = interface_secondary_subdomain
    value_type = max
  []
  [contact_pressure_average]
    type = AverageNodalVariableValue
    variable = interface_normal_lm
    block = interface_secondary_subdomain
  []
  [contact_pressure_min]
    type = NodalExtremeValue
    variable = interface_normal_lm
    block = interface_secondary_subdomain
    value_type = min
  []
  [interface_temperature_max]
    type = NodalExtremeValue
    variable = temperature
    block = interface_secondary_subdomain
    value_type = max
  []
  [interface_temperature_average]
    type = AverageNodalVariableValue
    variable = temperature
    block = interface_secondary_subdomain
  []
  [interface_temperature_min]
    type = NodalExtremeValue
    variable = temperature
    block = interface_secondary_subdomain
    value_type = min
  []
[]


[Executioner]
  type = Transient
  solve_type = NEWTON
  automatic_scaling = false
  line_search = 'none'

  # mortar contact solver options
  petsc_options = '-snes_converged_reason -pc_svd_monitor'
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_type'
  petsc_options_value = ' lu       superlu_dist'
  snesmf_reuse_base = false

  nl_rel_tol = 1e-7
  nl_max_its = 20
  l_max_its = 50

  dt = 0.125
  end_time = 1
[]

[Outputs]
  csv = true
  perf_graph = true
[]

(test/tests/misc/check_error/double_restrict_uo.i)

[Mesh]
  file = sq-2blk.e
[]

[Variables]
  [./u]
    order = FIRST
    family = LAGRANGE
    block = 1
  [../]
  [./v]
    order = FIRST
    family = LAGRANGE
  [../]
[]

[Kernels]
  [./diff_u]
    type = Diffusion
    variable = u
  [../]
  [./diff_v]
    type = Diffusion
    variable = v
  [../]
[]

[BCs]
  [./left_u]
    type = DirichletBC
    variable = u
    boundary = 6
    value = 0
  [../]
  [./right_u]
    type = NeumannBC
    variable = u
    boundary = 8
    value = 4
  [../]
  [./left_v]
    type = DirichletBC
    variable = v
    boundary = 6
    value = 1
  [../]
  [./right_v]
    type = DirichletBC
    variable = v
    boundary = 3
    value = 6
  [../]
[]

[Postprocessors]
  # This test demonstrates that you can have a block restricted NodalPostprocessor
  [./restricted_max]
    type = NodalExtremeValue
    variable = v
    block = 1   # Block restricted
    boundary = 1 # Boundary restricted
  [../]
[]

[Executioner]
  type = Steady

  solve_type = 'PJFNK'
[]

[Outputs]
  exodus = true
[]

(modules/stochastic_tools/examples/surrogates/polynomial_regression/sub.i)

[Mesh]
  type = GeneratedMesh
  dim = 1
  nx = 100
  xmax = 1
  elem_type = EDGE3
[]

[Variables]
  [T]
    order = SECOND
    family = LAGRANGE
  []
[]

[Kernels]
  [diffusion]
    type = MatDiffusion
    variable = T
    diffusivity = k
  []
  [source]
    type = BodyForce
    variable = T
    value = 1.0
  []
[]

[Materials]
  [conductivity]
    type = GenericConstantMaterial
    prop_names = k
    prop_values = 2.0
  []
[]

[BCs]
  [right]
    type = DirichletBC
    variable = T
    boundary = right
    value = 300
  []
[]

[Executioner]
  type = Steady
  solve_type = PJFNK
  petsc_options_iname = '-pc_type -pc_hypre_type'
  petsc_options_value = 'hypre boomeramg'
[]

[Postprocessors]
  [max]
    type = NodalExtremeValue
    variable = T
    value_type = max
  []
[]

[Outputs]
[]

(modules/stochastic_tools/test/tests/vectorpostprocessors/multiple_stochastic_results/sub.i)

[Mesh]
  type = GeneratedMesh
  dim = 1
  nx = 10
[]

[Variables]
  [u]
  []
[]

[Kernels]
  [diff]
    type = Diffusion
    variable = u
  []
  [time]
    type = TimeDerivative
    variable = u
  []
[]

[BCs]
  [left]
    type = DirichletBC
    variable = u
    boundary = left
    value = 0
  []
  [right]
    type = DirichletBC
    variable = u
    boundary = right
    value = 1
  []
[]

[Executioner]
  type = Transient
  num_steps = 5
  dt = 0.01
  solve_type = PJFNK
  petsc_options_iname = '-pc_type -pc_hypre_type'
  petsc_options_value = 'hypre boomeramg'
[]

[Controls]
  [stochastic]
    type = SamplerReceiver
  []
[]

[Postprocessors]
  [avg]
    type = AverageNodalVariableValue
    variable = u
  []
  [max]
    type = NodalExtremeValue
    value_type = MAX
    variable = u
  []
[]

[Outputs]
[]

(modules/tensor_mechanics/test/tests/plane_stress/ad_weak_plane_stress_incremental.i)

[GlobalParams]
  displacements = 'disp_x disp_y'
  temperature = temp
  out_of_plane_strain = strain_zz
[]

[Mesh]
  [./square]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 2
    ny = 2
  [../]
[]

[Variables]
  [./disp_x]
  [../]
  [./disp_y]
  [../]
  [./strain_zz]
  [../]
[]

[AuxVariables]
  [./temp]
  [../]
  [./nl_strain_zz]
    order = CONSTANT
    family = MONOMIAL
  [../]
[]

[Postprocessors]
  [./react_z]
    type = ADMaterialTensorIntegral
    rank_two_tensor = stress
    index_i = 2
    index_j = 2
  [../]
  [./min_strain_zz]
    type = NodalExtremeValue
    variable = strain_zz
    value_type = min
  [../]
  [./max_strain_zz]
    type = NodalExtremeValue
    variable = strain_zz
    value_type = max
  [../]
[]

[Modules/TensorMechanics/Master]
  [./plane_stress]
    planar_formulation = WEAK_PLANE_STRESS
    strain = SMALL
    incremental = true
    generate_output = 'stress_xx stress_xy stress_yy stress_zz strain_xx strain_xy strain_yy'
    eigenstrain_names = eigenstrain
    use_automatic_differentiation = true
  [../]
[]

[AuxKernels]
  [./tempfuncaux]
    type = FunctionAux
    variable = temp
    function = tempfunc
    use_displaced_mesh = false
  [../]
  [./strain_zz]
    type = ADRankTwoAux
    rank_two_tensor = total_strain
    variable = nl_strain_zz
    index_i = 2
    index_j = 2
  [../]
[]

[Functions]
  [./pull]
    type = PiecewiseLinear
    x='0     1   100'
    y='0  0.00  0.00'
  [../]
  [./tempfunc]
    type = ParsedFunction
    value = '(1 - x) * t'
  [../]
[]

[BCs]
  [./bottomx]
    type = DirichletBC
    boundary = 0
    variable = disp_x
    value = 0.0
  [../]
  [./bottomy]
    type = DirichletBC
    boundary = 0
    variable = disp_y
    value = 0.0
  [../]
[]

[Materials]
  [./elasticity_tensor]
    type = ADComputeIsotropicElasticityTensor
    poissons_ratio = 0.3
    youngs_modulus = 1e6
  [../]
  [./thermal_strain]
    type = ADComputeThermalExpansionEigenstrain
    thermal_expansion_coeff = 0.02
    stress_free_temperature = 0.5
    eigenstrain_name = eigenstrain
  [../]
  [./stress]
    type = ADComputeStrainIncrementBasedStress
  [../]
[]

[Executioner]
  type = Transient

  solve_type = 'NEWTON'

# controls for nonlinear iterations
  nl_max_its = 15
  nl_rel_tol = 1e-14
  nl_abs_tol = 1e-12

# time control
  start_time = 0.0
  dt = 1.0
  dtmin = 1.0
  end_time = 2.0
[]

[Outputs]
  file_base = 'weak_plane_stress_incremental_out'
  exodus = true
[]

(modules/tensor_mechanics/test/tests/plane_stress/weak_plane_stress_small.i)

[GlobalParams]
  order = FIRST
  family = LAGRANGE
  displacements = 'disp_x disp_y'
  temperature = temp
  out_of_plane_strain = strain_zz
[]

[Mesh]
  [./square]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 2
    ny = 2
  [../]
[]

[Variables]
  [./disp_x]
  [../]
  [./disp_y]
  [../]
  [./strain_zz]
  [../]
[]

[AuxVariables]
  [./temp]
  [../]
  [./nl_strain_zz]
    order = CONSTANT
    family = MONOMIAL
  [../]
[]

[Postprocessors]
  [./react_z]
    type = MaterialTensorIntegral
    rank_two_tensor = stress
    index_i = 2
    index_j = 2
  [../]
  [./min_strain_zz]
    type = NodalExtremeValue
    variable = strain_zz
    value_type = min
  [../]
  [./max_strain_zz]
    type = NodalExtremeValue
    variable = strain_zz
    value_type = max
  [../]
[]

[Modules/TensorMechanics/Master]
  [plane_stress]
    planar_formulation = WEAK_PLANE_STRESS
    strain = SMALL
    generate_output = 'stress_xx stress_xy stress_yy stress_zz strain_xx strain_xy strain_yy'
    eigenstrain_names = eigenstrain
  []
[]

[AuxKernels]
  [./tempfuncaux]
    type = FunctionAux
    variable = temp
    function = tempfunc
    use_displaced_mesh = false
  [../]
  [./strain_zz]
    type = RankTwoAux
    rank_two_tensor = total_strain
    variable = nl_strain_zz
    index_i = 2
    index_j = 2
  [../]
[]

[Functions]
  [./pull]
    type = PiecewiseLinear
    x='0     1   100'
    y='0  0.00  0.00'
  [../]
  [./tempfunc]
    type = ParsedFunction
    value = '(1 - x) * t'
  [../]
[]

[BCs]
  [./bottomx]
    type = DirichletBC
    boundary = 0
    variable = disp_x
    value = 0.0
  [../]
  [./bottomy]
    type = DirichletBC
    boundary = 0
    variable = disp_y
    value = 0.0
  [../]
[]

[Materials]
  [./elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    poissons_ratio = 0.3
    youngs_modulus = 1e6
  [../]
  [./thermal_strain]
    type = ComputeThermalExpansionEigenstrain
    thermal_expansion_coeff = 0.02
    stress_free_temperature = 0.5
    eigenstrain_name = eigenstrain
  [../]
  [./stress]
    type = ComputeLinearElasticStress
  [../]
[]

[Executioner]
  type = Transient

  solve_type = PJFNK
  line_search = none

# controls for linear iterations
  l_max_its = 100
  l_tol = 1e-06

# controls for nonlinear iterations
  nl_max_its = 15
  nl_rel_tol = 1e-14
  nl_abs_tol = 1e-12

# time control
  start_time = 0.0
  dt = 1.0
  dtmin = 1.0
  end_time = 2.0
[]

[Outputs]
  exodus = true
[]

(modules/tensor_mechanics/test/tests/thermal_expansion_function/small_const.i)

# This tests the thermal expansion coefficient function using both
# options to specify that function: mean and instantaneous.  There
# two blocks, each containing a single element, and these use the
# two variants of the function.

# In this test, the instantaneous CTE function has a constant value,
# while the mean CTE function is an analytic function designed to
# give the same response.  If \bar{alpha}(T) is the mean CTE function,
# and \alpha(T) is the instantaneous CTE function,

# \bar{\alpha}(T) = 1/(T-Tref) \intA^{T}_{Tsf} \alpha(T) dT

# where Tref is the reference temperature used to define the mean CTE
# function, and Tsf is the stress-free temperature.

# This version of the test uses small deformation theory.  The results
# from the two models are identical.

[Mesh]
  file = 'blocks.e'
[]

[GlobalParams]
  displacements = 'disp_x disp_y disp_z'
[]

[AuxVariables]
  [./temp]
    order = FIRST
    family = LAGRANGE
  [../]
[]

[Modules/TensorMechanics/Master]
  [./all]
    strain = SMALL
    add_variables = true
    eigenstrain_names = eigenstrain
    generate_output = 'strain_xx strain_yy strain_zz'
  [../]
[]

[BCs]
  [./left]
    type = DirichletBC
    variable = disp_x
    boundary = 3
    value = 0.0
  [../]

  [./bottom]
    type = DirichletBC
    variable = disp_y
    boundary = 2
    value = 0.0
  [../]

  [./back]
    type = DirichletBC
    variable = disp_z
    boundary = 1
    value = 0.0
  [../]
[]

[AuxKernels]
  [./temp]
    type = FunctionAux
    variable = temp
    block = '1 2'
    function = temp_func
  [../]
[]

[Materials]
  [./elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 1e6
    poissons_ratio = 0.3
  [../]
  [./small_stress]
    type = ComputeLinearElasticStress
  [../]
  [./thermal_expansion_strain1]
    type = ComputeMeanThermalExpansionFunctionEigenstrain
    block = 1
    thermal_expansion_function = cte_func_mean
    thermal_expansion_function_reference_temperature = 0.5
    stress_free_temperature = 0.0
    temperature = temp
    eigenstrain_name = eigenstrain
  [../]
  [./thermal_expansion_strain2]
    type = ComputeInstantaneousThermalExpansionFunctionEigenstrain
    block = 2
    thermal_expansion_function = cte_func_inst
    stress_free_temperature = 0.0
    temperature = temp
    eigenstrain_name = eigenstrain
  [../]
[]

[Functions]
  [./cte_func_mean]
    type = ParsedFunction
    vars = 'tsf tref scale' #stress free temp, reference temp, scale factor
    vals = '0.0 0.5  1e-4'
    value = 'scale * (t - tsf) / (t - tref)'
  [../]
  [./cte_func_inst]
    type = PiecewiseLinear
    xy_data = '0 1.0
               2 1.0'
    scale_factor = 1e-4
  [../]

  [./temp_func]
    type = PiecewiseLinear
    xy_data = '0 1
               1 2'
  [../]
[]

[Postprocessors]
  [./disp_1]
    type = NodalExtremeValue
    variable = disp_x
    boundary = 101
  [../]

  [./disp_2]
    type = NodalExtremeValue
    variable = disp_x
    boundary = 102
  [../]
[]

[Executioner]
  type = Transient

  solve_type = PJFNK
  l_max_its = 100
  l_tol = 1e-4
  nl_abs_tol = 1e-8
  nl_rel_tol = 1e-12

  start_time = 0.0
  end_time = 1.0
  dt = 0.1
[]

[Outputs]
  csv = true
[]

(modules/tensor_mechanics/test/tests/action/action_multi_eigenstrain_same_conditions.i)

# This tests a thermal expansion coefficient function using defined on both
# blocks.  There two blocks, each containing a single element, and these use
# automatic_eigenstrain_names function of the TensorMechanicsAction to ensure
# the names are passed correctly.

# In this test, the instantaneous CTE function has a constant value,
# while the mean CTE function is an analytic function designed to
# give the same response.  If \bar{alpha}(T) is the mean CTE function,
# and \alpha(T) is the instantaneous CTE function,

# \bar{\alpha}(T) = 1/(T-Tref) \intA^{T}_{Tsf} \alpha(T) dT

# where Tref is the reference temperature used to define the mean CTE
# function, and Tsf is the stress-free temperature.

# This version of the test uses finite deformation theory.
# The two models produce very similar results.  There are slight
# differences due to the large deformation treatment.

[Mesh]
  file = 'blocks.e'
[]

[GlobalParams]
  displacements = 'disp_x disp_y disp_z'
[]

[AuxVariables]
  [./temp]
    order = FIRST
    family = LAGRANGE
  [../]
[]

[Problem]
  solve = false
[]

[Modules/TensorMechanics/Master]
  [./block1]
    block = 1
    strain = FINITE
    add_variables = true
    automatic_eigenstrain_names = true
    generate_output = 'strain_xx strain_yy strain_zz'
  [../]
  [./block2]
    block = 2
    strain = FINITE
    add_variables = true
    automatic_eigenstrain_names = true
    generate_output = 'strain_xx strain_yy strain_zz'
  [../]
[]

[BCs]
  [./left]
    type = DirichletBC
    variable = disp_x
    boundary = 3
    value = 0.0
  [../]

  [./bottom]
    type = DirichletBC
    variable = disp_y
    boundary = 2
    value = 0.0
  [../]

  [./back]
    type = DirichletBC
    variable = disp_z
    boundary = 1
    value = 0.0
  [../]
[]

[AuxKernels]
  [./temp]
    type = FunctionAux
    variable = temp
    block = '1 2'
    function = temp_func
  [../]
[]

[Materials]
  [./elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 1e6
    poissons_ratio = 0.3
  [../]
  [./small_stress]
    type = ComputeFiniteStrainElasticStress
  [../]
  [./thermal_expansion_strain1]
    type = ComputeMeanThermalExpansionFunctionEigenstrain
    block = '1 2'
    thermal_expansion_function = cte_func_mean
    thermal_expansion_function_reference_temperature = 0.5
    stress_free_temperature = 0.0
    temperature = temp
    eigenstrain_name = eigenstrain
  [../]
[]

[Functions]
  [./cte_func_mean]
    type = ParsedFunction
    vars = 'tsf tref scale' #stress free temp, reference temp, scale factor
    vals = '0.0 0.5  1e-4'
    value = 'scale * (t - tsf) / (t - tref)'
  [../]
  [./cte_func_inst]
    type = PiecewiseLinear
    xy_data = '0 1.0
               2 1.0'
    scale_factor = 1e-4
  [../]

  [./temp_func]
    type = PiecewiseLinear
    xy_data = '0 1
               1 2'
  [../]
[]

[Postprocessors]
  [./disp_1]
    type = NodalExtremeValue
    variable = disp_x
    boundary = 101
  [../]

  [./disp_2]
    type = NodalExtremeValue
    variable = disp_x
    boundary = 102
  [../]
[]

[Executioner]
  type = Transient

  solve_type = PJFNK
  l_max_its = 100
  l_tol = 1e-4
  nl_abs_tol = 1e-8
  nl_rel_tol = 1e-12

  start_time = 0.0
  end_time = 1.0
  dt = 0.1
[]

[Outputs]
  csv = true
[]

(modules/tensor_mechanics/test/tests/action/ad_converter_action_multi_eigenstrain.i)

# This tests the thermal expansion coefficient function using both
# options to specify that function: mean and instantaneous.  There
# two blocks, each containing a single element, and these use the
# two variants of the function.

# In this test, the instantaneous CTE function has a constant value,
# while the mean CTE function is an analytic function designed to
# give the same response.  If \bar{alpha}(T) is the mean CTE function,
# and \alpha(T) is the instantaneous CTE function,

# \bar{\alpha}(T) = 1/(T-Tref) \intA^{T}_{Tsf} \alpha(T) dT

# where Tref is the reference temperature used to define the mean CTE
# function, and Tsf is the stress-free temperature.

# This version of the test uses finite deformation theory.
# The two models produce very similar results.  There are slight
# differences due to the large deformation treatment.

[Mesh]
  file = 'blocks.e'
[]

[GlobalParams]
  displacements = 'disp_x disp_y disp_z'
[]

[AuxVariables]
  [./temp]
    order = FIRST
    family = LAGRANGE
  [../]
[]

[Problem]
  solve = false
[]

[Modules/TensorMechanics/Master]
  [./block1]
    block = 1
    strain = FINITE
    add_variables = true
    automatic_eigenstrain_names = true
    generate_output = 'strain_xx strain_yy strain_zz'
    use_automatic_differentiation = true
  [../]
  [./block2]
    block = 2
    strain = FINITE
    add_variables = true
    automatic_eigenstrain_names = true
    generate_output = 'strain_xx strain_yy strain_zz'
    use_automatic_differentiation = true
  [../]
[]

[BCs]
  [./left]
    type = DirichletBC
    variable = disp_x
    boundary = 3
    value = 0.0
  [../]

  [./bottom]
    type = DirichletBC
    variable = disp_y
    boundary = 2
    value = 0.0
  [../]

  [./back]
    type = DirichletBC
    variable = disp_z
    boundary = 1
    value = 0.0
  [../]
[]

[AuxKernels]
  [./temp]
    type = FunctionAux
    variable = temp
    block = '1 2'
    function = temp_func
  [../]
[]

[Materials]
  [./elasticity_tensor]
    type = ADComputeIsotropicElasticityTensor
    youngs_modulus = 1e6
    poissons_ratio = 0.3
  [../]
  [./small_stress]
    type = ADComputeFiniteStrainElasticStress
  [../]
  [./thermal_expansion_strain1]
    type = ComputeMeanThermalExpansionFunctionEigenstrain
    block = 1
    thermal_expansion_function = cte_func_mean
    thermal_expansion_function_reference_temperature = 0.5
    stress_free_temperature = 0.0
    temperature = temp
    eigenstrain_name = reg_eigenstrain1
  [../]
  [./converter1]
    type = RankTwoTensorMaterialADConverter
    block = 1
    reg_props_in = 'reg_eigenstrain1'
    ad_props_out = 'eigenstrain1'

  [../]
  [./thermal_expansion_strain2]
    type = ADComputeInstantaneousThermalExpansionFunctionEigenstrain
    block = 2
    thermal_expansion_function = cte_func_inst
    stress_free_temperature = 0.0
    temperature = temp
    eigenstrain_name = eigenstrain2
  [../]
[]

[Functions]
  [./cte_func_mean]
    type = ParsedFunction
    vars = 'tsf tref scale' #stress free temp, reference temp, scale factor
    vals = '0.0 0.5  1e-4'
    value = 'scale * (t - tsf) / (t - tref)'
  [../]
  [./cte_func_inst]
    type = PiecewiseLinear
    xy_data = '0 1.0
               2 1.0'
    scale_factor = 1e-4
  [../]

  [./temp_func]
    type = PiecewiseLinear
    xy_data = '0 1
               1 2'
  [../]
[]

[Postprocessors]
  [./disp_1]
    type = NodalExtremeValue
    variable = disp_x
    boundary = 101
  [../]

  [./disp_2]
    type = NodalExtremeValue
    variable = disp_x
    boundary = 102
  [../]
[]

[Executioner]
  type = Transient

  solve_type = PJFNK
  l_max_its = 100
  l_tol = 1e-4
  nl_abs_tol = 1e-8
  nl_rel_tol = 1e-12

  start_time = 0.0
  end_time = 1.0
  dt = 0.1
[]

[Outputs]
  csv = true
[]

(modules/tensor_mechanics/test/tests/action/action_multi_eigenstrain.i)

# This tests the thermal expansion coefficient function using both
# options to specify that function: mean and instantaneous.  There
# two blocks, each containing a single element, and these use the
# two variants of the function.

# In this test, the instantaneous CTE function has a constant value,
# while the mean CTE function is an analytic function designed to
# give the same response.  If \bar{alpha}(T) is the mean CTE function,
# and \alpha(T) is the instantaneous CTE function,

# \bar{\alpha}(T) = 1/(T-Tref) \intA^{T}_{Tsf} \alpha(T) dT

# where Tref is the reference temperature used to define the mean CTE
# function, and Tsf is the stress-free temperature.

# This version of the test uses finite deformation theory.
# The two models produce very similar results.  There are slight
# differences due to the large deformation treatment.

[Mesh]
  file = 'blocks.e'
[]

[GlobalParams]
  displacements = 'disp_x disp_y disp_z'
[]

[AuxVariables]
  [./temp]
    order = FIRST
    family = LAGRANGE
  [../]
[]

[Problem]
  solve = false
[]

[Modules/TensorMechanics/Master]
  [./block1]
    block = 1
    strain = FINITE
    add_variables = true
    automatic_eigenstrain_names = true
    generate_output = 'strain_xx strain_yy strain_zz'
  [../]
  [./block2]
    block = 2
    strain = FINITE
    add_variables = true
    automatic_eigenstrain_names = true
    generate_output = 'strain_xx strain_yy strain_zz'
  [../]
[]

[BCs]
  [./left]
    type = DirichletBC
    variable = disp_x
    boundary = 3
    value = 0.0
  [../]

  [./bottom]
    type = DirichletBC
    variable = disp_y
    boundary = 2
    value = 0.0
  [../]

  [./back]
    type = DirichletBC
    variable = disp_z
    boundary = 1
    value = 0.0
  [../]
[]

[AuxKernels]
  [./temp]
    type = FunctionAux
    variable = temp
    block = '1 2'
    function = temp_func
  [../]
[]

[Materials]
  [./elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 1e6
    poissons_ratio = 0.3
  [../]
  [./small_stress]
    type = ComputeFiniteStrainElasticStress
  [../]
  [./thermal_expansion_strain1]
    type = ComputeMeanThermalExpansionFunctionEigenstrain
    block = 1
    thermal_expansion_function = cte_func_mean
    thermal_expansion_function_reference_temperature = 0.5
    stress_free_temperature = 0.0
    temperature = temp
    eigenstrain_name = eigenstrain1
  [../]
  [./thermal_expansion_strain2]
    type = ComputeInstantaneousThermalExpansionFunctionEigenstrain
    block = 2
    thermal_expansion_function = cte_func_inst
    stress_free_temperature = 0.0
    temperature = temp
    eigenstrain_name = eigenstrain2
  [../]
[]

[Functions]
  [./cte_func_mean]
    type = ParsedFunction
    vars = 'tsf tref scale' #stress free temp, reference temp, scale factor
    vals = '0.0 0.5  1e-4'
    value = 'scale * (t - tsf) / (t - tref)'
  [../]
  [./cte_func_inst]
    type = PiecewiseLinear
    xy_data = '0 1.0
               2 1.0'
    scale_factor = 1e-4
  [../]

  [./temp_func]
    type = PiecewiseLinear
    xy_data = '0 1
               1 2'
  [../]
[]

[Postprocessors]
  [./disp_1]
    type = NodalExtremeValue
    variable = disp_x
    boundary = 101
  [../]

  [./disp_2]
    type = NodalExtremeValue
    variable = disp_x
    boundary = 102
  [../]
[]

[Executioner]
  type = Transient

  solve_type = PJFNK
  l_max_its = 100
  l_tol = 1e-4
  nl_abs_tol = 1e-8
  nl_rel_tol = 1e-12

  start_time = 0.0
  end_time = 1.0
  dt = 0.1
[]

[Outputs]
  csv = true
[]

(modules/contact/examples/3d_berkovich/indenter_berkovich_friction.i)

[Mesh]
  file = indenter.e
[]

[GlobalParams]
  displacements = 'disp_x disp_y disp_z'
  volumetric_locking_correction = true
  order = FIRST
  family = LAGRANGE
[]

[Variables]
  [./disp_x]
  [../]
  [./disp_y]
  [../]
  [./disp_z]
  [../]
[]

[AuxVariables]
  [./saved_x]
  [../]
  [./saved_y]
  [../]
  [./saved_z]
  [../]
[]

[AuxKernels]
[]

[Functions]
  [./push_down]
    type = ParsedFunction
    value = 'if(t < 1.5, -t, t-3.0)'
  [../]
[]


[Modules/TensorMechanics/Master]
  [./all]
    add_variables = true
    strain = FINITE
    block = '1 2'
    use_automatic_differentiation = false
    generate_output = 'stress_xx stress_xy stress_xz stress_yy stress_zz'
    save_in = 'saved_x saved_y saved_z'
    use_finite_deform_jacobian = true
  [../]
[]

[BCs]
  [./botz]
    type = DirichletBC
    variable = disp_z
    boundary = 101
    value = 0.0
  [../]
  [./boty]
    type = DirichletBC
    variable = disp_y
    boundary = 101
    value = 0.0
  [../]
  [./botx]
    type = DirichletBC
    variable = disp_x
    boundary = 101
    value = 0.0
  [../]

  [./boty111]
    type = DirichletBC
    variable = disp_y
    boundary = 111
    value = 0.0
  [../]
  [./botx111]
    type = DirichletBC
    variable = disp_x
    boundary = 111
    value = 0.0
  [../]

  [./topz]
    type = FunctionDirichletBC
    variable = disp_z
    boundary = '201'
    function = push_down
  [../]
  [./topy]
    type = DirichletBC
    variable = disp_y
    boundary = 201
    value = 0.0
  [../]
  [./topx]
    type = DirichletBC
    variable = disp_x
    boundary = 201
    value = 0.0
  [../]
[]


[UserObjects]
  [./slip_rate_gss]
    type = CrystalPlasticitySlipRateGSS
    variable_size = 48
    slip_sys_file_name = input_slip_sys_bcc48.txt
    num_slip_sys_flowrate_props = 2
    flowprops = '1 48 0.0001 0.01'
    uo_state_var_name = state_var_gss
    slip_incr_tol = 10.0
    block = 1
  [../]
  [./slip_resistance_gss]
    type = CrystalPlasticitySlipResistanceGSS
    variable_size = 48
    uo_state_var_name = state_var_gss
    block = 1
  [../]
  [./state_var_gss]
    type = CrystalPlasticityStateVariable
    variable_size = 48
    groups = '0 24 48'
    group_values = '900 1000' #120
    uo_state_var_evol_rate_comp_name = state_var_evol_rate_comp_gss
    scale_factor = 1.0
    block = 1
  [../]
  [./state_var_evol_rate_comp_gss]
    type = CrystalPlasticityStateVarRateComponentGSS
    variable_size = 48
    hprops = '1.4 1000 1200 2.5'
    uo_slip_rate_name = slip_rate_gss
    uo_state_var_name = state_var_gss
    block = 1
  [../]
[]

[Materials]
  [./crysp]
    type = FiniteStrainUObasedCP
    block = 1
    stol = 1e-2
    tan_mod_type = exact
    uo_slip_rates = 'slip_rate_gss'
    uo_slip_resistances = 'slip_resistance_gss'
    uo_state_vars = 'state_var_gss'
    uo_state_var_evol_rate_comps = 'state_var_evol_rate_comp_gss'
    maximum_substep_iteration = 25
  [../]
  [./elasticity_tensor]
    type = ComputeElasticityTensorCP
    block = 1
    C_ijkl = '265190 113650 113650 265190 113650 265190 75769 75769 75760'
    fill_method = symmetric9
  [../]
  [./elasticity_tensor_indenter]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 1000000.0
    poissons_ratio = 0.3
    block = 2
  [../]
  [./stress_indenter]
    type = ComputeFiniteStrainElasticStress
    block = 2
  [../]
[]

[Postprocessors]
  [./stress_zz]
    type = ElementAverageValue
    variable = stress_zz
    block = 1
  [../]
  [./resid_z]
    type = NodalSum
    variable = saved_z
    boundary = 201
  [../]
  [./disp_z]
    type = NodalExtremeValue
    variable = disp_z
    boundary = 201
  [../]
[]

[Executioner]
  type = Transient
  solve_type = 'PJFNK'

  petsc_options_iname = '-pc_type -pc_factor_mat_solver_type'
  petsc_options_value = 'lu    superlu_dist'
  line_search = 'none'

  l_max_its = 60
  nl_max_its = 50
  dt = 0.004
  dtmin = 0.00001
  end_time = 1.8
  nl_rel_tol = 1e-8
  nl_abs_tol = 1e-6 # 6 if no friction
  l_tol = 1e-3
  automatic_scaling = true

[]

[Outputs]
  [./my_checkpoint]
    type = Checkpoint
    interval = 50
  [../]
  exodus = true
  csv = true
  print_linear_residuals = true
  print_perf_log = true
  [./console]
    type = Console
    max_rows = 5
  [../]
[]

[Preconditioning]
  [./smp]
    type = SMP
    full = true
  [../]
[]

[Dampers]
  [./contact_slip]
    type = ContactSlipDamper
    primary = '202'
    secondary = '102'
  [../]
[]

[Contact]
  [./ind_base]
    primary = 202
    secondary = 102
    model = coulomb
    friction_coefficient = 0.4
    normalize_penalty = true
    formulation = tangential_penalty
    penalty = 1e7
    capture_tolerance = 0.0001
  [../]
[]

(modules/stochastic_tools/test/tests/surrogates/poly_chaos/sub.i)

[Mesh]
  type = GeneratedMesh
  dim = 1
  nx = 100
  xmax = 10
[]

[Variables]
  [u]
  []
[]

[Kernels]
  [diffusion]
    type = MatDiffusion
    variable = u
    diffusivity = D
  []
  [absorption]
    type = MaterialReaction
    variable = u
    coefficient = sig
  []
  [source]
    type = BodyForce
    variable = u
    value = 1.0
  []
[]

[Materials]
  [diffusivity]
    type = GenericConstantMaterial
    prop_names = D
    prop_values = 2.0
  []
  [xs]
    type = GenericConstantMaterial
    prop_names = sig
    prop_values = 2.0
  []
[]

[BCs]
  [left]
    type = DirichletBC
    variable = u
    boundary = left
    value = 0
  []
  [right]
    type = DirichletBC
    variable = u
    boundary = right
    value = 0
  []
[]

[Executioner]
  type = Steady
  solve_type = PJFNK
  petsc_options_iname = '-pc_type -pc_hypre_type'
  petsc_options_value = 'hypre boomeramg'
[]

[Controls]
  [stochastic]
    type = SamplerReceiver
  []
[]

[Postprocessors]
  [avg]
    type = AverageNodalVariableValue
    variable = u
  []
  [max]
    type = NodalExtremeValue
    variable = u
    value_type = max
  []
[]

[Outputs]
[]

(modules/tensor_mechanics/test/tests/plane_stress/ad_weak_plane_stress_small.i)

[GlobalParams]
  order = FIRST
  family = LAGRANGE
  displacements = 'disp_x disp_y'
  temperature = temp
  out_of_plane_strain = strain_zz
[]

[Mesh]
  [./square]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 2
    ny = 2
  [../]
[]

[Variables]
  [./disp_x]
  [../]
  [./disp_y]
  [../]
  [./strain_zz]
  [../]
[]

[AuxVariables]
  [./temp]
  [../]
  [./nl_strain_zz]
    order = CONSTANT
    family = MONOMIAL
  [../]
[]

[Postprocessors]
  [./react_z]
    type = ADMaterialTensorIntegral
    rank_two_tensor = stress
    index_i = 2
    index_j = 2
  [../]
  [./min_strain_zz]
    type = NodalExtremeValue
    variable = strain_zz
    value_type = min
  [../]
  [./max_strain_zz]
    type = NodalExtremeValue
    variable = strain_zz
    value_type = max
  [../]
[]

[Modules/TensorMechanics/Master]
  [./plane_stress]
    planar_formulation = WEAK_PLANE_STRESS
    strain = SMALL
    generate_output = 'stress_xx stress_xy stress_yy stress_zz strain_xx strain_xy strain_yy'
    eigenstrain_names = eigenstrain
    use_automatic_differentiation = true
  [../]
[]

[AuxKernels]
  [./tempfuncaux]
    type = FunctionAux
    variable = temp
    function = tempfunc
    use_displaced_mesh = false
  [../]
  [./strain_zz]
    type = ADRankTwoAux
    rank_two_tensor = total_strain
    variable = nl_strain_zz
    index_i = 2
    index_j = 2
  [../]
[]

[Functions]
  [./pull]
    type = PiecewiseLinear
    x='0     1   100'
    y='0  0.00  0.00'
  [../]
  [./tempfunc]
    type = ParsedFunction
    value = '(1 - x) * t'
  [../]
[]

[BCs]
  [./bottomx]
    type = DirichletBC
    boundary = 0
    variable = disp_x
    value = 0.0
  [../]
  [./bottomy]
    type = DirichletBC
    boundary = 0
    variable = disp_y
    value = 0.0
  [../]
[]

[Materials]
  [./elasticity_tensor]
    type = ADComputeIsotropicElasticityTensor
    poissons_ratio = 0.3
    youngs_modulus = 1e6
  [../]
  [./thermal_strain]
    type = ADComputeThermalExpansionEigenstrain
    thermal_expansion_coeff = 0.02
    stress_free_temperature = 0.5
    eigenstrain_name = eigenstrain
  [../]
  [./stress]
    type = ADComputeLinearElasticStress
  [../]
[]

[Executioner]
  type = Transient

  solve_type = 'NEWTON'

  # controls for nonlinear iterations
  nl_max_its = 15
  nl_rel_tol = 1e-14
  nl_abs_tol = 1e-12

  # time control
  start_time = 0.0
  dt = 1.0
  dtmin = 1.0
  end_time = 2.0
[]

[Outputs]
  file_base = 'weak_plane_stress_small_out'
  exodus = true
[]

(modules/stochastic_tools/examples/surrogates/gaussian_process/sub.i)

[Mesh]
  type = GeneratedMesh
  dim = 1
  nx = 100
  xmax = 0.03
  elem_type = EDGE3
[]

[Variables]
  [T]
    order = SECOND
    family = LAGRANGE
  []
[]

[Kernels]
  [diffusion]
    type = MatDiffusion
    variable = T
    diffusivity = k
  []
  [source]
    type = BodyForce
    variable = T
    value = 10000
  []
[]

[Materials]
  [conductivity]
    type = GenericConstantMaterial
    prop_names = k
    prop_values = 5.0
  []
[]

[BCs]
  [right]
    type = DirichletBC
    variable = T
    boundary = right
    value = 300
  []
[]

[Executioner]
  type = Steady
  solve_type = PJFNK
  petsc_options_iname = '-pc_type -pc_hypre_type'
  petsc_options_value = 'hypre boomeramg'
[]

[Postprocessors]
  [avg]
    type = AverageNodalVariableValue
    variable = T
  []
  [max]
    type = NodalExtremeValue
    variable = T
    value_type = max
  []
[]

[Outputs]
[]

(modules/combined/test/tests/restart-transient-from-ss-with-stateful/sub_ss.i)

[GlobalParams]
  displacements = 'disp_x disp_y'
[]

[Mesh]
  [gen]
    type = GeneratedMeshGenerator
    nx = 8
    ny = 8
    xmin = -82.627
    xmax = 82.627
    ymin = -82.627
    ymax = 82.627
    dim = 2
  []
  [./extra_nodes_x]
    type = ExtraNodesetGenerator
    input = 'gen'
    new_boundary = 'no_x'
    coord = '0 82.627 0'
  [../]
  [./extra_nodes_y]
    type = ExtraNodesetGenerator
    input = 'extra_nodes_x'
    new_boundary = 'no_y'
    coord = '-82.627 0 0'
  [../]
[]

[Problem]
  type = ReferenceResidualProblem
  reference_vector = 'ref'
  extra_tag_vectors = 'ref'
[]

[AuxVariables]
  [./temp]
  [../]
[]

[Modules/TensorMechanics/Master]
  # FINITE strain when strain is large, i.e., visible movement.
  # SMALL strain when things are stressed, but may not move.
  [./fuel]
    add_variables = true
    strain = FINITE
    temperature = temp
    eigenstrain_names = 'thermal_eigenstrain'
    generate_output = 'vonmises_stress stress_xx stress_yy hydrostatic_stress max_principal_stress strain_xy elastic_strain_xx stress_xy'
    extra_vector_tags = 'ref'
    use_finite_deform_jacobian = true
    incremental = true
  [../]
[]

[BCs]
  [./no_x]
    type = DirichletBC
    variable = disp_x
    boundary = 'no_x'
    value = 0.0
    preset = true
  [../]
  [./no_y]
    type = DirichletBC
    preset = true
    variable = disp_y
    boundary = 'no_y'
    value = 0.0
  [../]
[]

[Materials]
  [./elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 3e10   # Pa
    poissons_ratio = 0.33    # unitless
  [../]
  [./thermal_strains]
    type = ComputeThermalExpansionEigenstrain
    temperature = temp
    thermal_expansion_coeff = 2e-6 # 1/K
    stress_free_temperature = 500 # K
    eigenstrain_name = 'thermal_eigenstrain'
  [../]
  [./stress_finite] # goes with FINITE strain formulation
    type = ComputeFiniteStrainElasticStress
  [../]
[]

[Postprocessors]
  [./avg_temp]
    type = ElementAverageValue
    variable = temp
  [../]
  [./disp_x_max_element]
    type = ElementExtremeValue
    value_type = max
    variable = disp_x
    execute_on = 'initial timestep_end'
  [../]
  [./disp_y_max_element]
    type = ElementExtremeValue
    value_type = max
    variable = disp_y
    execute_on = 'initial timestep_end'
  [../]
  [./disp_x_max_nodal]
    type = NodalExtremeValue
    value_type = max
    variable = disp_x
    execute_on = 'initial timestep_end'
  [../]
  [./disp_y_max_nodal]
    type = NodalExtremeValue
    value_type = max
    variable = disp_y
    execute_on = 'initial timestep_end'
  [../]
[]

[Executioner]
  type = Steady
  petsc_options_iname = '-pc_type -pc_hypre_type -ksp_gmres_restart'
  petsc_options_value = 'hypre boomeramg 300'
  line_search = 'none'

  l_tol = 1e-02
  nl_rel_tol = 5e-04
  nl_abs_tol = 1e-2

  l_max_its = 50
  nl_max_its = 25
[]

[Outputs]
  exodus = true
  print_linear_residuals = false
  perf_graph = true
[]

(modules/tensor_mechanics/tutorials/introduction/mech_step04a.i)

#
# We study the effects of volumetric locking
# https://mooseframework.inl.gov/modules/tensor_mechanics/tutorials/introduction/answer04b.html
#

[GlobalParams]
  displacements = 'disp_x disp_y'
  # elem_type applies to the GeneratedMeshGenerator blocks
  elem_type = QUAD4
  # volumetric_locking_correction applies to the TensorMechanics Master action
  volumetric_locking_correction = false
  # uniform_refine applies to the final mesh
  uniform_refine = 0
[]

[Mesh]
  [generated1]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 5
    ny = 15
    xmin = -0.6
    xmax = -0.1
    ymax = 5
    bias_y = 0.9
    boundary_name_prefix = pillar1
  []
  [generated2]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 4
    ny = 15
    xmin = 0.1
    xmax = 0.6
    ymax = 5
    bias_y = 0.9
    boundary_name_prefix = pillar2
    boundary_id_offset = 4
  []
  [collect_meshes]
    type = MeshCollectionGenerator
    inputs = 'generated1 generated2'
  []
[]

[Modules/TensorMechanics/Master]
  [all]
    add_variables = true
    strain = FINITE
  []
[]

[BCs]
  [bottom_x]
    type = DirichletBC
    variable = disp_x
    boundary = 'pillar1_bottom pillar2_bottom'
    value = 0
  []
  [bottom_y]
    type = DirichletBC
    variable = disp_y
    boundary = 'pillar1_bottom pillar2_bottom'
    value = 0
  []
  [Pressure]
    [sides]
      boundary = 'pillar1_left pillar2_right'
      function = 1e4*t
    []
  []
[]

[Materials]
  [elasticity]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 1e9
    # near incopmpressible material
    poissons_ratio = 0.49
  []
  [stress]
    type = ComputeFiniteStrainElasticStress
  []
[]

[Postprocessors]
  [x_deflection]
    type = NodalExtremeValue
    value_type = max
    variable = disp_x
  []
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  line_search = none
  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'
  end_time = 5
  dt = 0.5
  [Predictor]
    type = SimplePredictor
    scale = 1
  []
[]

[Outputs]
  exodus = true
  csv = true
[]

(modules/tensor_mechanics/test/tests/dynamics/time_integration/hht_test_ti.i)

# Test for  HHT time integration

# The test is for an 1D bar element of  unit length fixed on one end
# with a ramped pressure boundary condition applied to the other end.
# alpha, beta and gamma are HHT time integration parameters
# The equation of motion in terms of matrices is:
#
# M*accel + alpha*(K*disp - K*disp_old) + K*disp = P(t+alpha dt)*Area
#
# Here M is the mass matrix, K is the stiffness matrix, P is the applied pressure
#
# This equation is equivalent to:
#
# density*accel + alpha*(Div stress - Div stress_old) +Div Stress= P(t+alpha dt)
#
# The first term on the left is evaluated using the Inertial force kernel
# The next two terms on the left involving alpha are evaluated using the
# DynamicStressDivergenceTensors Kernel
# The residual due to Pressure is evaluated using Pressure boundary condition
#
# The system will come to steady state slowly after the pressure becomes constant.
# Alpha equal to zero will result in Newmark integration.

[Mesh]
  type = GeneratedMesh
  dim = 3
  nx = 1
  ny = 1
  nz = 1
  xmin = 0.0
  xmax = 0.1
  ymin = 0.0
  ymax = 1.0
  zmin = 0.0
  zmax = 0.1
[]


[Variables]
  [./disp_x]
  [../]
  [./disp_y]
  [../]
  [./disp_z]
  [../]
[]

[AuxVariables]
  [./vel_x]
  [../]
  [./accel_x]
  [../]
  [./vel_y]
  [../]
  [./accel_y]
  [../]
  [./vel_z]
  [../]
  [./accel_z]
  [../]
  [./stress_yy]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./strain_yy]
    order = CONSTANT
    family = MONOMIAL
  [../]

[]

[Kernels]
  [./DynamicTensorMechanics]
    displacements = 'disp_x disp_y disp_z'
    hht_alpha = 0.11
  [../]
  [./inertia_x]
    type = InertialForce
    variable = disp_x
  [../]
  [./inertia_y]
    type = InertialForce
    variable = disp_y
  [../]
  [./inertia_z]
    type = InertialForce
    variable = disp_z
  [../]

[]

[AuxKernels]
  [./accel_x] # These auxkernls are only for checking output
    type = TestNewmarkTI
    displacement = disp_x
    variable = accel_x
    first = false
  [../]
  [./accel_y]
    type = TestNewmarkTI
    displacement = disp_y
    variable = accel_y
    first = false
  [../]
  [./accel_z]
    type = TestNewmarkTI
    displacement = disp_z
    variable = accel_z
    first = false
  [../]
  [./vel_x]
    type = TestNewmarkTI
    displacement = disp_x
    variable = vel_x
  [../]
  [./vel_y]
    type = TestNewmarkTI
    displacement = disp_y
    variable = vel_y
  [../]
  [./vel_z]
    type = TestNewmarkTI
    displacement = disp_z
    variable = vel_z
  [../]
  [./stress_yy]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_yy
    index_i = 0
    index_j = 1
  [../]
  [./strain_yy]
    type = RankTwoAux
    rank_two_tensor = total_strain
    variable = strain_yy
    index_i = 0
    index_j = 1
  [../]
[]


[BCs]
  [./top_y]
    type = DirichletBC
    variable = disp_y
    boundary = top
    value=0.0
  [../]
  [./top_x]
    type = DirichletBC
    variable = disp_x
    boundary = top
    value=0.0
  [../]
  [./top_z]
    type = DirichletBC
    variable = disp_z
    boundary = top
    value=0.0
  [../]
  [./bottom_x]
    type = DirichletBC
    variable = disp_x
    boundary = bottom
    value=0.0
  [../]
  [./bottom_z]
    type = DirichletBC
    variable = disp_z
    boundary = bottom
    value=0.0
  [../]
  [./Pressure]
    [./Side1]
      boundary = bottom
      function = pressure
      displacements = 'disp_x disp_y disp_z'
      factor = 1
      alpha = 0.11
    [../]
  [../]
[]

[Materials]
  [./Elasticity_tensor]
    type = ComputeElasticityTensor
    block = 0
    fill_method = symmetric_isotropic
    C_ijkl = '210e9 0'
  [../]

  [./strain]
    type = ComputeSmallStrain
    block = 0
    displacements = 'disp_x disp_y disp_z'
  [../]

  [./stress]
    type = ComputeLinearElasticStress
    block = 0
  [../]
  [./density]
    type = GenericConstantMaterial
    block = 0
    prop_names = 'density'
    prop_values = '7750'
  [../]

[]

[Executioner]
  type = Transient
  start_time = 0
  end_time = 2
  dt = 0.1

  # Time integration scheme
  scheme = 'newmark-beta'
[]


[Functions]
  [./pressure]
    type = PiecewiseLinear
    x = '0.0 0.1 0.2 1.0 2.0 5.0'
    y = '0.0 0.1 0.2 1.0 1.0 1.0'
    scale_factor = 1e9
  [../]
[]

[Postprocessors]
  [./_dt]
    type = TimestepSize
  [../]
  [./disp]
    type = NodalExtremeValue
    variable = disp_y
    boundary = bottom
  [../]
  [./vel]
    type = NodalExtremeValue
    variable = vel_y
    boundary = bottom
  [../]
  [./accel]
    type = NodalExtremeValue
    variable = accel_y
    boundary = bottom
  [../]
  [./stress_yy]
    type = ElementAverageValue
    variable = stress_yy
  [../]
  [./strain_yy]
    type = ElementAverageValue
    variable = strain_yy
  [../]
[]

[Outputs]
  file_base = 'hht_test_out'
  exodus = true
  perf_graph = true
[]

Example Input File Syntax
Input Parameters
Input Files

Install MOOSE

New Users

Examples and Tutorials

Application Usage

Physics and Syntax

Application Development

Framework Development

MOOSEDocs

Infrastructure

Questions

Information and Tools

INL Applications and Remote Access