{ "cells": [ { "cell_type": "markdown", "id": "3bc51e72", "metadata": {}, "source": [ "# Shell analysis walkthrough\n", "\n", "This tutorial builds one shell model end to end: a mesh, a domain, boundary\n", "conditions, a material, loads, the solve, and the outputs — then makes the thickness a\n", "design variable and takes a total derivative.\n", "\n", "It assumes the `hermit` conda environment (DOLFINx 0.9 or 0.11) is active. Every\n", "number Hermit returns is a `csdl.Variable`, so the whole thing has to run inside a\n", "CSDL recorder." ] }, { "cell_type": "code", "execution_count": 1, "id": "d39ef462", "metadata": { "execution": { "iopub.execute_input": "2026-09-09T21:36:39.486064Z", "iopub.status.busy": "2026-09-09T21:36:39.485975Z", "iopub.status.idle": "2026-09-09T21:36:40.148951Z", "shell.execute_reply": "2026-09-09T21:36:40.148490Z" } }, "outputs": [], "source": [ "import numpy as np\n", "import csdl_alpha as csdl\n", "\n", "import hermit as hm\n", "\n", "recorder = csdl.Recorder(inline=True)\n", "recorder.start()" ] }, { "cell_type": "markdown", "id": "9667f1cf", "metadata": {}, "source": [ "## 1. A mesh\n", "\n", "A 2 × 10 cantilever plate of quadrilaterals. Hermit takes an ordinary DOLFINx mesh, so\n", "in a real workflow this comes from your geometry pipeline — `hm.read_mesh(\"plate.xdmf\")`\n", "reads one from XDMF. Here we build it in memory to keep the tutorial self-contained.\n", "\n", "Two traps are worth knowing about, and both are silent:\n", "\n", "- `create_mesh` **swaps its last two positional arguments** between DOLFINx 0.9 and\n", " 0.11, so dispatch on the signature rather than a version string.\n", "- basix orders quadrilateral vertices by tensor product — `(0,0), (1,0), (0,1), (1,1)`\n", " — **not** cyclically around the cell. Hand `create_mesh` a cyclic quad and you get\n", " bowtie cells, a solve that runs happily, and answers off by orders of magnitude." ] }, { "cell_type": "code", "execution_count": 2, "id": "5226633c", "metadata": { "execution": { "iopub.execute_input": "2026-09-09T21:36:40.150341Z", "iopub.status.busy": "2026-09-09T21:36:40.150164Z", "iopub.status.idle": "2026-09-09T21:36:40.154940Z", "shell.execute_reply": "2026-09-09T21:36:40.154606Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "105 vertices, 80 cells\n" ] } ], "source": [ "import inspect\n", "\n", "import basix.ufl\n", "import ufl\n", "from mpi4py import MPI\n", "from dolfinx.mesh import create_mesh\n", "\n", "LENGTH, WIDTH = 10.0, 2.0\n", "nx, ny = 20, 4\n", "\n", "xs = np.linspace(0.0, LENGTH, nx + 1)\n", "ys = np.linspace(0.0, WIDTH, ny + 1)\n", "points = np.array([[x, y, 0.0] for y in ys for x in xs])\n", "# tensor-product vertex order per cell, not cyclic\n", "cells = np.array([[j * (nx + 1) + i, j * (nx + 1) + i + 1,\n", " (j + 1) * (nx + 1) + i, (j + 1) * (nx + 1) + i + 1]\n", " for j in range(ny) for i in range(nx)], dtype=np.int64)\n", "\n", "el = ufl.Mesh(basix.ufl.element(\"Lagrange\", \"quadrilateral\", 1, shape=(3,)))\n", "if list(inspect.signature(create_mesh).parameters)[2] == \"e\": # DOLFINx 0.11\n", " mesh = create_mesh(MPI.COMM_WORLD, cells, el, points)\n", "else: # DOLFINx 0.9\n", " mesh = create_mesh(MPI.COMM_WORLD, cells, points, el)\n", "\n", "print(f\"{points.shape[0]} vertices, {cells.shape[0]} cells\")" ] }, { "cell_type": "markdown", "id": "4d3ce5e7", "metadata": {}, "source": [ "## 2. The domain\n", "\n", "`ShellDomain` is the one-time setup: the mixed state space, the mesh tags, and the\n", "maps between mesh-file ordering and FE dof ordering. It carries no CSDL and no\n", "per-solve cost, so build it **once** and reuse it across every load case, sweep and\n", "optimizer iteration.\n", "\n", "The default element is `CG2CG1` — quadratic displacement, linear rotation. `CG1CG1` is\n", "cheaper; `CG2CR1` uses a Crouzeix–Raviart rotation that resists shear locking on thin\n", "shells, and is triangle-only.\n", "\n", "The mesh must be serial: `ShellDomain` requires its vertex numbering to be a\n", "permutation of the file order, and raises on an MPI-partitioned mesh." ] }, { "cell_type": "code", "execution_count": 3, "id": "7232c5e6", "metadata": { "execution": { "iopub.execute_input": "2026-09-09T21:36:40.156068Z", "iopub.status.busy": "2026-09-09T21:36:40.155942Z", "iopub.status.idle": "2026-09-09T21:36:40.211813Z", "shell.execute_reply": "2026-09-09T21:36:40.211337Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "105 nodes, 80 cells, quadrature degree 4\n" ] } ], "source": [ "domain = hm.ShellDomain(mesh, element=\"CG2CG1\")\n", "print(f\"{domain.n_nodes} nodes, {domain.n_cells} cells, quadrature degree {domain.quadrature_degree}\")" ] }, { "cell_type": "markdown", "id": "c9c47817", "metadata": {}, "source": [ "## 3. Boundary conditions\n", "\n", "A BC is a *region predicate* plus a dof mask. The predicate is the raw DOLFINx entity\n", "locator — it receives a `(3, N)` coordinate array and returns an `(N,)` boolean mask —\n", "and `hm.near` / `hm.on_plane` build the common ones.\n", "\n", "`hm.clamp` fixes all six generalized dofs; `hm.pin` takes a named subset of\n", "`(\"ux\", \"uy\", \"uz\", \"rx\", \"ry\", \"rz\")`; `hm.symmetry` applies a plane mask; `hm.gauge`\n", "pins a single point to remove leftover rigid-body modes. All of them compose with `+`,\n", "and each takes `value=` for a non-zero prescribed displacement.\n", "\n", "Like the domain, boundary conditions are compiled once and reused." ] }, { "cell_type": "code", "execution_count": 4, "id": "4bec68fc", "metadata": { "execution": { "iopub.execute_input": "2026-09-09T21:36:40.212998Z", "iopub.status.busy": "2026-09-09T21:36:40.212894Z", "iopub.status.idle": "2026-09-09T21:36:40.215747Z", "shell.execute_reply": "2026-09-09T21:36:40.215427Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1 penalty region(s), 0 strong constraint(s)\n" ] } ], "source": [ "bcs = hm.clamp(domain, where=hm.near(\"x\", 0.0))\n", "print(f\"{len(bcs.penalty_terms)} penalty region(s), {len(bcs.strong)} strong constraint(s)\")" ] }, { "cell_type": "markdown", "id": "1d0002ec", "metadata": {}, "source": [ "## 4. A material\n", "\n", "Every material input may be a scalar (broadcast to the mesh), a per-vertex or per-cell\n", "array, or a `Field` on its own FE space — they need not agree with each other. Here\n", "they are all constants." ] }, { "cell_type": "code", "execution_count": 5, "id": "0952b285", "metadata": { "execution": { "iopub.execute_input": "2026-09-09T21:36:40.216821Z", "iopub.status.busy": "2026-09-09T21:36:40.216718Z", "iopub.status.idle": "2026-09-09T21:36:40.227923Z", "shell.execute_reply": "2026-09-09T21:36:40.227323Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "ABD assembled on ('DG', 1, (3, 3))\n" ] } ], "source": [ "E_VAL, NU_VAL, H_VAL, RHO_VAL = 4.32e8, 0.0, 0.2, 1.0\n", "\n", "material = hm.isotropic(domain, E=E_VAL, nu=NU_VAL, thickness=H_VAL, density=RHO_VAL)\n", "print(\"ABD assembled on\", material.A.space)" ] }, { "cell_type": "markdown", "id": "5f2fc2c5", "metadata": {}, "source": [ "## 5. Loads\n", "\n", "`hm.pressure` is a scalar load along the shell normal. `hm.traction` and `hm.moment`\n", "take explicit global vectors per unit area — on a curved surface those are very\n", "different loads, and self weight is a traction, not a pressure.\n", "\n", "`hm.edge_traction` / `hm.edge_moment` / `hm.edge_pressure` apply the same loads per\n", "unit length on the exterior facets a `where=` predicate selects, and `hm.point_load`\n", "places a consistent point force or moment at any physical coordinate. `Loads` compose\n", "with `+`, and each term keeps its own FE space all the way into the residual." ] }, { "cell_type": "code", "execution_count": 6, "id": "ad0adbcf", "metadata": { "execution": { "iopub.execute_input": "2026-09-09T21:36:40.229096Z", "iopub.status.busy": "2026-09-09T21:36:40.228999Z", "iopub.status.idle": "2026-09-09T21:36:40.244604Z", "shell.execute_reply": "2026-09-09T21:36:40.244287Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1 pressure term(s)\n" ] } ], "source": [ "PRESSURE_Z = 2.0\n", "loads = hm.pressure(domain, PRESSURE_Z)\n", "\n", "# composing is just addition -- e.g. adding a downward line load at the free tip:\n", "# loads = loads + hm.edge_traction(domain, [0.0, 0.0, -1.0], where=hm.near(\"x\", LENGTH))\n", "print(f\"{len(loads.pressure_terms)} pressure term(s)\")" ] }, { "cell_type": "markdown", "id": "83ae7061", "metadata": {}, "source": [ "## 6. Solve\n", "\n", "`hm.solve` is a CSDL implicit custom operation wrapping the assembly, the direct MUMPS\n", "factorization, and the adjoint. It returns a `ShellState` that back-references\n", "everything it was given, so the output functions need only the state.\n", "\n", "`material`, `loads` and `bcs` must all be built against this exact domain — DOLFINx\n", "silently drops a boundary condition located against a second, structurally identical\n", "function space, so Hermit checks rather than assumes." ] }, { "cell_type": "code", "execution_count": 7, "id": "61ba759a", "metadata": { "execution": { "iopub.execute_input": "2026-09-09T21:36:40.245533Z", "iopub.status.busy": "2026-09-09T21:36:40.245423Z", "iopub.status.idle": "2026-09-09T21:36:40.468550Z", "shell.execute_reply": "2026-09-09T21:36:40.468052Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "state.disp_solid: (1422,) mixed-space dofs\n" ] } ], "source": [ "state = hm.solve(domain, material, loads, bcs)\n", "print(\"state.disp_solid:\", state.disp_solid.shape, \"mixed-space dofs\")" ] }, { "cell_type": "markdown", "id": "dfb33e93", "metadata": {}, "source": [ "## 7. Scalar outputs\n", "\n", "Each output is a free function of the state. They are all `csdl.Variable`s, ready to\n", "become an objective or a constraint." ] }, { "cell_type": "code", "execution_count": 8, "id": "0c8531cc", "metadata": { "execution": { "iopub.execute_input": "2026-09-09T21:36:40.469484Z", "iopub.status.busy": "2026-09-09T21:36:40.469372Z", "iopub.status.idle": "2026-09-09T21:36:40.507835Z", "shell.execute_reply": "2026-09-09T21:36:40.507436Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "compliance : 1.387702e-01\n", "elastic energy : 6.938509e-02 (= compliance / 2)\n", "mass : 4.000000 (exact 4.0)\n", "centre of gravity: [5. 1. 0.]\n", "tip deflection : 8.676101e-03 (Euler-Bernoulli 8.680556e-03)\n" ] } ], "source": [ "compliance = hm.compliance(state)\n", "mass = hm.mass(state)\n", "cg = hm.center_of_gravity(state)\n", "energy = hm.elastic_energy(state)\n", "\n", "Ix = WIDTH * H_VAL**3 / 12.0\n", "eb_tip = PRESSURE_Z * WIDTH * LENGTH**4 / (8.0 * E_VAL * Ix)\n", "\n", "print(f\"compliance : {float(np.ravel(compliance.value)[0]):.6e}\")\n", "print(f\"elastic energy : {float(np.ravel(energy.value)[0]):.6e} (= compliance / 2)\")\n", "print(f\"mass : {float(np.ravel(mass.value)[0]):.6f} (exact {RHO_VAL * H_VAL * WIDTH * LENGTH})\")\n", "print(f\"centre of gravity: {np.round(np.ravel(cg.value), 6)}\")\n", "print(f\"tip deflection : {np.abs(state.disp_solid.value).max():.6e} (Euler-Bernoulli {eb_tip:.6e})\")" ] }, { "cell_type": "markdown", "id": "911f9f9e", "metadata": {}, "source": [ "## 8. Field outputs, and the frame they live in\n", "\n", "Field outputs come back as `Field` objects, which you can ask for on any DG or CG\n", "space and recover by L2 `project`, `interpolate`, or DG0 `average`.\n", "\n", "Frames matter here. A **DG** strain field defaults to the element-local in-plane\n", "frame, which follows the mesh parametrisation and is therefore a per-cell artefact: on\n", "an unstructured mesh, plotting the raw `xx` component gives a patchy picture of a\n", "smooth physical field. A **CG** field defaults to global Cartesian instead, because\n", "per-cell frames are ambiguous at a shared node — and because the components are then a\n", "full 3-D tensor, it carries six Voigt components rather than three. Use `to_global()`\n", "or `to_frame(...)` to move between representations.\n", "\n", "This flat, axis-aligned plate is the one case where the distinction does not show:\n", "every element frame already coincides with global $x$/$y$. `ex_orientation_fields.py`\n", "runs the same comparison on a triangle mesh, where the element frames alternate cell\n", "to cell and the three pictures are visibly different.\n" ] }, { "cell_type": "code", "execution_count": 9, "id": "3f019569", "metadata": { "execution": { "iopub.execute_input": "2026-09-09T21:36:40.508629Z", "iopub.status.busy": "2026-09-09T21:36:40.508525Z", "iopub.status.idle": "2026-09-09T21:36:40.599998Z", "shell.execute_reply": "2026-09-09T21:36:40.599647Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "curvature: (80, 3) (n_cells, [xx, yy, 2xy]), kind=strain2\n", " root cell, element frame : [ 0.00033015 -0. 0. ]\n", " root cell, global frame : [ 0.00033015 -0. 0. ]\n", " (identical here: the plate is flat and axis-aligned)\n", "\n", " CG1 curvature field: ('Lagrange', 1, (6,)), kind=tensor3 (global Voigt-6)\n", "\n", "nodal displacement: (105, 3); mean w at the free tip = 8.676101e-03\n", " (positive: pressure acts along +n, and the plate normal is +z)\n" ] } ], "source": [ "eps, kappa, gamma = hm.strain_fields(state) # DG2, element-local frame\n", "print(f\"curvature: {kappa.values.shape} (n_cells, [xx, yy, 2xy]), kind={kappa.kind}\")\n", "print(f\" root cell, element frame : {np.round(kappa.values[0], 8)}\")\n", "print(f\" root cell, global frame : {np.round(kappa.to_global().values[0], 8)}\")\n", "print(\" (identical here: the plate is flat and axis-aligned)\")\n", "\n", "_, kappa_cg, _ = hm.strain_fields(state, space=(\"Lagrange\", 1)) # smooth nodal\n", "print(f\"\\n CG1 curvature field: {kappa_cg.space}, kind={kappa_cg.kind} (global Voigt-6)\")\n", "\n", "u_nodal = hm.nodal_displacement(state).value # (n_nodes, 3), mesh-file order\n", "tip = domain.node_coords[:, 0] > LENGTH - 1e-9\n", "print(f\"\\nnodal displacement: {u_nodal.shape}; mean w at the free tip = {u_nodal[tip, 2].mean():.6e}\")\n", "print(\" (positive: pressure acts along +n, and the plate normal is +z)\")\n" ] }, { "cell_type": "markdown", "id": "1561a9e4", "metadata": {}, "source": [ "## 9. Fields as inputs: spaces and ordering\n", "\n", "`Field` is also how you give Hermit a *spatially varying* input. The builder you pick\n", "declares which ordering your array is in — that is the whole point of having several:\n", "\n", "| builder | input | space |\n", "|---|---|---|\n", "| `hm.constant` | one value | any |\n", "| `hm.from_nodal` | `(n_nodes, ...)`, **mesh-file** vertex order | CG1 |\n", "| `hm.from_cells` | `(n_cells, ...)`, **mesh-file** cell order | DG0 |\n", "| `hm.from_function` | a callable of the dof coordinates | any |\n", "| `hm.from_coeffs` | raw coefficients, **FE dof** order | any |\n", "\n", "Getting this wrong is silent — the two orderings coincide only when the mesh\n", "permutation happens to be the identity — so say which one you mean." ] }, { "cell_type": "code", "execution_count": 10, "id": "0beb18b0", "metadata": { "execution": { "iopub.execute_input": "2026-09-09T21:36:40.600973Z", "iopub.status.busy": "2026-09-09T21:36:40.600871Z", "iopub.status.idle": "2026-09-09T21:36:40.802850Z", "shell.execute_reply": "2026-09-09T21:36:40.802289Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "uniform t=0.2 : compliance 1.387702e-01, mass 4.0000\n", "tapered 0.3->0.1: compliance 6.413965e-02, mass 4.0000\n" ] } ], "source": [ "# a thickness that tapers from root to tip, given per mesh vertex\n", "taper = 0.3 - 0.02 * domain.node_coords[:, 0]\n", "tapered = hm.isotropic(domain, E=E_VAL, nu=NU_VAL,\n", " thickness=hm.from_nodal(domain, taper), density=RHO_VAL)\n", "tapered_state = hm.solve(domain, tapered, loads, bcs)\n", "\n", "print(f\"uniform t=0.2 : compliance {float(np.ravel(compliance.value)[0]):.6e}, \"\n", " f\"mass {float(np.ravel(mass.value)[0]):.4f}\")\n", "print(f\"tapered 0.3->0.1: compliance {float(np.ravel(hm.compliance(tapered_state).value)[0]):.6e}, \"\n", " f\"mass {float(np.ravel(hm.mass(tapered_state).value)[0]):.4f}\")" ] }, { "cell_type": "markdown", "id": "06ec48c7", "metadata": {}, "source": [ "## 10. Derivatives\n", "\n", "Everything above is one CSDL graph, so a total derivative is a `compute_totals` call.\n", "Reverse-mode sensitivities are available with respect to thickness, density, ABD, $E$,\n", "$\\nu$, ply angles and heights, every load coefficient, the fibre orientation, and the\n", "mesh coordinates.\n", "\n", "Note that the design variable has to be a `csdl.Variable` *before* it is built into a\n", "field, so that the graph records the dependence." ] }, { "cell_type": "code", "execution_count": 11, "id": "b3264009", "metadata": { "execution": { "iopub.execute_input": "2026-09-09T21:36:40.804114Z", "iopub.status.busy": "2026-09-09T21:36:40.804008Z", "iopub.status.idle": "2026-09-09T21:36:41.659378Z", "shell.execute_reply": "2026-09-09T21:36:41.658954Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "dC/dt: (105,), all negative = True\n", " most sensitive node at x = 0.50 (the root)\n", " least sensitive node at x = 10.00 (the tip)\n" ] } ], "source": [ "t = csdl.Variable(value=H_VAL * np.ones(domain.n_nodes), name=\"thickness\")\n", "dv_material = hm.isotropic(domain, E=E_VAL, nu=NU_VAL,\n", " thickness=hm.from_nodal(domain, t), density=RHO_VAL)\n", "dv_state = hm.solve(domain, dv_material, loads, bcs)\n", "dv_compliance = hm.compliance(dv_state)\n", "\n", "sim = csdl.experimental.PySimulator(recorder)\n", "dC_dt = np.asarray(sim.compute_totals([dv_compliance], [t])[dv_compliance, t]).ravel()\n", "\n", "print(f\"dC/dt: {dC_dt.shape}, all negative = {bool((dC_dt < 0).all())}\")\n", "print(f\" most sensitive node at x = {domain.node_coords[np.argmin(dC_dt), 0]:.2f} (the root)\")\n", "print(f\" least sensitive node at x = {domain.node_coords[np.argmax(dC_dt), 0]:.2f} (the tip)\")" ] }, { "cell_type": "code", "execution_count": 12, "id": "2c6192cf", "metadata": { "execution": { "iopub.execute_input": "2026-09-09T21:36:41.660728Z", "iopub.status.busy": "2026-09-09T21:36:41.660579Z", "iopub.status.idle": "2026-09-09T21:36:41.662538Z", "shell.execute_reply": "2026-09-09T21:36:41.662201Z" } }, "outputs": [], "source": [ "recorder.stop()" ] }, { "cell_type": "markdown", "id": "e63231fa", "metadata": {}, "source": [ "## Where to go next\n", "\n", "- **{doc}`../advanced_tutorials/optimization_and_shape`** — driving an optimizer, and\n", " differentiating with respect to the mesh coordinates.\n", "- **Examples** — `ex_composite_plate.py` (ply-angle sweep and its derivatives),\n", " `ex_fiber_orientation.py` and `ex_orientation_fields.py` (laminate orientation and\n", " strain-field frames), `ex_thickness_opt.py` (a full SLSQP run).\n", "- **Verification benchmarks** — 28 cases with published or closed-form references, each\n", " runnable on its own." ] } ], "metadata": { "kernelspec": { "display_name": "hermit", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.12.13" } }, "nbformat": 4, "nbformat_minor": 5 }