{
  "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 plane reference geometry \u2014 bilinear shape functions + 2\u00d72 Gauss\n\nThe 4-node bilinear plane element maps to the natural\ncoordinates $(\\xi, \\eta) \\in [-1, 1]^2$.  Two\ntranslational DOFs per node; 8 DOFs per element.\n\nShape functions (bilinear Lagrange on the reference square):\n\n\\begin{align}N_i(\\xi, \\eta) = \\tfrac{1}{4}(1 + \\xi_i\\xi)(1 + \\eta_i\\eta),\n    \\qquad i = 1, \\ldots, 4,\\end{align}\n\nwith $(\\xi_i, \\eta_i) \\in \\{-1, +1\\}^2$ at the corresponding\ncorners.  Stiffness and consistent mass use 2\u00d72 Gauss-Legendre\n(4 points), exact for the bilinear $B^T D B$ and\n$N^T N$ integrands.\n\nPlane stress vs plane strain selected by the ``_PLANE_MODE``\nmaterial flag (Cook 2002 \u00a73.5\u2013\u00a73.6).\n\n## References\n* Zienkiewicz, O. C. and Taylor, R. L. (2013) *The Finite Element\n  Method*, 7th ed., \u00a76 + App. G.\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, \u00a75, \u00a711.3.\n* Hughes, T. J. R. (2000) *The Finite Element Method*, Dover, \u00a73.1.\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 numpy as np\nimport pyvista as pv"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Reference quad\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "corners = np.array([[-1.0, -1.0, 0.0], [+1.0, -1.0, 0.0], [+1.0, +1.0, 0.0], [-1.0, +1.0, 0.0]])\ncells = np.hstack([[4], np.arange(4, dtype=np.int64)])\ncell_types = np.array([pv.CellType.QUAD], dtype=np.uint8)\nref_quad = pv.UnstructuredGrid(cells, cell_types, corners)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "2\u00d72 Gauss-Legendre on the reference square\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "g = 1.0 / np.sqrt(3.0)\ngauss = np.array([[a, b, 0.0] for a in (-g, +g) for b in (-g, +g)])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Render\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plotter = pv.Plotter(off_screen=True, window_size=(560, 480))\nplotter.add_mesh(\n    ref_quad, style=\"wireframe\", color=\"black\", line_width=2.5, label=\"reference QUAD4\"\n)\nplotter.add_points(\n    corners, render_points_as_spheres=True, point_size=18, color=\"black\", label=\"corner nodes (4)\"\n)\nplotter.add_points(\n    gauss, render_points_as_spheres=True, point_size=16, color=\"#d62728\", label=\"2\u00d72 Gauss\"\n)\nplotter.view_xy()\nplotter.camera.zoom(1.4)\nplotter.add_legend(face=None, size=(0.22, 0.10), bcolor=\"white\")\nplotter.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Verify partition of unity at the Gauss points\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "xi_i = np.array([-1.0, +1.0, +1.0, -1.0])\neta_i = np.array([-1.0, -1.0, +1.0, +1.0])\nsums = []\nfor q in gauss:\n    s = sum(0.25 * (1 + xi_i[i] * q[0]) * (1 + eta_i[i] * q[1]) for i in range(4))\n    sums.append(s)\nsums = np.array(sums)\nnp.testing.assert_allclose(sums, 1.0, atol=1e-14)\nprint(f\"4 Gauss points; \u03a3 N_i(\u03be_q, \u03b7_q) = {sums.round(15)} \u2014 partition of unity verified.\")"
      ]
    }
  ],
  "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
}