{
  "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# Modal participation factors + effective modal mass\n\nFor base-excitation analysis (seismic, shock, hard-mounted\nrandom vibration), the **modal participation factor**\n$\\Gamma_{i}$ and **effective modal mass**\n$M_{\\mathrm{eff},i}$ of each mode tell you which modes\nwill respond when the structure is shaken in a given direction \u2014\nand which can safely be ignored.  Together they answer the\npractical engineering question:\n\n    *How many modes do I need to capture 90 % of the response\n    in the y direction?*\n\nDefinitions (Bathe 2014 \u00a79.3.4; Chopra 2017 \u00a713.1):\n\n\\begin{align}:label: modal-participation\n\n    \\Gamma_{i}^{(d)}\n    \\;=\\;\n    \\frac{\\phi_{i}^{T}\\, M\\, r^{(d)}}{\\phi_{i}^{T}\\, M\\, \\phi_{i}},\n    \\qquad\n    M_{\\mathrm{eff},i}^{(d)}\n    \\;=\\;\n    \\bigl(\\Gamma_{i}^{(d)}\\bigr)^{2}\\,\n    \\phi_{i}^{T}\\, M\\, \\phi_{i},\\end{align}\n\nwith $\\phi_{i}$ the i-th mode shape (``mode_shapes[:, i]``),\n$M$ the consistent mass matrix, and $r^{(d)}$ the\n**influence vector** for direction $d$ \u2014 a vector of ones\non every DOF in direction $d$ and zeros elsewhere.  When\nmode shapes are mass-normalised ($\\phi^{T} M \\phi = I$),\nthe formulas simplify to\n$\\Gamma_{i}^{(d)} = \\phi_{i}^{T} M r^{(d)}$ and\n$M_{\\mathrm{eff},i}^{(d)} = \\Gamma_{i}^{(d) 2}$.\n\nTotal-mass property:\n\n\\begin{align}\\sum_{i=1}^{N_{\\mathrm{dof}}} M_{\\mathrm{eff},i}^{(d)}\n    \\;=\\;\n    M_{\\mathrm{total}}\\quad\n    \\text{(every direction; sum over all modes)}.\\end{align}\n\nSo the running cumulative sum of $M_{\\mathrm{eff},i}$\ngives a clean diagnostic \u2014 the first $k$ modes capture\n$\\sum_{i \\le k} M_{\\mathrm{eff},i}^{(d)} / M_{\\mathrm{total}}$\nof the structure's response in direction $d$.\n\n## Implementation\n\nA slender HEX8-EAS prismatic cantilever (the\n:class:`~femorph_solver.validation.problems.CantileverHigherModes`\ngeometry).  Solve 12 modes, get $M$ from\n:meth:`Model.mass_matrix`, build the X / Y / Z influence vectors,\nand compute $\\Gamma$ + $M_{\\mathrm{eff}}$ for every\nmode in every direction.  Print a participation-factor table\nsorted by mode number plus the cumulative-mass curve, and render\nthe cumulative bar chart.\n\n## References\n\n* Bathe, K.-J. (2014) *Finite Element Procedures*, 2nd ed.,\n  Prentice Hall, \u00a79.3.4 (mode-superposition + participation\n  factors).\n* Chopra, A. K. (2017) *Dynamics of Structures*, 5th ed.,\n  Pearson, \u00a713.1 \u2014 modal effective mass.\n* Cook, R. D., Malkus, D. S., Plesha, M. E., Witt, R. J.\n  (2002) *Concepts and Applications of Finite Element\n  Analysis*, 4th ed., Wiley, \u00a711.6 \u2014 base excitation.\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\n\nimport femorph_solver\nfrom femorph_solver import ELEMENTS"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Build a slender cantilever beam (HEX8 EAS)\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "E = 2.0e11\nNU = 0.30\nRHO = 7850.0\nL = 4.0\nWIDTH = 0.05\nHEIGHT = 0.05\n\nNX, NY, NZ = 60, 3, 3\nxs = np.linspace(0.0, L, NX + 1)\nys = np.linspace(0.0, WIDTH, NY + 1)\nzs = np.linspace(0.0, HEIGHT, NZ + 1)\ngrid = pv.StructuredGrid(*np.meshgrid(xs, ys, zs, indexing=\"ij\")).cast_to_unstructured_grid()\n\nm = femorph_solver.Model.from_grid(grid)\nm.assign(\n    ELEMENTS.HEX8(integration=\"enhanced_strain\"),\n    material={\"EX\": E, \"PRXY\": NU, \"DENS\": RHO},\n)\n\n# Full clamp at x = 0.\npts = np.asarray(m.grid.points)\nclamped = np.where(pts[:, 0] < 1e-9)[0]\nm.fix(nodes=(clamped + 1).tolist(), dof=\"ALL\")\n\n# Total mass for the running-sum normalization.\ntotal_mass = RHO * L * WIDTH * HEIGHT\nprint(\"Cantilever beam \u2014 modal participation factors\")\nprint(f\"  L = {L} m, cross = {WIDTH} \u00d7 {HEIGHT} m, \u03c1 = {RHO} kg/m\u00b3\")\nprint(f\"  total mass M_total = \u03c1\u00b7L\u00b7A = {total_mass:.3f} kg\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Modal solve + grab the consistent mass matrix\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "N_MODES = 12\nres = m.solve_modal(n_modes=N_MODES)\nfreqs = np.asarray(res.frequency, dtype=np.float64)\nshapes_dof = np.asarray(res.mode_shapes, dtype=np.float64)\n# ``mode_shapes`` is shaped ``(n_dof, n_modes)``.  Each column is a\n# mass-normalised eigenvector \u03c6_i with \u03c6_i^T \u00b7 M \u00b7 \u03c6_i = 1.\n\nM = m.mass_matrix()\nn_dof = M.shape[0]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Build the X / Y / Z influence vectors\n\nr^(d) is 1 on every DOF aligned with direction d, 0 elsewhere.\nDOFs are interleaved (UX, UY, UZ) per node \u2014 read them off the\nDOF map.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "dof_map = m.dof_map()\nr_x = np.zeros(n_dof)\nr_y = np.zeros(n_dof)\nr_z = np.zeros(n_dof)\nfor row, (_node, dof_idx) in enumerate(dof_map.tolist()):\n    if dof_idx == 0:  # UX\n        r_x[row] = 1.0\n    elif dof_idx == 1:  # UY\n        r_y[row] = 1.0\n    elif dof_idx == 2:  # UZ\n        r_z[row] = 1.0"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Compute \u0393 and M_eff for every mode \u00d7 every direction\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "directions = ((\"X\", r_x), (\"Y\", r_y), (\"Z\", r_z))\ngamma = np.zeros((N_MODES, 3), dtype=np.float64)\nm_eff = np.zeros((N_MODES, 3), dtype=np.float64)\nfor i in range(N_MODES):\n    phi = shapes_dof[:, i]\n    # Mass-normalised: \u03c6^T M \u03c6 = 1, so \u0393 = \u03c6^T M r and M_eff = \u0393\u00b2\n    M_phi = M @ phi\n    modal_mass_i = float(phi @ M_phi)  # \u2248 1 for normalised modes\n    for d_idx, (_, r) in enumerate(directions):\n        L_i = float(phi @ M @ r)  # = modal-load coefficient\n        gamma_i = L_i / modal_mass_i\n        gamma[i, d_idx] = gamma_i\n        m_eff[i, d_idx] = gamma_i**2 * modal_mass_i\n\ncum_eff = np.cumsum(m_eff, axis=0)\ncum_pct = cum_eff / total_mass * 100.0\n\nprint()\nprint(\n    f\"{'mode':>4}  {'f [Hz]':>10}  {'\u0393_X':>10}  {'\u0393_Y':>10}  {'\u0393_Z':>10}  \"\n    f\"{'M_eff,X [kg]':>12}  {'M_eff,Y':>9}  {'M_eff,Z':>9}  \"\n    f\"{'cum_X %':>9}  {'cum_Y %':>9}  {'cum_Z %':>9}\"\n)\nprint(\"  \" + \"-\" * 130)\nfor i in range(N_MODES):\n    print(\n        f\"{i + 1:>4}  {freqs[i]:>10.3f}  \"\n        f\"{gamma[i, 0]:>+10.4f}  {gamma[i, 1]:>+10.4f}  {gamma[i, 2]:>+10.4f}  \"\n        f\"{m_eff[i, 0]:>11.3f}   {m_eff[i, 1]:>8.3f}  {m_eff[i, 2]:>8.3f}  \"\n        f\"{cum_pct[i, 0]:>8.2f}%  {cum_pct[i, 1]:>8.2f}%  {cum_pct[i, 2]:>8.2f}%\"\n    )\n\nprint()\nprint(\n    f\"  total M_eff after {N_MODES} modes: \"\n    f\"X = {cum_pct[-1, 0]:.1f} %, \"\n    f\"Y = {cum_pct[-1, 1]:.1f} %, \"\n    f\"Z = {cum_pct[-1, 2]:.1f} %  (of M_total = {total_mass:.3f} kg)\"\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Verify the partition-of-unity property\n\nSanity check: a sufficient mode count must capture nearly the\nfull mass in every direction.  Twelve modes on a slender beam\nusually clear 90 % in the bending directions.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "assert cum_pct[-1, 1] > 80.0, (\n    f\"Y-direction effective mass {cum_pct[-1, 1]:.1f} % below 80 % \u2014 \"\n    \"increase n_modes or check the mass matrix.\"\n)\nassert cum_pct[-1, 2] > 80.0, f\"Z-direction effective mass {cum_pct[-1, 2]:.1f} % below 80 %.\""
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Render the cumulative effective-mass curve\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt  # noqa: E402\n\nfig, ax = plt.subplots(1, 1, figsize=(7.0, 4.0))\nmode_idx = np.arange(1, N_MODES + 1)\nax.plot(mode_idx, cum_pct[:, 0], \"o-\", color=\"#1f77b4\", lw=2, label=\"cum. M_eff,X\")\nax.plot(mode_idx, cum_pct[:, 1], \"s-\", color=\"#d62728\", lw=2, label=\"cum. M_eff,Y\")\nax.plot(mode_idx, cum_pct[:, 2], \"^-\", color=\"#2ca02c\", lw=2, label=\"cum. M_eff,Z\")\nax.axhline(90.0, color=\"black\", ls=\"--\", lw=1.0, label=\"90 % design target\")\nax.set_xlabel(\"number of modes included\")\nax.set_ylabel(\"cumulative effective mass  [% of M_total]\")\nax.set_title(\"Cantilever \u2014 modal effective mass per direction\", fontsize=11)\nax.set_xticks(mode_idx)\nax.set_ylim(0.0, 105.0)\nax.legend(loc=\"lower right\", fontsize=9)\nax.grid(True, ls=\":\", alpha=0.5)\nfig.tight_layout()\nfig.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Take-aways\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "print()\nprint(\"Take-aways:\")\nprint(\n    \"  \u2022 Modal participation factor \u0393_i = (\u03c6_i^T M r) / (\u03c6_i^T M \u03c6_i) \"\n    \"answers 'how strongly does mode i respond to base excitation in direction d?'.\"\n)\nprint(\n    \"  \u2022 Effective modal mass M_eff,i = \u0393_i\u00b2 \u00b7 (\u03c6_i^T M \u03c6_i) sums to the \"\n    \"total mass of the structure when summed over ALL modes.\"\n)\nprint(\n    \"  \u2022 The first transverse-bending mode dominates one direction (~60-70 %) \"\n    \"and barely contributes to the others \u2014 a single number tells you which \"\n    \"modes can be skipped per excitation axis.\"\n)\nprint(\n    \"  \u2022 Design codes (ASCE 7, Eurocode 8) usually require \u2265 90 % cumulative \"\n    \"effective mass per direction \u2014 read directly off the table above to \"\n    \"decide n_modes for a response-spectrum analysis.\"\n)"
      ]
    }
  ],
  "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
}