Ordered SIMP: Multimaterial

In addition to establish a maximum volume fraction and optimize the distribution of material over the domain, one can optimize multiple materials within the same domain and assign, additionally, individual cost to each of them (see Tavakoli and Mohseni (2014)). This allows to enforce another constraint: The total cost over the domain.

This approach uses the object DensityUpdateTwoConstraints, which allows multiple material optimization to be applied to blocks individually. The materials in this example have pseudo-densities of 0.4, 0.7, and 1.0; in addition to the void material of 0.0.

Cost and Young's modulus are defined as interpolated quantities across the three materials:

Listing 1: MBB Material interpolation for multimaterial optimization

[Materials<<<{"href": "../../../../syntax/Materials/index.html"}>>>]
  [elasticity_tensor]
    type = ComputeVariableIsotropicElasticityTensor<<<{"description": "Compute an isotropic elasticity tensor for elastic constants that change as a function of material properties", "href": "../../../../source/materials/ComputeVariableIsotropicElasticityTensor.html"}>>>
    youngs_modulus<<<{"description": "Name of material property defining the Young's Modulus"}>>> = E_phys
    poissons_ratio<<<{"description": "Name of material property defining the Poisson's Ratio"}>>> = poissons_ratio
    args<<<{"description": "Variable dependence for the Young's Modulus and Poisson's Ratio materials"}>>> = 'mat_den'
  []
  [E_phys]
    type = DerivativeParsedMaterial<<<{"description": "Parsed Function Material with automatic derivatives.", "href": "../../../../source/materials/DerivativeParsedMaterial.html"}>>>
    # ordered multimaterial simp
    expression<<<{"description": "Parsed function (see FParser) expression for the parsed material"}>>> = "A1:=(${E0}-${E1})/(${rho0}^${power}-${rho1}^${power}); B1:=${E0}-A1*${rho0}^${power}; E1:=A1*mat_den^${power}+B1; A2:=(${E1}-${E2})/(${rho1}^${power}-${rho2}^${power}); B2:=${E1}-A2*${rho1}^${power}; E2:=A2*mat_den^${power}+B2; A3:=(${E2}-${E3})/(${rho2}^${power}-${rho3}^${power}); B3:=${E2}-A3*${rho2}^${power}; E3:=A3*mat_den^${power}+B3; if(mat_den<${rho1},E1,if(mat_den<${rho2},E2,E3))"
    coupled_variables<<<{"description": "Vector of variables used in the parsed function"}>>> = 'mat_den'
    property_name<<<{"description": "Name of the parsed material property"}>>> = E_phys
  []

  [Cost_mat]
    type = DerivativeParsedMaterial<<<{"description": "Parsed Function Material with automatic derivatives.", "href": "../../../../source/materials/DerivativeParsedMaterial.html"}>>>
    # ordered multimaterial simp
    expression<<<{"description": "Parsed function (see FParser) expression for the parsed material"}>>> = "A1:=(${C0}-${C1})/(${rho0}^(1/${power})-${rho1}^(1/${power})); B1:=${C0}-A1*${rho0}^(1/${power}); C1:=A1*mat_den^(1/${power})+B1; A2:=(${C1}-${C2})/(${rho1}^(1/${power})-${rho2}^(1/${power})); B2:=${C1}-A2*${rho1}^(1/${power}); C2:=A2*mat_den^(1/${power})+B2; A3:=(${C2}-${C3})/(${rho2}^(1/${power})-${rho3}^(1/${power})); B3:=${C2}-A3*${rho2}^(1/${power}); C3:=A3*mat_den^(1/${power})+B3; if(mat_den<${rho1},C1,if(mat_den<${rho2},C2,C3))"
    coupled_variables<<<{"description": "Vector of variables used in the parsed function"}>>> = 'mat_den'
    property_name<<<{"description": "Name of the parsed material property"}>>> = Cost_mat
  []
  [CostDensity]
    type = ParsedMaterial<<<{"description": "Parsed expression Material.", "href": "../../../../source/materials/ParsedMaterial.html"}>>>
    property_name<<<{"description": "Name of the parsed material property"}>>> = CostDensity
    coupled_variables<<<{"description": "Vector of variables used in the parsed function"}>>> = 'mat_den Cost'
    expression<<<{"description": "Parsed function (see FParser) expression for the parsed material"}>>> = 'mat_den*Cost'
  []
  [poissons_ratio]
    type = GenericConstantMaterial<<<{"description": "Declares material properties based on names and values prescribed by input parameters.", "href": "../../../../source/materials/GenericConstantMaterial.html"}>>>
    prop_names<<<{"description": "The names of the properties this material will have"}>>> = poissons_ratio
    prop_values<<<{"description": "The values associated with the named properties"}>>> = 0.3
  []
  [stress]
    type = ComputeLinearElasticStress<<<{"description": "Compute stress using elasticity for small strains", "href": "../../../../source/materials/ComputeLinearElasticStress.html"}>>>
  []
  [dc]
    type = ComplianceSensitivity<<<{"description": "Computes compliance sensitivity needed for SIMP method.", "href": "../../../../source/materials/ComplianceSensitivity.html"}>>>
    design_density<<<{"description": "Design density variable name."}>>> = mat_den
    youngs_modulus<<<{"description": "DerivativeParsedMaterial for Youngs modulus."}>>> = E_phys
  []
  [cc]
    type = CostSensitivity<<<{"description": "Computes cost sensitivity needed for multimaterial SIMP method.", "href": "../../../../source/materials/CostSensitivity.html"}>>>
    design_density<<<{"description": "Design density variable name."}>>> = mat_den
    cost<<<{"description": "DerivativeParsedMaterial for cost of materials."}>>> = Cost_mat
    outputs<<<{"description": "Vector of output names where you would like to restrict the output of variables(s) associated with this object"}>>> = 'exodus'
  []
[]

The optimization process and the enforcement of volume and cost constraints are driven by the DensityUpdateTwoConstraints, radial average, and sensitivity filter user objects:

Listing 2: MBB User objects for multimaterial optimization

[UserObjects<<<{"href": "../../../../syntax/UserObjects/index.html"}>>>]
  [rad_avg]
    type = RadialAverage<<<{"description": "Perform a radial average of a material property", "href": "../../../../source/userobjects/RadialAverage.html"}>>>
    radius<<<{"description": "Cut-off radius for the averaging"}>>> = 3
    weights<<<{"description": "Distance based weight function"}>>> = linear
    prop_name<<<{"description": "Name of the material property to average"}>>> = sensitivity
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = TIMESTEP_END
    force_preaux<<<{"description": "Forces the UserObject to be executed in PREAUX"}>>> = true
  []
  [rad_avg_cost]
    type = RadialAverage<<<{"description": "Perform a radial average of a material property", "href": "../../../../source/userobjects/RadialAverage.html"}>>>
    radius<<<{"description": "Cut-off radius for the averaging"}>>> = 3
    weights<<<{"description": "Distance based weight function"}>>> = linear
    prop_name<<<{"description": "Name of the material property to average"}>>> = cost_sensitivity
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = TIMESTEP_END
    force_preaux<<<{"description": "Forces the UserObject to be executed in PREAUX"}>>> = true
  []
  [update]
    type = DensityUpdateTwoConstraints<<<{"description": "Compute updated densities based on sensitivities using an optimality criteria method to keep the volume and cost constraints satisified.", "href": "../../../../source/userobjects/DensityUpdateTwoConstraints.html"}>>>
    # This is
    density_sensitivity<<<{"description": "Name of the density_sensitivity variable."}>>> = Dc
    cost_density_sensitivity<<<{"description": "Name of the cost density sensitivity variable."}>>> = Cc
    cost<<<{"description": "Name of the cost variable."}>>> = Cost
    cost_fraction<<<{"description": "Cost fraction."}>>> = ${cost_frac}
    design_density<<<{"description": "Design density variable name."}>>> = mat_den
    volume_fraction<<<{"description": "Volume fraction."}>>> = ${vol_frac}
    bisection_lower_bound<<<{"description": "Lower bound for the bisection algorithm."}>>> = 0
    bisection_upper_bound<<<{"description": "Upper bound for the bisection algorithm."}>>> = 1.0e16 # 100
    bisection_move<<<{"description": "Bisection move for the updated solution."}>>> = 0.05
    adaptive_move<<<{"description": "Whether incremental moves in the bisection algorithm will be reduced as the number of iterations increases. Note that the time must correspond to iteration number for better results."}>>> = true
    relative_tolerance<<<{"description": "Relative tolerance on both compliance and cost to end the bisection method."}>>> = 1.0e-3
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = TIMESTEP_BEGIN
  []
  # Provides Dc
  [calc_sense]
    type = SensitivityFilter<<<{"description": "Computes the filtered sensitivity using a radial average user object.", "href": "../../../../source/userobjects/SensitivityFilter.html"}>>>
    density_sensitivity<<<{"description": "Name of the density_sensitivity variable."}>>> = Dc
    design_density<<<{"description": "Design density variable name."}>>> = mat_den
    filter_UO<<<{"description": "Radial Average user object"}>>> = rad_avg
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = TIMESTEP_END
    force_postaux<<<{"description": "Forces the UserObject to be executed in POSTAUX"}>>> = true
  []
  # Provides Cc
  [calc_sense_cost]
    type = SensitivityFilter<<<{"description": "Computes the filtered sensitivity using a radial average user object.", "href": "../../../../source/userobjects/SensitivityFilter.html"}>>>
    density_sensitivity<<<{"description": "Name of the density_sensitivity variable."}>>> = Cc
    design_density<<<{"description": "Design density variable name."}>>> = mat_den
    filter_UO<<<{"description": "Radial Average user object"}>>> = rad_avg_cost
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = TIMESTEP_END
    force_postaux<<<{"description": "Forces the UserObject to be executed in POSTAUX"}>>> = true
  []
[]

The final distribution of materials for the bridge structure (note that symmetry is being used) looks as follows:

References

Rouhollah Tavakoli and Seyyed Mohammad Mohseni. Alternating active-phase algorithm for multimaterial topology optimization problems: a 115-line matlab implementation. Structural and Multidisciplinary Optimization, 49:621–642, 2014.

@article{tavakoli2014alternating,
    author = "Tavakoli, Rouhollah and Mohseni, Seyyed Mohammad",
    title = "Alternating active-phase algorithm for multimaterial topology optimization problems: a 115-line MATLAB implementation",
    journal = "Structural and Multidisciplinary Optimization",
    volume = "49",
    pages = "621--642",
    year = "2014",
    publisher = "Springer"
}

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

vol_frac = 0.4
cost_frac = 0.3

power = 4

E0 = 1.0e-6
E1 = 0.2
E2 = 0.6
E3 = 1.0

rho0 = 1.0e-6
rho1 = 0.4
rho2 = 0.7
rho3 = 1.0

C0 = 1.0e-6
C1 = 0.5
C2 = 0.8
C3 = 1.0

[GlobalParams]
  displacements = 'disp_x disp_y'
[]

[Mesh]
  [MeshGenerator]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 50
    ny = 50
    xmin = 0
    xmax = 50
    ymin = 0
    ymax = 50
  []
  [node]
    type = ExtraNodesetGenerator
    input = MeshGenerator
    new_boundary = hold
    nodes = 0
  []
  [push_left]
    type = ExtraNodesetGenerator
    input = node
    new_boundary = push_left
    coord = '25 0 0'
  []
  [push_center]
    type = ExtraNodesetGenerator
    input = push_left
    new_boundary = push_center
    coord = '50 0 0'
  []
[]

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

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

[AuxKernels]
  [Cost]
    type = MaterialRealAux
    variable = Cost
    property = Cost_mat
  []
[]

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

[BCs]
  [no_y]
    type = DirichletBC
    variable = disp_y
    boundary = hold
    value = 0.0
  []
  [no_x_symm]
    type = DirichletBC
    variable = disp_x
    boundary = right
    value = 0.0
  []

[]
[NodalKernels]
  [push_left]
    type = NodalGravity
    variable = disp_y
    boundary = push_left
    gravity_value = -1e-3
    mass = 1
  []
  [push_center]
    type = NodalGravity
    variable = disp_y
    boundary = push_center
    gravity_value = -1e-3
    mass = 1
  []
[]
[Materials]
  [elasticity_tensor]
    type = ComputeVariableIsotropicElasticityTensor
    youngs_modulus = E_phys
    poissons_ratio = poissons_ratio
    args = 'mat_den'
  []
  [E_phys]
    type = DerivativeParsedMaterial
    # ordered multimaterial simp
    expression = "A1:=(${E0}-${E1})/(${rho0}^${power}-${rho1}^${power}); "
                 "B1:=${E0}-A1*${rho0}^${power}; E1:=A1*mat_den^${power}+B1; "
                 "A2:=(${E1}-${E2})/(${rho1}^${power}-${rho2}^${power}); "
                 "B2:=${E1}-A2*${rho1}^${power}; E2:=A2*mat_den^${power}+B2; "
                 "A3:=(${E2}-${E3})/(${rho2}^${power}-${rho3}^${power}); "
                 "B3:=${E2}-A3*${rho2}^${power}; E3:=A3*mat_den^${power}+B3; "
                 "if(mat_den<${rho1},E1,if(mat_den<${rho2},E2,E3))"
    coupled_variables = 'mat_den'
    property_name = E_phys
  []

  [Cost_mat]
    type = DerivativeParsedMaterial
    # ordered multimaterial simp
    expression = "A1:=(${C0}-${C1})/(${rho0}^(1/${power})-${rho1}^(1/${power})); "
                 "B1:=${C0}-A1*${rho0}^(1/${power}); C1:=A1*mat_den^(1/${power})+B1; "
                 "A2:=(${C1}-${C2})/(${rho1}^(1/${power})-${rho2}^(1/${power})); "
                 "B2:=${C1}-A2*${rho1}^(1/${power}); C2:=A2*mat_den^(1/${power})+B2; "
                 "A3:=(${C2}-${C3})/(${rho2}^(1/${power})-${rho3}^(1/${power})); "
                 "B3:=${C2}-A3*${rho2}^(1/${power}); C3:=A3*mat_den^(1/${power})+B3; "
                 "if(mat_den<${rho1},C1,if(mat_den<${rho2},C2,C3))"
    coupled_variables = 'mat_den'
    property_name = Cost_mat
  []
  [CostDensity]
    type = ParsedMaterial
    property_name = CostDensity
    coupled_variables = 'mat_den Cost'
    expression = 'mat_den*Cost'
  []
  [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
  []
  [cc]
    type = CostSensitivity
    design_density = mat_den
    cost = Cost_mat
    outputs = 'exodus'
  []
[]

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

[UserObjects]
  [rad_avg]
    type = RadialAverage
    radius = 3
    weights = linear
    prop_name = sensitivity
    execute_on = TIMESTEP_END
    force_preaux = true
  []
  [rad_avg_cost]
    type = RadialAverage
    radius = 3
    weights = linear
    prop_name = cost_sensitivity
    execute_on = TIMESTEP_END
    force_preaux = true
  []
  [update]
    type = DensityUpdateTwoConstraints
    # This is
    density_sensitivity = Dc
    cost_density_sensitivity = Cc
    cost = Cost
    cost_fraction = ${cost_frac}
    design_density = mat_den
    volume_fraction = ${vol_frac}
    bisection_lower_bound = 0
    bisection_upper_bound = 1.0e16 # 100
    bisection_move = 0.05
    adaptive_move = true
    relative_tolerance = 1.0e-3
    execute_on = TIMESTEP_BEGIN
  []
  # Provides Dc
  [calc_sense]
    type = SensitivityFilter
    density_sensitivity = Dc
    design_density = mat_den
    filter_UO = rad_avg
    execute_on = TIMESTEP_END
    force_postaux = true
  []
  # Provides Cc
  [calc_sense_cost]
    type = SensitivityFilter
    density_sensitivity = Cc
    design_density = mat_den
    filter_UO = rad_avg_cost
    execute_on = TIMESTEP_END
    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'
  nl_abs_tol = 1e-10
  dt = 1.0
  num_steps = 40
[]

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

[Postprocessors]
  [total_vol]
    type = ElementIntegralVariablePostprocessor
    variable = mat_den
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [mesh_volume]
    type = VolumePostprocessor
    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
  []
  [cost_sensitivity]
    type = ElementIntegralMaterialProperty
    mat_prop = cost_sensitivity
  []
  [cost]
    type = ElementIntegralMaterialProperty
    mat_prop = CostDensity
  []
  [cost_frac]
    type = ParsedPostprocessor
    expression = 'cost / mesh_volume'
    pp_names = 'cost mesh_volume'
  []

[]

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

vol_frac = 0.4
cost_frac = 0.3

power = 4

E0 = 1.0e-6
E1 = 0.2
E2 = 0.6
E3 = 1.0

rho0 = 1.0e-6
rho1 = 0.4
rho2 = 0.7
rho3 = 1.0

C0 = 1.0e-6
C1 = 0.5
C2 = 0.8
C3 = 1.0

[GlobalParams]
  displacements = 'disp_x disp_y'
[]

[Mesh]
  [MeshGenerator]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 50
    ny = 50
    xmin = 0
    xmax = 50
    ymin = 0
    ymax = 50
  []
  [node]
    type = ExtraNodesetGenerator
    input = MeshGenerator
    new_boundary = hold
    nodes = 0
  []
  [push_left]
    type = ExtraNodesetGenerator
    input = node
    new_boundary = push_left
    coord = '25 0 0'
  []
  [push_center]
    type = ExtraNodesetGenerator
    input = push_left
    new_boundary = push_center
    coord = '50 0 0'
  []
[]

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

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

[AuxKernels]
  [Cost]
    type = MaterialRealAux
    variable = Cost
    property = Cost_mat
  []
[]

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

[BCs]
  [no_y]
    type = DirichletBC
    variable = disp_y
    boundary = hold
    value = 0.0
  []
  [no_x_symm]
    type = DirichletBC
    variable = disp_x
    boundary = right
    value = 0.0
  []

[]
[NodalKernels]
  [push_left]
    type = NodalGravity
    variable = disp_y
    boundary = push_left
    gravity_value = -1e-3
    mass = 1
  []
  [push_center]
    type = NodalGravity
    variable = disp_y
    boundary = push_center
    gravity_value = -1e-3
    mass = 1
  []
[]
[Materials]
  [elasticity_tensor]
    type = ComputeVariableIsotropicElasticityTensor
    youngs_modulus = E_phys
    poissons_ratio = poissons_ratio
    args = 'mat_den'
  []
  [E_phys]
    type = DerivativeParsedMaterial
    # ordered multimaterial simp
    expression = "A1:=(${E0}-${E1})/(${rho0}^${power}-${rho1}^${power}); "
                 "B1:=${E0}-A1*${rho0}^${power}; E1:=A1*mat_den^${power}+B1; "
                 "A2:=(${E1}-${E2})/(${rho1}^${power}-${rho2}^${power}); "
                 "B2:=${E1}-A2*${rho1}^${power}; E2:=A2*mat_den^${power}+B2; "
                 "A3:=(${E2}-${E3})/(${rho2}^${power}-${rho3}^${power}); "
                 "B3:=${E2}-A3*${rho2}^${power}; E3:=A3*mat_den^${power}+B3; "
                 "if(mat_den<${rho1},E1,if(mat_den<${rho2},E2,E3))"
    coupled_variables = 'mat_den'
    property_name = E_phys
  []

  [Cost_mat]
    type = DerivativeParsedMaterial
    # ordered multimaterial simp
    expression = "A1:=(${C0}-${C1})/(${rho0}^(1/${power})-${rho1}^(1/${power})); "
                 "B1:=${C0}-A1*${rho0}^(1/${power}); C1:=A1*mat_den^(1/${power})+B1; "
                 "A2:=(${C1}-${C2})/(${rho1}^(1/${power})-${rho2}^(1/${power})); "
                 "B2:=${C1}-A2*${rho1}^(1/${power}); C2:=A2*mat_den^(1/${power})+B2; "
                 "A3:=(${C2}-${C3})/(${rho2}^(1/${power})-${rho3}^(1/${power})); "
                 "B3:=${C2}-A3*${rho2}^(1/${power}); C3:=A3*mat_den^(1/${power})+B3; "
                 "if(mat_den<${rho1},C1,if(mat_den<${rho2},C2,C3))"
    coupled_variables = 'mat_den'
    property_name = Cost_mat
  []
  [CostDensity]
    type = ParsedMaterial
    property_name = CostDensity
    coupled_variables = 'mat_den Cost'
    expression = 'mat_den*Cost'
  []
  [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
  []
  [cc]
    type = CostSensitivity
    design_density = mat_den
    cost = Cost_mat
    outputs = 'exodus'
  []
[]

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

[UserObjects]
  [rad_avg]
    type = RadialAverage
    radius = 3
    weights = linear
    prop_name = sensitivity
    execute_on = TIMESTEP_END
    force_preaux = true
  []
  [rad_avg_cost]
    type = RadialAverage
    radius = 3
    weights = linear
    prop_name = cost_sensitivity
    execute_on = TIMESTEP_END
    force_preaux = true
  []
  [update]
    type = DensityUpdateTwoConstraints
    # This is
    density_sensitivity = Dc
    cost_density_sensitivity = Cc
    cost = Cost
    cost_fraction = ${cost_frac}
    design_density = mat_den
    volume_fraction = ${vol_frac}
    bisection_lower_bound = 0
    bisection_upper_bound = 1.0e16 # 100
    bisection_move = 0.05
    adaptive_move = true
    relative_tolerance = 1.0e-3
    execute_on = TIMESTEP_BEGIN
  []
  # Provides Dc
  [calc_sense]
    type = SensitivityFilter
    density_sensitivity = Dc
    design_density = mat_den
    filter_UO = rad_avg
    execute_on = TIMESTEP_END
    force_postaux = true
  []
  # Provides Cc
  [calc_sense_cost]
    type = SensitivityFilter
    density_sensitivity = Cc
    design_density = mat_den
    filter_UO = rad_avg_cost
    execute_on = TIMESTEP_END
    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'
  nl_abs_tol = 1e-10
  dt = 1.0
  num_steps = 40
[]

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

[Postprocessors]
  [total_vol]
    type = ElementIntegralVariablePostprocessor
    variable = mat_den
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [mesh_volume]
    type = VolumePostprocessor
    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
  []
  [cost_sensitivity]
    type = ElementIntegralMaterialProperty
    mat_prop = cost_sensitivity
  []
  [cost]
    type = ElementIntegralMaterialProperty
    mat_prop = CostDensity
  []
  [cost_frac]
    type = ParsedPostprocessor
    expression = 'cost / mesh_volume'
    pp_names = 'cost mesh_volume'
  []

[]

References