strainsplit.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.
# ===========================================================================

# ===========================================================================
#  Tension–compression split of the elastic energy (Miehe SPECTRAL split)
# ===========================================================================
#
#  WHY DO WE NEED THIS?  (read this first — it is the whole point of the file)
#  ------------------------------------------------------------------------
#  In the simplest phase-field model the *entire* strain energy drives the
#  crack and the *entire* stiffness is degraded by g(φ).  That has two
#  unphysical consequences:
#
#    1. a crack can nucleate and grow under pure COMPRESSION, and
#    2. the two faces of an open crack can freely interpenetrate, because the
#       material has lost all stiffness even against closing (compressive) load.
#
#  Real brittle cracks only open and grow under TENSION.  To capture this we
#  split the elastic energy into a "tensile" part ψ⁺ and a "compressive" part
#  ψ⁻, and let the damage degrade (and be driven by) ONLY the tensile part:
#
#        ψ(ε, φ)  =  g(φ) ψ⁺(ε)  +  ψ⁻(ε) ,        g(φ) = (1-φ)² + k .
#
#  A fully broken point (φ = 1, g ≈ 0) then loses its tensile stiffness but
#  keeps its full compressive stiffness — crack faces can reopen but cannot
#  pass through each other, and compression never drives damage.
#
#  This file implements the SPECTRAL (principal-strain) split of
#  C. Miehe, F. Welschinger, M. Hofacker, IJNME 83 (2010) 1273–1311, and
#  C. Miehe, M. Hofacker, F. Welschinger, CMAME 199 (2010) 2765–2778.
#
#  ------------------------------------------------------------------------
#  THE MATHEMATICS
#  ------------------------------------------------------------------------
#  Write the isotropic elastic energy with the Lamé constants λ, μ:
#
#        ψ(ε) = ½ λ (tr ε)²  +  μ (ε : ε) ,        σ = λ (tr ε) I + 2μ ε .
#
#  Diagonalise the (symmetric) strain tensor into principal strains εₐ and
#  principal directions nₐ  (its "spectral decomposition"):
#
#        ε = Σₐ εₐ (nₐ ⊗ nₐ) .
#
#  Define the Macaulay brackets (positive/negative parts)
#
#        ⟨x⟩₊ = max(x, 0) ,      ⟨x⟩₋ = min(x, 0) ,      ⟨x⟩₊ + ⟨x⟩₋ = x ,
#
#  and split the strain by keeping only positive / only negative principal
#  strains:
#
#        ε⁺ = Σₐ ⟨εₐ⟩₊ (nₐ ⊗ nₐ) ,     ε⁻ = Σₐ ⟨εₐ⟩₋ (nₐ ⊗ nₐ) .
#
#  Miehe's energy split (this is the exact definition we implement):
#
#        ψ⁺(ε) = ½ λ ⟨tr ε⟩₊²  +  μ (ε⁺ : ε⁺)  =  ½ λ ⟨tr ε⟩₊² + μ Σₐ ⟨εₐ⟩₊² ,
#        ψ⁻(ε) = ½ λ ⟨tr ε⟩₋²  +  μ (ε⁻ : ε⁻)  =  ½ λ ⟨tr ε⟩₋² + μ Σₐ ⟨εₐ⟩₋² .
#
#  Because ⟨x⟩₊² + ⟨x⟩₋² = x², these sum back to the total energy: ψ⁺+ψ⁻ = ψ.
#
#  Differentiating the energy gives the split stress σ = ∂ψ/∂ε:
#
#        σ⁺ = ∂ψ⁺/∂ε = λ ⟨tr ε⟩₊ I + 2μ ε⁺ ,
#        σ⁻ = ∂ψ⁻/∂ε = λ ⟨tr ε⟩₋ I + 2μ ε⁻ ,        σ⁺ + σ⁻ = σ  (full stress).
#
#  Differentiating once more gives the consistent (algorithmic) TANGENT
#  C = ∂σ/∂ε, needed because the split makes the stress a NONLINEAR function
#  of the strain (so the mechanics is solved by Newton iteration, see
#  `solve_displacement` in elasticity.jl):
#
#        C⁺ = ∂σ⁺/∂ε = λ H(tr ε) (I ⊗ I) + 2μ ℙ⁺ ,
#        C⁻ = ∂σ⁻/∂ε = λ H(-tr ε) (I ⊗ I) + 2μ (𝕀 - ℙ⁺) ,
#
#  where H is the Heaviside step (H = d⟨·⟩₊/dx) and ℙ⁺ = ∂ε⁺/∂ε is the
#  fourth-order "positive projection" tensor of the strain (Miehe 1998).
#
#  ------------------------------------------------------------------------
#  PLANE STRAIN.  Out of plane ε_zz = 0, so 0 is automatically a third
#  principal strain.  Since ⟨0⟩₊ = ⟨0⟩₋ = 0 it contributes nothing to ε⁺, ε⁻
#  or to tr ε.  Hence we may work entirely with the in-plane 2×2 strain
#  tensor and its two principal strains — no special out-of-plane term.
#
#  ------------------------------------------------------------------------
#  "MANDEL" NOTATION — the one bookkeeping trick used below.
#  A symmetric 2×2 tensor has 3 independent numbers.  We store it as the
#  3-vector
#
#        â = [a₁₁ ,  a₂₂ ,  √2 · a₁₂]                     (Mandel form)
#
#  The √2 is chosen so that the tensor double-contraction becomes the ORDINARY
#  dot product,  a : b = â · b̂,  and consequently every fourth-order tensor
#  becomes an ordinary SYMMETRIC 3×3 matrix acting on these 3-vectors.  This
#  makes the projection ℙ⁺ and the tangents C⁺, C⁻ plain 3×3 matrices that we
#  can build and read directly.  At the very end we convert Mandel → the
#  "engineering" convention [εₓₓ, ε_yy, γₓy] used by the B-matrix and the
#  assembler, via the diagonal map T = diag(1, 1, 1/√2):
#
#        σ_eng = T σ̂ ,        C_eng = T Ĉ T .
#
#  (Engineering strain uses γₓy = 2ε₁₂, and stress keeps σₓy once; T carries
#  exactly those factors of 2.  As a check, for the undamaged solid this
#  reproduces the standard plane-strain matrix D — see the unit tests.)

# --- small constants used by the Mandel bookkeeping -------------------------
const SQRT2 = sqrt(2.0)
const IVEC  = [1.0, 1.0, 0.0]         # the 2×2 identity tensor, in Mandel form
const T_M2E = [1.0, 1.0, 1.0/SQRT2]   # diagonal of the Mandel→engineering map T
const IMAT  = Matrix(1.0I, 3, 3)      # 3×3 identity = the 4th-order identity 𝕀 (Mandel)

# --- Macaulay brackets (positive / negative parts) and the Heaviside step ---
"Positive part (Macaulay bracket)  ⟨x⟩₊ = max(x, 0)."
pos(x) = max(x, 0.0)
"Negative part  ⟨x⟩₋ = min(x, 0);  note ⟨x⟩₊ + ⟨x⟩₋ = x."
neg(x) = min(x, 0.0)
"Heaviside step  H(x) = d⟨x⟩₊/dx  (=1 for x>0, else 0); the slope of the ramp ⟨x⟩₊."
heaviside(x) = x > 0.0 ? 1.0 : 0.0

"""
    principal_strains(ε) -> (ε1, ε2)

Principal (eigen)values of the in-plane 2×2 strain tensor, given the strain in
ENGINEERING form `ε = [εₓₓ, ε_yy, γₓy]` with the tensor shear `ε₁₂ = γₓy/2`.

For a symmetric 2×2 tensor the eigenvalues have the closed form

    ε₁,₂ = ½ trε ± R ,   trε = εₓₓ + ε_yy ,   R = √( (½(εₓₓ−ε_yy))² + ε₁₂² ) ,

so `ε1 ≥ ε2` (the major/minor principal strains).  `R` is the "radius" of the
strain Mohr circle; `R = 0` means the two principal strains coincide.
"""
function principal_strains(ε)
    trε = ε[1] + ε[2]
    a   = 0.5 * (ε[1] - ε[2])     # half the in-plane strain difference
    ε12 = 0.5 * ε[3]              # tensor shear strain = γₓy / 2
    R   = hypot(a, ε12)          # √(a² + ε₁₂²) ≥ 0, computed without overflow
    return 0.5 * trε + R, 0.5 * trε - R
end

"""
    driving_energy_spectral(λ, μ, ε) -> ψ⁺

The **positive (tensile) elastic energy density** ψ⁺ — the crack driving force
of Miehe's spectral model,

    ψ⁺ = ½ λ ⟨trε⟩₊²  +  μ ( ⟨ε₁⟩₊² + ⟨ε₂⟩₊² ) .

This is the quantity whose running maximum in time becomes the irreversible
history field H (see `update_history!`).  It is `0` whenever the point is in
pure compression (all principal strains ≤ 0 and trε ≤ 0), so compression never
drives damage.
"""
function driving_energy_spectral(λ, μ, ε)
    ε1, ε2 = principal_strains(ε)
    trε    = ε1 + ε2
    return 0.5 * λ * pos(trε)^2 + μ * (pos(ε1)^2 + pos(ε2)^2)
end

"""
    spectral_stress_tangent(λ, μ, ε) -> (σ⁺, σ⁻, C⁺, C⁻)

Split stresses and consistent tangents at a strain `ε` (ENGINEERING form
`[εₓₓ, ε_yy, γₓy]`), for the Miehe spectral decomposition:

    σ⁺ = λ⟨trε⟩₊ I + 2μ ε⁺ ,     C⁺ = λ H(trε) I⊗I + 2μ ℙ⁺ ,
    σ⁻ = λ⟨trε⟩₋ I + 2μ ε⁻ ,     C⁻ = λ H(−trε) I⊗I + 2μ (𝕀 − ℙ⁺) .

Returned in engineering convention so they drop straight into the assembler:
`σ` is `[σₓₓ, σ_yy, σₓy]` (length 3) and `C` is a symmetric 3×3 matrix.

The degraded stress and tangent are then formed by the caller as
`σ = g(φ) σ⁺ + σ⁻` and `C = g(φ) C⁺ + C⁻`.

Method (all in Mandel form, then mapped to engineering — see the file header):
the fourth-order positive projection ℙ⁺ = ∂ε⁺/∂ε is built from the eigen-
projections m̂₁, m̂₂ of the strain and the difference quotient β of ⟨·⟩₊,

    ℙ⁺ = H(ε₁) m̂₁m̂₁ᵀ + H(ε₂) m̂₂m̂₂ᵀ + β (𝕀 − m̂₁m̂₁ᵀ − m̂₂m̂₂ᵀ) ,
    β  = (⟨ε₁⟩₊ − ⟨ε₂⟩₊)/(ε₁ − ε₂)   ∈ [0, 1]      (→ H(ε₁) as ε₂ → ε₁) .

The third term is the "shear" contribution that appears when the eigenvectors
rotate; β is bounded in [0,1] so ℙ⁺ (and hence C⁺, C⁻) is symmetric positive
semi-definite — the mechanical tangent stays well-posed.
"""
function spectral_stress_tangent(λ, μ, ε)
    εxx, εyy, γxy = ε[1], ε[2], ε[3]
    ε12 = 0.5 * γxy                 # tensor shear
    trε = εxx + εyy
    ctr = 0.5 * trε                 # centre of the strain Mohr circle
    a   = 0.5 * (εxx - εyy)
    R   = hypot(a, ε12)             # radius of the Mohr circle
    ε1  = ctr + R                   # major principal strain
    ε2  = ctr - R                   # minor principal strain

    # Mandel strain vector  ê = [εₓₓ, ε_yy, √2 ε₁₂].
    ε_m = [εxx, εyy, SQRT2 * ε12]

    # Eigenprojections of the 2×2 strain, as Mandel 3-vectors m̂₁, m̂₂.
    # For distinct eigenvalues,  M₁ = (ε − ε₂I)/(ε₁ − ε₂),  M₂ = I − M₁,
    # giving the closed forms below (all entries are bounded because |a|,|ε₁₂| ≤ R).
    # When ε₁ = ε₂ (R = 0) the strain is spherical: use m̂₁ = m̂₂ = ½ I and
    # β = H(ctr) — plugging these into the ℙ⁺ formula correctly yields H(ctr) 𝕀.
    if R > 0.0
        M1 = [(a + R) / (2R), (R - a) / (2R), SQRT2 * ε12 / (2R)]
        M2 = IVEC .- M1
        β  = (pos(ε1) - pos(ε2)) / (2R)          # difference quotient of ⟨·⟩₊, in [0,1]
    else
        M1 = 0.5 .* IVEC
        M2 = 0.5 .* IVEC
        β  = heaviside(ctr)
    end

    # Positive / negative strain in Mandel form:  ε⁺ = ⟨ε₁⟩₊ M₁ + ⟨ε₂⟩₊ M₂.
    ε_pos = pos(ε1) .* M1 .+ pos(ε2) .* M2
    ε_neg = ε_m .- ε_pos

    # Positive projection tensor ℙ⁺ (3×3, symmetric) — see the docstring.
    P1   = M1 * M1'
    P2   = M2 * M2'
    Ppos = heaviside(ε1) .* P1 .+ heaviside(ε2) .* P2 .+ β .* (IMAT .- P1 .- P2)

    # Split stresses in Mandel form:  σ± = λ⟨trε⟩± I + 2μ ε±.
    s_pos = (λ * pos(trε)) .* IVEC .+ (2μ) .* ε_pos
    s_neg = (λ * neg(trε)) .* IVEC .+ (2μ) .* ε_neg

    # Split tangents in Mandel form:  C± = λ H(±trε) I⊗I + 2μ ℙ± ,  ℙ⁻ = 𝕀 − ℙ⁺.
    IxI  = IVEC * IVEC'
    C_pos_m = (λ * heaviside(trε)) .* IxI .+ (2μ) .* Ppos
    C_neg_m = (λ * heaviside(-trε)) .* IxI .+ (2μ) .* (IMAT .- Ppos)

    # Convert Mandel → engineering:  σ_eng = T σ̂,  C_eng = T Ĉ T,  T = diag(1,1,1/√2).
    σ_pos = T_M2E .* s_pos
    σ_neg = T_M2E .* s_neg
    C_pos = (T_M2E .* C_pos_m) .* T_M2E'
    C_neg = (T_M2E .* C_neg_m) .* T_M2E'
    return σ_pos, σ_neg, C_pos, C_neg
end