| |
| """Publication figures for the hierarchical fracture study (PDF + 300 dpi PNG) and captions.json.""" |
| import os, json, glob, numpy as np, pandas as pd |
| import matplotlib; matplotlib.use('Agg'); import matplotlib.pyplot as plt |
| from matplotlib.lines import Line2D |
| from PIL import Image |
| from load_runs import load_all, ROOT |
| from render_utils import load_run, load_frames, draw_structure, curve_xy, geometry_arrays, member_polys, LEVEL_COLORS, frame_at_displacement |
| from matplotlib.collections import PolyCollection |
| FIG = os.path.join(ROOT, 'figures'); os.makedirs(FIG, exist_ok=True) |
| plt.rcParams.update({'font.size': 8, 'axes.labelsize': 8, 'axes.titlesize': 9, 'legend.fontsize': 6.5, 'xtick.labelsize': 7, 'ytick.labelsize': 7, 'figure.dpi': 110, 'axes.grid': True, 'grid.alpha': 0.25, 'axes.spines.top': False, 'axes.spines.right': False, 'pdf.fonttype': 42}) |
| CAP = {} |
| def save(fig, name, caption): |
| fig.savefig(os.path.join(FIG, name + '.pdf'), bbox_inches='tight'); fig.savefig(os.path.join(FIG, name + '.png'), dpi=300, bbox_inches='tight'); plt.close(fig); CAP[name] = caption; print('saved', name) |
| df = load_all() |
| BASE = (df.grip=='clamped')&(df.beam=='euler-bernoulli')&(df.alpha==1)&(df.beta==1)&(df.cascade=='sequential') |
| base = df[BASE] |
| def fam(d, **kw): |
| m = np.ones(len(d), bool) |
| for k, v in kw.items(): m &= (d[k].isin(v) if isinstance(v, (list, tuple, set)) else (d[k] == v)).values |
| return d[m] |
| R0F = dict(depth=4, branching='4x4x8x8', pattern='rect', defectType='none', profile='uniform', kappa=1, pruning=0, bracing=0) |
| def mean_sd(sub, col): return sub[col].mean(), (sub[col].std(ddof=1) if len(sub) > 1 else 0.0), len(sub) |
| PAL = {'flat': '#444444', 1.0: '#7f7f7f', 1.25: '#9ecae1', 1.35: '#6baed6', 1.45: '#4292c6', 1.6: '#2171b5', 1.8171: '#08519c', 2.0: '#e6550d', 2.5: '#a63603', 3.0: '#7a0177'} |
| def gcol(g): return PAL.get(round(float(g), 4), PAL.get(round(float(g), 2), '#2ca02c')) |
|
|
| |
| def fig1(): |
| fig = plt.figure(figsize=(7.2, 5.0)); gs = fig.add_gridspec(2, 5, height_ratios=[1, 2.0], hspace=0.22, wspace=0.08) |
| refs = [('ref4', '(a) Ref. 1: leaf venation'), ('ref2', '(b) Ref. 2: 2-level honeycomb'), ('ref3', '(c) Ref. 3: fibrous net'), ('ref5', '(d) Ref. 4: triangulated lattice'), ('ref1', '(e) Ref. 5: polar web')] |
| for i, (rid, lab) in enumerate(refs): |
| f = glob.glob(os.path.join(ROOT, 'reference_images', rid + '_*.jpg'))[0]; ax = fig.add_subplot(gs[0, i]); ax.imshow(Image.open(f)); ax.set_title(lab, fontsize=6); ax.axis('off') |
| run = load_run('S2_R0_s1'); ax = fig.add_subplot(gs[1, 0:2]); draw_structure(ax, run, color='level', grips=True); ax.set_title('(f) generated design R0: D = 4, γ = 1.6, η = 0.3\n(colour = hierarchy level, dark = coarsest)', fontsize=7) |
| |
| ax = fig.add_subplot(gs[1, 2:4]); nodes, n1, n2, lev, w = geometry_arrays(run); g = run['geometry']['members'] |
| polys = member_polys(nodes[n1], nodes[n2], w); cols = [LEVEL_COLORS[(l - 1) % 12] for l in lev]; ax.add_collection(PolyCollection(polys, facecolors=cols, edgecolors='none')) |
| seen = set() |
| for m in g: |
| for flag, n in ((m[4], m[0]), (m[5], m[1])): |
| if flag and n not in seen: seen.add(n); ax.add_patch(plt.Circle((nodes[n, 0], nodes[n, 1]), 0.35, fill=False, ec='#e5484d', lw=0.8)) |
| ax.set_xlim(18, 40); ax.set_ylim(18, 40); ax.set_aspect('equal'); ax.set_xticks([]); ax.set_yticks([]); ax.grid(False); ax.set_title('(g) zoom: interlevel (T-junction)\ninterfaces marked by red circles', fontsize=7) |
| ax = fig.add_subplot(gs[1, 4]); W, H = 40, 40; ax.add_patch(plt.Rectangle((0, 0), W, H, fc='#eef2f6', ec='#333', lw=0.8)); ax.add_patch(plt.Rectangle((-2, -4), W + 4, 4, fc='#777')); ax.add_patch(plt.Rectangle((-2, H), W + 4, 4, fc='#777')) |
| ax.annotate('', xy=(W / 2, H + 12), xytext=(W / 2, H + 4), arrowprops=dict(arrowstyle='-|>', lw=1.5, color='#c0392b')); ax.text(W / 2 + 2, H + 8, 'δ', color='#c0392b', fontsize=9) |
| ax.text(W / 2, -9, 'fixed grip', ha='center', fontsize=6.5); ax.text(W / 2, H / 2, 'W × H × t\n40 × 40 × 3 mm\nρ = 0.25 fixed\nE = 3 GPa\nσ_t = 40 MPa', ha='center', va='center', fontsize=6.5) |
| ax.set_xlim(-6, W + 6); ax.set_ylim(-12, H + 14); ax.set_aspect('equal'); ax.axis('off'); ax.set_title('(h) tensile test', fontsize=7) |
| save(fig, 'fig1_concept', 'Reference images and the adopted interpretation. (a-e) The five supplied images, read as one planar hierarchical reticulate-network family: nested cells whose walls are load-bearing members of decreasing width, fine members attaching to coarser ones at T-junctions (interlevel interfaces), closed loops (redundancy), pruned blind endings (Ref. 2), triangulation (Ref. 4), polar ordering and a defect (Ref. 5). (f) A generated specimen of the reference design R0 (depth D = 4, branching 4x4x8x8, width ratio γ = 1.6, geometric disorder η = 0.3, relative density ρ = 0.25); colour = hierarchy level (dark blue = coarsest). (g) Zoom of the upper-right quadrant with interlevel interfaces marked. (h) Displacement-controlled tension between a fixed and a moving clamped grip.') |
|
|
| |
| def fig2(): |
| d = df.copy(); fig, ax = plt.subplots(2, 3, figsize=(7.2, 4.6)); ax = ax.flat |
| st = {'stage1-baseline': ('o', '#7f7f7f'), 'stage2-reconnaissance': ('s', '#2171b5'), 'stage3-hypotheses': ('^', '#e6550d'), 'stage4-discriminating': ('D', '#31a354'), 'stage4b-mechanism': ('v', '#756bb1'), 'stage5-holdout': ('*', '#d62728'), 'stage6-representative': ('p', '#8c564b')} |
| for s, (mk, c) in st.items(): |
| sub = d[d.stage == s]; ax[0].scatter(sub.gamma + np.random.default_rng(1).normal(0, 0.012, len(sub)), sub.disorder + np.random.default_rng(2).normal(0, 0.004, len(sub)), marker=mk, c=c, s=14, alpha=.7, label=s.split('-')[0] + ' ' + s.split('-')[1][:6]) |
| ax[1].scatter(sub.depth + np.random.default_rng(3).normal(0, 0.05, len(sub)), sub.Gamma, marker=mk, c=c, s=14, alpha=.7); ax[2].scatter(sub.gamma, sub.kappa, marker=mk, c=c, s=14, alpha=.7) |
| ax[3].scatter(sub.gamma, sub.bracing - sub.pruning, marker=mk, c=c, s=14, alpha=.7) |
| ax[0].set_xlabel('width ratio γ'); ax[0].set_ylabel('geometric disorder η'); ax[0].set_title('(a) γ – η') |
| ax[1].set_xlabel('hierarchy depth D'); ax[1].set_ylabel('total width contrast Γ = w₁/w_D'); ax[1].set_yscale('log'); ax[1].set_title('(b) D – Γ'); ax[1].legend(*ax[0].get_legend_handles_labels(), ncol=1, loc='upper left', fontsize=5.5) |
| ax[2].set_xlabel('γ'); ax[2].set_ylabel('interface strength κ'); ax[2].set_yscale('log'); ax[2].set_title('(c) γ – κ') |
| ax[3].set_xlabel('γ'); ax[3].set_ylabel('bracing b (>0) / −pruning p (<0)'); ax[3].set_title('(d) connectivity controls') |
| cnt = d.groupby('stage').size(); ax[4].barh([s.replace('-', '\n') for s in cnt.index], cnt.values, color=[st[s][1] for s in cnt.index]); ax[4].set_xlabel('simulations'); ax[4].set_title('(e) runs per stage'); ax[4].tick_params(axis='y', labelsize=5.5) |
| cats = {'rect': (d.pattern == 'rect').sum(), 'hex': (d.pattern == 'hex').sum(), 'random': (d.pattern == 'random').sum(), 'polar': (d.pattern == 'polar').sum(), 'notched': (d.defectType == 'notch').sum(), 'fine-straight': (d.profile == 'fineStraight').sum(), 'coarse-straight': (d.profile == 'coarseStraight').sum(), 'β≠1 / α≠1': ((d.beta != 1) | (d.alpha != 1)).sum(), 'roller / Timoshenko / batch': ((d.grip == 'roller') | (d.beam == 'timoshenko') | (d.cascade == 'batch')).sum()} |
| ax[5].barh(list(cats.keys()), list(cats.values()), color='#6baed6'); ax[5].set_xlabel('simulations'); ax[5].set_title('(f) design languages & variants'); ax[5].tick_params(axis='y', labelsize=6) |
| fig.tight_layout(); n = len(d); nd = d.design.nunique() |
| save(fig, 'fig2_design_space', f'Coverage of the investigated design space: {n} simulations of {nd} unique designs (points jittered slightly for visibility). (a) Width ratio γ versus geometric disorder η; (b) hierarchy depth versus total width contrast Γ; (c) interlevel interface strength κ; (d) connectivity controls (bracing fraction positive, pruning fraction negative); (e) simulations per study stage; (f) design languages and model variants. All runs use W = H = 40 mm, t = 3 mm, ρ = 0.25 and 1024 finest cells except the three 4096-cell resolution checks.') |
|
|
| |
| def pareto_front(pts): |
| idx = np.argsort(-pts[:, 0]); front = []; best = -np.inf |
| for i in idx: |
| if pts[i, 1] > best: front.append(i); best = pts[i, 1] |
| return front |
| PRINT_SET = [('P01', 'S1_flat_rect_dis0.3', 'flat disordered grid (baseline; stiff, strong, abrupt)'), ('P02', 'S2_gamma1', 'hierarchical topology, uniform widths γ = 1'), ('P03', 'S2_R0', 'reference R0: γ = 1.6'), ('P04', 'S2_gamma2.5', 'γ = 2.5: two-tier progressive failure'), ('P05', 'S2_gamma3', 'γ = 3: toughest of the γ sweep'), ('P06', 'S4_fineStraight_g2.5', 'ordered fine level + γ = 2.5 (Pareto: stiff, strong, tough)'), ('P07', 'S4_coarseStraight_g2.5', 'straight skeleton + γ = 2.5 (strongest, abrupt)'), ('P08', 'S2_brace1', 'fully braced (triangulated): localized brittle crack'), ('P09', 'S4_random_g2.5', 'organic leaf-like Voronoi, γ = 2.5'), ('P10', 'S2_depth8_g1.6', 'deep hierarchy D = 8'), ('P11', 'S2_kappa0.5', 'weak interlevel interfaces κ = 0.5'), ('P12', 'S3_notch0.25_g2.5', 'notched γ = 2.5 (flaw tolerance)'), ('P13', 'S4_gamma0.7', 'inverted hierarchy γ = 0.7 (negative control)')] |
| def fig3(): |
| d = base[(base.defectType == 'none')]; d = d[d.disorder == 0.3] |
| g = d.groupby('design').agg(k0=('k0Norm', 'mean'), k0sd=('k0Norm', 'std'), Wf=('WfNorm', 'mean'), Wfsd=('WfNorm', 'std'), Fp=('FpeakNorm', 'mean'), Fpsd=('FpeakNorm', 'std'), n=('id', 'count'), gamma=('gamma', 'first'), pattern=('pattern', 'first'), profile=('profile', 'first'), depth=('depth', 'first'), bracing=('bracing', 'first'), pruning=('pruning', 'first'), kappa=('kappa', 'first')).reset_index() |
| fig, ax = plt.subplots(1, 3, figsize=(7.2, 2.7)) |
| for a, (xc, yc, xl) in zip(ax[:2], [('k0', 'Wf', 'normalised stiffness k₀H/(EWt)'), ('Fp', 'Wf', 'normalised strength F_peak/(σ_t W t)')]): |
| mk = np.where(g.pattern == 'rect', 'o', np.where(g.pattern == 'hex', 'h', np.where(g.pattern == 'random', 'P', 'X'))) |
| for i, r in g.iterrows(): |
| a.errorbar(r[xc], r.Wf, xerr=r[xc + 'sd'] if r.n > 1 else None, yerr=r.Wfsd if r.n > 1 else None, fmt=mk[i], ms=4 if r.n > 1 else 3, mfc=gcol(r.gamma) if r.depth > 1 else '#444', mec='k' if r.profile != 'uniform' else 'none', mew=0.6, ecolor='#bbb', elinewidth=0.5, capsize=0, alpha=0.9) |
| front = pareto_front(g[[xc, 'Wf']].values); fr = g.iloc[front].sort_values(xc); a.plot(fr[xc], fr.Wf, '-', color='#e5484d', lw=1, alpha=.8, zorder=0) |
| for _, r in fr.iterrows(): a.annotate(r.design.replace('S1_', '').replace('S2_', '').replace('S3_', '').replace('S4_', '').replace('S5_', '').replace('_g2.5', ' γ2.5').replace('flat_rect_dis0.3', 'flat').replace('coarseStraight', 'straight skeleton').replace('fineStraight', 'straight fine').replace('HO10', 'HO10 (straight fine, κ1.5)'), (r[xc], r.Wf), fontsize=5, xytext=(-4, 4), textcoords='offset points', ha='right') |
| a.set_xlabel(xl); a.set_ylabel('normalised work to failure W_f E/(σ_t² W H t)') |
| ax[0].set_title('(a) stiffness – toughness (η = 0.3)', fontsize=8); ax[1].set_title('(b) strength – toughness (η = 0.3)', fontsize=8) |
| |
| d2 = base[(base.defectType == 'none') & (base.profile == 'uniform') & (base.pattern == 'rect')] |
| g2 = d2.groupby(['design']).agg(k0=('k0Norm', 'mean'), Wf=('WfNorm', 'mean'), eta=('disorder', 'first'), gamma=('gamma', 'first'), depth=('depth', 'first')).reset_index(); g2 = g2[g2.eta > 0] |
| sc = ax[2].scatter(g2.k0, g2.Wf, c=g2.eta, cmap='magma_r', s=14, vmin=0.05, vmax=1.0, edgecolors=np.where(g2.depth == 1, 'k', 'none'), linewidths=0.5) |
| for eta in [0.1, 0.2, 0.3, 0.5]: |
| s = g2[g2.eta == eta]; fr = s.iloc[pareto_front(s[['k0', 'Wf']].values)].sort_values('k0'); ax[2].plot(fr.k0, fr.Wf, '-', color=plt.get_cmap('magma_r')((eta - 0.05) / 0.95), lw=0.8) |
| plt.colorbar(sc, ax=ax[2], label='disorder η'); ax[2].set_xlabel('normalised stiffness'); ax[2].set_ylabel('normalised work to failure'); ax[2].set_title('(c) rect family: one front per η\n(black edge = flat D = 1)', fontsize=8) |
| handles = [Line2D([], [], marker='o', ls='', mfc=gcol(g_), mec='none', label=f'γ = {g_}') for g_ in [1.0, 1.6, 2.0, 2.5, 3.0]] + [Line2D([], [], marker='o', ls='', mfc='#444', mec='none', label='flat (D = 1)'), Line2D([], [], marker='o', ls='', mfc='w', mec='k', label='ordered-level profile'), Line2D([], [], marker='h', ls='', mfc='#2ca02c', mec='none', label='hex'), Line2D([], [], marker='P', ls='', mfc='#2ca02c', mec='none', label='random'), Line2D([], [], color='#e5484d', label='Pareto front')] |
| ax[0].legend(handles=handles, ncol=2, fontsize=5, loc='upper left'); ax[0].set_ylim(top=ax[0].get_ylim()[1] * 1.25); fig.tight_layout() |
| save(fig, 'fig3_pareto', 'Stiffness–toughness and strength–toughness maps at fixed material (ρ = 0.25) and disorder η = 0.3 (a, b): each symbol is a design (mean over seeds, error bars = ±1 s.d. where n ≥ 3); fill colour = width ratio γ (grey = flat baseline), symbol = design language, black edge = designs with an ordered (straight) level; the red line joins the Pareto-optimal designs. (c) Rect family across disorder levels: colour = η, one Pareto front per η, black-edged points are flat (non-hierarchical) designs. Normalisations: k₀H/(EWt), F_peak/(σ_t W t) and W_f E/(σ_t² W H t).') |
| g.to_csv(os.path.join(ROOT, 'data', 'processed', 'pareto_designs_eta0.3.csv'), index=False) |
| return g |
|
|
| |
| def fig4(): |
| panels = [('(a) flat vs hierarchical (η = 0.3)', [('S1_flat_rect_dis0.3_s1', 'flat D = 1', '#444'), ('S2_gamma1', 'D = 4, γ = 1', gcol(1)), ('S2_R0_s1', 'D = 4, γ = 1.6', gcol(1.6)), ('S2_gamma2.5', 'D = 4, γ = 2.5', gcol(2.5))]), |
| ('(b) width ratio γ (D = 4)', [('S2_gamma1', 'γ = 1', gcol(1)), ('S2_gamma1.25', 'γ = 1.25', gcol(1.25)), ('S2_R0_s1', 'γ = 1.6', gcol(1.6)), ('S2_gamma2', 'γ = 2', gcol(2)), ('S2_gamma3', 'γ = 3', gcol(3)), ('S4_gamma0.7_s1', 'γ = 0.7 (inverted)', '#2ca02c')]), |
| ('(c) depth at fixed scales (Γ = 6)', [('S3_fixedscales_D2_G6_s1', 'D = 2', '#c6dbef'), ('S3_fixedscales_D4_G6_s1', 'D = 4', '#6baed6'), ('S3_fixedscales_D6_G6_s1', 'D = 6', '#2171b5'), ('S3_fixedscales_D9_G6_s1', 'D = 9', '#08306b')]), |
| ('(d) interface strength κ (γ = 1.6)', [('S2_kappa0.3', 'κ = 0.3', '#fdae6b'), ('S2_kappa0.5', 'κ = 0.5', '#f16913'), ('S2_R0_s1', 'κ = 1', gcol(1.6)), ('S2_kappa2', 'κ = 2', '#54278f')]), |
| ('(e) bracing / pruning (γ = 1.6)', [('S3_prune0.4_g1.6_s1', 'pruned 40 %', '#31a354'), ('S2_R0_s1', 'reference', gcol(1.6)), ('S2_brace0.5', 'braced 50 %', '#e6550d'), ('S2_brace1', 'braced 100 %', '#a50f15')]), |
| ('(f) disorder profile (γ = 2.5)', [('S2_gamma2.5', 'uniform η = 0.3', gcol(2.5)), ('S4_coarseStraight_g2.5_s1', 'straight skeleton', '#7a0177'), ('S4_fineStraight_g2.5_s1', 'straight fine level', '#238b45')])] |
| fig, ax = plt.subplots(2, 3, figsize=(7.2, 4.4)); ax = ax.flat |
| for a, (title, items) in zip(ax, panels): |
| for rid, lab, c in items: |
| r = load_run(rid); d, F = curve_xy(r); m = r['metrics']; a.plot(d, F, color=c, lw=0.9, label=f'{lab} (W_f = {m["workToFailure"]:.1f})') |
| a.set_title(title, fontsize=8); a.set_xlabel('grip displacement δ [mm]'); a.set_ylabel('reaction force F [N]'); a.legend(loc='center right') |
| fig.tight_layout(); save(fig, 'fig4_curves', 'Representative load–displacement responses (seed 1, physical units; δ/H = applied strain, F/(σ_t W t) = normalised force, 1 N corresponds to 2.1e-4 in normalised units). Vertical drops are rupture avalanches resolved exactly by the event-driven solver. (a) Flat disordered grid versus hierarchical designs of increasing width ratio; (b) γ sweep including the inverted hierarchy γ = 0.7; (c) hierarchy depth at fixed coarsest (2x2) and finest (1024 cells) scales and fixed total contrast Γ = 6; (d) interlevel interface strength; (e) connectivity via pruning and bracing; (f) disorder profiles at γ = 2.5. Work to failure W_f in N·mm in the legends.') |
|
|
| |
| def fig5(): |
| rows = [('S1_flat_rect_dis0.3_s1', 'flat grid, η = 0.3: single crack'), ('REP_brace1_g1.6', 'braced 100 %, γ = 1.6: localized crack'), ('S2_R0_s1', 'R0 (γ = 1.6): fine level first, skeleton later'), ('S2_gamma3', 'γ = 3: distributed progressive failure')] |
| fig, ax = plt.subplots(len(rows), 5, figsize=(7.2, 1.75 * len(rows)), gridspec_kw={'width_ratios': [1.1, 1.15, 1.15, 1.15, 0.7], 'wspace': 0.08, 'hspace': 0.45}) |
| for i, (rid, lab) in enumerate(rows): |
| r = load_run(rid); fr = load_frames(r); d, F = curve_xy(r); m = r['metrics'] |
| dpk = m['dPeak']; Fpk = m['Fpeak']; ipk = int(np.argmax(F)); below = np.where(F[ipk:] < 0.5 * Fpk)[0]; d_half = d[ipk + below[0]] if len(below) else d[-1] |
| ds = np.array([f_['d'] for f_ in fr['index']]); Fs = np.array([f_['F'] for f_ in fr['index']]) |
| k_peak = int(np.argmax(Fs)); after = np.where((ds >= d_half * 0.999) & (np.arange(fr['n']) > k_peak))[0]; k_half = int(after[0]) if len(after) else min(k_peak + 1, fr['n'] - 1) |
| if fr['index'][k_half]['nEvents'] == fr['index'][k_peak]['nEvents'] and k_half + 1 < fr['n']: k_half += 1 |
| picks = [k_peak, k_half, fr['n'] - 1] |
| a = ax[i, 0]; a.plot(d, F, color='#0b7a6e', lw=0.9); a.set_xlabel('δ [mm]', labelpad=1); a.set_ylabel('F [N]', labelpad=1); a.text(0.03, 0.97, lab, transform=a.transAxes, va='top', ha='left', fontsize=6, fontweight='bold', bbox=dict(fc='white', ec='none', alpha=0.85, pad=1.5)) |
| for k, c in zip(picks, ['#f2a541', '#e6550d', '#7a0177']): a.plot([fr['index'][k]['d']], [fr['index'][k]['F']], 'o', color=c, ms=4, mec='k', mew=0.4) |
| deform = float(np.clip(round(0.02 * r['H'] / max(dpk, 1e-6)), 1, 20)) |
| for j, (k, c) in enumerate(zip(picks, ['#f2a541', '#e6550d', '#7a0177'])): |
| draw_structure(ax[i, 1 + j], r, frame=k, frames=fr, color='r', deform=deform, grips=False); fi = fr['index'][k] |
| ax[i, 1 + j].set_title(f'δ={fi["d"]:.2f} mm F={fi["F"]:.0f} N\nfailed {fi["nEvents"]}/{r["E"]}', fontsize=5.5, color=c, pad=2); ax[i, 1 + j].set_ylim(-0.02 * r['H'], 1.02 * r['H'] + deform * fi['d']) |
| bands = m['damageBands']; a = ax[i, 4]; a.barh(np.arange(10) + 0.5, bands, height=0.9, color='#e5484d'); a.set_ylim(0, 10); a.set_yticks([0, 5, 10]); a.set_yticklabels(['0', 'H/2', 'H']); a.yaxis.tick_right(); a.set_xlabel('damage fraction', labelpad=1); a.set_title(f'Λ = {m["localization"]:.2f}', fontsize=7) |
| save(fig, 'fig5_snapshots', 'Fracture sequences contrasting localized brittle and distributed progressive failure. Left: load–displacement curve with the three snapshot states (at peak load, when the force first falls below half the peak, and at the final state). Middle: deformed geometry (displacements magnified 1–20x as needed) with intact members coloured by utilisation σ_max/σ_c (viridis, yellow = at rupture) and failed members drawn as red lines. Right: material-weighted fraction of failed members in ten horizontal bands with the damage-localisation index Λ (1 = single band). Rows: flat disordered grid; fully braced hierarchical lattice; reference R0; high width-contrast design.') |
|
|
| |
| def pm(ax, x, sub, col, color, label=None, marker='o', norm=1.0, ls='-'): |
| g = sub.groupby(x)[col]; m = g.mean() / norm; s = g.std(ddof=1).fillna(0) / norm |
| ax.errorbar(m.index, m.values, yerr=s.values, fmt=marker + ls, color=color, ms=3.5, lw=1, capsize=2, label=label) |
| ax.scatter(sub[x], sub[col] / norm, s=6, color=color, alpha=0.35, zorder=0) |
| def fig6(): |
| fig, ax = plt.subplots(2, 3, figsize=(7.2, 4.6)); ax = ax.flat |
| r0 = fam(base, **R0F); r0 = r0[r0.disorder == 0.3] |
| ref = r0[r0.gamma == 1] |
| for col, c, lab in [('k0', '#2171b5', 'stiffness k₀'), ('Fpeak', '#e6550d', 'peak force F_peak'), ('Wf', '#31a354', 'work to failure W_f')]: pm(ax[0], 'gamma', r0[r0.gamma >= 0.7], col, c, lab, norm=ref[col].mean()) |
| ax[0].axhline(1, color='#999', lw=0.6); ax[0].set_xlabel('width ratio γ'); ax[0].set_ylabel('value / value at γ = 1'); ax[0].set_title('(a) width contrast (D = 4, η = 0.3)'); ax[0].legend() |
| h5 = base[base.id.str.contains('fixedscales_D\\d+_G6')] |
| for col, c, lab in [('k0', '#2171b5', 'k₀ [N/mm]'), ('Fpeak', '#e6550d', 'F_peak [N] ×5'), ('Wf', '#31a354', 'W_f [N mm] ×40')]: pm(ax[1], 'depth', h5, col, c, lab, norm={'k0': 1, 'Fpeak': 0.2, 'Wf': 1 / 40}[col]) |
| ax[1].set_xlabel('hierarchy depth D (coarsest 2×2, finest 1024 cells, Γ = 6)'); ax[1].set_title('(b) depth per se'); ax[1].legend() |
| fbase = fam(base, pattern='rect', defectType='none', profile='uniform', kappa=1, pruning=0, bracing=0) |
| for g_, c in [(1.0, gcol(1)), (1.6, gcol(1.6)), (2.5, gcol(2.5))]: |
| sub = fbase[(fbase.depth == 4) & (fbase.branching == '4x4x8x8') & (fbase.gamma == g_) & (fbase.disorder.isin([0.1, 0.2, 0.3, 0.5]))]; pm(ax[2], 'disorder', sub, 'Wf', c, f'D = 4, γ = {g_}') |
| sub = fbase[(fbase.depth == 1) & (fbase.disorder.isin([0.1, 0.2, 0.3, 0.5]))]; pm(ax[2], 'disorder', sub, 'Wf', '#444', 'flat D = 1', marker='s', ls='--') |
| ax[2].set_xlabel('geometric disorder η'); ax[2].set_ylabel('W_f [N mm]'); ax[2].set_title('(c) disorder × width contrast'); ax[2].legend() |
| k = fam(base, depth=4, branching='4x4x8x8', pattern='rect', defectType='none', profile='uniform', pruning=0, bracing=0); k = k[k.disorder == 0.3] |
| for g_, c in [(1.6, gcol(1.6)), (2.5, gcol(2.5))]: |
| sub = k[k.gamma == g_]; pm(ax[3], 'kappa', sub, 'Fpeak', c, f'F_peak, γ = {g_}'); pm(ax[3], 'kappa', sub, 'Wf', c, f'W_f ×5, γ = {g_}', marker='s', norm=0.2, ls='--') |
| ax[3].set_xscale('log'); ax[3].set_xticks([0.3, 0.5, 1, 2, 3, 5]); ax[3].set_xticklabels(['0.3', '0.5', '1', '2', '3', '5']); ax[3].xaxis.set_minor_formatter(plt.NullFormatter()); ax[3].set_xlabel('interlevel interface strength κ'); ax[3].set_ylabel('F_peak [N], 5 W_f [N mm]'); ax[3].set_title('(d) interface strength'); ax[3].legend() |
| con = fam(base, depth=4, branching='4x4x8x8', pattern='rect', defectType='none', profile='uniform', kappa=1); con = con[(con.disorder == 0.3) & (con.gamma.isin([1.0, 1.6]))] |
| for g_, c in [(1.0, gcol(1)), (1.6, gcol(1.6))]: |
| sub = con[con.gamma == g_].copy(); sub['zr'] = sub.z.round(2) |
| br = sub[sub.pruning == 0]; pr_ = sub[sub.bracing == 0] |
| pm(ax[4], 'zr', br, 'Fpeak', c, f'F_peak, γ = {g_}'); pm(ax[4], 'zr', pr_, 'Fpeak', c, None); pm(ax[4], 'zr', br, 'localization', c, f'Λ ×100, γ = {g_}', marker='s', norm=0.01, ls='--'); pm(ax[4], 'zr', pr_, 'localization', c, None, marker='s', norm=0.01, ls='--') |
| ax[4].set_xlabel('mean coordination number z (pruning ← 3.1 → bracing)'); ax[4].set_ylabel('F_peak [N], 100 Λ'); ax[4].set_title('(e) connectivity'); ax[4].legend() |
| |
| n6 = base[(base.defectType == 'notch') | (base.defectType == 'none')]; n6 = n6[(n6.pattern == 'rect') & (n6.profile == 'uniform') & (n6.kappa == 1) & (n6.pruning == 0) & (n6.bracing == 0) & (n6.disorder == 0.3) & ((n6.depth == 1) | ((n6.depth == 4) & (n6.branching == '4x4x8x8') & (n6.gamma.isin([1.6, 2.5]))))] |
| n6 = n6.assign(famc=np.where(n6.depth == 1, 'flat', np.where(n6.gamma == 2.5, 'D = 4, γ = 2.5', 'D = 4, γ = 1.6'))) |
| refm = n6[n6.defectType == 'none'].groupby(['famc', 'seed']).Fpeak.mean(); rows = [] |
| for _, r in n6[n6.defectType == 'notch'].iterrows(): |
| if (r.famc, r.seed) in refm.index: rows.append({'famc': r.famc, 'a': r.defectSize, 'ratio': r.Fpeak / refm[(r.famc, r.seed)]}) |
| rr = pd.DataFrame(rows) |
| for f_, c in [('flat', '#444'), ('D = 4, γ = 1.6', gcol(1.6)), ('D = 4, γ = 2.5', gcol(2.5))]: pm(ax[5], 'a', rr[rr.famc == f_], 'ratio', c, f_) |
| ax[5].plot([0, 0.5], [1, 0.5], 'k--', lw=0.8, label='net section 1 − a/W'); ax[5].set_xlabel('edge-notch length a/W'); ax[5].set_ylabel('F_peak(notched) / F_peak(unnotched)'); ax[5].set_title('(f) notch sensitivity'); ax[5].legend() |
| fig.tight_layout(); save(fig, 'fig6_effects', 'Effects of the design parameters at fixed material (small dots = individual simulations, symbols = mean ± 1 s.d. over seeds; n = 3–8 per condition). (a) Width ratio γ at D = 4, η = 0.3, normalised by the γ = 1 means. (b) Hierarchy depth at fixed coarsest and finest scales and fixed total contrast Γ = 6. (c) Work to failure versus geometric disorder for the flat grid and hierarchical designs of increasing γ. (d) Peak force and work to failure versus interlevel interface strength κ (log axis). (e) Peak force and damage-localisation index versus mean coordination number z, varied by pruning (z < 3.1) and bracing (z > 3.1). (f) Notch sensitivity: peak force of edge-notched specimens relative to the unnotched specimen with the same seed; dashed line = net-section prediction.') |
| def fig6b(): |
| fig, ax = plt.subplots(1, 4, figsize=(7.2, 2.3)) |
| fbase = fam(base, depth=4, branching='4x4x8x8', pattern='rect', defectType='none', profile='uniform', kappa=1, pruning=0, bracing=0) |
| xs, ys, es = [], [], [] |
| for eta in [0.1, 0.2, 0.3, 0.5]: |
| a1 = fbase[(fbase.disorder == eta) & (fbase.gamma == 1)]; a2 = fbase[(fbase.disorder == eta) & (fbase.gamma == 2.5)] |
| if len(a1) and len(a2): xs.append(eta); ys.append(a2.k0.mean() / a1.k0.mean()); es.append(ys[-1] * np.sqrt((a1.k0.std(ddof=1) / a1.k0.mean() / np.sqrt(len(a1)))**2 + (a2.k0.std(ddof=1) / a2.k0.mean() / np.sqrt(len(a2)))**2)) |
| ax[0].errorbar(xs, ys, yerr=es, fmt='o-', color='#2171b5', label='β = 1 (physical)', capsize=2) |
| mech = df[df.id.str.startswith('S4b')] |
| for cond, c, lab in [('beta10', '#a63603', 'β = 10 (bending stiffened)'), ('beta0.1', '#31a354', 'β = 0.1 (bending softened)')]: |
| a1 = mech[(mech.id.str.contains(cond + '_g1_'))]; a2 = mech[(mech.id.str.contains(cond + '_g2.5_'))]; ax[0].plot([0.3], [a2.k0.mean() / a1.k0.mean()], 's', color=c, label=lab, ms=5) |
| e0 = df[df.id.isin(['S4b_eta0_g1', 'S4b_eta0_g2.5'])]; ax[0].plot([0.0], [e0[e0.gamma == 2.5].k0.iloc[0] / e0[e0.gamma == 1].k0.iloc[0]], 'D', color='#444', label='η = 0 (ordered)', ms=5) |
| ax[0].axhline(1, color='#999', lw=0.6); ax[0].set_xlabel('disorder η'); ax[0].set_ylabel('k₀(γ = 2.5) / k₀(γ = 1)'); ax[0].set_title('(a) stiffness gain of\nwidth contrast', fontsize=8); ax[0].set_ylim(0.8, 2.15); ax[0].legend(fontsize=5, loc='lower right') |
| prof = fam(base, depth=4, branching='4x4x8x8', pattern='rect', defectType='none', kappa=1, pruning=0, bracing=0); prof = prof[prof.disorder == 0.3] |
| labs = ['uniform', 'coarseStraight', 'fineStraight']; x = np.arange(3); wdt = 0.38 |
| for j, (g_, c) in enumerate([(1.0, gcol(1)), (2.5, gcol(2.5))]): |
| for k_, (col, a) in enumerate([('k0', ax[1]), ('Fpeak', ax[2])]): |
| m = [prof[(prof.gamma == g_) & (prof.profile == p)][col].mean() for p in labs]; s = [prof[(prof.gamma == g_) & (prof.profile == p)][col].std(ddof=1) for p in labs] |
| a.bar(x + (j - 0.5) * wdt, m, wdt, yerr=s, color=c, capsize=2, label=f'γ = {g_}'); a.set_xticks(x); a.set_xticklabels(['uniform η', 'straight\nskeleton', 'straight\nfine level'], fontsize=6) |
| ax[1].set_ylabel('k₀ [N/mm]'); ax[1].set_title('(b) straightened level:\nstiffness', fontsize=8); ax[1].legend(loc='upper left'); ax[2].set_ylabel('F_peak [N]'); ax[2].set_title('(c) straightened level:\nstrength', fontsize=8) |
| |
| pr = fam(base, depth=4, branching='4x4x8x8', pattern='rect', defectType='none', profile='uniform', kappa=1, bracing=0); pr = pr[pr.disorder == 0.3] |
| for g_, c in [(1.0, gcol(1)), (1.6, gcol(1.6)), (2.5, gcol(2.5))]: |
| sub = pr[pr.gamma == g_]; ref = sub[sub.pruning == 0] |
| for col, mk, ls, nm in [('Wf', 'o', '-', 'W_f'), ('Fpeak', 's', '--', 'F_peak')]: pm(ax[3], 'pruning', sub, col, c, f'{nm}, γ = {g_}', marker=mk, norm=ref[col].mean(), ls=ls) |
| ax[3].axhline(1, color='#999', lw=0.6); ax[3].set_xlabel('pruning fraction p'); ax[3].set_ylabel('value / unpruned'); ax[3].set_title('(d) pruning the finest level:\nW_f (solid), F_peak (dashed)', fontsize=8); ax[3].set_ylim(0.3, 1.8); ax[3].legend(fontsize=4.5, ncol=2, loc='upper left') |
| fig.tight_layout(); save(fig, 'fig6b_mechanism', 'Mechanism tests. (a) Stiffness gain of width contrast (γ = 2.5 over γ = 1, same topology and material) versus disorder for the physical bending stiffness (β = 1), for numerically stiffened (β = 10) and softened (β = 0.1) bending, and for the perfectly ordered lattice (η = 0): the gain vanishes without bending compliance. (b, c) Stiffness and peak force when the geometric disorder η = 0.3 is applied uniformly to all levels, only to the fine levels (straight skeleton) or only to the coarse levels (straight fine level), at γ = 1 and γ = 2.5 (mean ± s.d., n = 3–8). (d) Retention of work to failure and peak force under pruning of the finest level.') |
|
|
| |
| def fig7(): |
| hs = json.load(open(os.path.join(ROOT, 'data', 'processed', 'holdout_summary.json'))); rows = pd.DataFrame(hs['rows']) |
| fig, ax = plt.subplots(2, 3, figsize=(7.2, 4.6)); ax = ax.flat |
| for a, (k, lab) in zip(ax[:3], [('k0', 'initial stiffness k₀ [N/mm]'), ('Fpeak', 'peak force F_peak [N]'), ('Wf', 'work to failure W_f [N mm]')]): |
| a.errorbar(rows[f'{k}_model'], rows[f'{k}_obs'], yerr=rows[f'{k}_sd'], fmt='o', color='#2171b5', ms=4, capsize=2, label='regression model') |
| a.plot(rows[f'{k}_expert'], rows[f'{k}_obs'], '^', color='#e6550d', ms=4, mfc='none', label='mechanistic estimate') |
| lo, hi = min(rows[f'{k}_obs'].min(), rows[f'{k}_model'].min()) * 0.7, max(rows[f'{k}_obs'].max(), rows[f'{k}_model'].max()) * 1.4; xx = np.array([lo, hi]) |
| sd = json.load(open(os.path.join(ROOT, 'data', 'holdout_predictions.json')))['models'][k]['residual_sd'] |
| a.plot(xx, xx, 'k-', lw=0.8); a.plot(xx, xx * np.exp(2 * sd), 'k--', lw=0.5); a.plot(xx, xx * np.exp(-2 * sd), 'k--', lw=0.5); a.set_xscale('log'); a.set_yscale('log'); a.set_xlim(lo, hi); a.set_ylim(lo, hi) |
| for _, r in rows.iterrows(): a.annotate(r.id.replace('HO', ''), (r[f'{k}_model'], r[f'{k}_obs']), fontsize=5, xytext=(2, 2), textcoords='offset points') |
| a.set_xlabel('predicted'); a.set_ylabel('observed (mean ± s.d., 3 seeds)'); a.set_title(lab, fontsize=8) |
| a.xaxis.set_major_locator(plt.LogLocator(base=10, numticks=4)); a.xaxis.set_minor_formatter(plt.NullFormatter()); a.yaxis.set_minor_formatter(plt.NullFormatter()) |
| ax[0].legend(loc='lower right') |
| for a, (k, lab, lim) in zip(ax[3:5], [('ppw', 'post-peak work fraction', (0, 1)), ('localization', 'damage-localisation index Λ', (0, 0.6))]): |
| a.errorbar(rows[f'{k}_model'], rows[f'{k}_obs'], yerr=rows[f'{k}_sd'], fmt='o', color='#2171b5', ms=4, capsize=2); a.plot(rows[f'{k}_expert'], rows[f'{k}_obs'], '^', color='#e6550d', ms=4, mfc='none'); a.plot(lim, lim, 'k-', lw=0.8); a.set_xlim(lim); a.set_ylim(lim) |
| for _, r in rows.iterrows(): a.annotate(r.id.replace('HO', ''), (r[f'{k}_model'], r[f'{k}_obs']), fontsize=5, xytext=(2, 2), textcoords='offset points') |
| a.set_xlabel('predicted'); a.set_ylabel('observed'); a.set_title(lab, fontsize=8) |
| a = ax[5]; a.axis('off'); s = hs['summary'] |
| txt = ['holdout summary (12 designs × 3 seeds)', '', 'median |error| model / mechanistic:', f" k₀ {s['k0']['model_median_abs_pct_error']:.0f} % / {s['k0']['expert_median_abs_pct_error']:.0f} %", f" F_peak {s['Fpeak']['model_median_abs_pct_error']:.0f} % / {s['Fpeak']['expert_median_abs_pct_error']:.0f} %", f" W_f {s['Wf']['model_median_abs_pct_error']:.0f} % / {s['Wf']['expert_median_abs_pct_error']:.0f} %", '', f"95 % band coverage (model): k₀ {100*s['k0']['model_coverage_95']:.0f} %, F_peak {100*s['Fpeak']['model_coverage_95']:.0f} %, W_f {100*s['Wf']['model_coverage_95']:.0f} %", f"mode agreement: model {100*s['mode']['model_agreement']:.0f} %, mechanistic {100*s['mode']['expert_agreement']:.0f} %", '', 'labels = holdout design number'] |
| a.text(0, 1, '\n'.join(txt), va='top', fontsize=6.5, family='monospace') |
| fig.tight_layout(); save(fig, 'fig7_holdout', 'Holdout test: predictions recorded before running twelve unseen designs (three seeds each) versus observations. Circles: log-linear regression fitted to stages 1–4 (dashed lines: ±2 residual s.d. of the fit); triangles: hand-written mechanistic estimates. Panels: initial stiffness, peak force, work to failure (log axes), post-peak work fraction, damage-localisation index, and a summary of median absolute errors, interval coverage and qualitative fracture-mode agreement.') |
|
|
| |
| def fig8(): |
| v = json.load(open(os.path.join(ROOT, 'data', 'validation', 'validation_results.json')))['results'] |
| fig, ax = plt.subplots(figsize=(7.2, 3.4)); ids = [r['id'] for r in v]; err = np.array([max(r['error'], 1e-17) if isinstance(r['error'], (int, float)) else 1 for r in v]); tol = np.array([max(r['tolerance'], 1e-17) for r in v]) |
| cols = ['#31a354' if r['pass'] else '#e5484d' for r in v]; y = np.arange(len(v)) |
| ax.barh(y, err, color=cols, height=0.7); ax.plot(tol, y, 'k|', ms=8, label='tolerance'); ax.set_xscale('log'); ax.set_yticks(y); ax.set_yticklabels([f"{r['id']} {r['name'][:58]}" for r in v], fontsize=5.5); ax.invert_yaxis(); ax.set_xlabel('error (relative unless noted; 1e-17 shown for exactly zero)'); ax.legend(loc='upper right'); ax.set_xlim(1e-17, 10) |
| fig.tight_layout(); save(fig, 'fig8_validation', 'Validation suite (31 tests): observed error (bars, log scale) against the acceptance tolerance (ticks) for element benchmarks against analytical beam theory, structural benchmarks, reproducibility, scaling laws, symmetry, equilibrium residuals, energy conservation, step-size sensitivity, export consistency and generator checks. All tests pass; zero errors are plotted at 1e-17.') |
|
|
| if __name__ == '__main__': |
| fig1(); fig2(); g = fig3(); fig4(); fig5(); fig6(); fig6b(); fig7(); fig8() |
| json.dump(CAP, open(os.path.join(FIG, 'captions.json'), 'w'), indent=1); json.dump(PRINT_SET, open(os.path.join(ROOT, 'data', 'processed', 'print_set.json'), 'w'), indent=1) |
| print('all figures done') |
|
|
