{ "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# TET10 reference geometry \u2014 corner + mid-edge nodes, Keast points\n\nRenders the 10-node quadratic tetrahedron's reference element\n$\\{(\\xi_1, \\xi_2, \\xi_3) : \\xi_i \\ge 0, \\sum_i \\xi_i \\le 1\\}$:\n\n* 4 **corner nodes** at the unit-tet vertices.\n* 6 **mid-edge nodes** at the centre of every edge.\n* **Keast (1986) 4-point** quadrature rule used for both\n $K^e$ and $M^e$ \u2014 degree-2 exact, weights $1/24$.\n\nShape functions (volume coordinates $L_i$):\n\n\\begin{align}N^c_i = L_i (2 L_i - 1), \\quad i = 1, \\ldots, 4\n \\qquad\\text{and}\\qquad\n N^m_{ij} = 4 L_i L_j\\end{align}\n\nfor the four corners and six mid-edges respectively (Cook 2002\nTable 6.5-1).\n\nKeast 4-point rule on the unit tet:\n\n\\begin{align}(\\xi_1, \\xi_2, \\xi_3) = \\left(\\tfrac{1 \\pm 3\\alpha}{4},\n \\tfrac{1 \\mp \\alpha}{4}, \\tfrac{1 \\mp \\alpha}{4}\\right),\n \\quad \\alpha = \\tfrac{1}{\\sqrt{5}},\n \\quad w_q = \\tfrac{1}{24}.\\end{align}\n\n## References\n* Keast, P. (1986) \"Moderate-degree tetrahedral quadrature\n formulas,\" *Computer Methods in Applied Mechanics and\n Engineering* 55, 339\u2013348.\n* Cook, R. D., Malkus, D. S., Plesha, M. E., Witt, R. J. (2002)\n *Concepts and Applications of Finite Element Analysis*, 4th\n ed., Wiley, \u00a79.2\u2013\u00a79.3.\n* Hughes, T. J. R. (2000) *The Finite Element Method*, Dover, \u00a73.8.\n\nImplementation: :class:`femorph_solver.elements.tet10.Tet10`.\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" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Node layout \u2014 4 corners + 6 mid-edges\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "corners = np.array(\n [\n [0.0, 0.0, 0.0],\n [1.0, 0.0, 0.0],\n [0.0, 1.0, 0.0],\n [0.0, 0.0, 1.0],\n ]\n)\nmid_edges = np.array(\n [\n # ordering matches VTK_QUADRATIC_TETRA: edges 0-1, 1-2, 2-0,\n # 0-3, 1-3, 2-3.\n 0.5 * (corners[0] + corners[1]),\n 0.5 * (corners[1] + corners[2]),\n 0.5 * (corners[2] + corners[0]),\n 0.5 * (corners[0] + corners[3]),\n 0.5 * (corners[1] + corners[3]),\n 0.5 * (corners[2] + corners[3]),\n ]\n)\nall_nodes = np.vstack([corners, mid_edges])\n\ncells = np.hstack([[10], np.arange(10, dtype=np.int64)])\ncell_types = np.array([pv.CellType.QUADRATIC_TETRA], dtype=np.uint8)\nref_tet = pv.UnstructuredGrid(cells, cell_types, all_nodes)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Keast 4-point integration points\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "alpha = 1.0 / np.sqrt(5.0)\nplus = (1.0 + 3.0 * alpha) / 4.0\nminus = (1.0 - alpha) / 4.0\nkeast = np.array(\n [\n [plus, minus, minus],\n [minus, plus, minus],\n [minus, minus, plus],\n [minus, minus, minus],\n ]\n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Render the figure\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "plotter = pv.Plotter(off_screen=True, window_size=(640, 480))\nplotter.add_mesh(\n ref_tet,\n style=\"wireframe\",\n color=\"black\",\n line_width=2,\n label=\"reference TET10\",\n)\nplotter.add_points(\n corners,\n render_points_as_spheres=True,\n point_size=18,\n color=\"black\",\n label=\"corner nodes (4)\",\n)\nplotter.add_points(\n mid_edges,\n render_points_as_spheres=True,\n point_size=12,\n color=\"#1f77b4\",\n label=\"mid-edge nodes (6)\",\n)\nplotter.add_points(\n keast,\n render_points_as_spheres=True,\n point_size=16,\n color=\"#d62728\",\n label=\"Keast 4-pt (K + M)\",\n)\nplotter.add_axes(line_width=4, color=\"black\")\nplotter.view_isometric()\nplotter.camera.zoom(1.1)\nplotter.add_legend(face=None, size=(0.22, 0.14), bcolor=\"white\")\nplotter.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Sanity: Keast weights sum to the unit-tet volume = 1/6\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "w_q = 1.0 / 24.0\nsum_w = 4 * w_q\nassert abs(sum_w - 1.0 / 6.0) < 1e-15\nprint(f\"Keast 4-pt: \u03a3 w_q = {sum_w:.6f}; unit-tet volume = {1 / 6:.6f}. Match.\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Verify the partition of unity at the Keast points\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "L = np.column_stack([1.0 - keast.sum(axis=1), keast[:, 0], keast[:, 1], keast[:, 2]])\n# Corner shape functions L_i (2 L_i - 1).\nN_corner = L * (2.0 * L - 1.0)\n# Mid-edge shape functions 4 L_i L_j over the 6 edges.\nedge_pairs = [(0, 1), (1, 2), (2, 0), (0, 3), (1, 3), (2, 3)]\nN_mid = np.column_stack([4.0 * L[:, i] * L[:, j] for i, j in edge_pairs])\nsums = N_corner.sum(axis=1) + N_mid.sum(axis=1)\nprint(\"partition-of-unity \u03a3 N_i at each Keast point:\")\nfor q, s in zip(keast, sums, strict=True):\n print(f\" \u03be\u2081={q[0]:+.4f} \u03be\u2082={q[1]:+.4f} \u03be\u2083={q[2]:+.4f} \u03a3={s:.15f}\")\nnp.testing.assert_allclose(sums, 1.0, atol=1.0e-14)\nprint(\"OK \u2014 all 10 shape functions sum to 1 at every Keast point.\")" ] } ], "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 }