{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Modal analysis\n",
    "\n",
    "***\n",
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## General formulation of the problem\n",
    "\n",
    "---\n",
    "\n",
    "The goal of modal analysis is to determine the natural mode shapes and frequencies of an object or structure during free vibration. Let us again consider the problem of elastodynamics (neglecting the volumetric force)\n",
    "\\begin{equation}\n",
    "    \\rho \\ddot{\\mathbf{u}}\n",
    "    =\n",
    "    \\operatorname{div} \\sigma,\n",
    "    \\quad\n",
    "    \\text{in $\\Omega$}\n",
    "\\end{equation}\n",
    "where $\\boldsymbol{\\sigma}$ denotes the stress tensor which is, for an isotropic material, given by\n",
    "\\begin{equation}\n",
    "    \\boldsymbol{\\sigma}=\\lambda\\mathrm{tr}(\\boldsymbol{\\varepsilon})\\mathbf{I} + 2\\mu\\boldsymbol{\\varepsilon}, \\end{equation}\n",
    "where the strain tensor $\\boldsymbol{\\varepsilon}$ takes the form\n",
    "\\begin{equation}\n",
    "    \\boldsymbol{\\varepsilon}=\\frac{1}{2}\\left(\\nabla\\mathbf{u}+ \\left( \\nabla\\mathbf{u} \\right)^{\\top} \\right),\n",
    "\\end{equation}\n",
    "while $\\lambda, \\mu$ are Lame's constants which can be expressed in terms of the Young's modulus $E$ and Poisson's ratio $\\nu$ as\n",
    "\\begin{equation}\n",
    "    \\lambda = \\frac{E \\nu}{(1 + \\nu) (1 - 2 \\nu)},\n",
    "    \\qquad\n",
    "    \\mu = \\frac{E}{2(1 + \\nu)}.\n",
    "\\end{equation}\n",
    "\n",
    "Testing the governing equation by $\\delta\\mathbf{u}$, integrating over the spatial domain $\\Omega$, and integrating by parts yields the weak formulation of our problem. That is, find $\\mathbf{u}$ such that\n",
    "\\begin{equation}\n",
    "    \\int_\\Omega\n",
    "        \\boldsymbol{\\sigma}:\\delta\\boldsymbol{\\varepsilon}\n",
    "    \\,\\mathrm{dV}\n",
    "    + \n",
    "    \\int_\\Omega\n",
    "        \\rho\\ddot{\\mathbf{u}}\\cdot\\delta\\mathbf{u}\n",
    "    \\,\\mathrm{dV}\n",
    "     = \n",
    "     0, \n",
    "     \\quad\n",
    "     \\forall \\delta\\mathbf{u},\n",
    "\\end{equation}\n",
    "where $\\delta\\boldsymbol{\\varepsilon}$ is the symmetric gradient of $\\delta\\mathbf{u}$. To perform the modal analysis we first assume that the displacement field takes the following form (referred to as *normal mode*)\n",
    "\\begin{equation}\n",
    "    \\mathbf{u}(\\mathbf{x}, t)=\\mathbf{u}_A(\\mathbf{x})\\mathrm{e}^{\\mathrm{i}\\omega t},\n",
    "\\end{equation}\n",
    "where $\\mathbf{u}_A$ denotes the amplitude of the oscillating mode and $\\omega$ denotes the frequency of the mode. As a consequence we get\n",
    "\\begin{equation}\n",
    "    \\ddot{\\mathbf{u}}(\\mathbf{x}, t)\n",
    "    =\n",
    "    -\\omega^2\\mathbf{u}_A(\\mathbf{x})\\mathrm{e}^{\\mathrm{i}\\omega t}\n",
    "    =\n",
    "    -\\omega^2 \\mathbf{u}(\\mathbf{x}, t),\n",
    "\\end{equation}\n",
    "and the weak formulation thus reduces to the following problem\n",
    "\\begin{equation}\n",
    "    \\int_\\Omega\n",
    "    \\boldsymbol{\\sigma}:\\delta\\boldsymbol{\\varepsilon}\n",
    "    \\, \\mathrm{dV}\n",
    "    - \n",
    "    \\omega^2\n",
    "    \\int_\\Omega\n",
    "        \\rho\\mathbf{u}\\cdot\\delta\\mathbf{u}\n",
    "    \\,\\mathrm{dV}\n",
    "    = 0, \n",
    "    \\quad\n",
    "    \\forall \\delta\\mathbf{u}.\n",
    "\\end{equation}\n",
    "This formulation represents an eigenvalue problem which cannot be solved by the default FEniCS solver. The problem needs to be reformulated into an algebraic (generalized) eigenvalue problem\n",
    "\\begin{equation}\n",
    "    \\mathbf{A}\\mathbf{u} - \\omega^2\\mathbf{M}\\mathbf{u} = \\mathbf{0},\n",
    "\\end{equation}\n",
    "where $\\mathbf{A}$, $\\mathbf{M}$ are matrices assembled from the governing equations. This reformulation can be done easily in FEniCS, as well as solving the algebraic problem itself."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Implementation\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As our model example we consider a 3D plate whose mesh is generated by the `BoxMesh` function. We do not prescribe any Dirichlet boundary condition. Otherwise, the first blocks of code are practically the same as in the previous examples so we do not comment them further."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "import fenics as fe\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\n",
    "# --------------------\n",
    "# Functions and classes\n",
    "# --------------------\n",
    "# Strain function\n",
    "def epsilon(u):\n",
    "    return 0.5*(fe.nabla_grad(u) + fe.nabla_grad(u).T)\n",
    "\n",
    "# Stress function\n",
    "def sigma(u):\n",
    "    return lmbda*fe.div(u)*fe.Identity(3) + 2*mu*epsilon(u)\n",
    "\n",
    "# --------------------\n",
    "# Parameters\n",
    "# --------------------\n",
    "# Young modulus, poisson number and density\n",
    "E, nu = 70.0E9, 0.23\n",
    "rho = 2500.0\n",
    "\n",
    "# Lame's constants\n",
    "mu = E/2./(1+nu)\n",
    "lmbda = E*nu/(1+nu)/(1-2*nu)\n",
    "\n",
    "l_x, l_y, l_z = 1.0, 1.0, 0.01  # Domain dimensions\n",
    "n_x, n_y, n_z = 20, 20, 2  # Number of elements\n",
    "\n",
    "# --------------------\n",
    "# Geometry\n",
    "# --------------------\n",
    "mesh = fe.BoxMesh(fe.Point(0.0, 0.0, 0.0), fe.Point(l_x, l_y, l_z), n_x, n_y, n_z)\n",
    "\n",
    "# --------------------\n",
    "# Function spaces\n",
    "# --------------------\n",
    "V = fe.VectorFunctionSpace(mesh, \"CG\", 2)\n",
    "u_tr = fe.TrialFunction(V)\n",
    "u_test = fe.TestFunction(V)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we define the weak formulation of our problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "# --------------------\n",
    "# Forms & matrices\n",
    "# --------------------\n",
    "a_form = fe.inner(sigma(u_tr), epsilon(u_test))*fe.dx\n",
    "m_form = rho*fe.inner(u_tr, u_test)*fe.dx"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we create two `PETScMatrix()` objects that will hold the matrices assembled from the weak formulation. (PETSc is a widely used library for sparse matrix computations.) By calling `fe.assemble()`, FEniCS performs integration and localization and outputs the result as a matrix or a vector. The last two lines assemble the standard stiffness matrix `A` and mass matrix `M`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = fe.PETScMatrix()\n",
    "M = fe.PETScMatrix()\n",
    "A = fe.assemble(a_form, tensor=A)\n",
    "M = fe.assemble(m_form, tensor=M)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we want apply Dirichlet boundary condition on matrices `A` and `M`, we call `.zero` method of `DirichletBC` object. For example for `bc` object as `bc.zero(A)` and `bc.zero(M)`.\n",
    "\n",
    "The first step for solving the algebraic eigenvalue problem lies in creating an instance of the `SLEPcEigenSolver` constructor with two arguments: the stiffness matrix `A` and the mass matrix `M`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "# --------------------\n",
    "# Eigen-solver\n",
    "# --------------------\n",
    "fe.PETScOptions.set(\"eps_view\")\n",
    "eigensolver = fe.SLEPcEigenSolver(A, M)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Various parameters of the instance `eigensolver` are then set to control the solver."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "eigensolver.parameters[\"problem_type\"] = \"gen_hermitian\"\n",
    "eigensolver.parameters[\"spectrum\"] = \"smallest real\"\n",
    "eigensolver.parameters[\"spectral_transform\"] = \"shift-and-invert\"\n",
    "eigensolver.parameters[\"spectral_shift\"] = 100.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The possible choices for these parameters are summarized below.\n",
    "\n",
    "| \"problem_type\" | Description |\n",
    "| -- | -- |\n",
    "| “hermitian” | Hermitian |\n",
    "| “non_hermitian” | Non-Hermitian |\n",
    "| “gen_hermitian” | Generalized Hermitian |\n",
    "| “gen_non_hermitian” | Generalized Non-Hermitian |\n",
    "| “pos_gen_non_hermitian” | Generalized Non-Hermitian with positive semidefinite M |\n",
    "\n",
    "| \"spectrum\" | Description |\n",
    "| -- | -- |\n",
    "| “largest magnitude” | eigenvalues with largest magnitude |\n",
    "| “smallest magnitude” | eigenvalues with smallest magnitude |\n",
    "| “largest real” | eigenvalues with largest double part |\n",
    "| “smallest real” | eigenvalues with smallest double part |\n",
    "| “largest imaginary” | eigenvalues with largest imaginary part |\n",
    "| “smallest imaginary” | eigenvalues with smallest imaginary part |\n",
    "| “target magnitude” | eigenvalues closest to target in magnitude (versions >= 3.1) |\n",
    "| “target real” | eigenvalues closest to target in real part (versions >= 3.1) |\n",
    "| “target imaginary” | eigenvalues closest to target in imaginary part (versions >= 3.1) |\n",
    "\n",
    "| \"spectral_transform\" | Description |\n",
    "| -- | -- |\n",
    "| “shift-and-invert” | A shift-and-invert transform |\n",
    "| \" \" | No spectral transform |\n",
    "\n",
    "More information about the solver's parameters can be found in the Dolfin [documentation](https://fenicsproject.org/docs/dolfin/1.4.0/python/programmers-reference/cpp/la/SLEPcEigenSolver.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The eigenvalue problem is then solved by the `solve()` method which takes one parameter: the number of eigenvalues to be calculated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "N_eig = 12   # number of eigenvalues\n",
    "eigensolver.solve(N_eig)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The solution is contained in the `eigensolver` object. We shall extract it and save it to an `.xdmf` file. For this purpose we create an instance of the `XDMFFile()` constructor and set some of its parameters. If true, the parameter `flush_output` allows one to inspect the results while the process is still running. Similarly, if true, the parameter `functions_share_mesh` optimizes the saving process by considering that all functions have the same mesh."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "# --------------------\n",
    "# Post-process\n",
    "# --------------------\n",
    "# Set up file for exporting results\n",
    "file_results = fe.XDMFFile(\"MA.xdmf\")\n",
    "file_results.parameters[\"flush_output\"] = True\n",
    "file_results.parameters[\"functions_share_mesh\"] = True"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, for each eigenvalue we calculate the corresponding eigenfrequency and save the corresponding eigenvector as a FeniCS `Function`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Eigenfrequency:  0.00814 [Hz]\n",
      "Eigenfrequency:  0.00835 [Hz]\n",
      "Eigenfrequency:  0.00857 [Hz]\n",
      "Eigenfrequency:  0.01140 [Hz]\n",
      "Eigenfrequency:  0.01184 [Hz]\n",
      "Eigenfrequency:  0.01208 [Hz]\n",
      "Eigenfrequency: 35.16121 [Hz]\n",
      "Eigenfrequency: 50.83763 [Hz]\n",
      "Eigenfrequency: 59.78057 [Hz]\n",
      "Eigenfrequency: 89.80385 [Hz]\n",
      "Eigenfrequency: 90.36883 [Hz]\n",
      "Eigenfrequency: 153.76596 [Hz]\n"
     ]
    }
   ],
   "source": [
    "# Eigenfrequencies\n",
    "for i in range(0, N_eig):\n",
    "    # Get i-th eigenvalue and eigenvector\n",
    "    # r - real part of eigenvalue\n",
    "    # c - imaginary part of eigenvalue\n",
    "    # rx - real part of eigenvector\n",
    "    # cx - imaginary part of eigenvector\n",
    "    r, c, rx, cx = eigensolver.get_eigenpair(i)\n",
    "    \n",
    "    # Calculation of eigenfrequency from real part of eigenvalue\n",
    "    freq_3D = fe.sqrt(r)/2/fe.pi\n",
    "    print(\"Eigenfrequency: {0:8.5f} [Hz]\".format(freq_3D))\n",
    "    \n",
    "    # Initialize function and assign eigenvector\n",
    "    eigenmode = fe.Function(V, name=\"Eigenvector \" + str(i))\n",
    "    eigenmode.vector()[:] = rx\n",
    "    \n",
    "    # Write i-th eigenfunction to xdmf file\n",
    "    file_results.write(eigenmode, i)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here are two examples of the computed eigenmodes which are visualized using the external tool [ParaView](https://www.paraview.org/)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First eigenmode. | Fifth eigenmode.\n",
    "-|-  \n",
    "![alt](figures/mode1.png) | ![alt](figures/mode5.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Complete code\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Eigenfrequency:  0.00814 [Hz]\n",
      "Eigenfrequency:  0.00835 [Hz]\n",
      "Eigenfrequency:  0.00857 [Hz]\n",
      "Eigenfrequency:  0.01140 [Hz]\n",
      "Eigenfrequency:  0.01184 [Hz]\n",
      "Eigenfrequency:  0.01208 [Hz]\n",
      "Eigenfrequency: 35.16121 [Hz]\n",
      "Eigenfrequency: 50.83763 [Hz]\n",
      "Eigenfrequency: 59.78057 [Hz]\n",
      "Eigenfrequency: 89.80385 [Hz]\n",
      "Eigenfrequency: 90.36883 [Hz]\n",
      "Eigenfrequency: 153.76596 [Hz]\n"
     ]
    }
   ],
   "source": [
    "import fenics as fe\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\n",
    "# --------------------\n",
    "# Functions and classes\n",
    "# --------------------\n",
    "# Strain function\n",
    "def epsilon(u):\n",
    "    return 0.5*(fe.nabla_grad(u) + fe.nabla_grad(u).T)\n",
    "\n",
    "# Stress function\n",
    "def sigma(u):\n",
    "    return lmbda*fe.div(u)*fe.Identity(3) + 2*mu*epsilon(u)\n",
    "\n",
    "# --------------------\n",
    "# Parameters\n",
    "# --------------------\n",
    "# Young modulus, poisson number and density\n",
    "E, nu = 70.0E9, 0.23\n",
    "rho = 2500.0\n",
    "\n",
    "# Lame's constants\n",
    "mu = E/2./(1+nu)\n",
    "lmbda = E*nu/(1+nu)/(1-2*nu)\n",
    "\n",
    "l_x, l_y, l_z = 1.0, 1.0, 0.01  # Domain dimensions\n",
    "n_x, n_y, n_z = 20, 20, 2  # Number of elements\n",
    "\n",
    "# --------------------\n",
    "# Geometry\n",
    "# --------------------\n",
    "mesh = fe.BoxMesh(fe.Point(0.0, 0.0, 0.0), fe.Point(l_x, l_y, l_z), n_x, n_y, n_z)\n",
    "\n",
    "# --------------------\n",
    "# Function spaces\n",
    "# --------------------\n",
    "V = fe.VectorFunctionSpace(mesh, \"CG\", 2)\n",
    "u_tr = fe.TrialFunction(V)\n",
    "u_test = fe.TestFunction(V)\n",
    "\n",
    "# --------------------\n",
    "# Forms & matrices\n",
    "# --------------------\n",
    "a_form = fe.inner(sigma(u_tr), epsilon(u_test))*fe.dx\n",
    "m_form = rho*fe.inner(u_tr, u_test)*fe.dx\n",
    "\n",
    "A = fe.PETScMatrix()\n",
    "M = fe.PETScMatrix()\n",
    "A = fe.assemble(a_form, tensor=A)\n",
    "M = fe.assemble(m_form, tensor=M)\n",
    "\n",
    "# --------------------\n",
    "# Eigen-solver\n",
    "# --------------------\n",
    "eigensolver = fe.SLEPcEigenSolver(A, M)\n",
    "eigensolver.parameters['problem_type'] = 'gen_hermitian'\n",
    "eigensolver.parameters[\"spectrum\"] = \"smallest real\"\n",
    "eigensolver.parameters['spectral_transform'] = 'shift-and-invert'\n",
    "eigensolver.parameters['spectral_shift'] = 100.0\n",
    "\n",
    "N_eig = 12   # number of eigenvalues\n",
    "eigensolver.solve(N_eig)\n",
    "\n",
    "# --------------------\n",
    "# Post-process\n",
    "# --------------------\n",
    "# Set up file for exporting results\n",
    "file_results = fe.XDMFFile(\"MA.xdmf\")\n",
    "file_results.parameters[\"flush_output\"] = True\n",
    "file_results.parameters[\"functions_share_mesh\"] = True\n",
    "\n",
    "# Eigenfrequencies\n",
    "for i in range(0, N_eig):\n",
    "    # Get i-th eigenvalue and eigenvector\n",
    "    # r - real part of eigenvalue\n",
    "    # c - imaginary part of eigenvalue\n",
    "    # rx - real part of eigenvector\n",
    "    # cx - imaginary part of eigenvector\n",
    "    r, c, rx, cx = eigensolver.get_eigenpair(i)\n",
    "    \n",
    "    # Calculation of eigenfrequency from real part of eigenvalue\n",
    "    freq_3D = fe.sqrt(r)/2/fe.pi\n",
    "    print(\"Eigenfrequency: {0:8.5f} [Hz]\".format(freq_3D))\n",
    "    \n",
    "    # Initialize function and assign eigenvector\n",
    "    eigenmode = fe.Function(V, name=\"Eigenvector \" + str(i))\n",
    "    eigenmode.vector()[:] = rx\n",
    "    \n",
    "    # Write i-th eigenfunction to xdmf file\n",
    "    file_results.write(eigenmode, i)\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.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}