<jupyter_start><jupyter_text>1D Elastic Wave Equation with a Discontinuous Galerkin-FEM in the Time Domain with Free Surface Boundary Conditions  Notebook setup<jupyter_code>import numpy
from matplotlib import pyplot
%matplotlib inline<jupyter_output><empty_output><jupyter_text>Governing Equation  Displacement FormulationConsider the elastic wave equation in 1D written in its displacement form:$$\rho\partial_{tt}u = \partial_{x}(\mu\partial_{x}u) + f $$with $x \in \Omega = [x_{min},x_{max}]$ and $t \in [0,t_{max}]$, and where $u(x,t)$ is the displacement field, $\rho(x)$ is the material density, $\mu(x)$ is the material shear modulus, and $f(x,t)$ is a source term. The displacement field is subjected to free surface (or stress free) boundary conditions\begin{align}\left.\mu\partial_{x}u(x,t)\right|_{x=x_{min}} &= 0, \\ \left.\mu\partial_{x}u(x,t)\right|_{x=x_{max}} &= 0,\end{align} and initial conditions$$\left.\partial_{t}u(x,t)\right|_{t=0} = \left.\partial_{tt}u(x,t)\right|_{t=0} = 0.$$  Velocity-Stress FormulationFor a Discontinuous Galerkin (DG) approach, it will be convenient to rewrite the displacement form of the elastic wave equation in a velocity-stress formulation in order to reduce the time and space derivatives to first order. For that, we define two other variables namely velocity $v$ and stress $\sigma$:\begin{align}&\partial_{t}u = v, \\ &\sigma = \mu\partial_{x}u.\end{align}Upon substitution into the displacement form, these variables will give rise to the velocity-stress formulation of the elastic wave equation\begin{align}\partial_{t}v &= 1/\rho\:\partial_{x}\sigma + f/\rho, \\\partial_{t}\sigma &= \mu\partial_{x}v.\end{align}or, in vector notation$$\frac{\partial}{\partial t} \left[\begin{array}{c}v \\\sigma\end{array}\right] = \left[\begin{array}{cc}	  0 & 1/\rho \\	\mu &   0\end{array}\right] \frac{\partial}{\partial x} \left[\begin{array}{c}	v \\	\sigma\end{array}\right] + \left[\begin{array}{c}	f/\rho \\	0\end{array}\right]$$ As before, this system of PDEs will be endowed with free surface boundary conditions\begin{align}\left.\partial_{x}v(x,t)\right|_{x=x_{min}} =& \ 0, \ \ \left.\sigma(x,t)\right|_{x=x_{min}} = 0, \\ \left.\partial_{x}v(x,t)\right|_{x=x_{max}} =& \ 0, \ \ \left.\sigma(x,t)\right|_{x=x_{max}} = 0,\end{align} as well as initial conditions$$\left.v(x,t)\right|_{t=0} = \left.\sigma(x,t)\right|_{t=0} = 0.$$  Discontinuous Galerkin-FEM DiscretizationWe approximate $\Omega$ by $K$ nonoverlapping elements, $x \in [x^{k}_{l},x^{k}_{r}] = \Omega^{k}$. On each of these elements we express the local displacement solution as a polynomial of order N$$x \in \Omega^{k}: u^{k}_{h}(x,t) = \sum^{N+1}_{j=1}\hat{u}^{k}_{j}(t)\psi_{j}(x) = \sum^{N+1}_{i=1}u^{k}_{h}(x^{k}_{i},t)\ell^{k}_{i}(x).$$In the last statement, we have introduced two complementary expressions for the local solution known as the modal form, which uses the local polynomial basis, $\psi_{j}(x)$, and the nodal form, which is defined through $N_{p} = N+1$ local grid points, $x^{k}_{i} \in \Omega^{k}$, and uses the interpolating Lagrange polynomials, $\ell^{k}_{i}(x)$, as basis functions. These two representations are connected by the definition of the expansion coefficients, $\hat{u}^{k}_{j}$.From now on, we will choose the modal representation for the local solution but we should keep in mind that any of the two are mathematically equivalent (although not computationally!) and can be consistently used to discretize the elastic wave equation.The global solution $u(x,t)$ is assumed to be approximated by the piecewise N-th order polynomial approximation $u_{h}(x,t)$,$$u(x,t) \approx u_{h}(x,t) = \bigoplus^{K}_{k=1}u^{k}_{h}(x,t),$$defined as the direct sum of the $K$ local polynomial solutions $u^{k}_{h}(x,t)$. In exactly the same way, for the velocity-stress formulation of the elastic wave equation, we express the local velocity and local stress solutions as\begin{align}v^{k}_{h}(x,t) &= \sum^{N+1}_{j=1}\hat{v}^{k}_{j}(t)\psi_{j}(x), \\\sigma^{k}_{h}(x,t) &= \sum^{N+1}_{j=1}\hat{\sigma}^{k}_{j}(t)\psi_{j}(x).\end{align}Going back to the differential form of the velocity-stress formulation, we proceed as in the Galerkin projection forming the residual  $$	R_{h}(x,t) = \frac{\partial}{\partial t} \left[	\begin{array}{c}		v_{h} \\		\sigma_{h}	\end{array}	\right] - \left[	\begin{array}{cc}		  0 & 1/\rho \\		\mu &   0	\end{array}	\right] \frac{\partial}{\partial x} \left[	\begin{array}{c}		v_{h} \\		\sigma_{h}	\end{array}	\right] - \left[	\begin{array}{c}		f/\rho \\		0	\end{array}	\right]$$and requiring it to be orthogonal to the space of test functions $V_{h} = \oplus^{K}_{k=1}V^{k}_{h}$ where $V^{k}_{h} = span\{\psi_{j}(\Omega^{k})\}^{N+1}_{j=1}$$$\int_{\Omega^{k}} R_{h}(x,t)\psi_{i}(x) dx = 0, \ 1 \leq i \leq N+1,\label{orthog-condition}$$on all $K$ elements. This yields exactly $N+1$ equations for the $N+1$ unknowns on each element. Integrating by parts once for the velocity field and twice for the stress field we get the local weak form of the elastic wave equation\begin{align}\int_{\Omega^{k}} \left(\frac{\partial v^{k}_{h}}{\partial t}\psi_{i} + \frac{1}{\rho}\sigma^{k}_{h}\frac{\partial \psi_{i}}{\partial x} - \frac{f^{k}}{\rho}\psi_{i}\right) dx &= \int_{\partial \Omega^{k}} \mathbf{\hat{n}} \cdot \left(\frac{1}{\rho}\sigma^{k}_{h}\psi_{i}\right) dx = \left[\frac{1}{\rho}\sigma^{k}_{h}\psi_{i}\right]^{x^{k}_{r}}_{x^{k}_{l}},  \ 1 \leq i \leq N+1, \\\int_{\Omega^{k}} \left(\frac{\partial \sigma^{k}_{h}}{\partial t}\psi_{i} + \mu v^{k}_{h}\frac{\partial \psi_{i}}{\partial x} + \frac{\partial \mu}{\partial x} v^{k}_{h}\psi_{i} \right) dx &= \int_{\partial \Omega^{k}} \mathbf{\hat{n}} \cdot (\mu v^{k}_{h}\psi_{i}) dx = \left[\mu v^{k}_{h}\psi_{i}\right]^{x^{k}_{r}}_{x^{k}_{l}},  \ 1 \leq i \leq N+1.\end{align}At this point, using the standard Discontinuous Galerkin (DG) approach, we replace the physical fluxes represented by the line integrals computed along the interfaces between elements by some numerical fluxes in order to recover the global solution from the $K$ local solutions. This leads to the following semidiscrete scheme\begin{align}\int_{\Omega^{k}} \left(\frac{\partial v^{k}_{h}}{\partial t}\psi_{i} + \frac{1}{\rho}\sigma^{k}_{h}\frac{\partial \psi_{i}}{\partial x} - \frac{f^{k}}{\rho}\psi_{i}\right) dx &= \left[\frac{1}{\rho}(\sigma_{h})^{*}\psi_{i}\right]^{x^{k}_{r}}_{x^{k}_{l}},  \ 1 \leq i \leq N+1,\\\int_{\Omega^{k}} \left(\frac{\partial \sigma^{k}_{h}}{\partial t}\psi_{i} + \mu v^{k}_{h}\frac{\partial \psi_{i}}{\partial x} + \frac{\partial \mu}{\partial x} v^{k}_{h}\psi_{i} \right) dx &= \left[\mu (v_{h})^{*}\psi_{i}\right]^{x^{k}_{r}}_{x^{k}_{l}},  \ 1 \leq i \leq N+1.\end{align}Now, we have a total of $K\times(N+1)$ equations for the same number of unkowns. It is worth noting that the numerical fluxes are responsible for recovering the global solution from the local solutions as well as imposing (weakly) the boundary conditions.The expressions above are known as the DG weak form of the velocity-stress formulation of the 1D elastic wave equations.  To get the matrix form of the semidiscrete scheme, we continue by substituting the approximations for the velocity and the stress into the DG weak form\begin{align}\sum^{N+1}_{j=1} \frac{d}{dt}\hat{v}^{k}_{j} \int_{\Omega^{k}} \psi_{j}\psi_{i} dx &+ \frac{1}{\rho} \sum^{N+1}_{j=1} \hat{\sigma}^{k}_{j} \int_{\Omega^{k}} \psi_{j}\frac{\partial \psi_{i}}{\partial x} dx - \frac{1}{\rho}\int_{\Omega^{k}} f^{k}\psi_{i} dx = \frac{1}{\rho}\left[(\sigma_{h})^{*}\psi_{i}(x^{k}_{r})-(\sigma_{h})^{*}\psi_{i}(x^{k}_{l})\right],  \\\sum^{N+1}_{j=1} \frac{d}{dt}\hat{\sigma}^{k}_{j} \int_{\Omega^{k}} \psi_{j}\psi_{i} dx &+ \sum^{N+1}_{j=1} \hat{v}^{k}_{j} \int_{\Omega^{k}}\mu \psi_{j}\frac{\partial \psi_{i}}{\partial x} dx + \sum^{N+1}_{j=1} \hat{v}^{k}_{j} \int_{\Omega^{k}} \frac{\partial \mu}{\partial x} \psi_{j}\psi_{i} dx = \mu(x^{k}_{r})(v_{h})^{*}\psi_{i}(x^{k}_{r})-\mu(x^{k}_{l})(v_{h})^{*}\psi_{i}(x^{k}_{l}).\end{align}where $\rho$ is assumed to be constant. These equations can be expressed in matrix form in the following way\begin{align}\hat{\mathcal{M}}^{k}\frac{d}{dt}\hat{\boldsymbol{v}}^{k}_{h} &+ \frac{1}{\rho}\hat{\mathcal{S}}^{k}_{v}\hat{\boldsymbol{\sigma}}^{k}_{h} - \frac{1}{\rho}\hat{\mathbf{f}}^{k} = \frac{1}{\rho}\left[(\sigma_{h})^{*}\boldsymbol{\psi}(x^{k}_{r})-(\sigma_{h})^{*}\boldsymbol{\psi}(x^{k}_{l})\right],\\\hat{\mathcal{M}}^{k}\frac{d}{dt}\hat{\boldsymbol{\sigma}}^{k}_{h} &+ \hat{\mathcal{S}}^{k}_{s}\hat{\boldsymbol{v}}^{k}_{h} + \hat{\mathcal{S}}^{k}_{m}\hat{\boldsymbol{v}}^{k}_{h} = \mu(x^{k}_{r})(v_{h})^{*}\boldsymbol{\psi}(x^{k}_{r})-\mu(x^{k}_{l})(v_{h})^{*}\boldsymbol{\psi}(x^{k}_{l}).\end{align}where, for $1 \leq i,j \leq N+1$,\begin{align*}\hat{\mathcal{M}}^{k} &= \int_{\Omega^{k}} \psi_{i}\psi_{j} dx, \\\hat{\mathcal{S}}^{k}_{v} &= \int_{\Omega^{k}} \psi_{j}\frac{\partial \psi_{i}}{\partial x} dx, \\\hat{\mathcal{S}}^{k}_{s} &= \int_{\Omega^{k}}\mu \psi_{j}\frac{\partial \psi_{i}}{\partial x} dx, \\\hat{\mathcal{S}}^{k}_{m} &=\int_{\Omega^{k}} \frac{\partial \mu}{\partial x} \psi_{j}\psi_{i} dx, \\\hat{\mathbf{f}}^{k} &= \int_{\Omega^{k}} f^{k}\psi_{i} dx.\end{align*}are the local mass, velocity and stress stiffness matrices and source vector, respectively. Furthermore, we have\begin{align*}\hat{\boldsymbol{v}}^{k}_{h} &= [\hat{v}^{k}_{1},\ldots,\hat{v}^{k}_{N+1}]^{T}, \\\hat{\boldsymbol{\sigma}}^{k}_{h} &= [\hat{\sigma}^{k}_{1},\ldots,\hat{\sigma}^{k}_{N+1}]^{T}, \\\boldsymbol{\psi} &= [\psi_{1}(x),\ldots,\psi_{N+1}(x)]^{T}.\end{align*}as the vectors of local velocity and stress solutions and the local test functions, respectively. In order to perform the integrals over each element numerically, we transform coordinates by introducing the affine mapping$$x \in \Omega^{k}: x(\xi) = x^{k}_{l} + \frac{1+\xi}{2}h^{k}, h^{k} = x^{k}_{r} - x^{k}_{l},$$with the reference coordinate $\xi \in \hat{\Omega} = \left[-1,1\right]$. Now, the local polynomial representations can be rewritten as$$x \in \Omega^{k}: u^{k}_{h}(x(\xi),t) = \sum^{N+1}_{j=1}\hat{u}^{k}_{j}(t)\hat{\psi}_{j}(\xi) = \sum^{N+1}_{i=1}u^{k}_{h}(x^{k}_{i},t)l^{k}_{i}(\xi).$$and in exactly the same way for the approximations $v^{k}_{h}$ and $\sigma^{k}_{h}$. We should note that the Jacobian of the previous mapping is given by $$J = \left|\frac{dx}{d\xi}\right| = \frac{h^{k}}{2}$$and the derivatives of the local test functions are computed using the chain rule$$\frac{d\psi(x(\xi))}{dx} = \frac{d\psi}{d\xi}\frac{d\xi}{dx} = \frac{2}{h^{k}}\frac{d\hat{\psi}}{d\xi}.$$ After the development of a suitable local representation of the solution using the reference element $\hat{\Omega}$, we return to the discrete equations and rewrite them as follows\begin{align}\sum^{N+1}_{j=1} \frac{d}{dt}\hat{v}^{k}_{j} \ \frac{h^{k}}{2}\int_{\hat{\Omega}} \hat{\psi}_{j}\hat{\psi}_{i} d\xi &+ \frac{1}{\rho} \sum^{N+1}_{j=1} \hat{\sigma}^{k}_{j} \int_{\hat{\Omega}} \hat{\psi}_{j}\frac{\partial \hat{\psi}_{i}}{\partial \xi} d\xi - \frac{1}{\rho}\frac{h^{k}}{2}\int_{\hat{\Omega}} \hat{f}\hat{\psi}_{i} d\xi = \frac{1}{\rho}\left[(\sigma_{h})^{*}\hat{\psi}_{i}(1)-(\sigma_{h})^{*}\hat{\psi}_{i}(-1)\right],  \\\sum^{N+1}_{j=1} \frac{d}{dt}\hat{\sigma}^{k}_{j} \ \frac{h^{k}}{2}\int_{\hat{\Omega}} \hat{\psi}_{j}\hat{\psi}_{i} d\xi &+ \sum^{N+1}_{j=1} \hat{v}^{k}_{j} \int_{\hat{\Omega}}\mu \hat{\psi}_{j}\frac{\partial \hat{\psi}_{i}}{\partial \xi} d\xi + \sum^{N+1}_{j=1} \hat{v}^{k}_{j} \ \frac{h^{k}}{2}\int_{\hat{\Omega}} \frac{\partial \mu}{\partial x} \hat{\psi}_{j}\hat{\psi}_{i} d\xi = \mu(x^{k}_{r})(v_{h})^{*}\hat{\psi}_{i}(1)-\mu(x^{k}_{l})(v_{h})^{*}\hat{\psi}_{i}(-1).\end{align}or in matrix form\begin{align}\frac{h^{k}}{2}\mathcal{M}\frac{d}{dt}\hat{\boldsymbol{v}}^{k}_{h} &+ \frac{1}{\rho}\mathcal{S}_{v}\hat{\boldsymbol{\sigma}}^{k}_{h} - \frac{1}{\rho}\frac{h^{k}}{2}\mathbf{f} = \frac{1}{\rho}\left[(\sigma_{h})^{*}\boldsymbol{\hat{\psi}}(1)-(\sigma_{h})^{*}\boldsymbol{\hat{\psi}}(-1)\right], \\\frac{h^{k}}{2}\mathcal{M}\frac{d}{dt}\hat{\boldsymbol{\sigma}}^{k}_{h} &+ \mathcal{S}_{s}\hat{\boldsymbol{v}}^{k}_{h} + \frac{h^{k}}{2}\mathcal{S}_{m}\hat{\boldsymbol{v}}^{k}_{h} = \mu(x^{k}_{r})(v_{h})^{*}\boldsymbol{\hat{\psi}}(1)-\mu(x^{k}_{l})(v_{h})^{*}\boldsymbol{\hat{\psi}}(-1).\end{align}where, for $1 \leq i,j \leq N+1$,\begin{align*}\mathcal{M} &= \int_{\hat{\Omega}} \hat{\psi}_{i}\hat{\psi}_{j} d\xi, \\\mathcal{S}_{v} &= \int_{\hat{\Omega}} \hat{\psi}_{j}\frac{\partial \hat{\psi}_{i}}{\partial \xi} d\xi, \\\mathcal{S}_{s} &= \int_{\hat{\Omega}}\mu \hat{\psi}_{j}\frac{\partial \hat{\psi}_{i}}{\partial \xi} d\xi, \\\mathcal{S}_{m} &=\int_{\hat{\Omega}} \frac{\partial \mu}{\partial x} \hat{\psi}_{j}\hat{\psi}_{i} d\xi, \\\mathbf{f} &= \int_{\hat{\Omega}} \hat{f}\hat{\psi}_{i} d\xi.\end{align*}are the elemental mass, velocity and stress stiffness matrices and source vector defined over the reference element, respectively. The vectors $\hat{\boldsymbol{v}}^{k}_{h}$ and $\hat{\boldsymbol{\sigma}}^{k}_{h}$ are defined exactly as before. The vector of elemental test functions is defined as $$\boldsymbol{\hat{\psi}} = [\hat{\psi}_{1}(\xi),\ldots,\hat{\psi}_{N+1}(\xi)]^{T}.$$As usual, integration on the reference element $\hat{\Omega} = [-1,1]$ is performed using Gaussian quadrature $$\int^{1}_{-1} f(\xi) d\xi = \sum^{N_{i}}_{l=1} w_{l} f(\xi_{l}),$$where $w_{l}$ and $\xi_{l}$ are the quadrature weights and points, respectively, and  $N_{i}$ is total number of integration points needed to exactly integrate a polynomial of order at most $2N$ which is the polynomial order associated to the mass matrix coefficients. In the present DG approach, we considered two options for the distribution of quadrature points: standard Gauss-Legendre (GL) and Gauss-Legendre-Lobatto (GLL) points. The first option considers only internal points on the reference element while the second includes its boundaries $-1$ and $1$. Now we are in a position to choose the numerical fluxes $(v_{h})^{*}$ and $(\sigma_{h})^{*}$. The role of these fluxes is to guarantee stability of the formulation by mimicking the flow of information in the underlying partial differential equation. A simple choice that preserves stability in the case of the elastic wave equation is the Lax-Friedrichs flux\begin{align}(v_{h})^{*} &= \{\{v_{h}\}\} + \frac{Z}{2}[[\sigma_{h}]], \ \ \text{where} \ \ Z = -(\rho \: \mu)^{-1/2}, \\(\sigma_{h})^{*} &= \{\{\sigma_{h}\}\} + \frac{Y}{2}[[v_{h}]], \ \ \text{where} \ \ Y = -(\rho \: \mu)^{1/2}.\end{align}In the expressions above we have used the average and jump operators which are usually defined as follows$$\{\{u\}\} = \frac{u^{-}+u^{+}}{2}$$where $u$ can be both a scalar or a vector quantity and\begin{align}[[u]] &= \hat{\mathbf{n}}^{-}u^{-} + \hat{\mathbf{n}}^{+}u^{+}, \\[[\mathbf{u}]] &= \hat{\mathbf{n}}^{-}\cdot\mathbf{u}^{-} + \hat{\mathbf{n}}^{+}\cdot\mathbf{u}^{+}\end{align}where $\hat{\mathbf{n}}$ is defined as the normal vector along the left or right boundaries of the $k$ element. We observe that the jump operator is defined differently depending on whether $u$ is a scalar or a vector, $\mathbf{u}$. The superscripts "$-$" and "$+$" stand for the interior and exterior information of the element, respectively. It is worth noting that $$-\mu Z = -\frac{Y}{\rho} = c$$where $c$ is the propagation wave speed. Upon substitution of the expressions of the Lax-Friedrichs flux at the right hand side of the matrix equations, we come up with the final matrix form of the proposed DG scheme. Let us consider the numerical flux, $(\sigma_{h})^{*}$, first\begin{align*}\frac{1}{\rho}\left[(\sigma_{h})^{*}\boldsymbol{\hat{\psi}}(1)-(\sigma_{h})^{*}\boldsymbol{\hat{\psi}}(-1)\right] =& \frac{1}{\rho} \left[\left(\{\{\sigma^{k}_{h}\}\} + \frac{Y}{2}[[v^{k}_{h}]]\right)\boldsymbol{\hat{\psi}}(1) - \left(\{\{\sigma^{k}_{h}\}\} + \frac{Y}{2}[[v^{k}_{h}]]\right)\boldsymbol{\hat{\psi}}(-1)\right] \\=& \frac{1}{\rho}\left[\left(\frac{1}{2}(\sigma^{-}_{h}+\sigma^{+}_{h})+\frac{Y}{2}(\hat{\mathbf{n}}^{-}_{r}v^{-}_{h} + \hat{\mathbf{n}}^{+}_{r}v^{+}_{h})\right)\boldsymbol{\hat{\psi}}(1)\right. \\	&-\left.\left(\frac{1}{2}(\sigma^{-}_{h}+\sigma^{+}_{h})+\frac{Y}{2}(\hat{\mathbf{n}}^{-}_{l}v^{-}_{h} + \hat{\mathbf{n}}^{+}_{l}v^{+}_{h})\right)\boldsymbol{\hat{\psi}}(-1)\right]\end{align*}where the subscripts $r$ or $l$ for the normal vector $\hat{\mathbf{n}}$ indicates that we are taking the normal vector at the right or left side of the $k$ element, respectively. Substituting the expressions for the approximations $\sigma_{h}$ and $v_{h}$ into the equation above and remembering that $\hat{\mathbf{n}}^{-} = -\hat{\mathbf{n}}^{+}$ yields\begin{align}\frac{1}{\rho}\left[(\sigma_{h})^{*}\boldsymbol{\hat{\psi}}(1)-(\sigma_{h})^{*}\boldsymbol{\hat{\psi}}(-1)\right] =& \: \frac{1}{\rho}\left[\frac{1}{2}\left((\hat{\boldsymbol{\sigma}}^{k}_{h})^{T}\boldsymbol{\hat{\psi}}(1)+(\hat{\boldsymbol{\sigma}}^{k+1}_{h})^{T}\boldsymbol{\hat{\psi}}(-1)\right)\boldsymbol{\hat{\psi}}(1)\right. \\&+\frac{Y}{2}\hat{\mathbf{n}}^{-}_{r}\left((\hat{\boldsymbol{v}}^{k}_{h})^{T}\boldsymbol{\hat{\psi}}(1)-(\hat{\boldsymbol{v}}^{k+1}_{h})^{T}\boldsymbol{\hat{\psi}}(-1)\right)\boldsymbol{\hat{\psi}}(1) \nonumber \\&-\frac{1}{2}\left((\hat{\boldsymbol{\sigma}}^{k}_{h})^{T}\boldsymbol{\hat{\psi}}(-1)+(\hat{\boldsymbol{\sigma}}^{k-1}_{h})^{T}\boldsymbol{\hat{\psi}}(1)\right)\boldsymbol{\hat{\psi}}(-1) \nonumber \\&-\left.\frac{Y}{2}\hat{\mathbf{n}}^{-}_{l}\left((\hat{\boldsymbol{v}}^{k}_{h})^{T}\boldsymbol{\hat{\psi}}(-1)-(\hat{\boldsymbol{v}}^{k-1}_{h})^{T}\boldsymbol{\hat{\psi}}(1)\right)\boldsymbol{\hat{\psi}}(-1)\right] \\\\=& \: \frac{1}{\rho}\left[\frac{1}{2}\left([\boldsymbol{\hat{\psi}}(1)\otimes\boldsymbol{\hat{\psi}}(1)]\:\hat{\boldsymbol{\sigma}}^{k}_{h}+[\boldsymbol{\hat{\psi}}(1)\otimes\boldsymbol{\hat{\psi}}(-1)]\:\hat{\boldsymbol{\sigma}}^{k+1}_{h}\right)\right. \\&+\frac{Y}{2}\hat{\mathbf{n}}^{-}_{r}\left([\boldsymbol{\hat{\psi}}(1)\otimes\boldsymbol{\hat{\psi}}(1)]\:\hat{\boldsymbol{v}}^{k}_{h}-[\boldsymbol{\hat{\psi}}(1)\otimes\boldsymbol{\hat{\psi}}(-1)]\:\hat{\boldsymbol{v}}^{k+1}_{h}\right) \\&-\frac{1}{2}\left([\boldsymbol{\hat{\psi}}(-1)\otimes\boldsymbol{\hat{\psi}}(-1)]\:\hat{\boldsymbol{\sigma}}^{k}_{h}+[\boldsymbol{\hat{\psi}}(-1)\otimes\boldsymbol{\hat{\psi}}(1)]\:\hat{\boldsymbol{\sigma}}^{k-1}_{h}\right) \nonumber \\&-\left.\frac{Y}{2}\hat{\mathbf{n}}^{-}_{l}\left([\boldsymbol{\hat{\psi}}(-1)\otimes\boldsymbol{\hat{\psi}}(-1)]\:\hat{\boldsymbol{v}}^{k}_{h}-[\boldsymbol{\hat{\psi}}(-1)\otimes\boldsymbol{\hat{\psi}}(1)]\:\hat{\boldsymbol{v}}^{k-1}_{h}\right)\right]\end{align}Tensor multiplication of vectors of elemental test functions gives rise to matrices which we will call flux matrices. They are defined as follows\begin{align*}\mathcal{F}^{k}_{r} &= [\boldsymbol{\hat{\psi}}(1)\otimes\boldsymbol{\hat{\psi}}(1)], \\\mathcal{F}^{k+1}_{l} &= [\boldsymbol{\hat{\psi}}(1)\otimes\boldsymbol{\hat{\psi}}(-1)], \\\mathcal{F}^{k}_{l} &= [\boldsymbol{\hat{\psi}}(-1)\otimes\boldsymbol{\hat{\psi}}(-1)], \\\mathcal{F}^{k-1}_{r} &= [\boldsymbol{\hat{\psi}}(-1)\otimes\boldsymbol{\hat{\psi}}(1)].\end{align*} We now use the flux matrices definitions and rearrange terms to write the final matricial form for the velocity and stress equations\begin{align}\frac{1}{\rho}\left[(\sigma_{h})^{*}\boldsymbol{\hat{\psi}}(1)-(\sigma_{h})^{*}\boldsymbol{\hat{\psi}}(-1)\right] = \: \frac{1}{\rho}\left[\frac{1}{2}\mathcal{F}^{k+1}_{l}\left(\hat{\boldsymbol{\sigma}}^{k+1}_{h}-(\hat{\mathbf{n}}^{-}_{r}\cdot Y)\hat{\boldsymbol{v}}^{k+1}_{h}\right) \right. &+ \frac{1}{2}\mathcal{F}^{k}_{r}\left(\hat{\boldsymbol{\sigma}}^{k}_{h}+(\hat{\mathbf{n}}^{-}_{r}\cdot Y)\hat{\boldsymbol{v}}^{k}_{h}\right) \nonumber \\-\frac{1}{2}\mathcal{F}^{k-1}_{r}\left(\hat{\boldsymbol{\sigma}}^{k-1}_{h}-(\hat{\mathbf{n}}^{-}_{l}\cdot Y)\hat{\boldsymbol{v}}^{k-1}_{h}\right) &- \left.\frac{1}{2}\mathcal{F}^{k}_{l}\left(\hat{\boldsymbol{\sigma}}^{k}_{h}+(\hat{\mathbf{n}}^{-}_{l}\cdot Y)\hat{\boldsymbol{v}}^{k}_{h}\right)\right].\end{align} Repeating exactly the same steps for the numerical flux, $(v_{h})^{*}$, yields\begin{align}\mu(x^{k}_{r})(v_{h})^{*}\boldsymbol{\hat{\psi}}(1)-\mu(x^{k}_{l})(v_{h})^{*}\boldsymbol{\hat{\psi}}(-1) = \: &\mu(x^{k}_{r})\left[\frac{1}{2}\mathcal{F}^{k+1}_{l}\left(\hat{\boldsymbol{v}}^{k+1}_{h}-(\hat{\mathbf{n}}^{-}_{r}\cdot Z)\hat{\boldsymbol{\sigma}}^{k+1}_{h}\right) + \frac{1}{2}\mathcal{F}^{k}_{r}\left(\hat{\boldsymbol{v}}^{k}_{h}+(\hat{\mathbf{n}}^{-}_{r}\cdot Z)\hat{\boldsymbol{\sigma}}^{k}_{h}\right)\right] \\-&\mu(x^{k}_{l}) \left[\frac{1}{2}\mathcal{F}^{k-1}_{r}\left(\hat{\boldsymbol{v}}^{k-1}_{h}-(\hat{\mathbf{n}}^{-}_{l}\cdot Z)\hat{\boldsymbol{\sigma}}^{k-1}_{h}\right)+\frac{1}{2}\mathcal{F}^{k}_{l}\left(\hat{\boldsymbol{v}}^{k}_{h}+(\hat{\mathbf{n}}^{-}_{l}\cdot Z)\hat{\boldsymbol{\sigma}}^{k}_{h}\right)\right].\end{align} Noting that the exterior normal vector, $\hat{\mathbf{n}}^{-}$, assumes only two values, $\hat{\mathbf{n}}^{-}_{r} = 1$ and $\hat{\mathbf{n}}^{-}_{l} = -1$, at the right and the left boundaries, respectively, of any $k$ element and the variables, $Z$, and, $Y$, can be computed from the expressions for the fluxes, the matrix form of the proposed DG scheme for the velocity-stress formulation of the elastic wave equation can finally be written as\begin{align}\frac{h^{k}}{2}\mathcal{M}\frac{d}{dt}\hat{\boldsymbol{v}}^{k}_{h} + \frac{1}{\rho}\mathcal{S}_{v}\hat{\boldsymbol{\sigma}}^{k}_{h} - \frac{1}{\rho}\frac{h^{k}}{2}\mathbf{f} =& \: \frac{1}{\rho}\left[\frac{1}{2}\mathcal{F}^{k+1}_{l}\left(\hat{\boldsymbol{\sigma}}^{k+1}_{h}+(\rho \: \mu(x^{k}_{r}))^{1/2}\hat{\boldsymbol{v}}^{k+1}_{h}\right) + \frac{1}{2}\mathcal{F}^{k}_{r}\left(\hat{\boldsymbol{\sigma}}^{k}_{h}-(\rho \: \mu(x^{k}_{r}))^{1/2}\hat{\boldsymbol{v}}^{k}_{h}\right)\right. \\&- \ \frac{1}{2}\mathcal{F}^{k-1}_{r}\left(\hat{\boldsymbol{\sigma}}^{k-1}_{h}-(\rho \: \mu(x^{k}_{l}))^{1/2}\hat{\boldsymbol{v}}^{k-1}_{h}\right) - \left.\frac{1}{2}\mathcal{F}^{k}_{l}\left(\hat{\boldsymbol{\sigma}}^{k}_{h}+(\rho \: \mu(x^{k}_{l}))^{1/2}\hat{\boldsymbol{v}}^{k}_{h}\right)\right], \\\frac{h^{k}}{2}\mathcal{M}\frac{d}{dt}\hat{\boldsymbol{\sigma}}^{k}_{h} + \mathcal{S}_{s}\hat{\boldsymbol{v}}^{k}_{h} + \frac{h^{k}}{2}\mathcal{S}_{m}\hat{\boldsymbol{v}}^{k}_{h} =& \: \mu(x^{k}_{r})\left[\frac{1}{2}\mathcal{F}^{k+1}_{l}\left(\hat{\boldsymbol{v}}^{k+1}_{h}+(\rho \: \mu(x^{k}_{r}))^{-1/2}\hat{\boldsymbol{\sigma}}^{k+1}_{h}\right) + \frac{1}{2}\mathcal{F}^{k}_{r}\left(\hat{\boldsymbol{v}}^{k}_{h}-(\rho \: \mu(x^{k}_{r}))^{-1/2}\hat{\boldsymbol{\sigma}}^{k}_{h}\right)\right] \\-& \: \mu(x^{k}_{l}) \left[\frac{1}{2}\mathcal{F}^{k-1}_{r}\left(\hat{\boldsymbol{v}}^{k-1}_{h}-(\rho \: \mu(x^{k}_{l}))^{-1/2}\hat{\boldsymbol{\sigma}}^{k-1}_{h}\right)+\frac{1}{2}\mathcal{F}^{k}_{l}\left(\hat{\boldsymbol{v}}^{k}_{h}+(\rho \: \mu(x^{k}_{l}))^{-1/2}\hat{\boldsymbol{\sigma}}^{k}_{h}\right)\right].\end{align} It is important to observe some important properties of the DG formulation:- The scheme is completely local (there is no global assembly).- Boundary conditions and interface continuity are enforced weakly via numerical fluxes.- The solutions are piecewise smooth, polynomial in general, but discontinuous between elements.- The scheme is well suited for order and element sizes adaption. - Parallelism is very natural because of the strictly local formulation. The source term $f$ in the original displacement form of the elastic wave equation is defined as a point source$$f(x,t) = \theta(t)\delta(x-x_{\text{source}})$$where $x_{\text{source}}$ is chosen to be any quadrature point inside any $k$ element within the domain $\Omega$. We do not allow $x_{\text{source}}$ to be chosen at the element interfaces, which is a possible situation when using GLL quadrature points, in order to avoid the definition of a double valued Dirac delta function. To consider this possibility we should include some upwind feature in the definition of the Dirac delta function but we postponed this tricky theoretical approach to the future. Having defined the source term $f^{k}$ on the $k$-th element, we can write the final form of the vector source term $\hat{\mathbf{f}}^{k}$ as\begin{align}\hat{\mathbf{f}}^{k} &= \int_{\hat{\Omega}}\theta(t)\delta(x-x_{\text{source}})\psi_{i}(x) \: dx, \ i = 1,\ldots,N+1, \\&= \theta(t)\psi_{i}(x_{\text{source}}), \ i = 1,\ldots,N+1.\end{align}We observe that the vector form of the source term as defined above will have non zero entries on one single element, $k$. It is also worth noting that when using the Lagrange basis functions, the resulting source vector on the $k$-th element, will have just a single non zero entry $\mathbf{f}^{k}_{i}$ corresponding to the coefficient associated to the $i^{th}$-Lagrange polynomial $\ell^{k}_{i}(x_{\text{source}})$ defined at the collocation point $x_{\text{source}}$ which must be one of the GLL nodes. In this case the source term simplifies to $\mathbf{f}^{k}_{i} = \theta(t)$ because of the cardinality property of the Lagrange polynomial.  Time MarchingThe DG-FEM spatial discretization of the elastic wave equation produces a system of first order ODEs in terms of the velocity and stress fields of the problem.  Before discretizing the semidiscrete problem given, we define the right hand side operator $\mathcal{L}_{h}$ for the variables $\boldsymbol{v}_{h}$ and $\boldsymbol{\sigma}_{h}$ as\begin{align}\mathcal{L}_{h}(\boldsymbol{v}_{h},t) =& \mathcal{M}^{-1} \! \left\{\frac{1}{\rho}\mathbf{f} - \frac{2}{h^{k}}\frac{1}{\rho}\mathcal{S}_{v}\hat{\boldsymbol{\sigma}}^{k}_{h} + \frac{2}{h^{k}}\frac{1}{\rho}\left[\frac{1}{2}\mathcal{F}^{k+1}_{l}\left(\hat{\boldsymbol{\sigma}}^{k+1}_{h}+(\rho\mu(x^{k}_{r}))^{1/2}\hat{\boldsymbol{v}}^{k+1}_{h}\right) + \frac{1}{2}\mathcal{F}^{k}_{r}\left(\hat{\boldsymbol{\sigma}}^{k}_{h}-(\rho\mu(x^{k}_{r}))^{1/2}\hat{\boldsymbol{v}}^{k}_{h}\right)\right.\right. \\-& \left.\ \frac{1}{2}\mathcal{F}^{k-1}_{r}\left(\hat{\boldsymbol{\sigma}}^{k-1}_{h}-(\rho\mu(x^{k}_{l}))^{1/2}\hat{\boldsymbol{v}}^{k-1}_{h}\right) - \left.\frac{1}{2}\mathcal{F}^{k}_{l}\left(\hat{\boldsymbol{\sigma}}^{k}_{h}+(\rho\mu(x^{k}_{l}))^{1/2}\hat{\boldsymbol{v}}^{k}_{h}\right)\right]\right\}, \\\mathcal{L}_{h}(\boldsymbol{\sigma}_{h},t) \! =& \mathcal{M}^{-1}\!\left\{ \! -\frac{2}{h^{k}}\mathcal{S}_{s}\hat{\boldsymbol{v}}^{k}_{h} - \mathcal{S}_{m}\hat{\boldsymbol{v}}^{k}_{h} + \frac{2}{h^{k}}\mu(x^{k}_{r}) \! \left[\frac{1}{2}\mathcal{F}^{k+1}_{l}\left(\hat{\boldsymbol{v}}^{k+1}_{h}+(\rho\mu(x^{k}_{r}))^{-1/2}\hat{\boldsymbol{\sigma}}^{k+1}_{h}\right) + \frac{1}{2}\mathcal{F}^{k}_{r}\left(\hat{\boldsymbol{v}}^{k}_{h}-(\rho\mu(x^{k}_{r}))^{-1/2}\hat{\boldsymbol{\sigma}}^{k}_{h}\right)\right]\right. \\-& \left.\: \frac{2}{h^{k}}\mu(x^{k}_{l}) \left[\frac{1}{2}\mathcal{F}^{k-1}_{r}\left(\hat{\boldsymbol{v}}^{k-1}_{h}-(\rho\mu(x^{k}_{l}))^{-1/2}\hat{\boldsymbol{\sigma}}^{k-1}_{h}\right)+\frac{1}{2}\mathcal{F}^{k}_{l}\left(\hat{\boldsymbol{v}}^{k}_{h}+(\rho\mu(x^{k}_{l}))^{-1/2}\hat{\boldsymbol{\sigma}}^{k}_{h}\right)\right]\right\}.\end{align} As we are also interested in the displacement field of the problem, we can add this field to the system and write it in a new vector variable, $\boldsymbol{\phi}$,$$\boldsymbol{\phi}_{h} = \left\{\begin{array}{c}\boldsymbol{u}_{h} \\\boldsymbol{v}_{h} \\\boldsymbol{\sigma}_{h}\end{array}\right\}.$$Now, the semidiscrete problem can be rewritten in the form\begin{align}\frac{d\boldsymbol{\phi}_{h}}{dt} = \mathcal{L}_{h}(\boldsymbol{\phi}_{h},t),\end{align}or, in vector notation$$\frac{d\boldsymbol{\phi}_{h}}{dt} = \frac{d}{dt}\left\{\begin{array}{c}\boldsymbol{u}_{h} \\\boldsymbol{v}_{h} \\\boldsymbol{\sigma}_{h}\end{array}\right\} = \left\{\begin{array}{c}\boldsymbol{v}_{h} \\\mathcal{L}_{h}(\boldsymbol{v}_{h},t) \\\mathcal{L}_{h}(\boldsymbol{\sigma}_{h},t)\end{array}\right\}.$$In the above expression, we are assuming that the right hand side operator $\mathcal{L}_{h}$ evaluated at the displacement variable $\mathbf{u}_{h}$ simply returns the variable $\boldsymbol{v}_{h}$, i.e., $\mathcal{L}_{h}(\boldsymbol{u}_{h},t) = \boldsymbol{v}_{h}$. To march this system of equations in time, we can use the classical and reliable forth order four stage explicit Runge-Kutta method (ERK)\begin{align}&\boldsymbol{k}^{(1)} = \mathcal{L}_{h}(\boldsymbol{\phi}^{n}_{h},t^{n}), \\&\boldsymbol{k}^{(2)} = \mathcal{L}_{h}\left(\boldsymbol{\phi}^{n}_{h}+\frac{1}{2}\mathit{\Delta} t\boldsymbol{k}^{(1)},t^{n}+\frac{1}{2}\mathit{\Delta} t\right), \\&\boldsymbol{k}^{(3)} = \mathcal{L}_{h}\left(\boldsymbol{\phi}^{n}_{h}+\frac{1}{2}\mathit{\Delta} t\boldsymbol{k}^{(2)},t^{n}+\frac{1}{2}\mathit{\Delta} t\right), \\&\boldsymbol{k}^{(4)} = \mathcal{L}_{h}\left(\boldsymbol{\phi}^{n}_{h}+\mathit{\Delta} t\boldsymbol{k}^{(3)},t^{n}+\mathit{\Delta} t\right), \\&\boldsymbol{\phi}^{n+1}_{h} = \boldsymbol{\phi}^{n}_{h}+\frac{1}{6}\mathit{\Delta} t\left(\boldsymbol{k}^{(1)}+2\boldsymbol{k}^{(2)}+2\boldsymbol{k}^{(3)}+\boldsymbol{k}^{(4)}\right),\end{align}to advance from $\boldsymbol{\phi}^{n}_{h}$ to $\boldsymbol{\phi}^{n+1}_{h}$ which are separated by the timestep, $\mathit{\Delta} t$. It is important to note that the ERK method is conditionally stable and requires the timestep size $\mathit{\Delta} t$ be estimated using the $CFL$ (Courant-Friedrichs-Lewy) condition. We define the CFL as$$CFL := \max_{k}\frac{c^{k}\mathit{\Delta} t}{\mathit{\Delta} x^{k}},$$where $c^{k}$ is the local propagation wave speed on element $k$ and $\mathit{\Delta} x^{k} = h^{k}$ is the local spatial grid resolution. The specific value of the $CFL$ (constant) number depends on the method used for the space and time discretizations but, in general, its order of magnitude is about 1. For the present case, we consider $CFL = 0.4$. As a last remark, we observe that, for high order methods such as the DG scheme described in this text, we have to consider the minimum grid space as the minimum distance between the mapped quadrature points (GL or GLL) on the mesh elements.  Python implementation  Simulation parametersThe first step to use the Python DG framework is the definition of the simulation parameters. There is a Python class specifically designed for this, called _SimulationData_, implemented on `simulation_data.py`. Also on the same file are two service classes, _QuadratureNodes_ and _BasisType_ used to represent minemonically the types of quadrature points and basis available in the framework. QuadratureNodesThere are two variables defined in _QuadratureNodes_, both related to Gaussian integration:- `QuadratureNodes.GL`, corresponding to the Legendre integration points, and- `QuadratureNodes.GLL`, associated with the Legendre-Lobatto points. BasisTypeHere we have three variables:- `BasisType.NODAL`, used to indicate we are using a nodal basis;- `BasisType.MODAL`, the same as above, for modal basis, and > Remember to use `%load simulation_data.py` to bring the file to the notebook! Remember also to __run__ the cell!<jupyter_code># %load simulation_data.py
#
# Common simulation data
# ======================
# It presents a struct to store simulation data common to all problemas, like
# number of elements, polynomial order and so on.
#
# by <NAME>. (Matlab and Python versions)
#    <NAME> (Python version)
#
from math import ceil, floor

#
# Used to specify the chosen basis.
#
class BasisType(object):
    """The available basis types."""
    NODAL = 0
    MODAL = 1
    MODAL_ORTHOGONAL = 2
    MODAL_ORTHONORMAL = 3

#
# Used to specify the chosen type of points distribution for Gaussian quadrature
#
class QuadratureNodes(object):
    """The available quadrature nodes distribution."""
    GL  = 0     # Gauss-Legendre
    GLL = 1     # Gauss-Legendre-Lobatto

class SimulationData:
    """Creates a simulation data structure."""
    def __init__(self,
                 n_elements      = 1,
                 poly_order      = 0,
                 spatial_domain  = (0.0, 1.0),
                 temporal_domain = (0.0, 1.0),
                 basis_type      = BasisType.NODAL,
                 node_dist       = QuadratureNodes.GL,
                 non_linear      = False):
        """
        Parameters
        ----------
        n_elements
            Number of elements
        poly_order
            Polynomial order of approximation
        spatial_domain
            Spatial domain
        temporal_domain
            Temporal domain
        basis_type
            The basis used in the approximation
        node_dist
            The quadrature points distribution
        """
        self.n_elements      = n_elements
        self.poly_order      = poly_order
        self.spatial_domain  = spatial_domain
        self.temporal_domain = temporal_domain
        self.basis_type      = basis_type
        self.node_dist       = node_dist
        self.non_linear      = non_linear

    def mass_order(self):
        """Number of integration points needed for mass matrix assembly."""
        return self.poly_order*2

    def nonlinear_order(self):
        """Number of integration points needed for nonlinear terms."""
        return self.poly_order*3

    def spatial_domain_length(self):
        """Length of the spatial domain."""
        return self.spatial_domain[1] - self.spatial_domain[0]

    def n_local_dof(self):
        """Local degrees of freedom.
        
        In 1D we have the polynomial order plus 1, to account for the constant
        term."""
        return self.poly_order+1

    def n_global_dof(self):
        """Global degrees of freedom."""
        return self.n_local_dof()*self.n_elements

    def nip(self):
        order = self.mass_order() if self.non_linear == False else self.nonlinear_order()

        if self.node_dist == QuadratureNodes.GL:
            return int(ceil(0.5*(order+2)))
        elif self.node_dist == QuadratureNodes.GLL:
            # floor(nip) will keep the nodal spectral approach Np = N+1
            return int(floor((order+3)/2))
        else:
            raise AssertionError( "Wrong quadrature node distribution!\n" )

    def t_min(self):
        return self.temporal_domain[0]

    def t_max(self):
        return self.temporal_domain[1]

    # Boilerplate code
    def __str__(self):
        return ("SimulationData:\n"
                "  Number of elements: %s\n"
                "  Polynomial order  : %s\n"
                "  Spatial domain    : %s\n"
                "  Temporal domain   : %s") % (self.n_elements,
                                               self.poly_order,
                                               self.spatial_domain,
                                               self.temporal_domain)

#-- simulation_data.py ---------------------------------------------------------

# - Elements: 10
# - Polynomial order: 8
# - Spatial domain: [0, 2]
# - Temporal domain: [0, 10]
# - Gaussian integration with Legendre-Lobatto points
# - Modal basis (orthonormal)
sim_data = SimulationData(n_elements=20, poly_order=8, spatial_domain=[0.0, 9000.0],
                          temporal_domain=[0, 10], node_dist=QuadratureNodes.GL,
                          basis_type=BasisType.MODAL)<jupyter_output><empty_output><jupyter_text>The finite element meshGiven that we are on a 1D setting, a very simple mesh generator is easy to implement. Of course, on 2D or 3D we would resort to a third-party application. Note that the mesh is created based on the simulation data previously defined. The mesh is also responsible for evaluating the Jacobian.The mesh class is _Mesh1D_, defined in `mesh_gen_1d.py`, consists -- at this point -- of two data members:- `coord`, $x$ coordinate of every node on the mesh, and- `conn`, connectivity matrix, associating elements and nodes<jupyter_code># %load mesh_gen_1d.py
# 
# 1D mesh generation
# ==================
# Generate simple equidistant grid with a given number of elements elements
#
# by <NAME> Jr. (Matlab and Python versions)
#    <NAME> (Python version)
#
from simulation_data import SimulationData
import numpy as np

class Mesh1D:
    """Simple one-dimensional mesh.
    
    Arguments:
        sim_data - SimulationData object
    """
    def __init__(self, sim_data):
        # Generate node coordinates
        # - Note that we have n_elements+1 points on the domain.
        # - The name 'coord' is classic in FEM codes.
        self.coord = numpy.linspace( sim_data.spatial_domain[0],
                                  sim_data.spatial_domain[1],
                                  sim_data.n_elements+1 )

        # Element to node connectivity
        # - The name 'conn' is classic in FEM codes.
        self.conn = numpy.array([numpy.arange(0, sim_data.n_elements, 1),
                              numpy.arange(1, sim_data.n_elements+1, 1)]).T

    def jacobian(self, sim_data, quad_data):
        nip = sim_data.nip()

        # The Jacobian
        self.J = 0.5*(self.coord[1:]-self.coord[0:-1])

        # Elements' mid-points
        B = 0.5*(self.coord[1:]+self.coord[0:-1])
        # We repeat and reshape the array. Why?
        B = numpy.repeat(B, nip).reshape(sim_data.n_elements, nip).T

        # The transformed points. The construct 'numpy.newaxis' is needed here
        # because we want a de facto matricial product (rank-1 update).
        #
        # x = xi^T * J
        self.x = quad_data.xi[numpy.newaxis].T.dot(self.J[numpy.newaxis]) + B

#-- mesh_gen_1d.py -------------------------------------------------------------

# Generate simple mesh
mesh = Mesh1D(sim_data)
# Lets plot the mesh
pyplot.figure(1)
pyplot.clf()
pyplot.xlim(sim_data.spatial_domain[0]-0.1, sim_data.spatial_domain[1]+0.1)
pyplot.plot(mesh.coord, numpy.zeros(mesh.coord.shape), "ro")
pyplot.title("FEM mesh")
pyplot.show()<jupyter_output><empty_output><jupyter_text>Gaussian quadrature points and weightsHaving defined the mesh, now we proceed with the creation of the Gaussian quadrature points and integration weights. For this purpose we can resort to the _JacobiGaussQuad_ class, defined in `jacobi_gauss_quad.py`. Objects of this class have two data members:- `xi`, the integration points, and- `w`, the integration weights<jupyter_code># %load jacobi_gauss_quad.py
#
# Points and weights for Gaussian quadrature
#
# by <NAME> Jr. (Matlab and Python versions)
#    <NAME> (Python version)
#
import numpy
from math import gamma
from jacobi_p import jacobi_p
from djacobi_p import djacobi_p
from jacobi_roots import jacobi_roots
from simulation_data import QuadratureNodes

class JacobiGaussQuad:
    """Points and weights for Gaussian quadrature based on Jacobi polynomials"""
    def __init__(self, sim_data):
        # Sets parameters to obtain quadrature points from Legendre polynomials
        # Legendre == Jacobi(0,0)
        alpha = 0.0
        beta  = 0.0
        nip = sim_data.nip()

        # Case 1: Gauss-Legendre quadrature
        if sim_data.node_dist == QuadratureNodes.GL:
            self.xi, self.w = self._gl(nip, alpha, beta)
        # Case 2: Gauss-Legendre-Lobato quadrature
        elif sim_data.node_dist == QuadratureNodes.GLL:
            self.xi, self.w = self._gll(nip, alpha, beta)
        else:
            raise AssertionError("Unknown quadrature type!")

    # Case 1: Gauss-Legendre quadrature
    def _gl(self, nip, alpha, beta):
        xi = jacobi_roots(nip, alpha, beta)

        C1 = (2.0**(alpha+beta+1.0))*gamma(alpha+nip+1.0)*gamma(beta+nip+1.0)
        C2 = gamma(nip+1.0)*gamma(alpha+beta+nip+1.0)*(1.0-xi**2 )

        DPm = djacobi_p(xi, nip, alpha, beta)

        w = C1*DPm**(-2)/C2

        return xi, w

    # Case 2: Gauss-Legendre-Lobato quadrature
    def _gll(self, nip, alpha, beta):
        r = jacobi_roots(nip-2, alpha+1.0, beta+1.0)
        xi = numpy.empty(r.shape[0]+2)
        xi[0] = -1.0
        xi[1:-1] = r
        xi[-1] = 1.0

        C1 = (2.0**(alpha+beta+1.0))*gamma(alpha+nip)*gamma(beta+nip)
        C2 = (nip-1)*gamma(nip)*gamma(alpha+beta+nip+1.0)

        Pm = jacobi_p(xi, nip-1, alpha, beta)

        w = C1*Pm**(-2)/C2

        w[ 0] = w[ 0]*(beta+1.0)
        w[-1] = w[-1]*(alpha+1.0)

        return xi, w

    def n(self):
        """Number of integration points / weights."""
        return self.xi.shape[0]

#-- jacobi_gauss_quad.py -------------------------------------------------------

# Creating the quadrature points and weights
quad_data = JacobiGaussQuad(sim_data)
# Lets check the distribution of points and respective weights
pyplot.figure(2)
pyplot.clf()
pyplot.xlim(-1.1, 1.1)
pyplot.plot(quad_data.xi, quad_data.w, "bo")
pyplot.plot(quad_data.xi, numpy.zeros(quad_data.xi.shape), "ro")
pyplot.show()<jupyter_output><empty_output><jupyter_text>Basis functionsThe _BasisFunctions_ class implements the three possible choices of basis functions (nodal, modal e modal-orthogonal) in a convenient matricial form accessible through the data member `psi`. Basis functions derivatives are also evaluated (data member `Dpsi`). Beyond the _BasisFunctions_ class itself, the file `basis_functions.py` brings the individual implementations of Lagrange, Legendre and Bubble functions.<jupyter_code># %load basis_functions.py
#
# Basis functions and respective derivatives
#
# by <NAME>. (Matlab and Python versions)
#    <NAME> (Python version)
#
import numpy
from sys import float_info
from math import fabs, sqrt
from jacobi_p import jacobi_p
from djacobi_p import djacobi_p, d_pq, d2_pq
from jacobi_roots import jacobi_roots
from simulation_data import QuadratureNodes, BasisType

def gll_nodes(n_local_dof):
    """Gauss-Legendre-Lobatto auxiliar points.
    
    Parameters
    ----------
    n_local_dof
        Per-element number of degrees of freedom
    """
    if n_local_dof < 2:
        raise AssertionError("Local degress of freedom must be >= 2!")

    # Setting GLL nodes
    alpha = 0.0
    beta  = 0.0

    xj = numpy.empty(n_local_dof)

    # GLL points are the roots of a Jacobi polynomial with alpha = beta = 1,
    # plus the interval extrema, -1 and 1.
    xj[0] = -1.0
    xj[1:-1] = jacobi_roots(n_local_dof-2, alpha+1.0, beta+1.0)
    xj[-1] = 1.0

    return xj

def lagrange(n_local_dof, xi, xj):
    """Lagrange basis function.

    Parameters
    ----------
    n_local_dof
        Number of local degrees of freedom
    xi
        Distribution of integration points on the element defining the
        Lagrangian basis
    xj
        Points where we evaluate the basis
    """
    # Number of Lagrange basis functions
    Q   = n_local_dof
    QQ1 = Q*(Q-1)

    # Number of evaluation points
    npts = xi.shape[0]

    # Allocates matrix to store Lagrange basis evaluated at npts points
    Hj = numpy.ones([npts, Q])
    Hj[:, 0] = ((-1.0)**(Q-1)/QQ1)*(xi-1.0)*(0.5*Q*jacobi_p(xi, Q-2, 1.0, 1.0))
    Hj[:,-1] = (1.0/QQ1)*(xi+1.0)*(0.5*Q*jacobi_p(xi, Q-2, 1.0, 1.0))
    
    Lp  =  jacobi_p(xj, Q-1, 0.0, 0.0)
    DLp = djacobi_p(xi, Q-1, 0.0, 0.0)

    for j in range(1, Q-1):   # from the second to the second-last column
        for i in range(0, npts):
            if fabs(xi[i]-xj[j]) > float_info.epsilon:
                Hj[i, j] = (xi[i]**2-1.0)*DLp[i]/(QQ1*Lp[j]*(xi[i]-xj[j]))

    return Hj

def bubble_basis(n_local_dof, xi):
    """Modal bubble functions.
    n_local_dof
        Number of local degrees of freedom
    xi
        Distribution of integration points on the element
    """
    # Number of bubble basis functions
    Q = n_local_dof

    # Number of evaluation points
    npts = xi.shape[0]

    # Allocates matrix to store Lagrange basis evaluated at npts points
    Bj = numpy.empty([npts, Q])
    
    # Boundary nodal basis functions (necessary to enforce C^(0) continuity)
    Bj[:, 0] = 0.5*(1.0-xi)
    Bj[:,-1] = 0.5*(1.0+xi)

    # Bubble modes
    for j in range(1, Q-1):
        Lj   = jacobi_p(xi, j+1, 0.0, 0.0)
        Ljm2 = jacobi_p(xi, j-1, 0.0, 0.0)
        Bj[:, j] = (Lj-Ljm2)/sqrt(2.0*(2.0*j+1.0))

    return Bj
    
def legendre_basis(n_local_dof, xi):
    """Legendre basis functions.
    n_local_dof
        Number of local degrees of freedom
    xi
        Distribution of integration points on the element
    """
    # Number of basis functions
    Q = n_local_dof

    # Number of evaluation points
    npts = xi.shape[0]

    # Allocates matrix to store Lagrange basis evaluated at npts points
    Lj = numpy.empty([npts, Q])
    
    # Legendre basis
    for j in range(0, Q):
        Lj[:, j] = jacobi_p(xi, j, 0.0, 0.0)

    return Lj
    
def orthonormal_legendre_basis(n_local_dof, xi):
    """Legendre basis functions.
    n_local_dof
        Number of local degrees of freedom
    xi
        Distribution of integration points on the element
    """
    # Number of basis functions
    Q = n_local_dof

    # Number of evaluation points
    npts = xi.shape[0]

    # Allocates matrix to store Lagrange basis evaluated at npts points
    Lj = numpy.empty([npts, Q])
    
    # Legendre basis with the orthonormalization factor
    for j in range(0, Q):
        Lj[:, j] = jacobi_p(xi, j, 0.0, 0.0)/numpy.sqrt(2.0/(2.0*j+1.0))

    return Lj
    
class BasisFunctions:
    """Basis functions and derivatives on quadrature points."""
    def __init__(self, sim_data, gauss_quad):
        """
        Parameters
        ----------
        sim_data
            SimulationData object
        """
        n_local_dof = sim_data.n_local_dof()

        # Constructing the basis functions PHI
        if sim_data.basis_type == BasisType.NODAL:
            xj = gll_nodes(n_local_dof)
            self.psi = lagrange(n_local_dof, gauss_quad.xi, xj)
        elif sim_data.basis_type == BasisType.MODAL:
            self.psi = bubble_basis(n_local_dof, gauss_quad.xi)
        elif sim_data.basis_type == BasisType.MODAL_ORTHOGONAL:
            self.psi = legendre_basis(n_local_dof, gauss_quad.xi)
        elif sim_data.basis_type == BasisType.MODAL_ORTHONORMAL:
            self.psi = orthonormal_legendre_basis(n_local_dof, gauss_quad.xi)
        else:
            raise AssertionError( "Wrong basis type!\n" )

        # Now we construct the derivative of the basis functions, DPHI/DX
        nip = sim_data.nip()

        self.D = numpy.zeros([nip, nip])

        if nip > 1:
            dp  =  d_pq(gauss_quad.xi, nip, 0.0, 0.0, sim_data.node_dist)
            dp2 = d2_pq(gauss_quad.xi, nip, 0.0, 0.0, sim_data.node_dist)
        
            for i in range(0, nip):
                for j in range(0, nip):
                    # Never do "ifs" inside loops! Just for pedagogic reasons.
                    if i == j:
                        self.D[i, j] = dp2[i]/(2.0*dp[i])
                    else:
                        dx = gauss_quad.xi[i] - gauss_quad.xi[j]
                        self.D[i, j] = dp[i]/(dp[j]*dx)

        self.dpsi = self.D.dot(self.psi)

#-- basis_functions.py ---------------------------------------------------------

# Basis functions and derivatives evaluated at quadrature points
basis = BasisFunctions(sim_data, quad_data)<jupyter_output><empty_output><jupyter_text>JacobianAs we have seen in the theoretical part, there is a Jacobian involved in the numerical integration due to the transformation of all integrals from the real elements to the standard one. The Jacobian is evaluated by the same _Mesh1D_ class used before. After Jacobian evaluation two more data members become available:- `J`, the Jacobian matrix, and- `x`, a mapping of the integration points (remember, they are all in $[-1,1]$) to the real elements<jupyter_code>mesh.jacobian(sim_data, quad_data)<jupyter_output><empty_output><jupyter_text>Material propertiesIn order to be able to run the simulation, we need to evaluate the properties associated with the medium. They are- Density, $\rho(x)$ in $kg/m^3$- Shear modulus, $\mu(x)$ in $Pa$ and its derivative $d\mu/dx$- Wavefield velocity, $c$ in $m/s$The class _MaterialProperties_ is responsible for the evaluation and book keeping of all these properties.<jupyter_code># %load material_properties.py
#
# Sets several material properties for the elastic wave equation
#
# by <NAME>. (Matlab and Python versions)
#    <NAME> (Python version)
#
import numpy
from simulation_data import QuadratureNodes
from math import ceil

class MaterialProperties:
    """Material properties for the elastic wave equation"""
    def __init__(self, sim_data, quad_data, mesh, x_source):
        """MaterialProperties constructor

        Arguments:
            sim_data - SimulationData object
            quad_data - Quadrature points and weights
            mesh - 1D Finite Elements mesh
            x_source - Ricker wavelet point source location
        """
        nip   = sim_data.nip()
        nldof = sim_data.n_local_dof()
        dx    = sim_data.spatial_domain_length()

        self.rho   = numpy.zeros(sim_data.n_elements)         # Material density [kg/m^3]
        self.mu    = numpy.zeros([nip, sim_data.n_elements])  # Material shear modulus [Pa]
        self.c     = numpy.zeros([nip, sim_data.n_elements])  # Material wavefield velocity [m/s]
        self.dmudx = numpy.zeros([nip, sim_data.n_elements])  # Derivative of shear modulus

        # mu evaluated at the boundaries of each element including ghosts
        self.mub = numpy.zeros([nldof, sim_data.n_elements+2]) 
    
        # We could have set that above, but here is clearer
        self.rho.fill(1000.0)

        distance = 0.5*dx-mesh.x

        # Signal changes
        sg = numpy.sign(distance)

        distance = numpy.exp(-(numpy.fabs(distance)-500.0)/100.0)

        self.mu    = 1.2e10-5.96e9/(1.0+distance)
        self.dmudx = 5.96E9*(distance*sg)/(100.0*(1.0+distance)**2)
        self.c     = numpy.sqrt(self.mu/numpy.reshape(numpy.repeat(self.rho, nip), (nip, sim_data.n_elements)))

        distance1 = numpy.fabs(mesh.coord[0:-1]-0.5*dx)
        distance2 = numpy.fabs(mesh.coord[1:]-0.5*dx)

        self.mub[ 0, 1:-1] = 1.2e10-5.96e9/(1.0+numpy.exp(-(distance1-500.0)/100.0))
        self.mub[-1, 1:-1] = 1.2e10-5.96e9/(1.0+numpy.exp(-(distance2-500.0)/100.0))

        self.mub[-1,  0] = self.mub[0, 1]
        self.mub[ 0, -1] = self.mub[-1, -2]

        d = numpy.fabs(mesh.x.T.flatten()-x_source)
        dist = numpy.amin(d)         # The mininum
        globnode = numpy.argmin(d)   # The index
        self.el = int(ceil(globnode/nldof))
        self.locid = globnode-nldof*self.el

        if sim_data.node_dist == QuadratureNodes.GLL:
            if self.locid == 0:
                self.globid = globnode+1
                self.locid = self.locid+1
            elif locid == nldof-1:
                self.globid = globnode-1
                self.locid = self.locid-1
            else:
                self.globid = globnode
        elif sim_data.node_dist == QuadratureNodes.GL:
            self.globid = globnode

#-- material_properties.py -----------------------------------------------------

# Setting material properties
x_source = 3000.0
mat_prop = MaterialProperties(sim_data, quad_data, mesh, x_source)<jupyter_output><empty_output><jupyter_text>The source termThe `generate_source` function creates the source term $f(x)$, as well as the time step and the number of simulation steps.<jupyter_code># %load generate_source.py
#
# Sets several material properties for the elastic wave equation
#
# by <NAME>. (Matlab and Python versions)
#    <NAME> (Python version)
#
import numpy
from math import floor

def generate_source(sim_data, mesh, mat_prop):
    v_max = numpy.amax(mat_prop.c)
    Delta_x_min = numpy.amin(numpy.fabs(mesh.x[0, :]-mesh.x[1, :]))
    CFL = 0.4
    
    dt = CFL*Delta_x_min/v_max
    Nsteps = int(floor(sim_data.t_max()/dt))
    dt = sim_data.t_max()/float(Nsteps)
    
    t_source = 0.1*sim_data.t_max()
    sigma = t_source
    
    # Ricker wavelet over the domain [-5,+5]
    t_eval = numpy.linspace(-5.0, 5.0, floor(t_source/dt))
    Ricker_wavelet = numpy.zeros(Nsteps)
    Ricker_wavelet[0:t_eval.shape[0]] = (2.0/(numpy.sqrt(3.0*sigma)*numpy.pi**0.25))*\
            (1.0-(t_eval/sigma)**2)*numpy.exp(-0.5*(t_eval/sigma)**2)
    
    # Note: the 6e7 is an arbitrary reescaling
    return 6e7*Ricker_wavelet, dt, Nsteps

#-- generate_source.py ---------------------------------------------------------

# Calculates dt and number of time steps
theta, dt, Nsteps = generate_source(sim_data, mesh, mat_prop)<jupyter_output><empty_output><jupyter_text>Making room for the solutionWe will record every step of the solution so we can make a movie showing the temporal evolution.<jupyter_code># We will store the solution at every time step, so we need to make room for it
u = numpy.zeros([sim_data.n_global_dof(), Nsteps+1])

# Projection of the initial condition onto the elements
phi = numpy.zeros(3*sim_data.n_global_dof())

# Adjusts source term for modal basis functions
s = numpy.zeros(sim_data.n_global_dof())
if sim_data.basis_type == BasisType.MODAL or sim_data.basis_type == BasisType.MODAL_ORTHOGONAL:
    idx = sim_data.n_local_dof()*mat_prop.el
    s[idx:idx+sim_data.n_local_dof()] = basis.psi[mat_prop.locid, :]<jupyter_output><empty_output><jupyter_text>Auxiliary structures<jupyter_code># %load local_matrices.py
# 
# Local mass and stiffness matrices
# =================================
#
# by <NAME>. (Matlab and Python versions)
#    <NAME> (Python version)
#
import numpy

def local_mass(quad_data, basis):
    """Constructs the elemental mass matrix

    Arguments:
        quad_data - Quadrature points and weights
        basis - Basis and respective derivatives

    Returns:
        Mass matrix M, where m_ij = \int_k psi_i psi_j
    """
    return numpy.dot(quad_data.w*basis.psi.T, basis.psi)

def local_mass_diagonal(quad_data, basis):
    """Constructs the elemental mass matrix, diagonal version

    Arguments:
        quad_data - Quadrature points and weights
        basis - Basis and respective derivatives

    Returns:
        Mass matrix M, where m_ii = \int_k psi_i psi_i
    """
    return numpy.sum(quad_data.w*basis.psi.T**2, axis=1)

def local_stiffness(quad_data, basis):
    """Constructs the elemental stiffness matrix

    Arguments:
        quad_data - Quadrature points and weights
        basis - Basis and respective derivatives

    Returns:
        Stiffness matrix S, where s_ij = \int_k psi_i Dpsi_j/dx
    """
    return numpy.dot(quad_data.w*basis.psi.T, basis.dpsi)

def stress_stiffness(quad_data, basis, mat_prop, n):
    """Constructs the elemental stress stiffness matrix

    Arguments:
        quad_data - Quadrature points and weights
        basis - Basis and respective derivatives
        n - the element number
    """
    return (numpy.dot(quad_data.w*mat_prop.mu[:, n]*basis.psi.T, basis.dpsi).T, \
            numpy.dot(quad_data.w*mat_prop.dmudx[:, n]*basis.psi.T, basis.psi).T)

#-- local_matrices.py ----------------------------------------------------------

# %load eval_k_wave.py
# 
# Runge-Kutta 4 elemental evaluation (wave equation version)
# ==========================================================
#
# by <NAME>. (Matlab and Python versions)
#    <NAME> (Python version)
#
from simulation_data import BasisType
from local_matrices import local_mass, local_stiffness, stress_stiffness
from flux import Fluxes
import numpy

def eval_k_wave(sim_data, quad_data, basis, mesh, mat_prop, phi, s, t):
    ng = sim_data.n_global_dof()
    nl = sim_data.n_local_dof()
    K  = sim_data.n_elements   # saving some typing...

    # Allocates local vectors v and sigma
    v_t = numpy.zeros([(K+2)*nl, 2])     # size: K elements + 1 ghost
    sigma_t = numpy.zeros(v_t.shape) # size: K elements + 1 ghost

    # Gets v and sigma from previous stage of RK4
    v_t[nl:(K+1)*nl, 0] = phi[ng:2*ng]
    sigma_t[nl:(K+1)*nl,0] = phi[2*ng:3*ng]

    # Applies Neumann BC on v (v_ghost = v_element1)
    # d/dx v(x_bound) = 0 -> (v_element1 - v_ghost)/delta_x = 0 -> v_ghost = v_element1
    v_t[0:nl, 0] = v_t[nl:2*nl, 0]            # left side of the domain
    v_t[ng+nl:ng+2*nl, 0] = v_t[ng:ng+nl, 0]  # right side of the domain
 
    # Gets Mass, Stiffness and Flux matrices
    # Note that we multiply S by 2*a
    M = local_mass(quad_data, basis)
    S = local_stiffness(quad_data, basis)
    flux = Fluxes(sim_data)

    # Loop over total number of elements
    for i in range(1, sim_data.n_elements+1):
        Ss, Ssm = stress_stiffness(quad_data, basis, mat_prop, i-1)

        idx = numpy.arange(nl*i, nl*(i+1))
        m1 = 0.5*flux.FRp1.dot(sigma_t[idx+nl, 0]+numpy.sqrt(mat_prop.rho[i-1]*mat_prop.mub[nl-1, i])*v_t[idx+nl, 0] )
        m2 = 0.5*flux.FR.dot(sigma_t[idx, 0]-numpy.sqrt(mat_prop.rho[i-1]*mat_prop.mub[nl-1, i])*v_t[idx, 0])
        m3 = 0.5*flux.FL.dot(sigma_t[idx, 0]+numpy.sqrt(mat_prop.rho[i-1]*mat_prop.mub[0, i])*v_t[idx, 0])
        m4 = 0.5*flux.FLm1.dot(sigma_t[idx-nl, 0]-numpy.sqrt(mat_prop.rho[i-1]*mat_prop.mub[0, i])*v_t[idx-nl, 0])

        tmp = (s[idx-nl]-S.T.dot(sigma_t[idx, 0])+m1+m2-m3-m4)/(mesh.J[i-1]*mat_prop.rho[i-1])

        v_t[idx, 1] = numpy.linalg.solve(M, \
           tmp )

        m1 = 0.5*flux.FRp1.dot(v_t[idx+nl, 0]+(mat_prop.rho[i-1]*mat_prop.mub[nl-1, i])**(-0.5)*sigma_t[idx+nl, 0] )
        m2 = 0.5*flux.FR.dot(v_t[idx, 0]-(mat_prop.rho[i-1]*mat_prop.mub[nl-1, i])**(-0.5)*sigma_t[idx, 0])
        m3 = 0.5*flux.FL.dot(v_t[idx, 0]+(mat_prop.rho[i-1]*mat_prop.mub[0, i])**(-0.5)*sigma_t[idx, 0])
        m4 = 0.5*flux.FLm1.dot(v_t[idx-nl, 0]-(mat_prop.rho[i-1]*mat_prop.mub[0, i])**(-0.5)*sigma_t[idx-nl, 0])

        tmp = (-Ss.dot(v_t[idx, 0])-mesh.J[i-1]*Ssm.dot(v_t[idx, 0])+ \
                mat_prop.mub[nl-1, i]*(m1+m2)+mat_prop.mub[0, i]*(-m3-m4))/mesh.J[i-1]
        sigma_t[idx, 1] = numpy.linalg.solve(M,  tmp)

    # Assigns local vectors u, v and sigma to the new global vector phi
    ki = numpy.zeros(3*ng)

    ki[0:ng]      = v_t[nl:(K+1)*nl, 0]
    ki[ng:2*ng]   = v_t[nl:(K+1)*nl, 1]
    ki[2*ng:3*ng] = sigma_t[nl:(K+1)*nl, 1]

    return ki

#-- eval_k_wave.py -------------------------------------------------------------<jupyter_output><empty_output><jupyter_text>Finally, the time iteration!Load the `rk4_wave.py` file and be happy!<jupyter_code># %load rk4_wave.py
# 
# Runge-Kutta 4 steps method (wave equation version)
# ==================================================
#
# by <NAME>. (Matlab and Python versions)
#    <NAME> (Python version)
#
from simulation_data import BasisType
from eval_k_wave import eval_k_wave

def RK4_wave(sim_data, quad_data, basis, mesh, mat_prop, tstep, t, dt, u, src, theta, phi):
    """Runge-Kutta 4 steps.

    Arguments:
        sim_data - SimulationData object
        quad_data - Quadrature points and weights
        basis - Basis and respectivederivatives
        mesh - 1D Finite Elements mesh
        tstep - present time step (integer)
        t - simulation time
        dt - delta t
        u - solution
        s -
        phi -

    Returns:
        A step of the RK4 method
    """
    ngdof = sim_data.n_global_dof()
    nldof = sim_data.n_local_dof()

    # Sets source term
    s = src.copy()

    if sim_data.basis_type == BasisType.NODAL:
        s[mat_prop.globid] = theta[tstep-1]
    else:
        s[nldof*mat_prop.el:nldof*(mat_prop.el+1)] = s[nldof*mat_prop.el:nldof*(mat_prop.el+1)]*theta[tstep-1]

    # RK4
    k1 = eval_k_wave(sim_data, quad_data, basis, mesh, mat_prop, phi, s, t)
    k2 = eval_k_wave(sim_data, quad_data, basis, mesh, mat_prop, phi+0.5*dt*k1, s, t)
    k3 = eval_k_wave(sim_data, quad_data, basis, mesh, mat_prop, phi+0.5*dt*k2, s, t)
    k4 = eval_k_wave(sim_data, quad_data, basis, mesh, mat_prop, phi+dt*k3, s, t)
    
    phi = phi+dt*(k1+2.0*k2+2.0*k3+k4)/6.0

    # Recovers displacement solution u for modal type basis
    if sim_data.basis_type == BasisType.NODAL:
        u[:, tstep] = phi[0:ngdof]
    else:
        nip = sim_data.nip()
        nldof = sim_data.n_local_dof()

        u_hat = phi[0:ngdof]

        for k in range(0, sim_data.n_elements):
            u[nip*k:nip*(k+1), tstep] = \
                u_hat[nldof*k:nldof*(k+1)].dot(basis.psi.T)

    return (u, phi)

#-- rk4_wave.py ----------------------------------------------------------------

# Time marching (RK4 time integrator)
t = sim_data.t_min()

for tstep in range(1, Nsteps):
    u, phi = RK4_wave(sim_data, quad_data, basis, mesh, mat_prop, tstep, t, dt, u, s, theta, phi)
    t = t+dt<jupyter_output><empty_output><jupyter_text>Viewing the solutionIf you have installed the JSAnimation library you can generate a small movie showing the evolution of the solution along the time. If you get a message __`ImportError: No module named JSAnimation.IPython_display`__ do not dispair! We can still resort to matplotlib to plot a few steps.  Creating a movie<jupyter_code>from matplotlib import animation
from JSAnimation.IPython_display import display_animation
# Lets do some movies!
# From https://jakevdp.github.io/blog/2012/08/18/matplotlib-animation-tutorial/

# First set up the figure, the axis, and the plot element we want to animate
fig = pyplot.figure()
ax = pyplot.axes(xlim=(sim_data.spatial_domain[0], sim_data.spatial_domain[1]),ylim=(numpy.min(u), numpy.max(u)))
ax.grid(True)
line, = ax.plot([], [], lw=2)

# Initialization function: plot the background of each frame
def init():
    line.set_data([], [])
    return line,

# Animation function.  This is called sequentially
def animate(i):
    global u, mesh
    x = mesh.x.T.flat
    y = u[:, i]
    line.set_data(x, y)
    return line,

# call the animator.  blit=True means only re-draw the parts that have changed.
anim = animation.FuncAnimation(fig, animate, init_func=init,frames=Nsteps, interval=20, blit=True)
display_animation(anim, default_mode='once')

#pyplot.show()<jupyter_output><empty_output><jupyter_text>Plotting some steps of the solutionThe solution _u_ has _Nsteps_ columns, precisely one for each time step. In order to view some is simple, as shown in the code below.<jupyter_code># Plotting a few samples of the solution along the time
pyplot.figure(5)
pyplot.clf()
x = mesh.x.T.flat
pyplot.plot(x, u[:,150])  # tstep = 150
pyplot.plot(x, u[:,300])  # tstep = 300
pyplot.plot(x, u[:,600])  # tstep = 600
pyplot.plot(x, u[:,2100])  # tstep = 900
pyplot.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.20100107487205815, 'fraction_numerical': 0.016861052732902317, 'mean_word_length': 3.60267296758473, 'pattern_counts': {'":': 0, '<': 72, '<?xml version=': 0, '>': 77, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '503427', 'license': 'CC-BY-3.0', 'n_tokens_mistral': 22007, 'n_tokens_neox': 21051, 'n_words': 5149, 'path': 'lessons/08_dg/08_03_Elastic_Wave_1D.ipynb', 'repo_name': 'albertonogueira/numerical-mooc', 'size': '62561'}