{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Modal analysis of an elastic structure\n",
    "\n",
    "This program performs a dynamic modal analysis of an elastic cantilever beam\n",
    "represented by a 3D solid continuum. The eigenmodes are computed using the\n",
    "**SLEPcEigensolver** and compared against an analytical solution of beam theory. We also discuss the computation of modal participation factors.\n",
    "\n",
    "\n",
    "The first four eigenmodes of this demo will look as follows:"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    ".. image:: vibration_modes.gif\n",
    "   :scale: 80%"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first two fundamental modes are on top with bending along the weak axis (left) and along\n",
    "the strong axis (right), the next two modes are at the bottom.\n",
    "\n",
    "## Implementation\n",
    "\n",
    "After importing the relevant modules, the geometry of a beam of length $L=20$\n",
    "and rectangular section of size $B\\times H$ with $B=0.5, H=1$ is first defined:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from dolfin import *\n",
    "import numpy as np\n",
    "\n",
    "L, B, H = 20., 0.5, 1.\n",
    "\n",
    "Nx = 200\n",
    "Ny = int(B/L*Nx)+1\n",
    "Nz = int(H/L*Nx)+1\n",
    "\n",
    "mesh = BoxMesh(Point(0.,0.,0.),Point(L,B,H), Nx, Ny, Nz)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Material parameters and elastic constitutive relations are classical (here we take $\\nu=0$) and we also introduce the material density $\\rho$ for later definition of the mass matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "E, nu = Constant(1e5), Constant(0.)\n",
    "rho = Constant(1e-3)\n",
    "\n",
    "# Lame coefficient for constitutive relation\n",
    "mu = E/2./(1+nu)\n",
    "lmbda = E*nu/(1+nu)/(1-2*nu)\n",
    "\n",
    "def eps(v):\n",
    "    return sym(grad(v))\n",
    "def sigma(v):\n",
    "    dim = v.geometric_dimension()\n",
    "    return 2.0*mu*eps(v) + lmbda*tr(eps(v))*Identity(dim)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Standard `FunctionSpace` is defined and boundary conditions correspond to a fully clamped support at $x=0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "V = VectorFunctionSpace(mesh, 'Lagrange', degree=1)\n",
    "u_ = TrialFunction(V)\n",
    "du = TestFunction(V)\n",
    "\n",
    "\n",
    "def left(x, on_boundary):\n",
    "    return near(x[0],0.)\n",
    "\n",
    "bc = DirichletBC(V, Constant((0.,0.,0.)), left)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The system stiffness matrix $[K]$ and mass matrix $[M]$ are respectively obtained from assembling the corresponding variational forms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<dolfin.cpp.la.PETScMatrix at 0x7fa5da8c0a40>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "k_form = inner(sigma(du),eps(u_))*dx\n",
    "l_form = Constant(1.)*u_[0]*dx\n",
    "K = PETScMatrix()\n",
    "b = PETScVector()\n",
    "assemble_system(k_form, l_form, bc, A_tensor=K, b_tensor=b)\n",
    "\n",
    "m_form = rho*dot(du,u_)*dx\n",
    "M = PETScMatrix()\n",
    "assemble(m_form, tensor=M)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Matrices $[K]$ and $[M]$ are first defined as `PETScMatrix` and forms are assembled into it to ensure that they have the right type.\n",
    "Note that boundary conditions have been applied to the stiffness matrix using `assemble_system` so as to preserve symmetry (a dummy `l_form` and right-hand side vector have been introduced to call this function).\n",
    "\n",
    "Modal dynamic analysis consists in solving the following generalized eigenvalue problem $[K]\\{U\\}=\\lambda[M]\\{U\\}$ where the eigenvalue is related to the eigenfrequency $\\lambda=\\omega^2$. This problem can be solved using the `SLEPcEigenSolver`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "eigensolver = SLEPcEigenSolver(K, M)\n",
    "eigensolver.parameters['problem_type'] = 'gen_hermitian'\n",
    "eigensolver.parameters['spectral_transform'] = 'shift-and-invert'\n",
    "eigensolver.parameters['spectral_shift'] = 0."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The problem type is specified to be a generalized eigenvalue problem with Hermitian matrices. By default, SLEPc computes the largest eigenvalues. Here we instead look for the smallest eigenvalues (they should all be real). A  spectral transform is therefore performed using the keyword `shift-invert` i.e. the original problem is transformed into an equivalent problem with eigenvalues given by $\\dfrac{1}{\\lambda - \\sigma}$ instead of $\\lambda$ where $\\sigma$ is the value of the spectral shift. It is therefore much easier to compute eigenvalues close to $\\sigma$ i.e. close to $\\sigma = 0$ in the present case. Eigenvalues are then transformed back by SLEPc to their original value $\\lambda$.\n",
    "\n",
    "\n",
    "We now ask SLEPc to extract the first 6 eigenvalues by calling its solve function and extract the corresponding eigenpair (first two arguments of `get_eigenpair` correspond to the real and complex part of the eigenvalue, the last two to the real and complex part of the eigenvector)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Computing 6 first eigenvalues...\n",
      "Solid FE:  2.13092 [Hz]   Beam theory:  2.01925 [Hz]\n",
      "Solid FE:  4.09493 [Hz]   Beam theory:  4.03850 [Hz]\n",
      "Solid FE: 13.32102 [Hz]   Beam theory: 12.65443 [Hz]\n",
      "Solid FE: 25.41440 [Hz]   Beam theory: 25.30886 [Hz]\n",
      "Solid FE: 37.15137 [Hz]   Beam theory: 35.43277 [Hz]\n",
      "Solid FE: 69.97631 [Hz]   Beam theory: 70.86554 [Hz]\n"
     ]
    }
   ],
   "source": [
    "N_eig = 6   # number of eigenvalues\n",
    "print(\"Computing {} first eigenvalues...\".format(N_eig))\n",
    "eigensolver.solve(N_eig)\n",
    "\n",
    "# Exact solution computation\n",
    "from scipy.optimize import root\n",
    "from math import cos, cosh\n",
    "falpha = lambda x: cos(x)*cosh(x)+1\n",
    "alpha = lambda n: root(falpha, (2*n+1)*pi/2.)['x'][0]\n",
    "\n",
    "# Set up file for exporting results\n",
    "file_results = XDMFFile(\"modal_analysis.xdmf\")\n",
    "file_results.parameters[\"flush_output\"] = True\n",
    "file_results.parameters[\"functions_share_mesh\"] = True\n",
    "\n",
    "eigenmodes = []\n",
    "# Extraction\n",
    "for i in range(N_eig):\n",
    "    # Extract eigenpair\n",
    "    r, c, rx, cx = eigensolver.get_eigenpair(i)\n",
    "\n",
    "    # 3D eigenfrequency\n",
    "    freq_3D = sqrt(r)/2/pi\n",
    "\n",
    "    # Beam eigenfrequency\n",
    "    if i % 2 == 0: # exact solution should correspond to weak axis bending\n",
    "        I_bend = H*B**3/12.\n",
    "    else:          #exact solution should correspond to strong axis bending\n",
    "        I_bend = B*H**3/12.\n",
    "    freq_beam = alpha(i/2)**2*sqrt(float(E)*I_bend/(float(rho)*B*H*L**4))/2/pi\n",
    "\n",
    "    print(\"Solid FE: {:8.5f} [Hz]   Beam theory: {:8.5f} [Hz]\".format(freq_3D, freq_beam))\n",
    "\n",
    "    # Initialize function and assign eigenvector\n",
    "    eigenmode = Function(V,name=\"Eigenvector \"+str(i))\n",
    "    eigenmode.vector()[:] = rx\n",
    "    \n",
    "    eigenmodes.append(eigenmode)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The beam analytical solution is obtained using the eigenfrequencies of a clamped beam in bending given by $\\omega_n = \\alpha_n^2\\sqrt{\\dfrac{EI}{\\rho S L^4}}$ where :math:`S=BH` is the beam section, :math:`I` the bending inertia and $\\alpha_n$ is the solution of the following nonlinear equation:\n",
    "\n",
    "\\begin{equation}\n",
    "\\cos(\\alpha)\\cosh(\\alpha)+1 = 0\n",
    "\\end{equation}\n",
    "\n",
    "the solution of which can be well approximated by $(2n+1)\\pi/2$ for $n\\geq 3$.\n",
    "\n",
    "Since the beam possesses two bending axis, each solution to the previous equation is\n",
    "associated with two frequencies, one with bending along the weak axis ($I=I_{\\text{weak}} = HB^3/12$)\n",
    "and the other along the strong axis ($I=I_{\\text{strong}} = BH^3/12$). Since $I_{\\text{strong}} = 4I_{\\text{weak}}$ for the considered numerical values, the strong axis bending frequency will be twice that corresponding to bending along the weak axis. The solution $\\alpha_n$ are computed using the\n",
    "`scipy.optimize.root` function with initial guess given by $(2n+1)\\pi/2$.\n",
    "\n",
    "With `Nx=400`, we obtain the following comparison between the FE eigenfrequencies and the beam theory eigenfrequencies :\n",
    "\n",
    "\n",
    "| Mode | Solid FE [Hz] | Beam theory [Hz] |\n",
    "| --- | ------ | ------- |\n",
    "| 1 |  2.04991 |  2.01925|\n",
    "| 2 |  4.04854 |  4.03850|\n",
    "| 3 | 12.81504 | 12.65443|\n",
    "| 4 | 25.12717 | 25.30886|\n",
    "| 5 | 35.74168 | 35.43277|\n",
    "| 6 | 66.94816 | 70.86554|\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Modal participation factors\n",
    "\n",
    "In this section we show how to compute modal participation factors for a lateral displacement in the $Y$ direction. Modal participation factors are defined as:\n",
    "\n",
    "\\begin{equation}\n",
    "q_i = \\{\\xi_i\\}[M]\\{U\\}\n",
    "\\end{equation}\n",
    "\n",
    "where $\\{\\xi_i\\}$ is the i-th eigenmode and $\\{U\\}$ is a vector of unit displacement in the considered direction. The corresponding effective mass is given by:\n",
    "\n",
    "\\begin{equation}\n",
    "m_{\\text{eff},i} = \\left(\\dfrac{\\{\\xi_i\\}[M]\\{U\\}}{\\{\\xi_i\\}[M]\\{\\xi_i\\}}\\right)^2 = \\left(\\dfrac{q_i}{m_i}\\right)^2\n",
    "\\end{equation}\n",
    "\n",
    "where $m_i$ is the modal mass which is in general equal to 1 for eigensolvers which adopt the mass matrix normalization convention.\n",
    "\n",
    "With `FEniCS`, the modal participation factor can be easily computed by taking the `action` of the mass form with both the mode and a unit displacement function. Let us now print the corresponding effective mass of the 6 first modes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--------------------------------------------------\n",
      "Mode 1:\n",
      "  Modal participation factor: -7.83e-02\n",
      "  Modal mass: 1.0000\n",
      "  Effective mass: 6.13e-03\n",
      "  Relative contribution: 61.30 %\n",
      "--------------------------------------------------\n",
      "Mode 2:\n",
      "  Modal participation factor: 7.88e-04\n",
      "  Modal mass: 1.0000\n",
      "  Effective mass: 6.21e-07\n",
      "  Relative contribution: 0.01 %\n",
      "--------------------------------------------------\n",
      "Mode 3:\n",
      "  Modal participation factor: 4.34e-02\n",
      "  Modal mass: 1.0000\n",
      "  Effective mass: 1.89e-03\n",
      "  Relative contribution: 18.85 %\n",
      "--------------------------------------------------\n",
      "Mode 4:\n",
      "  Modal participation factor: -4.40e-04\n",
      "  Modal mass: 1.0000\n",
      "  Effective mass: 1.94e-07\n",
      "  Relative contribution: 0.00 %\n",
      "--------------------------------------------------\n",
      "Mode 5:\n",
      "  Modal participation factor: 2.55e-02\n",
      "  Modal mass: 1.0000\n",
      "  Effective mass: 6.49e-04\n",
      "  Relative contribution: 6.49 %\n",
      "--------------------------------------------------\n",
      "Mode 6:\n",
      "  Modal participation factor: -3.71e-04\n",
      "  Modal mass: 1.0000\n",
      "  Effective mass: 1.38e-07\n",
      "  Relative contribution: 0.00 %\n",
      "\n",
      "Total relative mass of the first 6 modes: 86.65 %\n"
     ]
    }
   ],
   "source": [
    "u = Function(V, name=\"Unit displacement\")\n",
    "u.interpolate(Constant((0, 1, 0)))\n",
    "combined_mass = 0\n",
    "for i, xi in enumerate(eigenmodes):\n",
    "    qi = assemble(action(action(m_form, u), xi))\n",
    "    mi = assemble(action(action(m_form, xi), xi))\n",
    "    meff_i = (qi / mi) ** 2\n",
    "    total_mass = assemble(rho * dx(domain=mesh))\n",
    "\n",
    "    print(\"-\" * 50)\n",
    "    print(\"Mode {}:\".format(i + 1))\n",
    "    print(\"  Modal participation factor: {:.2e}\".format(qi))\n",
    "    print(\"  Modal mass: {:.4f}\".format(mi))\n",
    "    print(\"  Effective mass: {:.2e}\".format(meff_i))\n",
    "    print(\"  Relative contribution: {:.2f} %\".format(100 * meff_i / total_mass))\n",
    "\n",
    "    combined_mass += meff_i\n",
    "\n",
    "print(\n",
    "    \"\\nTotal relative mass of the first {} modes: {:.2f} %\".format(\n",
    "        N_eig, 100 * combined_mass / total_mass\n",
    "    )\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that the first, third and fifth mode are those having the larger participation which is consistent with the fact that they correspond to horizontal vibrations of the beam."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Format de la Cellule Texte Brut",
  "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.8.10"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": false,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}