- argsVariable dependence for the Young's Modulus and Poisson's Ratio materials
C++ Type:std::vector<VariableName>
Controllable:No
Description:Variable dependence for the Young's Modulus and Poisson's Ratio materials
- poissons_ratioName of material property defining the Poisson's Ratio
C++ Type:MaterialPropertyName
Controllable:No
Description:Name of material property defining the Poisson's Ratio
- youngs_modulusName of material property defining the Young's Modulus
C++ Type:MaterialPropertyName
Controllable:No
Description:Name of material property defining the Young's Modulus
Compute Variable Isotropic Elasticity Tensor
Compute an isotropic elasticity tensor for elastic constants that change as a function of material properties
Description
This model computes an elasticity tensor for which the elastic moduli, prescribed in terms of the Young's modulus and Poisson's ratio, vary as defined by user-specified material properties. Every time this material is evaluated, the full tensor is updated to reflect the current values of those elastic constants.
Example Input File Syntax
[./elasticity_tensor]
type = ComputeVariableIsotropicElasticityTensor
args = temp
youngs_modulus = youngs_modulus
poissons_ratio = poissons_ratio
[../]
(modules/combined/test/tests/thermal_elastic/thermal_elastic.i)Input Parameters
- base_nameOptional parameter that allows the user to define multiple mechanics material systems on the same block, i.e. for multiple phases
C++ Type:std::string
Controllable:No
Description:Optional parameter that allows the user to define multiple mechanics material systems on the same block, i.e. for multiple phases
- blockThe list of blocks (ids or names) that this object will be applied
C++ Type:std::vector<SubdomainName>
Controllable:No
Description:The list of blocks (ids or names) that this object will be applied
- boundaryThe list of boundaries (ids or names) from the mesh where this object applies
C++ Type:std::vector<BoundaryName>
Controllable:No
Description:The list of boundaries (ids or names) from the mesh where this object applies
- computeTrueWhen false, MOOSE will not call compute methods on this material. The user must call computeProperties() after retrieving the MaterialBase via MaterialBasePropertyInterface::getMaterialBase(). Non-computed MaterialBases are not sorted for dependencies.
Default:True
C++ Type:bool
Controllable:No
Description:When false, MOOSE will not call compute methods on this material. The user must call computeProperties() after retrieving the MaterialBase via MaterialBasePropertyInterface::getMaterialBase(). Non-computed MaterialBases are not sorted for dependencies.
- constant_onNONEWhen ELEMENT, MOOSE will only call computeQpProperties() for the 0th quadrature point, and then copy that value to the other qps.When SUBDOMAIN, MOOSE will only call computeQpProperties() for the 0th quadrature point, and then copy that value to the other qps. Evaluations on element qps will be skipped
Default:NONE
C++ Type:MooseEnum
Options:NONE, ELEMENT, SUBDOMAIN
Controllable:No
Description:When ELEMENT, MOOSE will only call computeQpProperties() for the 0th quadrature point, and then copy that value to the other qps.When SUBDOMAIN, MOOSE will only call computeQpProperties() for the 0th quadrature point, and then copy that value to the other qps. Evaluations on element qps will be skipped
- declare_suffixAn optional suffix parameter that can be appended to any declared properties. The suffix will be prepended with a '_' character.
C++ Type:MaterialPropertyName
Controllable:No
Description:An optional suffix parameter that can be appended to any declared properties. The suffix will be prepended with a '_' character.
- elasticity_tensor_prefactorOptional function to use as a scalar prefactor on the elasticity tensor.
C++ Type:FunctionName
Controllable:No
Description:Optional function to use as a scalar prefactor on the elasticity tensor.
- 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
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
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.
Optional Parameters
- control_tagsAdds user-defined labels for accessing object parameters via control logic.
C++ Type:std::vector<std::string>
Controllable:No
Description:Adds user-defined labels for accessing object parameters via control logic.
- enableTrueSet the enabled status of the MooseObject.
Default:True
C++ Type:bool
Controllable:Yes
Description:Set the enabled status of the MooseObject.
- implicitTrueDetermines whether this object is calculated using an implicit or explicit form
Default:True
C++ Type:bool
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
Controllable:No
Description:The seed for the master random number generator
- use_displaced_meshFalseWhether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used.
Default:False
C++ Type:bool
Controllable:No
Description:Whether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used.
Advanced Parameters
- output_propertiesList of material properties, from this material, to output (outputs must also be defined to an output type)
C++ Type:std::vector<std::string>
Controllable:No
Description:List of material properties, from this material, to output (outputs must also be defined to an output type)
- outputsnone Vector of output names where you would like to restrict the output of variables(s) associated with this object
Default:none
C++ Type:std::vector<OutputName>
Controllable:No
Description:Vector of output names where you would like to restrict the output of variables(s) associated with this object
Outputs Parameters
Input Files
- (modules/combined/examples/phase_field-mechanics/interface_stress.i)
- (modules/combined/test/tests/optimization/compliance_sensitivity/2d_mmb_2material.i)
- (modules/combined/test/tests/optimization/compliance_sensitivity/2d_mbb.i)
- (modules/combined/examples/optimization/multi-load/single_subapp_one.i)
- (modules/combined/test/tests/thermo_mech/youngs_modulus_function_temp.i)
- (modules/combined/examples/optimization/three_materials.i)
- (modules/combined/test/tests/thermal_elastic/thermal_elastic.i)
- (modules/combined/examples/optimization/thermomechanical/thermal_sub.i)
- (modules/combined/test/tests/optimization/optimization_density_update/top_opt_3d_pde_filter.i)
- (modules/solid_mechanics/test/tests/j_integral_vtest/fgm_5.i)
- (modules/combined/test/tests/optimization/compliance_sensitivity/2d_mmb_2material_cost_initial.i)
- (modules/combined/test/tests/optimization/optimization_density_update/top_opt_2d_pde_filter.i)
- (modules/combined/examples/optimization/multi-load/square_subapp_two.i)
- (modules/combined/examples/optimization/multi-load/single_main.i)
- (modules/solid_mechanics/test/tests/j_integral_vtest/j_int_fgm_sif.i)
- (modules/combined/test/tests/optimization/compliance_sensitivity/2d_mbb_pde.i)
- (modules/combined/test/tests/optimization/compliance_sensitivity/2d_mbb_pde_amr.i)
- (modules/combined/examples/optimization/multi-load/square_main.i)
- (modules/combined/test/tests/optimization/optimization_density_update/top_opt_3d.i)
- (modules/combined/examples/optimization/2d_mbb_pde_amr.i)
- (modules/combined/test/tests/optimization/optimization_density_update/top_opt_2d.i)
- (modules/combined/examples/optimization/2d_mbb.i)
- (modules/combined/examples/optimization/multi-load/square_subapp_one.i)
- (modules/combined/examples/optimization/3d_mbb.i)
- (modules/combined/examples/optimization/2d_mbb_pde.i)
- (modules/combined/test/tests/thermal_elastic/derivatives.i)
- (modules/combined/examples/optimization/thermomechanical/thermomechanical_main.i)
- (modules/combined/examples/optimization/thermomechanical/structural_sub.i)
- (modules/combined/test/tests/optimization/compliance_sensitivity/2d_mmb_2material_cost.i)
- (modules/combined/examples/optimization/helmholtz_multimat_nostrip.i)
- (modules/combined/test/tests/optimization/compliance_sensitivity/three_materials_thermal.i)
- (modules/combined/examples/optimization/helmholtz_multimat_strip.i)
- (modules/combined/test/tests/optimization/compliance_sensitivity/thermal_test.i)
- (modules/combined/test/tests/optimization/compliance_sensitivity/3d_mbb.i)
- (modules/combined/test/tests/optimization/compliance_sensitivity/paper_three_materials_test.i)
- (modules/combined/examples/optimization/multi-load/single_subapp_two.i)
(modules/combined/test/tests/thermal_elastic/thermal_elastic.i)
# Patch Test
# This test is designed to compute constant xx, yy, zz, xy, yz, and xz
# stress on a set of irregular hexes. The mesh is composed of one
# block with seven elements. The elements form a unit cube with one
# internal element. There is a nodeset for each exterior node.
# The cube is displaced by 1e-6 units in x, 2e-6 in y, and 3e-6 in z.
# The faces are sheared as well (1e-6, 2e-6, and 3e-6 for xy, yz, and
# zx). This gives a uniform strain/stress state for all six unique
# tensor components. This displacement is again applied in the second
# step.
# With Young's modulus at 1e6 and Poisson's ratio at 0, the shear
# modulus is 5e5 (G=E/2/(1+nu)). Therefore, for the mechanical strain,
#
# stress xx = 1e6 * 1e-6 = 1
# stress yy = 1e6 * 2e-6 = 2
# stress zz = 1e6 * 3e-6 = 3
# stress xy = 2 * 5e5 * 1e-6 / 2 = 0.5
# (2 * G * gamma_xy / 2 = 2 * G * epsilon_xy)
# stress yz = 2 * 5e5 * 2e-6 / 2 = 1
# stress zx = 2 * 5e5 * 3e-6 / 2 = 1.5
# Young's modulus is a function of temperature for this test. The
# temperature changes from 100 to 500. The Young's modulus drops
# due to that temperature change from 1e6 to 6e5.
# Poisson's ratio also is a function of temperature and changes from
# 0 to 0.25.
# At the end of the temperature ramp, E=6e5 and nu=0.25. This gives
# G=2.4e=5. lambda=E*nu/(1+nu)/(1-2*nu)=2.4E5. The final stress
# is therefore
# stress xx = 2.4e5 * 12e-6 + 2*2.4e5*2e-6 = 3.84
# stress yy = 2.4e5 * 12e-6 + 2*2.4e5*4e-6 = 4.80
# stress zz = 2.4e5 * 12e-6 + 2*2.4e5*6e-6 = 5.76
# stress xy = 2 * 2.4e5 * 2e-6 / 2 = 0.48
# (2 * G * gamma_xy / 2 = 2 * G * epsilon_xy)
# stress yz = 2 * 2.4e5 * 4e-6 / 2 = 0.96
# stress xz = 2 * 2.4e5 * 6e-6 / 2 = 1.44
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
[]
[Mesh]
file = thermal_elastic.e
[]
[Functions]
[./ramp1]
type = PiecewiseLinear
x = '0. 1. 2.'
y = '0. 1. 2.'
scale_factor = 1e-6
[../]
[./ramp2]
type = PiecewiseLinear
x = '0. 1. 2.'
y = '0. 1. 2.'
scale_factor = 2e-6
[../]
[./ramp3]
type = PiecewiseLinear
x = '0. 1. 2.'
y = '0. 1. 2.'
scale_factor = 3e-6
[../]
[./ramp4]
type = PiecewiseLinear
x = '0. 1. 2.'
y = '0. 1. 2.'
scale_factor = 4e-6
[../]
[./ramp6]
type = PiecewiseLinear
x = '0. 1. 2.'
y = '0. 1. 2.'
scale_factor = 6e-6
[../]
[./tempFunc]
type = PiecewiseLinear
x = '0 1 2'
y = '100.0 100.0 500.0'
[../]
[]
[Variables]
[./temp]
initial_condition = 100.0
[../]
[]
[Modules/TensorMechanics/Master]
[./all]
add_variables = true
generate_output = 'stress_xx stress_yy stress_zz stress_xy stress_xz stress_yz'
strain = FINITE
[../]
[]
[Kernels]
[./heat]
type = Diffusion
variable = temp
[../]
[]
[BCs]
[./node1_x]
type = DirichletBC
variable = disp_x
boundary = 1
value = 0.0
[../]
[./node1_y]
type = FunctionDirichletBC
variable = disp_y
boundary = 1
function = ramp2
[../]
[./node1_z]
type = FunctionDirichletBC
variable = disp_z
boundary = 1
function = ramp3
[../]
[./node2_x]
type = FunctionDirichletBC
variable = disp_x
boundary = 2
function = ramp1
[../]
[./node2_y]
type = FunctionDirichletBC
variable = disp_y
boundary = 2
function = ramp2
[../]
[./node2_z]
type = FunctionDirichletBC
variable = disp_z
boundary = 2
function = ramp6
[../]
[./node3_x]
type = FunctionDirichletBC
variable = disp_x
boundary = 3
function = ramp1
[../]
[./node3_y]
type = DirichletBC
variable = disp_y
boundary = 3
value = 0.0
[../]
[./node3_z]
type = FunctionDirichletBC
variable = disp_z
boundary = 3
function = ramp3
[../]
[./node4_x]
type = DirichletBC
variable = disp_x
boundary = 4
value = 0.0
[../]
[./node4_y]
type = DirichletBC
variable = disp_y
boundary = 4
value = 0.0
[../]
[./node4_z]
type = DirichletBC
variable = disp_z
boundary = 4
value = 0.0
[../]
[./node5_x]
type = FunctionDirichletBC
variable = disp_x
boundary = 5
function = ramp1
[../]
[./node5_y]
type = FunctionDirichletBC
variable = disp_y
boundary = 5
function = ramp4
[../]
[./node5_z]
type = FunctionDirichletBC
variable = disp_z
boundary = 5
function = ramp3
[../]
[./node6_x]
type = FunctionDirichletBC
variable = disp_x
boundary = 6
function = ramp2
[../]
[./node6_y]
type = FunctionDirichletBC
variable = disp_y
boundary = 6
function = ramp4
[../]
[./node6_z]
type = FunctionDirichletBC
variable = disp_z
boundary = 6
function = ramp6
[../]
[./node7_x]
type = FunctionDirichletBC
variable = disp_x
boundary = 7
function = ramp2
[../]
[./node7_y]
type = FunctionDirichletBC
variable = disp_y
boundary = 7
function = ramp2
[../]
[./node7_z]
type = FunctionDirichletBC
variable = disp_z
boundary = 7
function = ramp3
[../]
[./node8_x]
type = FunctionDirichletBC
variable = disp_x
boundary = 8
function = ramp1
[../]
[./node8_y]
type = FunctionDirichletBC
variable = disp_y
boundary = 8
function = ramp2
[../]
[./node8_z]
type = DirichletBC
variable = disp_z
boundary = 8
value = 0.0
[../]
[./temp]
type = FunctionDirichletBC
variable = temp
boundary = '10 12'
function = tempFunc
[../]
[]
[Materials]
[./youngs_modulus]
type = PiecewiseLinearInterpolationMaterial
x = '100 500'
y = '1e6 6e5'
property = youngs_modulus
variable = temp
[../]
[./poissons_ratio]
type = PiecewiseLinearInterpolationMaterial
x = '100 500'
y = '0 0.25'
property = poissons_ratio
variable = temp
[../]
[./elasticity_tensor]
type = ComputeVariableIsotropicElasticityTensor
args = temp
youngs_modulus = youngs_modulus
poissons_ratio = poissons_ratio
[../]
[./stress]
type = ComputeFiniteStrainElasticStress
[../]
[]
[Executioner]
type = Transient
solve_type = 'PJFNK'
nl_rel_tol = 1e-9
nl_abs_tol = 1e-9
l_max_its = 20
start_time = 0.0
dt = 1.0
end_time = 2.0
[]
[Outputs]
exodus = true
[]
(modules/combined/examples/phase_field-mechanics/interface_stress.i)
[Mesh]
type = GeneratedMesh
dim = 3
nx = 50
ny = 50
nz = 50
xmax = 10
ymax = 10
zmax = 10
xmin = -10
ymin = -10
zmin = -10
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
[]
[Functions]
[./sphere]
type = ParsedFunction
expression = 'r:=sqrt(x^2+y^2+z^2); R:=(4.0-r)/2.0; if(R>1,1,if(R<0,0,3*R^2-2*R^3))'
[../]
[]
[AuxVariables]
[./eta]
[./InitialCondition]
type = FunctionIC
function = sphere
[../]
[../]
[]
[Modules/TensorMechanics/Master]
[./all]
add_variables = true
generate_output = 'hydrostatic_stress stress_xx'
[../]
[]
[Materials]
[./ym]
type = DerivativeParsedMaterial
property_name = ym
expression = (1-eta)*7+0.5
coupled_variables = eta
[../]
[./elasticity]
type = ComputeVariableIsotropicElasticityTensor
poissons_ratio = 0.45
youngs_modulus = ym
args = eta
[../]
[./stress]
type = ComputeLinearElasticStress
[../]
[./interface]
type = ComputeInterfaceStress
v = eta
stress = 1.0
[../]
[]
[VectorPostprocessors]
[./line]
type = SphericalAverage
variable = 'hydrostatic_stress'
radius = 10
bin_number = 40
[../]
[]
[Executioner]
type = Steady
[]
[Outputs]
exodus = true
csv = true
[]
(modules/combined/test/tests/optimization/compliance_sensitivity/2d_mmb_2material.i)
vol_frac = 0.5
power = 1
E0 = 1e-5
E1 = 0.6
E2 = 1.0
rho0 = 0.0
rho1 = 0.4
rho2 = 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
[]
[mat_den]
family = MONOMIAL
order = CONSTANT
# initial_condition = ${vol_frac}
[]
[]
[ICs]
[mat_den]
type = RandomIC
seed = 5
variable = mat_den
max = '${fparse vol_frac+0.15}'
min = '${fparse vol_frac-0.15}'
[]
[]
[Modules/TensorMechanics/Master]
[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
[]
[poissons_ratio]
type = GenericConstantMaterial
prop_names = poissons_ratio
prop_values = 0.3
[]
[stress]
type = ComputeLinearElasticStress
[]
[dc]
type = ComplianceSensitivity2
design_density = mat_den
youngs_modulus = E_phys
[]
[]
[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
[]
[update]
type = DensityUpdate
density_sensitivity = Dc
design_density = mat_den
volume_fraction = ${vol_frac}
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
[]
[]
[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 = 70
[]
[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
[]
[]
(modules/combined/test/tests/optimization/compliance_sensitivity/2d_mbb.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
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
[]
[mat_den]
family = MONOMIAL
order = CONSTANT
initial_condition = ${vol_frac}
[]
[]
[Modules/TensorMechanics/Master]
[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
# 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]
[rad_avg]
type = RadialAverage
radius = 1.2
weights = linear
prop_name = sensitivity
execute_on = TIMESTEP_END
force_preaux = true
[]
[update]
type = DensityUpdate
density_sensitivity = Dc
design_density = mat_den
volume_fraction = ${vol_frac}
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
[]
[]
[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 = 70
[]
[Outputs]
[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
[]
[]
(modules/combined/examples/optimization/multi-load/single_subapp_one.i)
power = 2
E0 = 1.0
Emin = 1.0e-6
[GlobalParams]
displacements = 'disp_x disp_y'
[]
[Mesh]
# final_generator = 'MoveRight'
[Bottom]
type = GeneratedMeshGenerator
dim = 2
nx = 80
ny = 40
xmin = 0
xmax = 150
ymin = 0
ymax = 75
[]
[left_load]
type = ExtraNodesetGenerator
input = Bottom
new_boundary = left_load
coord = '37.5 75 0'
[]
[right_load]
type = ExtraNodesetGenerator
input = left_load
new_boundary = right_load
coord = '112.5 75 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]
[]
[]
[AuxVariables]
[Dc]
family = MONOMIAL
order = SECOND
initial_condition = -1.0
[]
[Cc]
family = MONOMIAL
order = CONSTANT
initial_condition = -1.0
[]
[mat_den]
family = MONOMIAL
order = CONSTANT
initial_condition = 0.1
[]
[sensitivity_var]
family = MONOMIAL
order = SECOND
initial_condition = -1.0
[]
[]
[AuxKernels]
[sensitivity_kernel]
type = MaterialRealAux
property = sensitivity
variable = sensitivity_var
check_boundary_restricted = false
execute_on = 'TIMESTEP_END'
[]
[]
[Modules/TensorMechanics/Master]
[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
[]
[]
[NodalKernels]
[push_left]
type = NodalGravity
variable = disp_y
boundary = left_load
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
# 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.0
[]
[stress]
type = ComputeLinearElasticStress
[]
[dc]
type = ComplianceSensitivity
design_density = mat_den
youngs_modulus = E_phys
[]
[]
[Preconditioning]
[smp]
type = SMP
full = true
[]
[]
[UserObjects]
[rad_avg]
type = RadialAverage
radius = 3
weights = linear
prop_name = sensitivity
force_preaux = true
execute_on = 'TIMESTEP_END'
[]
# No SIMP optimization in subapp
[calc_sense]
type = SensitivityFilter
density_sensitivity = Dc
design_density = mat_den
filter_UO = rad_avg
force_postaux = true
execute_on = 'TIMESTEP_END'
[]
[]
[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
function = 'total_vol / mesh_volume'
pp_names = 'total_vol mesh_volume'
[]
[sensitivity]
type = ElementIntegralMaterialProperty
mat_prop = sensitivity
execute_on = 'TIMESTEP_BEGIN TIMESTEP_END NONLINEAR'
[]
[objective]
type = ElementIntegralMaterialProperty
mat_prop = strain_energy_density
execute_on = 'INITIAL TIMESTEP_END'
[]
[]
(modules/combined/test/tests/thermo_mech/youngs_modulus_function_temp.i)
# ---------------------------------------------------------------------------
# This test is designed to verify the variable elasticity tensor functionality in the
# ComputeFiniteStrainElasticStress class with the elasticity_tensor_has_changed flag
# by varying the young's modulus with temperature. A constant strain is applied
# to the mesh in this case, and the stress varies with the changing elastic constants.
#
# Geometry: A single element cube in symmetry boundary conditions and pulled
# at a constant displacement to create a constant strain in the x-direction.
#
# Temperature: The temperature varies from 400K to 700K in this simulation by
# 100K each time step. The temperature is held constant in the last
# timestep to ensure that the elasticity tensor components are constant
# under constant temperature.
#
# Results: Because Poisson's ratio is set to zero, only the stress along the x
# axis is non-zero. The stress changes with temperature.
#
# Temperature(K) strain_{xx}(m/m) Young's Modulus(Pa) stress_{xx}(Pa)
# 400 0.001 10.0e6 1.0e4
# 500 0.001 10.0e6 1.0e4
# 600 0.001 9.94e6 9.94e3
# 700 0.001 9.93e6 9.93e3
#
# The tensor mechanics results align exactly with the analytical results above
# when this test is run with ComputeIncrementalSmallStrain. When the test is
# run with ComputeFiniteStrain, a 0.05% discrepancy between the analytical
# strains and the simulation strain results is observed, and this discrepancy
# is carried over into the calculation of the elastic stress.
#-------------------------------------------------------------------------
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
[]
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 1
nz = 1
[]
[Variables]
[./disp_x]
[../]
[./disp_y]
[../]
[./disp_z]
[../]
[./temp]
initial_condition = 400
[../]
[]
[AuxVariables]
[./stress_xx]
order = CONSTANT
family = MONOMIAL
[../]
[./elastic_strain_xx]
order = CONSTANT
family = MONOMIAL
[../]
[]
[Functions]
[./temperature_function]
type = PiecewiseLinear
x = '1 4'
y = '400 700'
[../]
[]
[Kernels]
[./heat]
type = Diffusion
variable = temp
[../]
[./TensorMechanics]
use_displaced_mesh = true
[../]
[]
[AuxKernels]
[./stress_xx]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_xx
index_i = 0
index_j = 0
[../]
[./elastic_strain_xx]
type = RankTwoAux
rank_two_tensor = elastic_strain
variable = elastic_strain_xx
index_i = 0
index_j = 0
execute_on = timestep_end
[../]
[]
[BCs]
[./u_left_fix]
type = DirichletBC
variable = disp_x
boundary = left
value = 0.0
[../]
[./u_bottom_fix]
type = DirichletBC
variable = disp_y
boundary = bottom
value = 0.0
[../]
[./u_back_fix]
type = DirichletBC
variable = disp_z
boundary = back
value = 0.0
[../]
[./u_pull_right]
type = DirichletBC
variable = disp_x
boundary = right
value = 0.001
[../]
[./temp_bc_1]
type = FunctionDirichletBC
variable = temp
preset = false
boundary = '1 2 3 4'
function = temperature_function
[../]
[]
[Materials]
[./youngs_modulus]
type = PiecewiseLinearInterpolationMaterial
xy_data = '0 10e+6
599.9999 10e+6
600 9.94e+6
99900 10e3'
property = youngs_modulus
variable = temp
[../]
[./elasticity_tensor]
type = ComputeVariableIsotropicElasticityTensor
args = temp
youngs_modulus = youngs_modulus
poissons_ratio = 0.0
[../]
[./strain]
type = ComputeIncrementalSmallStrain
[../]
[./stress]
type = ComputeFiniteStrainElasticStress
[../]
[]
[Executioner]
type = Transient
end_time = 5
[]
[Postprocessors]
[./elastic_strain_xx]
type = ElementAverageValue
variable = elastic_strain_xx
[../]
[./elastic_stress_xx]
type = ElementAverageValue
variable = stress_xx
[../]
[./temp]
type = AverageNodalVariableValue
variable = temp
[../]
[]
[Outputs]
exodus = true
[]
(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
[]
[]
[Modules/TensorMechanics/Master]
[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
function = '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
function = 'cost / mesh_volume'
pp_names = 'cost mesh_volume'
[]
[]
(modules/combined/test/tests/thermal_elastic/thermal_elastic.i)
# Patch Test
# This test is designed to compute constant xx, yy, zz, xy, yz, and xz
# stress on a set of irregular hexes. The mesh is composed of one
# block with seven elements. The elements form a unit cube with one
# internal element. There is a nodeset for each exterior node.
# The cube is displaced by 1e-6 units in x, 2e-6 in y, and 3e-6 in z.
# The faces are sheared as well (1e-6, 2e-6, and 3e-6 for xy, yz, and
# zx). This gives a uniform strain/stress state for all six unique
# tensor components. This displacement is again applied in the second
# step.
# With Young's modulus at 1e6 and Poisson's ratio at 0, the shear
# modulus is 5e5 (G=E/2/(1+nu)). Therefore, for the mechanical strain,
#
# stress xx = 1e6 * 1e-6 = 1
# stress yy = 1e6 * 2e-6 = 2
# stress zz = 1e6 * 3e-6 = 3
# stress xy = 2 * 5e5 * 1e-6 / 2 = 0.5
# (2 * G * gamma_xy / 2 = 2 * G * epsilon_xy)
# stress yz = 2 * 5e5 * 2e-6 / 2 = 1
# stress zx = 2 * 5e5 * 3e-6 / 2 = 1.5
# Young's modulus is a function of temperature for this test. The
# temperature changes from 100 to 500. The Young's modulus drops
# due to that temperature change from 1e6 to 6e5.
# Poisson's ratio also is a function of temperature and changes from
# 0 to 0.25.
# At the end of the temperature ramp, E=6e5 and nu=0.25. This gives
# G=2.4e=5. lambda=E*nu/(1+nu)/(1-2*nu)=2.4E5. The final stress
# is therefore
# stress xx = 2.4e5 * 12e-6 + 2*2.4e5*2e-6 = 3.84
# stress yy = 2.4e5 * 12e-6 + 2*2.4e5*4e-6 = 4.80
# stress zz = 2.4e5 * 12e-6 + 2*2.4e5*6e-6 = 5.76
# stress xy = 2 * 2.4e5 * 2e-6 / 2 = 0.48
# (2 * G * gamma_xy / 2 = 2 * G * epsilon_xy)
# stress yz = 2 * 2.4e5 * 4e-6 / 2 = 0.96
# stress xz = 2 * 2.4e5 * 6e-6 / 2 = 1.44
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
[]
[Mesh]
file = thermal_elastic.e
[]
[Functions]
[./ramp1]
type = PiecewiseLinear
x = '0. 1. 2.'
y = '0. 1. 2.'
scale_factor = 1e-6
[../]
[./ramp2]
type = PiecewiseLinear
x = '0. 1. 2.'
y = '0. 1. 2.'
scale_factor = 2e-6
[../]
[./ramp3]
type = PiecewiseLinear
x = '0. 1. 2.'
y = '0. 1. 2.'
scale_factor = 3e-6
[../]
[./ramp4]
type = PiecewiseLinear
x = '0. 1. 2.'
y = '0. 1. 2.'
scale_factor = 4e-6
[../]
[./ramp6]
type = PiecewiseLinear
x = '0. 1. 2.'
y = '0. 1. 2.'
scale_factor = 6e-6
[../]
[./tempFunc]
type = PiecewiseLinear
x = '0 1 2'
y = '100.0 100.0 500.0'
[../]
[]
[Variables]
[./temp]
initial_condition = 100.0
[../]
[]
[Modules/TensorMechanics/Master]
[./all]
add_variables = true
generate_output = 'stress_xx stress_yy stress_zz stress_xy stress_xz stress_yz'
strain = FINITE
[../]
[]
[Kernels]
[./heat]
type = Diffusion
variable = temp
[../]
[]
[BCs]
[./node1_x]
type = DirichletBC
variable = disp_x
boundary = 1
value = 0.0
[../]
[./node1_y]
type = FunctionDirichletBC
variable = disp_y
boundary = 1
function = ramp2
[../]
[./node1_z]
type = FunctionDirichletBC
variable = disp_z
boundary = 1
function = ramp3
[../]
[./node2_x]
type = FunctionDirichletBC
variable = disp_x
boundary = 2
function = ramp1
[../]
[./node2_y]
type = FunctionDirichletBC
variable = disp_y
boundary = 2
function = ramp2
[../]
[./node2_z]
type = FunctionDirichletBC
variable = disp_z
boundary = 2
function = ramp6
[../]
[./node3_x]
type = FunctionDirichletBC
variable = disp_x
boundary = 3
function = ramp1
[../]
[./node3_y]
type = DirichletBC
variable = disp_y
boundary = 3
value = 0.0
[../]
[./node3_z]
type = FunctionDirichletBC
variable = disp_z
boundary = 3
function = ramp3
[../]
[./node4_x]
type = DirichletBC
variable = disp_x
boundary = 4
value = 0.0
[../]
[./node4_y]
type = DirichletBC
variable = disp_y
boundary = 4
value = 0.0
[../]
[./node4_z]
type = DirichletBC
variable = disp_z
boundary = 4
value = 0.0
[../]
[./node5_x]
type = FunctionDirichletBC
variable = disp_x
boundary = 5
function = ramp1
[../]
[./node5_y]
type = FunctionDirichletBC
variable = disp_y
boundary = 5
function = ramp4
[../]
[./node5_z]
type = FunctionDirichletBC
variable = disp_z
boundary = 5
function = ramp3
[../]
[./node6_x]
type = FunctionDirichletBC
variable = disp_x
boundary = 6
function = ramp2
[../]
[./node6_y]
type = FunctionDirichletBC
variable = disp_y
boundary = 6
function = ramp4
[../]
[./node6_z]
type = FunctionDirichletBC
variable = disp_z
boundary = 6
function = ramp6
[../]
[./node7_x]
type = FunctionDirichletBC
variable = disp_x
boundary = 7
function = ramp2
[../]
[./node7_y]
type = FunctionDirichletBC
variable = disp_y
boundary = 7
function = ramp2
[../]
[./node7_z]
type = FunctionDirichletBC
variable = disp_z
boundary = 7
function = ramp3
[../]
[./node8_x]
type = FunctionDirichletBC
variable = disp_x
boundary = 8
function = ramp1
[../]
[./node8_y]
type = FunctionDirichletBC
variable = disp_y
boundary = 8
function = ramp2
[../]
[./node8_z]
type = DirichletBC
variable = disp_z
boundary = 8
value = 0.0
[../]
[./temp]
type = FunctionDirichletBC
variable = temp
boundary = '10 12'
function = tempFunc
[../]
[]
[Materials]
[./youngs_modulus]
type = PiecewiseLinearInterpolationMaterial
x = '100 500'
y = '1e6 6e5'
property = youngs_modulus
variable = temp
[../]
[./poissons_ratio]
type = PiecewiseLinearInterpolationMaterial
x = '100 500'
y = '0 0.25'
property = poissons_ratio
variable = temp
[../]
[./elasticity_tensor]
type = ComputeVariableIsotropicElasticityTensor
args = temp
youngs_modulus = youngs_modulus
poissons_ratio = poissons_ratio
[../]
[./stress]
type = ComputeFiniteStrainElasticStress
[../]
[]
[Executioner]
type = Transient
solve_type = 'PJFNK'
nl_rel_tol = 1e-9
nl_abs_tol = 1e-9
l_max_its = 20
start_time = 0.0
dt = 1.0
end_time = 2.0
[]
[Outputs]
exodus = true
[]
(modules/combined/examples/optimization/thermomechanical/thermal_sub.i)
vol_frac = 0.4
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 = 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 = '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
boundaries_old = 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
[]
[]
[Modules/TensorMechanics/Master]
[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
[]
# 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 = 500
[]
[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
function = '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
function = '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/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'
[]
[]
[]
[Modules/TensorMechanics/Master]
[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/solid_mechanics/test/tests/j_integral_vtest/fgm_5.i)
[GlobalParams]
order = FIRST
family = LAGRANGE
displacements = 'disp_x disp_y'
volumetric_locking_correction = true
[]
[Mesh]
[file]
type = FileMeshGenerator
file = 2d_paulino.e
[]
# uniform_refine = 3
[]
[AuxVariables]
[react_z]
[]
[]
[DomainIntegral]
integrals = 'JIntegral InteractionIntegralKI'
boundary = 1001
radius_inner = '0.01 0.04 0.1 0.2'
radius_outer = '0.01 0.04 0.1 0.2'
crack_direction_method = CrackDirectionVector
crack_direction_vector = '1 0 0' # is it +?
2d = true
axis_2d = 2
incremental = true
symmetry_plane = 1
functionally_graded_youngs_modulus = elastic_mod_material_mat
functionally_graded_youngs_modulus_crack_dir_gradient = elastic_mod_material_der_mat
youngs_modulus = 2e6
poissons_ratio = 0.3
block = '1'
[]
[Physics/SolidMechanics/QuasiStatic]
[master]
strain = FINITE
add_variables = true
incremental = true
generate_output = 'stress_xx stress_yy stress_zz vonmises_stress strain_xx strain_yy'
decomposition_method = EigenSolution
planar_formulation = PLANE_STRAIN
[]
[]
[Functions]
[parsed_load]
type = ParsedFunction
symbol_names = 'E1 E2 beta'
symbol_values = '1e3 3e3 5'
expression = '-1.0*((E1 + E2) / 2 + (E1 - E2)/2 * tanh(beta*(x+0.1)))'
[]
[elastic_mod_material_der]
type = ParsedFunction
symbol_names = 'E1 E2 beta'
symbol_values = '1e6 3e6 5'
expression = '(E1 - E2) / 2 * beta * (1.0 - tanh(beta*(x+0.1)) * tanh(beta*(x+0.1)))'
[]
[elastic_mod_material]
type = ParsedFunction
symbol_names = 'E1 E2 beta'
symbol_values = '1e6 3e6 5'
expression = '(E1 + E2) / 2 + (E1 - E2)/2 * tanh(beta*(x+0.1))'
[]
[]
[BCs]
[plane_1_x]
type = DirichletBC
variable = disp_x
boundary = 10001
value = 0.0
[]
[plane_y]
type = DirichletBC
variable = disp_y
boundary = '10005 6 1' #10001
value = 0.0
[]
[Pressure]
[Side1]
boundary = 4
function = parsed_load # BCs
[]
[]
[]
[Materials]
[generic_materials]
type = GenericFunctionMaterial
prop_names = 'elastic_mod_material_mat elastic_mod_material_der_mat'
prop_values = 'elastic_mod_material elastic_mod_material_der'
[]
[elasticity_tensor]
type = ComputeVariableIsotropicElasticityTensor
youngs_modulus = elastic_mod_material_mat
poissons_ratio = 0.3
args = ''
[]
[elastic_stress]
type = ComputeFiniteStrainElasticStress
[]
[]
[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'
l_max_its = 50
nl_max_its = 20
nl_abs_tol = 1e-5
nl_rel_tol = 1e-8
l_tol = 1e-6
start_time = 0.0
dt = 1.0
end_time = 1
num_steps = 1
[]
[Postprocessors]
[_dt]
type = TimestepSize
[]
[nl_its]
type = NumNonlinearIterations
[]
[lin_its]
type = NumLinearIterations
[]
[react_z]
type = NodalSum
variable = react_z
boundary = '10005 6 1'
[]
[]
[Outputs]
execute_on = 'timestep_end'
csv = true
exodus = true
[]
(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
[]
[]
[Modules/TensorMechanics/Master]
[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/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'
[]
[]
[]
[Modules/TensorMechanics/Master]
[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/examples/optimization/multi-load/square_subapp_two.i)
power = 1.0
E0 = 1.0
Emin = 1.0e-6
[GlobalParams]
displacements = 'disp_x disp_y'
[]
[Mesh]
[Bottom]
type = GeneratedMeshGenerator
dim = 2
nx = 100
ny = 100
xmin = 0
xmax = 150
ymin = 0
ymax = 150
[]
[left_load]
type = ExtraNodesetGenerator
input = Bottom
new_boundary = left_load
coord = '0 150 0'
[]
[right_load]
type = ExtraNodesetGenerator
input = left_load
new_boundary = right_load
coord = '150 150 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]
[]
[]
[AuxVariables]
[Dc]
family = MONOMIAL
order = SECOND
initial_condition = -1.0
[]
[Cc]
family = MONOMIAL
order = CONSTANT
initial_condition = -1.0
[]
[mat_den]
family = MONOMIAL
order = CONSTANT
initial_condition = 0.25
[]
[sensitivity_var]
family = MONOMIAL
order = SECOND
initial_condition = -1.0
[]
[]
[AuxKernels]
[sensitivity_kernel]
type = MaterialRealAux
check_boundary_restricted = false
property = sensitivity
variable = sensitivity_var
execute_on = 'TIMESTEP_END'
[]
[]
[Modules/TensorMechanics/Master]
[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
[]
[no_x_right]
type = DirichletBC
variable = disp_x
boundary = right_support
value = 0.0
[]
[]
[NodalKernels]
[push_right]
type = NodalGravity
variable = disp_y
boundary = right_load
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
# 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
[]
[]
[Preconditioning]
[smp]
type = SMP
full = true
[]
[]
[UserObjects]
# We do averaging in subapps
[rad_avg]
type = RadialAverage
radius = 8
weights = linear
prop_name = sensitivity
execute_on = TIMESTEP_END
force_preaux = true
[]
[calc_sense]
type = SensitivityFilter
density_sensitivity = Dc
design_density = mat_den
filter_UO = rad_avg
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 = 10
[]
[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
function = 'total_vol / mesh_volume'
pp_names = 'total_vol mesh_volume'
[]
[sensitivity]
type = ElementIntegralMaterialProperty
mat_prop = sensitivity
execute_on = 'TIMESTEP_BEGIN TIMESTEP_END NONLINEAR'
[]
[objective]
type = ElementIntegralMaterialProperty
mat_prop = strain_energy_density
execute_on = 'INITIAL TIMESTEP_END'
[]
[]
(modules/combined/examples/optimization/multi-load/single_main.i)
vol_frac = 0.3
power = 1.1
E0 = 1.0
Emin = 1.0e-6
[GlobalParams]
displacements = 'disp_x disp_y'
[]
[Mesh]
# final_generator = 'MoveRight'
[Bottom]
type = GeneratedMeshGenerator
dim = 2
nx = 80
ny = 40
xmin = 0
xmax = 150
ymin = 0
ymax = 75
[]
[left_load]
type = ExtraNodesetGenerator
input = Bottom
new_boundary = left_load
coord = '37.5 75 0'
[]
[right_load]
type = ExtraNodesetGenerator
input = left_load
new_boundary = right_load
coord = '112.5 75 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]
[]
[]
[AuxVariables]
[mat_den]
family = MONOMIAL
order = CONSTANT
initial_condition = 0.02
[]
[sensitivity_one]
family = MONOMIAL
order = SECOND
initial_condition = -1.0
[]
[sensitivity_two]
family = MONOMIAL
order = SECOND
initial_condition = -1.0
[]
[total_sensitivity]
family = MONOMIAL
order = SECOND
initial_condition = -1.0
[]
[]
[AuxKernels]
[total_sensitivity]
type = ParsedAux
variable = total_sensitivity
expression = '0.5*sensitivity_one + 0.5*sensitivity_two'
coupled_variables = 'sensitivity_one sensitivity_two'
execute_on = 'LINEAR TIMESTEP_END'
[]
[]
[Modules/TensorMechanics/Master]
[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
[]
[]
[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.0
[]
[stress]
type = ComputeLinearElasticStress
[]
[]
[Preconditioning]
[smp]
type = SMP
full = true
[]
[]
[UserObjects]
# We do filtering in the subapps
[update]
type = DensityUpdate
density_sensitivity = total_sensitivity
design_density = mat_den
volume_fraction = ${vol_frac}
execute_on = MULTIAPP_FIXED_POINT_BEGIN
[]
[]
[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
function = 'total_vol / mesh_volume'
pp_names = 'total_vol mesh_volume'
[]
[sensitivity]
type = ElementIntegralVariablePostprocessor
variable = total_sensitivity
[]
[]
[MultiApps]
[sub_app_one]
type = TransientMultiApp
input_files = single_subapp_one.i
[]
[sub_app_two]
type = TransientMultiApp
input_files = single_subapp_two.i
[]
[]
[Transfers]
# First SUB-APP
# To subapp densities
[subapp_one_density]
type = MultiAppCopyTransfer
to_multi_app = sub_app_one
source_variable = mat_den # Here
variable = mat_den
[]
# From subapp sensitivity
[subapp_one_sensitivity]
type = MultiAppCopyTransfer
from_multi_app = sub_app_one
source_variable = Dc # sensitivity_var
variable = sensitivity_one # Here
[]
# Second SUB-APP
# To subapp densities
[subapp_two_density]
type = MultiAppCopyTransfer
to_multi_app = sub_app_two
source_variable = mat_den # Here
variable = mat_den
[]
# From subapp sensitivity
[subapp_two_sensitivity]
type = MultiAppCopyTransfer
from_multi_app = sub_app_two
source_variable = Dc # sensitivity_var
variable = sensitivity_two # Here
[]
[]
(modules/solid_mechanics/test/tests/j_integral_vtest/j_int_fgm_sif.i)
[GlobalParams]
order = FIRST
family = LAGRANGE
displacements = 'disp_x disp_y disp_z'
volumetric_locking_correction = true
[]
[Mesh]
[file]
type = FileMeshGenerator
file = multiple_blocks_bimaterial.e
[]
partitioner = centroid
centroid_partitioner_direction = z
[]
[AuxVariables]
[SED]
order = CONSTANT
family = MONOMIAL
[]
[resid_z]
[]
[]
[Functions]
[rampConstantUp]
type = PiecewiseLinear
x = '0. 1.'
y = '0. 1'
scale_factor = -68.95 #MPa
[]
[elastic_mod_material_der]
type = ParsedFunction
expression = 'if(y < 229, 0.0, if(y>279, 0, 20680*0460517019*exp(0.0460517019*(y-229))))'
[]
[elastic_mod_material]
type = ParsedFunction
expression = 'if(y < 229, 20680, if(y>279, 206800, 20680*exp(0.0460517019*(y-229))))'
[]
[]
[DomainIntegral]
integrals = 'JIntegral InteractionIntegralKI'
boundary = 1001
crack_direction_method = CurvedCrackFront
crack_end_direction_method = CrackDirectionVector
crack_direction_vector_end_1 = '0.0 1.0 0.0'
crack_direction_vector_end_2 = '1.0 0.0 0.0'
radius_inner = '12.5 25.0 100'
radius_outer = '25.0 37.5 150.0'
intersecting_boundary = '1 2'
symmetry_plane = 2
incremental = true
functionally_graded_youngs_modulus = elastic_mod_material_mat
youngs_modulus = 20680
poissons_ratio = 0.3
block = '1 2'
[]
[Physics/SolidMechanics/QuasiStatic]
[master]
strain = FINITE
add_variables = true
incremental = true
generate_output = 'stress_xx stress_yy stress_zz vonmises_stress'
block = '1 2'
[]
[]
[AuxKernels]
[SED]
type = MaterialRealAux
variable = SED
property = strain_energy_density
execute_on = timestep_end
[]
[]
[BCs]
[crack_y]
type = DirichletBC
variable = disp_z
boundary = 6
value = 0.0
[]
[no_y]
type = DirichletBC
variable = disp_y
boundary = 12
value = 0.0
[]
[no_x]
type = DirichletBC
variable = disp_x
boundary = 1
value = 0.0
[]
[Pressure]
[Side1]
boundary = 5
function = rampConstantUp # BCs
[]
[]
[]
[Materials]
[generic_materials]
type = GenericFunctionMaterial
prop_names = 'elastic_mod_material_mat elastic_mod_material_der_mat'
prop_values = 'elastic_mod_material elastic_mod_material_der'
[]
[elasticity_tensor]
type = ComputeVariableIsotropicElasticityTensor
youngs_modulus = elastic_mod_material_mat
poissons_ratio = 0.3
args = ''
[]
[elastic_stress]
type = ComputeFiniteStrainElasticStress
block = '1 2'
[]
[]
[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'
l_max_its = 50
nl_max_its = 20
nl_abs_tol = 1e-5
nl_rel_tol = 1e-8
l_tol = 1e-6
start_time = 0.0
dt = 1
end_time = 1
num_steps = 1
[]
[Postprocessors]
[_dt]
type = TimestepSize
[]
[nl_its]
type = NumNonlinearIterations
[]
[lin_its]
type = NumLinearIterations
[]
[react_z]
type = NodalSum
variable = resid_z
boundary = 5
[]
[]
[Outputs]
execute_on = 'timestep_end'
csv = true
exodus = true
[]
(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'
[]
[]
[]
[Modules/TensorMechanics/Master]
[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'
[]
[]
(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
[]
[]
[]
[Modules/TensorMechanics/Master]
[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/combined/examples/optimization/multi-load/square_main.i)
# This example is intended to reproduce a 2D example with opposing horizontal
# loads (see [1]). This test has an undefined solution if reguar SIMP is applied.
# Using multi-loads SIMP, on the other hand, generates a structure that optimizes
# the response to both loads individually,
# [1]. Lat. Am. j. solids struct. 12 (5), May 2015
# Topological derivative-based topology optimization of structures subject to multiple load-cases
vol_frac = 0.5
power = 1.0
E0 = 1.0
Emin = 1.0e-6
[GlobalParams]
displacements = 'disp_x disp_y'
[]
[Mesh]
[Bottom]
type = GeneratedMeshGenerator
dim = 2
nx = 100
ny = 100
xmin = 0
xmax = 150
ymin = 0
ymax = 150
[]
[left_load]
type = ExtraNodesetGenerator
input = Bottom
new_boundary = left_load
coord = '0 150 0'
[]
[right_load]
type = ExtraNodesetGenerator
input = left_load
new_boundary = right_load
coord = '150 150 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]
[]
[]
[AuxVariables]
[mat_den]
family = MONOMIAL
order = CONSTANT
[]
[sensitivity_one]
family = MONOMIAL
order = SECOND
initial_condition = -1.0
[]
[sensitivity_two]
family = MONOMIAL
order = SECOND
initial_condition = -1.0
[]
[total_sensitivity]
family = MONOMIAL
order = SECOND
initial_condition = -1.0
[]
[]
[ICs]
[mat_den]
type = RandomIC
seed = 7
variable = mat_den
max = '${fparse vol_frac+0.35}'
min = '${fparse vol_frac-0.35}'
[]
[]
[AuxKernels]
[total_sensitivity]
type = ParsedAux
variable = total_sensitivity
expression = '0.5*sensitivity_one + 0.5*sensitivity_two'
coupled_variables = 'sensitivity_one sensitivity_two'
execute_on = 'LINEAR TIMESTEP_END'
[]
[]
[Modules/TensorMechanics/Master]
[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
[]
[no_x_right]
type = DirichletBC
variable = disp_x
boundary = right_support
value = 0.0
[]
[]
[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
[]
[]
[Preconditioning]
[smp]
type = SMP
full = true
[]
[]
[UserObjects]
[update]
type = DensityUpdate
density_sensitivity = total_sensitivity
design_density = mat_den
volume_fraction = ${vol_frac}
execute_on = MULTIAPP_FIXED_POINT_BEGIN
[]
[]
[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 = 10
[]
[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
function = 'total_vol / mesh_volume'
pp_names = 'total_vol mesh_volume'
[]
[sensitivity]
type = ElementIntegralVariablePostprocessor
variable = total_sensitivity
[]
[]
[MultiApps]
[sub_app_one]
type = TransientMultiApp
input_files = square_subapp_one.i
[]
[sub_app_two]
type = TransientMultiApp
input_files = square_subapp_two.i
[]
[]
[Transfers]
# First SUB-APP
# To subapp densities
[subapp_one_density]
type = MultiAppCopyTransfer
to_multi_app = sub_app_one
source_variable = mat_den # Here
variable = mat_den
[]
# From subapp sensitivity
[subapp_one_sensitivity]
type = MultiAppCopyTransfer
from_multi_app = sub_app_one
source_variable = Dc # sensitivity_var
variable = sensitivity_one # Here
[]
# Second SUB-APP
# To subapp densities
[subapp_two_density]
type = MultiAppCopyTransfer
to_multi_app = sub_app_two
source_variable = mat_den # Here
variable = mat_den
[]
# From subapp sensitivity
[subapp_two_sensitivity]
type = MultiAppCopyTransfer
from_multi_app = sub_app_two
source_variable = Dc # sensitivity_var
variable = sensitivity_two # Here
[]
[]
(modules/combined/test/tests/optimization/optimization_density_update/top_opt_3d.i)
vol_frac = 0.5
E0 = 1e5
Emin = 1e-2
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'
[]
[]
[AuxVariables]
[compliance]
family = MONOMIAL
order = CONSTANT
[]
[Dc]
family = MONOMIAL
order = CONSTANT
initial_condition = -1.0
[]
[mat_den]
family = MONOMIAL
order = CONSTANT
initial_condition = ${vol_frac}
[]
[]
[Modules/TensorMechanics/Master]
[all]
strain = SMALL
add_variables = true
incremental = false
[]
[]
[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
[]
[]
[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]
[rad_avg]
type = RadialAverage
radius = 0.5
weights = constant
prop_name = sensitivity
execute_on = TIMESTEP_END
execution_order_group = -1
[]
[update]
type = DensityUpdate
density_sensitivity = Dc
design_density = mat_den
volume_fraction = ${vol_frac}
execute_on = TIMESTEP_BEGIN
[]
[calc_sense]
type = SensitivityFilter
density_sensitivity = Dc
design_density = mat_den
filter_UO = rad_avg
execute_on = TIMESTEP_END
[]
[]
[Executioner]
type = Transient
solve_type = NEWTON
petsc_options_iname = '-pc_type'
petsc_options_value = 'lu '
nl_abs_tol = 1e-10
l_max_its = 200
start_time = 0.0
dt = 1.0
num_steps = 10
[]
[Outputs]
[out]
type = Exodus
time_step_interval = 10
[]
[]
(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
[]
[]
[]
[Modules/TensorMechanics/Master]
[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/test/tests/optimization/optimization_density_update/top_opt_2d.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
[]
[]
[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
[]
[Dc]
family = MONOMIAL
order = CONSTANT
initial_condition = -1.0
[]
[mat_den]
family = MONOMIAL
order = CONSTANT
initial_condition = ${vol_frac}
[]
[]
[Modules/TensorMechanics/Master]
[all]
strain = SMALL
add_variables = true
incremental = false
[]
[]
[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
[]
[]
[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]
[rad_avg]
type = RadialAverage
radius = 0.5
weights = constant
prop_name = sensitivity
execute_on = TIMESTEP_END
force_preaux = true
execution_order_group = -1
[]
[update]
type = DensityUpdate
density_sensitivity = Dc
design_density = mat_den
volume_fraction = ${vol_frac}
execute_on = TIMESTEP_BEGIN
[]
[calc_sense]
type = SensitivityFilter
density_sensitivity = Dc
design_density = mat_den
filter_UO = rad_avg
execute_on = TIMESTEP_END
[]
[]
[Executioner]
type = Transient
solve_type = NEWTON
petsc_options_iname = '-pc_type '
petsc_options_value = 'lu'
nl_abs_tol = 1e-10
start_time = 0.0
dt = 1.0
num_steps = 50
[]
[Outputs]
[out]
type = Exodus
time_step_interval = 10
[]
[]
(modules/combined/examples/optimization/2d_mbb.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 = pull
nodes = 0
[]
[push]
type = ExtraNodesetGenerator
input = node
new_boundary = push
coord = '30 10 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}
[]
[Dc]
family = MONOMIAL
order = CONSTANT
initial_condition = -1.0
[]
[mat_den]
family = MONOMIAL
order = CONSTANT
initial_condition = ${vol_frac}
[]
[]
[Modules/TensorMechanics/Master]
[all]
strain = SMALL
add_variables = true
incremental = false
[]
[]
[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
[]
[]
[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]
[rad_avg]
type = RadialAverage
radius = 1.2
weights = linear
prop_name = sensitivity
execute_on = TIMESTEP_END
force_preaux = true
[]
[update]
type = DensityUpdate
density_sensitivity = Dc
design_density = mat_den
volume_fraction = ${vol_frac}
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
[]
[]
[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 = 70
[]
[Outputs]
[out]
type = Exodus
execute_on = 'TIMESTEP_END'
[]
print_linear_residuals = false
[]
[Postprocessors]
[total_vol]
type = ElementIntegralVariablePostprocessor
variable = mat_den
execute_on = 'INITIAL TIMESTEP_END'
[]
[]
(modules/combined/examples/optimization/multi-load/square_subapp_one.i)
power = 1.0
E0 = 1.0
Emin = 1.0e-6
[GlobalParams]
displacements = 'disp_x disp_y'
[]
[Mesh]
[Bottom]
type = GeneratedMeshGenerator
dim = 2
nx = 100
ny = 100
xmin = 0
xmax = 150
ymin = 0
ymax = 150
[]
[left_load]
type = ExtraNodesetGenerator
input = Bottom
new_boundary = left_load
coord = '0 150 0'
[]
[right_load]
type = ExtraNodesetGenerator
input = left_load
new_boundary = right_load
coord = '150 150 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]
[]
[]
[AuxVariables]
[Dc]
family = MONOMIAL
order = SECOND
initial_condition = -1.0
[]
[Cc]
family = MONOMIAL
order = CONSTANT
initial_condition = -1.0
[]
[mat_den]
family = MONOMIAL
order = CONSTANT
initial_condition = 0.25
[]
[sensitivity_var]
family = MONOMIAL
order = SECOND
initial_condition = -1.0
[]
[]
[AuxKernels]
[sensitivity_kernel]
type = MaterialRealAux
property = sensitivity
variable = sensitivity_var
check_boundary_restricted = false
execute_on = 'TIMESTEP_END'
[]
[]
[Modules/TensorMechanics/Master]
[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
[]
[no_x_right]
type = DirichletBC
variable = disp_x
boundary = right_support
value = 0.0
[]
[]
[NodalKernels]
[push_left]
type = NodalGravity
variable = disp_y
boundary = left_load
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
# 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
[]
[]
[Preconditioning]
[smp]
type = SMP
full = true
[]
[]
[UserObjects]
# We do averaging in subapps
[rad_avg]
type = RadialAverage
radius = 8
weights = linear
prop_name = sensitivity
force_preaux = true
execute_on = 'TIMESTEP_END'
[]
# Provides Dc
[calc_sense]
type = SensitivityFilter
density_sensitivity = Dc
design_density = mat_den
filter_UO = rad_avg
force_postaux = true
execute_on = 'TIMESTEP_END'
[]
[]
[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 = 10
[]
[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
function = 'total_vol / mesh_volume'
pp_names = 'total_vol mesh_volume'
[]
[sensitivity]
type = ElementIntegralMaterialProperty
mat_prop = sensitivity
execute_on = 'TIMESTEP_BEGIN TIMESTEP_END NONLINEAR'
[]
[objective]
type = ElementIntegralMaterialProperty
mat_prop = strain_energy_density
execute_on = 'INITIAL TIMESTEP_END'
[]
[]
(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
[]
[]
[]
[Modules/TensorMechanics/Master]
[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'
[]
[]
(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'
[]
[]
[]
[Modules/TensorMechanics/Master]
[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/combined/test/tests/thermal_elastic/derivatives.i)
[Mesh]
type = GeneratedMesh
dim = 1
nx = 2
[]
[Variables]
[./a]
[./InitialCondition]
type = RandomIC
min = -1
max = 1
[../]
[../]
[./b]
[./InitialCondition]
type = RandomIC
min = -1
max = 1
[../]
[../]
[]
[Debug]
[./MaterialDerivativeTest]
[./elastic]
prop_name = elasticity_tensor
prop_type = RankFourTensor
derivative_order = 1
args = 'a b'
[../]
[../]
[]
[Problem]
kernel_coverage_check = false
[]
[Materials]
[./youngs_modulus]
type = DerivativeParsedMaterial
property_name = youngs_modulus
expression = '23.1 * a^4 + 10.7 * b^2'
coupled_variables = 'a b'
[../]
[./poissons_ratio]
type = DerivativeParsedMaterial
property_name = poissons_ratio
expression = '0.2 * a^2 + 0.29 * b^3'
coupled_variables = 'a b'
[../]
[./elasticity_tensor]
type = ComputeVariableIsotropicElasticityTensor
args = 'a b'
youngs_modulus = youngs_modulus
poissons_ratio = poissons_ratio
[../]
[]
[Executioner]
type = Steady
[]
(modules/combined/examples/optimization/thermomechanical/thermomechanical_main.i)
vol_frac = 0.4
power = 2.0
E0 = 1.0e-6
E1 = 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 = '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
boundaries_old = left
[]
[dirichlet_bc]
type = SideSetsFromNodeSetsGenerator
input = extra
[]
[]
[Variables]
[disp_x]
[]
[disp_y]
[]
[]
[AuxVariables]
[mat_den]
family = MONOMIAL
order = FIRST
initial_condition = 0.02
[]
[sensitivity_one]
family = MONOMIAL
order = FIRST
initial_condition = -1.0
[]
[sensitivity_two]
family = MONOMIAL
order = FIRST
initial_condition = -1.0
[]
[total_sensitivity]
family = MONOMIAL
order = FIRST
initial_condition = -1.0
[]
[]
[AuxKernels]
[total_sensitivity]
type = ParsedAux
variable = total_sensitivity
expression = '(1-1.0e-7)*sensitivity_one + 1.0e-7*sensitivity_two'
coupled_variables = 'sensitivity_one sensitivity_two'
execute_on = 'LINEAR TIMESTEP_END'
[]
[]
[Modules/TensorMechanics/Master]
[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
[]
[]
[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 = '${E1} + (mat_den ^ ${power}) * (${E1}-${E0})'
coupled_variables = 'mat_den'
property_name = E_phys
[]
[poissons_ratio]
type = GenericConstantMaterial
prop_names = poissons_ratio
prop_values = 0.3
[]
[stress]
type = ComputeLinearElasticStress
[]
[]
[Preconditioning]
[smp]
type = SMP
full = true
[]
[]
[UserObjects]
# We do filtering in the subapps
[update]
type = DensityUpdate
density_sensitivity = total_sensitivity
design_density = mat_den
volume_fraction = ${vol_frac}
execute_on = MULTIAPP_FIXED_POINT_BEGIN
[]
[]
[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 = 2
[]
[Outputs]
[out]
type = CSV
execute_on = 'TIMESTEP_END'
[]
print_linear_residuals = false
exodus = true
[]
[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
function = 'total_vol / mesh_volume'
pp_names = 'total_vol mesh_volume'
[]
[sensitivity]
type = ElementIntegralVariablePostprocessor
variable = total_sensitivity
[]
[]
[MultiApps]
[sub_app_one]
type = TransientMultiApp
input_files = structural_sub.i
[]
[sub_app_two]
type = TransientMultiApp
input_files = thermal_sub.i
[]
[]
[Transfers]
# First SUB-APP: STRUCTURAL
# To subapp densities
[subapp_one_density]
type = MultiAppCopyTransfer
to_multi_app = sub_app_one
source_variable = mat_den # Here
variable = mat_den
[]
# From subapp sensitivity
[subapp_one_sensitivity]
type = MultiAppCopyTransfer
from_multi_app = sub_app_one
source_variable = Dc # sensitivity_var
variable = sensitivity_one # Here
[]
# Second SUB-APP: HEAT CONDUCTIVITY
# To subapp densities
[subapp_two_density]
type = MultiAppCopyTransfer
to_multi_app = sub_app_two
source_variable = mat_den # Here
variable = mat_den
[]
# From subapp sensitivity
[subapp_two_sensitivity]
type = MultiAppCopyTransfer
from_multi_app = sub_app_two
source_variable = Tc # sensitivity_var
variable = sensitivity_two # Here
[]
[]
(modules/combined/examples/optimization/thermomechanical/structural_sub.i)
vol_frac = 0.4
power = 2.0
E0 = 1.0e-6
E1 = 1.0
rho0 = 0.0
rho1 = 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 = '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
boundaries_old = left
[]
[dirichlet_bc]
type = SideSetsFromNodeSetsGenerator
input = extra
[]
[]
[Variables]
[disp_x]
[]
[disp_y]
[]
[]
[AuxVariables]
[Dc]
family = MONOMIAL
order = FIRST
initial_condition = -1.0
[]
[mat_den]
family = MONOMIAL
order = FIRST
initial_condition = ${vol_frac}
[]
[]
[Modules/TensorMechanics/Master]
[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 = -1.0e-3
mass = 1
[]
[push_center]
type = NodalGravity
variable = disp_y
boundary = push_center
gravity_value = -1.0e-3
mass = 1
[]
[]
[Materials]
[elasticity_tensor]
type = ComputeVariableIsotropicElasticityTensor
youngs_modulus = E_phys
poissons_ratio = poissons_ratio
args = 'mat_den'
[]
[E_phys]
type = DerivativeParsedMaterial
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
[]
[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
[]
[]
[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
[]
[calc_sense]
type = SensitivityFilter
density_sensitivity = Dc
design_density = mat_den
filter_UO = rad_avg
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 = 500
[]
[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
function = 'total_vol / mesh_volume'
pp_names = 'total_vol mesh_volume'
[]
[sensitivity]
type = ElementIntegralMaterialProperty
mat_prop = sensitivity
[]
[objective]
type = ElementIntegralMaterialProperty
mat_prop = strain_energy_density
execute_on = 'INITIAL TIMESTEP_END'
[]
[]
(modules/combined/test/tests/optimization/compliance_sensitivity/2d_mmb_2material_cost.i)
vol_frac = 0.5
power = 3
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
[]
[mat_den]
family = MONOMIAL
order = CONSTANT
initial_condition = ${vol_frac}
[]
[]
[Modules/TensorMechanics/Master]
[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]
type = DerivativeParsedMaterial
# ordered multimaterial simp
expression = "A1:=(${C0}-${C1})/(${rho0}^${power}-${rho1}^${power}); "
"B1:=${C0}-A1*${rho0}^${power}; C1:=A1*mat_den^${power}+B1; "
"A2:=(${C1}-${C2})/(${rho1}^${power}-${rho2}^${power}); "
"B2:=${C1}-A2*${rho1}^${power}; C2:=A2*mat_den^${power}+B2; "
"if(mat_den<${rho1},C1,C2)"
coupled_variables = 'mat_den'
property_name = 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
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
[]
[update]
type = DensityUpdate
density_sensitivity = Dc
design_density = mat_den
volume_fraction = ${vol_frac}
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
[]
[]
[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 = 70
[]
[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
[]
[]
(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
[]
[]
[Modules/TensorMechanics/Master]
[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
function = '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/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
[]
[]
[Modules/TensorMechanics/Master]
[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
function = '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
function = '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/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
[]
[]
[Modules/TensorMechanics/Master]
[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
function = '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/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
[]
[]
[Modules/TensorMechanics/Master]
[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
function = '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
function = '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/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
[]
[]
[]
[Modules/TensorMechanics/Master]
[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'
[]
[]
(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
[]
[]
[Modules/TensorMechanics/Master]
[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
function = '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
function = 'cost / mesh_volume'
pp_names = 'cost mesh_volume'
[]
[objective]
type = ElementIntegralMaterialProperty
mat_prop = strain_energy_density
execute_on = 'INITIAL TIMESTEP_END'
[]
[]
(modules/combined/examples/optimization/multi-load/single_subapp_two.i)
power = 2
E0 = 1.0
Emin = 1.0e-6
[GlobalParams]
displacements = 'disp_x disp_y'
[]
[Mesh]
[Bottom]
type = GeneratedMeshGenerator
dim = 2
nx = 80
ny = 40
xmin = 0
xmax = 150
ymin = 0
ymax = 75
[]
[left_load]
type = ExtraNodesetGenerator
input = Bottom
new_boundary = left_load
coord = '37.5 75 0'
[]
[right_load]
type = ExtraNodesetGenerator
input = left_load
new_boundary = right_load
coord = '112.5 75 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]
[]
[]
[AuxVariables]
[Dc]
family = MONOMIAL
order = SECOND
initial_condition = -1.0
[]
[Cc]
family = MONOMIAL
order = CONSTANT
initial_condition = -1.0
[]
[mat_den]
family = MONOMIAL
order = CONSTANT
initial_condition = 0.1
[]
[sensitivity_var]
family = MONOMIAL
order = SECOND
initial_condition = -1.0
[]
[]
[AuxKernels]
[sensitivity_kernel]
type = MaterialRealAux
check_boundary_restricted = false
property = sensitivity
variable = sensitivity_var
execute_on = 'TIMESTEP_END'
[]
[]
[Modules/TensorMechanics/Master]
[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
[]
[]
[NodalKernels]
[push_right]
type = NodalGravity
variable = disp_y
boundary = right_load
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
# 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.0
[]
[stress]
type = ComputeLinearElasticStress
[]
[dc]
type = ComplianceSensitivity
design_density = mat_den
youngs_modulus = E_phys
[]
[]
[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
[]
# No SIMP optimization in subapp
[calc_sense]
type = SensitivityFilter
density_sensitivity = Dc
design_density = mat_den
filter_UO = rad_avg
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
function = 'total_vol / mesh_volume'
pp_names = 'total_vol mesh_volume'
[]
[sensitivity]
type = ElementIntegralMaterialProperty
mat_prop = sensitivity
execute_on = 'TIMESTEP_BEGIN TIMESTEP_END NONLINEAR'
[]
[objective]
type = ElementIntegralMaterialProperty
mat_prop = strain_energy_density
execute_on = 'INITIAL TIMESTEP_END'
[]
[]