Many builders, one pipeline: visualizations

import base64
from pathlib import Path

import numpy as np
import pandas as pd
import plotly.graph_objects as go
from plotly.subplots import make_subplots

summary = pd.read_csv(Path("results") / "summary.csv")
MODELS = {  # builder: label used in the figures
    "abaqus_honeycomb": "Honeycomb",
    "gmshmodel_udfrp": "UD FRP",
    "microgen_tpms": "Schwarz-P TPMS",
    "texgen_weave": "Plain weave",
}


def stiffness(row):
    """Symmetric 6x6 Voigt stiffness (11, 22, 33, 23, 13, 12) from the C11..C66 columns."""
    C = np.zeros((6, 6))
    for i in range(6):
        for j in range(i, 6):
            C[i, j] = C[j, i] = row[f"C{i + 1}{j + 1}"]
    return C

MODEL_COLORS = dict(zip(MODELS, ["#4C78A8", "#F58518", "#54A24B", "#B279A2"]))
SYMBOLS = {"set1": "circle", "set2": "diamond-open"}
fig = go.Figure()
for builder, label in MODELS.items():
    for material in ("set1", "set2"):
        rows = summary[(summary.builder == builder) & (summary.material == material)].sort_values("density")
        fig.add_scatter(
            x=rows.density, y=rows.A_energy_ratio, mode="lines+markers",
            name=f"{label}, material set {material[-1]}", legendgroup=builder,
            line=dict(color=MODEL_COLORS[builder], dash="solid" if material == "set1" else "dot"),
            marker=dict(symbol=SYMBOLS[material], size=np.where(rows.baseline == 1, 13, 8),
                        line=dict(width=2, color=MODEL_COLORS[builder])),
            customdata=rows[["case"]], hovertemplate="%{customdata[0]}<br>ρ = %{x:.0f} kg/m³<br>A = %{y:.3g}<extra></extra>")
fig.update_layout(xaxis=dict(type="log", title="effective density (kg/m³)"),
                  yaxis=dict(type="log", title="A<sub>energy ratio</sub>"),
                  legend=dict(orientation="h", yanchor="bottom", y=1.02, xanchor="left", x=0),
                  template="plotly_white", height=600, margin=dict(t=20))
fig.show()

def directional_young(C, n_theta=40, n_phi=80):
    """Points r = E(n) n on the unit sphere, E(n) = 1 / (n n : S : n n)."""
    S_voigt = np.linalg.inv(C)
    pairs = [(0, 0), (1, 1), (2, 2), (1, 2), (0, 2), (0, 1)]
    S = np.zeros((3, 3, 3, 3))
    for I, (i, j) in enumerate(pairs):
        for J, (k, l) in enumerate(pairs):
            value = S_voigt[I, J] / ((1 if I < 3 else 2) * (1 if J < 3 else 2))
            for a, b in {(i, j), (j, i)}:
                for c, d in {(k, l), (l, k)}:
                    S[a, b, c, d] = value
    theta, phi = np.meshgrid(np.linspace(0, np.pi, n_theta), np.linspace(0, 2 * np.pi, n_phi))
    n = np.stack([np.sin(theta) * np.cos(phi), np.sin(theta) * np.sin(phi), np.cos(theta)])
    E = 1 / np.einsum("i...,j...,k...,l...,ijkl->...", n, n, n, n, S)
    return n * E, E


# axis triad with labels at the tips (plotly's axis titles mark box edges, not directions)
TRIAD = go.Scatter3d(x=[0, 1.2, None, 0, 0, None, 0, 0], y=[0, 0, None, 0, 1.2, None, 0, 0],
                     z=[0, 0, None, 0, 0, None, 0, 1.2], text=["", "x", "", "", "y", "", "", "z"],
                     mode="lines+text", line=dict(color="black", width=3), textfont=dict(size=14),
                     hoverinfo="skip", showlegend=False)

representative = summary[summary.baseline == 1]
fig = make_subplots(rows=4, cols=3, specs=[[{"type": "xy"}, {"type": "scene"}, {"type": "scene"}]] * 4,
                    column_titles=["SG", "Material set 1", "Material set 2"], row_titles=list(MODELS.values()),
                    column_widths=[0.3, 0.35, 0.35], horizontal_spacing=0.01, vertical_spacing=0.02)
for row, builder in enumerate(MODELS, start=1):
    for col, material in enumerate(("set1", "set2"), start=2):
        record = representative[(representative.builder == builder) & (representative.material == material)].iloc[0]
        (x, y, z), E = directional_young(stiffness(record))
        scale = E.max()
        fig.add_surface(x=x / scale, y=y / scale, z=z / scale, surfacecolor=E / scale,
                        cmin=0, cmax=1, colorscale="Viridis", showscale=row == 1 and col == 2,
                        colorbar=dict(title=dict(text="E / E<sub>max</sub>", side="top"), orientation="h",
                                      x=0.65, xanchor="center", y=-0.02, yanchor="top", len=0.5, thickness=15),
                        hovertemplate=f"E = %{{customdata:.3g}} GPa<extra>{record.case}</extra>",
                        customdata=E / 1e9, row=row, col=col)
        fig.add_trace(TRIAD, row=row, col=col)
    # first column: the SG image of material set 1 (same as in the Preprocess figure)
    image = Path("images") / f"{representative[(representative.builder == builder) & (representative.material == 'set1')].case.iloc[0]}.png"
    source = "data:image/png;base64," + base64.b64encode(image.read_bytes()).decode()
    fig.add_trace(go.Image(source=source, hoverinfo="skip"), row=row, col=1)
    fig.update_xaxes(visible=False, row=row, col=1)
    fig.update_yaxes(visible=False, row=row, col=1)
axis = dict(range=[-1.3, 1.3], showticklabels=False, title="")
scene = dict(xaxis=axis, yaxis=axis, zaxis=axis, aspectmode="cube")
fig.update_scenes(scene)
fig.update_layout(height=1250, margin=dict(l=0, r=30, t=40, b=80), template="plotly_white")
fig.show()