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

# ===========================================================================
#  Structured Q4 mesh of a rectangle  Ω = [x0,x1] × [y0,y1]
# ===========================================================================
#
#  Grid points are numbered row by row, x running fastest:
#
#        (ny)(nx+1)+1  ...............  (ny+1)(nx+1)
#             :                              :
#          nx+2   nx+3   ...............   2(nx+1)
#             1      2    ...............     nx+1
#
#  A cell (Q4 quadrilateral) stores its four corner nodes counter-clockwise:
#
#        4 ---- 3
#        |      |
#        1 ---- 2

struct Mesh
    coords::Matrix{Float64}   # node coordinates; row n = (x, y) of node n. size (nnode, 2)
    cells::Matrix{Int}        # connectivity; row e = the 4 node indices of cell e. size (ncell, 4)
    nx::Int                   # number of cells along x
    ny::Int                   # number of cells along y
    top::Vector{Int}          # node indices on the top edge    (y = y1)
    bottom::Vector{Int}       # node indices on the bottom edge (y = y0)
    left::Vector{Int}         # node indices on the left edge   (x = x0)
    right::Vector{Int}        # node indices on the right edge  (x = x1)
end

nnodes(m::Mesh) = size(m.coords, 1)     # total number of nodes
ncells(m::Mesh) = size(m.cells, 1)      # total number of cells

# Displacement degrees of freedom.  The global vector u interleaves the two
# components per node, so node n owns dofs (2n-1, 2n) = (uₓ, u_y).
xdof(n) = 2n - 1        # index of the x-displacement of node n in the global u
ydof(n) = 2n           # index of the y-displacement of node n in the global u

"Global node index of grid point (i,j),  i ∈ 1:nx+1 (column, x),  j ∈ 1:ny+1 (row, y)."
grid_node(i, j, nx) = (j - 1) * (nx + 1) + i

"""
    rectangle_mesh(nx, ny; x0, x1, y0, y1) -> Mesh

Build a structured mesh of `nx`×`ny` quadrilateral cells covering the rectangle
`[x0,x1] × [y0,y1]`.  Node coordinates are placed on a regular grid and the four
boundary node sets are collected by comparing coordinates to the edges.
"""
function rectangle_mesh(nx, ny; x0 = -0.5, x1 = 0.5, y0 = -0.5, y1 = 0.5)
    xs = range(x0, x1, length = nx + 1)      # the nx+1 grid x-coordinates
    ys = range(y0, y1, length = ny + 1)      # the ny+1 grid y-coordinates

    # Place every node on the (nx+1)×(ny+1) grid.
    coords = Matrix{Float64}(undef, (nx + 1) * (ny + 1), 2)
    for j in 1:ny+1, i in 1:nx+1
        n = grid_node(i, j, nx)
        coords[n, 1] = xs[i]
        coords[n, 2] = ys[j]
    end

    # Connect each cell to its four corner nodes, counter-clockwise.
    cells = Matrix{Int}(undef, nx * ny, 4)
    e = 0                                     # running cell index
    for j in 1:ny, i in 1:nx
        e += 1
        cells[e, 1] = grid_node(i,     j,     nx)   # lower-left
        cells[e, 2] = grid_node(i + 1, j,     nx)   # lower-right
        cells[e, 3] = grid_node(i + 1, j + 1, nx)   # upper-right
        cells[e, 4] = grid_node(i,     j + 1, nx)   # upper-left
    end

    # Collect boundary nodes by testing each coordinate against the four edges.
    tol = 1e-9 * max(x1 - x0, y1 - y0)        # geometric tolerance for "on the edge"
    onx(n, x) = abs(coords[n, 1] - x) ≤ tol
    ony(n, y) = abs(coords[n, 2] - y) ≤ tol
    N = size(coords, 1)
    top    = [n for n in 1:N if ony(n, y1)]
    bottom = [n for n in 1:N if ony(n, y0)]
    left   = [n for n in 1:N if onx(n, x0)]
    right  = [n for n in 1:N if onx(n, x1)]

    return Mesh(coords, cells, nx, ny, top, bottom, left, right)
end

"""
Nodes of a straight horizontal pre-crack on the line y = `crack_y`, running from
the left edge up to x = `tip_x`.  These nodes are held at φ = 1 (fully broken).
"""
function crack_line_nodes(mesh::Mesh; crack_y = 0.0, tip_x = 0.0)
    dy  = (maximum(mesh.coords[:, 2]) - minimum(mesh.coords[:, 2])) / mesh.ny  # vertical node spacing
    tol = 0.25 * dy    # accept nodes within a quarter-spacing of the crack line
    # A node is on the pre-crack if it is left of the tip and near the crack line.
    nodes = [n for n in 1:nnodes(mesh)
             if mesh.coords[n, 1] ≤ tip_x + 1e-12 && abs(mesh.coords[n, 2] - crack_y) ≤ tol]
    isempty(nodes) && error("No crack nodes found — use an even ny so nodes lie on y = crack_y.")
    return sort(nodes)
end