Shifted Cohesive Zone Physics System

Action to create Shifted Cohesive Zone Method (SCZM) objects for non-interface-fitted meshes.

The Shifted Cohesive Zone Physics creates the interface kernels, displacement-jump and global-traction materials, and optional geometry objects needed by the shifted cohesive-zone method (SCZM). Only small-strain kinematics are currently supported.

Setting up an SCZM problem involves two independent choices:

Select how the Physics obtains the distance and true normal for each interface.
Define a constitutive traction material on the same surrogate boundaries.

Select the interface geometry

The boundary parameter always lists the surrogate boundaries on which SCZM is applied. Exactly one of the following three approaches supplies their corresponding true-interface geometry.

Use an existing distance user object

Set sbm_distance_uo to an existing BoundaryShortestDistanceToSurface for full control over the distance calculation. The following example defines the object from a level-set function:

[UserObjects<<<{"href": "../../../UserObjects/index.html"}>>>]
  [sbm_distance_uo]
    type = BoundaryShortestDistanceToSurface<<<{"description": "Provides, at each quadrature point on a (surrogate) boundary, the vector from that point to its projection on the true surface, together with the true surface normal there. The true surface may be described by level-set/distance functions or a mesh-based representation.", "href": "../../../../source/userobjects/BoundaryShortestDistanceToSurface.html"}>>>
    surfaces<<<{"description": "Level-set or mesh-based functions that define the boundary distances."}>>> = signed_dist_func
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = 'Block1_Block2'
    execution_order_group<<<{"description": "Execution order groups are executed in increasing order (e.g., the lowest number is executed first). Note that negative group numbers may be used to execute groups before the default (0) group. Please refer to the user object documentation for ordering of user object execution within a group."}>>> = 0
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = INITIAL
  []
[]

The Physics can receive the object name directly or inherit it from [GlobalParams], as this example does:

[GlobalParams<<<{"href": "../../../GlobalParams/index.html"}>>>]
  displacements = 'disp_x disp_y'
  use_displaced_mesh = false
  sbm_distance_uo = sbm_distance_uo
[]

This approach is the most flexible, but the user is responsible for constructing the distance user object and any functions or surface data it needs.

Use one surface mesh per interface

Set generate_sbm_distance = true when a separate saved surface mesh is available for every interface. The Physics creates the surface-mesh builders, distance functions, and aggregate boundary distance user object. If each saved mesh has the same name as its corresponding surrogate boundary, only the following Physics configuration is needed:

[Physics<<<{"href": "../../index.html"}>>>/SolidMechanics<<<{"href": "../index.html"}>>>/ShiftedCohesiveZone<<<{"href": "index.html"}>>>]
  [czm_ik]
    use_automatic_differentiation<<<{"description": "Whether to use automatic differentiation to compute the Jacobian"}>>> = false
    boundary<<<{"description": "The list of boundary IDs from the mesh where the cohesive zone will be applied"}>>> = ${interfaces}
    generate_sbm_distance<<<{"description": "Automatically create surface mesh builders, distance functions, and the boundary distance user object from the supplied boundaries."}>>> = true
    check_surface_watertightness<<<{"description": "Whether generated surface mesh builders check mesh watertightness."}>>> = true
  []
[]

Set surface_meshes only when the saved mesh names differ from the boundary names. Entries in surface_meshes and boundary correspond by position. For example:


boundary = 'grain1_grain2 grain1_grain4'
surface_meshes = 'surface_12 surface_14'

Here, surface_12 supplies the geometry for grain1_grain2, and surface_14 supplies the geometry for grain1_grain4.

Use one complete interface mesh

Set complete_interface_mesh to one saved mesh containing the complete outward-oriented boundary of every subdomain. The Physics creates one SBMInterfaceManager and one boundary distance user object:

[Physics<<<{"href": "../../index.html"}>>>/SolidMechanics<<<{"href": "../index.html"}>>>/ShiftedCohesiveZone<<<{"href": "index.html"}>>>]
  [czm_ik]
    use_automatic_differentiation<<<{"description": "Whether to use automatic differentiation to compute the Jacobian"}>>> = false
    boundary<<<{"description": "The list of boundary IDs from the mesh where the cohesive zone will be applied"}>>> = ${interfaces}
    complete_interface_mesh<<<{"description": "Saved surface mesh containing the complete boundary of each subdomain. Setting this parameter automatically creates one interface manager and distance user object."}>>> = boundary_mesh
  []
[]

The mesh generator that supplies complete_interface_mesh must retain the mesh with save_with_name. Boundary names of the form grainX_grainY or blockX_blockY identify the two interface subdomains automatically. For other boundary names, set interface_subdomain_pairs with one pair per boundary in the same order. The pair selects an interface; its order does not control the final true-normal orientation.

Shifted integration is enabled by default and is the intended mode for SBM simulations. The no_shifted = true option is retained primarily for verification and comparison: it disables the shifted correction terms while using the same interface kernel and cohesive-zone model, providing a non-shifted baseline that isolates their effect. This option does not require a separate geometry source.

Define the traction law

The Physics does not create the constitutive traction law. Define an interface material such as BiLinearMixedModeTraction or SalehaniIrani3DCTraction and restrict it to the same surrogate interfaces/ boundaries as the SCZM Physics. The Action also creates the displacement-jump and global-traction materials. The displacement-jump material computes the local displacement jump used by the constitutive material to determine the local traction. The global-traction material then uses the local traction to compute the traction required by the equilibrium equations.

For multiple interfaces, define the boundary list once at the top of the input file and use ${interfaces} wherever it is needed. The ${interfaces} expression is replaced with the boundary list assigned to interfaces. Reusing this list in both the Physics and the traction material keeps their boundary specifications consistent. The complete-interface example uses this pattern:

interfaces = 'grain1_grain2 grain1_grain4 grain1_grain5 grain2_grain4 grain3_grain4 grain3_grain5 grain4_grain5'

[Physics<<<{"href": "../../index.html"}>>>/SolidMechanics<<<{"href": "../index.html"}>>>/ShiftedCohesiveZone<<<{"href": "index.html"}>>>]
  [czm_ik]
    use_automatic_differentiation<<<{"description": "Whether to use automatic differentiation to compute the Jacobian"}>>> = false
    boundary<<<{"description": "The list of boundary IDs from the mesh where the cohesive zone will be applied"}>>> = ${interfaces}
    complete_interface_mesh<<<{"description": "Saved surface mesh containing the complete boundary of each subdomain. Setting this parameter automatically creates one interface manager and distance user object."}>>> = boundary_mesh
  []
[]

[Materials<<<{"href": "../../../Materials/index.html"}>>>]
  [elastic_stress]
    type = ComputeLinearElasticStress<<<{"description": "Compute stress using elasticity for small strains", "href": "../../../../source/materials/ComputeLinearElasticStress.html"}>>>
  []
  [elasticity_tensor]
    type = ComputeIsotropicElasticityTensor<<<{"description": "Compute a constant isotropic elasticity tensor.", "href": "../../../../source/materials/ComputeIsotropicElasticityTensor.html"}>>>
    poissons_ratio<<<{"description": "Poisson's ratio for the material."}>>> = ${poi}
    youngs_modulus<<<{"description": "Young's modulus of the material."}>>> = ${E}
  []
  [interface_traction]
    type = BiLinearMixedModeTraction<<<{"description": "Mixed mode bilinear traction separation law.", "href": "../../../../source/materials/cohesive_zone_model/BiLinearMixedModeTraction.html"}>>>
    boundary<<<{"description": "The list of boundaries (ids or names) from the mesh where this object applies"}>>> = ${interfaces}
    penalty_stiffness<<<{"description": "Penalty stiffness."}>>> = 2e3
    GI_c<<<{"description": "Critical energy release rate in normal direction."}>>> = 1e3
    GII_c<<<{"description": "Critical energy release rate in shear direction."}>>> = 1e2
    normal_strength<<<{"description": "Tensile strength in normal direction."}>>> = 500
    shear_strength<<<{"description": "Tensile strength in shear direction."}>>> = 300
    displacements<<<{"description": "The string of displacements suitable for the problem statement"}>>> = 'disp_x disp_y'
    eta<<<{"description": "The power law parameter."}>>> = 2.2
    viscosity<<<{"description": "Viscosity."}>>> = 1e-3
  []
[]

Keeping the traction material explicit allows different constitutive laws to be assigned to different subsets of the SCZM boundaries. The Physics sets boundary, displacements, and base_name only on the kernels and materials that it creates; it does not modify user-defined materials in [Materials].

Non-AD tangent definition

The non-AD interface kernel uses stress and tangent material properties for its directional correction and Jacobian. When tangent_definition = auto, the standard Jacobian_mult property is interpreted as $var element = document.getElementById("moose-equation-59dde0ea-c243-46c4-bf11-8ddc1b0e5a16");katex.render("\\partial\\boldsymbol{\\sigma}/\\partial\\boldsymbol{\\varepsilon}", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ and pk1_jacobian is interpreted as $var element = document.getElementById("moose-equation-e6b59bab-6bea-40a9-b76a-223c61b9aaa3");katex.render("\\partial\\mathbf{P}/\\partial\\mathbf{F}", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ . A custom tangent property name requires an explicit tangent_definition:

Use stress_wrt_strain only for a derivative of stress with respect to symmetric strain ( $var element = document.getElementById("moose-equation-003e52f5-2e0a-4228-a609-ed3c0f57498e");katex.render("\\partial\\boldsymbol{\\sigma}/\\partial\\boldsymbol{\\varepsilon}", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ ).
Use pk1_wrt_deformation_gradient only for a derivative of first Piola-Kirchhoff stress with respect to the full deformation gradient ( $var element = document.getElementById("moose-equation-b80b1f01-cf25-490c-8813-b79f08e4046c");katex.render("\\partial\\mathbf{P}/\\partial\\mathbf{F}", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ ).

The AD interface kernel differentiates its residual automatically and does not use tangent or tangent_definition; setting either parameter with AD enabled is an error.

Input Parameters

boundaryThe list of boundary IDs from the mesh where the cohesive zone will be applied
C++ Type:std::vector<BoundaryName>
Controllable:No
Description:The list of boundary IDs from the mesh where the cohesive zone will be applied
displacementsThe nonlinear displacement variables for the problem
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:The nonlinear displacement variables for the problem

Required Parameters

active__all__ If specified only the blocks named will be visited and made active
Default:__all__
C++ Type:std::vector<std::string>
Controllable:No
Description:If specified only the blocks named will be visited and made active
base_nameMaterial property base name
C++ Type:std::string
Controllable:No
Description:Material property base name
check_surface_watertightnessFalseWhether generated surface mesh builders check mesh watertightness.
Default:False
C++ Type:bool
Controllable:No
Description:Whether generated surface mesh builders check mesh watertightness.
complete_interface_meshSaved surface mesh containing the complete boundary of each subdomain. Setting this parameter automatically creates one interface manager and distance user object.
C++ Type:MeshGeneratorName
Controllable:No
Description:Saved surface mesh containing the complete boundary of each subdomain. Setting this parameter automatically creates one interface manager and distance user object.
directional_correctionTrueAdd the directional correction terms.
Default:True
C++ Type:bool
Controllable:No
Description:Add the directional correction terms.
generate_sbm_distanceFalseAutomatically create surface mesh builders, distance functions, and the boundary distance user object from the supplied boundaries.
Default:False
C++ Type:bool
Controllable:No
Description:Automatically create surface mesh builders, distance functions, and the boundary distance user object from the supplied boundaries.
inactiveIf specified blocks matching these identifiers will be skipped.
C++ Type:std::vector<std::string>
Controllable:No
Description:If specified blocks matching these identifiers will be skipped.
interface_matching_tolerance1e-06Relative face-matching tolerance for the single-mesh interface manager.
Default:1e-06
C++ Type:Real
Unit:(no unit assumed)
Controllable:No
Description:Relative face-matching tolerance for the single-mesh interface manager.
interface_subdomain_pairsOrdered subdomain pair for each boundary when it cannot be inferred from boundary names.
C++ Type:std::vector<std::vector<unsigned short>>
Controllable:No
Description:Ordered subdomain pair for each boundary when it cannot be inferred from boundary names.
no_shiftedFalseDisable shifted terms.
Default:False
C++ Type:bool
Controllable:No
Description:Disable shifted terms.
sbm_distance_uoUserObject that provides signed distance and normal vector calculations.
C++ Type:UserObjectName
Controllable:No
Description:UserObject that provides signed distance and normal vector calculations.
strainSMALLStrain formulation
Default:SMALL
C++ Type:MooseEnum
Options:SMALL, FINITE
Controllable:No
Description:Strain formulation
stressstressName of the stress tensor material property for small strain SCZM.
Default:stress
C++ Type:MaterialPropertyName
Unit:(no unit assumed)
Controllable:No
Description:Name of the stress tensor material property for small strain SCZM.
surface_meshesNames of saved surface meshes corresponding one-to-one with boundary. If omitted, the boundary names are used.
C++ Type:std::vector<MeshGeneratorName>
Controllable:No
Description:Names of saved surface meshes corresponding one-to-one with boundary. If omitted, the boundary names are used.
tangentJacobian_multName of the material Jacobian tensor property for small strain SCZM.
Default:Jacobian_mult
C++ Type:MaterialPropertyName
Unit:(no unit assumed)
Controllable:No
Description:Name of the material Jacobian tensor property for small strain SCZM.
tangent_definitionautoMathematical definition of the tangent material property. 'auto' recognizes the standard Jacobian_mult property as d(stress)/d(strain) and pk1_jacobian as d(PK1)/d(F). Select 'stress_wrt_strain' or 'pk1_wrt_deformation_gradient' explicitly for a custom tangent property; use the latter only when the property is known to be d(PK1)/d(F).
Default:auto
C++ Type:MooseEnum
Options:auto, stress_wrt_strain, pk1_wrt_deformation_gradient
Controllable:No
Description:Mathematical definition of the tangent material property. 'auto' recognizes the standard Jacobian_mult property as d(stress)/d(strain) and pk1_jacobian as d(PK1)/d(F). Select 'stress_wrt_strain' or 'pk1_wrt_deformation_gradient' explicitly for a custom tangent property; use the latter only when the property is known to be d(PK1)/d(F).
use_automatic_differentiationFalseWhether to use automatic differentiation to compute the Jacobian
Default:False
C++ Type:bool
Controllable:No
Description:Whether to use automatic differentiation to compute the Jacobian
verboseFalseDisplay extra information.
Default:False
C++ Type:bool
Controllable:No
Description:Display extra information.
volumetric_locking_correctionFalseWhether to apply volumetric locking correction to the directional correction Jacobian term in SCZM. If not specified, this value is automatically inherited from the solid mechanics physics block.
Default:False
C++ Type:bool
Controllable:No
Description:Whether to apply volumetric locking correction to the directional correction Jacobian term in SCZM. If not specified, this value is automatically inherited from the solid mechanics physics block.

Optional Parameters

additional_generate_outputAdd scalar quantity output for stress and/or strain (will be appended to the list in `generate_output`)
C++ Type:MultiMooseEnum
Options:jump_x, jump_y, jump_z, pk1_traction_x, pk1_traction_y, pk1_traction_z, traction_x, traction_y, traction_z, normal_traction, normal_jump, tangent_traction, tangent_jump
Controllable:No
Description:Add scalar quantity output for stress and/or strain (will be appended to the list in `generate_output`)
additional_material_output_familySpecifies the family of FE shape functions to use for this variable.
C++ Type:MultiMooseEnum
Options:MONOMIAL
Controllable:No
Description:Specifies the family of FE shape functions to use for this variable.
additional_material_output_orderSpecifies the order of the FE shape function to use for this variable.
C++ Type:MultiMooseEnum
Options:CONSTANT, FIRST, SECOND, THIRD, FOURTH, FIFTH, SIXTH, SEVENTH, EIGHTH, NINTH
Controllable:No
Description:Specifies the order of the FE shape function to use for this variable.
generate_outputAdd scalar quantity output for stress and/or strain
C++ Type:MultiMooseEnum
Options:jump_x, jump_y, jump_z, pk1_traction_x, pk1_traction_y, pk1_traction_z, traction_x, traction_y, traction_z, normal_traction, normal_jump, tangent_traction, tangent_jump
Controllable:No
Description:Add scalar quantity output for stress and/or strain
material_output_familySpecifies the family of FE shape functions to use for this variable.
C++ Type:MultiMooseEnum
Options:MONOMIAL
Controllable:No
Description:Specifies the family of FE shape functions to use for this variable.
material_output_orderSpecifies the order of the FE shape function to use for this variable.
C++ Type:MultiMooseEnum
Options:CONSTANT, FIRST, SECOND, THIRD, FOURTH, FIFTH, SIXTH, SEVENTH, EIGHTH, NINTH
Controllable:No
Description:Specifies the order of the FE shape function to use for this variable.

Output 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.
diag_save_in_primaryThe displacement diagonal preconditioner terms on the primary side
C++ Type:std::vector<AuxVariableName>
Unit:(no unit assumed)
Controllable:No
Description:The displacement diagonal preconditioner terms on the primary side
diag_save_in_secondaryThe displacement diagonal preconditioner terms on the secondary side
C++ Type:std::vector<AuxVariableName>
Unit:(no unit assumed)
Controllable:No
Description:The displacement diagonal preconditioner terms on the secondary side
save_in_primaryThe displacement residuals on the primary side
C++ Type:std::vector<AuxVariableName>
Unit:(no unit assumed)
Controllable:No
Description:The displacement residuals on the primary side
save_in_secondaryThe displacement residuals on the secondary side
C++ Type:std::vector<AuxVariableName>
Unit:(no unit assumed)
Controllable:No
Description:The displacement residuals on the secondary side

Advanced Parameters

Associated Actions

Available Actions

Shifted Boundary Method App
ShiftedCohesiveZoneActionAction to create Shifted Cohesive Zone Method (SCZM) objects for non-interface-fitted meshes.

(modules/shifted_boundary_method/test/tests/shifted_cohesive_zone/angled_sczm_nonad_salehani.i)

nx = 8
nx_half = '${fparse nx / 2}'
E = 1e3
poi = 0.3
theta = 60
theta_rad = '${fparse theta * pi / 180}'
x0 = 0.5

[GlobalParams]
  displacements = 'disp_x disp_y'
  use_displaced_mesh = false
  sbm_distance_uo = sbm_distance_uo
[]

[Problem]
  extra_tag_vectors = 'ref'
[]

[Mesh]
  [gen]
    type = CartesianMeshGenerator
    dim = 2
    dx = '${x0} ${x0}'
    dy = 1
    ix = '${nx_half} ${nx_half}'
    iy = '${nx}'
    subdomain_id = '1 2'
  []
  [intercepted]
    type = SubdomainInterceptedGenerator
    input = gen
    subdomain_id_inside = 2
    subdomain_id_outside = 1
    lambda = 0.5
    is_domain_inside_surface = true
    signed_dist_function = '(x-${x0})*sin(${theta_rad}) - cos(${theta_rad})*(y-${x0})'
  []
  [break]
    type = BreakMeshByBlockGenerator
    input = intercepted
    split_interface = true
    add_interface_on_two_sides = true
  []
[]

[UserObjects]
  [sbm_distance_uo]
    type = BoundaryShortestDistanceToSurface
    surfaces = signed_dist_func
    boundary = 'Block1_Block2'
    execution_order_group = 0
    execute_on = INITIAL
  []
[]

[Functions]
  [displacement_with_time]
    type = ParsedFunction
    expression = '1e-2*t'
  []
  [signed_dist_func]
    type = ParsedFunction
    expression = '(x-${x0})*sin(${theta_rad}) - cos(${theta_rad})*(y-${x0})'
  []
[]

[AuxVariables]
  [react_x]
  []
[]

[AuxKernels]
  [react_x]
    type = TagVectorAux
    vector_tag = ref
    v = disp_x
    variable = react_x
    remove_variable_scaling = true
  []
[]

[BCs]
  [anchor_x]
    type = DirichletBC
    variable = disp_x
    boundary = left
    value = 0
  []
  [anchor_y]
    type = DirichletBC
    variable = disp_y
    boundary = bottom
    value = 0
  []
  [pull]
    type = FunctionDirichletBC
    variable = disp_x
    boundary = right
    function = displacement_with_time
    preset = false
  []
[]

[Postprocessors]
  [react_x]
    type = NodalSum
    variable = react_x
    boundary = right
  []
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_type'
  petsc_options_value = 'lu superlu_dist'
  automatic_scaling = true
  line_search = none
  nl_max_its = 50
  nl_abs_tol = 1e-12
  nl_rel_tol = 1e-10
  dt = 1
  end_time = 5
[]

[Outputs]
  [csv]
    type = CSV
    precision = 15
  []
[]

[Physics]
  [SolidMechanics]
    [QuasiStatic]
      [all]
        strain = SMALL
        extra_vector_tags = ref
        add_variables = true
        use_automatic_differentiation = false
      []
    []
  []
[]

[Physics/SolidMechanics/ShiftedCohesiveZone]
  [czm]
    use_automatic_differentiation = false
    boundary = Block1_Block2
  []
[]

[Materials]
  [elastic_stress]
    type = ComputeLinearElasticStress
  []
  [elasticity_tensor_1]
    type = ComputeIsotropicElasticityTensor
    poissons_ratio = ${poi}
    youngs_modulus = ${E}
    block = 1
  []
  [elasticity_tensor_2]
    type = ComputeIsotropicElasticityTensor
    poissons_ratio = ${poi}
    youngs_modulus = ${E}
    block = 2
  []
  [czm]
    type = SalehaniIrani3DCTraction
    normal_gap_at_maximum_normal_traction = 0.1
    tangential_gap_at_maximum_shear_traction = 0.1
    # Match the AD law's normal peak and initial tangential stiffness for comparison.
    maximum_normal_traction = 183.939720585721
    maximum_shear_traction = 303.265329856317
    boundary = Block1_Block2
  []
[]

(modules/shifted_boundary_method/test/tests/shifted_cohesive_zone/angled_sczm_nonad_salehani.i)

nx = 8
nx_half = '${fparse nx / 2}'
E = 1e3
poi = 0.3
theta = 60
theta_rad = '${fparse theta * pi / 180}'
x0 = 0.5

[GlobalParams]
  displacements = 'disp_x disp_y'
  use_displaced_mesh = false
  sbm_distance_uo = sbm_distance_uo
[]

[Problem]
  extra_tag_vectors = 'ref'
[]

[Mesh]
  [gen]
    type = CartesianMeshGenerator
    dim = 2
    dx = '${x0} ${x0}'
    dy = 1
    ix = '${nx_half} ${nx_half}'
    iy = '${nx}'
    subdomain_id = '1 2'
  []
  [intercepted]
    type = SubdomainInterceptedGenerator
    input = gen
    subdomain_id_inside = 2
    subdomain_id_outside = 1
    lambda = 0.5
    is_domain_inside_surface = true
    signed_dist_function = '(x-${x0})*sin(${theta_rad}) - cos(${theta_rad})*(y-${x0})'
  []
  [break]
    type = BreakMeshByBlockGenerator
    input = intercepted
    split_interface = true
    add_interface_on_two_sides = true
  []
[]

[UserObjects]
  [sbm_distance_uo]
    type = BoundaryShortestDistanceToSurface
    surfaces = signed_dist_func
    boundary = 'Block1_Block2'
    execution_order_group = 0
    execute_on = INITIAL
  []
[]

[Functions]
  [displacement_with_time]
    type = ParsedFunction
    expression = '1e-2*t'
  []
  [signed_dist_func]
    type = ParsedFunction
    expression = '(x-${x0})*sin(${theta_rad}) - cos(${theta_rad})*(y-${x0})'
  []
[]

[AuxVariables]
  [react_x]
  []
[]

[AuxKernels]
  [react_x]
    type = TagVectorAux
    vector_tag = ref
    v = disp_x
    variable = react_x
    remove_variable_scaling = true
  []
[]

[BCs]
  [anchor_x]
    type = DirichletBC
    variable = disp_x
    boundary = left
    value = 0
  []
  [anchor_y]
    type = DirichletBC
    variable = disp_y
    boundary = bottom
    value = 0
  []
  [pull]
    type = FunctionDirichletBC
    variable = disp_x
    boundary = right
    function = displacement_with_time
    preset = false
  []
[]

[Postprocessors]
  [react_x]
    type = NodalSum
    variable = react_x
    boundary = right
  []
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_type'
  petsc_options_value = 'lu superlu_dist'
  automatic_scaling = true
  line_search = none
  nl_max_its = 50
  nl_abs_tol = 1e-12
  nl_rel_tol = 1e-10
  dt = 1
  end_time = 5
[]

[Outputs]
  [csv]
    type = CSV
    precision = 15
  []
[]

[Physics]
  [SolidMechanics]
    [QuasiStatic]
      [all]
        strain = SMALL
        extra_vector_tags = ref
        add_variables = true
        use_automatic_differentiation = false
      []
    []
  []
[]

[Physics/SolidMechanics/ShiftedCohesiveZone]
  [czm]
    use_automatic_differentiation = false
    boundary = Block1_Block2
  []
[]

[Materials]
  [elastic_stress]
    type = ComputeLinearElasticStress
  []
  [elasticity_tensor_1]
    type = ComputeIsotropicElasticityTensor
    poissons_ratio = ${poi}
    youngs_modulus = ${E}
    block = 1
  []
  [elasticity_tensor_2]
    type = ComputeIsotropicElasticityTensor
    poissons_ratio = ${poi}
    youngs_modulus = ${E}
    block = 2
  []
  [czm]
    type = SalehaniIrani3DCTraction
    normal_gap_at_maximum_normal_traction = 0.1
    tangential_gap_at_maximum_shear_traction = 0.1
    # Match the AD law's normal peak and initial tangential stiffness for comparison.
    maximum_normal_traction = 183.939720585721
    maximum_shear_traction = 303.265329856317
    boundary = Block1_Block2
  []
[]

(modules/shifted_boundary_method/test/tests/grain_boundary_case/grain_boundary_sczm_separate_interfaces.i)

E = 1e3
poi = 0.3
interfaces = 'grain1_grain2 grain1_grain4 grain1_grain5 grain2_grain4 grain3_grain4 grain3_grain5 grain4_grain5'

ny = 6
nx = '${fparse ny*3/2}'

[GlobalParams]
  displacements = 'disp_x disp_y'
  use_displaced_mesh = false
[]

[Problem]
  extra_tag_vectors = 'ref'
[]

[Mesh]
  [boundary_mesh]
    type = FileMeshGenerator
    file = 'grain_boundary_only.msh'
  []

  [gen]
    type = GeneratedMeshGenerator
    dim = 2
    xmin = -0.6
    xmax = 0.6
    ymin = -0.4
    ymax = 0.4
    nx = ${nx}
    ny = ${ny}
  []

  [grain_gen]
    type = SubdomainGrainIDGenerator
    input = gen
    boundary_mesh = boundary_mesh
    lambda = 1
  []

  [break7]
    input = grain_gen
    type = BreakMeshByBlockGenerator
    split_interface = true
    add_interface_on_two_sides = true
  []

  [left_bottom]
    input = break7
    type = ExtraNodesetGenerator
    coord = '-0.6 -0.4'
    new_boundary = 'left_bottom'
  []

  [right_top]
    input = left_bottom
    type = ExtraNodesetGenerator
    coord = '0.6 0.4'
    new_boundary = 'right_top'
  []

  [grain1_grain2]
    type = FileMeshGenerator
    file = 'grain1_grain2.msh'
    save_with_name = 'grain1_grain2'
  []

  [grain1_grain4]
    type = FileMeshGenerator
    file = 'grain1_grain4.msh'
    save_with_name = 'grain1_grain4'
  []

  [grain1_grain5]
    type = FileMeshGenerator
    file = 'grain1_grain5.msh'
    save_with_name = 'grain1_grain5'
  []

  [grain2_grain4]
    type = FileMeshGenerator
    file = 'grain2_grain4.msh'
    save_with_name = 'grain2_grain4'
  []

  [grain3_grain4]
    type = FileMeshGenerator
    file = 'grain3_grain4.msh'
    save_with_name = 'grain3_grain4'
  []

  [grain3_grain5]
    type = FileMeshGenerator
    file = 'grain3_grain5.msh'
    save_with_name = 'grain3_grain5'
  []

  [grain4_grain5]
    type = FileMeshGenerator
    file = 'grain4_grain5.msh'
    save_with_name = 'grain4_grain5'
  []

  add_subdomain_ids = '10'
  add_subdomain_names = 'block_10'
  final_generator = 'right_top'
[]

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

[Physics]
  [SolidMechanics]
    [QuasiStatic]
      [all]
        strain = SMALL
        extra_vector_tags = 'ref'
        use_automatic_differentiation = false
        generate_output = 'stress_xx stress_xy stress_yy stress_zz strain_xx strain_xy strain_yy strain_zz'
      []
    []
  []
[]
[Physics/SolidMechanics/ShiftedCohesiveZone]
  [czm_ik]
    use_automatic_differentiation = false
    boundary = ${interfaces}
    generate_sbm_distance = true
    check_surface_watertightness = true
  []
[]

[Materials]
  [elastic_stress]
    type = ComputeLinearElasticStress
  []
  [elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    poissons_ratio = ${poi}
    youngs_modulus = ${E}
  []
  [interface_traction]
    type = BiLinearMixedModeTraction
    boundary = ${interfaces}
    penalty_stiffness = 2e3
    GI_c = 1e3
    GII_c = 1e2
    normal_strength = 500
    shear_strength = 300
    displacements = 'disp_x disp_y'
    eta = 2.2
    viscosity = 1e-3
  []
[]

[Executioner]
  type = Transient
  # solve_type = FD
  # automatic_scaling = false
  # residual_and_jacobian_together = false
  solve_type = NEWTON
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_type'
  petsc_options_value = 'lu superlu_dist'
  automatic_scaling = true
  # residual_and_jacobian_together = true
  line_search = 'none'

  nl_max_its = 500
  nl_abs_tol = 1e-12
  nl_rel_tol = 1e-50
  dt = 1
  end_time = 5
[]

[Functions]
  [displacement_with_time]
    type = ParsedFunction
    expression = '1e-2*t'
  []
[]

[Kernels]
[]

[AuxVariables]
  [react_x]
  []
[]

[AuxKernels]
  [react_x]
    type = TagVectorAux
    vector_tag = 'ref'
    v = 'disp_x'
    variable = 'react_x'
    remove_variable_scaling = true
  []
[]

[BCs]
  [anchor_x]
    type = DirichletBC
    variable = disp_x
    boundary = 'left'
    value = 0.0
  []
  [anchor_y]
    type = DirichletBC
    variable = disp_y
    boundary = 'bottom'
    value = 0.0
  []
  [displacement_x_right]
    #Anchors the left side against deformation in the x-direction
    type = FunctionDirichletBC
    variable = disp_x
    boundary = 'right'
    function = displacement_with_time
    preset = false
  []
[]

[Postprocessors]
  [react_x]
    type = NodalSum
    variable = 'react_x'
    boundary = 'right'
  []

  [area]
    type = ElementIntegralFunctorPostprocessor
    functor = 1
    block = '1 2'
    execute_on = 'final'
  []
[]

[Outputs]
  exodus = true
  # execute_on = 'final'
  [csv]
    type = CSV
    precision = 15
  []
[]

(modules/shifted_boundary_method/test/tests/grain_boundary_case/grain_boundary_sczm_coarse_lambda1.i)

E = 1e3
poi = 0.3
interfaces = 'grain1_grain2 grain1_grain4 grain1_grain5 grain2_grain4 grain3_grain4 grain3_grain5 grain4_grain5'

ny = 6
nx = '${fparse ny*3/2}'

[GlobalParams]
  displacements = 'disp_x disp_y'
  use_displaced_mesh = false
[]

[Problem]
  extra_tag_vectors = 'ref'
[]

[Mesh]
  [boundary_mesh]
    type = FileMeshGenerator
    file = 'grain_boundary_only.msh'
    save_with_name = boundary_mesh
  []

  [gen]
    type = GeneratedMeshGenerator
    dim = 2
    xmin = -0.6
    xmax = 0.6
    ymin = -0.4
    ymax = 0.4
    nx = ${nx}
    ny = ${ny}
    subdomain_ids = 10
  []

  [grain_gen]
    type = SubdomainGrainIDGenerator
    input = gen
    boundary_mesh = boundary_mesh
    lambda = 1
  []

  [break7]
    input = grain_gen
    type = BreakMeshByBlockGenerator
    split_interface = true
    add_interface_on_two_sides = true
  []

  [left_bottom]
    input = break7
    type = ExtraNodesetGenerator
    coord = '-0.6 -0.4'
    new_boundary = 'left_bottom'
  []

  [right_top]
    input = left_bottom
    type = ExtraNodesetGenerator
    coord = '0.6 0.4'
    new_boundary = 'right_top'
  []

  add_subdomain_ids = '10'
  add_subdomain_names = 'block_10'
  final_generator = 'right_top'
[]

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

[Physics]
  [SolidMechanics]
    [QuasiStatic]
      [all]
        strain = SMALL
        volumetric_locking_correction = true
        extra_vector_tags = 'ref'
        use_automatic_differentiation = false
        generate_output = 'stress_xx stress_xy stress_yy stress_zz strain_xx strain_xy strain_yy strain_zz'
      []
    []
  []
[]
[Physics/SolidMechanics/ShiftedCohesiveZone]
  [czm_ik]
    use_automatic_differentiation = false
    boundary = ${interfaces}
    complete_interface_mesh = boundary_mesh
  []
[]

[Materials]
  [elastic_stress]
    type = ComputeLinearElasticStress
  []
  [elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    poissons_ratio = ${poi}
    youngs_modulus = ${E}
  []
  [interface_traction]
    type = BiLinearMixedModeTraction
    boundary = ${interfaces}
    penalty_stiffness = 2e3
    GI_c = 1e3
    GII_c = 1e2
    normal_strength = 500
    shear_strength = 300
    displacements = 'disp_x disp_y'
    eta = 2.2
    viscosity = 1e-3
  []
[]

[Executioner]
  type = Transient
  # solve_type = FD
  # automatic_scaling = false
  # residual_and_jacobian_together = false
  solve_type = NEWTON
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_type'
  petsc_options_value = 'lu superlu_dist'
  automatic_scaling = true
  # residual_and_jacobian_together = true
  line_search = 'none'

  nl_max_its = 500
  nl_abs_tol = 1e-12
  nl_rel_tol = 1e-50
  dt = 1
  end_time = 5
[]

[Functions]
  [displacement_with_time]
    type = ParsedFunction
    expression = '1e-2*t'
  []
[]

[Kernels]
[]

[AuxVariables]
  [react_x]
  []
[]

[AuxKernels]
  [react_x]
    type = TagVectorAux
    vector_tag = 'ref'
    v = 'disp_x'
    variable = 'react_x'
    remove_variable_scaling = true
  []
[]

[BCs]
  [anchor_x]
    type = DirichletBC
    variable = disp_x
    boundary = 'left'
    value = 0.0
  []
  [anchor_y]
    type = DirichletBC
    variable = disp_y
    boundary = 'bottom'
    value = 0.0
  []
  [displacement_x_right]
    #Anchors the left side against deformation in the x-direction
    type = FunctionDirichletBC
    variable = disp_x
    boundary = 'right'
    function = displacement_with_time
    preset = false
  []
[]

[Postprocessors]
  [react_x]
    type = NodalSum
    variable = 'react_x'
    boundary = 'right'
  []

  [area]
    type = ElementIntegralFunctorPostprocessor
    functor = 1
    block = '1 2 3 4 5 10'
    execute_on = 'final'
  []
[]

[Outputs]
  exodus = true
  # execute_on = 'final'
  [csv]
    type = CSV
    precision = 15
  []
[]

(modules/shifted_boundary_method/test/tests/grain_boundary_case/grain_boundary_sczm_coarse_lambda1.i)

E = 1e3
poi = 0.3
interfaces = 'grain1_grain2 grain1_grain4 grain1_grain5 grain2_grain4 grain3_grain4 grain3_grain5 grain4_grain5'

ny = 6
nx = '${fparse ny*3/2}'

[GlobalParams]
  displacements = 'disp_x disp_y'
  use_displaced_mesh = false
[]

[Problem]
  extra_tag_vectors = 'ref'
[]

[Mesh]
  [boundary_mesh]
    type = FileMeshGenerator
    file = 'grain_boundary_only.msh'
    save_with_name = boundary_mesh
  []

  [gen]
    type = GeneratedMeshGenerator
    dim = 2
    xmin = -0.6
    xmax = 0.6
    ymin = -0.4
    ymax = 0.4
    nx = ${nx}
    ny = ${ny}
    subdomain_ids = 10
  []

  [grain_gen]
    type = SubdomainGrainIDGenerator
    input = gen
    boundary_mesh = boundary_mesh
    lambda = 1
  []

  [break7]
    input = grain_gen
    type = BreakMeshByBlockGenerator
    split_interface = true
    add_interface_on_two_sides = true
  []

  [left_bottom]
    input = break7
    type = ExtraNodesetGenerator
    coord = '-0.6 -0.4'
    new_boundary = 'left_bottom'
  []

  [right_top]
    input = left_bottom
    type = ExtraNodesetGenerator
    coord = '0.6 0.4'
    new_boundary = 'right_top'
  []

  add_subdomain_ids = '10'
  add_subdomain_names = 'block_10'
  final_generator = 'right_top'
[]

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

[Physics]
  [SolidMechanics]
    [QuasiStatic]
      [all]
        strain = SMALL
        volumetric_locking_correction = true
        extra_vector_tags = 'ref'
        use_automatic_differentiation = false
        generate_output = 'stress_xx stress_xy stress_yy stress_zz strain_xx strain_xy strain_yy strain_zz'
      []
    []
  []
[]
[Physics/SolidMechanics/ShiftedCohesiveZone]
  [czm_ik]
    use_automatic_differentiation = false
    boundary = ${interfaces}
    complete_interface_mesh = boundary_mesh
  []
[]

[Materials]
  [elastic_stress]
    type = ComputeLinearElasticStress
  []
  [elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    poissons_ratio = ${poi}
    youngs_modulus = ${E}
  []
  [interface_traction]
    type = BiLinearMixedModeTraction
    boundary = ${interfaces}
    penalty_stiffness = 2e3
    GI_c = 1e3
    GII_c = 1e2
    normal_strength = 500
    shear_strength = 300
    displacements = 'disp_x disp_y'
    eta = 2.2
    viscosity = 1e-3
  []
[]

[Executioner]
  type = Transient
  # solve_type = FD
  # automatic_scaling = false
  # residual_and_jacobian_together = false
  solve_type = NEWTON
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_type'
  petsc_options_value = 'lu superlu_dist'
  automatic_scaling = true
  # residual_and_jacobian_together = true
  line_search = 'none'

  nl_max_its = 500
  nl_abs_tol = 1e-12
  nl_rel_tol = 1e-50
  dt = 1
  end_time = 5
[]

[Functions]
  [displacement_with_time]
    type = ParsedFunction
    expression = '1e-2*t'
  []
[]

[Kernels]
[]

[AuxVariables]
  [react_x]
  []
[]

[AuxKernels]
  [react_x]
    type = TagVectorAux
    vector_tag = 'ref'
    v = 'disp_x'
    variable = 'react_x'
    remove_variable_scaling = true
  []
[]

[BCs]
  [anchor_x]
    type = DirichletBC
    variable = disp_x
    boundary = 'left'
    value = 0.0
  []
  [anchor_y]
    type = DirichletBC
    variable = disp_y
    boundary = 'bottom'
    value = 0.0
  []
  [displacement_x_right]
    #Anchors the left side against deformation in the x-direction
    type = FunctionDirichletBC
    variable = disp_x
    boundary = 'right'
    function = displacement_with_time
    preset = false
  []
[]

[Postprocessors]
  [react_x]
    type = NodalSum
    variable = 'react_x'
    boundary = 'right'
  []

  [area]
    type = ElementIntegralFunctorPostprocessor
    functor = 1
    block = '1 2 3 4 5 10'
    execute_on = 'final'
  []
[]

[Outputs]
  exodus = true
  # execute_on = 'final'
  [csv]
    type = CSV
    precision = 15
  []
[]

(modules/shifted_boundary_method/test/tests/grain_boundary_case/grain_boundary_sczm_coarse_lambda1.i)

E = 1e3
poi = 0.3
interfaces = 'grain1_grain2 grain1_grain4 grain1_grain5 grain2_grain4 grain3_grain4 grain3_grain5 grain4_grain5'

ny = 6
nx = '${fparse ny*3/2}'

[GlobalParams]
  displacements = 'disp_x disp_y'
  use_displaced_mesh = false
[]

[Problem]
  extra_tag_vectors = 'ref'
[]

[Mesh]
  [boundary_mesh]
    type = FileMeshGenerator
    file = 'grain_boundary_only.msh'
    save_with_name = boundary_mesh
  []

  [gen]
    type = GeneratedMeshGenerator
    dim = 2
    xmin = -0.6
    xmax = 0.6
    ymin = -0.4
    ymax = 0.4
    nx = ${nx}
    ny = ${ny}
    subdomain_ids = 10
  []

  [grain_gen]
    type = SubdomainGrainIDGenerator
    input = gen
    boundary_mesh = boundary_mesh
    lambda = 1
  []

  [break7]
    input = grain_gen
    type = BreakMeshByBlockGenerator
    split_interface = true
    add_interface_on_two_sides = true
  []

  [left_bottom]
    input = break7
    type = ExtraNodesetGenerator
    coord = '-0.6 -0.4'
    new_boundary = 'left_bottom'
  []

  [right_top]
    input = left_bottom
    type = ExtraNodesetGenerator
    coord = '0.6 0.4'
    new_boundary = 'right_top'
  []

  add_subdomain_ids = '10'
  add_subdomain_names = 'block_10'
  final_generator = 'right_top'
[]

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

[Physics]
  [SolidMechanics]
    [QuasiStatic]
      [all]
        strain = SMALL
        volumetric_locking_correction = true
        extra_vector_tags = 'ref'
        use_automatic_differentiation = false
        generate_output = 'stress_xx stress_xy stress_yy stress_zz strain_xx strain_xy strain_yy strain_zz'
      []
    []
  []
[]
[Physics/SolidMechanics/ShiftedCohesiveZone]
  [czm_ik]
    use_automatic_differentiation = false
    boundary = ${interfaces}
    complete_interface_mesh = boundary_mesh
  []
[]

[Materials]
  [elastic_stress]
    type = ComputeLinearElasticStress
  []
  [elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    poissons_ratio = ${poi}
    youngs_modulus = ${E}
  []
  [interface_traction]
    type = BiLinearMixedModeTraction
    boundary = ${interfaces}
    penalty_stiffness = 2e3
    GI_c = 1e3
    GII_c = 1e2
    normal_strength = 500
    shear_strength = 300
    displacements = 'disp_x disp_y'
    eta = 2.2
    viscosity = 1e-3
  []
[]

[Executioner]
  type = Transient
  # solve_type = FD
  # automatic_scaling = false
  # residual_and_jacobian_together = false
  solve_type = NEWTON
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_type'
  petsc_options_value = 'lu superlu_dist'
  automatic_scaling = true
  # residual_and_jacobian_together = true
  line_search = 'none'

  nl_max_its = 500
  nl_abs_tol = 1e-12
  nl_rel_tol = 1e-50
  dt = 1
  end_time = 5
[]

[Functions]
  [displacement_with_time]
    type = ParsedFunction
    expression = '1e-2*t'
  []
[]

[Kernels]
[]

[AuxVariables]
  [react_x]
  []
[]

[AuxKernels]
  [react_x]
    type = TagVectorAux
    vector_tag = 'ref'
    v = 'disp_x'
    variable = 'react_x'
    remove_variable_scaling = true
  []
[]

[BCs]
  [anchor_x]
    type = DirichletBC
    variable = disp_x
    boundary = 'left'
    value = 0.0
  []
  [anchor_y]
    type = DirichletBC
    variable = disp_y
    boundary = 'bottom'
    value = 0.0
  []
  [displacement_x_right]
    #Anchors the left side against deformation in the x-direction
    type = FunctionDirichletBC
    variable = disp_x
    boundary = 'right'
    function = displacement_with_time
    preset = false
  []
[]

[Postprocessors]
  [react_x]
    type = NodalSum
    variable = 'react_x'
    boundary = 'right'
  []

  [area]
    type = ElementIntegralFunctorPostprocessor
    functor = 1
    block = '1 2 3 4 5 10'
    execute_on = 'final'
  []
[]

[Outputs]
  exodus = true
  # execute_on = 'final'
  [csv]
    type = CSV
    precision = 15
  []
[]

Select the interface geometry
Define the traction law
Non-AD tangent definition
Input Parameters
Associated Actions
Available Actions