Density Update for Multimaterial Topology Optimization

Compute updated densities based on sensitivities using an optimality criteria method to keep the volume and cost constraints satisified.

Description

The DensityUpdateTwoConstraints class is an ElementUserObject that calculates and updates the design densities based on the filtered cost and compliance sensitivities. The density update is performed using the multi-material Solid Isotropic Material with Penalization (SIMP) method (see Zuo and Saitou (2017)). The class requires the design density variable, the filtered design density variable, the density sensitivity variable, the penalty power, the volume fraction, the cost fraction, the cost variable, and the cost sensitivity variable as inputs.

Example Input File

An example of how to use the DensityUpdateTwoConstraints class in an input file:

listing test/tests/materials/compliance_sensitivity/2d_mmb_2material_cost_initial.i block=UserObjects/update

Input Parameters

costName of the cost variable.
C++ Type:VariableName
Unit:(no unit assumed)
Controllable:No
Description:Name of the cost variable.
cost_density_sensitivityName of the cost density sensitivity variable.
C++ Type:VariableName
Unit:(no unit assumed)
Controllable:No
Description:Name of the cost density sensitivity variable.
cost_fractionCost fraction.
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Cost fraction.
density_sensitivityName of the density_sensitivity variable.
C++ Type:VariableName
Unit:(no unit assumed)
Controllable:No
Description:Name of the density_sensitivity variable.
design_densityDesign density variable name.
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:Design density variable name.
volume_fractionVolume fraction.
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Volume fraction.

Required Parameters

adaptive_moveFalseWhether 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.
Default:False
C++ Type:bool
Unit:(no unit assumed)
Controllable:No
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.
bisection_lower_bound0Lower bound for the bisection algorithm.
Default:0
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Lower bound for the bisection algorithm.
bisection_move0.01Bisection move for the updated solution.
Default:0.01
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Bisection move for the updated solution.
bisection_upper_bound1e+16Upper bound for the bisection algorithm.
Default:1e+16
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Upper bound for the bisection algorithm.
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
execute_onTIMESTEP_ENDThe 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.
Default:TIMESTEP_END
C++ Type:ExecFlagEnum
Unit:(no unit assumed)
Options:XFEM_MARK, FORWARD, ADJOINT, HOMOGENEOUS_FORWARD, ADJOINT_TIMESTEP_BEGIN, ADJOINT_TIMESTEP_END, NONE, INITIAL, LINEAR, NONLINEAR_CONVERGENCE, NONLINEAR, POSTCHECK, TIMESTEP_END, TIMESTEP_BEGIN, MULTIAPP_FIXED_POINT_END, MULTIAPP_FIXED_POINT_BEGIN, FINAL, CUSTOM
Controllable:No
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.
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.
relative_tolerance0.001Relative tolerance on both compliance and cost to end the bisection method.
Default:0.001
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Relative tolerance on both compliance and cost to end the bisection method.
thermal_sensitivityName of the thermal density sensitivity variable.
C++ Type:VariableName
Unit:(no unit assumed)
Controllable:No
Description:Name of the thermal density sensitivity variable.
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.
use_thermal_complianceFalseWhether to include the thermal compliance in the sensitivities to minimize in conjunction with stiffness compliance.
Default:False
C++ Type:bool
Unit:(no unit assumed)
Controllable:No
Description:Whether to include the thermal compliance in the sensitivities to minimize in conjunction with stiffness compliance.
weight_mechanical_thermalList of values between 0 and 1 to weight the stiffness and thermal sensitivities
C++ Type:std::vector<double>
Unit:(no unit assumed)
Controllable:No
Description:List of values between 0 and 1 to weight the stiffness and thermal sensitivities

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
Unit:(no unit assumed)
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>
Unit:(no unit assumed)
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
Unit:(no unit assumed)
Controllable:Yes
Description:Set the enabled status of the MooseObject.
execution_order_group0Execution order groups are executed in increasing order (e.g., the lowest number is executed first). Note that negative group numbers may be used to execute groups before the default (0) group. Please refer to the user object documentation for ordering of user object execution within a group.
Default:0
C++ Type:int
Unit:(no unit assumed)
Controllable:No
Description:Execution order groups are executed in increasing order (e.g., the lowest number is executed first). Note that negative group numbers may be used to execute groups before the default (0) group. Please refer to the user object documentation for ordering of user object execution within a group.
force_postauxFalseForces the UserObject to be executed in POSTAUX
Default:False
C++ Type:bool
Unit:(no unit assumed)
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
Unit:(no unit assumed)
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
Unit:(no unit assumed)
Controllable:No
Description:Forces the UserObject to be executed in PREIC during initial setup
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
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/compliance_sensitivity/thermal_test.i)
(modules/combined/test/tests/optimization/compliance_sensitivity/three_materials_thermal.i)
(modules/combined/test/tests/optimization/compliance_sensitivity/2d_mmb_2material_cost_initial.i)
(modules/combined/examples/optimization/three_materials.i)
(modules/combined/test/tests/optimization/compliance_sensitivity/paper_three_materials_test.i)
(modules/optimization/test/tests/simp/2d_twoconstraints.i)

References

Wenjie Zuo and Kazuhiro Saitou. Multi-material topology optimization using ordered simp interpolation. Structural and Multidisciplinary Optimization, 55:477–491, 2017.

@article{zuo2017multi,
    author = "Zuo, Wenjie and Saitou, Kazuhiro",
    title = "Multi-material topology optimization using ordered SIMP interpolation",
    journal = "Structural and Multidisciplinary Optimization",
    volume = "55",
    pages = "477--491",
    year = "2017",
    publisher = "Springer"
}

(modules/combined/test/tests/optimization/compliance_sensitivity/thermal_test.i)

vol_frac = 0.4
cost_frac = 10.0

power = 2.0

E0 = 1.0e-6
E1 = 1.0

rho0 = 0.0
rho1 = 1.0

C0 = 1.0e-6
C1 = 1.0

TC0 = 1.0e-16
TC1 = 1.0

[GlobalParams]
  displacements = 'disp_x disp_y'
[]

[Mesh]
  [MeshGenerator]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 10
    ny = 10
    xmin = 0
    xmax = 40
    ymin = 0
    ymax = 40
  []
  [node]
    type = ExtraNodesetGenerator
    input = MeshGenerator
    new_boundary = hold
    nodes = 0
  []
  [push_left]
    type = ExtraNodesetGenerator
    input = node
    new_boundary = push_left
    coord = '16 0 0'
  []
  [push_center]
    type = ExtraNodesetGenerator
    input = push_left
    new_boundary = push_center
    coord = '24 0 0'
  []
  [extra]
    type = SideSetsFromBoundingBoxGenerator
    input = push_center
    bottom_left = '-0.01 17.999  0'
    top_right = '5 22.001  0'
    boundary_new = n1
    included_boundaries = left
  []
  [dirichlet_bc]
    type = SideSetsFromNodeSetsGenerator
    input = extra
  []
[]

[Variables]
  [disp_x]
  []
  [disp_y]
  []
  [temp]
    initial_condition = 100.0
  []
[]

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

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

[Kernels]
  [heat_conduction]
    type = HeatConduction
    variable = temp
    diffusion_coefficient = thermal_cond
  []
  [heat_source]
    type = HeatSource
    value = 1e-2 # W/m^3
    variable = temp
  []
[]

[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
  []
  [left_n1]
    type = DirichletBC
    variable = temp
    boundary = n1
    value = 0.0
  []
  [top]
    type = NeumannBC
    variable = temp
    boundary = top
    value = 0
  []
  [bottom]
    type = NeumannBC
    variable = temp
    boundary = bottom
    value = 0
  []
  [right]
    type = NeumannBC
    variable = temp
    boundary = right
    value = 0
  []
  [left]
    type = NeumannBC
    variable = temp
    boundary = left
    value = 0
  []
[]

[NodalKernels]
  [push_left]
    type = NodalGravity
    variable = disp_y
    boundary = push_left
    gravity_value = 0.0 # -1e-8
    mass = 1
  []
  [push_center]
    type = NodalGravity
    variable = disp_y
    boundary = push_center
    gravity_value = 0.0 # -1e-8
    mass = 1
  []
[]

[Materials]
  [thermal_cond]
    type = DerivativeParsedMaterial
    # ordered multimaterial simp
    expression = "A1:=(${TC0}-${TC1})/(${rho0}^${power}-${rho1}^${power}); "
                 "B1:=${TC0}-A1*${rho0}^${power}; TC1:=A1*mat_den^${power}+B1; TC1"
    coupled_variables = 'mat_den'
    property_name = thermal_cond
    outputs = 'exodus'
  []
  [thermal_compliance]
    type = ThermalCompliance
    temperature = temp
    thermal_conductivity = thermal_cond
    outputs = 'exodus'
  []
  [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; E1"
    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; C1"
    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'
  []
  [tc]
    type = ThermalSensitivity
    design_density = mat_den
    thermal_conductivity = thermal_cond
    temperature = temp
    outputs = 'exodus'
  []
[]

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

[UserObjects]
  [rad_avg]
    type = RadialAverage
    radius = 0.1
    weights = linear
    prop_name = sensitivity
    execute_on = TIMESTEP_END
    force_preaux = true
  []
  [rad_avg_cost]
    type = RadialAverage
    radius = 0.1
    weights = linear
    prop_name = cost_sensitivity
    execute_on = TIMESTEP_END
    force_preaux = true
  []
  [rad_avg_thermal]
    type = RadialAverage
    radius = 0.1
    weights = linear
    prop_name = thermal_sensitivity
    execute_on = TIMESTEP_END
    force_preaux = true
  []
  [update]
    type = DensityUpdateTwoConstraints
    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.0e12 # 100

    use_thermal_compliance = true
    thermal_sensitivity = Tc
    weight_mechanical_thermal = '0 1'

    relative_tolerance = 1.0e-12
    bisection_move = 0.015
    adaptive_move = false
    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
  []
  # Provides Tc
  [calc_sense_thermal]
    type = SensitivityFilter
    density_sensitivity = Tc
    design_density = mat_den
    filter_UO = rad_avg_thermal
    execute_on = TIMESTEP_END
    force_postaux = true
  []
[]

[Executioner]
  type = Transient
  solve_type = PJFNK
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
  petsc_options_value = 'lu superlu_dist'
  nl_abs_tol = 1e-12
  dt = 1.0
  num_steps = 5
[]

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

[Postprocessors]
  [right_flux]
    type = SideDiffusiveFluxAverage
    variable = temp
    boundary = right
    diffusivity = 10
  []
  [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
  []
  [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'
  []
  [objective]
    type = ElementIntegralMaterialProperty
    mat_prop = strain_energy_density
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [objective_thermal]
    type = ElementIntegralMaterialProperty
    mat_prop = thermal_compliance
    execute_on = 'INITIAL TIMESTEP_END'
  []
[]

(modules/combined/test/tests/optimization/compliance_sensitivity/three_materials_thermal.i)

vol_frac = 0.4
cost_frac = 0.4

power = 4

# Stiffness (not optimized in this test)
E0 = 1.0e-6
E1 = 0.2
E2 = 0.6
E3 = 1.0

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

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

# Thermal conductivity
TC0 = 1.0e-6
TC1 = 0.2
TC2 = 0.6
TC3 = 1.0

[GlobalParams]
  displacements = 'disp_x disp_y'
[]

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

[Variables]
  [disp_x]
  []
  [disp_y]
  []
  [temp]
    initial_condition = 100.0
  []
[]

[AuxVariables]
  [Dc]
    family = MONOMIAL
    order = CONSTANT
    initial_condition = -1.0
  []
  [Cc]
    family = MONOMIAL
    order = CONSTANT
    initial_condition = -1.0
  []
  [Tc]
    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}
  []
[]

# [ICs]
#   [mat_den]
#     type = RandomIC
#     seed = 4
#     variable = mat_den
#     max = '${fparse vol_frac+0.25}'
#     min = '${fparse vol_frac-0.25}'
#   []
# []

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

[Kernels]
  [heat_conduction]
    type = HeatConduction
    variable = temp
    diffusion_coefficient = thermal_cond
  []
  [heat_source]
    type = HeatSource
    value = 1e-2 # W/m^3
    variable = temp
  []
[]

[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
  []
  [top]
    type = DirichletBC
    variable = temp
    boundary = top
    value = 0
  []
  [bottom]
    type = DirichletBC
    variable = temp
    boundary = bottom
    value = 0
  []
  [right]
    type = DirichletBC
    variable = temp
    boundary = right
    value = 0
  []
  [left]
    type = DirichletBC
    variable = temp
    boundary = left
    value = 0
  []
[]

[NodalKernels]
  [push_left]
    type = NodalGravity
    variable = disp_y
    boundary = push_left
    gravity_value = -1e-6 # -3
    mass = 1
  []
  [push_center]
    type = NodalGravity
    variable = disp_y
    boundary = push_center
    gravity_value = -1e-6 # -3
    mass = 1
  []
[]

[Materials]
  [thermal_cond]
    type = DerivativeParsedMaterial
    # ordered multimaterial simp
    expression = "A1:=(${TC0}-${TC1})/(${rho0}^${power}-${rho1}^${power}); "
                 "B1:=${TC0}-A1*${rho0}^${power}; TC1:=A1*mat_den^${power}+B1; "
                 "A2:=(${TC1}-${TC2})/(${rho1}^${power}-${rho2}^${power}); "
                 "B2:=${TC1}-A2*${rho1}^${power}; TC2:=A2*mat_den^${power}+B2; "
                 "A3:=(${TC2}-${TC3})/(${rho2}^${power}-${rho3}^${power}); "
                 "B3:=${TC2}-A3*${rho2}^${power}; TC3:=A3*mat_den^${power}+B3; "
                 "if(mat_den<${rho1},TC1,if(mat_den<${rho2},TC2,TC3))"
    coupled_variables = 'mat_den'
    property_name = thermal_cond
    outputs = 'exodus'
  []
  [thermal_compliance]
    type = ThermalCompliance
    temperature = temp
    thermal_conductivity = thermal_cond
    outputs = 'exodus'
  []
  [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'
  []
  [tc]
    type = ThermalSensitivity
    design_density = mat_den
    thermal_conductivity = thermal_cond
    temperature = temp
    outputs = 'exodus'
  []
[]

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

[UserObjects]
  [rad_avg]
    type = RadialAverage
    radius = 4
    weights = linear
    prop_name = sensitivity
    execute_on = TIMESTEP_END
    force_preaux = true
  []
  [rad_avg_cost]
    type = RadialAverage
    radius = 4
    weights = linear
    prop_name = cost_sensitivity
    execute_on = TIMESTEP_END
    force_preaux = true
  []
  [rad_avg_thermal]
    type = RadialAverage
    radius = 4
    weights = linear
    prop_name = thermal_sensitivity
    execute_on = TIMESTEP_END
    force_preaux = true
  []
  [update]
    type = DensityUpdateTwoConstraints
    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.0e12 # 100

    use_thermal_compliance = true
    thermal_sensitivity = Tc
    # Only account for thermal optimizxation
    weight_mechanical_thermal = '0 1'

    relative_tolerance = 1.0e-8
    bisection_move = 0.05
    adaptive_move = false
    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
  []
  # Provides Tc
  [calc_sense_thermal]
    type = SensitivityFilter
    density_sensitivity = Tc
    design_density = mat_den
    filter_UO = rad_avg_thermal
    execute_on = TIMESTEP_END
    force_postaux = true
  []
[]

[Executioner]
  type = Transient
  solve_type = PJFNK
  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 = 12
[]

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

[Postprocessors]
  [right_flux]
    type = SideDiffusiveFluxAverage
    variable = temp
    boundary = right
    diffusivity = 10
  []
  [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
  []
  [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'
  []
  [objective]
    type = ElementIntegralMaterialProperty
    mat_prop = strain_energy_density
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [objective_thermal]
    type = ElementIntegralMaterialProperty
    mat_prop = thermal_compliance
    execute_on = 'INITIAL TIMESTEP_END'
  []
[]

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

vol_frac = 0.4
cost_frac = 0.22 # Change back to 0.4

power = 2.0

E0 = 1.0e-6
E1 = 0.3
E2 = 1.0

rho0 = 1.0e-6
rho1 = 0.3
rho2 = 1.0

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

[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
    nodes = 0
  []
  [push]
    type = ExtraNodesetGenerator
    input = node
    new_boundary = push
    coord = '30 10 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_x]
    type = DirichletBC
    variable = disp_y
    boundary = hold
    value = 0.0
  []
  [no_y]
    type = DirichletBC
    variable = disp_x
    boundary = right
    value = 0.0
  []

[]
[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
    # 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; "
                 "if(mat_den<${rho1},E1,E2)"
    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; "
                 "if(mat_den<${rho1},C1,C2)"
    coupled_variables = 'mat_den'
    property_name = Cost_mat
  []

  [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 = 1.2
    weights = linear
    prop_name = sensitivity
    execute_on = TIMESTEP_END
    force_preaux = true
  []
  [rad_avg_cost]
    type = RadialAverage
    radius = 1.2
    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

    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-8
  dt = 1.0
  num_steps = 10 #50000
[]

[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'
  []
  [sensitivity]
    type = ElementIntegralMaterialProperty
    mat_prop = sensitivity
  []
  [cost_sensitivity]
    type = ElementIntegralMaterialProperty
    mat_prop = cost_sensitivity
  []
  [cost]
    type = ElementIntegralVariablePostprocessor
    variable = Cost
  []
[]

(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/test/tests/optimization/compliance_sensitivity/paper_three_materials_test.i)

vol_frac = 0.4
cost_frac = 0.2 #0.283 # Change back to 0.4

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 = 2
    weights = linear
    prop_name = sensitivity
    execute_on = TIMESTEP_END
    force_preaux = true
  []
  [rad_avg_cost]
    type = RadialAverage
    radius = 2
    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

    relative_tolerance = 1.0e-3
    bisection_move = 0.02
    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 = 25
[]

[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
  []
  [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'
  []
  [objective]
    type = ElementIntegralMaterialProperty
    mat_prop = strain_energy_density
    execute_on = 'INITIAL TIMESTEP_END'
  []
[]

(modules/optimization/test/tests/simp/2d_twoconstraints.i)

cost_frac = 0.3
vol_frac = 0.2

[Mesh]
  [planet]
    type = ConcentricCircleMeshGenerator
    has_outer_square = false
    radii = 1
    num_sectors = 10
    rings = 2
    preserve_volumes = false
  []

  [moon]
    type = ConcentricCircleMeshGenerator
    has_outer_square = false
    radii = 0.5
    num_sectors = 8
    rings = 2
    preserve_volumes = false
  []
  [combine]
    type = CombinerGenerator
    inputs = 'planet moon'
    positions = '0 0 0 -1.5 -0.5 0'
  []
[]

[AuxVariables]
  [mat_den]
    family = MONOMIAL
    order = CONSTANT
    initial_condition = 0.1
  []
  [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
  []
[]

[Variables]
  [u]
  []
  [v]
  []
[]

[Kernels]
  [diff_u]
    type = Diffusion
    variable = u
  []
  [dt_u]
    type = TimeDerivative
    variable = u
  []
  [diff_v]
    type = Diffusion
    variable = v
  []
  [dt_v]
    type = TimeDerivative
    variable = v
  []
[]

[Materials]
  [thermal_cond]
    type = GenericFunctionMaterial
    prop_values = '-1.4*abs(y)-2.7*abs(x)'
    prop_names = thermal_cond
    outputs = 'exodus'
  []
  [thermal_compliance_sensitivity]
    type = GenericFunctionMaterial
    prop_values = '-3*abs(y)-1.5*abs(x)'
    prop_names = thermal_sensitivity
    outputs = 'exodus'
  []
  [cost_sensitivity]
    type = GenericFunctionMaterial
    prop_values = '-0.3*y*y-0.5*abs(x*y)'
    prop_names = cost_sensitivity
    outputs = 'exodus'
  []
  [cost_sensitivity_parsed]
    type = DerivativeParsedMaterial
    expression = "if(mat_den<0.2,1.0,0.5)"
    coupled_variables = 'mat_den'
    property_name = cost_sensitivity_parsed
  []
  [cc]
    type = CostSensitivity
    design_density = mat_den
    cost = cost_sensitivity_parsed
    outputs = 'exodus'
    declare_suffix = 'for_testing'
  []
[]

[BCs]
  [flux_u]
    type = DirichletBC
    variable = u
    boundary = outer
    value = 3.0
  []

  [flux_v]
    type = DirichletBC
    variable = v
    boundary = outer
    value = 7.0
  []
[]

[UserObjects]
  [rad_avg]
    type = RadialAverage
    radius = 0.1
    weights = linear
    prop_name = thermal_sensitivity
    execute_on = TIMESTEP_END
    force_preaux = true
  []
  [rad_avg_cost]
    type = RadialAverage
    radius = 1.2
    weights = linear
    prop_name = cost_sensitivity
    execute_on = TIMESTEP_END
    force_preaux = true
  []
  [update]
    type = DensityUpdateTwoConstraints
    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
    relative_tolerance = 1.0e-3
    bisection_move = 0.15
    execute_on = TIMESTEP_BEGIN
  []
  [calc_sense]
    type = SensitivityFilter
    density_sensitivity = Dc
    design_density = mat_den
    filter_UO = rad_avg
    execute_on = TIMESTEP_END
    force_postaux = true
  []
  [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
  dt = 0.1
  num_steps = 3
  nl_rel_tol = 1e-04
[]

[Outputs]
  exodus = true
[]

Description
Example Input File
Input Parameters
Input Files
References