.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "gallery/elements/plane182/example_plane182.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_gallery_elements_plane182_example_plane182.py>`
        to download the full example code.

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_gallery_elements_plane182_example_plane182.py:


.. _ref_plane182_example:

QUAD4_PLANE — plane-stress uniaxial tension
===========================================

A square piece of plate is fixed on its left edge and pulled on its
right edge. We check the tip displacement against Hooke's law
``u_x = σ L / E`` and the Poisson contraction ``u_y = −ν u_x`` in
plane stress, then render the deformed mesh.

.. GENERATED FROM PYTHON SOURCE LINES 12-22

.. code-block:: Python


    from __future__ import annotations

    import numpy as np
    import pyvista as pv
    from vtkmodules.util.vtkConstants import VTK_QUAD

    import femorph_solver
    from femorph_solver import ELEMENTS








.. GENERATED FROM PYTHON SOURCE LINES 23-26

Problem set-up
--------------
Unit square of steel, 1 mm thick, pulled with 1 MN total along +x.

.. GENERATED FROM PYTHON SOURCE LINES 26-60

.. code-block:: Python


    E = 2.0e11  # Pa
    NU = 0.3
    RHO = 7850.0  # kg/m³
    THK = 1.0e-3  # m
    F_TOTAL = 1.0e6  # N

    coords = np.array(
        [(0.0, 0.0, 0.0), (1.0, 0.0, 0.0), (1.0, 1.0, 0.0), (0.0, 1.0, 0.0)],
        dtype=np.float64,
    )
    cells = np.array([4, 0, 1, 2, 3], dtype=np.int64)
    cell_types = np.array([VTK_QUAD], dtype=np.uint8)
    grid = pv.UnstructuredGrid(cells, cell_types, coords)

    m = femorph_solver.Model.from_grid(grid)
    m.assign(
        ELEMENTS.QUAD4_PLANE,
        material={"EX": E, "PRXY": NU, "DENS": RHO},
        real=(THK,),
    )
    # Plane-stress KEYOPT(3)=0 is the default; annotate it explicitly:
    m.materials[1]["_PLANE_MODE"] = "stress"

    # Clamp the x=0 edge in UX, and kill rigid-body UY without locking Poisson.
    for nn in (1, 4):
        m.fix(nodes=[nn], dof="UX")
    for nn in (1, 2):
        m.fix(nodes=[nn], dof="UY")

    F_each = F_TOTAL / 2.0
    for nn in (2, 3):
        m.apply_force(nn, fx=F_each)








.. GENERATED FROM PYTHON SOURCE LINES 61-66

Static solve + analytical comparison
------------------------------------
Out-of-plane stress is zero (plane stress), so the effective axial stress
on the x=1 edge is ``σ = F_TOTAL / (A · thickness)`` with edge length
1 m, giving strain ``σ / E`` and tip displacement ``σ · L / E``.

.. GENERATED FROM PYTHON SOURCE LINES 66-86

.. code-block:: Python


    res = m.solve()
    dof = m.dof_map()

    sigma = F_TOTAL / (1.0 * THK)
    ux_expected = sigma / E  # L = 1 m
    uy_expected = -NU * ux_expected

    print(f"Expected u_x (+x edge) = {ux_expected:.6e} m")
    for nn in (2, 3):
        row = np.where((dof[:, 0] == nn) & (dof[:, 1] == 0))[0][0]
        print(f"  node {nn} UX         = {res.displacement[row]:.6e}")
        assert np.isclose(res.displacement[row], ux_expected, rtol=1e-8)

    print(f"Expected u_y (+y edge) = {uy_expected:.6e} m")
    for nn in (3, 4):
        row = np.where((dof[:, 0] == nn) & (dof[:, 1] == 1))[0][0]
        print(f"  node {nn} UY         = {res.displacement[row]:.6e}")
        assert np.isclose(res.displacement[row], uy_expected, rtol=1e-8)





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    Expected u_x (+x edge) = 5.000000e-03 m
      node 2 UX         = 5.000000e-03
      node 3 UX         = 5.000000e-03
    Expected u_y (+y edge) = -1.500000e-03 m
      node 3 UY         = -1.500000e-03
      node 4 UY         = -1.500000e-03




.. GENERATED FROM PYTHON SOURCE LINES 87-89

Plot the deformed quad, coloured by displacement magnitude
----------------------------------------------------------

.. GENERATED FROM PYTHON SOURCE LINES 89-110

.. code-block:: Python


    grid = m.grid.copy()
    disp = np.zeros((grid.n_points, 3), dtype=np.float64)
    for i, nn in enumerate(grid.point_data["ansys_node_num"]):
        rows = np.where(dof[:, 0] == int(nn))[0]
        for r in rows:
            disp[i, int(dof[r, 1])] = res.displacement[r]
    grid.point_data["displacement"] = disp

    plotter = pv.Plotter(off_screen=True)
    warped = grid.warp_by_vector("displacement", factor=50.0)
    plotter.add_mesh(
        warped,
        scalars=np.linalg.norm(disp, axis=1),
        scalar_bar_args={"title": "disp [m]"},
        show_edges=True,
        cmap="viridis",
    )
    plotter.add_mesh(grid, style="wireframe", color="gray", opacity=0.4)
    plotter.view_xy()
    plotter.show()



.. image-sg:: /gallery/elements/plane182/images/sphx_glr_example_plane182_001.png
   :alt: example plane182
   :srcset: /gallery/elements/plane182/images/sphx_glr_example_plane182_001.png
   :class: sphx-glr-single-img






.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (0 minutes 0.229 seconds)


.. _sphx_glr_download_gallery_elements_plane182_example_plane182.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: example_plane182.ipynb <example_plane182.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: example_plane182.py <example_plane182.py>`

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: example_plane182.zip <example_plane182.zip>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_