% beamer version of chapter 10 "The Spectral Element Method" in "Computational Seismology, A Practical Introduction" (H. Igel)

\documentclass[aspectratio=169]{beamer}

\mode<presentation>
{
  %\usetheme{metropolis}
	\usetheme{metropolis}
  %\usecolortheme{beaver}
	\usecolortheme{crane}

  \setbeamercovered{transparent}
  \setbeamertemplate{navigation symbols}{}
}



\usepackage[english]{babel}
\usepackage[latin1]{inputenc}

\usepackage{times}
\usepackage[T1]{fontenc}
% Or whatever. Note that the encoding and the font should match. If T1
% does not look nice, try deleting the line with the fontenc.

% Math packages
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{mathtools}
\usepackage{graphicx}
%using a different letter style for capital letters
\usepackage{mathrsfs}

\usepackage{mdframed}
\usepackage{listings}
\usepackage[percent]{overpic}
\usepackage{capt-of}
\usepackage{tikz}
\usepackage{grffile}

% Own definitions ... 
\newcommand{\eqdef}{\mathrel{\mathop:}=} % for symbol :=
\newcommand*\Eval[3]{\left.#1\right\rvert_{#2}^{#3}} % for integration bound |
%\newcommand{\figpath}{../../../TheBook/Part2/Figs_se/}
%\newcommand{\figpath}{c:/mydata/Heiner/Text/Books/OUP/TheBook/Part2/Figs_se/}
%\newcommand{\figpath}{c:/mydata/Heiner/Text/Books/OUP/Submission/Igel-ComputationalSeismology/Part2/Figs_se/}
\newcommand{\figpath}{c:/mydata/Heiner/Text/Books/OUP/TheBook/Part2/Figs_se}



\newcommand{\fs}{\footnotesize}

\title[Computational Seismology] % (optional, use only with long paper titles)
{The Spectral-Element Method}

%\subtitle{Introduction} 

\author[Heiner Igel]{Heiner Igel}


\institute[] 
  {Department of Earth and Environmental Sciences \\ Ludwig-Maximilians-University Munich} 



\date[] % (optional)
{}

\subject{Talks}


%\beamerdefaultoverlayspecification{<+->}

\begin{document}

% ----------------------------------------------------------

\begin{frame}
  \titlepage
\end{frame}


% ---------------------------------------------------------

\section{Introduction}

% ----------------------------------------------------------

\subsection{Motivation}

\begin{frame}{Motivtion}
\footnotesize


\begin{itemize}
\item  High accuracy for wave propagation problems
\item  Flexibility with Earth model geometries
\item  Accurate implementation of boundary conditions
\item  Efficient parallelization possible
\end{itemize}

\begin{figure}
\includegraphics[scale=.15]{\figpath/fig_bridge.jpg}
\end{figure}

\end{frame}

% -----------------------------------------------------------

\subsection{History}


\begin{frame}{History}

\footnotesize

%\frametitle{History}
\begin{itemize}[<+->]
\item  Originated in fluid dynamics (Patera, 1984, Maday and Patera, 1989)
\item  First applications to seismic wave propagation by Priolo, Carcione and Seriani, 1994)
\item  Initial concepts around Chebyshev polynomials (Faccioli, Maggio, Quarteroni and Tagliani, 1996)
\item  Use of Lagrange polynomials (diagonal mass matrix) by Komatitsch and Vilotte (1998)
\item  Application to spherical geometry using the cubed sphere concept (Chaljub and Vilotte, 2004) 
\item Spectral elements in spherical coordinates by Fichtner et al. (2009)
\item  Community code "specfem" widely used for simulaton and inversion (e.g., Peter et al. 2011). 
\end{itemize}

\end{frame}

% -----------------------------------------------------------

\begin{frame}{Specfem 3D}
\begin{minipage}{0.35\textwidth}
\begin{figure}
\includegraphics[width=0.95\textwidth]{\figpath/fig_specfem3D.jpg}
\end{figure}
\end{minipage}
	\mbox{ }
    \hfill	
\begin{minipage}{0.6\textwidth}
Logo of the wellknown community code with a snapshot of global wave propagation for a simulation of the devastating M9.1 earthquake near Sumatra in December, 2004. \\
The open-source code is hosted by the CIG project \\
(www.geodynamics.org). 
\end{minipage}
\end{frame}


% -----------------------------------------------------------
\subsection{Spectral Elements in a Nutshell}


\begin{frame}{Spectral Elements in a Nutshell}

{
\footnotesize

\begin{columns}[t]
\begin{column}{5cm}
\begin{figure}
\includegraphics[scale=.25]{\figpath/fig_sem_element_lp.png}
\end{figure}

\end{column}
\begin{column}{5cm}


\begin{itemize}


\item Snapshot of the displacement field $u$ during a simulation in a strongly heterogeneous medium. 
\item Zoom into the displacement field inside an element discretised with order N=5 collocation points 
\item Exact interpolate using Lagrange polynomials. 
\item Gauss-Lobatto-Legendre Numerical integration scheme. 

\end{itemize}


\end{column}
\end{columns}

}

\end{frame}

% -----------------------------------------------------------

%

\begin{frame}
\begin{block}{1-D elastic wave equation}
\begin{equation*}
\rho(x) \ddot{u}(x,t)  =  \partial_x \left[ \mu(x) \partial_x u(x,t) \right] \ + \ f(x,t)
\end{equation*}
\end{block}
\vspace{1cm}
\begin{itemize}
\item[] $u$ displacement 
\item[] $f$ external force 
\item[]  $\rho$ mass density 
\item[] $\mu$ shear modulus 
\end{itemize}


\end{frame}

% --------------------------------------------


\begin{frame}{Boundary condition}
No traction perpendicular to the Earth's free surface
\begin{equation*}
\sigma_{ij} \; n_j = 0
\end{equation*}
normal vector  $n_j$, $\sigma_{ij}$ is the symmetric stress tensor
\begin{equation*}
\mu \partial_x u(x,t) |_{x=0,L} = 0
\end{equation*}
where our spatial boundaries are at $x=0,L$ and the stress-free condition applies at both ends 
\end{frame}

% --------------------------------------------


\begin{frame}{Spectral Element: Essentials}
\begin{itemize}
\item Weak formulation of the wave equation
\item  Transformation to the elemental level (Jacobian)
\item  Approximation of unknown function $u$ using Lagrange polynomials 
\item Evaluation of the 1st derivatives of the Lagrange polynomials 
\item Numerical integration scheme based on GLL quadrature 
\item Calculation of system matrices at elemental level
\item Assembly of global system of equations
\item Extrapolation in time using a simple finite-difference scheme
\end{itemize}
\end{frame}


% -----------------------------------------------

\begin{frame}{Galerkin Principle}
\begin{itemize}
\item The underlying principle of the finite-element method
\item Developed in context with structural engineering (Boris Galerkin, 1871-1945)
\item Also developed by Walther Ritz (1909) - variational principle
\item Conversion of a continuous operator problem (such as a differential equation) to a discrete problem 
\item Constraints on a finite set of basis functions
\end{itemize}
\end{frame}




% ----------------------------------------------------------

\begin{frame}{Weak Formulations}
Multiplication of pde with test function $w(x)$ on both sides.\\
G is here the complete computational domain defined with $x \in G = \left[ 0, L \right]$. 
\begin{equation*}
\int\limits_{G} \ w \ \rho \ \ddot u \ \mathrm{d}x \ -  \int\limits_{G} w \ \partial_x \ (\mu \ \partial_x u) \ \mathrm{d}x \ = \
\int\limits_{G} w \ f \ \mathrm{d}x 
\end{equation*}
\end{frame}

% -------------------------------------------------------------------

\begin{frame}{Integration by parts}
\begin{equation*}
\int\limits_{G} \ w \ \rho \ \ddot u \ \mathrm{d}x \ +  \int_{G} \mu \ \partial_x w \ \partial_x \  u \ \mathrm{d}x =
\int\limits_{G} w \ f \ \mathrm{d}x
\end{equation*}
where we made use of the boundary condition 
\begin{equation*}
\partial_x u(x,t) |_{x=0} = \partial_x u(x,t) |_{x=L} = 0 \ \
\end{equation*}
\end{frame}

% -------------------------------------------------------------------

\begin{frame}
\begin{block}{The approximate displacement field}
\begin{equation*}
u(x,t) \; \approx \; \overline{u}(x,t) \; = \; \sum_{i=1}^{n} u_i(t) \; \psi_i(x).
\end{equation*}
\end{block}
\vspace{0.5cm}
\begin{itemize}
\item Discretization of space introduced with this step
\item Specific basis function not yet specified
\item In principle basis functions defined on entire domain \\ (-> PS method)
\item Locality of basis functions lead to finite-element type method
\end{itemize}
\end{frame}

% -------------------------------------------------------------------------


\subsection{Global system}

\begin{frame}{Global system after discretization}

Use approximation of $u$ in weak form and (!) use the same basis function as test function
\begin{equation*}
\int\limits_{G} \psi_i \rho \ddot{\overline{u}} \mathrm{d}x +  \int\limits_{G} \mu  \partial_x \psi_i  \partial_x  \overline{u} \mathrm{d}x =
\int\limits_{G} \psi_i f \mathrm{d}x
\end{equation*}
with the requirement that the medium is at rest a t=0. 
\end{frame}

% ------------------------------------------------------------------------------


\begin{frame}{Including the function approximation for $u(x,t)$}
... leads to an equation for the unknown coefficients $u_i(t)$
\begin{equation*}
\begin{split}
&\sum_{i=1}^n \left[   \ddot{u}_i(t) \int_{G_e} \rho(x) \ \psi_j(x) \ \psi_i(x)  \ \mathrm{d}x \right]  \\
+ & \sum_{i=1}^n \left[ u_i(t) \int_{G_e} \mu(x) \ \partial_x \psi_j(x) \ \partial_x \psi_i(x) \ \mathrm{d}x  \right] \\ 
= & \int_{G} \psi_i \ f(x,t) \ \mathrm{d}x
\end{split}
\end{equation*}
for all basis functions $\psi_j$ with $j=1,...,n$. This is the classical algebro-differential equation for {\bf finite-element} type problems. 
\end{frame}

% -----------------------------------------------------------------------------

\begin{frame}{Matrix notation}

This system of equations with  the coefficients of the basis functions (meaning?!) as unknowns can be written in matrix notation
\begin{block}{}
\begin{equation*}
\textbf{M} \cdot \ddot{\textbf {u}} (t)  \ + \textbf{K} \cdot \textbf{u} (t)  =  \textbf{f} (t) 
\end{equation*}
\end{block}
\vspace{0.5cm}
\begin{itemize}
\item[] \textbf{M} Mass matrix
\item[] \textbf{K} Stiffness matrix
\item[] terminology originates from structural engineering problems
\end{itemize}
\end{frame}

% -------------------------------------------------------------------------------------------------------


\begin{frame}{Matrix system graphically}

The figure gives an actual graphical representation of the matrices for our 1D problem. 

\begin{figure}
\includegraphics[scale=0.325]{\figpath/fig_sem_symbolic.png}
\end{figure}
The unknown vector of coefficients {\bf u} is found by a simple finite-difference procedure. The solution requires the inversion of mass matrix {\bf M} which is trivial as it is diagonal. That is the key feature of the spectral element method.

\end{frame}

% ----------------------------------------------------------------------------------------

\begin{frame}{Mass and stiffness matrix}
Definition of the - at this  point - global mass matrix
\begin{equation*} 
M_{ji} \ = \ \int\limits_{G} \rho(x) \ \psi_j(x) \ \psi_i(x)  \ \mathrm{d}x
\end{equation*}
and the stiffness matrix
\begin{equation*} 
K_{ji} \ = \ \int\limits_{G} \mu(x) \ \partial_x \psi_j(x) \ \partial_x \psi_i(x) \ \mathrm{d}x
\end{equation*}

and the vector containing the volumetric forces $f(x,t)$ 
\begin{equation*}
f_j(t) \ = \ \int\limits_{G} \psi_i \ f(x,t) \ \mathrm{d}x \ \ . 
\end{equation*}
\end{frame}

% ------------------------------------------------------------

\begin{frame}{Mapping}
A simple centred finite-difference approximation of the 2nd derivative and the following mapping 
\begin{equation*}
\begin{split}
 \textbf{u}^{new}  &\rightarrow \ \textbf{u}(t + \mathrm{d}t) \\
 \textbf{u} \ \   &\rightarrow \ \textbf{u}(t) \\
 \textbf{u}^{old}   &\rightarrow \ \textbf{u}(t - \mathrm{d}t) \\
\end{split}
\end{equation*}
leads us to the solution for the coefficient vector $\textbf{u}(t + \mathrm{d}t)$ for the next time step as already well known from the other solution schemes in previous schemes. 
\end{frame}

% ------------------------------------------------------------

\subsection{Solution scheme}

\begin{frame}{Solution scheme}
\begin{block}{}
\begin{equation*}
\textbf{u}^{new} \ = \mathrm{d}t^2 \ \left[ \ \textbf{M}^{-1} \ (\textbf{f} - \textbf{K \ u}) \ \right] \ + 2 \textbf{u} \ - \textbf{u}^{old}
\end{equation*}
\end{block}
\vspace{0.5cm}
\begin{itemize}
\item[] General solution scheme for finite-element method (wave propagation)
\item[] Independent of choice of basis functions
\item[] Mass matrix needs to be inverted!
\item[] To Do: good choice of basis function, integration scheme for calculation of M and K
\end{itemize}
\end{frame}



% ----------------------------------------------------------

\begin{frame}{Element level}
\begin{figure}
\includegraphics[scale=.35]{\figpath/fig_element_mapping.png}
\end{figure}
\fs In order to facilitate the calculation of the space-dependent integrals we transform each element onto the standard interval [-1,1], illustrated here for $n_e=3$ elements. The elements share the boundary points. 

\end{frame}

% ------------------------------------------------------------

\begin{frame}{System at element level} 
\begin{equation*}
\begin{split}
& \sum_{i=1}^{n} \left[ \ddot{u}_i(t) \sum_{e=1}^{n_e} \int\limits_{G_e} \rho (x) \psi_j (x) \psi_i (x) \mathrm{d}x \right] \\
+& \sum_{i=1}^{n} \left[ u_i (t) \sum_{e=1}^{n_e} \int\limits_{G_e} \mu (x) \partial_x \psi_j (x) \partial_x \psi_i (x)\mathrm{d}x \right] \\
=& \sum_{e=1}^{n_e} \int\limits_{G_e} \psi_j (x) f(x, t) \mathrm{d}x  \ \ .
\end{split}
\end{equation*}
\small Because of the sum over $n_e$ there is a global dependence of the coefficients. How can we avoid this?  
\end{frame}

% ---------------------------------------------------

\begin{frame}{Local basis functions}
\begin{columns}[t]
\begin{column}{5cm}
\begin{figure}
\includegraphics[scale=.08]{\figpath/fig_localbasis.png}
\end{figure}
\end{column}

\begin{column}{5cm}
\vskip 0.5cm
\fs Illustration of local basis functions. By defining basis functions only inside elements the integrals can be evaluated in a local coordinate system. The graph assumes three elements of length two. A finite-element type linear basis function (dashed line) is shown along side a spectral-element type Lagrange polynomial basis function of degree N=5 (solid line).
\end{column}
\end{columns}

\end{frame}

% -----------------------------------------------------

\begin{frame}{$u$ inside element}
\begin{block}{}
\begin{equation*}
\Eval{\overline{u} (x,t)}{ x \in G_e}{} \ = \ \sum_{i=1}^{n} u_i^e (t) \psi_i^e (x)
\end{equation*}
\end{block}
\vspace{0.5cm}
\begin{itemize}
\item[] Now we can proceed with all calculations locally in $G_e$ (inside one element)
\item[] Here is the difference to pseudospectral methods
\item[] Sum is now over all basis functions inside one element (n turns out to be the order of the polynomial scheme)
\end{itemize}
\end{frame}

% -------------------------------------------------
 
\begin{frame}{Local system of equations}

\begin{equation*}
\begin{split}
& \sum_{i=1}^{n} \ddot{u}_i^e (t) \int\limits_{G_e} \rho (x) \psi_j^e (x) \psi_i^e (x) \mathrm{d}x \\
+ &\sum_{i=1}^{n}  u_i^e (t) \int\limits_{G_e} \mu (x) \partial_x \psi_j^e (x) \partial_x \psi_i^e (x)\mathrm{d}x  \\
= & \int\limits_{G_e} \psi_j^e (x) f(x, t) \mathrm{d}x \ \ .
\end{split}
\end{equation*}
\end{frame}

% -------------------------------------------

\begin{frame}{Matrix notation for local system}
\begin{block}{}
\begin{equation*}
\textbf{M}^e \cdot \ddot{\textbf{u}}^e (t) + \textbf{K}^e \cdot \textbf{u}^e (t) \ = \ \textbf{f}^e (t), 
\qquad e \ = \ 1, \dotsc, n_e 
\end{equation*}
\end{block}
\vspace{0.5cm}
Here $\textbf{u}^e$, $\textbf{K}^e$, $\textbf{M}^e$, and  $\textbf{f}^e$ are the coefficients of the unknown displacement inside the element, stiffness and mass matrices with information on the density and stiffnesses, and the forces, respectively.
\end{frame}

% -------------------------------------------


\begin{frame}{Coordinate transformation}
\begin{equation*}
F_e \ \colon \left[ -1, 1 \right] \ \rightarrow  \ G_e,  \qquad x \ = \ F_e(\xi), \qquad 
\end{equation*}
\begin{equation*}
 \xi \ = \ \xi(x) \ = \ F_e^{-1}(\xi), \qquad e \ = \ 1, \dotsc , n_e 
\end{equation*}
from our global system $x \in G$ to our local coordinates that we denote $\xi \in F_e$

\begin{equation*}
x(\xi) \ = \ F_e(\xi) \ = \ \Delta e \frac{(\xi + 1)}{2} + x_e
\end{equation*}
where $n_e$ is the number of elements, and $\xi \in [-1,1]$. Thus the physical coordinate $x$ can be related to the local coordinate $\xi$ via
\begin{equation*}
\xi(x) \ = \ \frac{2(x - x_e)}{\Delta e} - 1 
\end{equation*}

\end{frame}

% --------------------------------------------------

\begin{frame}{Integration of arbitrary function}
\begin{equation*}
\int\limits_{G_e}f(x) \mathrm{d}x \ = \ \int\limits_{-1}^1 f^e(\xi)\frac{\mathrm{d}x}{\mathrm{d}\xi}\mathrm{d}\xi \ 
\end{equation*}
the integrand has to be multiplied by the Jabobian $J$ before integration. 
\begin{equation*}
 \qquad  J \ = \ \frac{\mathrm{d}x}{\mathrm{d}\xi} \ = \ \frac{\Delta e}{2} \ \ .
\end{equation*}
we will also need
\begin{equation*}
\qquad J^{-1} \ = \ \frac{\mathrm{d}\xi}{\mathrm{d}x} \ = \ \frac{2}{\Delta e}  \ . 
\end{equation*}

\end{frame}

% -------------------------------------------------

\begin{frame}{Skewed 3D elements}

\begin{columns}[t]
\begin{column}{5cm}
\begin{figure}
\includegraphics[scale=.4]{\figpath/fig_Jacobian3D.png}
\end{figure}
\end{column}

\begin{column}{5cm}
\vskip 1cm
In 3D elements might be skewed and have curved boundaries. The calculation of the Jacobian is then carried out analytically by means of shape functions.
\end{column}
\end{columns}

\end{frame}

% ---------------------------------------------------

\begin{frame}{Assembly of system of equations}
\begin{equation*}
\begin{split}
   & \sum_{i=1}^{n}\ddot{u}_i^e(t) \int\limits_{-1}^1 \rho \left[x(\xi) \right] \psi_j^e\left[x(\xi)\right] \psi_i^e\left[x(\xi)\right] \frac{\mathrm{d}x}{\mathrm{d}\xi}\mathrm{d}\xi \\
+  &\sum_{i=1}^{n}u_i^e(t) \int\limits_{-1}^1 \mu \left[x(\xi)\right] \frac{\mathrm{d}\psi_j^e\left[x(\xi)\right]}{\mathrm{d}\xi} \frac{\mathrm{d}\psi_i^e\left[x(\xi)\right]}{\mathrm{d}\xi}  \left(\frac{\mathrm{d}\xi}{\mathrm{d}x} \right)^2 \frac{\mathrm{d}x}{\mathrm{d}\xi}\mathrm{d}\xi \\
=  &\int\limits_{-1}^1 \psi_j^e \left[x(\xi)\right] f\left[(x(\xi)), t\right] \frac{\mathrm{d}x}{\mathrm{d}\xi}\mathrm{d}\xi \ \ .
\end{split}
\end{equation*}
\fs System of $n$ equations for each index $j$ corresponding to one particular basis function. We need to find basis functions that allows efficient and accurate calculation of the required integrals. 

\end{frame} 

% -------------------------------------------------------

\subsection{Interpolation with Lagrange Polynomials} 

% -------------------------------------------------------

\begin{frame}{Lagrange polynomials}
Remember we seek to approximate $u(x,t)$ by a sum over space-dependent basis functions $\psi_i$ weighted by time-dependent coefficients $u_i(t)$.
\begin{equation*}
u(x,t) \ \approx \ \overline{u}(x,t) \ = \ \sum_{i=1}^{n} u_i(t) \ \psi_i(x) 
\end{equation*}
Our final choice: Lagrange polynomials:
\begin{equation*}
\psi_i \ \rightarrow \ \ell_i^{(N)} \ \eqdef \ \prod_{k = 1, \ k \neq i}^{N+1} \frac{\xi - \xi_k}{\xi_i-\xi_k}, \qquad   i = 1, 2, \dotsc , N + 1 
\end{equation*}
where $x_i$ are fixed points in the interval $[-1,1]$.
\end{frame}

% --------------------------------------------------------

\begin{frame}{Orthogonality of Lagrange polynomials}
\begin{equation*}
{\bf \ell_i^{(N)}(\xi_j) = \delta_{ij}}
\end{equation*} \\
\vspace{0.5cm}
They fulfill the condition that they {\it exactly} interpolate (or approximate) the function at  N+1 collocation points. Compare with discrete Fourier series on regular grids or Chebyshev polynomials on appropriate grid points (-> pseudospectral method). 
\end{frame}

% -----------------------------------------------------------

\begin{frame}{Lagrange polynomials graphically}
\begin{figure}
    \begin{minipage}{0.4\textwidth}
       \includegraphics[width=0.6\textwidth]{\figpath/fig_lp_2-6.png}
    \end{minipage} 
    	\mbox{ }
    \hfill	
    \begin{minipage}{0.55\textwidth}  
    \small
{\bf Top:} Family of $N+1$ Lagrange polynomials for $N=2$ defined in the interval $\xi \in [-1,1]$. Note their maximum value in the whole interval does not exceed unity. \\
{\bf Bottom:}  Same for $N=6$. The domain is divided into N intervals of uneven length. When using Lagrange polynomials for function interpolation the values are exactly recovered at the collocation points (squares).
    \end{minipage}
\end{figure}
\end{frame}


% ---------------------------------------------------------

\begin{frame}{Gauss-Lobatto-Legendre points}
\begin{figure}
\includegraphics[scale=.2, angle=0]{\figpath/fig_gll2.png}
\end{figure}
Illustration of the spatial distribution of Gauss-Lobatto-Legendre points in the interval [-1,1] from top to bottom for polynomial order 2 to 12 (from left to right). Note the increasing difference of largest to smallest interval between collocation points! Consequences?
\end{frame}

% ----------------------------------------------------------

\begin{frame}{Lagrange polynomials: some properties}

Mathematically the collocation property is expressed as 
\begin{equation*}
\ell_i^{(N)}(\xi_i)=1 \ \ \textrm{and}   \ \ \dot{\ell}_i^{(N)}(\xi_i) \ = \ 0
\end{equation*}
where the dot denotes a spatial derivative. The fact that
\begin{equation*}
\vert \ell_i^{(N)}(\xi) \vert \leq 1, \qquad \xi \in \left[ -1, 1 \right]  
\end{equation*}
minimizes the interpolation error in between the collocation points due to numerical inaccuracies. 
\end{frame}

% ----------------------------------------------------------

\begin{frame}{Function approximation}

This is the final mathematical description of the unknown field $u(x,t)$ for the spectral-element method based on Lagrange polynomials.  
\begin{equation*}
u^e(\xi) \ = \ \sum_{i=1}^{N+1}u^e(\xi_i)\ell_i(\xi)
\end{equation*}
Other options at this point are the Chebyhev  polynomials. They have equally good approximation properties (but ...) 

\end{frame}


% -----------------------------------------------------------

\begin{frame}{Interpolation with Lagrange Polynomials}

\begin{minipage}{0.4\textwidth}
\begin{figure}
       \includegraphics[width=0.55\textwidth]{\figpath/fig_lpint.png}
\end{figure}
\end{minipage}
    	\mbox{ }
    \hfill	
\begin{minipage}{0.55\textwidth}
The function to be approximated is given by the solid lines. The approximation is given by the dashed line exactly interpolating the function at the GLL points (squares). {\bf Top:} Order $N=2$ with  three grid points. {\bf Bottom:} Order $N=6$ with seven grid points.
\end{minipage}
\end{frame}

% --------------------------------------------------------


\begin{frame}{Spectral-element system with basis functions}
Including the Legendre polynomials in our local (element-based) system leads to 
\begin{equation*}
\begin{split}
&  \sum_{i =  1}^{N+1}\ddot{u}_i^e(t)  \int\limits_{-1}^{1} \rho (\xi) \ell_j(\xi) \ell_i(\xi) \frac{\mathrm{d}x}{\mathrm{d}\xi} \mathrm{d}\xi \\
+ & \sum_{i,k  =  1}^{N+1} u_i^e(t) \int\limits_{-1}^{1}    \mu (\xi) \dot{\ell}_j(\xi) \dot{\ell}_i(\xi) \left( \frac{\mathrm{d}\xi}{\mathrm{d}x} \right)^2 \frac{\mathrm{d}x}{\mathrm{d}\xi}  \mathrm{d}\xi \\
 = & \int\limits_{-1}^{1} \ell_j(\xi) f (\xi, t) \frac{\mathrm{d}x}{\mathrm{d}\xi}  \mathrm{d}\xi  
\end{split}
\end{equation*}
Because we want that $\rho$, $\mu$, and $f$ vary inside one element there is no way out carrying out the integration numerically. 

\end{frame}

% ------------------------------------------------------------------

\subsection[Numerical Integration]{Numerical Integration: The Gauss-Lobatto-Legendre Approach}

\begin{frame}{Integration scheme for an arbitrary function $f(x)$}
\begin{equation*}
\int\limits_{-1}^1 f(x) \mathrm{d}x \ \approx \ \int\limits_{-1}^1 P_N(x) \mathrm{d}x  \ = \ \sum_{i =  1}^{N+1} w_i f(x_i)
\end{equation*}
defined in the interval $x \in [-1,1]$ with 
\begin{equation*}
 P_N(x) \ = \ \sum_{i = 1}^{N+1}f(x_i)\ell_i^{(N)}(x)
\end{equation*}
\begin{equation*}
w_i \ = \ \int\limits_{-1}^1 \ell_i^{(N)} (x) \mathrm{d}x \ \ .
\end{equation*}
\end{frame}

% ---------------------------------------------------------------------

\begin{frame}{Numerical integration scheme}

\begin{minipage}{0.5\textwidth}
\begin{figure}
\includegraphics[width=0.4\textwidth]{\figpath/fig_integration.png}
\end{figure}
\end{minipage}
    	\mbox{ }
    \hfill	
\begin{minipage}{0.4\textwidth}
%\tiny
\begin{itemize}
\item Exact function (thick solid line)
\item Approximation by Lagrange polynomials (thin solid line)
\item  Difference between  true and approximate function (light gray)
\end{itemize}
\end{minipage}
\end{frame}

% --------------------------------------------------------------------

\begin{frame}{Collocation points and integration weights}

\begin{center}
\begin{tabular}[t]{lll}
$N$&$\xi_i$&$\omega_i$\\
\hline
2:&
0 & 4/3\\
&$\pm$ 1& 1/3 \\
&&\\
3:&
$\pm \sqrt{1/5}$ & 5/6  \\
&$\pm$ 1& 1/6 \\
&&\\
4:&
0&32/45\\
&$\pm \sqrt{3/7}$ & 49/90 \\
&$\pm$ 1& 1/10 \\
\end{tabular}\label{gll_weights}
\end{center}
Collocation points and integration weights of the Gauss-Lobatto-Legendre quadrature for order
$N=2,\dots,4$.
\end{frame}

% ------------------------------------------------------------------

\begin{frame}{After integration}

With the numerical integration scheme we obtain
\begin{equation*}
\begin{split}
& \Eval{\sum_{i,k  =  1}^{N+1}\ddot{u}_i^e(t)w_k \rho (\xi) \ell_j(\xi) \ell_i(\xi) \frac{\mathrm{d}x}{\mathrm{d}\xi}} {\xi = \xi_k}{} \\
+ & \Eval{\sum_{i,k  =  1}^{N+1}w_k u_i^e(t) \mu(\xi) \ell_j(\xi) \ell_i(\xi) \left( \frac{\mathrm{d}\xi}{\mathrm{d}x} \right)^2 \frac{\mathrm{d}x}{\mathrm{d}\xi}} {\xi = \xi_k}{} \\
\approx & \Eval{\sum_{i,k  =  1}^{N+1}w_k \ell_j(\xi) f (\xi, t) \frac{\mathrm{d}x}{\mathrm{d}\xi}} {\xi = \xi_k}{}
\end{split}
\end{equation*}
What is still missing is a formulation for the derivative of the Lagrange polynomials at the collocation points.
But: Major progress! We have replaced the integrals by sums!

\end{frame}

% ---------------------------------------------------------------

\begin{frame}{cont'd...}

\begin{block}{Solution equation for our spectral-element system at the element level}
\begin{equation*}
\sum_{i = 1}^{N+1} M_{ji}^e \ddot{u}_i^e(t) + \sum_{i = 1}^{N+1} K_{ji}^e u_i^e(t) \ = \ f_j^e(t), 
\qquad e \ = \ 1, \dotsc , n_e
\end{equation*}
\end{block}
\vfill
with
\begin{equation*}
\begin{split}
 M_{ji}^e& \ = \Eval{\ w_j \rho (\xi) \frac{\mathrm{d}x}{\mathrm{d}\xi} \delta_{ij} }{\xi = \xi_j}{} \\
 K_{ji}^e& \ = \ \Eval{\sum_{k = 1}^{N+1} w_k \mu (\xi) \ell_j (\xi)  \ell_i (\xi) \left(\frac{\mathrm{d}\xi}{\mathrm{d}x} \right)^2 \frac{\mathrm{d}x}{\mathrm{d}\xi} }{\xi = \xi_k}{} \\
 f_j^e& \ = \ \Eval{w_j f (\xi,t) \frac{\mathrm{d}x}{\mathrm{d}\xi} }{\xi = \xi_j}{}
\end{split} 
\end{equation*}

\end{frame}

% ----------------------------------------------------------------

\subsection{Derivatives of the Lagrange Polynomials}

% ---------------------------------------------------------------

\begin{frame}{Illustration of Legendre Polynomials}

\begin{minipage}{0.45\textwidth}
\begin{figure}
\includegraphics[width=0.95\textwidth]{\figpath/fig_legendre.png}
\end{figure}
\end{minipage}
    	\mbox{ }
    \hfill	
\begin{minipage}{0.5\textwidth}
The Legendre polynomials are used to calculate the first derivatives of the Legendre polynomials. They can also be used to calculate the integration weights of the GLL quadrature.

\begin{equation*}
\dot{\ell}(\xi_i) \ = \ \sum_{j=1}^{N} \tilde{d}_{ij}^{(1)} \ell(\xi_j)
\end{equation*}
\begin{equation*}
\partial_xu^e(\xi) \ = \ \sum_{i=1}^{N+1}u^e(\xi_i)\dot{\ell}_i(\xi)   \ \ 
\end{equation*}
\end{minipage}
\end{frame}

% --------------------------------------------------------------------------


\section{Global Assembly and Solution}
\subsection{Global Assembly}

% --------------------------------------------------------------------------

\begin{frame}{Global Assembly for the diagonal of the mass matrix}

\begin{equation*}
\textbf{M}_g = \left(\begin{array}{c}
M_{1,1}^{(1)}\\
M_{2,2}^{(1)}\\
M_{3,3}^{(1)}\\
\\
\\
\\
\\
\end{array}
\right)
+
\left(\begin{array}{c}
\\
\\
M_{1,1}^{(2)}\\
M_{2,2}^{(2)}\\
M_{3,3}^{(2)}\\
\\
\\
\end{array}
\right)
+
\left(\begin{array}{c}
\\
\\
\\
\\
M_{1,1}^{(3)}\\
M_{2,2}^{(3)}\\
M_{3,3}^{(3)}\\
\end{array}
\right)
=
\left(\begin{array}{c}
M_{1,1}^{(1)}\\
M_{2,2}^{(1)}\\
M_{3,3}^{(1)} + M_{1,1}^{(2)}\\
M_{2,2}^{(2)}\\
M_{3,3}^{(2)} + M_{1,1}^{(3)}\\
M_{2,2}^{(3)}\\
M_{3,3}^{(3)}\\
\end{array}
\right)
\end{equation*}
\end{frame}  



\begin{frame}{Global Assembly for the diagonal of the stiffness matrix}
\small
\begin{eqnarray*}
\textbf{K}_g & = &
\left(
\begin{array}{ccccccc}
K_{1,1}^{(1)}&K_{1,2}^{(1)}&K_{1,3}^{(1)} \ &&&&\\
K_{2,1}^{(1)}&K_{2,2}^{(1)}&K_{2,3}^{(1)} \  &&&\mathbf{0}&\\
K_{3,1}^{(1)}&K_{3,2}^{(1)}&K_{3,3}^{(1)}+K_{1,1}^{(2)}
&  K_{1,2}^{(2)}&  K_{1,3}^{(2)} \  &&\\
  &  &    \ K_{2,1}^{(2)} &K_{2,2}^{(2)} &K_{2,3}^{(2)} \  &&\\
&&  \ K_{3,1}^{(2)} &K_{3,2}^{(2)} &K_{3,3}^{(2)}+K  _{1,1}^{(3)} 
&  K_{1,2}^{(3)}&  K_{1,3}^{(3)}\\
&\mathbf{0}&&&  \ K_{2,1}^{(3)}&K_{2,2}^{(3)}&K_{2,3}^{(3)}\\
&&&&  \ K_{3,1}^{(2)}&K_{3,2}^{(2)}&K_{3,3}^{(2)}
\end{array}
\right)
\end{eqnarray*}

\end{frame}

% ----------------------------------------------------------------------------------------

\begin{frame}{Vector with information on the source}
\begin{equation*}
\textbf{f}_g = 
\left(\begin{array}{c}
f_{1}^{(1)}\\
f_{2}^{(1)}\\
f_{3}^{(1)} + f_{1}^{(2)}\\
f_{2}^{(2)}\\
f_{3}^{(2)} + f_{1}^{(3)}\\
f_{2}^{(3)}\\
f_{3}^{(3)}\\
\end{array}
\right)
\end{equation*}
\end{frame}

% -------------------------------------------------------------------------------------

\subsection{Global Solution}
\begin{frame}{Extrapolation for time-dependent coefficients $\textbf{u}_g$}

This is our final algorithm as it is implemented using Matlab or Python
\begin{equation*} 
\begin{split}
\textbf{u}_{g}(t+\mathrm{d}t) \ = \ & \mathrm{d}t^2 \ \left[ \ \textbf{M$_{g}$}^{-1} \ (\textbf{f$_{g}$}(t) - \textbf{K$_{g}$ \ u$_{g}$}(t)) \ \right] \ \\
+ &  \ 2 \textbf{u$_{g}$(t)} \ - \textbf{u$_{g}$}(t-\mathrm{d}t) 
\end{split}
\end{equation*}
Looks fairly simple, no?  
\end{frame}

% ---------------------------------------------------------------------------------------

\begin{frame}{Matlab Code: Matrices}

{\sffamily
{\footnotesize
\% Elemental Mass Matrix \\
for i = 1 : N + 1, \\
\indent  Me(i) = rho(i) * w(i) * J; \\
end   \\
(...)\\
\% Elemental Stiffness Matrix \\
for i = 1 : N + 1, \\
  \hskip 1cm  for j = 1 : N + 1, \\
   \hskip 2cm           sum = 0; \\
   \hskip 2cm          for k = 1 : N + 1, \\
   \hskip 2cm          sum = sum + mu(k) * w(k) * Ji\^{}2 * J * l1d(i,k) * l1d(j,k); \\
   \hskip 2cm          end \\ 
   \hskip 2cm           Ke(i,j) = sum; \\
   \hskip 1cm    end \\
end}}

\end{frame}

% ---------------------------------------------------------------------------------------

\begin{frame}{Matlab Code: Time extrapolation}

for it = 1 : nt, \\
(...) \\
  \hskip 1cm  \% Extrapolation \\
  \hskip 1cm    unew = dt\^{}2 * (Minv * (f - K * u)) + 2 * u - uold; \\
  \hskip 1cm    uold = u; \\
  \hskip 1cm   u = unew;\\
 (...)     \\
end 


\end{frame}

% ---------------------------------------------------------------------------------------

\begin{frame}{Spectral elements: work flow}

\fs

\begin{figure}
\includegraphics[scale=.35]{\figpath/fig_sem_workflow.png}
\end{figure}
\tiny A substantial part consists of preparing the interpolation and integration procedures required to initialize the global mass- and stiffness matrices. The final time-extrapolation is extremely compact and does not require the inversion of a global matrix as is the case in classic finite-element methods.

\end{frame}

% -------------------------------------------------------------------------------

\begin{frame}{Point source initialization}
\fs
\begin{columns}[t]
\begin{column}{5cm}

\begin{figure}
\includegraphics[scale=.12]{\figpath/fig_force.png}
\end{figure}

\end{column}

\begin{column}{5cm}
\vskip.5cm
{\bf Top:} A point-source polynomial representation (solid line) of a $\delta$-function (red bar). \\ {\bf Bottom:} A finite-source polynomial representation (solid line) as a superposition of point sources injected at some collocation points. For comparison with analytical solutions it is important to note that the spatial source function actually simulated is the polynomial representation of the (sum over) $\delta$-function(s).


\end{column}

\end{columns}

\end{frame}

% ------------------------------------------------------

\begin{frame}{sem1d: simulation examples}

\begin{figure}
\includegraphics[width=0.4\textwidth]{\figpath/fig_sem1d_comp.png}
\end{figure}
\end{frame}

% ------------------------------------------------------

\begin{frame}{Summary}
\begin{itemize}
\item The spectral element method combines the flexibility of finite-element methods concerning computational meshes with the spectral convergence of Lagrange basis functions used inside the elements.
\item The enormous success of the spectral element method is based upon the diagonal structure of the mass matrix that needs to be inverted to extrapolate the system in time. Due to this diagonality no matrix inversion techniques need to be employed allowing straight forward parallelisation of the algorithm. The diagonal mass matrix is possible through the coincidence of the collocation points of both interpolation and integration schemes (Gauss-Lobatto-Legendre).
\end{itemize}
\end{frame}

% ------------------------------------------------------

\begin{frame}{Summary}
\begin{itemize}
\item The errors of the spectral-element scheme accumulate from the (usually low-order finite-difference) time extrapolation scheme and the numerical integration using Gauss-Lobatto-Legendre quadrature.
\item The spectral-element method is particularly useful for simulation problems where the free-surface plays an important role, and/or in which surface waves need to be accurately modelled. Technically the reason is that the free-surface boundary is implicitly solved. 
\item A well engineered community code (SPECFEM3D, www.geodynamics.org) is available for Cartesian and spherical geometries including global Earth (or planetary scale) calculations. 
\end{itemize}
\end{frame}

% --------------------------------------------------------------------------

\end{document}



%\begin{frame}{Derivatives of the Lagrange Polynomials}
%\end{frame}
%
%
%\begin{frame}{Calculation of the derivatives}
%We need to know the derivatives only on the GLL collocation points by using Legendre polynomials: \\
%\begin{equation*}
%P_N(\xi) \ = \ \frac{1}{2^NN!} \frac{\mathrm{d}^N}{\mathrm{d}\xi^N}\left(\xi^2 - 1\right)^N
%\end{equation*}
%where $N$ denotes the polynomial degree. 
%\end{frame}
%
%
%\begin{frame}{Calculation with recursive formula}
%The Legendre polynomials can be calculated using the following recursive formula
%\begin{eqnarray*}
%    P_0(\xi) & = & 1  \\
%    P_1(\xi) & = & \xi \\
%    P_{n\geq2}(\xi) & = & \frac{1}{n}[(2n-1) \ \xi \ P_{n-1}(\xi)-(n-1)P_{n-2}(\xi)] 
%\end{eqnarray*}
%\end{frame} 
%
%
%\begin{frame}{Derivatives of Lagrange polynomials}
%\begin{equation*}
%\dot{\ell}(\xi_i) = \sum_{j=1}^{N} \tilde {\mathrm{d}}_{ij}^{ \ (1)} \ell(\xi_j)
%\end{equation*}
%with
%\begin{eqnarray*}
% \tilde {\mathrm{d}}_{ij}^{ \ (1)} & = & \left\{
%\begin{array}{ll}
%-\frac{1}{4} N(N+1) & \mbox{if}\ i  =  j  =  0\\
%&\\
%\frac{P_N(\xi_i)}{P_N(\xi_j)}\frac{1}{\xi_i - \xi_j} & \mbox{if} \ 0\leq i \leq N, \
% 0\leq j \leq N, \ i\neq j \\
%&\\
%0 & \mbox{if} \ 1 \leq i=j \leq N-1\\
%&\\
%\frac{1}{4} N(N+1) & \mbox{if} \ i=j=N 
%\end{array}
%\right . 
%\end{eqnarray*}
%\end{frame}
%
%
%\begin{frame}{Approximation of derivatives of an arbitrary function $u(\xi)$ known at the collocation points $\xi_i$}
%\end{frame}


