{
  "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# SDOF transient \u2014 step-load response\n\nA single-DOF mass-spring-damper subjected to a step load shows off the\nNewmark-\u03b2 transient integrator in\n:func:`~femorph_solver.solvers.transient.solve_transient`.  We integrate\nfor long enough to watch the underdamped oscillation settle toward the\nnew static equilibrium, then compare the result to the textbook\nsecond-order response.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from __future__ import annotations\n\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport scipy.sparse as sp\n\nfrom femorph_solver.solvers.transient import solve_transient"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Set up a 1-DOF oscillator\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "m = 1.0  # kg\nk = 4.0 * np.pi**2  # N / m \u2014 gives \u03c9n = 2\u03c0 rad/s (fn = 1 Hz)\nzeta = 0.05\nc = 2.0 * zeta * np.sqrt(k * m)\n\nM = sp.csr_array(np.array([[m]], dtype=float))\nK = sp.csr_array(np.array([[k]], dtype=float))\nC = sp.csr_array(np.array([[c]], dtype=float))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Step load at t = 0\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "F_step = 1.0  # N\nF = np.array([F_step])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Newmark-\u03b2 with the default unconditionally-stable parameters\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "dt = 1e-3\nn_steps = 5000\nresult = solve_transient(K, M, F, dt=dt, n_steps=n_steps, C=C)\nu_fs = result.displacement[:, 0]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Analytical comparison\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "omega_n = np.sqrt(k / m)\nomega_d = omega_n * np.sqrt(1 - zeta**2)\nu_static = F_step / k\nphi = np.arctan(zeta / np.sqrt(1 - zeta**2))\nu_exact = u_static * (\n    1\n    - np.exp(-zeta * omega_n * result.time)\n    / np.sqrt(1 - zeta**2)\n    * np.cos(omega_d * result.time - phi)\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Plot Newmark vs analytic\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig, ax = plt.subplots(figsize=(6, 3))\nax.plot(result.time, u_fs, label=\"Newmark-\u03b2\", color=\"#1f77b4\")\nax.plot(result.time, u_exact, \"--\", label=\"analytic\", color=\"black\", linewidth=0.8)\nax.axhline(\n    u_static,\n    color=\"red\",\n    linestyle=\":\",\n    linewidth=0.8,\n    label=f\"static limit u_s = {u_static:.4f}\",\n)\nax.set_xlabel(\"time [s]\")\nax.set_ylabel(\"displacement [m]\")\nax.set_title(\"SDOF step-load response (\u03b6 = 0.05)\")\nax.legend(loc=\"lower right\")\nax.grid(True, alpha=0.3)\nfig.tight_layout()\n\nerr = np.max(np.abs(u_fs - u_exact)) / u_static\nprint(f\"max relative error vs analytic: {err:.3e}\")"
      ]
    }
  ],
  "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
}