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

# ===========================================================================
#  Post-processing and output
# ===========================================================================
#
#  This file turns the raw solution vectors (`u`, `φ`) into things a user reads:
#
#    * scalar monitors per load step  → a CSV load history (`HistoryOutput.txt`)
#    * full nodal fields per snapshot  → a ParaView time series made of
#         - one `.vtu` file per saved step   (VTK "UnstructuredGrid", XML), and
#         - one `.pvd` file                  (VTK "Collection") that lists the
#           `.vtu` files and gives each of them a "time" value so ParaView can
#           animate the loading.
#
#  Every file is written here by hand as plain text — no external package.

# ---------------------------------------------------------------------------
#  Scalar monitors (the columns of the CSV load history)
# ---------------------------------------------------------------------------

"""
    reaction_top(mesh, f_int) -> Float64

Vertical reaction force carried by the top edge: the sum of the vertical
components of the internal force vector `f_int = Σₑ ∫ Bᵀ σ dΩ` over the
top-edge nodes.  At equilibrium the internal force vanishes on the free dofs,
so this sum is exactly the external force the supports must supply — the
"force" plotted against the imposed displacement in a force–displacement curve.

(We use the assembled `f_int` rather than `Kᵤ·u`, because with the
tension/compression split the stress is nonlinear in the strain and the
identity `f_int = Kᵤ·u` no longer holds.)
"""
reaction_top(mesh::Mesh, f_int) = sum(f_int[ydof(n)] for n in mesh.top)

"""
    crack_tip_x(mesh, φ; threshold=0.95) -> Float64

Estimate the crack-tip position as the right-most `x` coordinate of any node
whose damage `φ` exceeds `threshold` (nearly broken).  Returns `0.0` when no
node is damaged yet (the pre-crack reaches `x = 0`).
"""
function crack_tip_x(mesh::Mesh, φ; threshold = 0.95)
    tip = 0.0
    for n in 1:nnodes(mesh)
        if φ[n] ≥ threshold                 # this node is essentially broken
            tip = max(tip, mesh.coords[n, 1])
        end
    end
    return tip
end

"""
    max_free_phi(φ, crack_nodes) -> Float64

Largest damage value **outside** the prescribed initial crack.  The pre-crack
nodes are pinned to `φ = 1`, so they are excluded; what remains is the useful
monitor of how far damage has *grown* into the intact material.
"""
function max_free_phi(φ, crack_nodes)
    iscrack = falses(length(φ))              # mark the prescribed pre-crack nodes
    for n in crack_nodes
        iscrack[n] = true
    end
    return maximum((φ[n] for n in 1:length(φ) if !iscrack[n]); init = 0.0)
end

"Write the comment line naming the columns of the CSV load history."
function write_history_header(path)
    open(path, "w") do io
        println(io, "# step, imposed_top_uy, reaction_y_top, crack_tip_x, max_phi, max_free_phi, staggered_iterations, converged")
    end
end

"""
    log_step(problem, step, top_uy, u, φ, f_int, iters, converged)

Append one row of scalar monitors for the current load step to the CSV history:
step index, imposed top displacement, top reaction, crack tip, max damage,
max damage outside the pre-crack, staggered-iteration count, and whether the
staggered loop converged.  `f_int` is the assembled internal force at the
converged displacement (used for the reaction).
"""
function log_step(problem::Problem, step, top_uy, u, φ, f_int, iters, converged)
    R   = reaction_top(problem.mesh, f_int)                  # top-edge reaction force
    tip = crack_tip_x(problem.mesh, φ)                       # crack-tip x position
    mfp = max_free_phi(φ, problem.crack_nodes)               # grown damage monitor
    open(problem.histfile, "a") do io
        @printf(io, "%d, %.16e, %.16e, %.16e, %.16e, %.16e, %d, %s\n",
                step, top_uy, R, tip, maximum(φ), mfp, iters, converged)
    end
end

# ---------------------------------------------------------------------------
#  VTU: one snapshot of the mesh together with its nodal fields
# ---------------------------------------------------------------------------
#
#  A `.vtu` file is XML with this skeleton (indentation only for readability):
#
#    <VTKFile type="UnstructuredGrid">
#      <UnstructuredGrid>
#        <Piece NumberOfPoints=.. NumberOfCells=..>
#          <Points>   … node coordinates (x y z), z = 0 here …           </Points>
#          <Cells>    … connectivity, offsets, types …                   </Cells>
#          <PointData>… phase_field (scalar) and displacement (vector) … </PointData>
#        </Piece>
#      </UnstructuredGrid>
#    </VTKFile>
#
#  * connectivity  lists the node indices of every cell, **0-based** for VTK.
#  * offsets       are the running totals of nodes per cell: 4, 8, 12, …
#  * types         is the VTK cell code; 9 is VTK_QUAD (a 4-node quadrilateral).

"""
    write_vtu(path, mesh, u, φ)

Write one ParaView `.vtu` snapshot (XML, ASCII) holding the current nodal
damage `φ` and displacement `u` on `mesh`.
"""
function write_vtu(path, mesh::Mesh, u, φ)
    nn = nnodes(mesh)          # number of nodes  (points)
    ne = ncells(mesh)          # number of cells  (quads)
    open(path, "w") do io
        println(io, """<?xml version="1.0"?>""")
        println(io, """<VTKFile type="UnstructuredGrid" version="1.0" byte_order="LittleEndian">""")
        println(io, "  <UnstructuredGrid>")
        println(io, """    <Piece NumberOfPoints="$nn" NumberOfCells="$ne">""")

        # --- node coordinates:  x y z   (z = 0 for this 2D problem) ---
        println(io, "      <Points>")
        println(io, """        <DataArray type="Float64" NumberOfComponents="3" format="ascii">""")
        for n in 1:nn
            @printf(io, "          %.16e %.16e 0.0\n", mesh.coords[n, 1], mesh.coords[n, 2])
        end
        println(io, "        </DataArray>")
        println(io, "      </Points>")

        # --- cells: connectivity (0-based), offsets (cumulative), types (9=quad) ---
        println(io, "      <Cells>")
        println(io, """        <DataArray type="Int64" Name="connectivity" format="ascii">""")
        for e in 1:ne
            c = @view mesh.cells[e, :]
            @printf(io, "          %d %d %d %d\n", c[1] - 1, c[2] - 1, c[3] - 1, c[4] - 1)
        end
        println(io, "        </DataArray>")
        println(io, """        <DataArray type="Int64" Name="offsets" format="ascii">""")
        for e in 1:ne
            println(io, "          ", 4 * e)          # 4 nodes per quad → 4, 8, 12, …
        end
        println(io, "        </DataArray>")
        println(io, """        <DataArray type="UInt8" Name="types" format="ascii">""")
        for _ in 1:ne
            println(io, "          9")                # VTK_QUAD
        end
        println(io, "        </DataArray>")
        println(io, "      </Cells>")

        # --- nodal fields: damage (scalar) and displacement (3-vector, w=0) ---
        println(io, """      <PointData Scalars="phase_field" Vectors="displacement">""")
        println(io, """        <DataArray type="Float64" Name="phase_field" NumberOfComponents="1" format="ascii">""")
        for n in 1:nn
            @printf(io, "          %.16e\n", φ[n])
        end
        println(io, "        </DataArray>")
        println(io, """        <DataArray type="Float64" Name="displacement" NumberOfComponents="3" format="ascii">""")
        for n in 1:nn
            @printf(io, "          %.16e %.16e 0.0\n", u[xdof(n)], u[ydof(n)])
        end
        println(io, "        </DataArray>")
        println(io, "      </PointData>")

        println(io, "    </Piece>")
        println(io, "  </UnstructuredGrid>")
        println(io, "</VTKFile>")
    end
    return path
end

# ---------------------------------------------------------------------------
#  PVD: the collection that turns the .vtu snapshots into a time series
# ---------------------------------------------------------------------------

"""
    write_pvd(path, entries)

Write a ParaView `.pvd` collection.  `entries` is a list of `(time, filename)`
pairs: `filename` is the name of a `.vtu` file **relative to this `.pvd`
file's own folder**, and `time` is the value ParaView uses on its time axis
(here the imposed top displacement, i.e. the quasi-static "loading time").
Opening the `.pvd` in ParaView loads the whole animation at once.
"""
function write_pvd(path, entries)
    open(path, "w") do io
        println(io, """<?xml version="1.0"?>""")
        println(io, """<VTKFile type="Collection" version="1.0" byte_order="LittleEndian">""")
        println(io, "  <Collection>")
        for (time, file) in entries
            @printf(io, "    <DataSet timestep=\"%.16e\" group=\"\" part=\"0\" file=\"%s\"/>\n",
                    time, file)
        end
        println(io, "  </Collection>")
        println(io, "</VTKFile>")
    end
    return path
end