Fibre orientation and strain-field frames (with plots)
A carbon/epoxy cantilever plate on a triangle mesh – so the per-element in-plane
frame alternates cell to cell – with a curvilinear fibre path: the fibres fan from
0 deg at the clamped root to THETA_TIP at the free tip, set through a per-cell
fiber_angle field. The script then
reads the bending
curvatureback as aFieldand re-expresses the same component (kappa_xx) withto_frame– the raw element frame (a per-cell artefact -> patchy), the local fibre frame, and the global Cartesian frame (to_global()-> the smooth physical field),contrasts that DG (per-cell, discontinuous) field with a smooth CG-projected nodal field (
hm.strain_fields(..., space=("Lagrange", 1))-> a global-frameField),sweeps a uniform fibre angle and plots compliance and the Tsai-Wu failure index.
Writes ex_orientation_fields.png next to this file (needs matplotlib).
conda activate hermit
python examples/basic_examples/ex_orientation_fields.py
"""Fibre orientation and strain-field frames (with plots)
A carbon/epoxy cantilever plate on a **triangle** mesh -- so the per-element in-plane
frame alternates cell to cell -- with a *curvilinear* fibre path: the fibres fan from
0 deg at the clamped root to ``THETA_TIP`` at the free tip, set through a per-cell
``fiber_angle`` field. The script then
* reads the bending ``curvature`` back as a ``Field`` and re-expresses the same
component (``kappa_xx``) with ``to_frame`` -- the raw element frame (a per-cell
artefact -> patchy), the local fibre frame, and the global Cartesian frame
(``to_global()`` -> the smooth physical field),
* contrasts that DG (per-cell, discontinuous) field with a smooth CG-projected nodal
field (``hm.strain_fields(..., space=("Lagrange", 1))`` -> a global-frame ``Field``),
* sweeps a *uniform* fibre angle and plots compliance and the Tsai-Wu failure index.
Writes ``ex_orientation_fields.png`` next to this file (needs ``matplotlib``).
conda activate hermit
python examples/basic_examples/ex_orientation_fields.py
"""
import pathlib
import numpy as np
import csdl_alpha as csdl
import hermit as hm
from hermit import Layup
from caddee_materials import TransverseMaterial
MESH = pathlib.Path(__file__).parents[2] / "tests" / "meshes" / "plate_2x10_tri_4x20.xdmf"
OUT_PNG = pathlib.Path(__file__).with_name("ex_orientation_fields.png")
LENGTH, WIDTH = 10.0, 2.0
THETA_TIP = 60.0 # curvilinear fibre angle at the tip (deg)
TOTAL_H, DENSITY = 0.02, 1.6e3
PRESSURE_Z = 1.0e3
def clamped_at_x0(x):
return np.less(x[0], 1e-12)
def main():
mesh = hm.read_mesh(MESH)
rec = csdl.Recorder(inline=True)
rec.start()
domain = hm.ShellDomain(mesh, element="CG2CG1")
ud = TransverseMaterial(name="ud", EA=138e9, ET=10e9, vA=0.34, vT=0.67, GA=7e9, density=DENSITY)
ud.set_strength(F1t=1500e6, F1c=1200e6, F2t=50e6, F2c=200e6, F12=70e6, F23=40e6)
layup = Layup(ud, csdl.Variable(value=np.zeros(3)), np.full(3, TOTAL_H / 3), num_plies=3)
loads = hm.pressure(domain, PRESSURE_Z)
bcs = hm.clamp(domain, where=clamped_at_x0)
# -- per-cell fibre frame (file cell order) ----------------------------
frames = domain.local_frames() # (n_cells, 3, 3), file order
mids = domain.cell_centroids
theta = np.radians(THETA_TIP) * mids[:, 0] / LENGTH # 0 at root -> THETA_TIP at tip
fibre_dir = (np.cos(theta)[:, None] * frames[:, 0]
+ np.sin(theta)[:, None] * frames[:, 1]) # (n_cells, 3), file order
material = hm.laminate(domain, layup=layup, density=DENSITY,
orientation=hm.fiber_angle(domain, csdl.Variable(value=theta)))
# -- solve once, post-process twice (DG default, then smooth CG) --------
state = hm.solve(domain, material, loads, bcs)
compliance, failure_index = hm.compliance(state), hm.failure_index(state)
_, kap_dg, _ = hm.strain_fields(state) # ("DG", 2), element frame
_, kap_cg, _ = hm.strain_fields(state, space=("Lagrange", 1)) # smooth nodal, global frame
print(f" triangle mesh: {domain.n_cells} cells; per-cell fibre 0 deg (root) -> {THETA_TIP:.0f} deg (tip)")
print(f" compliance : {float(np.ravel(compliance.value)[0]):.4e}")
print(f" Tsai-Wu failure index : {float(np.ravel(failure_index.value)[0]):.4e}")
print(f" curvature DG field (element frame) : {kap_dg.values.shape} kind={kap_dg.kind}")
print(f" curvature CG field (global frame) : "
f"{np.asarray(kap_cg.coeffs.value).reshape(-1, 6).shape} kind={kap_cg.kind}")
# kappa_xx in three frames (DG, per-cell, user cell order)
kxx_elem = kap_dg.values[:, 0]
kxx_fibre = kap_dg.to_frame(fibre_dir).values[:, 0] # resolved along the local fibre
kxx_glob = kap_dg.to_global().values[:, 0] # smooth physical field
kxx_cg = np.asarray(kap_cg.coeffs.value).reshape(-1, 6)[:, 0] # global frame, CG1 nodes
# -- uniform fibre-angle sweep ----------------------------------------
sweep_deg = np.arange(0, 91, 15)
comp, fi = [], []
for deg in sweep_deg:
d = [np.cos(np.radians(deg)), np.sin(np.radians(deg)), 0.0]
m = hm.laminate(domain, layup=layup, density=DENSITY,
orientation=hm.fiber_direction(domain, d))
s = hm.solve(domain, m, loads, bcs)
comp.append(float(np.ravel(hm.compliance(s).value)[0]))
fi.append(float(np.ravel(hm.failure_index(s).value)[0]))
rec.stop()
# -- figure ----------------------------------------------------------
try:
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from matplotlib.colors import TwoSlopeNorm
from matplotlib.tri import Triangulation
except ImportError:
print("\n matplotlib not available -- skipping the figure")
return
verts = mesh.geometry.x[:, :2]
try:
cell_nodes = np.asarray(mesh.geometry.dofmaps[0])
except (AttributeError, IndexError):
cell_nodes = np.asarray(mesh.geometry.dofmap)
tri_cells = Triangulation(verts[:, 0], verts[:, 1], cell_nodes[domain.reverse_cell_idx])
cg_xy = domain.dof_coords(kap_cg.space)[:, :2]
tri_nodes = Triangulation(cg_xy[:, 0], cg_xy[:, 1])
allk = np.concatenate([kxx_elem, kxx_fibre, kxx_glob, kxx_cg])
norm = TwoSlopeNorm(vcenter=0.0, vmin=min(allk.min(), -1e-9), vmax=max(allk.max(), 1e-9))
fig, axes = plt.subplots(2, 3, figsize=(15, 6), constrained_layout=True)
kw = dict(cmap="RdBu_r", norm=norm)
kpanels = [
(axes[0, 1], kxx_elem, False, r"$\kappa_{xx}$ -- raw element frame (DG2, per-cell)"),
(axes[0, 2], kxx_fibre, False, r"$\kappa$ along local fibre -- to_frame(fibre)"),
(axes[1, 0], kxx_glob, False, r"$\kappa_{xx}$ -- global frame -- to_global()"),
(axes[1, 1], kxx_cg, True, r"$\kappa_{xx}$ -- global frame, smooth (CG1 project)"),
]
for ax, vals, gouraud, title in kpanels:
m = (ax.tripcolor(tri_nodes, vals, shading="gouraud", **kw) if gouraud
else ax.tripcolor(tri_cells, facecolors=vals, **kw))
ax.set_xlim(0, LENGTH); ax.set_ylim(0, WIDTH)
ax.set_title(title, fontsize=10)
fig.colorbar(m, ax=ax, fraction=0.03, pad=0.02)
ax = axes[0, 0]
ax.triplot(tri_cells, color="0.88", lw=0.3)
gx, gy = np.meshgrid(np.linspace(0.5, LENGTH - 0.5, 20), np.linspace(0.3, WIDTH - 0.3, 4))
gth = np.radians(THETA_TIP) * gx / LENGTH
q = ax.quiver(gx, gy, np.cos(gth), np.sin(gth), np.degrees(gth),
cmap="viridis", pivot="mid", scale=22, width=0.006)
ax.set_xlim(0, LENGTH); ax.set_ylim(0, WIDTH)
ax.set_title(f"curvilinear fibre_angle (0 deg root -> {THETA_TIP:.0f} deg tip)", fontsize=10)
fig.colorbar(q, ax=ax, fraction=0.03, pad=0.02, label="deg")
ax = axes[1, 2]
ax.plot(sweep_deg, comp, "o-", color="C0")
ax.set_xlabel("uniform fibre angle from beam axis (deg)")
ax.set_ylabel("compliance", color="C0"); ax.tick_params(axis="y", labelcolor="C0")
ax2 = ax.twinx()
ax2.plot(sweep_deg, fi, "s--", color="C3")
ax2.set_ylabel("Tsai-Wu failure index", color="C3"); ax2.tick_params(axis="y", labelcolor="C3")
ax.set_title("uniform fibre-angle sweep", fontsize=10)
fig.savefig(OUT_PNG, dpi=130)
rel = OUT_PNG.relative_to(pathlib.Path.cwd()) if OUT_PNG.is_relative_to(pathlib.Path.cwd()) else OUT_PNG
print(f"\n wrote {rel}")
if __name__ == "__main__":
main()