{
  "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 shell reference geometry \u2014 Mindlin\u2013Reissner shape functions\n\nThe 4-node flat Mindlin\u2013Reissner shell maps to the natural\ncoordinates $(\\xi, \\eta) \\in [-1, 1]^2$ and carries six\nDOFs per node \u2014 three translations\n$(u_x, u_y, u_z)$ plus three rotations\n$(\\theta_x, \\theta_y, \\theta_z)$ \u2014 for 24 DOFs per element.\n\nIn-plane shape functions (bilinear Lagrange):\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\nThe first-order shear-deformation theory (Mindlin 1951;\nReissner 1945) assumes the kinematic split\n\n\\begin{align}u_x = u_x^0 + z\\, \\theta_y, \\quad\n    u_y = u_y^0 - z\\, \\theta_x, \\quad\n    u_z = u_z^0,\\end{align}\n\nso the bending-membrane coupling is governed by independent\nrotation DOFs and the through-thickness shear strain is constant\nacross $z$.\n\nIntegration:\n\n* **Membrane + bending** \u2014 2\u00d72 Gauss-Legendre (4 points), exact\n  for the bilinear $B^T D B$ integrand on the reference\n  square.\n* **Transverse shear** \u2014 1\u00d71 selective-reduced Gauss (1 point at\n  the centre), per Malkus & Hughes (CMAME 1978), to suppress\n  shear locking on thin shells.\n\nDrilling-DOF stabilisation: a per-node penalty\n$\\alpha\\, G\\, t\\, A$ (Allman 1984; Hughes-Brezzi 1989)\nremoves the local rank deficiency when $\\theta_z$ is free.\n\n## References\n* Mindlin, R. D. (1951) \"Influence of rotatory inertia and shear\n  on flexural motions of isotropic elastic plates,\" *J. Applied\n  Mechanics* 18, 31\u201338.\n* Reissner, E. (1945) \"The effect of transverse shear deformation\n  on the bending of elastic plates,\" *J. Applied Mechanics* 12,\n  A-69\u2013A-77.\n* Malkus, D. S. and Hughes, T. J. R. (1978) \"Mixed finite element\n  methods \u2014 Reduced and selective integration techniques: A\n  unification of concepts,\" *CMAME* 15 (1), 63\u201381.\n* Allman, D. J. (1984) \"A compatible triangular element including\n  vertex rotations for plane elasticity analysis,\" *Computers &\n  Structures* 19 (1-2), 1\u20138.\n\nImplementation: :class:`femorph_solver.elements.quad4_shell.Quad4Shell`.\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 with both quadrature rules\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)\n\n# 2x2 Gauss for membrane + bending\ng = 1.0 / np.sqrt(3.0)\ngauss_4 = np.array([[a, b, 0.0] for a in (-g, +g) for b in (-g, +g)])\n\n# 1x1 selective-reduced for transverse shear\ngauss_1 = np.array([[0.0, 0.0, 0.0]])"
      ]
    },
    {
      "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=(640, 480))\nplotter.add_mesh(\n    ref_quad, style=\"wireframe\", color=\"black\", line_width=2.5, label=\"reference QUAD4 shell\"\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_4,\n    render_points_as_spheres=True,\n    point_size=14,\n    color=\"#d62728\",\n    label=\"2\u00d72 membrane / bending\",\n)\nplotter.add_points(\n    gauss_1,\n    render_points_as_spheres=True,\n    point_size=22,\n    color=\"#9467bd\",\n    label=\"1\u00d71 transverse shear\",\n)\nplotter.view_xy()\nplotter.camera.zoom(1.3)\nplotter.add_legend(face=None, size=(0.27, 0.14), bcolor=\"white\")\nplotter.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Sanity \u2014 N_i at all 5 quadrature locations sums to 1\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])\nall_pts = np.vstack([gauss_4, gauss_1])\nsums = []\nfor q in all_pts:\n    sums.append(sum(0.25 * (1 + xi_i[i] * q[0]) * (1 + eta_i[i] * q[1]) for i in range(4)))\nsums = np.array(sums)\nnp.testing.assert_allclose(sums, 1.0, atol=1e-14)\nprint(f\"5 quadrature locations; max |\u03a3 N_i - 1| = {np.max(np.abs(sums - 1)):.2e}\")\nprint(\"OK \u2014 partition of unity verified at every Gauss 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
}