{ "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# QUAD4 bilinear shape functions on the reference square\n\nRenders the four bilinear (4-node) shape functions\n$N_i(\\xi, \\eta)$ on the reference quadrilateral\n$\\hat\\Omega = [-1, +1]^2$:\n\n\\begin{align}N_i(\\xi, \\eta) = \\tfrac{1}{4}(1 + \\xi_i\\xi)(1 + \\eta_i\\eta),\n \\qquad i = 1, \\dots, 4,\\end{align}\n\nwith corner-node signs\n$(\\xi_i, \\eta_i) \\in \\{-1, +1\\}^2$. Each function equals 1\nat its own node and 0 at every other node \u2014 the Kronecker-delta\nproperty at the heart of standard Lagrange interpolation. The\nfour basis functions sum to 1 everywhere on the square (partition\nof unity), so any constant or linear field is reproduced exactly.\n\n## References\n\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, \u00a76.3 (4-node quad), Table 6.3-1.\n* Hughes, T. J. R. (2000) *The Finite Element Method \u2014\n Linear Static and Dynamic Finite Element Analysis*, Dover,\n \u00a73.5 + \u00a73.6.\n* Zienkiewicz, O. C., Taylor, R. L. (2013) *The Finite Element\n Method: Its Basis and Fundamentals*, 7th ed., \u00a76.4.\n\nImplementation: :class:`femorph_solver.elements.quad4_plane.Quad4Plane`.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from __future__ import annotations\n\nimport matplotlib.pyplot as plt\nimport numpy as np" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Corner signs in VTK_QUAD order: (-1,-1), (+1,-1), (+1,+1), (-1,+1).\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "corner_signs = np.array(\n [[-1, -1], [+1, -1], [+1, +1], [-1, +1]],\n dtype=float,\n)\n\n\ndef n_corner(xi, eta, signs):\n \"\"\"Bilinear corner shape function for a given (xi_i, eta_i).\"\"\"\n xi_i, eta_i = signs\n return 0.25 * (1 + xi_i * xi) * (1 + eta_i * eta)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Sample the reference square on a (xi, eta) grid\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "n = 81\nxi_grid, eta_grid = np.meshgrid(np.linspace(-1, 1, n), np.linspace(-1, 1, n), indexing=\"ij\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Render \u2014 2\u00d72 grid of contour plots, one per corner.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fig, axes = plt.subplots(2, 2, figsize=(9, 7), constrained_layout=True)\nfig.suptitle(\n \"QUAD4 bilinear shape functions on $\\\\hat\\\\Omega = [-1, +1]^2$\",\n fontsize=13,\n)\n\nfor idx, signs in enumerate(corner_signs):\n row, col = divmod(idx, 2)\n ax = axes[row, col]\n N = n_corner(xi_grid, eta_grid, signs)\n img = ax.contourf(xi_grid, eta_grid, N, levels=21, cmap=\"plasma\")\n ax.contour(\n xi_grid,\n eta_grid,\n N,\n levels=[0.0, 0.25, 0.5, 0.75, 1.0],\n colors=\"black\",\n linewidths=0.4,\n )\n sx, sy = signs.astype(int)\n ax.scatter(corner_signs[:, 0], corner_signs[:, 1], s=30, color=\"black\", zorder=5)\n ax.scatter([sx], [sy], s=110, edgecolor=\"black\", facecolor=\"white\", zorder=6)\n ax.set_aspect(\"equal\")\n ax.set_title(f\"$N_{{{idx + 1}}}$ $(\\\\xi_i, \\\\eta_i) = ({sx:+d}, {sy:+d})$\", fontsize=11)\n ax.set_xlabel(r\"$\\xi$\")\n ax.set_ylabel(r\"$\\eta$\")\n fig.colorbar(img, ax=ax, shrink=0.85)\n\nfig.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Verify partition of unity at every reference-element node.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "total = sum(n_corner(corner_signs[:, 0], corner_signs[:, 1], s) for s in corner_signs)\nnp.testing.assert_allclose(total, 1.0, atol=1e-14)\nprint(\"OK \u2014 QUAD4 partition-of-unity holds at every corner node.\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Verify Kronecker-delta interpolation at every reference-element node.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "basis = np.zeros((4, 4))\nfor k, (xi_k, eta_k) in enumerate(corner_signs):\n for i, signs in enumerate(corner_signs):\n basis[k, i] = n_corner(xi_k, eta_k, signs)\nnp.testing.assert_allclose(basis, np.eye(4), atol=1e-14)\nprint(\"OK \u2014 QUAD4 Kronecker-delta interpolation verified at every node.\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Verify partition of unity on a dense interior grid as well.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "dense_total = sum(n_corner(xi_grid, eta_grid, s) for s in corner_signs)\nnp.testing.assert_allclose(dense_total, 1.0, atol=1e-14)\nprint(\"OK \u2014 QUAD4 partition-of-unity holds across the full reference square.\")" ] } ], "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 }