Private API

InterfaceError

A type for a customized error that is thrown when a material model is not implemented correctly.

Fields

• type::DataType: Type that is used.
• func::String: Function that is used.
• hint::String: Optional explanation of how the interface method can be defined.

NaNError

A type for a customized error that is thrown when NaN values are detected in the internal force density field after the force density evaluation.

Fields

• time::Float64: Simulation time when NaNs were detected.
• step::Int: Simulation step when NaNs were detected.

HaloExchange

A type used for communication between body chunks. This type is used with both MPI and multithreaded simulations. Note that the tag is only used with MPI and has the value 0 in multithreaded simulations.

Fields

• tag::Int: Tag used for the MPI sending and receiving commands.
• src_chunk_id::Int: Index of the chunk that sends information.
• dest_chunk_id::Int: Index of the chunk that receives information.
• src_idxs::Vector{Int}: Indices of the points in the source chunk that send information.
• dest_idxs::Vector{Int}: Indices of the points in the destination chunk that receive information.

JobOptions{F,V}

A type that contains the options of a job.

Type Parameters

• F: Type for fields of simulation.
• V: Type for basename of vtk-files.

Fields

• export_allowed::Bool: Specify if data is exported for the job.
• root::String: Path of the folder where all data is saved.
• vtk::String: Path of the folder where vtk-files are saved.
• logfile::String: Complete path of the logfile.
• freq::Int: Frequency of exported time steps.
• fields::F: Exported fields of the job.
• vtk_filebase::V: Basename of exported vtk-files.

MPIHaloInfo

A type providing information about the halo points of body chunks.

Fields

• point_ids::Dict{Int,Vector{Int}}: Indices of local and halo points of each chunk.
• halos_points::Dict{Int,Vector{Int}}: Indices of halo points of each chunk.
• localizers::Dict{Int,Dict{Int,Int}}: Localizes global indices to local indices in the corresponding chunks.
• hidxs_by_src::Dict{Int,Dict{Int,Vector{Int}}}: Dict specifying for each chunk the indices of halo points but depending on the chunk they belong to.

MPIBodyDataHandler

A type for handling information and communication for a body chunk in MPI.

Type Parameters

• Sys: Type of system of the body.
• M: Type of the material model.
• P: Material parameter type.
• S: Storage type of the system.
• Bufs: Type for buffers.

Fields

• chunk::BodyChunk{Sys,M,P,S}: Body chunk of the body.
• n_halo_fields::Int: Number of chunks in communication with this chunk due to existing halo points.
• lth_exs_send::Vector{HaloExchange}: Local-to-halo-exchanges of the body chunk that send information.
• lth_exs_recv::Vector{HaloExchange}: Local-to-halo-exchanges of the body chunk that receive information.
• htl_exs_send::Vector{HaloExchange}: Halo-to-local-exchanges of the body chunk that send information.
• htl_exs_recv::Vector{HaloExchange}: Halo-to-local-exchanges of the body chunk that receive information.
• lth_send_bufs::Vector{Bufs}: Buffers for local-to-halo-exchanges that send information.
• lth_recv_bufs::Vector{Bufs}: Buffers for local-to-halo-exchanges that receive information.
• htl_send_bufs::Vector{Bufs}: Buffers for halo-to-local-exchanges that send information.
• htl_recv_bufs::Vector{Bufs}: Buffers for halo-to-local-exchanges that receive information.
• field_to_buf::Dict{Symbol,Int}: Dict specifying the index in the buffer of each field of the simulation.
• lth_reqs::Vector{Vector{MPI.Request}}: Pre-allocated buffer for requests for local-to-halo-exchanges.
• htl_reqs::Vector{Vector{MPI.Request}}: Pre-allocated buffer for requests for halo-to-local-exchanges.

SingleParamChunk

Type for a body chunk of a body with only one material parameter set.

MultiParamChunk

Type for a body chunk of a body with multiple material parameter sets.

ParameterHandler

A type used to manage multiple point parameters defined for the same body. It is used to assign different point parameters to the points of a body.

Type Parameters

• P<:AbstractPointParameters: Point parameter type.

Fields

• parameters::Vector{P}: All parameter sets defined in the simulation.
• point_mapping::Vector{Int}: Vector assigning the related parameter set to each material point.

ThreadsBodyDataHandler

A type for handling all data of a body in multithreading simulations.

Type Parameters

• Sys: Type of system of the body.
• M: Type of the material model.
• P: Material parameter type.
• S: Storage type of the system.

Fields

• n_chunks::Int: Number of chunks of the body.
• chunks::Vector{BodyChunk{Sys,M,P,S}}: All body chunks of the body.
• lth_exs::Vector{Vector{HaloExchange}}: All local-to-halo-exchanges of each body chunk of the body.
• htl_exs::Vector{Vector{HaloExchange}}: All halo-to-local-exchanges of each body chunk of the body.

ThreadsMultibodyDataHandler

A type for handling all data of multiple bodies in multithreading simulations.

Type Parameters

• BDH: Body data handler type.
• PC: Position cache type.
• VC: Volume cache type.

Fields

• n_bodies::Int: Number of bodies in the simulation.
• body_dhs::BDH: Tuple containing all body data handlers.
• body_names::Vector{Symbol}: Names of the bodies.
• body_idxs::Dict{Symbol,Int}: Names of bodies assigned to their indices.
• srf_contacts::Vector{ShortRangeForceContact}: All short range force contacts of this simulation.
• position_caches::PC: Positions of all points of all bodies (should be of type Vector{Matrix{Float64}}).
• volume_caches::VC: Volumes of all points of all bodies (should be of type Vector{Vector{Float64}}).

BodyChunk{System, Material, Params, Storage}

A type that contains all data of a body chunk. For parallel simulations, the body is divided into multiple chunks. Each BodyChunk instance contains all necessary information for the simulation on this specific chunk. This type is used for multithreading and MPI.

Type Parameters

• System<:AbstractSystem: Type of the system.
• Material<:AbstractMaterial: Type of the material model of the system.
• Params<:AbstractParameterSetup: Material parameters of the points in the body chunk.
• Storage<:AbstractStorage: Storage of all information that changes during the simulation.

Fields

• body_name::Symbol: Name of the body in multibody simulations.
• system::System: System with all information that is known before the simulation.
• mat::Material: Material model of the system.
• paramsetup::Params: Material parameters of the points in the body chunk.
• storage::Storage: Storage of all information that changes during the simulation.
• psets::Dict{Symbol,Vector{Int}}: Point sets of the chunk with local indices.
• sdbcs::Vector{SingleDimBC}: Single dimension boundary conditions.
• pdsdbcs::Vector{PosDepSingleDimBC}: Position dependent single dimension boundary conditions.
• cells::Vector{MeshCell{VTKCellType,Tuple{Int64}}}: Cells for vtk export.

Bond

Type that describes a bond of two points in a peridynamics body.

Fields

• neighbor::Int: The index of the neighbor point with which the bond is formed.
• length::Float64: The length of the bond.
• fail_permit::Bool: Describes whether failure is allowed or not for this bond.

BondSystem{Correction}

A type for a system for all peridynamic formulations that work with just bonds of two points.

Type Parameters

• Correction<:AbstractCorrection: Applied surface correction.

Fields

• position::Matrix{Float64}: Positions of all points of the system.
• volume::Vector{Float64}: Volumes of the points of the system.
• bonds::Vector{Bond}: Vector containing all bonds of the bond system.
• n_neighbors::Vector{Int}: Number of neighbors for each point of the system.
• bond_ids::Vector{UnitRange{Int}}: Range of the bonds vector containing bonds of considered point.
• correction::Correction: Applied surface correction.
• chunk_handler::ChunkHandler: Type to handle the chunks for the simulation. See ChunkHandler.

ChunkHandler

A type to handle a body chunk and its communication to other chunks.

Fields

• n_loc_points::Int: Number of local points that belong to the body chunk.
• point_ids::Vector{Int}: Indices of all local and halo points of the chunk.
• loc_points::UnitRange{Int}: Indices of local points of the chunk.
• halo_points::Vector{Int}: Indices of halo points of the chunk.
• hidxs_by_src::Dict{Int,UnitRange{Int}}: Dict specifying the indices of halo Points depending on the body chunk they belong to. So body_chunk => indices, with indices being the indices of the halo points in point_ids.
• localizer::Dict{Int,Int}: Localizes global indices to local indices in this chunk.

PointDecomposition

A type that describes how a body is divided into multiple body chunks.

Fields

• n_chunks::Int: Number of body chunks.
• decomp::Vector{UnitRange{Int}}: Indices of the points belonging to each chunk.
• point_src::Dict{Int,Int}: Dict that assigns all point indices to the chunk they belong to.

TwoNeighborInteraction

Type for two-neighbor interactions.

Fields

• oni_j::Int: One-neighbor interaction of considered point with point j.
• oni_k::Int: One-neighbor interaction of considered point with point k.
• surface::Float64: Surface spread by this two-neighbor interaction.

ThreeNeighborInteraction

Type for three-neighbor interactions.

Fields

• oni_j::Int: One-neighbor interaction of considered point with point j.
• oni_k::Int: One-neighbor interaction of considered point with point k.
• oni_l::Int: One-neighbor interaction of considered point with point l.
• volume::Float64: Volume spread by this three-neighbor interaction.

InteractionSystem

A peridynamic system type that is mainly designed for continuum-kinematics-inspired peridynamics [JMS19].

Fields

• position::Matrix{Float64}: Positions of all points of the system.
• one_nis::Vector{Bond}: Vector containing all one-neighbor interactions (bonds) of the system.
• two_nis::Vector{TwoNeighborInteraction}: Vector containing all two-neighbor interactions of the system.
• three_nis::Vector{ThreeNeighborInteraction}: Vector containing all three-neighbor interactions of the system.
• volume::Vector{Float64}: Volumes of the points of the system.
• volume_one_nis::Vector{Float64}: Effective volumes of one-neighbor interactions.
• volume_two_nis::Vector{Float64}: Effective volumes of two-neighbor interactions.
• volume_three_nis::Vector{Float64}: Effective volumes of three-neighbor interactions.
• n_one_nis::Vector{Int}: Number of one-neighbor interactions for each point of the system.
• n_two_nis::Vector{Int}: Number of two-neighbor interactions for each point of the system.
• n_three_nis::Vector{Int}: Number of three-neighbor interactions for each point of the system.
• one_ni_idxs::Vector{UnitRange{Int}}: Range of the one-neighbor interactions vector containing interactions of considered point.
• two_ni_idxs::Vector{UnitRange{Int}}: Range of the two-neighbor interactions vector containing interactions of considered point.
• three_ni_idxs::Vector{UnitRange{Int}}: Range of the three-neighbor interactions vector containing interactions of considered point.
• chunk_handler::ChunkHandler: Type to handle the chunks for the simulation. See ChunkHandler.

PointSetsPreCrack

Type describing a predefined crack in a peridynamic body.

Fields

• set_a::Symbol: Point set containing points on one side of the crack.
• set_b::Symbol: Point set with points on other side of the crack.
• filter_bonds::Bool: If true, the involved bonds are filtered out so no damage is present at the beginning of the simulation. Else, all involved bonds are marked broken from the beginning.

SingleDimBC{F}

Type for a boundary condition in a single dimension for a peridynamic simulation.

Type Parameters

• F<:Function: Time dependent function which describes the boundary condition.

Fields

• fun::F: Time dependent function which describes the boundary condition.
• field::Symbol: Field of the condition (e.g. velocity, force density).
• point_set::Symbol: Point set on which the condition is applied.
• dim::UInt8: Dimension in which the condition is applied.

PosSingleDimBC{F}

Type for a position dependent boundary condition in a single dimension for a peridynamic simulation.

Type Parameters

• F<:Function: A position dependent function which describes the boundary condition, not time dependent.

Fields

• fun::F: Position dependent function which describes the boundary condition.
• field::Symbol: Field of the condition (e.g. velocity, force density).
• point_set::Symbol: Point set on which the condition is applied.
• dim::UInt8: Dimension in which the condition is applied.

PosDepSingleDimBC{F}

Type for a position dependent boundary condition in a single dimension for a peridynamic simulation.

Type Parameters

• F<:Function: Position and time dependent function which describes the boundary condition.

Fields

• fun::F: Position and time dependent function which describes the boundary condition.
• field::Symbol: Field of the condition (e.g. velocity, force density).
• point_set::Symbol: Point set on which the condition is applied.
• dim::UInt8: Dimension in which the condition is applied.

BBPointParameters

Point parameters of the bond-based material: the discretization parameters, the elastic parameters with the Poisson's ratio of bond-based peridynamics, the bond constant bc and the parameters of the damage model. GBBMaterial uses them as they are, DHBBMaterial inherits them and halves the bond constant.

Computed together:

• (; δ, rho) = get_discretization_params(; horizon, rho)
• (; E, nu, G, K, λ, μ) = get_elastic_params_bb(; E, nu, G, K, lambda, mu)

Keywords of material!: horizon, rho, E, nu, G, K, lambda, mu.

DHBBPointParameters

Point parameters of the dual-horizon bond-based material: the BBPointParameters with half of the bond constant, because every bond is visited from both of its points.

Computed together:

• (; δ, rho) = get_discretization_params(; horizon, rho)
• (; E, nu, G, K, λ, μ) = get_elastic_params_bb(; E, nu, G, K, lambda, mu)

Keywords of material!: horizon, rho, E, nu, G, K, lambda, mu.

OSBPointParameters

Point parameters of the ordinary state-based material: exactly the StandardParameters.

Computed together:

• (; δ, rho) = get_discretization_params(; horizon, rho)
• (; E, nu, G, K, λ, μ) = get_elastic_params(; E, nu, G, K, lambda, mu)

Keywords of material!: horizon, rho, E, nu, G, K, lambda, mu.

CKIPointParameters

Type containing the material parameters for a continuum-kinematics-inspired peridynamics model.

Fields

• δ::Float64: Horizon.
• rho::Float64: Density.
• E::Float64: Young's modulus.
• nu::Float64: Poisson's ratio.
• G::Float64: Shear modulus.
• K::Float64: Bulk modulus.
• λ::Float64: 1st Lamé parameter.
• μ::Float64: 2nd Lamé parameter.
• Gc::Float64: Critical energy release rate.
• εc::Float64: Critical strain.
• C1::Float64: Material constant for one-neighbor interactions.
• C2::Float64: Material constant for two-neighbor interactions.
• C3::Float64: Material constant for three-neighbor interactions.

BACPointParameters

Type containing the material parameters for a peridynamics model using the bond-associated correspondence formulation of Chen and Spencer.

Fields

• δ::Float64: Horizon.
• δb::Float64: Bond-associated horizon.
• rho::Float64: Density.
• E::Float64: Young's modulus.
• nu::Float64: Poisson's ratio.
• G::Float64: Shear modulus.
• K::Float64: Bulk modulus.
• λ::Float64: 1st Lamé parameter.
• μ::Float64: 2nd Lamé parameter.
• Gc::Float64: Critical energy release rate.
• εc::Float64: Critical strain.
• bc::Float64: Bond constant.

CPointParameters

Point parameters of the correspondence formulation. In addition to the standard parameters, they carry the parameters of the constitutive model and the elastic stiffness tensor C, which the ZEMWan stabilization needs. CRMaterial uses them as well.

Computed together:

• (; δ, rho) = get_discretization_params(; horizon, rho)
• (; E, nu, G, K, λ, μ) = get_elastic_params(; E, nu, G, K, lambda, mu)

Keywords of material!: horizon, rho, E, nu, G, K, lambda, mu.

RKCPointParameters

Point parameters of the reproducing-kernel correspondence family: the standard parameters plus the parameters of the constitutive model the material carries. RKCRMaterial uses them as well.

Computed together:

• (; δ, rho) = get_discretization_params(; horizon, rho)
• (; E, nu, G, K, λ, μ) = get_elastic_params(; E, nu, G, K, lambda, mu)

Keywords of material!: horizon, rho, E, nu, G, K, lambda, mu.

SingleDimIC

Type for an initial condition in a single dimension for a peridynamic simulation.

Fields

• value::Float64: Value of the condition.
• field::Symbol: Field of the condition (e.g. velocity, force density).
• point_set::Symbol: Point set on which the condition is applied.
• dim::UInt8: Dimension in which the condition is applied.

PosDepSingleDimIC{F}

Type for a position dependent initial condition in a single dimension for a peridynamic simulation.

Type Parameters

• F<:Function: Position dependent function which describes the initial condition.

Fields

• fun::F: Position dependent function which describes the initial condition.
• field::Symbol: Field of the condition (e.g. velocity, force density).
• point_set::Symbol: Point set on which the condition is applied.
• dim::UInt8: Dimension in which the condition is applied.

ShortRangeForceContact

A type for contact simulations with the short range forces algorithm.

Type Parameters

• N: Neighborhood search object used by PointNeighbors.jl.

Fields

• body_id_a::Symbol: Index of a body of the contact.
• body_id_b::Symbol: Index of a body of the contact.
• radius::Float64: Search radius for contact.
• penalty_factor::Float64: Penalty factor for the contact simulation.
• nhs::N: Neighborhood search object used by PointNeighbors.jl.

failure_permit!(body, set_name, fail_permit)

Set the failure permission for points of the set set_name of a body.

Arguments

• body::AbstractBody: Body where the failure permission will be set.
• set_name::Symbol: The name of a point set of this body.
• fail_permit::Bool: If true, failure is allowed, and if false then no bonds of this point are allowed to break during the simulation.

Throws

• Error if the body does not contain a set with set_name.

get_frac_params(dmgmodel, δ, K; kwargs...)

Read or calculate the fracture parameters of a damage model from the fracture keywords of material!. This function has to be defined when creating a new damage model. Otherwise, a default method returns an empty named tuple (; ).

Every fracture keyword of material! is passed to this function, and a keyword the user did not specify arrives as nothing — that is how a damage model decides which keywords it accepts and how it converts them into each other. A method declares only the keywords its model reads and collects everything else in kwargs..., so a model that has no notion of a critical strain never has to mention epsilon_c:

function Peridynamics.get_frac_params(::MyDamage, δ, K; Gc=nothing, kwargs...)
    isnothing(Gc) && return (; Gc=0.0, εc=0.0)
    return (; Gc, εc=sqrt(5.0 * Gc / (9.0 * K * δ)))
end

The fracture keywords material! accepts are currently fixed to Gc and epsilon_c; a damage model cannot register its own keywords yet, and a keyword its method does not read is ignored.

Arguments

• dmgmodel::AbstractDamageModel: The damage model
• δ::Float64: Horizon
• K::Float64: Bulk modulus

Keywords

• Gc: Critical energy release rate, or nothing if not specified
• epsilon_c: Critical strain, or nothing if not specified

set_failure_permissions!(body, set_name, params)

Grant or prohibit failure permission depending on the submitted fracture parameters by calling failure_permit!.

If fracture parameters are found, failure is allowed. If no fracture parameters are found, failure is not allowed.

has_fracture(mat, params)

Return true if at least one fracture parameter is set !=0 in params and the system therefore is supposed to have failure allowed or return false if not.

check_pos_and_vol(n_points, position, volume)

Check if the positions and volumes for the points are correctly specified in the fields of a Body.

pre_submission_check(body::Body, time_solver; body_in_multibody_setup::Bool=false)
pre_submission_check(ms::AbstractMultibodySetup, time_solver)

Check if necessary material parameters and conditions are defined when defining a Job, and if every condition of the body can be applied by time_solver.

check_conditions_applicable(body::Body, time_solver)

Check that every condition of body is one that time_solver actually applies.

Prescribed displacement conditions are applied by apply_incr_boundary_conditions!, which only the NewtonKrylov solver calls. With any other solver the simulation would run to completion with none of the prescribed displacement applied.

get_paramsetup(body::AbstractBody, ::AbstractChunkHandler, ::SingleParamChunk)
get_paramsetup(body::AbstractBody, ch::AbstractChunkHandler, ::MultiParamChunk)

Return the parameters of a BodyChunk if only one parameter set is defined for the corresponding Body or the parameter handler if multiple parameter sets are defined.

get_params(paramhandler::ParameterHandler, point_id::Int)
get_params(params::AbstractPointParameters, ::Int)
get_params(chunk::BodyChunk, point_id::Int)

Return parameters of a specific point with index point_id of a Body with parameters params or parameter handler paramhandler or of the body chunk.

check_if_sets_intersect(point_sets, key_a, key_b)

Throw error if two sets chosen for precrack! have common points.

check_if_set_is_defined(point_sets, set_name)

Throw an error if no point set set_name is found in the dictionary point_sets.

find_points(f, position)

Find all points whose positions meet function f.

The function f accepts only one positional argument and will be used in a findall call. Depending on the argument name, a different input will be processed. See point_set!.

apply_precracks!(chunk, body)

Apply all predefined cracks for chunk by calling apply_precrack! for each crack.

apply_precrack!(chunk, body, crack)

Apply the predefined crack of the body for the considered chunk by breaking the concerned bonds.

point_sets_intersect(point_sets, key_a, key_b)

Return true if point sets key_a and key_b out of the dictionary point_sets have at least one common point, return false if they do not.

invreg(M::StaticMatrix{N,N,T}, λ::Real, β::Real) where {N,T}

Computes the regularized pseudo-inverse of a square static matrix M using a combination of Tikhonov regularization in the SVD domain and truncated singular value regularization.

Arguments

• M::StaticMatrix{N,N,T}: The square N×N static matrix with element type T to be inverted.
• λ::Real: Relative Tikhonov regularization parameter (dimensionless, non-negative). Controls the smoothing strength applied as $\lambda_{\text{eff}} = \lambda \sigma_{\max}$, where $\sigma_{\max}$ is the largest singular value of M.
• β::Real: Relative SVD truncation parameter (dimensionless, non-negative). Defines the cutoff threshold as $\beta_{\text{eff}} = \beta \sigma_{\max}$ for excluding small singular values.

Returns

• Minv::StaticMatrix{N,N,T}: The regularized pseudo-inverse of the input matrix M.

Regularization Techniques

The function applies two complementary regularization strategies:

1. SVD-based Tikhonov Regularization: For each singular value $\sigma_i$, the inverse is computed as $\sigma_i/(\sigma_i^2 + \lambda_{\text{eff}}^2)$, which smoothly dampens the contribution of small singular values without completely removing them.

2. Truncated SVD: Singular values below the threshold $\beta_{\text{eff}}$ are completely excluded by setting their contribution to zero, preventing numerical instability from near-zero singular values.

Parameter Selection Guidelines

• $λ$ (Tikhonov parameter):

  • Well-conditioned matrices: $\lambda = 0$ (no Tikhonov regularization, recommended default)
  • Mild regularization: $\lambda \in [0, 10^{-12}]$ (scale-invariant gentle smoothing)
  • Moderate regularization: $\lambda \in [10^{-12}, 10^{-4}]$ (for moderately ill-conditioned problems)
  • Note: Values $\lambda > 10^{-4}$ may introduce noticeable bias in the solution

• $β$ (truncation parameter): Primary regularization mechanism, less sensitive than $\lambda$.

  • Well-conditioned matrices: $\beta \in [\sqrt{\epsilon}, 10^{-6}]$ (remove numerical noise, recommended default)
  • Moderately ill-conditioned: $\beta \in [10^{-6}, 10^{-4}]$ (moderate truncation)
  • Severely ill-conditioned: $\beta \in [10^{-4}, 10^{-2}]$ (aggressive truncation)

See also the single-parameter method invreg(M, ε), which regularizes each singular value on its own and is therefore free of bias for the singular values that do not need regularizing.

invreg(M::StaticMatrix{N,N,T}, ε::Real) where {N,T}

Computes the regularized inverse of a square static matrix M with an adaptive, per-singular-value Tikhonov damping that is only active where it is needed.

Arguments

• M::StaticMatrix{N,N,T}: The square N×N static matrix with element type T to be inverted.
• ε::Real: Relative singular value floor (dimensionless, non-negative). The smallest singular value that is still trusted, expressed as a fraction of the largest one: $\sigma_f = \varepsilon \sigma_{\max}$.

Returns

• Minv::StaticMatrix{N,N,T}: The regularized inverse of the input matrix M.

Regularization

Every singular value gets its own damping parameter

\[\lambda_i = \max(0, \sigma_f - \sigma_i) , \qquad \sigma_i^{\mathrm{inv}} = \frac{\sigma_i}{\sigma_i^2 + \lambda_i^2}\]

so a singular value above the floor is damped with $\lambda_i = 0$, which is to say it is inverted exactly, while one below the floor is damped smoothly towards zero:

This bounds the inverse: the damped branch peaks at $\sigma_i = \sigma_f/\sqrt{2}$, so

\[\|M^{-1}\|_2 \leq \frac{1}{(2\sqrt{2}-2)\,\sigma_f} = \frac{1}{(2\sqrt{2}-2)\,\varepsilon\,\sigma_{\max}} \approx \frac{1.21}{\varepsilon\,\sigma_{\max}}\]

no matter how rank deficient M becomes, and it stays bounded without a jump: the mapping is $C^1$ across $\sigma_i = \sigma_f$, because the derivative of $\sigma/(\sigma^2 + (\sigma_f - \sigma)^2)$ at $\sigma_f^-$ is $-1/\sigma_f^2$, which is the derivative of $1/\sigma$ at $\sigma_f^+$. A hard truncation, as used by the two-parameter method above, is only $C^0$: as a singular value crosses the cutoff — which is what happens when bonds drop out of a family — the inverse jumps, and a time integrator sees that jump as a force impulse.

Parameter Selection Guidelines

• $\varepsilon = 10^{-3}$ (recommended default): a 3D moment matrix of an intact point has $\sigma_{\min}/\sigma_{\max} \gtrsim 0.05$, so the regularization stays inactive in the bulk and only starts to act once enough bonds have failed to depress the conditioning by roughly 50x.
• $\varepsilon = 10^{-2}$: more conservative, and also acts on the intact matrices of thin structures with few points across the thickness.
• $\varepsilon = 0$: no regularization, the exact pseudo-inverse.

get_invreg_params(epsilon, lambda, beta)

Resolves the moment matrix regularization keywords of a material of the RKC family into the triple (ε, λ, β) that is stored in it. Exactly one of the two regularizations is active: ε > 0 selects the adaptive invreg with a singular value floor, ε == 0 selects the legacy Tikhonov-and-truncation invreg with λ and β. A keyword that was not given is nothing: without any keyword the floor is 1e-3, with a legacy keyword the missing legacy parameter takes its old default (λ = 0, β = sqrt(eps())).

update_sim_success_from_log!(study::Study)

Read the study.logfile and update study.sim_success flags according to the last recorded status for each job. This allows resuming processing or submission after an interrupted run.

get_vector(M, i)

Return column i of the storage field M as a SVector{3}, e.g. the position of point i of a PointVector field.

Example

u = Peridynamics.get_vector(storage.displacement, i)

See also update_vector!, get_vector_diff.

get_vector_diff(M, i, j)

Return column j - column i of the storage field M as a SVector{3}, without building the two columns first. This is the bond vector of a bond from point i to point j:

ΔXij = Peridynamics.get_vector_diff(system.position, i, j)  # initial bond vector
Δxij = Peridynamics.get_vector_diff(storage.position, i, j) # current bond vector

update_vector!(M, i, V)

Write the SVector{3} V into column i of the storage field M, overwriting what is there. The inverse of get_vector.

Note that a force density is accumulated over the bonds of a point, so it is written with update_add_vector! and not with this function.

update_add_vector!(M, i, V)

Add the SVector{3} V to column i of the storage field M. This is how a force density is accumulated inside force_density_point!, where every bond of a point contributes a share.

Example

Peridynamics.update_add_vector!(storage.b_int, i, b)

get_tensor(M, i)

Return column i of the storage field M as a SMatrix{3,3}, e.g. the deformation gradient of point or bond i of a PointTensor or BondTensor field. The column is read in column-major order, which is how update_tensor! writes it.

Example

F = Peridynamics.get_tensor(storage.defgrad, i)

See also update_tensor!, get_vector.

update_tensor!(M, i, T)

Write the SMatrix{3,3} T into column i of the storage field M, in column-major order. The inverse of get_tensor.

Example

Peridynamics.update_tensor!(storage.defgrad, i, F)

get_sym_tensor(M, i)

Return column i of the storage field M as a symmetric SMatrix{3,3}, e.g. the plastic strain of point or bond i of a PointSymTensor or BondSymTensor field. The column holds the six independent components in Voigt order, (11, 22, 33, 23, 13, 12), which is how update_sym_tensor! writes it.

Note that a symmetric field has six rows, not nine, so get_tensor must not be used on it.

Example

εᵖ = Peridynamics.get_sym_tensor(state.bond_plastic_strain, idx)

See also update_sym_tensor!, get_tensor.

update_sym_tensor!(M, i, T)

Write the symmetric SMatrix{3,3} T into column i of the storage field M, as the six independent components in Voigt order, (11, 22, 33, 23, 13, 12). The inverse of get_sym_tensor.

Only the upper triangle of T is read, and get_sym_tensor mirrors it back. A tensor that is symmetric only up to round-off therefore comes back changed by that round-off, and a tensor with a real skew part is not symmetrized — its lower triangle is silently discarded.

Example

Peridynamics.update_sym_tensor!(state.bond_plastic_strain, idx, εᵖ + Δεᵖ)

sym_eigvals(A)

Return the three eigenvalues of the symmetric SMatrix{3,3} A in descending order, as an SVector{3}. They are computed in closed form (Smith, 1961) rather than by an eigendecomposition, so this allocates nothing and needs no eigenvectors.

Only the upper triangle of A is read.

See also hencky_and_invstretch.

hencky_and_invstretch(C)

Return the material logarithmic (Hencky) strain $\boldsymbol{\varepsilon} = \ln \boldsymbol{U} = \frac{1}{2} \ln \boldsymbol{C}$ and the inverse right stretch $\boldsymbol{U}^{-1} = \boldsymbol{C}^{-1/2}$ of the right Cauchy-Green tensor $\boldsymbol{C} = \boldsymbol{F}^T \boldsymbol{F}$, as a tuple of two SMatrix{3,3}.

Both are isotropic functions of C and are evaluated in closed form, from the analytic eigenvalues of sym_eigvals and the spectral projectors built from them by Lagrange interpolation, so no eigenvectors are computed and nothing is allocated. Repeated eigenvalues are handled by the reduced interpolation they call for.

This is the pair a finite-strain constitutive model needs to work in logarithmic strain space: $\boldsymbol{\varepsilon}$ and the rotated Kirchhoff stress $\hat{\boldsymbol{\tau}}$ are work conjugate, and the pull-back to the first Piola-Kirchhoff stress is $\boldsymbol{P} = \boldsymbol{F} \, \boldsymbol{U}^{-1} \hat{\boldsymbol{\tau}} \, \boldsymbol{U}^{-1}$.

Example

ε, Uinv = Peridynamics.hencky_and_invstretch(F' * F)
τ = params.λ * tr(ε) * I + 2 * params.μ * ε
P = F * (Uinv * τ * Uinv)

surface_correction_factor(correction, bond_id)

Return the surface correction factor of bond bond_id, by which a material multiplies the force density of that bond. The correction of a body chunk is obtained with get_correction and depends on the Correction type parameter of the material, see AbstractBondSystemMaterial.

NoCorrection returns 1, so a material that supports corrections can always write the multiplication unconditionally and pays nothing when no correction is used.

Example

correction = Peridynamics.get_correction(mat, i, j, bond_id)
scfactor = Peridynamics.surface_correction_factor(correction, bond_id)
b = scfactor * params.bc * ε / l .* Δxij

monomial(mat::AbstractRKCMaterial)

The symbol of the monomial basis of a material of the reproducing kernel correspondence family (:C1, :RK1, :RK2 or :PD2). It is the last type parameter of the material, so this is a compile-time constant, which makes the size of the moment matrix in rkc_weights! known to the compiler.

rkc_weights!(storage, system, mat, params, t, Δt, i)

Computes the weighted volume and the gradient weights Φ of every bond of point i from the moment matrix of its family. The monomial basis is the type parameter of mat, see monomial, so the size of the moment matrix is a compile-time constant and nothing here allocates. The regularized inverse of the moment matrix is invreg, with the adaptive singular value floor mat.epsilon if it is positive and the legacy pair mat.lambda, mat.beta otherwise, see get_invreg_params.

The neighbor of every bond enters the moment matrix with its kinematic_weight, so a damage model that weights the motion of a neighbor as data of the least-squares fit is honored here. Bonds that are not active contribute nothing.

custom_field(::Type{<:MyStorage}, ::Val{field})

Return true for every field name that a storage exports but that is not one of its own field names, e.g. a quantity that export_field derives on the fly. Without this, naming the field in Job(...; fields=(...,)) is rejected as a typo.

Defaults to false, so a field that is a field of the storage needs nothing.

Example

Peridynamics.custom_field(::Type{<:MyStorage}, ::Val{:von_mises_stress}) = true

See also export_field.

export_field(::Val{field}, mat, system, storage, paramsetup, t)

Return the point data that is written to the VTK file for field. The default returns the local points of the storage field of that name, so a field of the storage is exported without any further work.

Specialize it to derive a quantity that is not a storage field, or to reduce a bond field to a point field. A derived name also has to be announced with custom_field.

The returned array must have one column per local point, i.e. get_n_loc_points(system) of them; halo entries are owned by another chunk and must not be exported twice.

Example

Peridynamics.custom_field(::Type{<:MyStorage}, ::Val{:bond_damage_avg}) = true

function Peridynamics.export_field(::Val{:bond_damage_avg}, mat, system, storage::MyStorage,
                                   paramsetup, t)
    n = Peridynamics.get_n_loc_points(system)
    out = zeros(n)
    for i in 1:n
        bond_ids = Peridynamics.each_bond_idx(system, i)
        out[i] = sum(@view storage.bond_damage[bond_ids]) / length(bond_ids)
    end
    return out
end

sphere_shape_coords(ΔX0, radius, center)

Coordinates of points filling a disk or a ball of a given radius in concentric layers of spacing ΔX0, without the edges a cut cubic grid would leave on the surface. This is the primitive round_sphere and round_cylinder are built from.

Whether a disk or a ball is generated follows from center: a two-dimensional center gives a disk in the plane, a three-dimensional one gives a ball.

Unlike the point generators above, this returns coordinates only and no volumes, because the volume of a point depends on the shape being filled. For a solid of revolution with radius r(x), stack one disk per cross section and give each of the n points of that cross section the volume π * r(x)^2 * ΔX0 / n.

Arguments

• ΔX0::Real: Spacing of the points, and the thickness of one layer.
• radius::Real: Radius of the disk or ball.
• center::StaticVector: Center coordinates, added to all point coordinates. Its length decides the dimension: SVector{2} for a disk, SVector{3} for a ball.

Returns

• coords::Matrix{Float64}: A length(center)×n_points matrix with the point coordinates.

Examples

julia> using Peridynamics.StaticArrays

julia> coords = Peridynamics.sphere_shape_coords(1.0, 2.0, SVector{2}(0.0, 0.0));

julia> size(coords)
(2, 12)

julia> round.(coords[:, 1:4], digits=4)
2×4 Matrix{Float64}:
 0.5  -0.25   -0.25   1.5
 0.0   0.433  -0.433  0.0

@dmg_storage dmgmodel state

Declare the state a damage model carries per bond or per point, e.g. the accumulated ductile damage of a Johnson-Cook model or the number of load cycles a bond has survived. This is @cm_storage for a damage model instead of a constitutive model: the body accepts the same field declarations, so the state is allocated, sized, moved to another array backend and inherited from exactly like a storage.

The generated state is reached inside calc_failure!, calc_damage!, kinematic_weight and bond_integrity with damage_state, and a material carries it by declaring dmg_state::DamageState, see DamageState. A material that declares that field then supports every damage model without knowing any of them, because a model that needs per-bond variables brings them itself.

Differences to @storage

• Every field has to be declared with a field shape. A damage model has no init_field hook, because it does not know the material its state belongs to.
• The halo annotations @lth and @htl are not allowed. The state of a damage model is chunk-local, which is what bond state and point-local state are.
• The state is not point data, so it is not exported to VTK files. A quantity that should be exported belongs into the storage of the material, or is reduced to point data by an export_field method.

Example

struct MyDamage <: Peridynamics.AbstractDamageModel end

Peridynamics.@dmg_storage MyDamage struct MyDamageState
    bond_damage::BondScalar
end

calc_failure!(storage, system, mat, dmgmodel, paramsetup, i)

Decide which bonds of point i have failed and update the fracture bookkeeping of the storage accordingly. This is the one method a damage model has to define; it is called once per local point and per time step, right before force_density_point!.

A method sets storage.bond_active[bond_id] = false for every bond that fails and keeps storage.n_active_bonds[i] in sync, because that count is what calc_damage! turns into the damage of the point. A bond whose bond.fail_permit is false must never fail, which is how no_failure! and the pre-cracks are honored.

Arguments

• storage: The storage of the body chunk.
• system: The system of the body chunk.
• mat: The material.
• dmgmodel: The damage model, i.e. what a new model dispatches on.
• paramsetup: The parameters of the body chunk; resolve them with get_params.
• i::Int: The index of the local point that is evaluated.

Example

struct MyDamage <: Peridynamics.AbstractDamageModel end

function Peridynamics.calc_failure!(storage, system, mat, ::MyDamage, paramsetup, i)
    params = Peridynamics.get_params(paramsetup, i)
    for bond_id in Peridynamics.each_bond_idx(system, i)
        bond = system.bonds[bond_id]
        ε = ... # the stretch of the bond
        if ε > params.εc && bond.fail_permit
            storage.bond_active[bond_id] = false
        end
        storage.n_active_bonds[i] += storage.bond_active[bond_id]
    end
    return nothing
end

See also get_frac_params, has_fracture, AbstractDamageModel.

calc_damage!(storage, system, mat, dmgmodel, paramsetup, i)

Reduce the bond-wise state of the damage model to the damage of point i, which is the scalar written to storage.damage and exported as the :damage field. It is called once per local point and per time step, directly after calc_failure!.

The default is the fraction of broken bonds, 1 - n_active_bonds[i] / n_neighbors[i], which is what a model that deletes bonds wants. A model that degrades a bond continuously instead of deleting it defines its own method, so that a partially damaged bond is counted with its degree of damage rather than as intact; see kinematic_weight and bond_integrity.

Arguments

• storage: The storage of the body chunk. A stateful damage model reaches its own bond fields with damage_state.
• system: The system of the body chunk.
• mat: The material.
• dmgmodel: The damage model, i.e. what a new model dispatches on.
• paramsetup: The parameters of the body chunk; resolve them with get_params.
• i::Int: The index of the local point that is evaluated.

See also calc_failure!, @dmg_storage, AbstractDamageModel.

failure_by_sets!(storage, system, dmgmodel, set_a, set_b)

Break every bond between a point of set_a and a point of set_b, which is how a predefined crack is applied to a body chunk, and keep the fracture bookkeeping of the storage in sync. The default marks the bonds inactive; a damage model that carries its own state, see @dmg_storage, defines a method that also writes the crack into that state:

function Peridynamics.failure_by_sets!(storage, system::Peridynamics.AbstractBondSystem,
                                       dmg::MyDamage, set_a, set_b)
    Peridynamics.failure_by_sets!(storage, system, CriticalStretch(), set_a, set_b)
    # ... mark the broken bonds in `Peridynamics.damage_state(storage)`
    return nothing
end

damage_storage_type(dmgmodel)
damage_storage_type(dmgmodel, ::Type{FT})

Return the type of the state of a damage model, instantiated for the float type FT of the simulation, or Nothing for a model without state. This method is generated by @dmg_storage and is the damage-model analogue of storage_type.

get_dmg_storage(dmgmodel, solver, system)

Allocate the state of a damage model for one body chunk, or return nothing for a model without state. This method is generated by @dmg_storage and is the damage-model analogue of get_storage.

init_damage_state(mat, solver, system)

Return the value of the storage field that was declared with dmg_state::DamageState, which by default is the state of the damage model of mat. This is the escape hatch for a material family that assembles the state of its damage model itself.

damage_state(storage)

Return the state of the damage model that is carried by a storage, i.e. the field declared with dmg_state::DamageState, or nothing for a storage that does not declare one. This is how calc_failure!, calc_damage!, kinematic_weight and bond_integrity reach the per-bond variables of a stateful damage model.

has_damage_state(::Type{Storage})

Return whether a storage type declares a field with the DamageState declaration and can therefore carry the state of a damage model.

bond_integrity(dmgmodel, storage, bond_id)

Return the integrity of bond bond_id: the fraction in [0, 1] of its undamaged load-carrying capacity that the bond retains. An intact bond has integrity 1, a bond that stores and transmits nothing has integrity 0. This is the continuity 1 - d of classical damage mechanics, evaluated per bond: the free energy of a damaged bond is its undamaged free energy scaled by the integrity, and since the stress follows from the free energy, the stress the bond transmits is scaled with it.

The default is 1.0 for every damage model. A model that deletes bonds knows only intact bonds, because a failed bond is excluded through bond_active before the integrity is ever asked, so the default costs nothing. A model that softens bonds defines a method that reads its own state, see @dmg_storage and damage_state, e.g.

@inline function Peridynamics.bond_integrity(::MyDamage,
                                             storage::Peridynamics.AbstractStorage,
                                             bond_id)
    @inbounds d = Peridynamics.damage_state(storage).bond_damage[bond_id]
    return (1 - d)^2
end

The integrity is honored by the materials that evaluate their constitutive model per bond, i.e. RKCMaterial and RKCRMaterial: the stress and the strain energy density of every bond are scaled by it. One scalar per bond is the isotropic damage of classical damage mechanics; a model that degrades anisotropically, e.g. only the tensile part of the stress, instead specializes the stress computation of the material family for its damage model type. Combining a model that defines this method with a material that ignores it fails once when the Job is created, see supports_bond_integrity.

See also kinematic_weight, supports_bond_integrity, calc_failure!, @dmg_storage.

kinematic_weight(dmgmodel, storage, bond_id)

Return the weight in [0, 1] with which bond bond_id enters the reconstruction of the deformation gradient, i.e. the moment matrix and the gradient weights of RKCMaterial and RKCRMaterial.

The reconstruction is a weighted least-squares fit over the family, and every bond contributes the motion of its neighbor as data. The kinematic weight states how much of that data survives the damage of the bond: while a crack forms between two points, the neighbor turns into a point on the other side of a discontinuity, and a deformation gradient fitted through the jump produces spurious deformation and stress. Taking the weight back smoothly keeps the moment matrix a continuous function of the damage, where deleting the bond is a jump — which is what keeps fragmentation stable.

The kinematic weight is deliberately not the bond_integrity: the integrity states how much load a bond carries, the kinematic weight whether its neighbor still moves with the point. A softened bond can remain perfectly valid data. The default is therefore 1.0 for every damage model — a failed bond is excluded from the fit through bond_active, so a model that deletes bonds needs nothing else. Combining a model that defines this method with a material that ignores it fails once when the Job is created, see supports_kinematic_weight.

See also bond_integrity, supports_kinematic_weight, calc_failure!, @dmg_storage.

supports_bond_integrity(mat)

Return whether the force path of a material scales the stress and the strain energy density of every bond with its bond_integrity. Defaults to false; RKCMaterial and RKCRMaterial declare true. A custom material that applies the integrity in its own force routines declares it the same way:

Peridynamics.supports_bond_integrity(::MyMaterial) = true

The declaration is checked once when a Job is created, see check_damage_model: a damage model that defines bond_integrity combined with a material that ignores it fails there instead of silently not softening.

See also supports_kinematic_weight, bond_integrity.

supports_kinematic_weight(mat)

Return whether a material weights the bonds with their kinematic_weight when it reconstructs the deformation gradient. Defaults to false; RKCMaterial and RKCRMaterial declare true. A custom material that applies the weight in its own gradient reconstruction declares it the same way:

Peridynamics.supports_kinematic_weight(::MyMaterial) = true

The declaration is checked once when a Job is created, see check_damage_model: a damage model that defines kinematic_weight combined with a material that ignores it fails there instead of silently not softening.

See also supports_bond_integrity, kinematic_weight.

check_damage_model(spatial_setup)

Check that a damage model that defines bond_integrity or kinematic_weight is combined with a material whose force path calls the hooks, see supports_bond_integrity and supports_kinematic_weight. Throws a SofteningSupportError that names the ignored hooks and how to fix the setup.

This check is done once when a Job is created, next to check_constitutive_model.

SofteningSupportError

A type for a customized error that is thrown when a damage model defines bond_integrity or kinematic_weight, but is combined with a material whose force path never calls the hooks, so the softening would be silently ignored.

Fields

• dmgmodel::DataType: Damage model type that defines the ignored hooks.
• material::DataType: Material type that ignores them.
• reason::String: Why the combination cannot work.
• fix::String: How to fix it.

log_dmgmodel(dmgmodel; indentation)

Generate a logging message for a damage model. By default, logs only the damage model type. A damage model with properties worth logging defines its own method.

kernel(system, bond_id)

Return the value of the influence function $\omega$ of bond bond_id. The kernel is evaluated once during setup from the kernel function of the material, e.g. linear_kernel or cubic_b_spline_kernel, and the initial bond length, so a material only has to look it up.

Example

for bond_id in Peridynamics.each_bond_idx(system, i)
    ωij = Peridynamics.kernel(system, bond_id)
end

@params material struct PointParameters ... end
@params material PointParameters

Declare the point parameters of a material, i.e. the values material! assigns to a point set. The first form defines a new point parameter type and makes it the point parameters of material. The second form makes material use the point parameters that were defined for another material, e.g. CRMaterial uses the CPointParameters of CMaterial. This is the parameter counterpart of @storage: one definition per material, reuse through @inherit and the second form.

The macro generates the struct, the constructor that reads the keywords of material!, the converting constructor between float types, point_param_type, get_point_params, allowed_material_kwargs and the log_param_property methods, so that these cannot drift apart from each other.

Declaring a parameter

Every right-hand side follows one rule:

This is one rule for a constant default, a value derived from other parameters and a call that supplies a whole group, and it is why the order of the declarations matters. See @derived.

Types

Inside a definition, FT stands for the float type of the simulation, so a parameter whose type is built from it stays generic. This is how the correspondence formulation declares its elastic stiffness tensor:

@derived C::SArray{NTuple{4,3},FT,4,81} = get_hooke_matrix(nu, λ, μ)

The generated struct is parametric in FT when any parameter follows it:

struct BBPointParameters{FT<:Real,DMP} <: AbstractPointParameters
    δ::FT
    rho::FT
    ...
    bc::FT
    dmg_params::DMP
end

Peridynamics.point_param_type(mat) returns the instantiation with the float type of the simulation and Peridynamics.point_param_type(mat, Float32) the Float32 one. Point parameters are isbits, which is what lets them be captured by value in a kernel, so an array-valued parameter is an SArray, never an Array.

Reusing parameters with @inherit

@inherit includes all parameters of a parameter block defined with @params_fields or of the point parameters of another material, with the merge rules of @storage: two @inherits may contribute the same parameter only if they declare it identically, and a declaration in the body overrides an inherited one in place. This is how the dual-horizon bond-based material halves the bond constant of the bond-based one:

Peridynamics.@params DHBBMaterial struct DHBBPointParameters
    @inherit BBPointParameters
    @derived bc = 0.5 * 18 * K / (π * δ^4)
end

The package ships the blocks DiscretizationParameters, ElasticParameters, BBElasticParameters, BondHorizonParameters, InteractionParameters and StandardParameters, each documented with the parameters it exposes. The fracture parameters belong to the damage model, see @dmg_params. Names inside the definition are resolved in the module of the definition first and in Peridynamics second, so they need no import. By convention a block is named ...Parameters and a point parameter type ...PointParameters, as a field block is named ...Fields and a storage ...Storage.

Parameters of the models

A material whose constitutive model or damage model declares parameters of its own with @cm_params / @dmg_params gives them a place with a marker field, cm_params::ConstitutiveParameters or dmg_params::DamageParameters, exactly like a storage gives the state of its models a place with cm_state::ConstitutiveState. The model's declarations may read every parameter declared above the marker, and the model parameters are read flat off the point parameters, e.g. params.Gc. See ConstitutiveParameters and DamageParameters.

Examples

A material with two parameters of its own on top of the standard ones:

Peridynamics.@params MyMaterial struct MyPointParameters
    @inherit StandardParameters
    @log "yield stress" sigma_y = Inf
    hardening = 0.0
end

material!(body; horizon, rho, E, nu, Gc, sigma_y=250.0) then works, hardening is optional, anything else is rejected, and the yield stress appears in the simulation log.

A material that uses the point parameters defined for another one:

Peridynamics.@params MyOtherMaterial MyPointParameters

The generated constructor reads the material as mat, so the same type serves every material it is used for. Only point parameters defined with @params can be shared this way; a hand-written subtype of AbstractPointParameters defines point_param_type, get_point_params and allowed_material_kwargs itself.

@params_fields name begin ... end

Define a reusable block of point parameter declarations that can be included into a @params definition with @inherit. The body accepts exactly the same declarations as @params, including @inherit itself.

The macro defines a marker type name and registers the flattened declarations with param_fields_expr. A block emits no method of its own, in particular no log_param_property, because several blocks may describe the same parameter, e.g. ElasticParameters and BBElasticParameters. The @log labels travel with the declarations and are emitted by the @params definition that defines the type.

Example

Peridynamics.@params_fields DiscretizationParameters begin
    @derived (; δ, rho) = get_discretization_params(; horizon, rho)
    @log "horizon" δ
    @log "density" rho
end

The blocks this package ships are DiscretizationParameters, ElasticParameters, BBElasticParameters, BondHorizonParameters, InteractionParameters and StandardParameters, each documented with the parameters it exposes; FractureParameters belongs to the damage model and is inherited inside a @dmg_params declaration. The point parameters of a material defined with @params can be inherited the same way, e.g. @inherit BBPointParameters.

@kwarg keyword parameter

Declare the material! keyword of a parameter when it differs from the name of the parameter itself, e.g. @kwarg gamma_c gammac = 1e-10.

@derived parameter = expression
@derived (; parameter, ...) = call(args...; keyword, ...)

Declare a parameter that is computed and is not a material! keyword of its own. The second form declares a group of parameters that one call supplies together, such as the six elastic parameters, which are resolved from any two of six keywords. The call has to return a NamedTuple with one entry per destructured name.

Every right-hand side follows one rule:

@derived bc = 18 * K / (π * δ^4)
@derived (; δ, rho) = get_discretization_params(; horizon, rho)

The keywords written after ; are exactly the material! keywords the block consumes, so allowed_material_kwargs cannot disagree with the call that reads them. A keyword the user did not pass is not forwarded, which is why the provider decides on its own whether a keyword is required:

get_discretization_params(; horizon, rho)        # both required
get_elastic_params(; E=nothing, nu=nothing, …)   # any two of six

A body may declare any number of groups, in any order. Without this annotation a parameter with a default expression is a keyword that defaults to that expression.

@log "label" parameter
@log "label" declaration

Write a parameter to the simulation log under label. The first form labels a parameter that is already declared above, which is how a member of a @derived group is logged; the second declares and labels in one line. Without this annotation the parameter does not appear in the log, which is the default of log_param_property.

@log "yield stress" sigma_y = Inf
@log "bond constant" @derived bc = 18 * K / (π * δ^4)

@derived (; δ, rho) = get_discretization_params(; horizon, rho)
@log "horizon" δ

ParamFieldDecl

One parameter declaration of a @params or @params_fields body.

Fields

• name::Symbol: Name of the parameter, i.e. the field of the generated struct.
• type::Any: Declared type of the parameter: a Type, SimFloat if the parameter is declared without a type and follows the float type of the simulation, or a type expression over FT, e.g. SArray{NTuple{4,3},FT,4,81}, which follows it too.
• kwarg::Symbol: material! keyword the parameter is read from, or :none if it is not read from a keyword of its own.
• default::Any: Expression of the default value, or nothing if the keyword is required.
• provider::Any: Expression of the call that supplies the parameter, or nothing if the parameter is not a member of a @derived group.
• source::String: The right-hand side as it was written, for block_table and the error messages. It is not part of ==, because two declarations mean the same thing when they generate the same code, no matter how they were spelled.
• label::String: Label of the simulation-log line of the parameter, or "" if the parameter is not logged.

ParamFieldsSpec

The result of parsing a @params or @params_fields body: the flattened ParamFieldDecls in declaration order, and every material! keyword they consume, including the keywords a provider call reads that have no parameter of their own.

param_fields_expr(::Type{T})

Return the ParamFieldsSpec of a point parameter type or of a parameter block defined with @params_fields. This is the table @inherit reads at macro expansion time.

get_param_decls(block_args, mod)

Parse the body of a @params or @params_fields definition into a ParamFieldsSpec, expanding every @inherit and applying the merge rules documented there.

Types are resolved in mod while the macro is expanded, so an inherited declaration always keeps the meaning it has in the module that declared it. Default and provider expressions are escaped instead, because they are evaluated in the generated constructor, where they see the material mat, the keyword dictionary p and the parameters declared before them.

material_kwargs(p, names)

Return the material! keywords names that are present in p as a NamedTuple. This is what a provider call of a @derived declaration is passed: a keyword the user did not give is not forwarded, so the provider decides whether it is required.

required_param(p, kwarg)

Return the value of a required material! keyword, or throw an UndefKeywordError naming it. This is how a point parameter that is declared without a default value is read.

StandardParameters

Parameter block of the parameters a standard peridynamics model needs: the discretization and the elastic parameters, the bond constant bc, and the marker field that holds whatever parameters the damage model of the material declares — with the standard CriticalStretch, the fracture parameters Gc and εc. Materials that carry a constitutive model add cm_params::ConstitutiveParameters themselves, e.g. RKCPointParameters. See @params_fields, ConstitutiveParameters, DamageParameters.

Computed together:

• (; δ, rho) = get_discretization_params(; horizon, rho)
• (; E, nu, G, K, λ, μ) = get_elastic_params(; E, nu, G, K, lambda, mu)

Keywords of material!: horizon, rho, E, nu, G, K, lambda, mu.

DiscretizationParameters

Parameter block of the horizon δ and the density rho, the two parameters every peridynamic material needs. See @params_fields.

Computed together:

• (; δ, rho) = get_discretization_params(; horizon, rho)

Keywords of material!: horizon, rho.

ElasticParameters

Parameter block of the six elastic parameters. They are resolved together, because any two of the six keywords E, nu, G, K, lambda and mu determine all of them. See @params_fields.

Computed together:

• (; E, nu, G, K, λ, μ) = get_elastic_params(; E, nu, G, K, lambda, mu)

Keywords of material!: E, nu, G, K, lambda, mu.

FractureParameters

Parameter block of the critical energy release rate Gc and the critical stretch εc, resolved by get_frac_params of the damage model, which decides which of the fracture keywords it reads and how it converts them into each other. The block belongs to the damage model, so it is inherited inside a @dmg_params declaration: this is how CriticalStretch declares its parameters, and a custom damage model that wants the standard fracture keywords inherits it the same way. It reads the horizon δ and the bulk modulus K of the material parameters declared above the dmg_params::DamageParameters marker. See @params_fields.

Computed together:

• (; Gc, εc) = get_frac_params(model, δ, K; Gc, epsilon_c)

Keywords of material!: Gc, epsilon_c.

BondHorizonParameters

Parameter block of the bond horizon δb of a bond-associated material. It defaults to the horizon δ, so the block has to follow the one that provides it. See @params_fields.

Computed together:

• (; δb) = get_bond_horizon(δ; bond_horizon)

Keywords of material!: bond_horizon.

InteractionParameters

Parameter block of the three material constants of the continuum-kinematics-inspired formulation. They are resolved together, because they are either all derived from the elastic parameters or all specified by hand. See @params_fields.

Computed together:

• (; C1, C2, C3) = get_interaction_parameters(mat, (; δ, μ, λ); C1, C2, C3)

Keywords of material!: C1, C2, C3.

BBElasticParameters

Parameter block of the six elastic parameters of a bond-based material. Bond-based peridynamics fixes the Poisson's ratio at nu = 0.25, so a single elastic keyword is enough and any combination that results in another value is rejected. See @params_fields.

Computed together:

• (; E, nu, G, K, λ, μ) = get_elastic_params_bb(; E, nu, G, K, lambda, mu)

Keywords of material!: E, nu, G, K, lambda, mu.

@cm_params model struct ModelParameters ... end

Declare the point parameters a constitutive model owns, with the declaration language of @params_fields. The macro generates the parameter struct (parametric in the float type, isbits), its constructor, the registration that resolves a cm_params::ConstitutiveParameters marker field of a material to this type, the material! keywords the declarations consume, and the simulation-log methods of the @log labels — so a model registers its own keywords without the material knowing them.

struct LinearHardening <: Peridynamics.AbstractConstitutiveModel end

Peridynamics.@cm_params LinearHardening struct LinearHardeningParameters
    @log "yield stress" sigma_y
    @log "hardening modulus" H = 0.0
end

Inside the block the model instance is available as model (the material mat is not), and a declaration can read every material parameter declared above the marker field, e.g. the horizon δ or the bulk modulus K. A material whose point parameters carry the marker cm_params::ConstitutiveParameters supports the model; its parameters are read flat, e.g. params.sigma_y, or through the marker field, e.g. params.cm_params.sigma_y.

See also @dmg_params, ConstitutiveParameters, @params.

@dmg_params dmgmodel struct ModelParameters ... end

This is @cm_params for a damage model: the generated parameter type resolves the dmg_params::DamageParameters marker field of a material. The standard fracture parameters are declared this way — CriticalStretch owns Gc and εc through the FractureParameters block:

Peridynamics.@dmg_params CriticalStretch struct CriticalStretchParameters
    @inherit FractureParameters
end

A custom damage model that wants the standard fracture keywords inherits the same block; one that brings its own keywords declares them, and material! accepts them exactly when the body's damage model reads them.

See also DamageParameters, get_frac_params, @params.

ConstitutiveParameters

Marker for one field of a @params definition or a @params_fields block, e.g. cm_params::ConstitutiveParameters. The field holds the parameters the constitutive model of the material declares with @cm_params — whatever model that is: the parameter type is resolved per model when the point parameter type is instantiated, and it is Nothing for a model that declares no parameters. A material whose point parameters carry this marker therefore supports every constitutive model, parameterized or not, without knowing any of them.

See also DamageParameters and the storage-side analogue ConstitutiveState.

DamageParameters

Marker for one field of a @params definition or a @params_fields block, e.g. dmg_params::DamageParameters. The field holds the parameters the damage model of the material declares with @dmg_params, resolved per model exactly like ConstitutiveParameters. The standard fracture parameters Gc and εc live here: they belong to CriticalStretch, not to the material.

CriticalStretchParameters

The point parameters of CriticalStretch: the FractureParameters block, declared with @dmg_params. They occupy the dmg_params::DamageParameters marker field of every material used with the standard damage model, and are read flat off the point parameters, e.g. params.Gc.

Computed together:

• (; Gc, εc) = get_frac_params(model, δ, K; Gc, epsilon_c)

Keywords of material!: Gc, epsilon_c.

constitutive_param_type(model, FT)
damage_param_type(dmgmodel, FT)

Return the point parameter type a model declares with @cm_params / @dmg_params, instantiated with the float type FT, or Nothing for a model without parameters. This is how point_param_type resolves the marker fields of a @params definition.

constitutive_param_type(model, FT)
damage_param_type(dmgmodel, FT)

Return the point parameter type a model declares with @cm_params / @dmg_params, instantiated with the float type FT, or Nothing for a model without parameters. This is how point_param_type resolves the marker fields of a @params definition.

get_cm_params(model, FT, mat_params, p)
get_dmg_params(dmgmodel, FT, mat_params, p)

Construct the model-owned parameters of one point parameter set, or nothing for a model without parameters. Called by the generated point parameter constructor at the position of the marker field, with mat_params the NamedTuple of every material parameter declared above the marker and p the keyword dictionary of material!. The methods are generated by @cm_params / @dmg_params.

get_cm_params(model, FT, mat_params, p)
get_dmg_params(dmgmodel, FT, mat_params, p)

Construct the model-owned parameters of one point parameter set, or nothing for a model without parameters. Called by the generated point parameter constructor at the position of the marker field, with mat_params the NamedTuple of every material parameter declared above the marker and p the keyword dictionary of material!. The methods are generated by @cm_params / @dmg_params.

constitutive_param_kwargs(model)
damage_param_kwargs(dmgmodel)

The material! keywords the parameters of a model consume, () for a model without parameters. material! accepts the union of the material's own keywords and the keywords of its models, so a keyword is accepted exactly when something reads it.

constitutive_param_kwargs(model)
damage_param_kwargs(dmgmodel)

The material! keywords the parameters of a model consume, () for a model without parameters. material! accepts the union of the material's own keywords and the keywords of its models, so a keyword is accepted exactly when something reads it.

check_model_params(mat)

Check that the point parameters of a material have a place for the parameters its models declare, before material! instantiates them: a model with a @cm_params / @dmg_params declaration combined with a material whose @params definition lacks the marker field is an error naming the missing marker, and so is a parameter name that exists in more than one place.

all_material_kwargs(mat)

Every keyword material! accepts for a material: the keywords its own parameter declarations consume, plus the keywords of its constitutive model and of its damage model. A keyword is accepted exactly when something reads it, and a keyword two of them declare is an error naming both owners.

force_density_point!(storage, system, mat, params, t, Δt, i)

Calculate the internal force density of point i and accumulate it into storage.b_int. This is the one method a material has to define; everything else about a material is declared with @params and @storage.

It is called once per local point and per time step, from inside a loop over each_point_idx that runs on every thread and on every MPI rank. A method must therefore only ever write to the columns of point i, and read every other point through the system; then the same material runs under multithreading and MPI without a change.

Arguments

• storage: The storage of the body chunk, holding every field of the simulation.
• system: The system of the body chunk, e.g. a BondSystem. Dispatch on it to say which discretization the material is written for.
• mat: The material.
• params: The point parameters of point i. When a body has several parameter sets, the fallback resolves them with get_params before calling this method, so a method only ever sees one set.
• t::Real: The current simulation time.
• Δt::Real: The current time step.
• i::Int: The index of the local point that is evaluated.

Example

A bond-based material, in full:

function Peridynamics.force_density_point!(storage, system::Peridynamics.BondSystem,
                                           mat::MyMaterial, params, t, Δt, i)
    for bond_id in Peridynamics.each_bond_idx(system, i)
        bond = system.bonds[bond_id]
        j, L = bond.neighbor, bond.length

        ΔXij = Peridynamics.get_vector_diff(system.position, i, j)
        Δxij = Peridynamics.get_vector_diff(storage.position, i, j)
        l = norm(Δxij)
        ε = (l - L) / L

        # the damage model decides whether the bond is still there
        stretch_based_failure!(storage, system, bond, params, ε, i, bond_id)

        b = storage.bond_active[bond_id] * params.bc * ε / l .* Δxij
        Peridynamics.update_add_vector!(storage.b_int, i, b .* system.volume[j])
    end
    return nothing
end

See also storage_type, get_params, each_bond_idx.

storage_type(mat)
storage_type(mat, ::Type{FT})

Return the storage type of a material, instantiated for the float type FT of the simulation. The storage generated by @storage is parametric in the array type of every field, and this method returns the instantiation with the arrays of the CPU, so that body_chunk_type and everything built from it stay concrete.

init_field(material, solver, system, ::Val{field})

Return the initial value of a storage field. Every field of a storage generated by @storage is allocated by this function, in this order:

1. the system, which owns position and the fields of the fracture bookkeeping,
2. the field shape of a field declared with @lth or @htl, because whether a field has halo entries is storage knowledge that a time solver cannot have,
3. the time solver, which says whether it needs the field and at which extent,
4. the field shape of every other declared field,
5. an InterfaceError that asks for an init_field method.

A method of this function for a specific material is more specific than the fallback above and therefore wins, which is the escape hatch for a field that no field shape describes and for a shaped field that has to be filled by hand.

init_field_storage(::Type{Storage}, system, field, extent)

Return the value of field as it follows from the field shape it was declared with in the @storage definition of Storage, or nothing if the field was declared with a plain container type or does not have the requested extent. Methods of this function are generated by @storage, so a field declared with a field shape never needs an init_field method.

@storage_fields name begin ... end

Define a reusable block of storage field declarations that can be included into a @storage definition with @inherit. The body accepts exactly the same field declarations as @storage, including @inherit itself.

The macro defines a marker type name and registers the (flattened) field declarations with storage_fields_expr.

Example

@storage_fields VelocityVerletFields begin
    @lth position::PointVector{Float64}
    displacement::PointVector
    velocity::PointVector
    velocity_half::PointVector
    acceleration::PointVector
    b_int::PointVector
    b_ext::PointVector
end

@inherit storage_or_field_block

Include all field declarations of another storage or of a field block defined with @storage_fields into a @storage or @storage_fields definition. The inherited fields keep their order and are spliced in at the position of the @inherit.

Two @inherits may contribute the same field, as long as they declare it identically. A field declared in the body itself overrides an inherited field of the same name and keeps its position, which is how a family changes the annotation of an inherited field, for example from b_int::PointVector to @htl b_int::PointVector.

The inherited storage or field block has to be defined by an earlier top-level statement, and it has to be resolvable in the module of the definition, so field blocks of this package are written as Peridynamics.RKCStorage.

Example

Peridynamics.@storage MyMaterial struct MyStorage
    @inherit Peridynamics.RKCStorage
    my_point_field::PointScalar
    my_bond_field::BondScalar
end

@lth field

Annotate a field of a @storage definition for local-to-halo exchange. It has one entry per local and halo point, and the halo entries are updated from the owning chunk during the simulation.

@htl field

Annotate a field of a @storage definition for halo-to-local exchange. It has one entry per local and halo point, and the halo entries are added to the local entries of the owning chunk during the simulation.

LocalPoints

Marker for fields that have one entry per local point of a chunk. This is the extent of every field that is not annotated for halo exchange.

HaloPoints

Marker for fields that have one entry per local and halo point of a chunk. This is the extent of every field annotated with @lth or @htl, because these fields are exchanged between chunks.

SimFloat

Placeholder element type of a field shape that is written without one, e.g. PointVector instead of PointVector{Float64}. It is resolved to float_type(system) when the field is allocated, so such a field follows the floating point type of the simulation, while a shape with an explicit element type keeps it for every simulation.

AbstractPointFieldShape{T}

Supertype of all field shapes with one column per point. Every field with such a shape is point data and can therefore be exported. The number of entries is determined by the field annotation: @lth and @htl give one entry per local and halo point, an unannotated field one entry per local point.

AbstractBondFieldShape{T}

Supertype of all field shapes with one column per bond.

AbstractDofFieldShape{T}

Supertype of all field shapes with one entry per degree of freedom. These fields are neither point nor bond data, so they cannot be exported.

PointScalar
PointScalar{T}

Field shape of a scalar quantity per point, stored as a Vector{T}. Without an element type it follows the float type of the simulation, see SimFloat. See @storage for how field shapes are used.

PointVector
PointVector{T}

Field shape of a vectorial quantity per point, stored as a Matrix{T} with one column per point and get_n_dim(system) rows. See @storage.

PointTensor
PointTensor{T}

Field shape of a second order tensor per point, stored as a Matrix{T} with one column per point and get_n_dim(system)^2 rows. See @storage.

PointSymTensor
PointSymTensor{T}

Field shape of a symmetric second order tensor per point in Voigt notation, stored as a Matrix{T} with one column per point and d * (d + 1) ÷ 2 rows, where d = get_n_dim(system). See @storage.

PointField{N}
PointField{N,T}

Field shape of a quantity per point with N rows, stored as a Matrix{T} with one column per point. In contrast to every other point shape the number of rows does not follow the spatial dimension, so this is the shape for a per-point quantity that is not a vector or a tensor. See @storage.

BondScalar
BondScalar{T}

Field shape of a scalar quantity per bond, stored as a Vector{T}. See @storage.

BondVector
BondVector{T}

Field shape of a vectorial quantity per bond, stored as a Matrix{T} with one column per bond and get_n_dim(system) rows. See @storage.

BondTensor
BondTensor{T}

Field shape of a second order tensor per bond, stored as a Matrix{T} with one column per bond and get_n_dim(system)^2 rows. See @storage.

BondSymTensor
BondSymTensor{T}

Field shape of a symmetric second order tensor per bond in Voigt notation, stored as a Matrix{T} with one column per bond and d * (d + 1) ÷ 2 rows, where d = get_n_dim(system). See @storage.

BondField{N}
BondField{N,T}

Field shape of a quantity per bond with N rows, stored as a Matrix{T} with one column per bond. The number of rows does not follow the spatial dimension, see PointField.

DofVector
DofVector{T}

Field shape of a quantity per degree of freedom, stored as a Vector{T} with get_n_loc_dof(system) entries, or get_n_dof(system) entries if the field is declared with @lth or @htl. See @storage.

ConstitutiveState

Declaration of the storage field that carries the state of the constitutive model:

Peridynamics.@storage MyMaterial struct MyStorage
    @inherit VelocityVerletFields BondFracFields
    cm_state::ConstitutiveState
end

This is not a field shape: the field is not an array, it is whatever the constitutive model of the material declares with @cm_storage, and nothing for a model without state. The generated storage gets one extra type parameter CMS for it, which storage_type fills with constitutive_storage_type(get_constitutive_model(mat), FT), so the storage stays concrete whichever model is used. The state is reached with constitutive_state.

A storage may declare at most one such field, and it cannot be annotated with @lth or @htl, because the halo exchange does not descend into a nested field.

DamageState

Declaration of the storage field that carries the state of the damage model:

Peridynamics.@storage MyMaterial struct MyStorage
    @inherit VelocityVerletFields BondFracFields
    dmg_state::DamageState
end

This is the damage-model twin of ConstitutiveState: the field is not an array, it is whatever the damage model of the material declares with @dmg_storage, and nothing for a model without state. The generated storage gets one extra type parameter DMS for it, which storage_type fills with damage_storage_type(get_dmgmodel(mat), FT), so the storage stays concrete whichever damage model is used. The state is reached with damage_state.

A material that declares this field supports every damage model, stateful or not, without knowing any of them: a model that needs per-bond variables brings them itself instead of the material having to allocate them for it.

A storage may declare at most one such field, and it cannot be annotated with @lth or @htl, because the halo exchange does not descend into a nested field.

FullField()
FullField(extent)

The answer of a time solver that needs a storage field: allocate it from its field shape. Without an argument the field gets the extent it is declared with, which is HaloPoints for a field annotated with @lth or @htl and LocalPoints otherwise. With an extent the declaration is overruled, which is how NewtonKrylov gets b_int with halo entries although the field is declared without a halo annotation.

EmptyField()

The answer of a time solver that does not need a storage field: allocate it as the zero-sized array of the type it has in the storage, see alloc_empty_field.

VelocityVerletFields

The storage fields required by the VelocityVerlet time solver, see @storage_fields. A material family that accumulates into halo points overrides b_int with @htl b_int::PointVector.

DynamicRelaxationFields

The storage fields required by the DynamicRelaxation time solver, see @storage_fields. The fields it shares with VelocityVerletFields are declared identically, so a storage can inherit both blocks.

NewtonKrylovFields

The storage fields required by the NewtonKrylov time solver, see @storage_fields. The working fields of the solver have one entry per degree of freedom and are declared with DofVector; every other solver gets them as empty arrays.

BondFracFields

The storage fields of the fracture bookkeeping of a bond system, see @storage_fields. All of them are allocated by init_field_system.

InteractionFracFields

The storage fields of the fracture bookkeeping of an interaction system, see @storage_fields. All of them are allocated by init_field_system.

RKCFields

The storage fields of the reproducing kernel correspondence family, see @storage_fields. They are exactly the fields declared by req_storage_fields(::AbstractRKCMaterial).

block_table(x)

Return the markdown table of what x exposes, where x is one of

These are exactly the things @inherit accepts, so this answers "which names do I get if I inherit this".

Example

julia> print(Peridynamics.block_table(Peridynamics.DiscretizationParameters))
| parameter | type | `material!` keyword | value | simulation log |
|:---|:---|:---|:---|:---|
| `δ` | simulation float | `horizon` | required | horizon |
| `rho` | simulation float | `rho` | required | density |

field_shape(x)

Return the field shape instance a storage field is declared with, or nothing if x is a plain container type such as Matrix{Float64}. This is how @storage tells the two kinds of field declaration apart without having to know the shapes by name, so a new shape only needs methods of field_shape, container_type and field_n_rows to be usable.

field_n_dims(shape)

Return the number of dimensions of the array that holds a field of the given field_shape.

container_type(shape)
container_type(shape, eltype_expr)

Return the array type that holds a field of the given field_shape. With one argument this is a concrete type, where a shape that does not pin its element type uses default_float_type. With two arguments it is the type expression used in the generated storage struct, where eltype_expr is either a type or the symbol of a type parameter of the storage.

field_n_rows(shape, system)

Return the number of rows of a field with the given field_shape, or nothing if the field is stored as a Vector. Everything that depends on the spatial dimension is resolved here through get_n_dim(system), so this is the single place that has to change when a system with a different dimension is added.

field_n_entries(shape, system, extent)

Return the number of entries (columns) of a field with the given field_shape. For point shapes this is the number of local points (LocalPoints) or the number of local and halo points (HaloPoints), depending on the field annotation. For bond shapes it is the number of bonds.

field_eltype(shape)
field_eltype(shape, system)

Return the element type of a field with the given field_shape. With one argument this is SimFloat for a shape that does not pin its element type, with a system it is the resolved type, i.e. float_type(system).

concrete_eltype(shape)

Return the element type a field with the given field_shape has in a simulation with the default_float_type.

alloc_field(shape, system, extent)
alloc_field(shape, system, extent, init)

Allocate the array of a storage field with the given field_shape. The element type is field_eltype(shape, system), the size follows from the shape, the system and the extent. The initial value init is

A UniformScaling requires a shape whose rows are the entries of a square tensor.

alloc_empty_field(::Type{Storage}, field)

Allocate the zero-sized array of the type that field has in Storage. This is what a time solver gets for a field it does not need, see EmptyField. The type comes from the storage and not from the field shape, so this works for a field declared with a container type as well, and gives arrays of the right backend for a storage that lives on a device.

alloc_solver_field(marker, shape, system, extent, init)

Allocate the array of a storage field that a time solver answered for with the FullField marker, where shape, extent and init are the field shape, the extent and the initial value of the declaration. Methods of this function are called by the init_field_storage methods that @storage generates.

storage_field_type(x)

Return the type of a storage field declared as x. A field_shape is resolved to its container_type, everything else is a plain container type and is returned unchanged. @storage wraps every declared field type with this function, so both kinds of declaration can be mixed freely inside one storage.

storage_fields_expr(::Type{T})

Return the field declarations of a storage type or of a field block defined with @storage_fields as a Vector of StorageFieldDecl. This is the table @inherit reads at macro expansion time.

StorageFieldDecl

One field declaration of a @storage or @storage_fields block: the field name, its annotation (:none, :lth or :htl), the declared type expression (a container type or a field_shape) and the init expression of the declaration, or nothing if it does not specify one.

get_field_decls(block_args, mod)

Parse the body of a @storage or @storage_fields definition into a Vector of StorageFieldDecl, expanding every @inherit and applying the merge rules documented there.

Field types and initial values are resolved in mod while the macro is expanded, so an inherited declaration always keeps the meaning it has in the module that declared it.

StorageTypeParam

One derived type parameter of a storage generated by @storage: its name, the element type of the fields it stands for, which is SimFloat if they follow the float type of the simulation, and the number of dimensions n_dims of their arrays.

derive_storage_type_params(decls)

Derive the type parameters of a storage from its field declarations. Returns the derived StorageTypeParams in the order they appear in the struct header, a Dict that maps the name of a field to the name of its type parameter, whether any field follows the float type of the simulation and the header therefore starts with FT, and the names of the fields declared with ConstitutiveState and DamageState, which contribute the trailing parameters CMS and DMS in that order.

storage_field_decl_type(decl, field_params)

Return the type of a field in the generated storage struct, which is the name of its derived type parameter, or the type it is declared with if it has none.

storage_param_bound(param)

Return the type expression a derived type parameter is bounded by in the generated storage struct, e.g. AbstractMatrix{FT} for M_FT. The bound is abstract, so that a storage can hold the arrays of any backend.

storage_param_default(param)

Return the type expression of the array a derived type parameter stands for in a simulation on the CPU, e.g. Matrix{FT} for M_FT. This is what storage_type instantiates the storage with.

get_storage_header(storage_expr)

Return the name and the supertype of the storage struct of a @storage definition. A definition that declares its own type parameters is rejected, because @storage derives them from the field declarations.

storage_header_expr(storage_type, supertype, params, uses_sim_float)

Assemble the header of the generated storage struct from the derived type parameters, see derive_storage_type_params. A storage without a derived parameter, i.e. one whose fields are all declared with a type that is not an Array, keeps a plain header.

required_fields(::Type{Material})

Return all storage fields that can be determined from the material type alone. This is the cheap check done by @storage at macro expansion time.

Everything that depends on the material instance, the damage model or the time solver cannot be known here and is therefore checked at setup time with check_storage_contract.

required_fields_timesolver(::Type{TS})

Return all storage fields required by the time solver type TS.

req_storage_fields(mat)
req_storage_fields(mat, dmgmodel)
req_storage_fields(solver)

Return the storage fields that are required by a material family, by a damage model in combination with a material, or by a time solver. All methods return a Tuple of Symbols and default to an empty tuple, so only the parts of the contract that really exist have to be declared.

These declarations are collected by check_storage_contract when a Job is created.

Example

Peridynamics.req_storage_fields(::MyMaterialFamily) = (:defgrad, :weighted_volume)
Peridynamics.req_storage_fields(::MyMaterialFamily, ::MyDamageModel) = (:bond_damage,)

storage_contract(mat::AbstractMaterial, solver::AbstractTimeSolver)

Return the complete storage contract of mat in combination with solver as a Tuple of (fields, reason) pairs, where reason describes why the fields are required.

check_storage_contract(spatial_setup, solver::AbstractTimeSolver)

Check if the storage of every material of spatial_setup contains all fields that are required by the material itself, by its damage model and by the time solver solver. Throws a StorageContractError that names the missing fields and the reason why they are required.

This check is done once when a Job is created. In contrast to the check inside @storage, the damage model and the time solver are known here, so the parts of the contract that depend on them can finally be verified.

StorageContractError

A type for a customized error that is thrown when the storage of a material does not contain all fields that are required by the material, its damage model and the time solver of the submitted job. See check_storage_contract.

Fields

• storage::DataType: Storage type that misses fields.
• material::DataType: Material type the storage belongs to.
• solver::DataType: Time solver of the job.
• missing_fields::Vector{Pair{Symbol,String}}: Missing fields and the reason why they are required.

default_float_type()

Return the floating point type of all simulation data. This is the single place that decides it, so nothing below hardcodes Float64 any more. It will become a property of Body and flow into the system, the point parameters and the storage; until then every float_type method returns this value.

float_type(x)

Return the floating point type of the simulation data of x. Every storage field declared with a field shape that does not pin its element type, e.g. PointVector instead of PointVector{Float64}, is allocated with this type, see @storage.

get_n_loc_points(system)

Return the number of local points of a chunk, i.e. the points this chunk is responsible for and whose equation of motion it integrates. Always ≤ get_n_points(system).

A field declared with the default extent LocalPoints has this many entries.

See also get_n_points, each_point_idx.

get_n_points(system)

Return the number of points of a chunk, local and halo. Halo points are the points of other chunks that this chunk needs in order to evaluate the force density of its own points; they are read from the neighboring chunks and never integrated here.

A field annotated with @lth or @htl has this many entries, i.e. the extent HaloPoints.

See also get_n_loc_points.

each_point_idx(system)

Return an iterator over the indices of the local points of a chunk. This is the loop the force density calculation of a body chunk runs over, and it is what makes a material work under multithreading and MPI without a change: every chunk only ever iterates its own points.

Example

for i in Peridynamics.each_point_idx(system)
    force_density_point!(storage, system, mat, params, t, Δt, i)
end

See also each_bond_idx, get_n_loc_points.

each_bond_idx(system, i)

Return an iterator over the bond indices of point i. A bond index addresses every bond field of the storage, i.e. every field declared with a Bond... field shape, and system.bonds[bond_id] gives the bond itself with its neighbor and its length.

Example

for bond_id in Peridynamics.each_bond_idx(system, i)
    bond = system.bonds[bond_id]
    j, L = bond.neighbor, bond.length
    storage.bond_active[bond_id] || continue
end

See also each_point_idx, get_n_bonds.

get_n_bonds(system)

Return the number of bonds of a chunk. Every storage field declared with a Bond... field shape has this many entries, and so does the state of a constitutive model that indexes per bond, see @cm_storage.

get_dmgmodel(mat::AbstractMaterial)

Return the damage model of mat, or nothing if the material has no damage model.

@cm_storage model state

Declare the state a history-dependent constitutive model integrates over time, e.g. the plastic strain of a plasticity model. This is @storage for a constitutive model instead of a material: the body accepts the same field declarations, so the state is allocated, sized, moved to another array backend and inherited from exactly like a storage.

The generated state is reached inside first_piola_kirchhoff with constitutive_state, and a material carries it by declaring cm_state::ConstitutiveState, see ConstitutiveState.

Differences to @storage

• Every field has to be declared with a field shape. A constitutive model has no init_field hook, because it does not know the material its state belongs to.
• The halo annotations @lth and @htl are not allowed. The state of a constitutive model is chunk-local, which is what bond state and point-local state are.
• The state is not point data, so it is not exported to VTK files. A quantity that should be exported belongs into the storage of the material.

Example

struct J2Plasticity <: Peridynamics.AbstractConstitutiveModel end

Peridynamics.@cm_storage J2Plasticity struct J2PlasticityState
    bond_plastic_strain::BondTensor
    bond_eqps::BondScalar
end

J2Plasticity is history-dependent by this declaration alone, see is_history_dependent.

get_constitutive_model(mat)

Return the constitutive model of a material, or nothing for a material family that does not have one, e.g. the bond-based materials. Every family that stores a constitutive model defines this method, which is how the setup checks and the storage of the constitutive state reach the model without knowing the material.

Example

Peridynamics.get_constitutive_model(mat::MyMaterial) = mat.constitutive_model

first_piola_kirchhoff(model, storage, params, F, idx, Δt)

Return the first Piola-Kirchhoff stress $\boldsymbol{P}$ of a constitutive model for the deformation gradient $\boldsymbol{F}$. This is the one method a constitutive model has to define.

Arguments

• model: The constitutive model.

• storage: The storage of the body chunk. The state of a history-dependent model is reached with constitutive_state.

• params: The point parameters of the point that is evaluated.

• F::SMatrix{3,3}: The deformation gradient.

• idx::Int: The index of the evaluated quantity, which is what a history-dependent model indexes its state with. Which index this is follows from the material family:

  | material family | idx | |:–-|:–-| | CMaterial | point index | | RKCMaterial, BACMaterial | bond index |

• Δt::Real: The time step, i.e. the increment a history-dependent model integrates over.

A model without state

A model that only needs the deformation gradient defines the four-argument form and ignores this one, which is bridged to it:

function Peridynamics.first_piola_kirchhoff(::MyModel, storage, params, F)
    return ...
end

A model with state

struct J2Plasticity <: Peridynamics.AbstractConstitutiveModel end

Peridynamics.@cm_storage J2Plasticity struct J2PlasticityState
    bond_plastic_strain::BondTensor
    bond_eqps::BondScalar
end

function Peridynamics.first_piola_kirchhoff(cm::J2Plasticity, storage, params, F, idx, Δt)
    state = Peridynamics.constitutive_state(storage)
    εᵖ = Peridynamics.get_tensor(state.bond_plastic_strain, idx)
    # ... radial return ...
    Peridynamics.update_tensor!(state.bond_plastic_strain, idx, εᵖ_new)
    return P
end

strain_energy_density(model, storage, params, F, idx)

Return the strain energy density $\Psi$ of a constitutive model. In contrast to first_piola_kirchhoff this method must not change the state of the model: it is also called when a field is exported, i.e. outside of the time integration. See first_piola_kirchhoff for the meaning of idx.

A model that only needs the deformation gradient defines the four-argument form strain_energy_density(model, storage, params, F), which this one is bridged to.

constitutive_storage_type(model)
constitutive_storage_type(model, ::Type{FT})

Return the type of the state of a constitutive model, instantiated for the float type FT of the simulation, or Nothing for a model without state. This method is generated by @cm_storage and is the constitutive-model analogue of storage_type.

get_cm_storage(model, solver, system)

Allocate the state of a constitutive model for one body chunk, or return nothing for a model without state. This method is generated by @cm_storage and is the constitutive-model analogue of get_storage.

init_constitutive_state(mat, solver, system)

Return the value of the storage field that was declared with cm_state::ConstitutiveState, which by default is the state of the constitutive model of mat. This is the escape hatch for a material family that assembles the state of its constitutive model itself.

constitutive_state(storage)

Return the state of the constitutive model that is carried by a storage, i.e. the field declared with cm_state::ConstitutiveState, or nothing for a storage that does not declare one. This is how first_piola_kirchhoff reaches the state of a history-dependent model.

has_constitutive_state(::Type{Storage})

Return whether a storage type declares a field with the ConstitutiveState declaration and can therefore carry the state of a history-dependent constitutive model.

is_history_dependent(model)

Return whether a constitutive model integrates an internal state over time, e.g. the plastic strain of a plasticity model. Defaults to whether the model declares a state with @cm_storage, so a model that keeps its history in the storage of the material has to declare it:

Peridynamics.is_history_dependent(::MyPlasticModel) = true

A history-dependent model may only be used with a time solver that evaluates the force density once per time step, and only with a material whose storage carries the state. Both are checked once when a Job is created.

supports_history_dependence(solver)
supports_history_dependence(mat)

Return whether a time solver or a material may be combined with a history-dependent constitutive model. Defaults to true for both.

A time solver answers false if it evaluates the force density more than once per time step, because the model would then integrate its history several times per step. NewtonKrylov does this for the Jacobian-vector products and the line search.

A material answers false if the index that first_piola_kirchhoff receives does not always index the same quantity, because the state of the model could then not be indexed at all.

check_constitutive_model(spatial_setup, solver)

Check that a history-dependent constitutive model is combined with a time solver that evaluates it once per time step, with a material whose evaluation index is well defined, and with a storage that carries the state of the model. Throws a HistoryDependenceError that names the reason and how to fix it.

This check is done once when a Job is created, next to check_storage_contract.

HistoryDependenceError

A type for a customized error that is thrown when a history-dependent constitutive model is combined with a time solver, a material or a storage that cannot integrate its history correctly. See check_constitutive_model.

Fields

• model::DataType: Constitutive model that integrates a history.
• material::DataType: Material type the model belongs to.
• solver::DataType: Time solver of the job.
• reason::String: Why the combination cannot work.
• fix::String: How to fix it.

@storage material storage
@storage material solver storage

A macro for automatic creation of a storage type. The macro then maps a specified material type with the storage structure. When no solver is specified, this storage is always used independently of the solver type. Specifying a solver allows the usage of customized storages for a solver.

Declaring a field

A field is declared with a field shape, which says what the field is. The shape determines the container type, the element type and the allocation, so a field declared with a shape needs no init_field method:

Everything that depends on the spatial dimension goes through get_n_dim(system), so a shaped field needs no change when a system of a different dimension is added.

A point shape makes the field point data, which is what can be exported to VTK files. A bond or dof shape is not point data.

Element type

A shape written without an element type follows the float type of the simulation, see SimFloat. A shape with an element type pins it for every simulation, which is why position is declared PointVector{Float64}: bond vectors are position differences and a smaller float type loses them over a large domain. PointScalar{Bool} and PointScalar{Int} are the non-float cases.

Annotations

Only the halo exchange is annotated, everything else follows from the shape. Both annotations give the field one entry per local and halo point:

• @lth: local-to-halo exchange during the simulation, i.e. the halo entries are updated from the chunk that owns the points.
• @htl: halo-to-local exchange, i.e. the halo entries are added to the local entries of the owning chunk.
• @inherit: includes all fields of one or more other storages or of field blocks defined with @storage_fields. See @inherit for the merge rules.

Initial value

The initial value of a shaped field is written = value:

A UniformScaling requires a shape whose rows are the entries of a square tensor, e.g. PointTensor.

Anything a shape does not describe

A field can also be declared with a plain container type, e.g. zem_stiffness::MArray{NTuple{4,3},Float64,4,81}. Such a field is not point data and it needs an init_field method that allocates it. A container type that is an Array, e.g. Vector{Float64}, still takes part in the type parameters described below, everything else keeps the type it is declared with.

A field that a time solver works with is the one case where a container type is not enough. A solver says only whether it needs the field and at which extent and leaves the number of rows and the element type to the shape, so it cannot size a container-typed field. Inheriting the field block of the solver, e.g. @inherit VelocityVerletFields, declares all of them with the right shape.

An init_field method is also the escape hatch for a shaped field: it is more specific than the generic fallback and therefore wins, so a field can keep its shape — and with it its type, its size and its export status — while being filled by hand:

function Peridynamics.init_field(::MyMaterial, ::Peridynamics.AbstractTimeSolver,
                                 system::Peridynamics.AbstractSystem, ::Val{:my_field})
    return ...
end

The generated type

The generated struct must not declare type parameters of its own, because the macro derives them from the declarations: one parameter per distinct combination of element type and number of dimensions, plus a leading FT for the float type of the simulation. The storage of the bond-based material below becomes

struct BBStorage{FT<:Real,M_F64<:AbstractMatrix{Float64},M_FT<:AbstractMatrix{FT},
                 V_FT<:AbstractVector{FT},V_Int<:AbstractVector{Int},
                 V_Bool<:AbstractVector{Bool}} <: AbstractStorage
    position::M_F64
    displacement::M_FT
    ...
end

so storage_type(mat) returns the instantiation with the arrays of the CPU, storage_type(mat, Float32) the one with Float32 arrays, and Adapt.adapt(backend, s) moves a whole storage to another array backend. Dispatch on the storage therefore has to be written ::Type{<:BBStorage} instead of ::Type{BBStorage}, while ::BBStorage is unchanged.

Example

Example definition of the storage for the bond-based material:

@storage BBMaterial struct BBStorage <: AbstractStorage
    @inherit VelocityVerletFields DynamicRelaxationFields NewtonKrylovFields BondFracFields
    strain_energy_density::PointScalar
    bond_length::BondScalar
end

@autoinfiltrate
@autoinfiltrate condition::Bool

Invoke the @infiltrate macro of the package Infiltrator.jl to create a breakpoint for ad-hoc interactive debugging in the REPL. If the optional argument condition is given, the breakpoint is only enabled if condition evaluates to true.

As opposed to using Infiltrator.@infiltrate directly, this macro does not require Infiltrator.jl to be added as a dependency to Peridynamics.jl. As a bonus, the macro will also attempt to load the Infiltrator module if it has not yet been loaded manually.

Note: For this macro to work, the Infiltrator.jl package needs to be installed in your current Julia environment stack.

See also: Infiltrator.jl

velocity_databc!(body, data, set_name, dims)

Specifies velocity boundary conditions for points of the set set_name in body. The value of the boundary condition is assigned by reading the corresponding positions in the matrix data. Multiple dimensions can be handled at once.

Arguments

• body::AbstractBody: Body the condition is specified on.
• data::Matrix: A matrix of size length(dims) x n_points that contains the values of the boundary condition for each point in the body. But only the conditions of points contained in the set set_name are applied during the simulation! It should be noted, that the value of the data matrix is constant and currently cannot be updated during the simulation. The data matrix is not checked for NaN values, since this is handled in the apply_bc! function. If it contains NaN values, then these values are ignored.
• set_name::Symbol: The name of a point set of this body. The condition applies only to the points in this set, even if the data matrix contains values for all points in the body.
• dims::Vector{Union{Integer,Symbol}}: Vector containing the directions of the condition that should be applied, either specified as Symbol or integer.
  • x-direction: :x or 1
  • y-direction: :y or 2
  • z-direction: :z or 3
  It should not contain more than 3 elements, and the elements should be unique. The order of the elements does not matter, however it must match the values in the data matrix. So if the first column of the data matrix contains the values for the x-direction, then the first element of dims should be 1 or :x, and so on.

Throws

• Errors if the body does not contain a set with set_name.
• Errors if the directions are not correctly specified.
• Errors if the dimensions of the data matrix are incorrect.

forcedensity_databc!(body, data, set_name, dims)

Specifies forcedensity boundary conditions for points of the set set_name in body. The value of the boundary condition is assigned by reading the corresponding positions in the matrix data. Multiple dimensions can be handled at once.

Arguments

• body::AbstractBody: Body the condition is specified on.
• data::Matrix: A matrix of size length(dims) x n_points that contains the values of the boundary condition for each point in the body. But only the conditions of points contained in the set set_name are applied during the simulation! It should be noted, that the value of the data matrix is constant and currently cannot be updated during the simulation. The data matrix is not checked for NaN values, since this is handled in the apply_bc! function. If it contains NaN values, then these values are ignored.
• set_name::Symbol: The name of a point set of this body. The condition applies only to the points in this set, even if the data matrix contains values for all points in the body.
• dims::Vector{Union{Integer,Symbol}}: Vector containing the directions of the condition that should be applied, either specified as Symbol or integer.
  • x-direction: :x or 1
  • y-direction: :y or 2
  • z-direction: :z or 3
  It should not contain more than 3 elements, and the elements should be unique. The order of the elements does not matter, however it must match the values in the data matrix. So if the first column of the data matrix contains the values for the x-direction, then the first element of dims should be 1 or :x, and so on.

Throws

• Errors if the body does not contain a set with set_name.
• Errors if the directions are not correctly specified.
• Errors if the dimensions of the data matrix are incorrect.