FunctionDiffusion

Diffusion with a function coefficient.

Overview

The FunctionDiffusion object represents a diffusion term with a function coefficient for scalar field variables. This term is

$var element = document.getElementById("moose-equation-dddeb565-3eae-45e3-b051-5de8f6dcf098");katex.render(" - \\nabla \\cdot f(\\mathbf{r}, t) \\nabla v", element, {displayMode:true,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$

where

$var element = document.getElementById("moose-equation-d2c113e7-d086-4a1c-b244-aadf5c4d9ed4");katex.render("v", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ is a scalar field solution variable, and
$var element = document.getElementById("moose-equation-232d29a5-725c-4240-a11d-0fe0660a2f7f");katex.render("f(\\mathbf{r}, t)", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ is a scalar function serving as a diffusion coefficient.

Note that the "v" parameter is optional. If no variable is entered there, then the kernel's nonlinear variable will be used in the operator as usual.

Example Input File Syntax

[Kernels]
  [diff]
    type = FunctionDiffusion
    variable = u
    function = 'x'
  []
[]

Input Parameters

variableThe name of the variable that this residual object operates on
C++ Type:NonlinearVariableName
Unit:(no unit assumed)
Controllable:No
Description:The name of the variable that this residual object operates on

Required Parameters

blockThe list of blocks (ids or names) that this object will be applied
C++ Type:std::vector<SubdomainName>
Unit:(no unit assumed)
Controllable:No
Description:The list of blocks (ids or names) that this object will be applied
displacementsThe displacements
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:The displacements
function1Function multiplier for diffusion term.
Default:1
C++ Type:FunctionName
Unit:(no unit assumed)
Controllable:No
Description:Function multiplier for diffusion term.
prop_getter_suffixAn optional suffix parameter that can be appended to any attempt to retrieve/get material properties. The suffix will be prepended with a '_' character.
C++ Type:MaterialPropertyName
Unit:(no unit assumed)
Controllable:No
Description:An optional suffix parameter that can be appended to any attempt to retrieve/get material properties. The suffix will be prepended with a '_' character.
use_interpolated_stateFalseFor the old and older state use projected material properties interpolated at the quadrature points. To set up projection use the ProjectedStatefulMaterialStorageAction.
Default:False
C++ Type:bool
Unit:(no unit assumed)
Controllable:No
Description:For the old and older state use projected material properties interpolated at the quadrature points. To set up projection use the ProjectedStatefulMaterialStorageAction.
vCoupled concentration variable for kernel to operate on; if this is not specified, the kernel's nonlinear variable will be used as usual
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:Coupled concentration variable for kernel to operate on; if this is not specified, the kernel's nonlinear variable will be used as usual

Optional Parameters

absolute_value_vector_tagsThe tags for the vectors this residual object should fill with the absolute value of the residual contribution
C++ Type:std::vector<TagName>
Unit:(no unit assumed)
Controllable:No
Description:The tags for the vectors this residual object should fill with the absolute value of the residual contribution
extra_matrix_tagsThe extra tags for the matrices this Kernel should fill
C++ Type:std::vector<TagName>
Unit:(no unit assumed)
Controllable:No
Description:The extra tags for the matrices this Kernel should fill
extra_vector_tagsThe extra tags for the vectors this Kernel should fill
C++ Type:std::vector<TagName>
Unit:(no unit assumed)
Controllable:No
Description:The extra tags for the vectors this Kernel should fill
matrix_tagssystemThe tag for the matrices this Kernel should fill
Default:system
C++ Type:MultiMooseEnum
Unit:(no unit assumed)
Options:nontime, system
Controllable:No
Description:The tag for the matrices this Kernel should fill
vector_tagsnontimeThe tag for the vectors this Kernel should fill
Default:nontime
C++ Type:MultiMooseEnum
Unit:(no unit assumed)
Options:nontime, time
Controllable:No
Description:The tag for the vectors this Kernel should fill

Tagging Parameters

control_tagsAdds user-defined labels for accessing object parameters via control logic.
C++ Type:std::vector<std::string>
Unit:(no unit assumed)
Controllable:No
Description:Adds user-defined labels for accessing object parameters via control logic.
diag_save_inThe name of auxiliary variables to save this Kernel's diagonal Jacobian contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)
C++ Type:std::vector<AuxVariableName>
Unit:(no unit assumed)
Controllable:No
Description:The name of auxiliary variables to save this Kernel's diagonal Jacobian contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)
enableTrueSet the enabled status of the MooseObject.
Default:True
C++ Type:bool
Unit:(no unit assumed)
Controllable:Yes
Description:Set the enabled status of the MooseObject.
implicitTrueDetermines whether this object is calculated using an implicit or explicit form
Default:True
C++ Type:bool
Unit:(no unit assumed)
Controllable:No
Description:Determines whether this object is calculated using an implicit or explicit form
save_inThe name of auxiliary variables to save this Kernel's residual contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)
C++ Type:std::vector<AuxVariableName>
Unit:(no unit assumed)
Controllable:No
Description:The name of auxiliary variables to save this Kernel's residual contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)
seed0The seed for the master random number generator
Default:0
C++ Type:unsigned int
Unit:(no unit assumed)
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
Unit:(no unit assumed)
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/combined/test/tests/optimization/optimization_density_update/top_opt_2d_pde_filter.i)
(modules/combined/test/tests/optimization/compliance_sensitivity/3d_mbb.i)
(test/tests/kernels/function_diffusion/function_diffusion.i)
(modules/combined/test/tests/optimization/optimization_density_update/top_opt_3d_pde_filter.i)
(modules/combined/test/tests/optimization/compliance_sensitivity/2d_mbb_pde_amr.i)
(modules/optimization/test/tests/executioners/constrained/shape_optimization/forward.i)
(modules/combined/examples/optimization/2d_mbb_pde.i)
(modules/electromagnetics/test/tests/kernels/scalar_complex_helmholtz/scalar_complex_helmholtz.i)
(modules/combined/examples/optimization/helmholtz_multimat_nostrip.i)
(modules/combined/test/tests/optimization/compliance_sensitivity/2d_mbb_pde.i)
(test/tests/preconditioners/reuse/convergence.i)
(modules/combined/examples/optimization/helmholtz_multimat_strip.i)
(modules/combined/examples/optimization/2d_mbb_pde_amr.i)
(modules/combined/examples/optimization/3d_mbb.i)

(test/tests/kernels/function_diffusion/function_diffusion.i)

[Mesh]
  type = GeneratedMesh
  dim = 2
  nx = 10
  ny = 10
[]

[Variables]
  [u]
  []
[]

[Kernels]
  [diff]
    type = FunctionDiffusion
    variable = u
    function = 'x'
  []
[]

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

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

[Outputs]
  exodus = true
[]

(modules/combined/test/tests/optimization/optimization_density_update/top_opt_2d_pde_filter.i)

vol_frac = 0.4
E0 = 1e5
Emin = 1e-4
power = 2
[GlobalParams]
  displacements = 'disp_x disp_y'
[]

[Mesh]
  [MeshGenerator]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 40
    ny = 20
    xmin = 0
    xmax = 20
    ymin = 0
    ymax = 10
  []
  [node]
    type = ExtraNodesetGenerator
    input = MeshGenerator
    new_boundary = pull
    nodes = 0
  []
[]

[Variables]
  [Dc]
    initial_condition = -1.0
  []
[]

[AuxVariables]
  [sensitivity]
    family = MONOMIAL
    order = FIRST
    initial_condition = -1.0
    [AuxKernel]
      type = MaterialRealAux
      variable = sensitivity
      property = sensitivity
      execute_on = LINEAR
    []
  []

  [compliance]
    family = MONOMIAL
    order = CONSTANT
  []

  [mat_den]
    family = MONOMIAL
    order = CONSTANT
    initial_condition = ${vol_frac}
  []
  [Dc_elem]
    family = MONOMIAL
    order = CONSTANT
    initial_condition = -1.0
    [AuxKernel]
      type = SelfAux
      variable = Dc_elem
      v = Dc
      execute_on = 'TIMESTEP_END'
    []
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [all]
    strain = SMALL
    add_variables = true
    incremental = false
  []
[]

[Kernels]
  [diffusion]
    type = FunctionDiffusion
    variable = Dc
    function = 0.05
  []
  [potential]
    type = Reaction
    variable = Dc
  []
  [source]
    type = CoupledForce
    variable = Dc
    v = sensitivity
  []
[]

[BCs]
  [no_x]
    type = DirichletBC
    variable = disp_x
    boundary = right
    value = 0.0
  []
  [no_y]
    type = DirichletBC
    variable = disp_y
    boundary = right
    value = 0.0
  []
  [boundary_penalty]
    type = ADRobinBC
    variable = Dc
    boundary = 'left top'
    coefficient = 10
  []
[]
[NodalKernels]
  [pull]
    type = NodalGravity
    variable = disp_y
    boundary = pull
    gravity_value = -1
    mass = 1
  []
[]
[Materials]
  [elasticity_tensor]
    type = ComputeVariableIsotropicElasticityTensor
    youngs_modulus = E_phys
    poissons_ratio = poissons_ratio
    args = 'mat_den'
  []
  [E_phys]
    type = DerivativeParsedMaterial
    # Emin + (density^penal) * (E0 - Emin)
    expression = '${Emin} + (mat_den ^ ${power}) * (${E0}-${Emin})'
    coupled_variables = 'mat_den'
    property_name = E_phys
  []
  [poissons_ratio]
    type = GenericConstantMaterial
    prop_names = poissons_ratio
    prop_values = 0.3
  []
  [stress]
    type = ComputeLinearElasticStress
  []
  [dc]
    type = ComplianceSensitivity
    design_density = mat_den
    youngs_modulus = E_phys
    incremental = false
  []
[]

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

[UserObjects]
  [update]
    type = DensityUpdate
    density_sensitivity = Dc_elem
    design_density = mat_den
    volume_fraction = ${vol_frac}
    execute_on = TIMESTEP_BEGIN
    force_postaux = true
  []
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  petsc_options_iname = '-pc_type '
  petsc_options_value = 'lu'
  nl_abs_tol = 1e-10
  line_search = none
  dt = 1.0
  num_steps = 30
[]

[Outputs]
  [out]
    type = Exodus
    time_step_interval = 10
  []
[]

(modules/combined/test/tests/optimization/compliance_sensitivity/3d_mbb.i)

vol_frac = 0.5
E0 = 1
Emin = 1e-8
power = 3
[GlobalParams]
  displacements = 'disp_x disp_y disp_z'
[]

[Mesh]
  [MeshGenerator]
    type = GeneratedMeshGenerator
    dim = 3
    nx = 30
    ny = 10
    nz = 10
    xmin = 0
    xmax = 30
    ymin = 0
    ymax = 10
    zmin = 0
    zmax = 10
  []
  [node]
    type = ExtraNodesetGenerator
    input = MeshGenerator
    new_boundary = hold_y
    coord = '0 0 0; 0 0 10'
  []
  [push]
    type = ExtraNodesetGenerator
    input = node
    new_boundary = push
    coord = '30 10 5'
  []
[]

[Variables]
  [disp_x]
  []
  [disp_y]
  []
  [disp_z]
  []
  [Dc]
    initial_condition = -1.0
  []
[]

[AuxVariables]
  [sensitivity]
    family = MONOMIAL
    order = FIRST
    initial_condition = -1.0
    [AuxKernel]
      type = MaterialRealAux
      variable = sensitivity
      property = sensitivity
      execute_on = LINEAR
    []
  []
  [mat_den]
    family = MONOMIAL
    order = CONSTANT
    initial_condition = ${vol_frac}
  []
  [Dc_elem]
    family = MONOMIAL
    order = CONSTANT
    initial_condition = -1.0
    [AuxKernel]
      type = SelfAux
      variable = Dc_elem
      v = Dc
      execute_on = 'TIMESTEP_END'
    []
  []
  [mat_den_nodal]
    family = L2_LAGRANGE
    order = FIRST
    initial_condition = ${vol_frac}
    [AuxKernel]
      type = SelfAux
      execute_on = TIMESTEP_END
      variable = mat_den_nodal
      v = mat_den
    []
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [all]
    strain = SMALL
    add_variables = true
    incremental = false
  []
[]
[Kernels]
  [diffusion]
    type = FunctionDiffusion
    variable = Dc
    function = 0.15 # radius coeff
  []
  [potential]
    type = Reaction
    variable = Dc
  []
  [source]
    type = CoupledForce
    variable = Dc
    v = sensitivity
  []
[]
[BCs]
  [no_x]
    type = DirichletBC
    variable = disp_y
    boundary = hold_y
    value = 0.0
  []
  [no_y]
    type = DirichletBC
    variable = disp_x
    boundary = right
    value = 0.0
  []
  [boundary_penalty]
    type = ADRobinBC
    variable = Dc
    boundary = 'left top front back'
    coefficient = 10
  []
  [boundary_penalty_right]
    type = ADRobinBC
    variable = Dc
    boundary = 'right'
    coefficient = 10
  []
[]
[NodalKernels]
  [push]
    type = NodalGravity
    variable = disp_y
    boundary = push
    gravity_value = -1
    mass = 1
  []
[]

[Materials]
  [elasticity_tensor]
    type = ComputeVariableIsotropicElasticityTensor
    youngs_modulus = E_phys
    poissons_ratio = poissons_ratio
    args = 'mat_den'
  []
  [E_phys]
    type = DerivativeParsedMaterial
    # Emin + (density^penal) * (E0 - Emin)
    expression = '${Emin} + (mat_den ^ ${power}) * (${E0}-${Emin})'
    coupled_variables = 'mat_den'
    property_name = E_phys
  []
  [poissons_ratio]
    type = GenericConstantMaterial
    prop_names = poissons_ratio
    prop_values = 0.3
  []
  [stress]
    type = ComputeLinearElasticStress
  []
  [dc]
    type = ComplianceSensitivity
    design_density = mat_den
    youngs_modulus = E_phys
    incremental = false
  []
[]

[Preconditioning]
  [smp]
    type = SMP
    full = true
  []
[]
[UserObjects]
  [update]
    type = DensityUpdate
    density_sensitivity = Dc_elem
    design_density = mat_den
    volume_fraction = ${vol_frac}
    execute_on = TIMESTEP_BEGIN
    force_postaux = true
  []
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
  petsc_options_value = 'lu superlu_dist'
  line_search = none
  nl_abs_tol = 1e-4
  l_max_its = 200
  start_time = 0.0
  dt = 1.0
  num_steps = 2
[]

[Outputs]
  [out]
    type = CSV
    execute_on = 'INITIAL TIMESTEP_END'
  []
  print_linear_residuals = false
[]

[Postprocessors]
  [total_vol]
    type = ElementIntegralVariablePostprocessor
    variable = mat_den
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [sensitivity]
    type = ElementIntegralMaterialProperty
    mat_prop = sensitivity
  []
[]

[Controls]
  [first_period]
    type = TimePeriod
    start_time = 0.0
    end_time = 10
    enable_objects = 'BCs::boundary_penalty_right'
    execute_on = 'initial timestep_begin'
  []
[]

(test/tests/kernels/function_diffusion/function_diffusion.i)

[Mesh]
  type = GeneratedMesh
  dim = 2
  nx = 10
  ny = 10
[]

[Variables]
  [u]
  []
[]

[Kernels]
  [diff]
    type = FunctionDiffusion
    variable = u
    function = 'x'
  []
[]

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

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

[Outputs]
  exodus = true
[]

(modules/combined/test/tests/optimization/optimization_density_update/top_opt_3d_pde_filter.i)

vol_frac = 0.4
E0 = 1e5
Emin = 1e-4
power = 2
[GlobalParams]
  displacements = 'disp_x disp_y disp_z'
[]

[Mesh]
  [MeshGenerator]
    type = GeneratedMeshGenerator
    dim = 3
    nx = 24
    ny = 12
    nz = 12
    xmin = 0
    xmax = 20
    ymin = 0
    ymax = 10
    zmin = 0
    zmax = 10
  []
  [middle_bottom_left_edge]
    type = ExtraNodesetGenerator
    input = MeshGenerator
    new_boundary = pull
    coord = '0 0 5'
  []
[]

[Variables]
  [Dc]
    initial_condition = -1.0
  []
[]

[AuxVariables]
  [sensitivity]
    family = MONOMIAL
    order = FIRST
    initial_condition = -1.0
    [AuxKernel]
      type = MaterialRealAux
      variable = sensitivity
      property = sensitivity
      execute_on = LINEAR
    []
  []

  [compliance]
    family = MONOMIAL
    order = CONSTANT
  []

  [mat_den]
    family = MONOMIAL
    order = CONSTANT
    initial_condition = ${vol_frac}
  []
  [Dc_elem]
    family = MONOMIAL
    order = CONSTANT
    initial_condition = -1.0
    [AuxKernel]
      type = SelfAux
      variable = Dc_elem
      v = Dc
      execute_on = 'TIMESTEP_END'
    []
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [all]
    strain = SMALL
    add_variables = true
    incremental = false
  []
[]

[Kernels]
  [diffusion]
    type = FunctionDiffusion
    variable = Dc
    function = 0.05
  []
  [potential]
    type = Reaction
    variable = Dc
  []
  [source]
    type = CoupledForce
    variable = Dc
    v = sensitivity
  []
[]

[BCs]
  [no_x]
    type = DirichletBC
    variable = disp_x
    boundary = right
    value = 0.0
  []
  [no_y]
    type = DirichletBC
    variable = disp_y
    boundary = right
    value = 0.0
  []
  [no_z]
    type = DirichletBC
    variable = disp_z
    boundary = right
    value = 0.0
  []
  [boundary_penalty]
    type = ADRobinBC
    variable = Dc
    boundary = 'left top front back'
    coefficient = 10
  []
[]
[NodalKernels]
  [pull]
    type = NodalGravity
    variable = disp_y
    boundary = pull
    gravity_value = -1
    mass = 1
  []
[]
[Materials]
  [elasticity_tensor]
    type = ComputeVariableIsotropicElasticityTensor
    youngs_modulus = E_phys
    poissons_ratio = poissons_ratio
    args = 'mat_den'
  []
  [E_phys]
    type = DerivativeParsedMaterial
    # Emin + (density^penal) * (E0 - Emin)
    expression = '${Emin} + (mat_den ^ ${power}) * (${E0}-${Emin})'
    coupled_variables = 'mat_den'
    property_name = E_phys
  []
  [poissons_ratio]
    type = GenericConstantMaterial
    prop_names = poissons_ratio
    prop_values = 0.3
  []
  [stress]
    type = ComputeLinearElasticStress
  []
  [dc]
    type = ComplianceSensitivity
    design_density = mat_den
    youngs_modulus = E_phys
    incremental = false
  []
[]

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

[UserObjects]
  [update]
    type = DensityUpdate
    density_sensitivity = Dc_elem
    design_density = mat_den
    volume_fraction = ${vol_frac}
    execute_on = TIMESTEP_BEGIN
  []
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  petsc_options_iname = '-pc_type '
  petsc_options_value = 'lu'
  nl_abs_tol = 1e-10
  line_search = none
  dt = 1.0
  num_steps = 10
[]

[Outputs]
  [out]
    type = Exodus
    time_step_interval = 10
  []
[]

(modules/combined/test/tests/optimization/compliance_sensitivity/2d_mbb_pde_amr.i)

vol_frac = 0.5
E0 = 1
Emin = 1e-8
power = 3
[GlobalParams]
  displacements = 'disp_x disp_y'
[]

[Mesh]
  [MeshGenerator]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 30
    ny = 10
    xmin = 0
    xmax = 30
    ymin = 0
    ymax = 10
  []
  [node]
    type = ExtraNodesetGenerator
    input = MeshGenerator
    new_boundary = hold
    nodes = 0
  []
  [push]
    type = ExtraNodesetGenerator
    input = node
    new_boundary = push
    coord = '30 10 0'
  []

[]

[Variables]
  [disp_x]
  []
  [disp_y]
  []
  [Dc]
    initial_condition = -1.0
  []
[]

[AuxVariables]

  [sensitivity]
    family = MONOMIAL
    order = FIRST
    initial_condition = -1.0
    [AuxKernel]
      type = MaterialRealAux
      variable = sensitivity
      property = sensitivity
      execute_on = LINEAR
    []
  []
  [mat_den]
    family = MONOMIAL
    order = CONSTANT
    initial_condition = ${vol_frac}
  []
  [Dc_elem]
    family = MONOMIAL
    order = CONSTANT
    initial_condition = -1.0
    [AuxKernel]
      type = SelfAux
      variable = Dc_elem
      v = Dc
      execute_on = 'TIMESTEP_END'
    []
  []
  [mat_den_nodal]
    family = L2_LAGRANGE
    order = FIRST
    initial_condition = ${vol_frac}
    [AuxKernel]
      type = SelfAux
      execute_on = TIMESTEP_END
      variable = mat_den_nodal
      v = mat_den
    []
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [all]
    strain = SMALL
    add_variables = true
    incremental = false
  []
[]
[Kernels]
  [diffusion]
    type = FunctionDiffusion
    variable = Dc
    function = 0.15 # radius coeff
  []
  [potential]
    type = Reaction
    variable = Dc
  []
  [source]
    type = CoupledForce
    variable = Dc
    v = sensitivity
  []
[]
[BCs]
  [no_x]
    type = DirichletBC
    variable = disp_y
    boundary = hold
    value = 0.0
  []
  [no_y]
    type = DirichletBC
    variable = disp_x
    boundary = right
    value = 0.0
  []
  [boundary_penalty]
    type = ADRobinBC
    variable = Dc
    boundary = 'left top'
    coefficient = 10
  []
  [boundary_penalty_right]
    type = ADRobinBC
    variable = Dc
    boundary = 'right'
    coefficient = 10
  []
[]
[NodalKernels]
  [push]
    type = NodalGravity
    variable = disp_y
    boundary = push
    gravity_value = -1
    mass = 1
  []
[]

[Materials]
  [elasticity_tensor]
    type = ComputeVariableIsotropicElasticityTensor
    youngs_modulus = E_phys
    poissons_ratio = poissons_ratio
    args = 'mat_den'
  []
  [E_phys]
    type = DerivativeParsedMaterial
    # Emin + (density^penal) * (E0 - Emin)
    expression = '${Emin} + (mat_den ^ ${power}) * (${E0}-${Emin})'
    coupled_variables = 'mat_den'
    property_name = E_phys
  []
  [poissons_ratio]
    type = GenericConstantMaterial
    prop_names = poissons_ratio
    prop_values = 0.3
  []
  [stress]
    type = ComputeLinearElasticStress
  []
  [dc]
    type = ComplianceSensitivity
    design_density = mat_den
    youngs_modulus = E_phys
    incremental = false
  []
[]

[Preconditioning]
  [smp]
    type = SMP
    full = true
  []
[]
[UserObjects]
  [update]
    type = DensityUpdate
    density_sensitivity = Dc_elem
    design_density = mat_den
    volume_fraction = ${vol_frac}
    execute_on = TIMESTEP_BEGIN
    force_postaux = true
  []
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
  petsc_options_value = 'lu superlu_dist'
  line_search = none
  nl_abs_tol = 1e-4
  l_max_its = 200
  start_time = 0.0
  dt = 1.0
  num_steps = 70
[]

[Outputs]
  [out]
    type = CSV
    execute_on = 'INITIAL TIMESTEP_END'
  []
  print_linear_residuals = false
[]

[Postprocessors]
  [total_vol]
    type = ElementIntegralVariablePostprocessor
    variable = mat_den
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [sensitivity]
    type = ElementIntegralMaterialProperty
    mat_prop = sensitivity
  []
[]

[Controls]
  [first_period]
    type = TimePeriod
    start_time = 0.0
    end_time = 40
    enable_objects = 'BCs::boundary_penalty_right'
    execute_on = 'initial timestep_begin'
  []
[]

[Adaptivity]
  max_h_level = 2
  recompute_markers_during_cycles = true
  interval = 1
  cycles_per_step = 1
  marker = density_marker
  [Indicators]
    [density_jump]
      type = ValueJumpIndicator
      variable = mat_den_nodal
    []
  []
  [Markers]
    [density_marker]
      type = ErrorToleranceMarker
      indicator = density_jump
      coarsen = 0.1
      refine = 0.1
    []
  []
[]

(modules/optimization/test/tests/executioners/constrained/shape_optimization/forward.i)

# This test is documented as an example for ConstrainedShapeOptimization. This
# test should not be changed without updating the documentation.

inner_radius = 6
outer_radius = 10
volume_constraint = 200

[GlobalParams]
  displacements = 'disp_x disp_y'
[]

[Mesh]
  [mesh]
    type = ConcentricCircleMeshGenerator
    has_outer_square = no
    num_sectors = 16
    radii = '${inner_radius} ${outer_radius}'
    rings = '16 16'
    preserve_volumes = false
  []
  [inner_radius]
    type = SideSetsBetweenSubdomainsGenerator
    input = mesh
    new_boundary = inner
    primary_block = 2
    paired_block = 1
  []

  [delete]
    type = BlockDeletionGenerator
    input = inner_radius
    block = 1
  []
  [gather_all]
    type = BoundingBoxNodeSetGenerator
    input = delete
    bottom_left = '-100 -100 -100'
    top_right = '100 100 100'
    new_boundary = total
  []
  [combine]
    type = SideSetsFromBoundingBoxGenerator
    input = gather_all
    bottom_left = '-100 -100 -100'
    top_right = '100 100 100'
    boundaries_old = 'inner outer'
    boundary_new = moving
  []

[]
[Variables]
  [disp_x]
  []
  [disp_y]
  []
  [T]
  []
[]

[AuxVariables]
  [dist_between]
    [AuxKernel]
      type = NearestNodeDistanceAux
      variable = dist_between
      paired_boundary = moving
      boundary = total
      block = 2
      use_displaced_mesh = false
      execute_on = "INITIAL"
    []
  []
[]

[Kernels]
  [disp_x]
    type = MatDiffusion
    variable = disp_x
    use_displaced_mesh = false
    diffusivity = diff_coef
  []
  [disp_y]
    type = MatDiffusion
    variable = disp_y
    use_displaced_mesh = false
    diffusivity = diff_coef
  []
  # run physics of interest on deformed mesh
  [Diffusion]
    type = FunctionDiffusion
    variable = T
    use_displaced_mesh = true
  []
  [Source]
    type = BodyForce
    variable = T
    value = 1
    use_displaced_mesh = true
  []
[]

[Materials]
  # perserve elements near the boundary
  [diff_coef]
    type = ParsedMaterial
    coupled_variables = 'dist_between'
    expression = '1/(dist_between+0.5)'
    property_name = 'diff_coef'
  []
  [h]
    type = ADGenericFunctionMaterial
    prop_names = h
    prop_values = h
  []
  # convection type boundary condition
  [convection_bc]
    type = ADParsedMaterial
    coupled_variables = "T"
    expression = "h*(100-T)"
    material_property_names = "h"
    property_name = convection
  []
[]

[Functions]
  [r1_x]
    type = ParsedOptimizationFunction
    expression = 'r1 * cos((atan(y/x)))'
    param_symbol_names = 'r0 r1'
    param_vector_name = 'params/radii'
  []
  [r1_y]
    type = ParsedOptimizationFunction
    expression = 'r1 * sin((atan(y/x)))'
    param_symbol_names = 'r0 r1'
    param_vector_name = 'params/radii'
  []
  [r0_x]
    type = ParsedOptimizationFunction
    expression = 'r0 * cos((atan(y/x)))'
    param_symbol_names = 'r0 r1'
    param_vector_name = 'params/radii'
  []
  [r0_y]
    type = ParsedOptimizationFunction
    expression = 'r0 * sin((atan(y/x)))'
    param_symbol_names = 'r0 r1'
    param_vector_name = 'params/radii'
  []
  [h]
    type = ParsedOptimizationFunction
    # r0+${inner_radius} is the true current inner radius
    expression = '10 /(pi * (r0+${inner_radius})^3)'
    param_symbol_names = 'r0 r1'
    param_vector_name = 'params/radii'
  []
  [eq_grad_r0]
    type = ParsedOptimizationFunction
    expression = '-2 * pi * (r0 + ${inner_radius})'
    param_symbol_names = 'r0 r1'
    param_vector_name = 'params/radii'
  []
  [eq_grad_r1]
    type = ParsedOptimizationFunction
    # r1+${outer_radius} is the true current outer radius
    expression = '2 * pi * (r1+${outer_radius})'
    param_symbol_names = 'r0 r1'
    param_vector_name = 'params/radii'
  []
[]

[BCs]
  [diffuse_r1_x]
    type = ADFunctionDirichletBC
    variable = disp_x
    boundary = 'outer'
    function = r1_x
    preset = false
  []
  [diffuse_r1_y]
    type = ADFunctionDirichletBC
    variable = disp_y
    boundary = 'outer'
    function = r1_y
    preset = false
  []
  [diffuse_r0_x]
    type = ADFunctionDirichletBC
    variable = disp_x
    boundary = 'inner'
    function = r0_x
    preset = false
  []
  [diffuse_r0_y]
    type = ADFunctionDirichletBC
    variable = disp_y
    boundary = 'inner'
    function = r0_y
    preset = false
  []
  # run physics on deformed mesh
  [convection]
    type = ADMatNeumannBC
    variable = T
    boundary = inner
    boundary_material = convection
    use_displaced_mesh = true
    value = 1
  []
[]

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

[Reporters]
  [params]
    type = ConstantReporter
    real_vector_names = 'radii'
    real_vector_values = '0 0'
    dof_id_type_vector_names = 'num_params'
    dof_id_type_vector_values = '2'
  []
[]

[Postprocessors]
  [current_volume]
    type = VolumePostprocessor
    use_displaced_mesh = true
  []
  # objective function
  [objective]
    type = NodalExtremeValue
    variable = T
  []
  [eq_constraint]
    type = ParsedPostprocessor
    pp_names = current_volume
    expression = 'current_volume - ${volume_constraint}'
  []
  [func_r0]
    type = FunctionValuePostprocessor
    function = eq_grad_r0
  []
  [func_r1]
    type = FunctionValuePostprocessor
    function = eq_grad_r1
  []
[]

[VectorPostprocessors]
  # convert "Real" postprocessors to vectors
  [vol_constraint]
    type = VectorOfPostprocessors
    postprocessors = 'eq_constraint'
    force_postaux = true
  []
  [eq_grad]
    type = VectorOfPostprocessors
    postprocessors = 'func_r0 func_r1'
    force_postaux = true
  []
[]

[Executioner]
  type = Steady
  solve_type = 'NEWTON'
  petsc_options_iname = '-pc_type -pc_factor_shift_type'
  petsc_options_value = 'lu       NONZERO'
  line_search = none
  nl_rel_tol = 1e-8
  nl_abs_tol = 1e-8
[]

[Outputs]
  console = false
[]

(modules/combined/examples/optimization/2d_mbb_pde.i)

vol_frac = 0.5
E0 = 1
Emin = 1e-8
power = 2
[GlobalParams]
  displacements = 'disp_x disp_y'
[]

[Mesh]
  [MeshGenerator]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 150
    ny = 50
    xmin = 0
    xmax = 30
    ymin = 0
    ymax = 10
  []
  [node]
    type = ExtraNodesetGenerator
    input = MeshGenerator
    new_boundary = hold_y
    nodes = 0
  []
  [push]
    type = ExtraNodesetGenerator
    input = node
    new_boundary = push
    coord = '30 10 0'
  []

[]

[Variables]
  [Dc]
    initial_condition = -1.0
  []
[]

[AuxVariables]
  [Emin]
    family = MONOMIAL
    order = CONSTANT
    initial_condition = ${Emin}
  []
  [power]
    family = MONOMIAL
    order = CONSTANT
    initial_condition = ${power}
  []
  [E0]
    family = MONOMIAL
    order = CONSTANT
    initial_condition = ${E0}
  []
  [sensitivity]
    family = MONOMIAL
    order = FIRST
    initial_condition = -1.0
    [AuxKernel]
      type = MaterialRealAux
      variable = sensitivity
      property = sensitivity
      execute_on = LINEAR
    []
  []
  [mat_den]
    family = MONOMIAL
    order = CONSTANT
    initial_condition = ${vol_frac}
  []
  [Dc_elem]
    family = MONOMIAL
    order = CONSTANT
    initial_condition = -1.0
    [AuxKernel]
      type = SelfAux
      variable = Dc_elem
      v = Dc
      execute_on = 'TIMESTEP_END'
    []
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [all]
    strain = SMALL
    add_variables = true
    incremental = false
  []
[]
[Kernels]
  [diffusion]
    type = FunctionDiffusion
    variable = Dc
    function = 0.15 # radius coeff
  []
  [potential]
    type = Reaction
    variable = Dc
  []
  [source]
    type = CoupledForce
    variable = Dc
    v = sensitivity
  []
[]
[BCs]
  [no_x]
    type = DirichletBC
    variable = disp_y
    boundary = hold_y
    value = 0.0
  []
  [no_y]
    type = DirichletBC
    variable = disp_x
    boundary = right
    value = 0.0
  []
  [boundary_penalty]
    type = ADRobinBC
    variable = Dc
    boundary = 'left top'
    coefficient = 10
  []
  [boundary_penalty_right]
    type = ADRobinBC
    variable = Dc
    boundary = 'right'
    coefficient = 10
  []
[]
[NodalKernels]
  [push]
    type = NodalGravity
    variable = disp_y
    boundary = push
    gravity_value = -1
    mass = 1
  []
[]

[Materials]
  [elasticity_tensor]
    type = ComputeVariableIsotropicElasticityTensor
    youngs_modulus = E_phys
    poissons_ratio = poissons_ratio
    args = 'Emin mat_den power E0'
  []
  [E_phys]
    type = DerivativeParsedMaterial
    # Emin + (density^penal) * (E0 - Emin)
    expression = '${Emin} + (mat_den ^ ${power}) * (${E0}-${Emin})'
    coupled_variables = 'mat_den'
    property_name = E_phys
  []
  [poissons_ratio]
    type = GenericConstantMaterial
    prop_names = poissons_ratio
    prop_values = 0.3
  []
  [stress]
    type = ComputeLinearElasticStress
  []
  [dc]
    type = ComplianceSensitivity
    design_density = mat_den
    youngs_modulus = E_phys
    incremental = false
  []
[]

[Preconditioning]
  [smp]
    type = SMP
    full = true
  []
[]
[UserObjects]
  [update]
    type = DensityUpdate
    density_sensitivity = Dc_elem
    design_density = mat_den
    volume_fraction = ${vol_frac}
    execute_on = TIMESTEP_BEGIN
    force_postaux = true
  []
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
  petsc_options_value = 'lu superlu_dist'
  line_search = none
  nl_abs_tol = 1e-4
  l_max_its = 200
  start_time = 0.0
  dt = 1.0
  num_steps = 70
[]

[Outputs]
  [out]
    type = Exodus
    execute_on = 'INITIAL TIMESTEP_END'
  []
  print_linear_residuals = false
[]

[Postprocessors]
  [total_vol]
    type = ElementIntegralVariablePostprocessor
    variable = mat_den
    execute_on = 'INITIAL TIMESTEP_END'
  []
[]

[Controls]
  [first_period]
    type = TimePeriod
    start_time = 0.0
    end_time = 10
    enable_objects = 'BCs::boundary_penalty_right'
    execute_on = 'initial timestep_begin'
  []
[]

(modules/electromagnetics/test/tests/kernels/scalar_complex_helmholtz/scalar_complex_helmholtz.i)

# problem: -(cu')' - k^2 * u = -F , 0 < x < L, u: R -> C
# u(x=0) = g0 , u(x=L) = gL
# k = a + jb
#     a = a(x) =  2 * (1 + x/L)
#     b = b(x) =      (1 + x/L)
# c = d + jh
#     d = d(x) = 12 * (1 + x/L)^2
#     h = h(x) =  4 * (1 + x/L)^2
# L = 10

[Mesh]
  [gmg]
    type = GeneratedMeshGenerator
    dim = 1
    xmin = 0
    xmax = 10
    nx = 100
  []
[]

[Variables]
  [u_real]
    order = FIRST
    family = LAGRANGE
  []
  [u_imag]
    order = FIRST
    family = LAGRANGE
  []
[]

[Functions]
  [k_real]
    type = ParsedFunction
    expression = '2*(1 + x/10)'
  []
  [k_imag]
    type = ParsedFunction
    expression = '(1 + x/10)'
  []
  [d_func]
    type = ParsedFunction
    expression = '12 * (1 + x/10)^2'
  []
  [h_func]
    type = ParsedFunction
    expression = '4 * (1 + x/10)^2'
  []
  [negative_h_func]
    type = ParsedFunction
    expression = '-4 * (1 + x/10)^2'
  []
  [RHS_real]
    type = MMSTestFunc
    L = 10
    g0_real = 1
    g0_imag = -1
    gL_real = 0
    gL_imag = 0
    component = real
  []
  [RHS_imag]
    type = MMSTestFunc
    L = 10
    g0_real = 1
    g0_imag = -1
    gL_real = 0
    gL_imag = 0
    component = imaginary
  []
[]

[Materials]
  [k_real_mat]
    type = ADGenericFunctionMaterial
    prop_names = k_real_mat
    prop_values = k_real
  []
  [k_imag_mat]
    type = ADGenericFunctionMaterial
    prop_names = k_imag_mat
    prop_values = k_imag
  []
  [wave_equation_coefficient]
    type = WaveEquationCoefficient
    k_real = k_real_mat
    k_imag = k_imag_mat
    eps_rel_real = 1
    eps_rel_imag = 0
    mu_rel_real = 1
    mu_rel_imag = 0
  []
  [negative_wave_equation_coefficient_imaginary]
    type = ADParsedMaterial
    property_name = negative_wave_equation_coefficient_imaginary
    material_property_names = wave_equation_coefficient_imaginary
    expression = '-1 * wave_equation_coefficient_imaginary'
  []
[]

[Kernels]
  [laplacian_real]
    type = FunctionDiffusion
    function = d_func
    variable = u_real
  []
  [coupledLaplacian_real]
    type = FunctionDiffusion
    function = negative_h_func
    v = u_imag
    variable = u_real
  []
  [coeffField_real]
    type = ADMatReaction
    reaction_rate = wave_equation_coefficient_real
    variable = u_real
  []
  [coupledField_real]
    type = ADMatCoupledForce
    v = u_imag
    mat_prop_coef = negative_wave_equation_coefficient_imaginary
    variable = u_real
  []
  [bodyForce_real]
    type = BodyForce
    function = RHS_real
    variable = u_real
  []
  [laplacian_imag]
    type = FunctionDiffusion
    function = d_func
    variable = u_imag
  []
  [coupledLaplacian_imag]
    type = FunctionDiffusion
    function = h_func
    v = u_real
    variable = u_imag
  []
  [coeffField_imag]
    type = ADMatReaction
    reaction_rate = wave_equation_coefficient_real
    variable = u_imag
  []
  [coupledField_imag]
    type = ADMatCoupledForce
    v = u_real
    mat_prop_coef = wave_equation_coefficient_imaginary
    variable = u_imag
  []
  [bodyForce_imag]
    type = BodyForce
    function = RHS_imag
    variable = u_imag
  []
[]

[BCs]
  [left_real]
    type = DirichletBC
    value = 1
    boundary = left
    variable = u_real
  []
  [left_imag]
    type = DirichletBC
    value = -1
    boundary = left
    variable = u_imag
  []
  [right_real]
    type = DirichletBC
    value = 0
    boundary = right
    variable = u_real
  []
  [right_imag]
    type = DirichletBC
    value = 0
    boundary = right
    variable = u_imag
  []
[]

[Executioner]
  type = Steady
  solve_type = 'PJFNK'
[]

[Outputs]
  exodus = true
[]

(modules/combined/examples/optimization/helmholtz_multimat_nostrip.i)

vol_frac = 0.35

power = 1.1

Emin = 1.0e-6
Ess = 0.475 # ss
Et = 1.0 # w

[GlobalParams]
  displacements = 'disp_x disp_y'
[]

[Mesh]
  # final_generator = 'MoveRight'
  [Bottom]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 320
    ny = 30
    xmin = 0
    xmax = 150
    ymin = 0
    ymax = 15
  []
  [RenameBottom]
    type = RenameBoundaryGenerator
    input = Bottom
    old_boundary = 'top bottom right left'
    new_boundary = 'top_bottom bottom_bottom right_bottom left_bottom'
  []
  [Top]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 320
    ny = 30
    xmin = 0
    xmax = 150
    ymin = 0
    ymax = 15
  []
  [MoveTop]
    type = TransformGenerator
    input = Top
    transform = TRANSLATE
    vector_value = '0 15 0'
  []
  [RenameTop]
    type = RenameBoundaryGenerator
    input = MoveTop
    old_boundary = 'top bottom right left'
    new_boundary = 'top_top bottom_top right_top left_top'
  []
  [bottom_gen]
    type = ParsedSubdomainMeshGenerator
    input = RenameBottom
    combinatorial_geometry = 'y <= 15'
    block_id = 1
  []
  [top_gen]
    type = ParsedSubdomainMeshGenerator
    input = RenameTop
    combinatorial_geometry = 'y > 15'
    block_id = 3
  []
  [stitch]
    type = StitchedMeshGenerator
    inputs = 'bottom_gen top_gen'
    stitch_boundaries_pairs = 'top_bottom bottom_top'
  []
  [left_load]
    type = ExtraNodesetGenerator
    input = stitch
    new_boundary = left_load
    coord = '37.5 30 0'
  []
  [right_load]
    type = ExtraNodesetGenerator
    input = left_load
    new_boundary = right_load
    coord = '112.5 30 0'
  []
  [left_support]
    type = ExtraNodesetGenerator
    input = right_load
    new_boundary = left_support
    coord = '0 0 0'
  []
  [right_support]
    type = ExtraNodesetGenerator
    input = left_support
    new_boundary = right_support
    coord = '150 0 0'
  []
[]

[Variables]
  [disp_x]
  []
  [disp_y]
  []
  [Dc]
    initial_condition = -1.0
  []
[]

[AuxVariables]
  [Cc]
    family = MONOMIAL
    order = CONSTANT
    initial_condition = -1.0
  []
  [mat_den]
    family = MONOMIAL
    order = CONSTANT
    initial_condition = ${vol_frac}
  []
  [sensitivity]
    family = MONOMIAL
    order = FIRST
    initial_condition = -1.0
    [AuxKernel]
      type = MaterialRealAux
      variable = sensitivity
      property = sensitivity
      execute_on = LINEAR
    []
  []
  [mat_den_nodal]
    family = L2_LAGRANGE
    order = FIRST
    initial_condition = ${vol_frac}
    [AuxKernel]
      type = SelfAux
      execute_on = TIMESTEP_END
      variable = mat_den_nodal
      v = mat_den
    []
  []
  [Dc_elem]
    family = MONOMIAL
    order = CONSTANT
    initial_condition = -1.0
    [AuxKernel]
      type = SelfAux
      variable = Dc_elem
      v = Dc
      execute_on = 'TIMESTEP_END'
    []
  []
[]

[Kernels]
  [diffusion]
    type = FunctionDiffusion
    variable = Dc
    function = 4.0
  []
  [potential]
    type = Reaction
    variable = Dc
  []
  [source]
    type = CoupledForce
    variable = Dc
    v = sensitivity
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [all]
    strain = SMALL
    add_variables = true
    incremental = false
  []
[]

[BCs]
  [no_y]
    type = DirichletBC
    variable = disp_y
    boundary = left_support
    value = 0.0
  []
  [no_x]
    type = DirichletBC
    variable = disp_x
    boundary = left_support
    value = 0.0
  []
  [no_y_right]
    type = DirichletBC
    variable = disp_y
    boundary = right_support
    value = 0.0
  []
  [boundary_penalty]
    type = ADRobinBC
    variable = Dc
    boundary = 'bottom_bottom right_bottom left_bottom top_top right_top left_top'
    coefficient = 10
  []
[]

[NodalKernels]
  [left_down]
    type = NodalGravity
    variable = disp_y
    boundary = left_load
    gravity_value = -1e-3
    mass = 1
  []
  [right_down]
    type = NodalGravity
    variable = disp_y
    boundary = right_load
    gravity_value = -1e-3
    mass = 1
  []
[]

[Materials]
  [elasticity_tensor_one]
    type = ComputeVariableIsotropicElasticityTensor
    youngs_modulus = E_phys_one
    poissons_ratio = poissons_ratio
    args = 'mat_den'
    block = '1'
  []
  [elasticity_tensor_three]
    type = ComputeVariableIsotropicElasticityTensor
    youngs_modulus = E_phys_three
    poissons_ratio = poissons_ratio
    args = 'mat_den'
    block = '3'
  []
  # One: Tungsten
  [E_phys_one]
    type = DerivativeParsedMaterial
    # Emin + (density^penal) * (E0 - Emin)
    expression = '${Emin} + (mat_den ^ ${power}) * (${Et}-${Emin})'
    coupled_variables = 'mat_den'
    property_name = E_phys_one
    block = '1'
    outputs = 'exodus'
  []
  # Three: SS316
  [E_phys_three]
    type = DerivativeParsedMaterial
    # Emin + (density^penal) * (E0 - Emin)
    expression = '${Emin} + (mat_den ^ ${power}) * (${Ess}-${Emin})'
    coupled_variables = 'mat_den'
    property_name = E_phys_three
    block = '3'
    outputs = 'exodus'
  []
  [poissons_ratio]
    type = GenericConstantMaterial
    prop_names = poissons_ratio
    prop_values = 0.3
  []
  [stress]
    type = ComputeLinearElasticStress
  []
  [dc_one]
    type = ComplianceSensitivity
    design_density = mat_den
    youngs_modulus = E_phys_one
    block = '1'
  []
  [dc_three]
    type = ComplianceSensitivity
    design_density = mat_den
    youngs_modulus = E_phys_three
    block = '3'
  []
[]

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

[UserObjects]
  [update_one]
    type = DensityUpdate
    density_sensitivity = Dc_elem
    design_density = mat_den
    volume_fraction = ${vol_frac}
    execute_on = TIMESTEP_BEGIN
    block = '1'
  []
  [update_three]
    type = DensityUpdate
    density_sensitivity = Dc_elem
    design_density = mat_den
    volume_fraction = ${vol_frac}
    execute_on = TIMESTEP_BEGIN
    block = '3'
  []
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
  petsc_options_value = 'lu superlu_dist'
  nl_abs_tol = 1e-10
  dt = 1.0
  num_steps = 90
[]

[Outputs]
  exodus = true
  [out]
    type = CSV
    execute_on = 'TIMESTEP_END'
  []
  print_linear_residuals = false
[]

[Postprocessors]
  [mesh_volume]
    type = VolumePostprocessor
    execute_on = 'initial timestep_end'
  []
  [total_vol]
    type = ElementIntegralVariablePostprocessor
    variable = mat_den
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [vol_frac]
    type = ParsedPostprocessor
    expression = 'total_vol / mesh_volume'
    pp_names = 'total_vol mesh_volume'
  []
  [sensitivity]
    type = ElementIntegralMaterialProperty
    mat_prop = sensitivity
    block = '1 3'
  []
  [objective_one]
    type = ElementIntegralMaterialProperty
    mat_prop = strain_energy_density
    execute_on = 'INITIAL TIMESTEP_END'
    block = '1'
  []
  [objective_three]
    type = ElementIntegralMaterialProperty
    mat_prop = strain_energy_density
    execute_on = 'INITIAL TIMESTEP_END'
    block = '3'
  []
[]

(modules/combined/test/tests/optimization/compliance_sensitivity/2d_mbb_pde.i)

vol_frac = 0.5
E0 = 1
Emin = 1e-8
power = 3
[GlobalParams]
  displacements = 'disp_x disp_y'
[]

[Mesh]
  [MeshGenerator]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 150
    ny = 50
    xmin = 0
    xmax = 30
    ymin = 0
    ymax = 10
  []
  [node]
    type = ExtraNodesetGenerator
    input = MeshGenerator
    new_boundary = hold_y
    nodes = 0
  []
  [push]
    type = ExtraNodesetGenerator
    input = node
    new_boundary = push
    coord = '30 10 0'
  []

[]

[Variables]
  [disp_x]
  []
  [disp_y]
  []
  [Dc]
    initial_condition = -1.0
  []
[]

[AuxVariables]

  [sensitivity]
    family = MONOMIAL
    order = FIRST
    initial_condition = -1.0
    [AuxKernel]
      type = MaterialRealAux
      variable = sensitivity
      property = sensitivity
      execute_on = LINEAR
    []
  []
  [mat_den]
    family = MONOMIAL
    order = CONSTANT
    initial_condition = ${vol_frac}
  []
  [Dc_elem]
    family = MONOMIAL
    order = CONSTANT
    initial_condition = -1.0
    [AuxKernel]
      type = SelfAux
      variable = Dc_elem
      v = Dc
      execute_on = 'TIMESTEP_END'
    []
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [all]
    strain = SMALL
    add_variables = true
    incremental = false
  []
[]
[Kernels]
  [diffusion]
    type = FunctionDiffusion
    variable = Dc
    function = 0.15 # radius coeff
  []
  [potential]
    type = Reaction
    variable = Dc
  []
  [source]
    type = CoupledForce
    variable = Dc
    v = sensitivity
  []
[]
[BCs]
  [no_x]
    type = DirichletBC
    variable = disp_y
    boundary = hold_y
    value = 0.0
  []
  [no_y]
    type = DirichletBC
    variable = disp_x
    boundary = right
    value = 0.0
  []
  [boundary_penalty]
    type = ADRobinBC
    variable = Dc
    boundary = 'left top'
    coefficient = 10
  []
  [boundary_penalty_right]
    type = ADRobinBC
    variable = Dc
    boundary = 'right'
    coefficient = 10
  []
[]
[NodalKernels]
  [push]
    type = NodalGravity
    variable = disp_y
    boundary = push
    gravity_value = -1
    mass = 1
  []
[]

[Materials]
  [elasticity_tensor]
    type = ComputeVariableIsotropicElasticityTensor
    youngs_modulus = E_phys
    poissons_ratio = poissons_ratio
    args = 'mat_den'
  []
  [E_phys]
    type = DerivativeParsedMaterial
    # Emin + (density^penal) * (E0 - Emin)
    expression = '${Emin} + (mat_den ^ ${power}) * (${E0}-${Emin})'
    coupled_variables = 'mat_den'
    property_name = E_phys
  []
  [poissons_ratio]
    type = GenericConstantMaterial
    prop_names = poissons_ratio
    prop_values = 0.3
  []
  [stress]
    type = ComputeLinearElasticStress
  []
  [dc]
    type = ComplianceSensitivity
    design_density = mat_den
    youngs_modulus = E_phys
    incremental = false
  []
[]

[Preconditioning]
  [smp]
    type = SMP
    full = true
  []
[]
[UserObjects]
  [update]
    type = DensityUpdate
    density_sensitivity = Dc_elem
    design_density = mat_den
    volume_fraction = ${vol_frac}
    execute_on = TIMESTEP_BEGIN
    force_postaux = true
  []
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
  petsc_options_value = 'lu superlu_dist'
  line_search = none
  nl_abs_tol = 1e-4
  l_max_its = 200
  start_time = 0.0
  dt = 1.0
  num_steps = 70
[]

[Outputs]
  [out]
    type = CSV
    execute_on = 'INITIAL TIMESTEP_END'
  []
  print_linear_residuals = false
[]

[Postprocessors]
  [total_vol]
    type = ElementIntegralVariablePostprocessor
    variable = mat_den
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [sensitivity]
    type = ElementIntegralMaterialProperty
    mat_prop = sensitivity
  []
[]

[Controls]
  [first_period]
    type = TimePeriod
    start_time = 0.0
    end_time = 10
    enable_objects = 'BCs::boundary_penalty_right'
    execute_on = 'initial timestep_begin'
  []
[]

(test/tests/preconditioners/reuse/convergence.i)

# Simple 3D test with diffusion, setup to make sure
# there is a sensible difference in the linear iteration
# counts with re-use versus without re-use

[Variables]
  [u]
  []
[]

[Mesh]
  [msh]
    type = GeneratedMeshGenerator
    dim = 3
    nx = 4
    ny = 4
    nz = 4
  []
[]

[Kernels]
  [diffusion]
    type = FunctionDiffusion
    variable = u
    function = 'arg'
  []
  [time]
    type = TimeDerivative
    variable = u
  []
  [body_force]
    type = BodyForce
    variable = u
    function = body
  []
[]

[Functions]
  [body]
    type = ParsedFunction
    expression = 100*sin(t)
  []
  [arg]
    type = ParsedFunction
    expression = 'x*y*z*cos(t)+1'
  []
[]

[BCs]
  [fix_concentration]
    type = DirichletBC
    preset = true
    boundary = left
    variable = u
    value = 0.0
  []
[]

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

[Executioner]
  type = Transient

  solve_type = 'newton'
  line_search = none

  petsc_options = ''
  petsc_options_iname = '-pc_type -ksp_type'
  petsc_options_value = 'lu gmres'
  l_tol = 1e-8
  l_max_its = 100

  reuse_preconditioner = false
  reuse_preconditioner_max_linear_its = 10

  nl_max_its = 10
  nl_rel_tol = 1e-8
  nl_abs_tol = 1e-10

  start_time = 0.0
  dt = 1.0
  dtmin = 1.0
  end_time = 10.0

  [./Adaptivity]
    interval = 5
    max_h_level = 1
    start_time = 11.0
    stop_time = 6.0
  [../]

[]

[Reporters/iteration_info]
  type = IterationInfo
[]

[Outputs]
  exodus = false
  [./csv]
    type = CSV
    file_base = base_case
  [../]
[]

(modules/combined/examples/optimization/helmholtz_multimat_strip.i)

vol_frac = 0.35

power = 1.1

Emin = 1.0e-6
Ess = 0.475 # ss
Et = 1.0 # w

[GlobalParams]
  displacements = 'disp_x disp_y'
[]

[Mesh]
  # final_generator = 'MoveRight'
  [Bottom]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 320
    ny = 30
    xmin = 0
    xmax = 150
    ymin = 0
    ymax = 15
  []
  [RenameBottom]
    type = RenameBoundaryGenerator
    input = Bottom
    old_boundary = 'top bottom right left'
    new_boundary = 'top_bottom bottom_bottom right_bottom left_bottom'
  []
  [Middle]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 320
    ny = 6
    xmin = 0
    xmax = 150
    ymin = 0
    ymax = 3
  []
  [MoveMiddle]
    type = TransformGenerator
    input = Middle
    transform = TRANSLATE
    vector_value = '0 15 0'
  []
  [RenameMiddle]
    type = RenameBoundaryGenerator
    input = MoveMiddle
    old_boundary = 'top bottom right left'
    new_boundary = 'top_middle bottom_middle right_middle left_middle'
  []
  [Top]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 320
    ny = 30
    xmin = 0
    xmax = 150
    ymin = 0
    ymax = 15
  []
  [MoveTop]
    type = TransformGenerator
    input = Top
    transform = TRANSLATE
    vector_value = '0 18 0'
  []
  [RenameTop]
    type = RenameBoundaryGenerator
    input = MoveTop
    old_boundary = 'top bottom right left'
    new_boundary = 'top_top bottom_top right_top left_top'
  []
  [bottom_gen]
    type = ParsedSubdomainMeshGenerator
    input = RenameBottom
    combinatorial_geometry = 'y <= 15'
    block_id = 1
  []
  [middle_gen]
    type = ParsedSubdomainMeshGenerator
    input = RenameMiddle
    combinatorial_geometry = 'y <= 18 & y > 15'
    block_id = 2
  []
  [top_gen]
    type = ParsedSubdomainMeshGenerator
    input = RenameTop
    combinatorial_geometry = 'y > 18'
    block_id = 3
  []
  [stitch]
    type = StitchedMeshGenerator
    inputs = 'bottom_gen middle_gen top_gen'
    stitch_boundaries_pairs = 'top_bottom bottom_middle; top_middle bottom_top'
  []
  [left_load]
    type = ExtraNodesetGenerator
    input = stitch
    new_boundary = left_load
    coord = '37.5 33 0'
  []
  [right_load]
    type = ExtraNodesetGenerator
    input = left_load
    new_boundary = right_load
    coord = '112.5 33 0'
  []
  [left_support]
    type = ExtraNodesetGenerator
    input = right_load
    new_boundary = left_support
    coord = '0 0 0'
  []
  [right_support]
    type = ExtraNodesetGenerator
    input = left_support
    new_boundary = right_support
    coord = '150 0 0'
  []
[]

[Variables]
  [disp_x]
  []
  [disp_y]
  []
  [Dc]
    initial_condition = -1.0
  []
[]

[AuxVariables]
  [Cc]
    family = MONOMIAL
    order = CONSTANT
    initial_condition = -1.0
  []
  [mat_den]
    family = MONOMIAL
    order = CONSTANT
    initial_condition = ${vol_frac}
  []

  [sensitivity]
    family = MONOMIAL
    order = FIRST
    initial_condition = -1.0
    [AuxKernel]
      type = MaterialRealAux
      variable = sensitivity
      property = sensitivity
      execute_on = LINEAR
    []
    block = '1 2 3'
  []
  [mat_den_nodal]
    family = L2_LAGRANGE
    order = FIRST
    initial_condition = ${vol_frac}
    [AuxKernel]
      type = SelfAux
      execute_on = TIMESTEP_END
      variable = mat_den_nodal
      v = mat_den
    []
  []
  [Dc_elem]
    family = MONOMIAL
    order = CONSTANT
    initial_condition = -1.0
    [AuxKernel]
      type = SelfAux
      variable = Dc_elem
      v = Dc
      execute_on = 'TIMESTEP_END'
    []
  []
[]

[Kernels]
  [diffusion]
    type = FunctionDiffusion
    variable = Dc
    function = 4.0
  []
  [potential]
    type = Reaction
    variable = Dc
  []
  [source]
    type = CoupledForce
    variable = Dc
    v = sensitivity
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [all]
    strain = SMALL
    add_variables = true
    incremental = false
  []
[]

[BCs]
  [no_y]
    type = DirichletBC
    variable = disp_y
    boundary = left_support
    value = 0.0
  []
  [no_x]
    type = DirichletBC
    variable = disp_x
    boundary = left_support
    value = 0.0
  []
  [no_y_right]
    type = DirichletBC
    variable = disp_y
    boundary = right_support
    value = 0.0
  []
  [boundary_penalty]
    type = ADRobinBC
    variable = Dc
    boundary = 'bottom_bottom right_bottom left_bottom top_top right_top left_top left_middle '
               'right_middle'
    coefficient = 10
  []
[]

[NodalKernels]
  [left_down]
    type = NodalGravity
    variable = disp_y
    boundary = left_load
    gravity_value = -1e-3
    mass = 1
  []
  [right_down]
    type = NodalGravity
    variable = disp_y
    boundary = right_load
    gravity_value = -1e-3
    mass = 1
  []
[]

[Materials]
  [sensitivity]
    type = ParsedMaterial
    property_name = 'sensitivity'
    block = '2'
    expression = '0'
  []
  [elasticity_tensor_one]
    type = ComputeVariableIsotropicElasticityTensor
    youngs_modulus = E_phys_one
    poissons_ratio = poissons_ratio
    args = 'mat_den'
    block = '1'
  []
  [elasticity_tensor_three]
    type = ComputeVariableIsotropicElasticityTensor
    youngs_modulus = E_phys_three
    poissons_ratio = poissons_ratio
    args = 'mat_den'
    block = '3'
  []
  [elasticity_tensor_two]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 1.0
    poissons_ratio = 0.3
    block = '2'
  []
  # One: Tungsten
  [E_phys_one]
    type = DerivativeParsedMaterial
    # Emin + (density^penal) * (E0 - Emin)
    expression = '${Emin} + (mat_den ^ ${power}) * (${Et}-${Emin})'
    coupled_variables = 'mat_den'
    property_name = E_phys_one
    block = '1'
    outputs = 'exodus'
  []
  # Three: SS316
  [E_phys_three]
    type = DerivativeParsedMaterial
    # Emin + (density^penal) * (E0 - Emin)
    expression = '${Emin} + (mat_den ^ ${power}) * (${Ess}-${Emin})'
    coupled_variables = 'mat_den'
    property_name = E_phys_three
    block = '3'
    outputs = 'exodus'
  []
  [poissons_ratio]
    type = GenericConstantMaterial
    prop_names = poissons_ratio
    prop_values = 0.3
  []
  [stress]
    type = ComputeLinearElasticStress
  []
  [dc_one]
    type = ComplianceSensitivity
    design_density = mat_den
    youngs_modulus = E_phys_one
    block = '1'
  []
  [dc_three]
    type = ComplianceSensitivity
    design_density = mat_den
    youngs_modulus = E_phys_three
    block = '3'
  []
[]

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

[UserObjects]
  [update_one]
    type = DensityUpdate
    density_sensitivity = Dc_elem
    design_density = mat_den
    volume_fraction = ${vol_frac}
    execute_on = TIMESTEP_BEGIN
    block = '1'
  []
  [update_three]
    type = DensityUpdate
    density_sensitivity = Dc_elem
    design_density = mat_den
    volume_fraction = ${vol_frac}
    execute_on = TIMESTEP_BEGIN
    block = '3'
  []
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
  petsc_options_value = 'lu superlu_dist'
  nl_abs_tol = 1e-10
  dt = 1.0
  num_steps = 90
[]

[Outputs]
  exodus = true
  [out]
    type = CSV
    execute_on = 'TIMESTEP_END'
  []
  print_linear_residuals = false
[]

[Postprocessors]
  [mesh_volume]
    type = VolumePostprocessor
    execute_on = 'initial timestep_end'
  []
  [total_vol]
    type = ElementIntegralVariablePostprocessor
    variable = mat_den
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [vol_frac]
    type = ParsedPostprocessor
    expression = 'total_vol / mesh_volume'
    pp_names = 'total_vol mesh_volume'
  []
  [sensitivity]
    type = ElementIntegralMaterialProperty
    mat_prop = sensitivity
    block = '1 3'
  []
  [objective_one]
    type = ElementIntegralMaterialProperty
    mat_prop = strain_energy_density
    execute_on = 'INITIAL TIMESTEP_END'
    block = '1'
  []
  [objective_three]
    type = ElementIntegralMaterialProperty
    mat_prop = strain_energy_density
    execute_on = 'INITIAL TIMESTEP_END'
    block = '3'
  []
[]

(modules/combined/examples/optimization/2d_mbb_pde_amr.i)

vol_frac = 0.5
E0 = 1
Emin = 1e-8
power = 2
[GlobalParams]
  displacements = 'disp_x disp_y'
[]

[Mesh]
  [MeshGenerator]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 30
    ny = 10
    xmin = 0
    xmax = 30
    ymin = 0
    ymax = 10
  []
  [node]
    type = ExtraNodesetGenerator
    input = MeshGenerator
    new_boundary = pull
    nodes = 0
  []
  [push]
    type = ExtraNodesetGenerator
    input = node
    new_boundary = push
    coord = '30 10 0'
  []

[]

[Variables]
  [Dc]
    initial_condition = -1.0
  []
[]

[AuxVariables]
  [Emin]
    family = MONOMIAL
    order = CONSTANT
    initial_condition = ${Emin}
  []
  [power]
    family = MONOMIAL
    order = CONSTANT
    initial_condition = ${power}
  []
  [E0]
    family = MONOMIAL
    order = CONSTANT
    initial_condition = ${E0}
  []
  [sensitivity]
    family = MONOMIAL
    order = FIRST
    initial_condition = -1.0
    [AuxKernel]
      type = MaterialRealAux
      variable = sensitivity
      property = sensitivity
      execute_on = LINEAR
    []
  []
  [mat_den]
    family = MONOMIAL
    order = CONSTANT
    initial_condition = ${vol_frac}
  []
  [Dc_elem]
    family = MONOMIAL
    order = CONSTANT
    initial_condition = -1.0
    [AuxKernel]
      type = SelfAux
      variable = Dc_elem
      v = Dc
      execute_on = 'TIMESTEP_END'
    []
  []
  [mat_den_nodal]
    family = L2_LAGRANGE
    order = FIRST
    initial_condition = ${vol_frac}
    [AuxKernel]
      type = SelfAux
      execute_on = TIMESTEP_END
      variable = mat_den_nodal
      v = mat_den
    []
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [all]
    strain = SMALL
    add_variables = true
    incremental = false
  []
[]
[Kernels]
  [diffusion]
    type = FunctionDiffusion
    variable = Dc
    function = 0.15 # radius coeff
  []
  [potential]
    type = Reaction
    variable = Dc
  []
  [source]
    type = CoupledForce
    variable = Dc
    v = sensitivity
  []
[]
[BCs]
  [no_x]
    type = DirichletBC
    variable = disp_y
    boundary = pull
    value = 0.0
  []
  [no_y]
    type = DirichletBC
    variable = disp_x
    boundary = right
    value = 0.0
  []
  [boundary_penalty]
    type = ADRobinBC
    variable = Dc
    boundary = 'left top'
    coefficient = 10
  []
  [boundary_penalty_right]
    type = ADRobinBC
    variable = Dc
    boundary = 'right'
    coefficient = 10
  []
[]
[NodalKernels]
  [pull]
    type = NodalGravity
    variable = disp_y
    boundary = push
    gravity_value = -1
    mass = 1
  []
[]

[Materials]
  [elasticity_tensor]
    type = ComputeVariableIsotropicElasticityTensor
    youngs_modulus = E_phys
    poissons_ratio = poissons_ratio
    args = 'Emin mat_den power E0'
  []
  [E_phys]
    type = DerivativeParsedMaterial
    # Emin + (density^penal) * (E0 - Emin)
    expression = '${Emin} + (mat_den ^ ${power}) * (${E0}-${Emin})'
    coupled_variables = 'mat_den'
    property_name = E_phys
  []
  [poissons_ratio]
    type = GenericConstantMaterial
    prop_names = poissons_ratio
    prop_values = 0.3
  []
  [stress]
    type = ComputeLinearElasticStress
  []
  [dc]
    type = ComplianceSensitivity
    design_density = mat_den
    youngs_modulus = E_phys
    incremental = false
  []
[]

[Preconditioning]
  [smp]
    type = SMP
    full = true
  []
[]
[UserObjects]
  [update]
    type = DensityUpdate
    density_sensitivity = Dc_elem
    design_density = mat_den
    volume_fraction = ${vol_frac}
    execute_on = TIMESTEP_BEGIN
    force_postaux = true
  []
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
  petsc_options_value = 'lu superlu_dist'
  line_search = none
  nl_abs_tol = 1e-4
  l_max_its = 200
  start_time = 0.0
  dt = 1.0
  num_steps = 70
[]

[Outputs]
  [out]
    type = Exodus
    execute_on = 'INITIAL TIMESTEP_END'
  []
  print_linear_residuals = false
[]

[Postprocessors]
  [total_vol]
    type = ElementIntegralVariablePostprocessor
    variable = mat_den
    execute_on = 'INITIAL TIMESTEP_END'
  []
[]

[Controls]
  [first_period]
    type = TimePeriod
    start_time = 0.0
    end_time = 40
    enable_objects = 'BCs::boundary_penalty_right'
    execute_on = 'initial timestep_begin'
  []
[]

[Adaptivity]
  max_h_level = 2
  recompute_markers_during_cycles = true
  interval = 1
  cycles_per_step = 1
  marker = density_marker
  [Indicators]
    [density_jump]
      type = ValueJumpIndicator
      variable = mat_den_nodal
    []
  []
  [Markers]
    [density_marker]
      type = ErrorToleranceMarker
      indicator = density_jump
      coarsen = 0.1
      refine = 0.1
    []
  []
[]

(modules/combined/examples/optimization/3d_mbb.i)

vol_frac = 0.5
E0 = 1
Emin = 1e-8
power = 2
[GlobalParams]
  displacements = 'disp_x disp_y disp_z'
[]

[Mesh]
  [MeshGenerator]
    type = GeneratedMeshGenerator
    dim = 3
    nx = 60
    ny = 20
    nz = 20
    xmin = 0
    xmax = 30
    ymin = 0
    ymax = 10
    zmin = 0
    zmax = 10
  []
  [node]
    type = ExtraNodesetGenerator
    input = MeshGenerator
    new_boundary = hold_y
    coord = '0 0 0; 0 0 10'
  []
  [push]
    type = ExtraNodesetGenerator
    input = node
    new_boundary = push
    coord = '30 10 5'
  []
[]

[Variables]
  [disp_z]
  []
  [Dc]
    initial_condition = -1.0
  []
[]

[AuxVariables]
  [Emin]
    family = MONOMIAL
    order = CONSTANT
    initial_condition = ${Emin}
  []
  [power]
    family = MONOMIAL
    order = CONSTANT
    initial_condition = ${power}
  []
  [E0]
    family = MONOMIAL
    order = CONSTANT
    initial_condition = ${E0}
  []
  [sensitivity]
    family = MONOMIAL
    order = FIRST
    initial_condition = -1.0
    [AuxKernel]
      type = MaterialRealAux
      variable = sensitivity
      property = sensitivity
      execute_on = LINEAR
    []
  []
  [mat_den]
    family = MONOMIAL
    order = CONSTANT
    initial_condition = ${vol_frac}
  []
  [Dc_elem]
    family = MONOMIAL
    order = CONSTANT
    initial_condition = -1.0
    [AuxKernel]
      type = SelfAux
      variable = Dc_elem
      v = Dc
      execute_on = 'TIMESTEP_END'
    []
  []
  [mat_den_nodal]
    family = L2_LAGRANGE
    order = FIRST
    initial_condition = ${vol_frac}
    [AuxKernel]
      type = SelfAux
      execute_on = TIMESTEP_END
      variable = mat_den_nodal
      v = mat_den
    []
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [all]
    strain = SMALL
    add_variables = true
    incremental = false
  []
[]
[Kernels]
  [diffusion]
    type = FunctionDiffusion
    variable = Dc
    function = 0.15 # radius coeff
  []
  [potential]
    type = Reaction
    variable = Dc
  []
  [source]
    type = CoupledForce
    variable = Dc
    v = sensitivity
  []
[]
[BCs]
  [no_x]
    type = DirichletBC
    variable = disp_y
    boundary = hold_y
    value = 0.0
  []
  [no_y]
    type = DirichletBC
    variable = disp_x
    boundary = right
    value = 0.0
  []
  [boundary_penalty]
    type = ADRobinBC
    variable = Dc
    boundary = 'left top front back'
    coefficient = 10
  []
  [boundary_penalty_right]
    type = ADRobinBC
    variable = Dc
    boundary = 'right'
    coefficient = 10
  []
[]
[NodalKernels]
  [push]
    type = NodalGravity
    variable = disp_y
    boundary = push
    gravity_value = -1
    mass = 1
  []
[]

[Materials]
  [elasticity_tensor]
    type = ComputeVariableIsotropicElasticityTensor
    youngs_modulus = E_phys
    poissons_ratio = poissons_ratio
    args = 'Emin mat_den power E0'
  []
  [E_phys]
    type = DerivativeParsedMaterial
    # Emin + (density^penal) * (E0 - Emin)
    expression = '${Emin} + (mat_den ^ ${power}) * (${E0}-${Emin})'
    coupled_variables = 'mat_den'
    property_name = E_phys
  []
  [poissons_ratio]
    type = GenericConstantMaterial
    prop_names = poissons_ratio
    prop_values = 0.3
  []
  [stress]
    type = ComputeLinearElasticStress
  []
  [dc]
    type = ComplianceSensitivity
    design_density = mat_den
    youngs_modulus = E_phys
    incremental = false
  []
[]

[Preconditioning]
  [smp]
    type = SMP
    full = true
  []
[]
[UserObjects]
  [update]
    type = DensityUpdate
    density_sensitivity = Dc_elem
    design_density = mat_den
    volume_fraction = ${vol_frac}
    execute_on = TIMESTEP_BEGIN
    force_postaux = true
  []
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
  petsc_options_value = 'lu superlu_dist'
  line_search = none
  nl_abs_tol = 1e-4
  l_max_its = 200
  start_time = 0.0
  dt = 1.0
  num_steps = 70
[]

[Outputs]
  [out]
    type = Exodus
    execute_on = 'INITIAL TIMESTEP_END'
  []
  print_linear_residuals = false
[]

[Postprocessors]
  [total_vol]
    type = ElementIntegralVariablePostprocessor
    variable = mat_den
    execute_on = 'INITIAL TIMESTEP_END'
  []
[]

[Controls]
  [first_period]
    type = TimePeriod
    start_time = 0.0
    end_time = 10
    enable_objects = 'BCs::boundary_penalty_right'
    execute_on = 'initial timestep_begin'
  []
[]

Overview
Example Input File Syntax
Input Parameters
Input Files