{
  "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# SPRING \u2014 2-node longitudinal spring stiffness\n\nThe 2-node 3D longitudinal spring carries an axial-only\nstiffness $k$ between two nodes connected by a unit\ndirection vector $\\hat d = (l, m, n)$.  No mass.  Three\ntranslational DOFs per node, 6 DOFs per element.\n\nLocal axial-only stiffness on the natural coordinate\n$s \\in [-1, +1]$:\n\n\\begin{align}K_\\text{local} = k \\begin{bmatrix} +1 & -1 \\\\ -1 & +1 \\end{bmatrix},\\end{align}\n\nembedded in 3D via the direction-cosine block\n$T = \\hat d\\, \\hat d^{\\!\\top}$:\n\n\\begin{align}K_e = k \\begin{bmatrix} +T & -T \\\\ -T & +T \\end{bmatrix}.\\end{align}\n\n## References\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, \u00a73.2.\n* Przemieniecki, J. S. (1968) *Theory of Matrix Structural\n  Analysis*, McGraw-Hill, \u00a73.\n\nImplementation: :class:`femorph_solver.elements.spring.Spring`.\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": [
        "## 2-node spring, drawn along an arbitrary direction\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Pick a non-axial direction so the figure shows the rotation\n# from local to global axis.\nend_a = np.array([0.0, 0.0, 0.0])\nend_b = np.array([1.0, 0.6, 0.4])\nlength = float(np.linalg.norm(end_b - end_a))\ndirection = (end_b - end_a) / length\n\n# Spring as a coil polyline along the I\u2192J axis.\nn_turns = 8\nn_pts = 240\nt = np.linspace(0.0, 1.0, n_pts)\naxial = end_a + np.outer(t, end_b - end_a)\n# Build two perpendicular vectors to ``direction`` for the coil radius.\nref = np.array([0.0, 0.0, 1.0]) if abs(direction[2]) < 0.9 else np.array([1.0, 0.0, 0.0])\ne1 = np.cross(direction, ref)\ne1 /= np.linalg.norm(e1)\ne2 = np.cross(direction, e1)\ncoil_radius = 0.05 * length\nphase = 2.0 * np.pi * n_turns * t\ncoil = axial + coil_radius * np.outer(np.cos(phase), e1) + coil_radius * np.outer(np.sin(phase), e2)\ncoil_pd = pv.lines_from_points(coil)"
      ]
    },
    {
      "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, 360))\nplotter.add_mesh(\n    coil_pd, color=\"#1f77b4\", line_width=2.5, label=\"k between nodes I and J (k = const)\"\n)\nplotter.add_points(\n    np.vstack([end_a, end_b]),\n    render_points_as_spheres=True,\n    point_size=18,\n    color=\"black\",\n    label=\"nodes I, J\",\n)\nplotter.add_axes(line_width=4, color=\"black\")\nplotter.view_isometric()\nplotter.camera.zoom(1.1)\nplotter.add_legend(face=None, size=(0.30, 0.10), bcolor=\"white\")\nplotter.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Verify the global stiffness assembly\n\nSymbolic check that ``K_e = k [[+T, -T], [-T, +T]]`` matches\nwhat the Cook \u00a73.2 derivation predicts on this direction\nvector.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "k = 1.0e3\nT = np.outer(direction, direction)\nK_e = np.block([[+k * T, -k * T], [-k * T, +k * T]])\n\n# Sanity \u2014 moving the I node by -d (away from J along the axis)\n# stretches the spring by 1 m.  ``K_e @ u`` returns the *external*\n# load that produces that displacement: pull on I in -d (force\n# magnitude k) and push on J in +d (force magnitude k).  The\n# internal restoring force is the negative of those.\nu = np.zeros(6)\nu[0:3] = -direction\nf_ext = K_e @ u\nnp.testing.assert_allclose(f_ext[0:3], -k * direction, atol=1e-12)\nnp.testing.assert_allclose(f_ext[3:6], +k * direction, atol=1e-12)\nprint(f\"K_e on direction d = {direction.round(4)}:\")\nprint(f\"  k = {k:.0f} N/m\")\nprint(\"  external load to displace I by -d: F_I = -k d, F_J = +k d  \u2713\")"
      ]
    }
  ],
  "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
}