% beamer version of chapter 4 "The Finite Volume 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}
\usepackage{mathrsfs}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%Implement Matlab and Python Codes%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\usepackage{listings}
\lstloadlanguages{bash,  Python,  Matlab}
\lstset{
alsolanguage=bash, %Setzen der Codesprache
alsolanguage=Python,
alsolanguage=Matlab,
basicstyle=\tiny \sffamily , % alle listings winzig drucken
keywordstyle=\color{blue} \bfseries , % Schlüsselwörter fett und schwarz drucken
identifierstyle=\color{black},
stringstyle=\color{orange}, % Strings werden orange gedruckt
%otherkeywords={*,+,-}=\color{green}, % andere Schlüsselwörter, die gekennzeichnet werden sollen
commentstyle=\color{purple!40!black}, % Kommentare blau drucken
%stringstyle=\ttfamily, % Strings im Code Schreibmaschinenähnlich - setzt sich etwas vom Code ab
backgroundcolor=\color{white!95!black},
showstringspaces=false, % Strings im Code ohne Kenntlichmachung von Leerzeichen
captionpos=b,   %Caption erscheint unterhalb des Codes
frame=single, %Code wird umrandet
numbers=none, % Nummerierung der Linien links vom Code
numberstyle=\footnotesize, % Art der Nummerierung, hier z. B.: winzige Schrift
%stepnumber=5, % Intervall der Nummern, hier z.B. 5
              % falls nicht gesetzt: Alle Nummern werden ausgegeben.
numbersep=9 pt % Abstand vom Code, hier z.B. 5 pt
} 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\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_fe}
\newcommand{\figpath}{c:/mydata/Heiner/Text/Books/OUP/TheBook/Part2/Figs_fe}

\newcommand{\fs}{\footnotesize}

\title[Computational Seismology] % (optional, use only with long paper titles)
{The Finite 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{The Finite-Element Method: 1D Elastic Wave Equation}

\begin{frame}{1D Elastic Wave Equation}
	Apply the Galerkin principle to the 1D elastic wave equation 
\begin{equation*}
	\rho \partial_t^2 u = \partial_x \; \mu \; \partial_x u + f \ . 
\end{equation*}
where again we omit space and time-dependencies. From now on we assume that the properties of the medium, density $\rho$, and shear modulus $\mu$ are both space-dependent. We obtain the weak form as
\begin{equation*}
	\int_D \rho \; \partial_t^2 \; u \; \varphi_j \; dx = \int_D \partial_x \; \mu \; \partial_x \; u \; \varphi_j \; dx + \int_D f \; \varphi_j \; dx \ . 
\end{equation*} 
Integration by parts of the term containing the space derivatives leads to 
\begin{equation*}
	\int_D  \partial_x \; \mu \; \partial_x \; u \; \varphi_j \; dx = \left[ \mu \partial_x \; u \; \varphi_j  \right] - \int_D  \mu \; \partial_x \; u \; \partial_x \; \varphi_j \; dx
\end{equation*}
\end{frame}

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

\begin{frame}{1D Elastic Wave Equation}
Assuming a stress-free boundary condition leads to
\begin{equation*}
	\int_D \rho \; \partial_t^2 \; u \; \varphi_j \; dx + \int_D  \mu \; \partial_x \; u \; \partial_x \; \varphi_j \; dx = \int_D f \; \varphi_j \; dx
\end{equation*}
where  $u$ is the continuous unknown displacement field. We replace the exact displacement field by an approximation $\overline{u}$ of the form
\begin{equation*}
	u(x,t) \rightarrow \overline{u} (x,t) = \sum_{i=1}^N u_i (t) \; \varphi_i (x)
\end{equation*}
where the coefficients $u_i$ are expected to correspond to a discrete representation of the solution field. The wave equation bevomes
\begin{equation*}
	\int_D \rho \; \partial_t^2 \; \overline{u} \; \varphi_j \; dx + \int_D  \mu \; \partial_x \;\overline{u} \; \partial_x \; \varphi_j \; dx = \int_D f \; \varphi_j \; dx
\end{equation*}
\end{frame}

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

\begin{frame}{1D Elastic Wave Equation}
Turning the continuous weak form into a system of linear equations
\begin{equation*}
\begin{split}
	&\int_D \rho \; \partial_t^2 \; \left( \sum_{i=1}^N u_i (t) \; \varphi_i (x) \right) \; \varphi_j \; dx \\ 
	+& \int_D  \mu \; \partial_x \; \left( \sum_{i=1}^N u_i (t) \; \varphi_i (x) \right) \; \partial_x \; \varphi_j \; dx \\
	=& \int_D f \; \varphi_j \; dx	\ .
\end{split}
\end{equation*}
Changing the order of integration and summation we obtain
\begin{equation*}
	\sum_{i=1}^N \partial_t^2 u_i \int_D \rho \; \varphi_i \; \varphi_j \; dx + \sum_{i=1}^N u_i \int_D \mu \; \partial_x \; \varphi_i \; \partial_x \; \varphi_j \; dx = \int_D f \; \varphi_j \; dx
\end{equation*}
using the fact that the unknown coefficients $u_i$ only depend on time. 
\end{frame}

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

\begin{frame}{1D Elastic Wave Equation}
Using matrix-vector notation with the following definitions for the time-dependent solution vector of displacement values $\textbf{u} (t)$, mass matrix $\textbf{M}$, the already well-known stiffness matrix $\textbf{K}$, and the source vector $\textbf{f}$:
\begin{equation*}
	\begin{split}
		\textbf{u} (t) & \rightarrow u_i (t) \\
		\textbf{M}& \rightarrow M_{ij} = \int_D \rho \; \varphi_i \; \varphi_j \; dx \\
		\textbf{K} & \rightarrow K_{ij} = \int_D \mu \; \partial_x \; \varphi_i \; \partial_x \; \varphi_j dx \\
		\textbf{f} & \rightarrow f_j = \int_D f \; \varphi_j \; dx \ .
	\end{split}
\end{equation*}
\end{frame}

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

\begin{frame}{1D Elastic Wave Equation}
Thus we can write the system of equations as
\begin{equation*}
	\ddot{\textbf{u}} \textbf{M} + \textbf{u} \textbf{K} = \textbf{f} 
\end{equation*}
or with transposed system matrices as
\begin{equation*}
	\textbf{M}^T \ddot{\textbf{u}} + \textbf{K}^T \textbf{u} = \textbf{f} \ . 
\end{equation*}
For the second time-derivative we use a standard finite-difference approximation 
\begin{equation*}
	\ddot{\textbf{u}} = \partial_t^2 \textbf{u} \approx \frac{\textbf{u} (t+dt) - 2 \textbf{u} (t) + \textbf{u} (t-dt)}{dt^2}
\end{equation*}
\end{frame}

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

\begin{frame}{1D Elastic Wave Equation}
Replacing the original partial derivative with respect to time to obtain
\begin{equation*}
	\textbf{M}^T \left[  \frac{\textbf{u} (t+dt) - 2 \textbf{u} (t) + \textbf{u} (t-dt)}{dt^2} \right] = \textbf{f} - \textbf{K}^T \textbf{u}  \ . 
\end{equation*}
Starting from an initial state $\textbf{u}(t=0)=0$ we can determine the displacement field at time $t+dt$ by 
\begin{equation*}
\textbf{u} (t+dt) = dt^2 (\textbf{M}^T)^{-1} \left[ \textbf{f} - \textbf{K}^T \textbf{u}   \right] + 2 \textbf{u} (t) - \textbf{u} (t - dt) \ . 
\end{equation*}
\end{frame}

% ---------------------------------------------------------
\subsection{The System Matrices}

\begin{frame}{The System Matrices}
	To calculate the entries of the system matrices we transform the space coordinate into a local system
\begin{equation*}
	\begin{split}	
	\xi & = x - x_i \\
	h_i & = x_{i+1} - x_i \ . \\
	\end{split}
\end{equation*}
However, now we allow the element size $h_i$ to vary. With the definition above element $i$ is defined in the interval $x \in  [x_i, x_{i+1}]$. In the local coordinate system the basis functions are defined by
\begin{equation*}
	\varphi_i (\xi) = 
	\begin{cases}
	\frac{\xi}{h_{i-1}} \;\; \; \; \text{for} \; - h_{i-1} < \xi \leq 0 \\
	1 - \frac{\xi}{h_i} \; \text{for} \; 0 < \xi < h_i \\
	0 \; \; \; \; \; \; \; \; \; \; \text{elsewhere}
	\end{cases}
\end{equation*}

\end{frame}

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

\begin{frame}{The System Matrices}
	\begin{minipage}{0.45\textwidth}
		\begin{figure}
			\includegraphics[width=0.9\textwidth]{\figpath/fig_fe_basis_h.png}
		\end{figure}			
	\end{minipage}
    \mbox{ }
    \hfill		
	\begin{minipage}{0.5\textwidth}
with the corresponding derivatives 
\begin{equation*}
	\partial_x \; \varphi_i (\xi) = 
	\begin{cases}
		\frac{1}{h_{i-1}} \; \text{for} \; -h_{i-1} < \xi \leq 0 \\
		- \frac{1}{h_i} \; \text{for} \; 0 < \xi < h_i \\
		\phantom{-} 0 \; \; \;  \; \; \text{elsewhere} \\
	\end{cases}
\end{equation*}
Example of a finite-element domain with irregular element sizes $h_i$. The basis functions (thick solid lines) are illustrated with the normalized derivatives (thin solid lines).
	\end{minipage}
\end{frame}

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

\begin{frame}{The Mass Matrix}
Looking at the global definition of the mass matrix $\textbf{M}$ with components
\begin{equation*}
 M_{ij} = \int_D \rho \; \varphi_i \; \varphi_j \; dx 
\end{equation*}
the only non-zero entries are around the diagonal and are of components $M_{i , i-1}$ $M_{i,i}$ and  $M_{i,i+1}$ for $i=2,...,N-1$. Elements $M_{11}$ and $M_{NN}$ have to be treated separately. For the diagonal elements we obtain
\begin{equation*}
	M_{ii} = \int_D \rho \; \varphi_i \; \varphi_i \; dx \ = \  \int_{D_{\xi}} \rho \; \varphi_i \; \varphi_i \; d \xi 
\end{equation*}
in the local coordinate system. 
\end{frame}

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

\begin{frame}{The Mass Matrix}
Integration has to be carried out over the elements to the left and right of the boundary points $x_i$. We thus obtain
\begin{equation*}
\begin{split}
M_{ii} &= \rho_{i-1} \int_{-h_{i-1}}^0 \left( \frac{\xi}{h_{i-1}} + 1  \right)^2 d \xi + \rho_i \int_0^{h_i} \left( 1 - \frac{\xi}{h_i} \right)^2 d \xi \\
	   &= \frac{1}{3} \left( \rho_{i-1} \; h_{i-1} + \rho_i \; h_i \right) \ 
\end{split}
\end{equation*}
For the off-diagonal elements the basis functions overlap only in one element
\begin{equation*}
	M_{i,i-1} = \rho_{i-1}  \int_{-h_{i-1}}^0 \left( \frac{\xi}{h_{i-1}} + 1  \right) \frac{-\xi}{h_{i-1}} d \xi = \frac{1}{6} \; \rho_{i-1} \; h_{i-1}
\end{equation*}
or
\begin{equation*}
M_{i,i+1} = \rho_i \int_0^{h_i} \frac{\xi}{h_i} \left( 1 - \frac{\xi}{h_i} \right) d \xi = \frac{1}{6} \; \rho_i \; h_i \ . 
\end{equation*}
\end{frame}

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

\begin{frame}{The Mass Matrix}
Just to illustrate the banded nature of the mass matrix, assuming constant element size $h$ and density $\rho$ the mass matrix is given by
\begin{equation*}
	\textbf{M} = \frac{\rho \; h}{6}
	\begin{pmatrix}
		 \ddots &  &   &   &  \smash{\makebox{\hspace{3mm} 0\hspace{-3mm}} }  & \\
		1 & 4  & 1  &   &   &   \\
		   & 1  & 4  & 1  &   &   \\
		   &   & 1  & 4  & 1  &   \\
		 \smash{\makebox{\hspace{3mm} 0\hspace{-3mm}} } &   &   &   & \ddots  &   \\
	\end{pmatrix}
\end{equation*}\\
\vspace{0.5cm}
\textbf{Note:} In the general case with varying element size the mass matrix is not symmetric. 
\end{frame}

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

\begin{frame}{Stiffness matrix}
The same concepts apply to the stiffness matrix. We move to the local coordinate system by
\begin{equation*}
	K_{ij} = \int_D \mu \; \partial_x \; \varphi_i \; \partial_x \; \varphi_j \; dx \ =\  \int_{D_{\xi}} \mu \; \partial_{\xi} \; \varphi_i \; \partial_{\xi} \; \varphi_j \; d \xi
\end{equation*}
to obtain for a diagonal element, assuming constant shear modulus $\mu$ inside each element
\begin{equation*}
	\begin{split}
	K_{ii} & = \mu_{i-1} \int_{-h_{i-1}}^0 \left( \frac{1}{h_{i-1}} \right)^2 d \xi + \mu_i \int_0^{h_i} \left( - \frac{1}{h_i} \right)^2 d \xi \\
	& = \frac{\mu_{i-1}}{h_{i-1}} + \frac{\mu_i}{h_i}
	\end{split}
\end{equation*}
\end{frame}

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

\begin{frame}{Stiffness matrix}
For the off-diagonal elements
\begin{equation*}
	\begin{split}
	K_{i,i+1} & = \mu_i \int_0^{h_i} \left( - \frac{1}{h_i} \right) \left( \frac{1}{h_i} \right) d \xi = - \frac{\mu_i}{h_i} \\
	K_{i,i-1} & = \mu_{i-1} \int_{-h_{i-1}}^0 \left( - \frac{1}{h_{i-1}} \right) \left( \frac{1}{h_{i-1}} \right) d \xi = - \frac{\mu_{i-1}}{h_{i-1}} \\
	\end{split}
\end{equation*}
while all other elements of the stiffness matrix are zero. For example, assuming constant shear modulus and element size the stiffness matrix reads
\begin{equation*}
	\textbf{K} = \frac{\mu}{h}
	\begin{pmatrix}
		 & \ddots  &   &   &  \smash{\makebox{\hspace{3mm} 0\hspace{-3mm}} }  &   \\
		 & -1  & 2  & -1  &   &   \\
		 &   & -1  & 2  & -1  &   \\
		 &   &   & -1  & 2  & -1  \\
		 &  \smash{\makebox{\hspace{3mm} 0\hspace{-3mm}} }  &   &   & \ddots  &   \\
	\end{pmatrix} \ . 
\end{equation*}
\end{frame}

% ---------------------------------------------------------
\subsection{Simulation Example: Comparison with Finite Differences}

\begin{frame}[fragile]{Simulation Example}
	\begin{minipage}{0.3\textwidth}
		\begin{table}
			\begin{tabular}[t]{ll}
Parameter & Value \\
\hline
$x_{max}$ & 10000m   \\
nx & 1000 \\
$v_s$ & 3000 m/s \\
$\rho$ & 2500 $kg/m^3$ \\
h & 10 m \\
eps & 0.5 \\
$f_0$ & 20 Hz
\end{tabular}
\end{table}				
	\end{minipage}
    \mbox{ }
    \hfill		
	\begin{minipage}{0.58\textwidth}
\begin{lstlisting}
% Time extrapolation
for i=1:nt,   
   % Finite-difference method
   unew=dt^2 * Minv * ((D*u)+f/h*src(i)) + 2*u - uold; 
   uold=u;
   u=unew;
   (...)
end
\end{lstlisting}
Matlab code segment for a finite-difference algorithm in matrix-vector form.
	\end{minipage}	
\end{frame}

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

\begin{frame}{Simulation Example}
	\begin{figure}
		\includegraphics[width=0.45\textwidth]{\figpath/fig_fe_matrices.png}
	\end{figure}
\footnotesize{The structure of the system matrices for the finite-element method are compared with the finite-difference method formulated with matrix-vector operations. {\bf Top row:} Stiffness and inverse mass matrix for the finite element method. {\bf Bottom row:} Stiffness (differential) matrix and diagonal mass matrix for the finite-difference method.}
\end{frame}

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

\begin{frame}[fragile]{Simulation Example}
	\begin{minipage}{0.5\textwidth}
	Mass matrix for the general case of varying element size
\begin{lstlisting}
% Mass matrix Mij
M=zeros(nx);
for i=2:nx-1,
   for j=2:nx-1,
      if i==j, 
         M(i,j)=1./3*(ro(i-1)*h(i-1)+ro(i)*h(i));
      elseif j==i+1
         M(i,j)=ro(i)*h(i)/6;
      elseif j==i-1
         M(i,j)=ro(i-1)*h(i-1)/6;
      else
         M(i,j)=0;
      end
   end
end
% Corner elements
M(1,1)=1/3*ro(1)*h(1);
M(nx,nx)=1/3*ro(nx)*h(nx-1);
(...)
\end{lstlisting}				
	\end{minipage}
    \mbox{ }
    \hfill		
	\begin{minipage}{0.45\textwidth}
	Finite element time extrapolation
\begin{lstlisting}
% Time extrapolation
for i=1:nt,   
   % Finite element method
   unew=dt^2 * Minv * (f*src(i)-K*u) 
   + 2*u - uold;                             
   uold=u;
   u=unew;
   (...)
end
\end{lstlisting}
	\end{minipage}
\end{frame}

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

\begin{frame}[fragile]{Simulation Example}
	\begin{minipage}{0.5\textwidth}
		\begin{figure}
			\includegraphics[width=0.3\textwidth]{\figpath/fig_fe_example.png}
\end{figure}				
	\end{minipage}
    \mbox{ }
    \hfill		
	\begin{minipage}{0.45\textwidth}
Snapshots of the displacement wavefield calculated with the finite-element method (solid line) are compared with the finite-difference method (dotted line) at various distances from the source using the same parameters. The length of the window is 500m.
	\end{minipage}
\end{frame}

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

\begin{frame}{h - adaptivity}
h - adaptivity: The simplicity with which the element size can vary due to geometrical features or the velocity model. \\
Think about an Earth model in which the seismic velocities have strong variations. 
\begin{itemize}
	\item[] P-velocity in the oceans (1.5 km/s)
	\item[] P-velocities at the core-mantle boundary (13 km/s)
	\item[] Any numerical scheme with globally constant element size has to be accurate for the shortest wavelength
	\item[] $\Longrightarrow$ Regions with higher velocities will be oversampled
\end{itemize}
\end{frame}

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

\begin{frame}{h - adaptivity}
	\begin{minipage}{0.5\textwidth}
		\begin{table}
			\begin{tabular}[t]{llll}
 & Left  & Middle & Right \\
\hline
x & 4600m & 1000m & 4600m \\
$v_s$ & 6000$\frac{m}{s}$ & 1500$\frac{m}{s}$ & 3000$\frac{m}{s}$\\
dx & 40m & 10m & 20m \\
$\rho$ & 2500$\frac{kg}{m^3}$ & 2500$\frac{kg}{m^3}$ &2500$\frac{kg}{m^3}$  \\
\\
\end{tabular}
\end{table}	
\begin{table}
\begin{tabular}[t]{ll}
Parameter & Value \\
\hline
nt & 18000 \\
dt & 3.3 ms \\
$f_0$ & 5 Hz \\
eps & 0.5 
\end{tabular}
\end{table}		
	\end{minipage}
    \mbox{ }
    \hfill		
	\begin{minipage}{0.45\textwidth}
		We demonstrate this in the 1D case with a strongly heterogeneous velocity model in which the number of grid points per wavelength is kept constant in the entire physical domain. The model mimics the situation in a fault zone with a central low-velocity zone (damage zone) with different material properties on the two sides of the fault.
	\end{minipage}
\end{frame}

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

\begin{frame}{h - adaptivity}
	\begin{figure}
		\includegraphics[width=0.4\textwidth]{\figpath/fig_fe_example_h.png}
	\end{figure}
	\footnotesize{Snapshots of displacement values are shown as a function time. Where displacement amplitudes are below a threshold, the velocity models is shown in gray scale. Note the polarity change of the reflections at the boundaries and the slope of the signals in the $x-t$ plane indicated their velocities.}
\end{frame}

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

\begin{frame}{h - adaptivity}
	\begin{figure}
		\includegraphics[width=0.4\textwidth]{\figpath/fig_fe_example_h_detail.png}
	\end{figure}
	\footnotesize{Detail of the finite-element simulation with varying element size at one of the domain boundaries. The crosses indicate the element boundaries and the changing element size. Note the continuous but non-differentiable behaviour of the displacement field at the interface.}
\end{frame}

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

\section{Shape Functions in 1D and 2D}


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

\begin{frame}{1D}
Let us recall how we replaced the originally continuous unknown field $u(x)$ by a sum over some basis functions $\varphi_i$
\begin{equation*}
u(x) = \sum_{i=1}^N c_i \varphi_i (x)
\end{equation*}
denoting the coefficients of the basis functions by $c_i$. Mapping all elements to a local coordinate system such
\begin{equation*}
\xi = \frac{x-x_i}{x_{i+1}-x_i}
\end{equation*}
where our reference element is defined with $\xi \in [0, 1]$.
\end{frame}

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

\begin{frame}{Linear Shape Functions}
We put ourselves at element level and assume that our unknown function $u(\xi)$ is linear
\begin{equation*}
u(\xi) = c_1 + c_2 \xi
\end{equation*}
where $c_i$ are real coefficients. Each element has two node points, namely the element boundaries at $\xi_{1,2}=0,1$. This leads to the following conditions and solutions for coefficients $c_i$
\begin{equation*}
\begin{split}
u_1 = c_1 &\rightarrow c_1 = u_1 \\
u_2 = c_1 + c_2  &\rightarrow c_2 = -u_1 + u_2 
\end{split} \ . 
\end{equation*}
\end{frame}

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

\begin{frame}{Linear Shape Functions}
This can also be written in matrix notation which will help us when dealing with high-order systems. We obtain
\begin{equation*}
\left[ \begin{array}{l} u_1 \\ u_2  \\ \end{array} \right] =
\left[ \begin{array}{ll} 1 & 0 \\ 1 & 1 \\ \end{array} \right]    
\left[ \begin{array}{l} c_1 \\ c_2  \\ \end{array} \right] 
\end{equation*}
and using matrix inversion
\begin{equation*}
\left[ \begin{array}{l} c_1 \\ c_2  \\ \end{array} \right] =
\left[ \begin{array}{ll} 1 & 0 \\ -1 & 1 \\ \end{array} \right] 
\left[ \begin{array}{l} u_1 \\ u_2  \\ \end{array} \right]     \ . 
\end{equation*}
\end{frame}


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

\begin{frame}{Linear Shape Functions}
With appropriate matrix and vector definitions this can be written as
\begin{equation*}
\textbf{u} = \textbf{A}\textbf{c} \rightarrow \textbf{c} = \textbf{A}^{-1}\textbf{u}
\end{equation*}
implying that to obtain coefficients \textbf{c} we need to calculate the inverse of \textbf{A}.
\begin{equation*}
\begin{split}
u(\xi) & = u_1 + (-u_1 + u_2) \xi \\
       & = u_1(1-\xi) + u_2 {\xi} \\
       & = u_1 N_1(\xi) + u_2 N_2(\xi) 
\end{split}
\end{equation*}
where we introduced a novel concept, the shape functions $N_i(\xi)$ with the following form
\begin{equation*}
N_1(\xi) = 1 - \xi \ \  , \ \  N_2(\xi) = \xi \ . 
\end{equation*}
\end{frame}

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

\begin{frame}{Linear Shape Functions}
	\begin{minipage}{0.4\textwidth}
		\begin{figure}
			\includegraphics[width=0.7\textwidth]{\figpath/fig_fe_shape1d.png}
\end{figure}				
	\end{minipage}
    \mbox{ }
    \hfill		
	\begin{minipage}{0.55\textwidth}
The sum over the weighted shape function of general order $N$ gives the approximate {\it continuous} representation of the solution field $u(\xi)$ inside the element.
\begin{equation*}
u(\xi) = \sum_{i=1}^N u_i N_i(\xi)
\end{equation*}
{\bf Top:} Linear shape functions as used in the development of the finite-element solution to static and dynamic elastic problems. Node points are indicated by crosses. {\bf Bottom:} Quadratic shape functions requiring one more node point at the center of the element.
	\end{minipage}
\end{frame}

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

\begin{frame}{Quadratic Shape Functions}
Describing our solution field by quadratic functions requires
\begin{equation*}
u(\xi) = c_1 + c_2 \xi + c_3 \xi^2
\end{equation*}
where we added one more node point at the center of the element $\xi_{1,2,3} = 0, \frac{1}{2}, 1$. With these node locations we obtain
\begin{equation*}
\begin{split}
u_1 & = c_1 \\
u_2 & = c_1  + 0.5 c_2 + 0.25 c_3 \\
u_3 & = c_1 + c_2 + c_3
\end{split}
\end{equation*}
and after inverting the resulting system matrix \textbf{A}
\begin{equation*}
\textbf(A)^{-1} = 
\left[
\begin{array}{lll}
1 & 0 & 0 \\
-3 & 4 & -1 \\
2 & -4 & 2 \\
\end{array}
\right]
\end{equation*}
\end{frame}

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

\begin{frame}{Quadratic Shape Functions}
We can represent the final quadratic solution field inside the element with
\begin{equation*}
\begin{split}
u(\xi) =& c_1 + c_2 \xi + c_3 \xi^2 \\
       =& u_1 (1-3\xi +2\xi^2) + \\
       =& u_2 (4\xi - 4 \xi^2) + \\
       =& u_3 (-xi + 2 \xi^2) 
\end{split}
\end{equation*}
resulting in the following shape functions
\begin{equation*}
\begin{split}
N_1(\xi) &= 1-3\xi +2\xi^2 \\
N_2(\xi) &= 4\xi - 4 \xi^2 \\
N_3(\xi) &= -\xi + 2 \xi^2
\end{split}
\end{equation*}
\end{frame}

% ---------------------------------------------------------
\subsection{2D}

\begin{frame}{Triangular Shape Functions}
	\begin{minipage}{0.4\textwidth}
		\begin{figure}
			\includegraphics[width=0.7\textwidth]{\figpath/fig_fe_tri1.png}
\end{figure}				
	\end{minipage}
    \mbox{ }
    \hfill		
	\begin{minipage}{0.55\textwidth}
Triangular elements. Mapping of physical coordinates ($x,y$, top) to a local reference frame (bottom) with coordinates $\xi,\eta$.
	\end{minipage}
\end{frame}

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

\begin{frame}{Triangular Shape Functions}
To perform the integration operations when calculating the system matrices we move to the local coordinate system $\xi,\eta \in [0, 1]$ (sometimes the reference space is chosen to be $[-1, 1]$) through 
\begin{equation*}
\begin{split}
x &= x_1 + (x_2 - x_1) \xi + (x_3 - x_1) \eta \\
y &= y_1 + (y_2 - y_1) \xi + (y_3 - y_1) \eta \ .
\end{split}
\end{equation*}
We seek to describe a linear function inside our triangle, therefore
\begin{equation*}
u(\xi,\eta) = c_1 + c_2 \xi + c_3 \eta \ . 
\end{equation*}
\end{frame}

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

\begin{frame}{Triangular Shape Functions}
We only know our function at the corners of the reference triangle, therefore the constraints for coefficients $c_i$ are
\begin{equation*}
\begin{split}
u_1 &= u(0,0) = c_1 \\
u_2 &= u(1,0) = c_1 + c_2 \\
u_3 &= u(0,1) = c_1 + c_3  \ \ . 
\end{split}
\end{equation*}
This leads - using the same matrix inversion approach described above to - the following shape functions for triangular elements
\begin{equation*}
\begin{split}
N_1(\xi,\eta) &= 1 - \xi - \eta \\
N_2(\xi,\eta) &= \xi \\
N_3(\xi,\eta) &= \eta \ . 
\end{split}
\end{equation*}
\end{frame}

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

\begin{frame}{Triangular Shape Functions}
	\begin{figure}
		\includegraphics[width=0.9\textwidth]{\figpath/fig_fe_tri2.png}
	\end{figure}
Triangular elements. The three corner nodes lead to an equivalent number of shape functions $N_i (\xi,\eta)$ with unit value at one of the corners.
\end{frame}

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

\begin{frame}{Rectangular Shape Functions}
	\begin{minipage}{0.4\textwidth}
		\begin{figure}
			\includegraphics[width=0.7\textwidth]{\figpath/fig_fe_rect1.png}
\end{figure}				
	\end{minipage}
    \mbox{ }
    \hfill		
	\begin{minipage}{0.55\textwidth}
Quadrilateral elements. Mapping of physical coordinates ($x,y$, top) to a local reference frame (bottom) with coordinates $\xi,\eta$.
	\end{minipage}
\end{frame}

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

\begin{frame}{Rectangular Shape Functions}
Accordingly, shape functions can be derived for general quadrilateral elements. We map space to a local coordinate system through
\begin{equation*}
\begin{split}
x &= x_1 + (x_2 - x_1) \xi + (x_4 - x_1) \eta + (x_3 - x_2) \xi \eta \\
y &= y_1 + (y_2 - y_1) \xi + (y_4 - y_1) \eta + (y_3 - y_2) \xi \eta  \ . 
\end{split}
\end{equation*}
Requiring linear behaviour of the function inside the element 
\begin{equation*}
u(\xi,\eta) = c_1 + c_2 \xi + c_3 \eta + c_4 \xi \eta
\end{equation*}
we obtain the following shape functions
\begin{equation*}
\begin{split}
N_1(\xi,\eta) &= (1 - \xi) (1 - \eta) \\
N_2(\xi,\eta) &= \xi ( 1  - \eta ) \\
N_3(\xi,\eta) &= \xi \eta \\
N_4(\xi,\eta) &= (1 - \xi) \eta \ \ . 
\end{split}
\end{equation*}
\end{frame}

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

\begin{frame}{Rectangular Shape Functions}
	\begin{figure}
		\includegraphics[width=0.5\textwidth]{\figpath/fig_fe_rect2.png}
	\end{figure}
Rectangular reference elements. The four shape functions $N_i (\xi,\eta)$ with unity value at one of the corners.
\end{frame}

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

\begin{frame}{Summary}
	\footnotesize{
	\begin{itemize}
		\item The finite-element method was originally developed  mostly for static structural engineering problems.
		\item The {\it element} concept relates to describing the solution field in an analogous way inside each element, thereby facilitating the required calculations of the system matrices.
		\item The finite-element approach can in principle be applied to elements of arbitrary shape. Most used shapes are triangles (tetrahedra) or quadrilateral (hexahedral) structures.
		\item The finite-element method is a series expansion method. The continuous solution field is replaced by a finite sum over (not necessarily orthogonal) basis functions.
		\item For static elastic problems or the elastic wave propagation problem finite-element analysis leads to a (large) system of linear equations. In general, the matrices are of size $N \times N$ where N is the number of degrees of freedom.
		\item Because of the specific interpolation properties of the basis functions, their coefficients take the meaning of the values of the solution field at specific node points.
	\end{itemize}
	}
\end{frame}

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

\begin{frame}{Summary}
	\footnotesize{
	\begin{itemize}
		\item In an initialization step the global stiffness and mass matrices have to be calculated. They depend on integrals over products of basis functions and their derivatives.
		\item If equation parameters (e.g., elastic parameters, density) vary inside elements, then numerical integration has to be performed.
		\item The stress-free surface condition can be implicitly solved. This is a major advantage for example for the simulation of surface waves.
		\item The classic finite-element method plays a minor role in seismology as its high-order sister, the spectral-element method, is more efficient. 
	\end{itemize}
	}
\end{frame}

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

\begin{frame}{Comprehension questions}
	\footnotesize{	
	\begin{enumerate}
\item[1] In which community was the finite-element method primarily developed? Give some typical problems. 
\item[2] What are weak and strong forms of partial differential equations? Give examples. 
\item[3] Discuss the pros and cons of the finite-element method vs. low-order finite-difference methods. 
\item[4] Present and discuss problem classes that can be handled well with the finite-volume method, compare to problems better handled with other methods. 
\item[5] Compare the spatial discretization strategies of finite-element and finite-difference methods. 
\item[6] Describe the concept of shape functions.
\item[7] Discuss qualitatively (use sketches) the use of basis functions. Compare with the pseudospectral method. 
\item[8] Is the finite-element method a global or a local scheme?
\item[9] Why does the finite-element method require the solution of a (possible huge) system of linear equations? What is the consequence for parallel computing?
\item[10] Why is the classic linear finite-element method not so much used for seismological research today?
\end{enumerate}
}
\end{frame}		

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

\begin{frame}{Theoretical questions}
	\footnotesize{	
\begin{enumerate}
\item[11] The scalar advection equation is simply 
\begin{equation*}
 \partial_t q(x,t) + c(x) \partial_x q(x,t) = 0
  \nonumber
\end{equation*}
where $q(x,t)$ is the scalar quantity to be advected and $c(x)$ is the advection velocity. Write down the weak form of this equation and perform integration by parts. 
What happens to the anti-derivative? Does it cancel out at the boundaries like in the 1D elastic wave equation? Note: This is the point of departure for the discontinuous Galerkin method).  
\item[12]  Are the linear basis functions 
\begin{equation*}
	\varphi_i (x) = 
	\begin{cases} \frac{x - x_{i-1}}{x_i - x_{i-1}} \;  \text{for} \;  x_{i-1} < x \leq x_i \\ 
	\frac{x_{i+1} - x}{x_{i+1} - x_i} \; \text{for} \; x_i < x < x_{i+1} \\
	\phantom{-} 0 \; \; \; \; \; \; \; \; \; \; \text{elsewhere}
	\end{cases}
	\nonumber
\end{equation*}
orthogonal?
\item[13] Calculate all entries of the stiffness matrix $ A_{i j} = \int_D \mu \partial_x \; \varphi_i \;  \partial_x \; \varphi_j$ for a static elastic problem with $\mu =70 GPa$ and $h = 1m$ for a problem with $n=5$ degrees of freedom. 
\end{enumerate}
}
\end{frame}	

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

\begin{frame}{Theoretical questions}
	\footnotesize{	
\begin{enumerate}
\item[14] A finite-element system has the following parameters: Element sizes $h=[1, 3, 0.5, 2, 4]$, density $\rho = [2, 3, 2, 3, 2] kg/m^3$. Calculate the entries of the mass matrix given by   $M_{ij} = \int_D \rho \; \varphi_i \; \varphi_j \; dx$ using linear basis functions.
\item[15] H-adaptivity. For the  simulation with varying velocities and element size, calculate the time step required for $\epsilon=0.5$ in each of the subdomains. Discuss the result.  
\item[16] Follow the approach of the derivation of shape functions and derive the cubic case in 1D: $u(x)=c_1 + c_2 \xi + c_3 \xi^2 + c_4 \xi^3$. What are key differences to quadratic and linear cases?   
\item[17] Derive the quadratic shape functions $N(\xi,\eta)$ for 2D triangles with the following node points: 
\begin{equation*}
\begin{split}
&P_1(0,0), P_2(1,0), P_3(0,1),\\
&P_4(1/2,0), P_5(1/2, 1/2), P_6(0, 1/2)
\end{split}
\nonumber
\end{equation*}
\item[18] Derive the quadratic shape functions $N(\xi,\eta)$ for 2D rectangles with the following node points: 
\begin{equation*}
\begin{split}
&P_1(0,0), P_2(1/2,0), P_3(1,0), P_4(1,1/2), \\
&P_5(1, 1), P_6(1/2, 1), P_7(0,1), P_8(0,1/2)
\end{split}
\nonumber
\end{equation*}
\end{enumerate}
}
\end{frame}		

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

\begin{frame}{Programming Exercise}
	\footnotesize{	
\begin{enumerate}
\item[19] Static elastic case. Extend the formulation given to arbitrary element sizes (calculation of stiffness matrix and source vector). Make examples and compare with the regular grid finite-difference solution (relaxation method). 
\item[20] 1D elastic case. Determine numerically the stability limit and compare with the finite-difference solution. 
\item[21] Initialize a strongly heterogeneous velocity model with spatially varying element size. Try to match the results with a regular grid finite-difference implementation of the same model. Discuss the two approaches in terms of time step, run time, memory usage. 
\item[22] Plot the high-order 2D shape functions derived in the problems above.   
\end{enumerate}
}
\end{frame}	

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

\end{document}
