{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Vibration of membranes under tension\n",
    "\n",
    "In the previous lecture, I discussed time-dependent diffusion problem. It was a parabolic PDE. Now, I will discuss vibration in Membranes. This problem, at the simplest level, takes the form:\n",
    "\n",
    "\\begin{align}\n",
    "\\rho {\\partial^2 \\over \\partial t^2} w(x, y, t) &= T \\nabla^2 w(x, y, t) + f(x, y, t)\\;\\;\\;\\;\\mbox{domain (circle of radius 1)}\\\\\n",
    "w(x, y, 0) &= 0\\;\\;\\;\\;\\;\\mbox{ at the boundary the membrane is fixed}\\\\\n",
    "w(x, y, 0) &= w_0(x, y)\\;\\;\\;\\mbox{initial displacement}\\\\\n",
    "{\\partial \\over \\partial t}w(x, y, t)\\left . \\right\\vert_{t = 0} &= v_0(x, y)\\;\\;\\;\\mbox{initial velocity}\n",
    "\\end{align}\n",
    "\n",
    "Note that we require two initial conditions, one each in displacement and velocity. \n",
    "In the problem that we solve below, we will use initial displacement and velocity profile as zero, i.e., $w_0(x, y) = v_0(x, y) = 0$. Note that $\\rho$ is density per unit area and $T$ is the tension in the membrane. The quantity $\\sqrt{{T \\over \\rho}} = c$,which is the wave speed in the membrane. We will hence, for simplicity, replace the main differential equation with a more compact version\n",
    "\n",
    "$$ {\\partial^2 \\over \\partial t^2} w(x, y, t) = c^2 \\nabla^2 w(x, y, t) + f(x, y, t),$$\n",
    "\n",
    "where, we apply oscillatory, loading of the form:\n",
    "$$f(x, y, t) = \\epsilon \\exp \\left(-\\frac{x^2 + y^2}{2\\sigma^2} \\right)\\cos \\omega t.$$\n",
    "The loading is centered about $x = y = 0$ and oscillates with frequency $\\omega$, and $\\epsilon$ is the magnitude of the load. \n",
    "\n",
    "Such PDEs with second order derivative in time are called as hyperbolic and the numerical techniques that are required need to be more sophisticated (look at the book by [Seshu](https://www.amazon.in/Textbook-Finite-Element-Analysis-Seshu/dp/8120323157/ref=sr_1_1?crid=24AHCDSKIUX56&dchild=1&keywords=finite+element+analysis+p+seshu&qid=1585123865&sprefix=Seshu+finite%2Caps%2C268&sr=8-1)). We will follow a simple rule -- extension of the Crank-Nicholson rule, and can be also thought of approximately as [mid-point Euler](https://en.wikipedia.org/wiki/Midpoint_method) scheme, or Trapezoidal rule. As described below, Newmark's scheme (or its variations) are demonstrably better (see the following [demo](https://comet-fenics.readthedocs.io/en/latest/demo/elastodynamics/demo_elastodynamics.py.html) or the book by [K. J. Bathe](https://www.amazon.in/Finite-Element-Procedures-United-States/dp/0133014584/ref=sr_1_1?dchild=1&keywords=finite+element+procedures&qid=1585124041&sr=8-1)). \n",
    "\n",
    "Our simple scheme is as follows. We use the notation $w$ and ${\\partial w \\over \\partial t} = v$ (`wdot` or velocity of the membrane.) With this notation, the scheme is:\n",
    "\\begin{align}\n",
    "{w^{t + \\Delta t} - w^t \\over \\Delta t} &= {1\\over 2}\\left(v^{t+\\Delta t} + v^t \\right),\\\\\n",
    "{v^{t + \\Delta t} - v^t \\over \\Delta t} &= {c^2 \\over 2}\\left(\\nabla^2 w^{t+\\Delta t} + \\nabla^2 w^t \\right) + f(x,y,t+{\\Delta t \\over 2}),\n",
    "\\end{align}\n",
    "where the forcing term is evaluated at $t + \\Delta t/2$. We can manipulate these two equations to provide us with two equations in terms of $w^{t + \\Delta t}$ and $v^{t + \\Delta t}$ as follows:\n",
    "\\begin{align}\n",
    "w^{t + \\Delta t} &= w^t + v^t\\Delta t  + {c^2 \\Delta t^2 \\over 4}\\left(\\underbrace{\\nabla^2 w^{t+\\Delta t}}_{\\ddot{u}^{t+\\Delta t}} + \\underbrace{\\nabla^2 w^t}_{\\ddot{u}^{t}} \\right) + {c^2 \\Delta t^2 \\over 2}f(x,y,t+{\\Delta t \\over 2}),\\\\\n",
    "v^{t + \\Delta t} &= 2 {w^{t + \\Delta t} - w^t \\over \\Delta t} + v^t\n",
    "\\end{align}\n",
    "Thus we can start with the initial conditions for $w^0$ and $v^0$. Solve the above equation using the following weak form (similar to what we did for heat diffusion problem)\n",
    "\n",
    "$$\\int_\\Omega w^{t+\\Delta t} \\delta w d\\Omega - \\int_\\Omega w^t \\delta \\Omega - \\Delta t \\int_\\Omega v^t \\delta w d\\Omega + {c^2 \\Delta t^2 \\over 4} \\int_{\\Omega} (\\nabla w^{t+\\Delta t} + \\nabla w^t)\\cdot \\nabla \\delta w d\\Omega - {c^2 \\Delta t^2 \\over 2} \\int_\\Omega f(x, y, t + {t \\over 2})\\delta w d\\Omega = 0,\n",
    "$$\n",
    "\n",
    "to obtain $w^{\\Delta t}$, which we use to update $v^{t + \\Delta t}$. The solution is thus propagated.  \n",
    "\n",
    "Note that after the finite element analysis you actually get a following equation:\n",
    "$$[M]\\{\\ddot{u}\\} + [K]\\{ u\\} = \\{F(t)\\},$$\n",
    "where:\n",
    "\\begin{align}\n",
    "[M] &= \\int_{\\Omega} \\{N\\}\\langle N \\rangle d\\Omega,\\\\\n",
    "[K] &=\\int_{\\Omega}   c^2 \\{B\\}\\langle B \\rangle d\\Omega,\\\\\n",
    "\\{F\\} &=  \\int_\\Omega \\{N\\} f(x, y, t) d\\Omega.\n",
    "\\end{align}\n",
    "We can explicity get these matrices in FEniCS and then potentially add damping terms, and also create more __efficient__ numerical schemes as shown in the [elastodynamics demo](https://comet-fenics.readthedocs.io/en/latest/demo/elastodynamics/demo_elastodynamics.py.html). However, to demonstrate the concept, we perform a simpler analysis as demonstrated below."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### FEniCS demo\n",
    "Load libraries, define the parameters of the simulations, and create a circular domain centered around the origin and having radius of 1. The domain has radius of 1. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import print_function\n",
    "from dolfin import *\n",
    "import numpy as np\n",
    "import mshr\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "\n",
    "T = 20.0            # final time\n",
    "num_steps = 1000     # number of time steps\n",
    "dt = T / num_steps # time step size\n",
    "c = 1.0 # wave speed in the material\n",
    "PI = np.pi # defining PI as a parameter\n",
    "epsilon = 0.05 # amplitude of the forcing\n",
    "omega = 2.0 # frequency of oscillation\n",
    "sig = 0.05 # width of Gaussian of the forcing\n",
    "\n",
    "\n",
    "Nmesh = 50\n",
    "domain = mshr.Circle(Point(0,0), 1)\n",
    "mesh = mshr.generate_mesh(domain, Nmesh)\n",
    "plot(mesh)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define relevant functions spaces on the mesh and the initial conditions for displacement `w` and velocity `wdot`. Define the boundary conditions. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define functionspace\n",
    "V = FunctionSpace(mesh, 'P', 1)\n",
    "\n",
    "# Define boundary condition\n",
    "w_D = Constant(0.)\n",
    "\n",
    "# Define initial value\n",
    "w_n = interpolate(Constant(0.0), V)\n",
    "wdot_n = interpolate(Constant(0.0), V) # initial velocity is given to be zero\n",
    "\n",
    "# Define boundary conditions\n",
    "def boundary(x, on_boundary):\n",
    "    return on_boundary\n",
    "\n",
    "bc = DirichletBC(V, w_D, boundary)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As in the case of the heat equation, define the variational problem. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define variational problem\n",
    "w = TrialFunction(V)\n",
    "wdot = TrialFunction(V)\n",
    "v = TestFunction(V)\n",
    "\n",
    "# External\n",
    "f_ext = Expression('epsilon*exp(-(x[0]*x[0] + x[1]*x[1])/(2*sig*sig))*cos(omega*t)', \\\n",
    "                   degree = 2, epsilon = epsilon, omega = omega, sig = sig, t = dt/2.)\n",
    "\n",
    "# using Crank-Nicholson scheme\n",
    "F = w*v*dx - w_n*v*dx - dt*wdot_n*v*dx - \\\n",
    "    c*c*dt*dt/2*f_ext*v*dx + \\\n",
    "    c*c*dt*dt/2./2.*(inner(grad(w),grad(v)) + inner(grad(w_n), grad(v)))*dx\n",
    "a, L = lhs(F), rhs(F)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Begin the update for the variational form. Save .pvd file for paraview. Also obtain elastic and kinetic energy. The initial variable `w_n` and `wdot_n` is now updated to carry the new value `w` and `wdot`. The process continues. The interplay is between the rate of oscillatory loading $\\omega$ and the speed of wave $c$, which depends on tension and density of the membrane. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# start the update\n",
    "# Time-stepping\n",
    "file = File(\"./Membrane/Displacement.pvd\", \"compressed\")\n",
    "\n",
    "w = Function(V)\n",
    "wdot = Function(V)\n",
    "t = 0\n",
    "for n in range(num_steps):\n",
    "\n",
    "    # Update current time\n",
    "    t += dt\n",
    "    f_ext.t = t\n",
    "\n",
    "    # Compute solution\n",
    "    solve(a == L, w, bc)\n",
    "    \n",
    "    # update velocity\n",
    "    wdot.vector()[:] = 2*(w.vector() - w_n.vector())/dt - wdot_n.vector()\n",
    "\n",
    "    E_kin = 0.5*c*c*assemble(wdot*wdot*dx)\n",
    "    E_elas = 0.5*assemble(inner(grad(w),grad(w))*dx)\n",
    "    E_tot = E_elas + E_kin\n",
    "    # Plot solution\n",
    "    #plot(w)\n",
    "    #plt.pause(0.5)\n",
    "\n",
    "    # Update previous solution\n",
    "    w_n.assign(w)\n",
    "    wdot_n.assign(wdot)\n",
    "    \n",
    "    if (np.mod(n,10) == 0):\n",
    "        file << (w, t)\n",
    "        print(\"Elastic energy is \", E_elas)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The __complete__ code is given here, and can be run as a python file. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQ0AAAD4CAYAAAD2OrMWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO2deXRV1dn/v09GMkBCQsAggTApiMzI5KqAFVqx1KkvS23VZV1iFV3Y5a9iq0DBLqVV+0q7VH5o9aetw8u71OKqTA5FUQoKCAHKDGFIUgYjYAiSaf/+uOfEk5t77j3DPnufYX/Wysq9Z9jnudNz9vA8z5cYY1AoFAqrpMk2QKFQBAvlNBQKhS2U01AoFLZQTkOhUNhCOQ2FQmGLDNkGOKFLly6svLxcthkKRWjZtGnTScZYSaJ9gXQa5eXl2Lhxo2wzFIrQQkSHzPap4YlCobCFchoKhcIWymkoFApbKKehUChsoZyGQqGwBRenQUQvEdFxItpusp+I6E9EtI+IKohohGHfD4lot7bvYR72KBQK7+DV0/h/AH6YZP/VAPprfzMAPA8ARJQO4Flt/yUAbiaiSzjZpFAoPIBLnAZj7BMiKk9yyLUAXmWxPPz1RFRIRKUAygHsY4wdAAAielM79t887FJ4x7fffovs7GxUVVWhvr4eb7zxBoYOHYqtW7eiU6dOOHPmjOm5GRkZaGpqwpAhQ1BRUYEf//jHAIBBgwYBALKysoS8BoUzRAV3XQjgiOH5UW1bou1jEjVARDMQ66WgZ8+e3lipaIUxBiJCfX099u/fj7fffjvlOVu3bgWApA4DAJqamgAAFRUVAIB33323zf94Jk+ejNLSUvTq1QtEBCKy/DoU/BHlNBJ9yizJ9vYbGVsCYAkAjBo1SlUO4gxjDPX19Vi5ciV2796NxsZG2Sa18v7777fblpmZiZkzZyInJ0f1TAQjymkcBVBmeN4DQDWALJPtCo9pbm7GN998g7/+9a+ora2VbU47CgsLcerUKdP9jY2NeOaZZ9psmzVrFrKzs5GTk+O1eZFGlNN4F8B92pzFGACnGWM1RHQCQH8i6g2gCsBNAG4RZFPkaGlpwY4dOywNNWSTzGGYsWjRotbH/fr1w0033YT09HSeZinAyWkQ0RsAJgLoQkRHAcwDkAkAjLHFAJYDmApgH4B6AHdo+5qI6D4AqwCkA3iJMbaDh02KGC0tLXjzzTexd+9eW+d17doVx48fd3xdIoLM+rP79u3D7373OwBAQUEBbrvtNhQVFUmzJ0xQEAsLjxo1iqks1+ScPHkSzz77bNJj0tPT0dzcLMgifzBx4kRcccUVajI1BUS0iTE2KtG+QKbGKxJz+PBh7NmzB5999pml46PmMABgzZo1WLNmDYqKinDVVVdhwIAByoHYRDmNELBlyxYsW7bM9nkFBQU4ffq0Bxb5n9raWixduhQAkJOTg4ceekiyRcFBDU8CzLp16xIuR9pBD7Tq37+/7XmPMDJ37lzV80Dy4YlyGgGDMYZNmzbhvffeS7jf7gRmx44d8c033/AyL+l1+vTp0xoA5mcKCgowa9asSDsP5TRCwuHDh/Hyyy/bOicnJwfnzp3zyKJwU1xcjHvvvRdpadFLBldOI8AwxlBdXY0XX3xRtimRZerUqRg2bBgyMzNlmyIM5TQCSm1tLf785z+3215cXIyvvvpKiA1WrxWFHs0dd9wRmbwn5TQCxvnz57F06VIcOHDA9rlRjL0QzaxZs1BYWCjbDE9J5jSiN1jzMfpQZOHChSkdhll+RRgdht9CwRctWoS33noLLS0tsk2RgorT8AmNjY14/PHHLR/vl6GAiGGJHx3h9u3bsX37djzwwAMoKCiQbY5QlNPwAR999BHWrl3rup3c3FzU19dzsMg6fnFestAzbaMU36GGJxJpaWnB/PnzuTgMAMIdRjzDhg2Ten2ZLFiwAAcPHpRthhDURKgkNm7caBqgJQo1aeoN8+bNk22Ca9REqI/QexepHEZeXh7X615wwQXttiVzGB07duR6/Sgxf/78UOf0KKchkPPnz2PhwoXttidaCTl79izXa//nP/+xdbyI0HK/0rt3b9dtPPPMMzh8+DAHa/yHchqCOHXqFBYuXJiw9qZXk4nJamcOGDDAk2uGAV5zEy+//HIokwB5iSUlFTwiol8R0RbtbzsRNRNRkbavkoi2afuCPVFhwpYtW9qUohOFcTY/PgR6165dyM3NxcSJEwVb5Q+ys7MTbs/Ly0NZWVnCfU54/fXXMX/+fG7t+QHXE6Ga4NEeAJMRKyD8BYCbGWMJtUuIaBqAXzLGrtSeVwIYxRg7afWaQZoIXbJkCWpqamSbofABc+bMCUzym9cToaOhCR4xxhoA6IJHZtwM4A0O1/U1zc3NmD9/vi8chsgvapSSuuzy2GOP+eL74BYe3yYzIaR2EFEuYvKNbxk2MwCriWiTJoiUECKaQUQbiWjjiRMnOJjtHY2NjbZT2L1EZLhznz59hF0rnuLiYmnXtsqSJUsCv7LCw2lYFjwCMA3AZ4wxo9DG5YyxEYjpuc4koisSncgYW8IYG8UYG1VSUuLOYg9hjOHxxx9HVVWVbFOksHv3bmnXFpX565ZnnnmG++qYSHg4DTMhpETchLihCWOsWvt/HMA7iA13AklzczMWLFgg2wyFD+ncuXOb50899RQ+//xzSda4g4fT+AKa4BERZSHmGNqJchJRAYAJAJYZtuURUUf9MYApALZzsEk4zc3NrTobfqdDhw4A2s8/+C2bNAzoS9tff/11u30rVqwIZCyHa6fBGGsCoAse7QSwlDG2g4h+QUS/MBx6PYDVjDFjv6wbgE+JaCuAzwG8xxhb6dYm0TDGAuMwACA/Px8A2sWMBCmk/JJLLjHd56dCObt27WrzvEuXLrjwwu+m/F5++WXU1dWJNssVKvfEJYyxQA9J0tLSIlsXQhQ9e/bE4cOHk77XP//5z7nGh7hF5Z54REtLi1CHwTsfBRC7shJV9CFIsvD0l156CZWVlYIscodyGg5hjOEvf/mL0GsGecZdgZQraq+88goaGhoEWeMc5TQcsmDBAlRXmy0SKcJAstwdr9p74oknpApnW0E5DQf8/e9/T7o/PgIzKhWdwgbvu/6ZM2csHbdgwQJfT0orp2GT2trapCphmZmZ7eYJ/H7nCAKdOnWSbYJQPv30U9kmmKKchg3q6uoS6pAYSZT6rnCP1bu0aPr27etJu7q6vR9RTsMijDE8/fTTss3ghhcrMX5C1JBw//79pvvMZCas8vHHH+OLL75w1YYXKKdhkcWLF8s2gStBXYnp2rWrpeNEDQkHDx5suo9HcaXly5f7bnirnIYFFi9ebEuJPRVu70Cp0MPEw8bYsWOTfg4y0vK3bduGiy66yNNr+C14UDmNFJw+fRrHjh3j2qYX5f2Mw41vv/2We/t28CqHJVUtisbGRuTm5npy7WSIkC5YudI/2RXKaaRAF8PxOyKHG2aRjd26dUN6erpnKx2HDh1KeYwM7RcRPbsNGzZI17XRUU4jCWGr7ZgIJz9wszvrsWPH0NzcjK+//hqjRiVMWwgVvXr1AiCucvuTTz6JpqYmIddKhnIaJgSlhL/boYDbpUw9yjE+2covCYU8ideCsdLz4Y0f3lflNEz44x//aOt43iHHVpEdOahHTR45ciTFkcFFLyWQSHDKC+LLFk6YMKH18apVqyC73KVyGgl4++23bZ8ThEQjWeg/uqCir3ZZ0TDReyM/+MEPHF0rNze3tWzhoEGDAMTiNYw899xzUm8WonRPJhLRaYP2yVyr54rm22+/xbZt21qfl5eXpzynS5cuHlrEh27dukm7dtCKzMRj586uD2tXrVrl6FrGyc4dO3aYHve3v/3NUfs8cO00NN2TZxErDHwJgJuJKFFZpbWMsWHa3wKb5wqBMYbf//73bbYlq3GgxwWcPBmTbDFWZPIbqZaNg94bCCPDhw9v83z8+PGtjysrK6WF1svQPeF1LneMeqdWtELi80yCXIGcR2+gR48eHCz5jssuu4xre34mkdP+8ssv2zxft25dm+ep8qC8QqTuyTgi2kpEK4hokM1zPaelpQVLlixp81zHTgRnENPgCwsLubRz9OhRLu3ouM27iK8A7mecOO2mpiZs2LDBA2uSI0r3ZDOAXoyxoQD+DEAvSGFZM8VrsaQPP/zQdF98BGeiZc6MjAwA1nMerMyViOLUqVOyTXBMMseQqAK4HzEOOwB74twyIkWF6J4wxs4wxuq0x8sBZBJRFyvnGtrwTCypubm5Xdcv1fHx2A264X1XjipBcQzJiP/uxVcwT8X69et5mpMSIbonRHQBaf12IhqtXfcrK+eK4PHHHxd9SV9E9oWB0tJS2SbYhndMz6pVq4TWcRGle/ITANs1fZM/AbiJxUh4rlub7NDc3OxZRe6hQ4d60q7iO4IoqOxFTM8nn3zCvU0zIq97EoX8EoU/GTt2LNavX4/i4mIuOrSPPvootwxjpXtiguwQbEW00ecizp8/z6U9ESn6QMSdhl0pxe7du3tkiSLK8IqYfe2114RU+Yqs03Ayj6F0TsQyZMgQS8cp4ervEKHSFlmn8dRTT8k2QZGCiooKS8epYeZ3vPrqq55fI5JOo7m52ZOSe2GGRxk9r2tpKmJ4nc4QSafhp3qLvDHrqusRq07hUWpuz549rttQpObFF1/0tP1IOg0/VD/yCrOuenwwWbIcmfz8fGRkZLSrVKUIDl7WE42c0+BZqbuoqKj1cVZWlvRkNauaIEDiHBm95mVdXR2ampoCU/JQ0R5j8iVvIuc04utluKG2thYAMGLECDQ0NEgXtXGrzSKj5qXCG06fPu3Z9zFyTsMIrxyAzZs3c2lHNiUlJUnriMgQI1I4x6tkvkg5je3bt7d5zisHYMSIEVzakc2JEyeSxq8ocWtrFBQUyDYBgHdFeiLlNN59l28CbXl5Obp27WqrpxFk4WWZPwYrldS8xM7q0+nTpy0dl52d7dQcqUTGaZw/f577nbKystL2PIIbJTS9JJysH5DVH4MXeJWJbBU7pQysFpq2k3PiVLVu2bJljs5LRmScRpCrU+noOQpmP6Agzjkk+oFdeeWV3K/DM9Tc2FbPnj3b7dcLTSdTlLeL0yLChw8f5maDTmScxuLFi7m15ZVWqVu8nnPgKa6sd/f1H5hRoe2jjz7idh0dnqHmxraS/SiNUhiy0Ff4eBIZp8ETWaXjZUFEKCws5BowFF99O16hzWqymhfYibcxvo6+fft6YY5reMfbiBJL+ikRVWh/64hoqGFfJRFt00SUPAnV9IPHjydIQwnGGLfhnR4QF99evOShnqwmIyrVjgq8Ma19//79XpjjGrsSo6kQJZZ0EMAExtgQAI8BiA9Xm6SJKHkiNW61JJwXY2kzjEMJN5GksqNQ7WLWXTZqzhhxM3HsFL8kMzoV38rNzfV0ZUaIWBJjbB1jTI80WY9Y1XFh/Otf/7J0nBdjaSu4idyTHYXqNbJXTWRSVVXlaAK3vr6+3coMz++JSLEknTsBrDA8ZwBWE9EmIpphdpLXuid24O3F7XSH/YTbzNlk2L3LZmVlhbIYD68JXJ5DJ1FiSbEDiSYh5jRmGzZfzhgbgdjwZiYRXZHoXKe6J16kwfOq6ajDM4lOJJdc4p3srt2aEA0NDWhubnbUpffralgi4ud+rPLaa69xs0GIWBIAENEQAC8CuJYx1lp6mTFWrf0/DuAdxIY73PCixL1f5f68mFwtLi423We1spZInBSgCcJqmP45mM39xKM7Fy8CAUWJJfUE8DaAWxljewzb84ioo/4YwBQAbRNEXOJFcIvVRCA91VxUxSov4jR4lNZXuMfu56A7Fy/mhFwPShljTUSkCx6lA3hJF0vS9i8GMBdAMYDntNn+Jm2lpBuAd7RtGQBeZ4yFpqyWfgcLc8WqjIwMEJH0ZLacnBzfrHr4laqqKscrMkZCLZbU0NCAJ554QoBFikSUlZW1CdrKz89HXV0dsrOzuc8LKawxb948S8dFVixp69atsk2INPFRnnogVFQdBk8NV5krbqF2GsuXL7d1vBdp6yNHjuTept8w+wJ369Yt4eOowlPDVeaKW6idRiKSJV15EX24adMmy8f6pXiLXcy+wMeOHWvzOGjRq2GER9xH5JyGl1WadfR1f7vr/zLrVYgg2fxZUAPcggaPHKLIOQ0R6HfeIKz/+wW/B7jFz0fYCTD0GjuRsAcOHHB9vdA6DZmTbUEt46YwxzgfkZ6eDtmpDEbsDDnszvMlIrROY+/evdKurfRCwk2iH6mXeTg84REhGlqnsXr1atkmKCyQTOBp+PDhAPiWzfMKOzVEZcIjQjS0TkNFBwYDY2FmY54LEeHLL78EAOzbt0+4XVZJlpsTVkLrNIyev0+fPsjKykJpaalEixSpMOZXGFda/Cz7EMTcHLdR4KF1GkYOHDiAhoYGTzJeFXxIlvKtFx/2K9nZ2cjMzJRa19QObuNlIuE0vIJnde4okKjcv45Zynf//v0dXat3795tBLq9RNfU8WOpgES4DfAKpdMQlYQXHyg2cOBALu2G1RnFlymwEvzmdBXs4MGDnpTvDwNu5/tC6TRkwSsFXkTUqiyMXWOrwW8inWifPn2EXUsWVgv5mBFKpyEj3T8rK0slZVlA/2zshI2LdKI8IibtIjrnyG26gnIanGhoaEB1dbXjGo4i8UMBXmPYuF8Do5LNwfBE/xFbCbzikaNz9OhRV+eLEksiIvqTtr+CiEZYPdcJMn8Ubrt+Ikg1ESa6BqqVwCgZoflelIpMhpXAKx45Om4djyixpKsB9Nf+ZgB43sa5tjHraVhV83aDiGvYwYkDsFoDVSSicol0mcXCwkIh1xOFUTLS7QSxELEk7fmrLMZ6AIVEVGrxXNuY3blErPf7LaaAtwPwc6AVD/TqYrxkKP2CUffE7Y1NlFiS2TGWhZbsiCUFse6pU0TPB8iQSVTwxa1DFCWWZHaMZaElO2JJfpjo8xr9jh+URCmFf3Bbq1SUWJLZMZaEluwSBacR5ju+nxTP/Lqy4wa3NxohYkna89u0VZSxAE4zxmosnmsbEcMTXtGfivYkCvqSdSNw8wMbO3YsR0v44Vb7RJRY0nIAUwHsA1AP4I5k57q1yesue+/evbFz505Pr6FoCy8hZJGsX79etgkJcbt8zaXvxRhbjphjMG5bbHjMAMy0eq5bvNA0NRKv5xFFCgsLQ7fCwJvhw4e31gTxEx07dnR1vooIdYDdnowfx8Vu8zmiMG/kFtEOw+oEp/TgLj/iN30NP65wuM3nCGLxmTCTl5dnWYypR48erq4VSqehiCZ5eXmBrwTvtPCvF+rwZiin4YJk3bx+/foJtEQBxJahg64T6/THb6dGhir3J5FkyUPGsF078Cgxr0iMnXwSImrXjXe7VGmVXr16CbmOU0L7DZWdoj506FDb5xQVFQlLx7aDnsQVZC644AJbqz2MsXYp5FVVVbzNSsihQ4c8bV/VCDVBdvDVli1bbJ9TW1uLyspKfP/73/fAIufoSVxBLkPopGSBE8cfBULrNIIgsGPGhx9+6Gn7TrU6rK64eB0nI4qtW7fKNoE7PD6b0DoNP94VzeYreJYJtFKB2+vl0sbGRk/bVziHR+8ptE7Dj0tvZjPjx44da7fNadSescCKnvg1bNiwNseEvSaGwpwxY8a4biO0TiPo8BCR1oN94udXRGfIuk3FVvCDR0Uy5TRCTDJxZR5YdQYyhithK9fHCx4pDaF2GoMGDZJtglS8LoxrNWyZRy5QfIm6VE5BJdN5R6idxpQpU2SboOBEfO1VWU7hoosuAuDPJERRhNppiAhKkh1EJgO3qdVuKSkpkSZMpavo+TEJMRVXX301l3ZC7S5FhGTz0DnJzMxEU1NTYAoi85ikdUOqwtJhoKioiLsW7ciRI7m04+pXRURFRPQ+Ee3V/rcT2SCiMiL6JxHtJKIdRDTLsO+3RFRFRFu0v6lu7DGxkXeT3GX0GhsbPXEYfisRIAsn74PIHmSi+RkvxKt53UTdtvIwgA8ZY/0BfKg9j6cJwIOMsYEAxgKYGSeI9N+MsWHaH9cKXgAwceJE3k060sJMtZJRWlrq1BxT3H5JwrJU6sQhWxWn5oGo+RleNxG3TuNaAK9oj18BcF38AYyxGsbYZu3xNwB2wkTbxAvGjRsn6lJJOX78eNL9NTU13K/ptq5mlDNu/RAcmEqqww7f+973uLXl9lvRTasqDu1/0tspEZUDGA5gg2HzfZq+60uJhjeGcy2LJRkJSx6EDHjohgYVXsp0TtIZBgwYgMLCQq5zN+PHj+fWVkqnQUQfENH2BH+25BOJKB/AWwAeYIzpfb/nAfQFMAxADYCnzc63I5YUT58+fWwdnwoeyt0KPsheyUlFoiS/VDeyXbt2cR+y8Ow5pVw9YYxdZbaPiI4RUSljrEbTZk3YByeiTMQcxmuMsbcNbR8zHPMCgH/YMd4qPXr0wIEDB7i1F5Q7cNeuXVMOi4KO7JUcJ6SqztW9e3dUV7vWDGsDz0lxt8OTdwHcrj2+HcCy+AMoZu1fAOxkjP0xbp9x9u96ANtd2pOQSZMmedFsSmQryIfdYQSVVHNNusPQe1FOJqRzcnJaH99zzz22z0+GW6exEMBkItoLYLL2HETUnYj0lZDLAdwK4MoES6t/IKJtRFQBYBKAX7q0x1f4TUFeESz0XpTVcH0jxpqhvG9eroK7GGNfAWhXZooxVo2YohoYY58isdAzGGO3urm+HaZPn46lS5eKupwigPTr1w/79u2TbQZXOnTowH0VLDJrar1795ZtgkICVsbynTvHFu3C5jAAb0IOIuM0OnTooCIkLeIkrby8vJy/IRywUh7ASfffLsYVPJExIDzjM3Qi4zQAb97AMBK/3GdlibmystIja9yRqCpaPCKKEh04cAAFBQVIT08Xqs3ixY0yUk7Di5DyKGBlidmPNVn9Rn19vesoXTvcf//9nrQbKachenjideUsL7FbR9StNmwUEF1NTK8Ry5tIOQ0AmDVrVuqDOCEiTsKrIsGi64iGnbS0tNawcDeJgP3797d0XH5+vmeFgiLnNKx4+/T0dAGW8EH9uIOBMQrUqaPPzc3F3r17LR17++23pz7IIZFzGkDq5B2e405R+p+K4GA3GW748OEA7OVQeRmNHEmnMXnyZGHXEqX/qbCO7PB+u3z55ZcAgO3brWVZ8CrrZ0YknUZYiV//v/TSSyVZ4m+ShfeHoWDw6NGjPW0/sk7jwQcfFHYtEXMk2dnZ7db/rd6ZRNgXlALMQSwYbGTqVO4VM9sRWadhp1J5fOy+3cI+XqzNxy8fuwkYEhE7wKMAs91Vh/iel7EOSxh6FIkYNWqU59eIrNMAkvc2ysrKWh/H1z/grRhmJX4kvtjM8OHDA1HDMzMzEz179uTSlt24A92R6oFnxkpYsnoUXpZQnDZtmpBYpEg7jWS9jSNHjgizw0rhW2OxmeLiYmzevFlIzoRbGhsbuSm9OS014KfAs1QFeOxirJshag4r0k4DAB5+OFEBdX+HRfOWUPCa8ePHSy1S7IciwVax+9nqdTPGjx8vrOfpue6JdlylVmxnCxFttHu+l5h9ofx0d4onvnThhAkTJFlijXXr1nG/w9pBZIKYW5zIYwBikzFF6J7oTNK0TYwzNXbO94zZs2fLuCw3Pv7449BO7ClSM2DAAKHFrj3XPfH4fC506NAh8FIHbib29ElWsyGZ3yt+R50bb7xR6PVE6Z4wAKuJaBMRzXBwvuc89NBDsi4tHX2S1WxIFsSK336F99zOpZdeKryXKUr35HLG2AgAVyMmy3iFXUOdiiVZJSMjQ7jHtopXmawK8fCe25HxnU3pNBhjVzHGLk3wtwzAMV2GIJnuiVZoGIyx4wDeAaDHuVo6XzvXsViSVfwadq0yWf3DxRdfnPKYZHEpRUVFlq/Vo0ePpPtFlnkwIkL3JI+IOuqPAUzBd/omKc8XjcjwcgU/3HT77Syv796923SfXnYhWR0VO2rwR48eNd1XXFwsvKiPjgjdk24APiWirQA+B/AeY2xlsvNlooYCwcRNt9+OqvyAAQNM9+mBVjwU+IxBW4m4++67XV/DKWTnDfMLo0aNYhs3bkx9oENaWlrw2GOPeda+kczMTHTt2tVXKfRDhgxBRUUFJk6ciDVr1sg2x1NycnKQlpbmyyFgp06dcObMmXbbr7vuOgwdOtTTaxPRprjwiFYiHxGaiLS0NPzsZz8Tcq3GxkZUVVX5SpeloqICQEyIOOycO3fOE4eRaj7CCrrDMOaTlJeXe+4wUqGchgl9+/YVer2DBw8KvZ4VeGSmRpVk8xF2MY4GbrnlFm7tOkU5jSTMnTvX0Xl67Ygg1RpViMVJnsgDDzzgiyBE5TSSQESOhin6Hbq5uRndunUD4K4CddgZPHiw1IQ2GaSa6IynZ8+evklUjNYn5QC3wxRd4UtGGntpaanwa9qhpKQEeXl52LZtm9SENrfEF76xIjNgNzHtjjvusHW8lyinYYF58+bJNsERNTU1bZ7LTGrTYyEGDRrUuu3EiRO+XLWwS/xKnlWZAav86le/4tqeW5TTsMhdd92VcPtFF10k2BLnyKx/qee17NixQ5oNXpGXl4devXp50nbfvn19V9tFOQ2LdO/eHffee2+bbfn5+dizZ0+7Y/V5DK9wkgbtpy+e34dNgL05qLNnz+LQoUPcbejTp4+wpX87KKdhg5KSEtxwww2tz+vq6hIel0qp3G3ujJOIQ78UFSosLGw3bPIjbueg3Nbq7NevH2699VZXbXiFcho2GTx4sOs29CxdrxLvkiFaBDueU6dOSb2+GTyCsYy4jbSePn06J0v4o5yGA3hNjHqR4h9PfKRpENMGRMAzGMstc+bM8UU8hhnKaThExIqKk+Cw+PV/0ZGmaWlpkYu54MmMGTN8//752zqf43W1LyciRnp1alm0tLSYxlzk5ORIS+cOArfeemsgJomV03BBTk4O5syZI9uMwHDu3Dnfzml4Tao6q9OmTbOlCi8T5TRckpaWpgr3KJKSk5OTtM7qT3/6U4wYMUKgRe5QToMD+fn5uPPOO2WbofApyYaME1mqdm0AAAlySURBVCZMQL9+/QRa4x7PxZKI6GJNJEn/O0NED2j7fktEVYZ93ktee0SPHj1CVdHcz7P3YWHgwIGYOHGibDNs47lYEmNstyaSNAzASAD1iBUX1vlvfT9jbHn8+UEiJycHjzzyiGwzXJObm9tG5DrRbL6VpCzFd8THx8ycOdPXsRjJEC2W9H0A+xlj/GNufUJGRgZ+/etfyzajHeXl5ZaPjY8eTbQawjspK+wY42NuvvlmdOnSRaI17hAllqRzE4A34rbdR0QVRPRSMi1Xr3VPeJKVleV6VSXZ8MAsqtMsv6RTp06orKx0ZY8bZEeh+ompU6cGKskxEaLEkkBEWQB+DOB/DZufB9AXwDAANQCeNjtfhO4JT9LS0kwDwMrKylKebxwexGMW1WnmaBIVp/WC+FBsfZkxSFGoBQUFntVrnT17Ni677DJP2hZJygILjLGrzPYR0TEiKmWM1aQSO0JMXW0zY6w1m8v4mIheAPAPa2YHh3nz5uGf//wnPvnkk9ZtXgVgJSvskpmZ2eqIiMiTH3J8KHYqOce8vDycPXsWhYWFXOM33Ly+06dPO1ZuT8bcuXND0+PyXCzJwM2IG5ro6moa1+M7EaVQMWnSJEybNq31+cmTJ4XbYOy5+OXOrxfgOXXqFLKzs7m165fXBwAjR47Eb37zm9A4DECMWBKIKFfb/3bc+X8gom1EVAFgEoBfurTHt4wYMQLXXmtrROeIoK3565w/f162CdwpLy/Hj370o9AtXyuxJMEwxrBgwQLZZig85sYbb/StNrAVlFiSjyAizJs3L9BLborkzJ07N9AOIxXyKs1GnJkzZ+L8+fNYuFC6fG2o0CdXZXDNNde0q0weRlRPQyLZ2dmYM2eOcDU3nbKyMqkVyr1AlsOYPXt2JBwGoHoa0tF1Y6urq/HCCy8IvfaRI0cAAKNHj8bnn38u9NoyyMjI4F6Rffr06Rg4cCDXNv2Omgj1EY2Njfj000/bxHQo/ElaWhruueee0M5NJZsIVT0NH5GZmYlJkyZh+PDhWLRokWxzFCZMnjwZ48ePl22GNFRPw8d88MEH+Oyzz2SbodAYOXIkrrnmmlAFapmhehoB5aqrrsKECROwfPlybNmyRbY5kaVr166YNm0ad5mDoKJ6GgGhsbERa9euxdq1a6XakZ+fbyoSJYry8nIhWbtEhNtvv90zyUU/k6ynoZxGwDh37hxWrFiBbdu2yTYl1Nx///0oKiqSbYY0lNMIIS0tLXj++eelJL+5wYtlT5489NBD7bRjooia0wghaWlpmDlzJgBgxYoVgYmzSOQwysrKWmNGZPHoo486EqeKIqqnESJWrlyJTZs2+fpO7jdEKOUFETU8iRgtLS3YsGEDVq9eLdsU31FYWIhbbrlFivh2kFDDk4iRlpaGcePGYdy4caivr8eTTz4p2yTp3HDDDRg8eLBsM0KBchohJzc3t7ULfvbsWbz++uuorq6WbJUY7r77bhQWFqJDhw6yTQkVrpwGEf0XgN8CGAhgNGMs4ZiBiH4IYBGAdAAvMsb0Cl9FAP4HQDmASgDTGWNfu7FJYU5eXh7uuusuADGZgvXr13sW91FaWoqamhpP2jZjypQp6NGjh6XCzQrnuJrTIKKBAFoA/F8A/yeR0yCidAB7ECv3dxTAFwBuZoz9m4j+AKCWMbaQiB4G0JkxNjvVddWcBn/q6uqwY8cOrFy5UrYplrn++uvRuXNn5SQ8wLM5DcbYTu0CyQ4bDWAfY+yAduybiIks/Vv7P1E77hUAawCkdBoK/uTn52PMmDEYM2YMAKC5uRktLS1YtmwZzpw5I3VJ9MILL0RVVRUefPBB5ObmJlR8U4hDxJzGhQCM37ijAMZoj9uILRGRqdgSEc0AMAMAevbs6ZGpCp309HSkp6fjJz/5SZvtzc3NSEtLw8mTJ1FXV4dXX30Vffv2xf79+9vIJNhh4MCB2LlzJ6ZPn47i4mKUlJSgpaVFxU34lJROg4g+AHBBgl2PMMaSSRa0NpFgm+0xEWNsCYAlQGx4Yvd8BR/0H3JJSQlKSko8i3NQDsO/uBJLsshRAMZBZw8A+vS9HbElhULhA0QMDr8A0J+IemvSjDchJrIE2BNbUigUPsCV0yCi64noKIBxAN4jolXa9laxJMZYE4D7AKwCsBPAUsbYDq2JhGJLCoXCv6gwcoVC0Q4llqRQKLihnIZCobCFchoKhcIWymkoFApbBHIilIhOADhk4dAuAIJVD689YXgNQDheR5ReQy/GWMKiI4F0GlYhoo1mM8BBIQyvAQjH61CvIYYanigUClsop6FQKGwRdqexRLYBHAjDawDC8TrUa0DI5zQUCgV/wt7TUCgUnFFOQ6FQ2CJUToOI/ouIdhBRCxGZLisR0Q+JaDcR7dNqk/oGIioioveJaK/2v7PJcZVEtI2IthCRL7L3Ur2vFONP2v4KIhohw85UWHgdE4notPbebyGiuTLsTAYRvUREx4lou8l+558FYyw0f4hVRb8YsVqjo0yOSQewH0AfAFkAtgK4RLbtBvv+AOBh7fHDAH5vclwlgC6y7bXzvgKYCmAFYtXcxgLYINtuh69jIoB/yLY1xeu4AsAIANtN9jv+LELV02CM7WSM7U5xWGuhY8ZYAwC90LFfuBaxIsvQ/l8n0RY7WHlfrwXwKouxHkChVrHNT/j9+2EJxtgnAGqTHOL4swiV07BIokLHF0qyJRFtii0DMCu2zACsJqJNWtFl2Vh5X/3+3gPWbRxHRFuJaAURDRJjGlccfxaBU1jzS6FjNyR7DTaauZwxVq1VcH+fiHZpdxdZWHlfpb/3FrBi42bEcjPqiGgqgL8D6O+5ZXxx/FkEzmkwbwsdCyHZayAiS8WWGWPV2v/jRPQOYt1qmU7Dyvsq/b23QEobGWNnDI+XE9FzRNSFMRakZDbHn0UUhyfJCh37gZTFlokoj4g66o8BTAGQcJZcIFbe13cB3KbN3I8FcFofivmIlK+DiC4gTSGMiEYj9jv6Sril7nD+Wcie5eU8Y3w9Yh70PIBjAFZp27sDWB43c7wHsVnyR2TbHfcaigF8CGCv9r8o/jUgNrO/Vfvb4ZfXkOh9BfALAL/QHhOAZ7X922CywiX7z8LruE9737cCWA9gvGybE7yGNwDUAGjUfhN38vosVBi5QqGwRRSHJwqFwgXKaSgUClsop6FQKGyhnIZCobCFchoKhcIWymkoFApbKKehUChs8f8B35NzOgalV0gAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from __future__ import print_function\n",
    "from dolfin import *\n",
    "import numpy as np\n",
    "import mshr\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "\n",
    "T = 20.0            # final time\n",
    "num_steps = 1000     # number of time steps\n",
    "dt = T / num_steps # time step size\n",
    "c = 0.2 # wave speed in the material\n",
    "PI = np.pi # defining PI as a parameter\n",
    "epsilon = 5 # amplitude of the forcing\n",
    "omega = 0.1 # frequency of oscillation\n",
    "sig = 0.05 # width of Gaussian of the forcing\n",
    "#eps2 = 0.1\n",
    "\n",
    "Nmesh = 50\n",
    "domain = mshr.Circle(Point(0,0), 1)\n",
    "mesh = mshr.generate_mesh(domain, Nmesh)\n",
    "plot(mesh)\n",
    "\n",
    "# Define functionspace\n",
    "V = FunctionSpace(mesh, 'P', 1)\n",
    "\n",
    "# Define boundary condition\n",
    "w_D = Constant(0.)\n",
    "\n",
    "# Define initial value\n",
    "#w_exp = Expression(\"cos(PI*sqrt(x[0]*x[0]+x[1]*x[1]))\", eps2 = eps2, PI = PI, degree = 2)\n",
    "#w_n = interpolate(w_exp, V)\n",
    "w_n = interpolate(Constant(0.0), V)\n",
    "wdot_n = interpolate(Constant(0.0), V) # initial velocity is given to be zero\n",
    "\n",
    "# Define boundary conditions\n",
    "def boundary(x, on_boundary):\n",
    "    return on_boundary\n",
    "\n",
    "bc = DirichletBC(V, w_D, boundary)\n",
    "\n",
    "# Define variational problem\n",
    "w = TrialFunction(V)\n",
    "wdot = TrialFunction(V)\n",
    "v = TestFunction(V)\n",
    "\n",
    "# External\n",
    "f_ext = Expression('epsilon*exp(-(x[0]*x[0] + x[1]*x[1])/(2*sig*sig))*cos(omega*t)', \\\n",
    "                   degree = 2, epsilon = epsilon, omega = omega, sig = sig, t = dt/2.)\n",
    "\n",
    "# using Crank-Nicholson scheme\n",
    "F = w*v*dx - w_n*v*dx - dt*wdot_n*v*dx - \\\n",
    "    c*c*dt*dt/2*f_ext*v*dx + \\\n",
    "    c*c*dt*dt/4.*(inner(grad(w),grad(v)) + inner(grad(w_n), grad(v)))*dx\n",
    "a, L = lhs(F), rhs(F)\n",
    "\n",
    "# start the update\n",
    "# Time-stepping\n",
    "file = File(\"./Membrane2/Displacement.pvd\", \"compressed\")\n",
    "filev = File(\"./Membrane2/Velocity.pvd\", \"compressed\")\n",
    "\n",
    "w = Function(V)\n",
    "wdot = Function(V)\n",
    "t = 0\n",
    "for n in range(num_steps):\n",
    "\n",
    "    # Update current time\n",
    "    t += dt\n",
    "    f_ext.t = t + dt/2\n",
    "\n",
    "    # Compute solution\n",
    "    solve(a == L, w, bc)\n",
    "    \n",
    "    # update velocity\n",
    "    wdot.vector()[:] = 2*(w.vector() - w_n.vector())/dt - wdot_n.vector()\n",
    "    \n",
    "\n",
    "    E_kin = 0.5*assemble(wdot*wdot*dx)\n",
    "    E_elas = 0.5*c*c*assemble(inner(grad(w),grad(w))*dx)\n",
    "    E_tot = E_elas + E_kin\n",
    "    # Plot solution\n",
    "    #plot(w)\n",
    "    #plt.pause(0.5)\n",
    "\n",
    "    # Update previous solution\n",
    "    w_n.assign(w)\n",
    "    wdot_n.assign(wdot)\n",
    "    if (np.mod(n,10) == 0):\n",
    "        file << (w, t)\n",
    "        filev << (wdot, t)\n",
    "        #print(\"Elastic energy is \", E_elas)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x118bc1198>"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVUAAAD8CAYAAADHaDe8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO2de7Bd113fP19dvSwlwXbkh/wIdkEwKMS46cUEQktC4tT2lCgBwtgwxIVQ1S2eoQwvpR4ChaEjkkKAwSQIGioewXgoxioWfqmFhNCA5SR2rNiuhWNiRcKy7BDjhx5X99c/zt7X++679t5r7b32ea7PzJl7Hmvvvc69537Ob/3WS2ZGIpFIJOKwatQVSCQSiWkiSTWRSCQikqSaSCQSEUlSTSQSiYgkqSYSiUREklQTiUQiIkmqiURiIpB0laRHJR2UtMPxuiT9Wvb6g5Je33SspLMl3SPpseznWV3rmaSaSCTGHklzwM3A1cBW4DpJW0vFrga2ZLftwIc8jt0B7DOzLcC+7HEnklQTicQkcAVw0MweN7OTwC3AtlKZbcDv2oBPAmdK2txw7DZgd3Z/N/COrhVd3fUEo2DTpk12ySWXjLoaM83iqc/2ev5Va17X6/kT9dx///3HzOycLue48s1n2DPPLnqV/fSDJw8AxwtP7TKzXYXHFwJPFh4fAr6pdBpXmQsbjj3PzI4AmNkRSed6VbiGiZTqJZdcwv79+0ddjannpSOX1rx6Qc9Xf6bylTM2f77nayck/X3Xczzz7CKfuNPvc7LhgieOm9l8XZUcz5Xn2FeV8Tk2GhMp1UR86gU6XrjqmkQ79RwCLi48vgg47Flmbc2xT0nanEWpm4GjXSuapDqDTJJAfSm/pyTZqeM+YIukS4EvAtcC31sqswe4UdItDJr3X85k+XTNsXuA64Gd2c/bu1Y0SXVG6FOkxxaPNxcCNq1a31sdyhTfbxLs5GNmC5JuBO4C5oCPmNkBSTdkr38Y2AtcAxwEXgR+oO7Y7NQ7gVslvQf4AvCurnVNUp1iYonUV5qxzxNLwimKnQ7MbC8DcRaf+3DhvgE/7Hts9vwzwFti1jNJdcroKtJYAo1BXV26CDf/HSW5JvogSXUK6CLS2BL9h9Prlj0+f+5E1PPnuOodKtqUIkj0QZLqBNNWpm1FWhZmX8dAOxkX31dbwSa5xmOBxbFq+QyLJNUJpI1MQz/cbWUYi7rr+wi3/H59JZui10RXklQnhL5FOmqJhtAmxZD/LkIi2BS9JtqQpDrm9CXTGBJ99vSGzucocvbci62OC5FsmxRBkmsihCTVMSVUpn2KNLY821wnRLjF9xlTsEmuCR+SVMcQX6H2IdJhCTQUV718RBsq2CTXRFeSVMeIUch0XCXqQ7nuTZL1EWySazxO2aqJytXHIkl1DIglU58PcEyJHl14VbRzlTl39XPBxxTfm69gY8k1iTWRk6Q6QmLItG+R9inONtf1la2vYGPJNUWtiZwk1RHhI9Rhy3RUAg2hXEcfyfoINqZck1hnmyTVIdM1Om2SaYhIY0n0WKTzbGrR5A+VbJNgfeSaotZEHUmqQ6RLdBpLpm1FGkucba/hK9zi+/MVbKhcU9SaqCNJdUg0CbWNTPsS6TAEGoqrTk2i9RVsn3JNUevskaTaM31Epz4yDRFpF4keO/WK1se62LTmef9rl+pdJ1kfwfYt11kT62lbNdFD9tqSpNojsaPTWDJtI9HY8mxznSbhFt+Xj2D7kGsSayJJtSeGLdSYMh2WQEMp16tOsj6C7UOuTVFrSgdMP1GkKukq4FcZ7P/y22a2s/T6TwDfV7jm1wHnmNmzkp4A/gk4DSw0bFM79rRt7vclUx+RdpHoM6c2tj7WxavXvOBd1ley+e+gSa7gFmyTXFPUOloknQ38EXAJ8ATwPWb2JUc5p6fqjpd0GfCbwKuAReAbzax2Fk5nqUqaA24GrmSwRex9kvaY2efyMmb2AeADWfnvAH7UzJ4tnObNZnasa11GTczotG+Zhoo0tjxDr+Mj2+J7cgnWN3oNjVzrotYk1qGwA9hnZjsl7cge/1SxQIOnnMdLWg38PvD9ZvaApFcDp5oqEyNSvQI4aGaPZ5W/BdgGfK6i/HXAH0a47ljRRqhtotM6mcYU6bAk6ku5Pk2S9RWsS64+aQHfqNUnHZDE2pltwJuy+7uBv6AkVeo9VXX824AHzewBWNoksJEYUr0QeLLw+BCDPbdXIGkDcBVwY+FpA+6WZMBvmtmuCHUaKrGE2odMfUUaW6JfOll9vrPW+jfxqwiRbP47iCnXmFFrEmtnzjOzIwBmdkTSuY4ydZ6qOv5rAJN0F3AOcIuZvb+pMjGkKsdzVlH2O4BPlJr+bzSzw9kbuUfSI2b2sRUXkbYD2wFe85rXdK1zNPoWap8ybSvSOmHGOj5UvMX3UiXYvuQaErXOUgfWKeZChvZtkrS/8HhXMcCSdC9wvuO4mzzPH+KpnNXAtwLfCLwI7JN0v5ntazqoK4eAiwuPLwIOV5S9llLT38wOZz+PSrqNQZi+QqrZL3gXwPz8fNMvYyjUCXVU0WmTTENF2lWgbXFd11e0TYLtItfQqDXlWb05VtdJbWZvrXpN0lOSNmdR5mbgqKNYnaeqjj8E/GXe3yNpL/B6oFaqq+pe9OQ+YIukSyWtZSDOPeVCkr4C+Dbg9sJzGyW9Mr/PIIfxUIQ69U6fQj268KpKoR5beJVTqMdOvaJSqM+c2rh0a+JLJzcuu40TbepW977rfmdVv+e6v43r7+n6ux9bPF67WE6XLcdnlD3A9dn96yk4pkCdp6qOvwu4TNKGrNPq26juK1qic6RqZguSbswqMAd8xMwOSLohe/3DWdF3AnebWTF0OA+4TVJel4+a2Z1d69Q3fQk1dmTqG5XGkOc/ngzbErrMmWvDtzIu17suks1/F20i19CoNeVZh85O4FZJ7wG+ALwLQNIFDIZOXVPlqbrjzexLkn6ZgZAN2GtmdzRVRmZj0ZIOYn5+3vbv399csAeGLdS+ZNpGpF3F2YY2soXmVEFdx1bVeNeqYVhVowRcuVbXmNa6dMCoxJrlDjuNGb/0da+wn/uTr/cq++6v+ZvO1xsX0oyqALoKNbQzKlSoTTINEekoBOqiXA9fyebvtUqubSLXWFFrilinmyRVT0KEOuzoNJZMu4j0uRPdJfyqdc3CDJVs8b27BNskV5dYYWXUWjVCIIl19khS9WBchVon0z5EGkOcoedvEm2x/r6CDZFrjKh1VsW6YKvGchnJvklSbWDUQg2NTn1k6ivSviXapg51kvUVbFu5hkStfYg1MRkkqdYwLKHGiE5jyLSNRF84sTb4mDo2rjtZ+3qxjl0F2yTXLlFrH2KdxGh1FklSrWCUQo0dndbJ1FekseUZcp0q0fpGsfn7D5VraNSaxJqAOIP/p45xE2rVAPamQfD/eHJ9pVCfO7G+UagvnFi7dBslxXrU1aXpPdX9PqD6y8n1u3f9nVwTBlyTBXwnCVSRJgeMNylSLTGOQnVRJYAuTfyu8jx+Yk2r49ava1xNbRnFerqi2KYUQV3kGhK1dkkH+ESs09ZxNSskqXpSN60wp61QY+VO2zbzQ2TaVpxtz9kkXF/BtpFrSK41iTWRk6RaoCpKbTuwf1hCbSNTX5H2IdEQytevk2z+ntrINTRqHRexjjMLNje2W/P0ScqpZoyzUKtyp21ypj450uMn1izdxg2futXlX6t+L3W/S9fv3pXnrsqzlvHJsa44T0VLKeVXx48kVcI+mG2E6urACBGqi9gyHWeRVuFT3zZydeHbiRVLrFWrW7lIYh0vUvO/hqaOKV+hrjhv6R8vpLlfJ1MXPlFpF04dj/cRWrN+odVxxfdQlR6oSg08d2L9ipRAVa7VNx3gmwooU04FpKFWk8nMR6q+zf5hCrWuuV+mTWTaNio9dXz1iltMXOcPvYZPaqBMjKi1TcTqShH5RKyJ8WampdpWqC5iCtVFlVDL+MjUlz4FGkKbOtS917qUQJlxFGtKA4w3qflfos3QKZ99eNoK1VemUN3UDxXpOFOsn0+6IH/vrrSAKyUQmg6IkQqo2xY7x3dEwDilARZs1djtyjsMxvs/qEd8v9Wbmv1tevljC7WrTNuIdPHEXPAxPqxad9q7bIhgm+RaFiusHH7lGnrVl1irNhRcdo4JHWo17cxk8z9mHnXFOSZIqCFN6sUTc8tufVG+ju+1fN9LVVqgSzrAlQP3SQWUaTMiwEVKA4yWmZSqiz46pvoUalVe0Cdv6iugULHp+FzQzZeQeoTItYzrdxqSZy3TJNYYQ61SfnX8iCJVSVdJelTSQUk7HK+/SdKXJX0mu73P99jYxPqwjVqoLmLI1EdgXSXZ5Ty+gvV5r75Rq2t0gGuyQJvOK59FnNuKdVaQdLakeyQ9lv08q6Kc0zWS3iXpgKRFSfOF56+UdL+kz2Y/v92nPp2lKmkOuBm4GtgKXCdpq6Pox83s8uz2c4HHRiFWs7+pYyqWUF3/zG2i0ybBNImqqzxD8b1eTLmWaZsOiCFWn45PH2YoWt0B7DOzLcC+7PEyGlzzEPCdwMdKhx0DvsPMXsdg6+rf86lMjEj1CuCgmT1uZieBW4BtQzi2F2LkUYt0EWqZKqHW4SNTF20kOndCQTdfmuriE702ydX1xdQ2HeAj1jJNYk3Rai3bgN3Z/d3AOxxlKl1jZg+b2aPlA8zs02Z2OHt4AFgvqTGxHUOqFwJPFh4fyp4r882SHpD055JeG3gskrZL2i9p/9NPPx1cyZC5/XWENvvLDEuodRLxkWkVMSTZdK46fAVbhY9cy8QSa5E2C420Eesoo9XTtmqpE6/pBmzK/7+z2/aAS51nZkcAsp/nOsp4u6aC7wI+bWYr9xgvEWNIleu/wEqPPwV8pZk9L+ka4E+BLZ7HDp402wXsApifn3eWCaVrsz80j9pGqG3yp03RqYsmkQ6T8vVOr3P/ufM62/qVw7AWT8zVDs86dXx15TCs4yfWrBh65Rp21bQpYXm4VdNQqzbjV30Yp7GrNRwzs/mqFyXdC5zveOkmz/N7u8Zx7dcCvwi8zad8jEj1EHBx4fFFwOFiATN7zsyez+7vBdZI2uRzbAzafFv7NPvr8Bn0HFuobaLTqqivSwQam6ZItuo9dIla20SsVcOtinTNr85qGsDM3mpmX++43Q48JWkzQPbzqOMUrVwj6SLgNuDdZvZ3PnWNIdX7gC2SLpW0FrgW2FOq2PmSlN2/IrvuMz7H9kXoVNSQKDVk+mlO1w6p2DIdZ/qQq4thiTWUNl/4U95ptYdBRxLZz9sdZYJdI+lM4A7gvWb2Cd/KdJaqmS0ANwJ3AQ8Dt5rZAUk3SLohK/bdwEOSHgB+DbjWBjiP7VqnIj4fptjN/jIhq01Bu/ypi1HJdO4lLbtVPdf5Oi3kWkWdWF0dWEXaiLVI19EAsxKt1rATuFLSY8CV2WMkXSBpL1R7Kiv3TkmHgG8G7pB0V3beG4GvBn66MBzUla9dhsyipCeHyvz8vO3fv9+rrEuqfeZSu+ZRYwi1SqYuYkalsWR5+ox2n0lX7tWVb4Xq6bB1013LedbyMoLlHGt5Smt5Omt594DyVNbyUoHl/GpxGmt5icCq6au+uVVJ99flOH04++vOsSs/8l1eZW/9lt/sfL1xYapnVA07So3RMVVmXIRajjRdt1i0Pbcrcg2NWkNGBoRGrLHTAEVStDo+TLVUXYR82EIH+TcR2jHVl1B9m/p9CLMtoXItE0usZZoWAQ8Ra+xOKxfDzK2eXly1NPOs6TZNTK1UY0SpZUIG+Zf/eWIJtarHuiwIV4QWKtNxxFf0vlFrqFibJls0tTwgLL/ahRStjoaplaqLug9Z383+OkKE6sIl1DK+Tf1xlakLX7mWiS3WNh1XRerSAKHRamL0TKVU20SpIcRs9vct1GmITptoil6rotYiscVaJmTGVZtFV3LGfZbVLDCVUnURM0ot0rXZ3wUfofowDJmuOtlcJgZ172WYYg39O4d0WqVodbyZGakW6StK7TrAPyRKjSHUWNHpqpPNt7pysWmKWot0FWsdXUYDpGh1cpk6qfqMS62jS5RaJmSA/zCF2lWmsaXYl2SHIdbQNECZtp1WKVodX2Z2j6qckA9jSJTaVx41hlDbECK7OY/vsNM1rePitRbDHLWyLi/JOZlg7oSWTRbQ8bllEwWqFmRxLcTiWoAlp7zwimufqyLlRVeKlBdcqaO8UaCLvhdaOW2roqa7JoWpilS7dlDFjFKLtPlgxRZq2+i0KXqcO77y5oPvcTEi2NgRaxNdJwUUiRWtpuFVw2OqpOoipIOq9jwBQ6i6Nvv7EGoodSILEejqE8tvdfgKtg0xxdpHGqBIyBCrOrr0HSTaM/VSLRLyIQv58NZRF6X6biE9TKFWiatOeGVx1knUt1zd9drKddhiLTKu0WrqsIrP1Eg1tIMqZN+pvqJUF029zH0LdcU5PEQaAx/Blmkj177FWqSvaDUx3kyNVGMyqii1qdnfl1BdcmqSaSghudcqwdbJNaguLcXqQ1/RapGQpQFTCmD4zIxUQzqoigwzSi1T11ESU6grji+Jq01O1EegIYKtq1/V+6itq6dYi4wyWg2ZxZdSAKNlKoZUdW3690GMKLVIMWrqS6hVkWnltSJ2KBfP5RputfoELKxbWb5YdtXJsCFYVcOtivgMtSoPsyoPsSrub9U0xKq8r1UVIcOrRsWiqdOXyqQyM5GqL8WmVZdxqUW6Rqm+zVBfoVY194u06UCKRUjqoVwuxgSCvtMAbQnpsCqSUgDDZSak2rbpXyTWuFSfKLVNsz9EqCuOdQjVeQ1PkZ5e//LN9bwvVdcry7VLOmAUaYCQVazadlilFMDoiCJVSVdJelTSQUk7HK9/n6QHs9tfS/qGwmtPSPpstv+L3x4pBUbR9O87Si3i0+xvK9RyRBja+16UZ5NIi8+HyLVpKFddHWOKtc3EgLbRatsOq1lF0tmS7pH0WPbzrIpyTk9JepekA5IWJc0Xnl8jaXfmp4clvdenPp2lKmkOuBm4GtgKXCdpa6nY54FvM7PLgJ8HdpVef7OZXT7qPWqqmv7jEqV2wae5X6aumR8ScVbRJGRXfcr4pAN88PliakoDdIlWfUkpACc7gH1mtgXYlz1eRoOnHgK+E/hY6bB3AevM7HXAvwD+vaRLmioTI1K9AjhoZo+b2UngFmBbsYCZ/bWZfSl7+EkGe24PhRhN/zomIUptK9QyvvLLWVj38s2Xpmv0KVbn9QLTAGV8o9W2KYC6aHWGFlnZBuzO7u8G3uEoU+kpM3vYzB51HGPARkmrgTOAk0Bj72AMqV4IPFl4fCh7ror3AH9eeGzA3ZLul7S96iBJ2yXtl7T/6aefblXRPpr+XegapbaZftpFqHUUBeoSqev1JtlWydUVQTflWX3EOk7RaszPmYth5FUXF7W0xXfTDdiU/39nt0oXODjPzI4AZD9d20iHegrgj4EXgCPAF4D/ZmbPNlUmxpAq1yfROU5F0psZSPVbC0+/0cwOZ/tp3yPpETMrh+GY2S6ytMH8/LxB92X+lh0XIWqNMSU1x3cIVRNFmTQJNVSmIRGo7zlckj+93l23ueMr61cceuV6vQnXMKvyilZFqlazKlK3ilWRphWsco6desWK7axzji68asVW1jk+K1eNkGN16T9J9wLnO166yfP83p4qcAVwGrgAOAv4uKR7zezxuoNiSPUQcHHh8UXA4XIhSZcBvw1cbWbP5M+b2eHs51FJtzF4Iyuk2oa2eSTffGqXpv8wotS+hOoTYdZRN4ogP3e5fvk5y8eGiDV0HGsV5bGrZVzLA3albknAWcDM3lr1mqSnJG02syOSNgNHHcW8PFXie4E7zewUcFTSJ4B5oFaqMZr/9wFbJF0qaS1wLbCnWEDSa4A/Ab7fzP5f4fmNkl6Z3wfexiBpHJ2Quf5VtG2SjSJKrWvuNgm1qtld1WQP6WzyLd90rbr619E2DdA1t1okdgog5VXZA1yf3b8euN1RptFTDr4AfLsGbATeADzSVJnOUjWzBeBG4C7gYeBWMzsg6QZJN2TF3ge8GviN0tCp84C/kvQA8LfAHWZ2Z9c6jYK6RajLDDuXWh42VfVaHb6Ca0OdYH1TDHWReB8TFUInBMSeDNB2FMCyc0zPeNWdwJWSHgOuzB4j6QJJe6HaU1m5d0o6BHwzcIeku7Lz3gy8gkGgdx/wO2b2YFNlokxTNbO9wN7Scx8u3P8h4Iccxz0OfEP5+bbEyKfG3HcdRh+lhgq1LLYqmfaFq5nvSgm48qzlVECXNEBobrVM2xRAXV511lMAVWTpxLc4nj8MXFN4vMJT2fO3Abc5nn+ewbCqICZ2RlXTN2rfs6ia9nKvos0GckVi7Xw6jkJtuk65PlUjA3yIveHgOKcAikz5eNWxYGKlGkLfQ6lCmv5lqub4u6LULp1TTc3hJqHGauqH4LpmqFi7pAGacqvjnAIoMqq8qpk4dXy1122amAmpjoo+FtVoQ9OSfW3ENUy61q/q/fexPXaRtrJo2wpKjAczL9U+86llunwjdx1C1fR8FbGE2mZzwLp6NInV9/03iTUk3dLXFOMivlOmu3RWJbox81Ktoo98apmQpn9bukSpbYXqs0h1G7n65HyriLXtC3RLARQZVV61iQkdATA2TI1Ufb9xY8z3L9IlnxqTNlFqTKG2jUJDj6urV1/RakzGISWUOqv6ZSKlunjqs7Wvx/7QtBn0H/Ofp22Pf0h01nbKadXygb638rl8KMpzGNFqlxEXfXfC9J2ySoQzkVKNxTiuR9mm6e81S6ghSvV9reqcrgVNfPasirkjK7SPVkOpSgHEyqvGXrh6RmZWjQXTNZbBge+HadidVL7bpXSZPeVLm2a/a9jSsudeqh4kv3CGVuwxlQ/Ud+095aI4+H9hXf3EgOLA//JeV6Ombv+qiWdxOJ1348ZMR6pVTMOe63XjM7v26Ofny6PMYvS5+iVbEmp+v3grl1l2bOCY0mEO9Yo16aItbVJQsfsPEn5MfaQam6pm2aTuGhkapZYFWJRp8WfO6uOLpTOs/B7PI9fT69tHkuVo1Yc+VrDywXcpwMRkkiLVDlQNgxlVD2+fO5wWz98k1NXHF5duc4Vb8bW8fDlqdZ0/Zt2L5+98zohD30ZJGqsal6mQ6qiGU8Wm/E/aZqsUqJdGVc956BqoLqG6RLr6xYWl27LnHXJddl5PsbYdCRCbus6qYU7DHMfO11kjNf8b6Htbi76JHb0WhV0W6lI0+uLy1ZlWvXCSxY1rS8+//NFbWL9q6VwLZzjm27dYwb9qt4BppG4ngMTwmYpIdZLpMiNnWLii1LJQc3Jxrnrh5NLN9TgvN1eIWPPzuTquxoFhd1bFHlZVJE0A6I+ZjVSH2Uwa91V4fKPAsuSKQs2b+sDLIj2+MjexuH7tssh1YcPqpQgXVi1FquXotC5a9Y1KXedoEwXHpjisaqpY1EQEDbEZ7//2KaWPsXtV8+v7uE65h78s1GUyfeGll+9vPGOFaFcDCxte/hiW0wChowF8RwGM23jVnKkbqzqDpOY/aapfHVVN8XKzH0pCfeGl5UItPVeUa96JtXL4VTdCI9BhrgEwCtKsquEQRaqSrpL0qKSDknY4XpekX8tef1DS632PTVTTRQKxorQ8Ss2b/MAymdqLLy7dXK8vO27FuePUsU+mZVjVJCPpbEn3SHos+3lWRTmnayR9QNIjmZtuk3Rm6bjXSHpe0o/71KezVCXNMdgg62pgK3CdpK2lYlcDW7LbduBDAccmAvAdTtWGUMktE6njsSvnCsuHWHW5fmJm2AHsM7MtwL7s8TIaXHMP8PVmdhnw/4D3lg7/IPDnvpWJEaleARw0s8fN7CRwC7CtVGYb8Ls24JPAmdn+3D7HJgIYZZ6wKMmyQFc8X04NsHIkQCLhyTZgd3Z/N/AOR5lK15jZ3dluqwCfBC7KD5L0DuBx4IBvZWJI9ULgycLjQ9lzPmV8jgVA0nZJ+yXtP/bM6c6VTowfp9enFP8Msyn//85u2wOOPc/MjgBkP891lPF1zQ+SRaWSNgI/BfyXgLpE6f13JZXKbbeqMj7HDp402wXsAnj9N3juE5wYGdqwwRmtakPWWbLxDOdxCzVibVquMKUHxgwLyjkfM7P5qhcl3Quc73jpJs/zN7pG0k3AAvAH2VP/BfigmT0v+efOY0j1EHBx4fFFwGHPMms9jk30QJshRaHiKot1SagZi+vdK5i4ZlUlZhsze2vVa5KekrTZzI5kacWjjmK1npJ0PfBvgLeYWS7bbwK+W9L7gTOBRUnHzezX6+oao711H7BF0qWS1gLXAntKZfYA785GAbwB+HIWpvscm2iJ74LNXVnYsJrFjWtflmQhCtWGDUu3JQqvL26sFmtIx9q4LRWYGCp7gOuz+9cDtzvKVLpG0lUMmvlvN7OlKMDM/qWZXWJmlwC/AvzXJqFCBKlmCd4bgbuAh4FbzeyApBsk3ZAV28sg2XsQ+C3gP9Yd27VOoUzqvOniUnVdhFEnpPy85ah24QyxsH7VsjxolViXkT1fjFIXNqye6XxqcTZVGvjfip3AlZIeA67MHiPpAkl7odE1vw68ErhH0mckfbhLZaLMqDKzvQzEWXzuw4X7Bvyw77HDYNPq59KKPi0YSFbLhjwNhLh6aVbV4vq1g5EAFWLNhbq4ce2y2VQL6wfTVIvN/75GM8Q67+kRpfdfveaFkVx3HDGzZ4C3OJ4/DFxTeOx0jZl9tcc1fta3PmmaaqI1A/mtWhoGtbBhtfcHqijU0+tXreigKkfeqenenbPn3MPcEnFJUk00kndQ5fPqi9FqLlYYjDNd2DCIWPNcaXHGVDF/6hJq2yi1ywLURVkPa+V/X85c65cKGNf0lRZHvw3NKJiKRNamVZMVxqxa9/I4W1tfPeb29BntmpZ1i09XCci3E6t47mUCzPKrp9evYmHD6qXb4sa1Szdg6flB3dxCjRWlzuIQq02rnxt1FWaeFKn2wPp1p5ZtqbJm/YLX8n+n11nwXPLFtS+vAdDnWM3yuZei12XDnwapgLzTKY9cl5/n5e/xOqG2iVL7pu2X3AmQdD0AABy4SURBVJr1C82Fhsz5c2O2WO0UMVNSPXf1c8Fbqpy19oWJX/0/Fq40QC7WciqgimLuNBdpnVBDolTfpn8fnV/FFkexJZKYPaai+R8b357V4vCXvhYZbhsd1dElBVAW39LjLNpcdsua98Wbq6zrvMXn2uCzLXfXfGrbnv82O6metTa8t//clAoYCTMVqdaxac3zU7WuanGx5r7SAss7rpanA16OXB3HFZ5vuxFhfr2ccdt6xZdxGKM6aX0S487UR6qxh5H49sg2UdVZ1XXco290FyNahYEIi1FrbeRaikyL0WmsCNX1HmLsuNpHiyExncx0pDrMCQC+nVVtKHZWFWkbrfrs21Q+X/laPrgE1/VLoQ1dh1LV5VPbdlLVfXkX01PjOpwKAJv+3RRcTGWkGrtns00+q03erIpYUZLv8KqQ87mi1nLkWaTq9fK56qj7cvDJpdY978uoZlLVUTecKg38Hx4TKdVVa1436io4idVZ1eYfNsY6AG3GrebXc12zLNkqkXZZOCXky2Acmv4xv2wT48lEStWFb7Ldt0e0zQiAUHwnAXTBJbIidZIK7dwqSrKqt72NSH3q6nrdJ0rtcxZVU9Pf90u4TUupjjRGtV+mRqox8M1Pte2sGuYg8BidM9Bt1EBXidbVI1SoMcemFlsSscanxu759w0eXMHIGZs/H7Uus8ZMSLUun9TntL4uTb1yCsCnCeobdYVEq/nro5zy2aXJD9UjC3x+X6Po9fftpEqMJzMh1RjUNcFiTQKInQKoi1ZDxeoqMwza1KvPFa2qolQfYudTfVtWqZNquEysVJuaKHV5o9h51RDKKYCQJmObaLWrWEcVtbqu00aow4pSYw2lKuObTx3HhVS0+PLfsek2TUysVPsidl41NDqJMRGgS+eL6wNeJbM+/hmqzhtTqONAnzOp6oKG1EnVP1Ml1boRAH03gbqkAGJHqyuOCYhWoVpqfcq1Tqaxm/xto9Q+OqjKDDOfmjqp+qGTVCWdLekeSY9lP89ylLlY0v+R9LCkA5J+pPDaz0r6YrYvzGckXVM+fhj4Np1iDW1pahr2Ea22EatvxFgsHyLYpmN8JV41Rtb1+jAWonb9fVM+tT98PJSVu0rSo5IOStpReP7nJT2YOehuSRcUXntvVv5RSf/apz5dI9UdwD4z2wLsyx6XWQB+zMy+DngD8MOSthZe/6CZXZ7dou5VFSOvWkddVFGMVn3+ofqIVusE4iNWCIsei8d0yaP5RqdV9e7a7O87Sq1q+pc/T5OcTx0yjR6SNAfcDFwNbAWuK3joA2Z2mZldDvwZ8L7smK0Mdl19LXAV8BvZeWrpKtVtwO7s/m7gHeUCZnbEzD6V3f8nBjsZXtjxukB/TZViFBDS5ArJjXWNVrumASBMrHVyjbFCVNO5Ygh1GEOoRr0g9YzmUxs9BFwBHDSzx83sJHBLdhxmVvylbQTyD8E24BYzO2Fmn2ewG/QVTZXpKtXzzOxIVrEjwLl1hSVdAvxz4G8KT9+Yhd4fqQrbs2O3S9ovaf/TTz9deY22edW23/YxO6zqop4+0gDgL1ZobqqXb77lfKLeMsMUatcotepvH/IlHHsRlTHMp27K/7+z2/aAY308dCHwZOHxIQrBnaRfkPQk8H1kkWrTMVU0SlXSvZIecty2NR1bOs8rgP8J/KfCN8OHgK8CLgeOAL9UdbyZ7TKzeTObP+ecc7yvGzsF0NQkC+mwColWXbRNA3QRK/jnTkPkGXKNUQm1CZ8oteoz0UfTf9T5VFn1l6jjc3Es///ObruWnau7h1yL+y79cc3sJjO7GPgD4EafY6poXIvOzN5a9ZqkpyRtNrMjkjYDRyvKrWEg1D8wsz8pnPupQpnfYpDPGAuKi1a/es0LPHOq+5Yq5b2rXKxad5rFE+60jWsPq9NnmNeOleXlAeuW7stfB/8VoWINWaq7XpNMXWViNvljRqkh+KagprnpH8FDh4CLC48vAg47yn0UuAP4mYBjltG1+b8HuD67fz1we7mAJAH/HXjYzH659NrmwsN3Ag+FVsDVZOk7BVCOInw7rNpQjla75FddEWvVgtPFMj64OqBcgmzTeVW3ClZdXbsKtc8ota7pXxelFpv+qYNqiUYPAfcBWyRdKmktgw6oPQCSthTKvR14pHDeayWtk3QpsAX426bKdJXqTuBKSY8BV2aPkXSBpLwn/43A9wPf7hg69X5Jn5X0IPBm4Ec71sfJsGdXde2wKkdBPtMhu4wI8EkHhC7Plwuyy8yZuuvGEGoIbRai9o1SY+0mUaQ2eBi/fGpXGj1kZgsMmvV3Megsv9XMDuTHZ6mEB4G3AT+SHXMAuBX4HHAn8MNm1vjP2GkpejN7BniL4/nDwDXZ/b/CnZvAzL6/y/Xr2LRqPccWu31Y6/atKu+yeuba4/zjSbcBNq47yQsnBv/lrhSAa1eA0DQAtE8FQHM6IC+T0+dW2HX47BbgvbCMZ5TaZk2GslDbRqm+HVSzvMmfj4eyx3uBFcM2zey7as79C8AvhNRnKmZUhX7Ltk0B9BWt+uCTBoCwiLVNOqBcNlYutelcvrsFxBZqmdhz/NtGqSFN/0nPp04aUyFVH9qmAOqig7a5VVezsE0aIMaWHr7pAF/Blo+te85HzFXXdcm0D6EOM0ot07aDasaa/mPHVEu1jw6rWNFqLLG6CB3A7hO15jQJtnhs8fjQyLZOpl3ypyG/m/Lv2idKrRNqmZBhVG07qEYapaZVqiYbn2/bth+wPqLVWHRNAxSpilrrOouaBOtL06aBVfWIKdS+m/19RKmJ8WNqpNqVchNqkqLV2GKtk6uPYKvkWFWmTd42pLkP4UKNMSa1jyi1TGr6jx9TL9VYywH2Fa2Oo1ihXlohTfnQaLbp3G1k2kWoLtqsQlW3cEqZui/tiWn6zzBTJdWuKYCQYSl9jwQYB7FCs8BCOp26HBsqU/B7301CbdvsD+mcGlWUmuiHqZJqFbE6rLpEq3VrAnSZyjgssRZvTbhk2Ua+Iddcdn2P6BTChVpFzGZ/rCjVh9T074eZkGqZkGg1Vm61TKw0AISJNcbuoKGSHeZ5fd9fG6G26e3vMl65S5Q6Dk1/Gax+ybxu08TUSbXvb99Y0SqMRqwQT645ZRm2vXUh5D31JdQy5b/3qKLU1EE1XKZOqlWEpAC6RKuhYvWhD7FCt5TAOBEi01hCdRHS7G+iy5TUcYhSZ5mplGqfY1Zh5Qe+7zQAdBPrMKPWYZHXu210Ct2E2rXZHzNKTR1U48VUSrWKmNHqsNMAUC1Wn1lXPlHrpMg1tJ6u6HSYQh23KDU1/ftlaqXad7RapmsaoItYwT3cKjRqhfGWa2jdfJr70E2oZXyEGrISVYpSJ4+plWoVw4pWfehbrNAuaoXlzetRS7ZNHXyb+12FGjr1uEuzP0aPf4pS+2eqpdpHtFon1tBo1YdYYm0TtS4rP0TJtr1W/p7aNvehm1D7bPaHMg5RqhbTkKqZoUu02kTsNAC0E2tI1NpmCcGy+HwlmL/e9njnOWveg29zH/oXaopSZ4Opl2qbaLVrGqBpNMAwxAr+UStUR3ltaBJmrEi3Saa+zX0YvVDrlvYL7Zwahyh1lukkVUlnS7pH0mPZz7Mqyj2R7UX1GUn7Q4/vg64fvC6jASC+WNtErVULh8SSa1+EyhTqo9NhC7VMaLO/zdbT0xylBnjoKkmPSjooaUfh+Z+X9GDmp7slXVB47TJJ/1fSgcxhjeLoGqnuAPaZ2RZgX/a4ijeb2eVmNt/y+Nb0Ea020ZQGcOESa8g6AU1Ra1u5joNgm+pSJ9PQ6BT6FWpIb39os39Go9RGj0iaA24Grga2AtdJ2pq9/AEzu8zMLgf+DHhfdsxq4PeBG8zstcCbgMaFOrpKdRuwO7u/G3jHkI/vRPkD2HcawKfjyjVoPFSsseUKoxGszzXbynQchFomdCqqr1CnOUrN8PHIFcBBM3vczE4Ct2THYWbFX/xGIP/AvQ140MweyMo947ObalepnmdmR7ILHgHOrShnwN2S7pe0vcXxSNouab+k/U8//XRwRfv6YI1SrG2jVmiWq69go+RgS+fzFWlVHetWl2r6vQxTqE3N/pCl/aqYIKFuyv+/s9v25kOW8PHIhcCThceHsucAkPQLkp4Evo8sUgW+BjBJd0n6lKSf9KlM4xbVku4Fzne8dJPPBTLeaGaHJZ0L3CPpETP7WMDxmNkuYBfA/Px8tHCpvJX1+XMn+IfTL6+ofPbcizx7esPS43NXP8fRhVctP8fq5zhWeK5ua2sXru2tX7XuOM+dWP5ccavrHNeW1zm5QMrbXxfJBeTaDrsoLR13b5cNfuNeu9K0cHTTMn1tZArDE+o0Nvu1aKw+vuhb/FgpNbj8XN095Nq7femDa2Y3ATdJei9wI/AzDPz4rcA3Ai8C+yTdb2b76i7UKFUze2tlLaWnJG02syOSNgNHK85xOPt5VNJtDELxjwFex8fijM2f56UjlzaWayPWOl695gWeObVx6fFZa1/gSyc3LitTJVZgmVxDxQrd5Qr+go2F7y6mXWUK4y/UMrPY7I/goUPAxYXHFwGHHeU+CtzBQKqHgL80s2PZdfYCr2eQt62ka/N/D3B9dv964PZyAUkbJb0yv88gT/GQ7/GxcX3Q2nzLx86vwuCf1CcdUDUyoGkKpY9g8rRAnayKze+mVEETrnP5RKVNdWzKm0L178zVQTgsobpo0+yfMXw8ch+wRdKlktYC12bHIWlLodzbgUey+3cBl0nakHVafRvwuabKdJXqTuBKSY8BV2aPkXRBZnWA84C/kvQA8LfAHWZ2Z93xoyC00wr6ESv45Vld//jQPDfdRzY5TeIqUiXHppsvPiLNaStTcH9hDSuHCvGa/dMUpXrQ6CEzW2DQrL8LeBi41cwO5MdLekjSgwyCvh/JjvkS8MsMhPwZ4FNmdkdTZWQ2+uEyoczPz9v+/fubC9bgSgMUc6s5xTQAsCwNADjTAMdKz5Xzq8VUALAiFZBTTgcAK/KswIp0QE5dSqBIXVrARVWKIBa+Ii/i+0URIlPoV6h95VFjCDXLHVbmOH145ZkX2T//Vz/iVfbj/+snO19vXJj6GVVV+KYB2oxf7Tti9UkHgP/eVyHRKyyPGsu3YZ4jpO5totM+hk3lhOZRE5PDzEq1Cp/8al9i7ZpnbZtrzckFFSLYMnWybCvOLvVskqlPcx+6CbVNx1Rq9k8uMy3Vqg9gU37VRQyxQrc8K9RHrSG7thbF1UWysQitS1uZDluoLqZGqIswd3zR6zZNzLRUob1YfTqunOftQay+USuEyzWnLNk+RdvlWk3vzzd3CtWthJhCnYbxqInlhPVQzDhN41dh5RjW8sQAWDk5oDyOFV7+x3WNZ4WVnVhVkwXA3ZFVFI9vh1aZcYhgc5q+KEJkCtVrocbqlIK420yPXZQ6w8x8pAr+0SrEGWoF7oi1j6gV6iNXaB+9jpq83k2RaWhTvyo6HaZQXUxMs3/GSVLNCBFrGa9/CA+xQnU6ILQTq4tcx1myvvULlSmER6exhOoiNfsnlyTVAl06rtqMCAB/sUJY1ArNcm3aX6ks2VGINuT6de+rSaZdmvvQTahd8qgpSh0/Uk7Vk6aFV6Dd4ivgXoAl/8ftmmsF9xoCOUUBVU0iKFIntra52abzNlH35VAlUqjfQ6rqCywJNdFEkmqJukVXYooVWLGyFaycfeXqxAL3oizQXq4QLtgyw4pkfXYwjS1TSEINRYvG6hfHpzNzWCSpOvBdzQr8xQorp7SGRK1VYgX3NFcfuYKfYKGdZGPhuw10nUiheXfTPqJTaCfUOiZBqLNMkmoFVWItR6vgJ1boJx0A7eUKfoIFt9j6EK2vQIv0JVMYnVBTT//kkqRaw6jFCv7pAKhOCcBysfgIFuolC+0EGIMmiebElil07+HvItTEZJCk2sAoxQpxo9YcH8GCfxTbJ74CzWkSKTRvwthHdArdhZqi1MkgSdWDYYoVVi4dWBe1Qnu5QjvB5sQWbahAc3xECu1kCkmoiTCSVDvSVazg14EF1Xtf+cgV4gk2p60EYxBLpBAmU+iWP4XZEaoWjVUvjCZFNEqSVD0JGWoFL//jlNcKgJULXceIWqE+3wr+0SuslJaPZPvCV6BFusgU+olOYXaEOsskqQbQJFZYuXtAjKgV4qQEckKi15xhSLaNPIt0FSnEiU4hCXWW6SRVSWcDfwRcAjwBfE+2r0uxzNdmZXL+GfA+M/sVST8L/Dvg6ey1/2xmexljmsawdkkHQPVOrXUpAaiWK8QXLHQXYCx8RAqjlSkkofaJj4eyclcBvwrMAb9tZjtLr/848AHgHDM7Jinf72otcBL4CTP730316Tr3fwewz8y2MNi2dUe5gJk9amaXm9nlwL9gsH/2bYUiH8xfH3eh5jR94H1Wt4LBP2DVKldVawdU/ZNvWvN8pRiqVsAqky/cUrWAy7jgW8f8fTc186ua+rGiU9dMqSTUqDR6SNIccDNwNbAVuE7S1sLrFzPYNPALhcOOAd9hZq9jsEvr7/lUpmvzfxvwpuz+buAvgJ+qKf8W4O/M7O87Xnfk5B/8rnlWaBe1wsqUADTnXHPqotecsrRCItlYtJF726gU6lfm7zs6hSTUlvh46ArgoJk9DiDpluy4fMvpDwI/SWF7azP7dOH4A8B6SevMrHb6W1epnmdmR7IKHJF0bkP5a4E/LD13o6R3A/uBH3OF7eNMaAcWVKcDwN2JBe5dW9vKFcIFC+0EV6Yo5pjRsNdme5FlCkmoY4KPhy4Eniw8PgR8E4CktwNfNLMHJFVd47uATzcJFTykKule4HzHSzc1HVs6z1rg7cB7C09/CPh5wLKfvwT8YMXx24HtAK95zWtCLt07sTqwoD5qhfZyhbiCbcs4iRSa94yKEZ3CjArVjFXHvYdUbZJU3Hd+l5ntyh9E8JDLliZpQ3aOt1UeKL0W+MW6MkUapWpmb6252FOSNmffDpuBozWnuhr4lJk9VTj30n1JvwX8WU09dgG7AObn562p3sOmbQcWuNMBsDJqheqUANTLFZqjV1gpqr4lG4KPRKFZpNBeppCi0544ZmbzVS9G8NAh4OLC44uAw8BXAZcCeZR6EfApSVeY2T9IuohBH9C7zezvfN5I1+b/HgYJ3J3Zz9tryl5Hqemf/yKyh+8EHupYn5HSRqwQN2qF5cJoG73mVImsT9n6yrNIDJFCPJlCEuoQ8fHQfcAWSZcCX2SQivxeMzsALKULJD0BzGe9/2cCdwDvNbNP+Famq1R3ArdKeg+DXrN3ZRW7gMGQhWuyxxsY9Kz9+9Lx75d0OYPm/xOO1yeOLmKF8Kg1p2v0mtMk2Zw24ouNj0ihP5lCEuqY0OghM1uQdCNwF4MhVR/JhFrHjcBXAz8t6aez595mZnUtcmQ2di3pRubn523//v3NBUdM05qsLrnCSrEWccm1SJVcl13Xo8xSWU/JDgNfiUJ3kebEjk5hMoQq6f665rgPX3HGZvuWS/6tV9k7H9nZ+XrjQppR1SOxo1aoj1whLHqFZsFWiaxP2YbIc9lxHpKEbjKFFJ0m6klS7RkfsYI7aq3KtUKzXKE59worReQbxbYVX0x8JZozKplCEuoskaQ6BJomCkD90CuoTgmEyBWa0wNtJds3oQLN8REptJMpJKHWcnoRXnhp1LUYOkmqQ8Rn76s2KQHwkyuECRaqZdaXbNvKs4ivSCHJNBGfJNUh0yVqBX+5QnzBLqtjBPnFJJZIobtMIQl1lklSHRG+USvUyxW6R6+wUkqhkh0WIfIsk2SaGAZJqiPEJ2qFerlC3Og1p0pew5JtF3kWaRIpNG8PnZr6iRCSVMeAYckVVkrGV7I5IbIrCjiWJH3wESnEkSkkoSaWk6Q6RvikBMBfrlAvWOgu2TqGIVJfgebEEikkmSbcJKmOGb5RKywXQJfotYhLUjFF24VQgeY0iRSSTHthcRF7sd3fbJJJUh1TQuQKYdEr+EsWqmXWl2zbyrOIj0ghyTQRnyTVMSe2XHO6SDYnhvxi4StRCBMpJJkmwkhSnRDayhWaBQtxJDssQgSaEypSSDJNtCNJdcIIlSusFEobyeYMW7ZtBJrTRqSQZJroRpLqhFL8xw8RLIRHsUV8JOcr3i7CdNFWopBEmohHkuoU0CZ6zWkTxTYRW5Yuugi0SJJpIjZJqlNEl+g1p0pWMWTbhVgShSTSobG4iD0/+h0ihk2S6pQSQ7BF6qQWQ7gxpekiiTQxLJJUZ4CyUGJItkjfQmxDkmhiVKwadQUSw+eMzZ9fdpsGpvE9JfyQdLakeyQ9lv08q6LcVZIelXRQ0g7H6z8uySRtyh6vkbRb0mclPSzpvT71SZFqwimh2NFsbJI4EwV2APvMbGcmyx3ATxULSJoDbmawq/Mh4D5Je8zsc9nrF2evfaFw2LuAdWb2umxH6M9J+kMze6KuMkmqCSdV0hqmbJM4E55sA96U3d8N/AUlqQJXAAfN7HEASbdkx30ue/2DwE8CtxeOMWCjpNXAGcBJoHGVoImU6v33339M0t+Puh49sAk4NupK9ESL96ZeKhKZaf2bfW3XEzxnz95194mPbvIsvl5Scd/5XWa2y/PY88zsCICZHZF0rqPMhcCThceHgG8CkPR24Itm9oC07DP3xwzEewTYAPyomT3bVJmJlKqZnTPqOvSBpP3Tsvd5mWl9b9P8vrqew8yuilEXAEn3Auc7XrrJ9xSO5yxr1t8EvM3x+hXAaeAC4Czg45LuzaPdKiZSqolEYrYws7dWvSbpKUmbsyh1M3DUUewQcHHh8UXAYeCrgEuBPEq9CPiUpCuA7wXuNLNTwFFJnwDmgVqppt7/RCIx6ewBrs/uX8/yvGjOfcAWSZdKWgtcC+wxs8+a2blmdomZXcJAvq83s39g0Gn17RqwEXgD8EhTZZJUxwvfHNIkMq3vLb2v0bMTuFLSYwx68HcCSLpA0l4AM1sAbgTuAh4GbjWzAw3nvRl4BfAQAyn/jpk92FQZmVnbN5JIJBKJEilSTSQSiYgkqSYSiUREklRHSMD0uieyqXKfiTHUpS88pgFK0q9lrz8o6fWjqGcbPN7bmyR9OfsbfUbS+0ZRz1AkfUTSUUkPVbw+sX+zUZGkOlry6XVbgH3Z4yrebGaXj+uYyMI0wKuBrcB1kraWil0NbMlu24EPDbWSLfF8bwAfz/5Gl5vZzw21ku35H0DdeNKJ/JuNkiTV0bKNwbQ6sp/vGGFdurI0DdDMTgL5NMAi24DftQGfBM7MxhWOOz7vbSIxs48BdbOEJvVvNjKSVEfLsul1gGt6HQzmIN8t6X5J24dWuzBc0wAvbFFmHPGt9zdLekDSn0t67XCq1juT+jcbGWlGVc9EmF4H8EYzO5zNab5H0iNZhDFOOKcBtigzjvjU+1PAV5rZ85KuAf6UQZN50pnUv9nISFLtmQjT6zCzw9nPo5JuY9AcHTepVk0DDC0zjjTW28yeK9zfK+k3JG0ys0lfbGVS/2YjIzX/R0vj9DpJGyW9Mr/PYOEHZ0/tiHFOAyyV2QO8O+tRfgPw5Tz9MeY0vjdJ5yubPJ7NG18FPDP0msZnUv9mIyNFqqNlJ3CrpPcwmGf8LhhMrwN+28yuAc4Dbsv+X1cDHzWzO0dU30rMbEFSPg1wDviImR2QdEP2+oeBvcA1wEHgReAHRlXfEDzf23cD/0HSAvAScK1NwHRFSX/IYC3STZIOAT8DrIHJ/puNkjRNNZFIJCKSmv+JRCIRkSTVRCKRiEiSaiKRSEQkSTWRSCQikqSaSCQSEUlSTSQSiYgkqSYSiURE/j9hGtDDJqxeAwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "c = plot(wdot)\n",
    "import matplotlib.pyplot as plt\n",
    "plt.colorbar(c)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The video generated from [Paraview](https://www.paraview.org) is also attached. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}