{
  "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# QUAD4_SHELL \u2014 uniaxial stretch of a flat square plate\n\nSingle QUAD4_SHELL element (a 1 m \u00d7 1 m flat plate of thickness ``t``) in\npure in-plane tension. Out-of-plane DOFs are suppressed so the problem\nreduces to plane-stress membrane behaviour: ``u_x = \u03c3 L / E`` and the\nPoisson contraction ``u_y = \u2212\u03bd u_x``.\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_QUAD\n\nimport femorph_solver\nfrom femorph_solver import ELEMENTS"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Problem data\nSteel plate, 10 mm thick, pulled with a total 100 kN along +x.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "E = 2.1e11\nNU = 0.3\nRHO = 7850.0\nTHK = 0.010  # 10 mm\nL = 1.0\nF_TOTAL = 1.0e5  # N spread over the +x edge (two corner nodes)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Build the model\nCorner ordering is VTK_QUAD: node 1 = (0,0), 2 = (1,0), 3 = (1,1),\n4 = (0,1). We clamp UX on the -x edge, anchor UY on a single corner\nto kill the in-plane translation mode, and suppress every\nout-of-plane DOF (UZ, ROTX, ROTY, ROTZ) so the element behaves like\na pure membrane.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "coords = np.array(\n    [\n        [0.0, 0.0, 0.0],\n        [L, 0.0, 0.0],\n        [L, L, 0.0],\n        [0.0, L, 0.0],\n    ],\n    dtype=np.float64,\n)\n\ncells = np.array([4, 0, 1, 2, 3], dtype=np.int64)\ncell_types = np.array([VTK_QUAD], dtype=np.uint8)\ngrid = pv.UnstructuredGrid(cells, cell_types, coords)\n\nm = femorph_solver.Model.from_grid(grid)\nm.assign(\n    ELEMENTS.QUAD4_SHELL,\n    material={\"EX\": E, \"PRXY\": NU, \"DENS\": RHO},\n    real=(THK,),\n)\n\n# -x edge clamped in UX\nm.fix(nodes=[1], dof=\"UX\")\nm.fix(nodes=[4], dof=\"UX\")\n# Anchor UY on one node to stop rigid-body in-plane translation\nm.fix(nodes=[1], dof=\"UY\")\n# Pure membrane: zero all out-of-plane DOFs everywhere\nfor nn in (1, 2, 3, 4):\n    m.fix(nodes=[nn], dof=\"UZ\")\n    m.fix(nodes=[nn], dof=\"ROTX\")\n    m.fix(nodes=[nn], dof=\"ROTY\")\n    m.fix(nodes=[nn], dof=\"ROTZ\")\n\nF_each = F_TOTAL / 2.0\nm.apply_force(2, fx=F_each)\nm.apply_force(3, fx=F_each)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Static solve + analytical comparison\nStress on the +x face is ``\u03c3 = F / (width \u00b7 t)``. On a 1 m width with\n10 mm thickness: ``\u03c3 = 10 MPa``. Expected ``u_x = \u03c3 L / E`` and\n``u_y = \u2212\u03bd u_x`` on the free corner nodes.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "res = m.solve_static()\ndof = m.dof_map()\n\nsigma = F_TOTAL / (L * THK)\nux_expected = sigma * L / E\nuy_expected = -NU * ux_expected\n\nprint(f\"In-plane stress \u03c3    = {sigma:.3e} Pa\")\nprint(f\"Expected u_x         = {ux_expected:.6e} m\")\nfor nn in (2, 3):\n    row = np.where((dof[:, 0] == nn) & (dof[:, 1] == 0))[0][0]\n    val = res.displacement[row]\n    print(f\"  node {nn} UX        = {val:.6e}\")\n    assert np.isclose(val, ux_expected, rtol=1e-8)\n\n# Node 3 is free in UY (node 1 is anchored, so check there)\nrow_uy3 = np.where((dof[:, 0] == 3) & (dof[:, 1] == 1))[0][0]\nuy3 = res.displacement[row_uy3]\nprint(f\"Expected u_y (node 3) = {uy_expected:.6e}\")\nprint(f\"Computed u_y (node 3) = {uy3:.6e}\")\nassert np.isclose(uy3, uy_expected, rtol=1e-6)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Plot the deformed plate\nScatter displacements onto point_data and warp the mesh. A modest\nmagnification (~1000\u00d7) makes the Poisson contraction visible.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "grid = m.grid.copy()\ndisp = np.zeros((grid.n_points, 3), dtype=np.float64)\nfor i, nn in enumerate(grid.point_data[\"ansys_node_num\"]):\n    rows = np.where(dof[:, 0] == int(nn))[0]\n    for r in rows:\n        d_idx = int(dof[r, 1])\n        if d_idx < 3:\n            disp[i, d_idx] = res.displacement[r]\ngrid.point_data[\"displacement\"] = disp\n\nplotter = pv.Plotter(off_screen=True)\nplotter.add_mesh(grid, style=\"wireframe\", color=\"gray\", line_width=2)\nplotter.add_mesh(\n    grid.warp_by_vector(\"displacement\", factor=2.0e3),\n    scalars=np.linalg.norm(disp, axis=1),\n    show_edges=True,\n    scalar_bar_args={\"title\": \"|u| [m]\"},\n)\nplotter.add_axes()\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
}