{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import pyvista\npyvista.OFF_SCREEN = True\npyvista.set_jupyter_backend('static')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n\n# HEX8 \u2014 elastic-strain post-processing\n\nSolve a HEX8 flat plate under uniaxial tension and recover the full\n6-component elastic-strain tensor on the mesh with\n:meth:`femorph_solver.result.StaticResult.elastic_strain` \u2014 element-nodal\nstrain averaged onto every grid point.\n\n``result.elastic_strain(model=m)`` returns the nodal-averaged Voigt\nstrain ``(n_points, 6)`` (the canonical post-processing output);\n:meth:`~femorph_solver.result.StaticResult.elastic_strain_per_element`\nreturns the per-element dict ``{elem_num: (n_nodes_in_elem, 6)}`` for\nverification workflows that need raw element-by-element values. Strain\nis computed at each element's own nodes as\n$\\varepsilon(\\xi_\\text{node}) = B(\\xi_\\text{node})\\cdot u_e$\n\u2014 no RST round-trip, no disk write.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from __future__ import annotations\n\nimport numpy as np\nimport pyvista as pv\nfrom vtkmodules.util.vtkConstants import VTK_HEXAHEDRON\n\nimport femorph_solver\nfrom femorph_solver import ELEMENTS" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Problem setup\nA 1 m \u00d7 0.4 m \u00d7 0.05 m steel plate meshed as a 20 \u00d7 8 \u00d7 1 HEX8\nbrick (160 elements). The ``x = 0`` face is held in ``UX`` (symmetry),\na single pin at the origin kills the ``UY`` / ``UZ`` rigid-body modes,\nand the ``x = LX`` face is pulled by a total force ``F`` split over\nits corner nodes.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "E = 2.1e11 # Pa\nNU = 0.30\nRHO = 7850.0\nLX, LY, LZ = 1.0, 0.4, 0.05\nNX, NY, NZ = 20, 8, 1\nF_TOTAL = 1.0e5 # N\n\nxs = np.linspace(0.0, LX, NX + 1)\nys = np.linspace(0.0, LY, NY + 1)\nzs = np.linspace(0.0, LZ, NZ + 1)\nxx, yy, zz = np.meshgrid(xs, ys, zs, indexing=\"ij\")\npoints = np.stack([xx.ravel(), yy.ravel(), zz.ravel()], axis=1)\n\n\n# Hex connectivity in VTK_HEXAHEDRON order (0-based VTK indices).\ndef _node_idx(i: int, j: int, k: int) -> int:\n return (i * (NY + 1) + j) * (NZ + 1) + k\n\n\ncells_flat: list[int] = []\nfor i in range(NX):\n for j in range(NY):\n for k in range(NZ):\n cells_flat.extend(\n [\n 8,\n _node_idx(i, j, k),\n _node_idx(i + 1, j, k),\n _node_idx(i + 1, j + 1, k),\n _node_idx(i, j + 1, k),\n _node_idx(i, j, k + 1),\n _node_idx(i + 1, j, k + 1),\n _node_idx(i + 1, j + 1, k + 1),\n _node_idx(i, j + 1, k + 1),\n ]\n )\n\nn_cells = NX * NY * NZ\ncell_types = np.full(n_cells, VTK_HEXAHEDRON, dtype=np.uint8)\ngrid = pv.UnstructuredGrid(np.asarray(cells_flat, dtype=np.int64), cell_types, points)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Build the femorph-solver model\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "m = femorph_solver.Model.from_grid(grid)\nm.assign(ELEMENTS.HEX8, material={\"EX\": E, \"PRXY\": NU, \"DENS\": RHO})\n\nnode_nums = np.asarray(m.grid.point_data[\"ansys_node_num\"])\npts = np.asarray(m.grid.points)\n\n# Symmetry BC: x=0 face clamped in UX; single pin at the origin in UY/UZ.\nx0_nodes = node_nums[pts[:, 0] < 1e-9].tolist()\nm.fix(nodes=x0_nodes, dof=\"UX\")\norigin_nodes = node_nums[(pts[:, 0] < 1e-9) & (pts[:, 1] < 1e-9) & (pts[:, 2] < 1e-9)].tolist()\nm.fix(nodes=origin_nodes, dof=\"UY\")\nm.fix(nodes=origin_nodes, dof=\"UZ\")\n\n# Traction on x=LX face: split F_TOTAL over its nodes.\nx_end_nodes = node_nums[pts[:, 0] > LX - 1e-9].tolist()\nfx_each = F_TOTAL / len(x_end_nodes)\nfor nn in x_end_nodes:\n m.apply_force(int(nn), fx=fx_each)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Static solve\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "res = m.solve_static()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Recover elastic strain\nDefault call returns nodal-averaged strain of shape ``(n_points, 6)``:\ncolumns are ``[\u03b5xx, \u03b5yy, \u03b5zz, \u03b3xy, \u03b3yz, \u03b3xz]`` with *engineering*\nshears (canonical Voigt strain-recovery output).\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "eps = res.elastic_strain(model=m)\nprint(f\"eps shape: {eps.shape}\")\n\n# Analytical: uniform \u03c3xx = F_TOTAL / (LY \u00b7 LZ), \u03b5xx = \u03c3 / E,\n# \u03b5yy = \u03b5zz = -\u03bd \u00b7 \u03b5xx.\nsigma_xx = F_TOTAL / (LY * LZ)\neps_xx_expected = sigma_xx / E\neps_yy_expected = -NU * eps_xx_expected\nprint(f\"\u03b5xx expected = {eps_xx_expected:.3e}\")\nprint(f\"\u03b5xx recovered (mean over nodes) = {eps[:, 0].mean():.3e}\")\nprint(f\"\u03b5yy recovered (mean) = {eps[:, 1].mean():.3e}\")\nprint(f\"\u03b5yy analytical = {eps_yy_expected:.3e}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Per-element arrays keyed by element number \u2014 useful when you want to\nsee jumps at element boundaries or compute element-wise strain norms.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "per_elem = res.elastic_strain_per_element(model=m)\nfirst_elem = next(iter(per_elem))\nprint(\n f\"per-element dict has {len(per_elem)} elements; \"\n f\"first key = {first_elem}, \"\n f\"strain block shape = {per_elem[first_elem].shape}\"\n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Visualise \u03b5xx on the deformed mesh\n:func:`femorph_solver.io.static_result_to_grid` scatters the per-node\ndisplacement onto ``(n_points, 3)`` UX/UY/UZ point data in one call \u2014\nno hand-rolled dof-map loop required. ``elastic_strain`` is already\nindexed in grid-point order so \u03b5xx drops straight onto the grid.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "grid = femorph_solver.io.static_result_to_grid(m, res)\ngrid.point_data[\"eps_xx\"] = eps[:, 0]\n\nwarped = grid.warp_by_vector(\"displacement\", factor=200.0)\n\nplotter = pv.Plotter(off_screen=True)\nplotter.add_mesh(\n grid,\n style=\"wireframe\",\n color=\"gray\",\n line_width=1,\n opacity=0.4,\n label=\"undeformed\",\n)\nplotter.add_mesh(\n warped,\n scalars=\"eps_xx\",\n show_edges=True,\n cmap=\"viridis\",\n scalar_bar_args={\"title\": \"\u03b5xx\"},\n label=\"\u03b5xx (deformed \u00d7200)\",\n)\nplotter.add_legend()\nplotter.add_axes()\nplotter.view_xy()\nplotter.show()" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "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.3" } }, "nbformat": 4, "nbformat_minor": 0 }