vem_phase_field_3d.py

"""
3D Phase-Field Fracture VEM on Polyhedral Meshes.

Extends vem_phase_field.py (2D) to full 3D polyhedral VEM.
Staggered (alternating minimization) with DI-dependent G_c.

World's first: 3D polyhedral VEM × phase-field × biofilm detachment.

Phase-field model:
  Momentum:  ∇·[g(d)·σ₀] = 0,  g(d) = (1-d)² + k
  Phase-field: G_c·l₀·Δd - G_c/l₀·d + 2(1-d)·ψ⁺ = 0

References:
  - Aldakheel, Hudobivnik, Hussein, Wriggers (2018) CMAME 341
  - Gain, Talischi, Paulino (2014) 3D VEM elasticity
"""

import numpy as np
import scipy.sparse as sp
from scipy.sparse.linalg import spsolve

from vem_3d import (isotropic_3d, traction_from_voigt, face_normal_area,
                    polyhedron_volume)


# ── Constitutive ─────────────────────────────────────────────────────────

def compute_Gc(DI, Gc_max=0.5, Gc_min=0.01, n=2):
    """G_c(DI) = G_c_min + (G_c_max - G_c_min)·(1 - DI)^n"""
    return Gc_min + (Gc_max - Gc_min) * (1.0 - np.clip(DI, 0, 1)) ** n

def compute_E_from_DI(DI, E_max=1000.0, E_min=30.0, n=2):
    """E(DI) = E_min + (E_max - E_min)·(1 - DI)^n"""
    return E_min + (E_max - E_min) * (1.0 - np.clip(DI, 0, 1)) ** n


# ── 3D Spectral Decomposition ────────────────────────────────────────────

def spectral_decomposition_3d(eps_voigt):
    """
    Spectral decomposition of 3D strain tensor.

    Parameters
    ----------
    eps_voigt : (6,) [ε_xx, ε_yy, ε_zz, γ_yz, γ_xz, γ_xy] (engineering shear)

    Returns
    -------
    eps_plus, eps_minus : (6,) Voigt arrays
    """
    eps_tensor = np.array([
        [eps_voigt[0],       0.5*eps_voigt[5], 0.5*eps_voigt[4]],
        [0.5*eps_voigt[5],   eps_voigt[1],     0.5*eps_voigt[3]],
        [0.5*eps_voigt[4],   0.5*eps_voigt[3], eps_voigt[2]],
    ])

    # Guard against NaN/Inf from singular solve
    if not np.all(np.isfinite(eps_tensor)):
        return np.zeros(6), np.zeros(6)

    eigvals, eigvecs = np.linalg.eigh(eps_tensor)

    eps_plus_t = np.zeros((3, 3))
    eps_minus_t = np.zeros((3, 3))
    for i in range(3):
        n = eigvecs[:, i:i+1]
        if eigvals[i] > 0:
            eps_plus_t += eigvals[i] * (n @ n.T)
        else:
            eps_minus_t += eigvals[i] * (n @ n.T)

    eps_plus = np.array([
        eps_plus_t[0,0], eps_plus_t[1,1], eps_plus_t[2,2],
        2.0*eps_plus_t[1,2], 2.0*eps_plus_t[0,2], 2.0*eps_plus_t[0,1]
    ])
    eps_minus = np.array([
        eps_minus_t[0,0], eps_minus_t[1,1], eps_minus_t[2,2],
        2.0*eps_minus_t[1,2], 2.0*eps_minus_t[0,2], 2.0*eps_minus_t[0,1]
    ])
    return eps_plus, eps_minus


def compute_psi_plus_3d(eps_voigt, E, nu):
    """
    Tensile elastic energy density ψ⁺ in 3D.
    ψ⁺ = λ/2·<tr(ε)>₊² + μ·tr(ε⁺·ε⁺)
    """
    mu = E / (2.0 * (1.0 + nu))
    lam = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu))

    eps_plus, _ = spectral_decomposition_3d(eps_voigt)
    tr_eps = eps_voigt[0] + eps_voigt[1] + eps_voigt[2]
    tr_plus = max(tr_eps, 0.0)

    eps_plus_t = np.array([
        [eps_plus[0],       0.5*eps_plus[5], 0.5*eps_plus[4]],
        [0.5*eps_plus[5],   eps_plus[1],     0.5*eps_plus[3]],
        [0.5*eps_plus[4],   0.5*eps_plus[3], eps_plus[2]],
    ])
    psi = lam / 2.0 * tr_plus**2 + mu * np.trace(eps_plus_t @ eps_plus_t)
    return max(psi, 0.0)


# ── 3D Degraded Elasticity VEM Assembly ──────────────────────────────────

def assemble_degraded_elasticity_3d(vertices, cells, cell_faces,
                                     E_field, nu, d_field, k_residual=1e-6):
    """
    3D VEM elasticity with degradation g(d) = (1-d)² + k.
    Returns K_global (sparse) and elem_strain_data for ψ⁺ computation.
    """
    n_nodes = len(vertices)
    n_dofs = 3 * n_nodes
    n_polys = 12

    rows_list, cols_list, vals_list = [], [], []
    elem_strain_data = []

    strain_ids = np.array([
        [1, 0, 0, 0, 0, 0],
        [0, 1, 0, 0, 0, 0],
        [0, 0, 1, 0, 0, 0],
        [0, 0, 0, 2, 0, 0],
        [0, 0, 0, 0, 2, 0],
        [0, 0, 0, 0, 0, 2],
    ], dtype=float)

    for el_id in range(len(cells)):
        vert_ids = cells[el_id].astype(int)
        coords = vertices[vert_ids]
        faces = cell_faces[el_id]
        n_v = len(vert_ids)
        n_el = 3 * n_v

        vmap = {int(g): loc for loc, g in enumerate(vert_ids)}

        E_el = E_field[el_id] if hasattr(E_field, '__len__') else E_field

        # Average d over element vertices
        d_avg = np.mean(d_field[vert_ids])
        g_d = (1.0 - d_avg)**2 + k_residual
        E_degraded = E_el * g_d

        C = isotropic_3d(E_degraded, nu)

        centroid = coords.mean(axis=0)
        h = max(np.linalg.norm(coords[i] - coords[j])
                for i in range(n_v) for j in range(i+1, n_v))
        vol = polyhedron_volume(vertices, faces)
        if vol < 1e-20:
            elem_strain_data.append(None)
            continue

        xc, yc, zc = centroid

        # D matrix (n_el × 12)
        D = np.zeros((n_el, n_polys))
        for i in range(n_v):
            dx = (coords[i, 0] - xc) / h
            dy = (coords[i, 1] - yc) / h
            dz = (coords[i, 2] - zc) / h
            r = 3 * i
            D[r,   :] = [1, 0, 0, 0,   dz, -dy, dx, 0,  0,  0,  dz, dy]
            D[r+1, :] = [0, 1, 0, -dz, 0,  dx,  0,  dy, 0,  dz, 0,  dx]
            D[r+2, :] = [0, 0, 1, dy,  -dx, 0,  0,  0,  dz, dy, dx, 0]

        # B matrix (12 × n_el)
        B = np.zeros((n_polys, n_el))
        for i in range(n_v):
            B[0, 3*i]     = 1.0 / n_v
            B[1, 3*i + 1] = 1.0 / n_v
            B[2, 3*i + 2] = 1.0 / n_v

        for face in faces:
            face_int = face.astype(int)
            pts = vertices[face_int]
            n_f, A_f = face_normal_area(pts)
            fc = pts.mean(axis=0)
            if np.dot(n_f, fc - centroid) < 0:
                n_f = -n_f
            k_f = len(face_int)

            for gv in face_int:
                if gv not in vmap:
                    continue
                li = vmap[gv]
                w = A_f / k_f

                wrot = w / (2.0 * vol)
                B[3, 3*li + 1] += -wrot * n_f[2]
                B[3, 3*li + 2] +=  wrot * n_f[1]
                B[4, 3*li + 0] +=  wrot * n_f[2]
                B[4, 3*li + 2] += -wrot * n_f[0]
                B[5, 3*li + 0] += -wrot * n_f[1]
                B[5, 3*li + 1] +=  wrot * n_f[0]

                for alpha in range(6):
                    eps_a = strain_ids[alpha] / h
                    sigma_a = C @ eps_a
                    t_f = traction_from_voigt(sigma_a, n_f)
                    B[6 + alpha, 3*li + 0] += w * t_f[0]
                    B[6 + alpha, 3*li + 1] += w * t_f[1]
                    B[6 + alpha, 3*li + 2] += w * t_f[2]

        G = B @ D
        try:
            projector = np.linalg.solve(G, B)
        except np.linalg.LinAlgError:
            elem_strain_data.append(None)
            continue

        G_tilde = G.copy()
        G_tilde[:6, :] = 0.0
        K_cons = projector.T @ G_tilde @ projector

        I_minus_PiD = np.eye(n_el) - D @ projector
        trace_k = np.trace(K_cons)
        stab_param = trace_k / n_el if trace_k > 0 else E_degraded
        K_stab = stab_param * (I_minus_PiD.T @ I_minus_PiD)
        K_local = K_cons + K_stab

        elem_strain_data.append({
            'vert_ids': vert_ids,
            'projector': projector,
            'h': h,
            'E_el': E_el,  # undegraded
        })

        gdofs = np.zeros(n_el, dtype=int)
        for i in range(n_v):
            gdofs[3*i]     = 3 * vert_ids[i]
            gdofs[3*i + 1] = 3 * vert_ids[i] + 1
            gdofs[3*i + 2] = 3 * vert_ids[i] + 2

        gi, gj = np.meshgrid(gdofs, gdofs, indexing='ij')
        rows_list.append(gi.ravel())
        cols_list.append(gj.ravel())
        vals_list.append(K_local.ravel())

    K_global = sp.coo_matrix(
        (np.concatenate(vals_list), (np.concatenate(rows_list), np.concatenate(cols_list))),
        shape=(n_dofs, n_dofs)
    ).tocsr()

    return K_global, elem_strain_data


def compute_element_strains_3d(u, elem_strain_data):
    """
    Recover element-average strains from 3D displacement via VEM projector.
    Returns list of (6,) Voigt strain arrays [ε_xx, ε_yy, ε_zz, γ_yz, γ_xz, γ_xy].
    """
    strains = []
    for ed in elem_strain_data:
        if ed is None:
            strains.append(np.zeros(6))
            continue
        vert_ids = ed['vert_ids']
        n_v = len(vert_ids)
        proj = ed['projector']
        h = ed['h']

        gdofs = np.zeros(3 * n_v, dtype=int)
        for i in range(n_v):
            gdofs[3*i]     = 3 * vert_ids[i]
            gdofs[3*i + 1] = 3 * vert_ids[i] + 1
            gdofs[3*i + 2] = 3 * vert_ids[i] + 2

        u_local = u[gdofs]
        c = proj @ u_local

        # c[6:12] correspond to strain modes scaled by 1/h
        eps_xx  = c[6]  / h
        eps_yy  = c[7]  / h
        eps_zz  = c[8]  / h
        gamma_yz = c[9]  / h  # engineering shear
        gamma_xz = c[10] / h
        gamma_xy = c[11] / h

        strains.append(np.array([eps_xx, eps_yy, eps_zz, gamma_yz, gamma_xz, gamma_xy]))
    return strains


# ── 3D Scalar VEM Assembly (Phase-Field Equation) ────────────────────────

def assemble_scalar_vem_3d(vertices, cells, cell_faces,
                            diffusion_field, reaction_field, source_field):
    """
    Assemble 3D scalar VEM system: -κ·Δd + r·d = s

    Scalar polynomial basis P_1: {1, (x-xc)/h, (y-yc)/h, (z-zc)/h} — 4 functions.
    B matrix rows 1-3 computed from face integrals (∇p · n).
    """
    n_nodes = len(vertices)
    n_polys = 4

    rows_list, cols_list, vals_list = [], [], []
    F_global = np.zeros(n_nodes)

    for el_id in range(len(cells)):
        vert_ids = cells[el_id].astype(int)
        coords = vertices[vert_ids]
        faces = cell_faces[el_id]
        n_v = len(vert_ids)

        vmap = {int(g): loc for loc, g in enumerate(vert_ids)}

        kappa = diffusion_field[el_id] if hasattr(diffusion_field, '__len__') else diffusion_field
        react = reaction_field[el_id] if hasattr(reaction_field, '__len__') else reaction_field
        source = source_field[el_id] if hasattr(source_field, '__len__') else source_field

        centroid = coords.mean(axis=0)
        h = max(np.linalg.norm(coords[i] - coords[j])
                for i in range(n_v) for j in range(i+1, n_v))
        vol = polyhedron_volume(vertices, faces)
        if vol < 1e-20:
            continue

        # D matrix (n_v × 4)
        D = np.zeros((n_v, n_polys))
        D[:, 0] = 1.0
        for i in range(n_v):
            D[i, 1] = (coords[i, 0] - centroid[0]) / h
            D[i, 2] = (coords[i, 1] - centroid[1]) / h
            D[i, 3] = (coords[i, 2] - centroid[2]) / h

        # B matrix (4 × n_v) using face integrals
        B = np.zeros((n_polys, n_v))
        B[0, :] = 1.0 / n_v  # constant mode: average

        # Gradient modes: B[α, i] = Σ_faces (∇p_α · n_f) · A_f / k_f
        for face in faces:
            face_int = face.astype(int)
            pts = vertices[face_int]
            n_f, A_f = face_normal_area(pts)
            fc = pts.mean(axis=0)
            if np.dot(n_f, fc - centroid) < 0:
                n_f = -n_f
            k_f = len(face_int)

            for gv in face_int:
                if gv not in vmap:
                    continue
                li = vmap[gv]
                w = A_f / k_f

                # ∇p_1 = (1/h, 0, 0), ∇p_2 = (0, 1/h, 0), ∇p_3 = (0, 0, 1/h)
                B[1, li] += w * n_f[0] / h
                B[2, li] += w * n_f[1] / h
                B[3, li] += w * n_f[2] / h

        G = B @ D
        try:
            G_inv = np.linalg.inv(G)
        except np.linalg.LinAlgError:
            continue
        projector = G_inv @ B

        # Diffusion stiffness
        G_tilde = G.copy()
        G_tilde[0, :] = 0.0
        K_cons = projector.T @ G_tilde @ projector * kappa

        # Stabilization
        I_minus_PiD = np.eye(n_v) - D @ projector
        trace_k = np.trace(K_cons)
        stab_param = 0.5 * abs(trace_k) / max(n_v, 1) if trace_k > 0 else kappa * 0.01
        K_stab = stab_param * (I_minus_PiD.T @ I_minus_PiD)
        K_diff = K_cons + K_stab

        # Lumped mass matrix
        M_lumped = (vol / n_v) * np.eye(n_v)

        # Element stiffness and source
        K_el = K_diff + react * M_lumped
        F_el = source * (vol / n_v) * np.ones(n_v)

        # Assemble
        ii, jj = np.meshgrid(vert_ids, vert_ids, indexing='ij')
        rows_list.append(ii.ravel())
        cols_list.append(jj.ravel())
        vals_list.append(K_el.ravel())
        F_global[vert_ids] += F_el

    K_global = sp.coo_matrix(
        (np.concatenate(vals_list), (np.concatenate(rows_list), np.concatenate(cols_list))),
        shape=(n_nodes, n_nodes)
    ).tocsr()
    return K_global, F_global


# ── 3D Phase-Field VEM Solver ────────────────────────────────────────────

class PhaseFieldVEM3D:
    """
    Staggered 3D phase-field fracture solver on polyhedral VEM mesh.

    Alternating minimization:
      1. Fix d → solve displacement u (degraded 3D VEM)
      2. Fix u → solve phase-field d (3D scalar VEM)
      3. Enforce irreversibility: d_new ≥ d_old
    """

    def __init__(self, vertices, cells, cell_faces, E_field, nu, Gc_field, l0=None):
        self.vertices = vertices
        self.cells = cells
        self.cell_faces = cell_faces
        self.E_field = np.asarray(E_field, dtype=float)
        self.nu = nu
        self.Gc_field = np.asarray(Gc_field, dtype=float)

        n_el = len(cells)
        if l0 is None:
            diams = []
            for el_id in range(n_el):
                coords = vertices[cells[el_id].astype(int)]
                n_v = len(coords)
                h = max(np.linalg.norm(coords[i] - coords[j])
                        for i in range(n_v) for j in range(i+1, n_v))
                diams.append(h)
            self.l0 = np.mean(diams) * 0.5
        else:
            self.l0 = l0

        self.n_nodes = len(vertices)
        self.n_dofs = 3 * self.n_nodes
        self.n_el = n_el

        self.d = np.zeros(self.n_nodes)
        self.u = np.zeros(self.n_dofs)
        self.psi_history = np.zeros(n_el)

    def solve_displacement(self, bc_fixed_dofs, bc_vals, load_dofs=None, load_vals=None):
        """Solve 3D displacement with degraded stiffness."""
        K, elem_data = assemble_degraded_elasticity_3d(
            self.vertices, self.cells, self.cell_faces,
            self.E_field, self.nu, self.d
        )

        F = np.zeros(self.n_dofs)
        if load_dofs is not None and load_vals is not None:
            F[load_dofs] += load_vals

        u = np.zeros(self.n_dofs)
        bc_set = set(bc_fixed_dofs.tolist())
        internal = np.array([i for i in range(self.n_dofs) if i not in bc_set])
        u[bc_fixed_dofs] = bc_vals
        F -= K[:, bc_fixed_dofs].toarray() @ bc_vals

        K_ii = K[np.ix_(internal, internal)]
        try:
            u[internal] = spsolve(K_ii, F[internal])
        except Exception:
            pass

        self.u = u
        self._elem_strain_data = elem_data
        return u

    def compute_psi_plus_field(self):
        """Compute tensile energy density ψ⁺ per element and update history."""
        strains = compute_element_strains_3d(self.u, self._elem_strain_data)

        psi_plus = np.zeros(self.n_el)
        for i, eps in enumerate(strains):
            E_el = self.E_field[i] if hasattr(self.E_field, '__len__') else self.E_field
            psi = compute_psi_plus_3d(eps, E_el, self.nu)
            psi_plus[i] = max(psi, self.psi_history[i])

        self.psi_history = psi_plus.copy()
        return psi_plus

    def solve_phase_field(self, psi_plus):
        """Solve 3D phase-field equation via scalar VEM."""
        diffusion = self.Gc_field * self.l0
        reaction = self.Gc_field / self.l0 + 2.0 * psi_plus
        source = 2.0 * psi_plus

        K, F = assemble_scalar_vem_3d(
            self.vertices, self.cells, self.cell_faces,
            diffusion, reaction, source
        )

        try:
            d_new = spsolve(K, F)
        except Exception:
            d_new = self.d.copy()

        d_new = np.clip(d_new, 0.0, 1.0)
        d_new = np.maximum(d_new, self.d)
        self.d = d_new
        return d_new

    def run(self, bc_fixed_dofs, bc_vals, load_schedule,
            max_stagger=30, tol=1e-4, verbose=False):
        """
        Incremental loading with staggered phase-field solve.

        Parameters
        ----------
        load_schedule: list of (load_dofs, load_vals) per step

        Returns: list of snapshots
        """
        snapshots = []

        for step, (l_dofs, l_vals) in enumerate(load_schedule):
            d_old = self.d.copy()

            for stag_iter in range(max_stagger):
                self.solve_displacement(bc_fixed_dofs, bc_vals, l_dofs, l_vals)
                psi_plus = self.compute_psi_plus_field()
                d_new = self.solve_phase_field(psi_plus)

                d_change = np.linalg.norm(d_new - d_old) / max(np.linalg.norm(d_new), 1e-10)

                if verbose and stag_iter % 5 == 0:
                    print(f"  Step {step+1}, stagger {stag_iter}: "
                          f"|Δd|/|d| = {d_change:.2e}, "
                          f"max(d) = {np.max(self.d):.4f}")

                if d_change < tol:
                    break
                d_old = d_new.copy()

            ux = self.u[0::3]
            uy = self.u[1::3]
            uz = self.u[2::3]
            mag = np.sqrt(ux**2 + uy**2 + uz**2)

            snapshots.append({
                'step': step,
                'u': self.u.copy(),
                'd': self.d.copy(),
                'u_max': np.max(mag),
                'd_max': np.max(self.d),
                'd_mean': np.mean(self.d),
                'psi_max': np.max(psi_plus),
                'stagger_iters': stag_iter + 1,
                'n_cracked': int(np.sum(self.d > 0.9)),
            })

            if verbose:
                s = snapshots[-1]
                print(f"  → Step {step+1}: |u|_max={s['u_max']:.4e}, "
                      f"d_max={s['d_max']:.4f}, cracked={s['n_cracked']}")

        return snapshots


# ── Pipeline: Real SAPA Data → 3D Phase-Field ────────────────────────────

def run_3d_phase_field_real(tiff_path, save_dir, downsample=3):
    """Run 3D phase-field on real SAPA dual-species data."""
    import os, time
    from pipeline_3d_real import (load_tiff_3d, segment_colonies_3d,
                                   build_voronoi_mesh_3d)

    os.makedirs(save_dir, exist_ok=True)
    print("Loading 3D TIFF...")
    channels = load_tiff_3d(tiff_path, downsample=downsample)
    voxel_xy = 0.5 * downsample

    print("Segmenting colonies...")
    centroids, species_fracs = segment_colonies_3d(channels, quantile=0.90,
                                                     min_volume_voxels=5, merge_distance_voxels=6)
    n_colonies = len(centroids)
    print(f"  {n_colonies} colonies detected")

    print("Building 3D Voronoi mesh...")
    # Normalize centroids to [0,1]³ for Voronoi
    cents_phys = centroids * voxel_xy
    c_min = cents_phys.min(axis=0)
    c_max = cents_phys.max(axis=0)
    c_range = c_max - c_min
    c_range[c_range < 1e-10] = 1.0
    cents_norm = (cents_phys - c_min) / c_range
    cents_norm = np.clip(cents_norm, 0.05, 0.95)
    vertices_raw, cells_raw, cell_faces_raw, seed_to_cell = build_voronoi_mesh_3d(cents_norm)

    # Vertex cleanup (critical!)
    used_verts = set()
    for cv in cells_raw:
        used_verts.update(cv.astype(int).tolist())
    for cf_list in cell_faces_raw:
        for face in cf_list:
            used_verts.update(face.astype(int).tolist())
    used_sorted = sorted(used_verts)
    old_to_new = {old: new for new, old in enumerate(used_sorted)}
    vertices = vertices_raw[used_sorted]
    cells = [np.array([old_to_new[int(v)] for v in cv]) for cv in cells_raw]
    cell_faces = []
    for cf_list in cell_faces_raw:
        new_faces = []
        for face in cf_list:
            new_face = np.array([old_to_new[int(v)] for v in face.astype(int)])
            new_faces.append(new_face)
        cell_faces.append(new_faces)

    n_el = len(cells)
    n_nodes = len(vertices)
    print(f"  Mesh: {n_el} cells, {n_nodes} vertices (cleaned)")

    # DI, E, Gc per cell — assign from nearest colony
    from pipeline_3d_real import compute_DI_shannon
    DI_per_cell = np.full(n_el, 0.5)
    for i, c_verts in enumerate(cells):
        cx = vertices[c_verts.astype(int)].mean(axis=0)
        # Nearest colony in normalized coords
        dists = np.linalg.norm(cents_norm - cx, axis=1)
        nearest = np.argmin(dists)
        if nearest < len(species_fracs):
            phi = species_fracs[nearest]
            di_sh = compute_DI_shannon(phi)
            pathogenicity = np.dot(phi, np.linspace(0.3, 0.9, len(phi)))
            DI_per_cell[i] = np.clip(0.4 * di_sh + 0.6 * pathogenicity, 0, 1)
    E_per_cell = compute_E_from_DI(DI_per_cell)
    Gc_per_cell = compute_Gc(DI_per_cell)

    print(f"  DI: [{DI_per_cell.min():.3f}, {DI_per_cell.max():.3f}]")
    print(f"  E:  [{E_per_cell.min():.0f}, {E_per_cell.max():.0f}] Pa")
    print(f"  Gc: [{Gc_per_cell.min():.4f}, {Gc_per_cell.max():.4f}] J/m²")

    # BCs: bottom fixed, top loaded
    zmin, zmax = vertices[:, 2].min(), vertices[:, 2].max()
    tol_bc = (zmax - zmin) * 0.05
    bottom = np.where(vertices[:, 2] < zmin + tol_bc)[0]
    top = np.where(vertices[:, 2] > zmax - tol_bc)[0]

    bc_dofs = np.concatenate([3*bottom, 3*bottom+1, 3*bottom+2])
    bc_vals = np.zeros(len(bc_dofs))

    # Load schedule: shear + compression on top
    # Real data: E=74-682 Pa, Gc=0.03-0.34 → need substantial load
    # Scale: lf ~ sqrt(2*Gc_min*E_min) * domain_size * factor
    domain_L = max(zmax - zmin,
                    vertices[:, 0].max() - vertices[:, 0].min(),
                    vertices[:, 1].max() - vertices[:, 1].min())
    lf_max = np.sqrt(2 * Gc_per_cell.min() * E_per_cell.min()) * domain_L * 2.5
    print(f"  Load factor max: {lf_max:.2f} (domain_L={domain_L:.3f})")
    n_steps = 30
    load_schedule = []
    for step in range(n_steps):
        lf = (step + 1) / n_steps * lf_max
        if len(top) > 0:
            l_dofs = np.concatenate([3*top, 3*top+1, 3*top+2])
            l_vals = np.concatenate([
                np.full(len(top), lf / len(top)),         # x-shear
                np.full(len(top), 0.0),                    # y = 0
                np.full(len(top), -lf * 0.3 / len(top))   # z-compression
            ])
        else:
            l_dofs, l_vals = None, None
        load_schedule.append((l_dofs, l_vals))

    # Solve
    print("Running 3D PhaseFieldVEM...")
    t0 = time.time()
    solver = PhaseFieldVEM3D(vertices, cells, cell_faces,
                              E_per_cell, 0.35, Gc_per_cell)
    snapshots = solver.run(bc_dofs, bc_vals, load_schedule, verbose=True)
    elapsed = time.time() - t0
    print(f"  Done in {elapsed:.1f}s")

    # Visualization
    _plot_3d_results(vertices, cells, cell_faces, DI_per_cell, E_per_cell,
                     Gc_per_cell, snapshots, save_dir,
                     filename='phase_field_3d_real.png',
                     title_prefix='3D Phase-Field Fracture on Real Biofilm (SAPA)')

    return snapshots


def _plot_3d_results(verts, cells, cell_faces, DI, E, Gc, snapshots, save_dir,
                     filename='phase_field_3d.png', title_prefix='3D Phase-Field Fracture'):
    """Generate 3D phase-field result figure."""
    import os
    import matplotlib
    matplotlib.use('Agg')
    import matplotlib.pyplot as plt

    fig, axes = plt.subplots(2, 3, figsize=(18, 11))
    n_el = len(cells)

    # Helper: get cell centroid projections for 2D scatter
    cx = np.array([verts[c.astype(int)].mean(axis=0) for c in cells])

    # (a) DI field (xy projection)
    ax = axes[0, 0]
    sc = ax.scatter(cx[:, 0], cx[:, 1], c=DI, cmap='RdYlGn_r',
                    s=40, edgecolors='k', linewidth=0.3, vmin=0, vmax=1)
    fig.colorbar(sc, ax=ax, label='DI', shrink=0.8)
    ax.set_title(f'(a) DI [{DI.min():.2f}, {DI.max():.2f}]')
    ax.set_aspect('equal')

    # (b) G_c field
    ax = axes[0, 1]
    sc = ax.scatter(cx[:, 0], cx[:, 1], c=Gc, cmap='YlOrRd_r',
                    s=40, edgecolors='k', linewidth=0.3)
    fig.colorbar(sc, ax=ax, label='G_c [J/m²]', shrink=0.8)
    ax.set_title(f'(b) G_c [{Gc.min():.3f}, {Gc.max():.3f}]')
    ax.set_aspect('equal')

    # (c) E field
    ax = axes[0, 2]
    sc = ax.scatter(cx[:, 0], cx[:, 1], c=E, cmap='viridis',
                    s=40, edgecolors='k', linewidth=0.3)
    fig.colorbar(sc, ax=ax, label='E [Pa]', shrink=0.8)
    ax.set_title(f'(c) E [{E.min():.0f}, {E.max():.0f}] Pa')
    ax.set_aspect('equal')

    if snapshots:
        # Find last pre-catastrophic snapshot (d_max < 0.99)
        pre_cat = [s for s in snapshots if s['d_max'] < 0.99]
        show_snap = pre_cat[-1] if pre_cat else snapshots[-1]
        final = snapshots[-1]

        # (d) Phase-field d (last pre-catastrophic)
        ax = axes[1, 0]
        d_per_cell = np.array([np.mean(show_snap['d'][c.astype(int)]) for c in cells])
        sc = ax.scatter(cx[:, 0], cx[:, 1], c=d_per_cell, cmap='hot',
                        s=40, edgecolors='k', linewidth=0.3, vmin=0, vmax=1)
        fig.colorbar(sc, ax=ax, label='d (damage)', shrink=0.8)
        cat_step = show_snap['step'] + 1
        ax.set_title(f'(d) d at step {cat_step} (max={show_snap["d_max"]:.3f})')
        ax.set_aspect('equal')

        # (e) Damage evolution
        ax = axes[1, 1]
        steps = [s['step']+1 for s in snapshots]
        d_maxs = [s['d_max'] for s in snapshots]
        u_maxs = [s['u_max'] for s in snapshots]
        ax.plot(steps, d_maxs, 'ro-', ms=4, label='d_max')
        # Mark catastrophic transition
        if pre_cat and len(pre_cat) < len(snapshots):
            cat_idx = len(pre_cat)
            ax.axvline(x=steps[cat_idx], color='gray', ls='--', alpha=0.7,
                       label=f'catastrophic (step {steps[cat_idx]})')
            ax.legend(fontsize=8)
        ax2 = ax.twinx()
        # Only plot u_max for pre-catastrophic steps (post-cat is meaningless)
        pre_steps = steps[:len(pre_cat)] if pre_cat else steps
        pre_u = u_maxs[:len(pre_cat)] if pre_cat else u_maxs
        ax2.plot(pre_steps, pre_u, 'b^-', ms=4, label='|u|_max')
        ax2.set_ylabel('|u|_max', color='blue')
        ax.set_xlabel('Load step')
        ax.set_ylabel('d_max', color='red')
        ax.set_title('(e) Damage & displacement evolution')
        ax.grid(True, alpha=0.3)

        # (f) DI vs d correlation (pre-catastrophic)
        ax = axes[1, 2]
        sc = ax.scatter(DI, d_per_cell, c=Gc, cmap='coolwarm', s=30,
                        edgecolors='k', linewidth=0.3)
        fig.colorbar(sc, ax=ax, label='G_c [J/m²]', shrink=0.8)
        ax.set_xlabel('DI')
        ax.set_ylabel(f'd (step {cat_step})')
        ax.set_title('(f) DI vs damage (color=G_c)')
        ax.grid(True, alpha=0.3)

    cat_step_info = f'catastrophic at step {len(pre_cat)+1}' if pre_cat and len(pre_cat) < len(snapshots) else ''
    plt.suptitle(f'{title_prefix}\n'
                 f'{n_el} polyhedral VEM elements, {len(snapshots)} steps'
                 f'{", " + cat_step_info if cat_step_info else ""}',
                 fontsize=13)
    plt.tight_layout()
    out = os.path.join(save_dir, filename)
    plt.savefig(out, dpi=200, bbox_inches='tight')
    plt.close()
    print(f"Saved: {out}")


# ── Demo: Synthetic 3D Phase-Field ───────────────────────────────────────

def demo_3d_phase_field():
    """Demo on synthetic 3D Voronoi mesh with DI gradient."""
    import os, time
    from vem_3d_advanced import make_voronoi_mesh_3d

    print("Building 3D Voronoi mesh...")
    vertices_raw, cells_raw, cell_faces_raw = make_voronoi_mesh_3d(n_seeds=30, seed=42)

    # Vertex cleanup — remove floating vertices
    used_verts = set()
    for cv in cells_raw:
        used_verts.update(cv.astype(int).tolist())
    for cf_list in cell_faces_raw:
        for face in cf_list:
            used_verts.update(face.astype(int).tolist())
    used_sorted = sorted(used_verts)
    old_to_new = {old: new for new, old in enumerate(used_sorted)}
    vertices = vertices_raw[used_sorted]
    cells = [np.array([old_to_new[int(v)] for v in cv]) for cv in cells_raw]
    cell_faces = []
    for cf_list in cell_faces_raw:
        new_faces = []
        for face in cf_list:
            new_face = np.array([old_to_new[int(v)] for v in face.astype(int)])
            new_faces.append(new_face)
        cell_faces.append(new_faces)

    n_el = len(cells)
    n_nodes = len(vertices)
    print(f"  {n_el} cells, {n_nodes} vertices (cleaned from {len(vertices_raw)})")

    # DI gradient: center = dysbiotic, edges = commensal
    center = np.array([0.5, 0.5, 0.5])
    DI_per_cell = np.zeros(n_el)
    for i, c in enumerate(cells):
        cx = vertices[c.astype(int)].mean(axis=0)
        r = np.linalg.norm(cx - center)
        DI_per_cell[i] = np.clip(0.8 - 1.2 * r, 0.05, 0.95)

    E_field = compute_E_from_DI(DI_per_cell)
    Gc_field = compute_Gc(DI_per_cell)
    nu = 0.35

    print(f"  DI: [{DI_per_cell.min():.3f}, {DI_per_cell.max():.3f}]")
    print(f"  Gc: [{Gc_field.min():.4f}, {Gc_field.max():.4f}]")

    # BCs
    tol_bc = 0.05
    bottom = np.where(vertices[:, 2] < tol_bc)[0]
    top = np.where(vertices[:, 2] > 1.0 - tol_bc)[0]

    bc_dofs = np.concatenate([3*bottom, 3*bottom+1, 3*bottom+2])
    bc_vals = np.zeros(len(bc_dofs))

    n_steps = 30
    load_schedule = []
    for step in range(n_steps):
        lf = (step + 1) / n_steps * 4.5  # tuned for progressive damage
        if len(top) > 0:
            l_dofs = np.concatenate([3*top, 3*top+1, 3*top+2])
            l_vals = np.concatenate([
                np.full(len(top), lf / len(top)),
                np.full(len(top), 0.0),
                np.full(len(top), -lf * 0.3 / len(top))
            ])
        else:
            l_dofs, l_vals = None, None
        load_schedule.append((l_dofs, l_vals))

    print("Running 3D PhaseFieldVEM (synthetic)...")
    t0 = time.time()
    solver = PhaseFieldVEM3D(vertices, cells, cell_faces, E_field, nu, Gc_field)
    snapshots = solver.run(bc_dofs, bc_vals, load_schedule, verbose=True)
    elapsed = time.time() - t0
    print(f"  Done in {elapsed:.1f}s")

    save_dir = os.path.join(os.path.dirname(__file__), 'results', '3d_phase_field')
    os.makedirs(save_dir, exist_ok=True)
    _plot_3d_results(vertices, cells, cell_faces, DI_per_cell, E_field,
                     Gc_field, snapshots, save_dir,
                     filename='phase_field_3d_synthetic.png',
                     title_prefix='3D Phase-Field Fracture (Synthetic Voronoi)')

    final = snapshots[-1]
    print(f"\nFinal: d_max={final['d_max']:.4f}, |u|_max={final['u_max']:.4e}, "
          f"cracked={final['n_cracked']}")
    return snapshots


if __name__ == '__main__':
    import sys
    if len(sys.argv) > 1 and sys.argv[1] == '--real':
        from pathlib import Path
        data_dir = Path(__file__).parent / '3d_data'
        save_dir = Path(__file__).parent / 'results' / '3d_phase_field'
        sapa_path = data_dir / 'SAPA_cluster2_3d.tif'
        run_3d_phase_field_real(str(sapa_path), str(save_dir))
    else:
        demo_3d_phase_field()