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

# ===========================================================================
#  Q4 finite element: shape functions, isoparametric map, quadrature
# ===========================================================================
#
#  Reference element is the square (ξ,η) ∈ [-1,1]².  Bilinear shape functions
#
#      N₁ = ¼(1-ξ)(1-η)      N₂ = ¼(1+ξ)(1-η)
#      N₄ = ¼(1-ξ)(1+η)      N₃ = ¼(1+ξ)(1+η)
#
#  map the reference element onto a physical cell,  x(ξ,η) = Σₐ Nₐ(ξ,η) xₐ.

# The four 2×2 Gauss points (ξ,η) on the reference square [-1,1]².  Integrating
# a function means summing its value at these points times detJ (each weight = 1).
const GAUSS = let a = 1 / √3
    ((-a, -a), (a, -a), (a, a), (-a, a))
end

"""
    shape(ξ, η) -> (N, ∂N∂ξ, ∂N∂η)

The four bilinear shape functions `N` and their derivatives with respect to the
reference coordinates `ξ`, `η`, all evaluated at the point `(ξ,η)`.  Each is a
length-4 vector, one entry per corner node.
"""
function shape(ξ, η)
    N    = 0.25 .* [(1 - ξ) * (1 - η), (1 + ξ) * (1 - η), (1 + ξ) * (1 + η), (1 - ξ) * (1 + η)]
    ∂N∂ξ = 0.25 .* [-(1 - η),           (1 - η),           (1 + η),          -(1 + η)]
    ∂N∂η = 0.25 .* [-(1 - ξ),          -(1 + ξ),           (1 + ξ),           (1 - ξ)]
    return N, ∂N∂ξ, ∂N∂η
end

"""
Evaluate the element at Gauss point (ξ,η) for a cell whose 4 node coordinates
are the rows of `X` (4×2).  Returns

    N    shape function values                              (length 4)
    ∇N   physical gradients, row a = [∂Nₐ/∂x  ∂Nₐ/∂y]        (4×2)
    detJ Jacobian determinant  (so that  dΩ = detJ · weight)

The chain rule gives  [∂N/∂ξ ∂N/∂η] = ∇N · J  with the Jacobian
J = ∂(x,y)/∂(ξ,η),  hence  ∇N = [∂N/∂ξ ∂N/∂η] · J⁻¹.
"""
function element(X, ξ, η)
    N, ∂N∂ξ, ∂N∂η = shape(ξ, η)
    ref∇N = [∂N∂ξ ∂N∂η]          # 4×2: row a = [∂Nₐ/∂ξ, ∂Nₐ/∂η]  (reference gradients)
    J     = X' * ref∇N           # 2×2 Jacobian J = ∂(x,y)/∂(ξ,η) of the isoparametric map
    detJ  = det(J)               # local area scaling; must be positive
    detJ > 0 || error("Non-positive Jacobian determinant — check node ordering.")
    ∇N = ref∇N / J               # 4×2 physical gradients (right-division solves ∇N·J = ref∇N)
    return N, ∇N, detJ
end

"""
Strain–displacement matrix B (3×8) for engineering strain

    [εₓₓ, ε_yy, γₓy] = B · uₑ ,      γₓy = ∂u/∂y + ∂v/∂x ,

where uₑ = [uₓ¹ u_y¹ uₓ² u_y² uₓ³ u_y³ uₓ⁴ u_y⁴].
"""
function strain_matrix(∇N)
    B = zeros(3, 8)
    for a in 1:4
        B[1, 2a-1] = ∇N[a, 1]     # εₓₓ = ∂u/∂x
        B[2, 2a  ] = ∇N[a, 2]     # ε_yy = ∂v/∂y
        B[3, 2a-1] = ∇N[a, 2]     # γₓy = ∂u/∂y + ∂v/∂x
        B[3, 2a  ] = ∇N[a, 1]
    end
    return B
end

"Global displacement dofs [uₓ¹,u_y¹,…,uₓ⁴,u_y⁴] for a cell's 4 `nodes`."
function cell_dofs(nodes)
    d = Vector{Int}(undef, 8)
    for a in 1:4
        d[2a-1] = xdof(nodes[a])
        d[2a  ] = ydof(nodes[a])
    end
    return d
end