phasefield.jl

# ===========================================================================
#  PhaseFieldFracture — a readable, pure-Julia phase-field fracture solver
#
#  Author:     Yang Bai — Materials Mechanics Laboratory (MMLab)
#  Contact:    yangbai90@outlook.com
#  Copyright:  © 2026 Materials Mechanics Laboratory (MMLab). All rights reserved.
# ===========================================================================

# ===========================================================================
#  AT2 phase-field damage
# ===========================================================================
#
#  Damage φ ∈ [0,1]:  φ = 0 intact,  φ = 1 fully broken.  The regularized
#  fracture energy (AT2 crack-surface density) is
#
#      Gc ∫_Ω ( φ²/(2ℓ) + (ℓ/2)|∇φ|² ) dΩ .
#
#  Minimizing the total energy g(φ)ψ⁺ + fracture energy over φ, with the
#  irreversible history field H replacing the tensile elastic energy ψ⁺, gives
#  the linear strong form
#
#      (2H + Gc/ℓ) φ − Gc ℓ Δφ = 2H .
#
#  (Only the TENSILE energy ψ⁺ drives the crack — see strainsplit.jl.  The
#  phase-field equation itself is unchanged; only the definition of H changes,
#  from max ψ to max ψ⁺, when the tension/compression split is switched on.)
#
#  Its weak form yields, on every cell,
#
#      Kₑ[a,b] = ∫ ( (2H + Gc/ℓ) Nₐ N_b + Gc ℓ ∇Nₐ·∇N_b ) dΩ
#      fₑ[a]   = ∫ ( 2H Nₐ ) dΩ .
#
#  Irreversibility (no crack healing) is imposed in two ways:
#    • the history field H can only grow in time, and
#    • the nodal solution is projected so that φ ≥ φ from the previous step.

"""
Assemble the global AT2 phase-field system  K φ = f  for a fixed history field
`H` (values per Gauss point per cell, size `length(GAUSS) × ncells`).
"""
function phasefield_system(mesh::Mesh, mat::Material, H)
    n     = nnodes(mesh)                 # one damage dof per node
    ncell = ncells(mesh)
    # COO triplets again; each cell contributes a 4×4 block = 16 stiffness entries.
    I = Vector{Int}(undef, 16 * ncell)   # row indices
    J = Vector{Int}(undef, 16 * ncell)   # column indices
    V = Vector{Float64}(undef, 16 * ncell)  # values
    f = zeros(n)                         # global right-hand side (source term)
    p = 0                                # running write position into I/J/V

    for e in 1:ncell
        nodes = @view mesh.cells[e, :]   # the 4 node indices of this cell
        X     = mesh.coords[nodes, :]    # their coordinates, a 4×2 matrix

        Kₑ = zeros(4, 4)                 # element phase-field matrix
        fₑ = zeros(4)                    # element source vector
        for (q, (ξ, η)) in enumerate(GAUSS)      # q indexes the Gauss point (for H)
            N, ∇N, detJ = element(X, ξ, η)
            H_q = H[q, e]                # history (driving force) at this Gauss point
            # N*N' is the 4×4 mass block [a,b]=NₐN_b; ∇N*∇N' is the 4×4 diffusion
            # block [a,b]=∇Nₐ·∇N_b.  Together they give Kₑ; fₑ is the 2H source.
            Kₑ .+= ((2H_q + mat.Gc / mat.ℓ) .* (N * N') .+
                    mat.Gc * mat.ℓ .* (∇N * ∇N')) .* detJ
            fₑ .+= (2H_q .* N) .* detJ
        end

        for a in 1:4                     # scatter the 4×4 block and 4-vector globally
            f[nodes[a]] += fₑ[a]
            for b in 1:4
                p += 1
                I[p] = nodes[a]; J[p] = nodes[b]; V[p] = Kₑ[a, b]
            end
        end
    end
    return sparse(I, J, V, n, n), f
end

"""
Update the trial history field from a displacement field `u`:

    H_trial = max( H_committed , ψ⁺(ε(u)) )   at every Gauss point,

where ψ⁺ is the crack driving energy `driving_energy(mat, ε)` — the TENSILE
part of the strain energy for the spectral split, or the total energy for the
`:none` model.  Taking the running maximum makes the driving force monotone in
time, which is what prevents the damage from reversing (no crack healing).
"""
function update_history!(H_trial, H_committed, mesh::Mesh, mat::Material, u)
    for e in 1:ncells(mesh)
        nodes = @view mesh.cells[e, :]
        X     = mesh.coords[nodes, :]
        uₑ    = u[cell_dofs(nodes)]
        for (q, (ξ, η)) in enumerate(GAUSS)
            _, ∇N, _ = element(X, ξ, η)
            ε = strain_matrix(∇N) * uₑ
            H_trial[q, e] = max(H_committed[q, e], driving_energy(mat, ε))
        end
    end
    return H_trial
end

"""
Solve the phase-field problem for the given history `H`, then enforce
irreversibility:  φ ≥ `φ_prev` everywhere, φ ∈ [0,1], and φ = 1 on the
pre-crack `crack_nodes`.
"""
function solve_phase(mesh::Mesh, mat::Material, H, φ_prev, crack_nodes)
    K, f = phasefield_system(mesh, mat, H)
    bc   = Dict(n => 1.0 for n in crack_nodes)      # φ = 1 on the notch
    φ    = solve_dirichlet(K, f, bc)
    @. φ = clamp(max(φ, φ_prev), 0.0, 1.0)          # no healing, keep in [0,1]
    for n in crack_nodes
        φ[n] = 1.0
    end
    return φ
end