% 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{Introduction}

% ---------------------------------------------------------
\subsection{Motivation}

	\begin{frame}{Motivation}
		\begin{minipage}{0.45\textwidth}
			\begin{figure}
				\includegraphics[width=0.95\textwidth]{\figpath/fig_fe_grenoble.png}
			\end{figure}
		\end{minipage}
    \mbox{ }
    \hfill		
		\begin{minipage}{0.5\textwidth}
			\begin{block}{}
				\begin{itemize}
					\item Static and time-dependent PDEs
					\item Implicit approach
    				\item Irregular grids     
				\end{itemize}	
			\end{block}	
		\end{minipage}
	\end{frame}

% ---------------------------------------------------------
\subsection{History}

\begin{frame}{History}
	\footnotesize
	\begin{itemize}[<+->]
		\item As the free-surface boundary condition is implicitly fulfilled study of surface wave propagation (Lysmer 1972, Schlue 1979). 
		\item Seismic scattering problems were simulated with the method in the dissertation by Day (1977). 
		\item In addition to the physical propagation modes, parasitic modes for high-order implementations were found Kelly (1990).
		\item Li (1994) presented parallel implementations on the legendary CM-2 massively parallel supercomputer.
		\item Finite-element principles are also the basis for the so-called direct solution method that was introduced by R. Geller and co-workers.
		\item Problems in seismic shaking hazards and engineering seismology were conducted by the group of J. Bielak and co-workers.
		\item The methods were later extended to the problem of full waveform inversion (Askan 2008). 
		\item Hybrid methods that make use of advantages of both finite-difference and finite-element methods were presented by Moczo (2010). 
	\end{itemize}
\end{frame}

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

\begin{frame}{Finite Elements in a Nutshell}
		\begin{minipage}{0.45\textwidth}
			\begin{figure}
				\includegraphics[width=0.5\textwidth]{\figpath/fig_fe_nutshell.png}
			\end{figure}
		\end{minipage}
    \mbox{ }
    \hfill		
		\begin{minipage}{0.5\textwidth}
			1D elastic wave equation with space-dependent density $\rho$,  shear modulus $\mu$, and forcing term $f(x,t$) 
\begin{equation*}
	\rho \partial_t^2 u = \partial_x \; \mu \; \partial_x u + f \ . 
\end{equation*}
We seek to find solutions to the displacement field $u(x,t)$, thus, we replace it by a finite sum over (here at first linear) basis functions $\varphi_i$. Our unknowns are the coefficients of the basis functions  $\varphi_i$ 
\begin{equation*}
	u(x) \approx \overline{u} (x) = \sum_{i = 1}^N u_i(t) \varphi_i(x) \ .
\end{equation*}
		\end{minipage}	
\end{frame}

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

\begin{frame}{Finite Elements in a Nutshell}
	Furthermore, we formulate a so-called {\it weak} form of the wave equation, multiplying the original strong form by a test function $\varphi_j$ of the same basis, followed by an integration over the entire physical domain. This leads to a linear system of equations of the form  
\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*}
where we seek to find the approximate displacement field $\overline{u}$. Given appropriate initial conditions, the solution at the next time step $\textbf{u} (t+dt)$ can be found by the following matrix-vector equation
\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*}
where \textbf{M} and \textbf{K} are the mass matrix and stiffness matrix, respectively. 
\end{frame}

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

\begin{frame}{Finite Elements in a Nutshell}
	\begin{itemize}
		\item Global matrices in the sense that if a physical domain is discretized with $N$ elements, then the matrix sizes is $N \times N$.
		\item One of the matrices has to be inverted.
		\item Mass matrix  \textbf{M} consists of elements of the form $\int_D \rho \varphi_i  \varphi_j dx$ and the stiffness matrix \textbf{K} is built up with elements of the form $\int_D \mu \nabla \varphi_i \nabla  \varphi_j dx$.
		\item These integrals can be computed in an elegant way for each element by mapping the physical space to a local reference space.
	\end{itemize}
\end{frame}


% ---------------------------------------------------------
\section{Static Elasticity}


\begin{frame}{Discretization}
Departing from the 1D elastic wave equation
\begin{equation*}
		\rho \partial_t^2 u(x,t) = \partial_x \mu(x) \partial_x u(x,t) + f(x,t) 
\end{equation*}
we assume the following:
\vspace{0.5cm}
\begin{itemize}
	\item[] Independent in time: $\partial_t^2 u(x,t) = 0$
	\item[] Elastic properties of our 1D medium are independent of space: $\mu(x) = \text{const.}$
\end{itemize}
\vspace{0.5cm}
that leads to the equation
\begin{equation*}
		- \mu \partial_x^2 u \ = f \ . 
\end{equation*}		
\end{frame}

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


	\begin{frame}{Illustration}
		\begin{minipage}{0.45\textwidth}
			\begin{figure}
				\includegraphics[width=0.95\textwidth]{\figpath/fig_fe_static.png}
			\end{figure}
		\end{minipage}
    \mbox{ }
    \hfill		
		\begin{minipage}{0.5\textwidth}
			\begin{block}{}
Static elasticity. A string with
homogeneous properties (density and
shear modulus) is pulled with a certain
force. The Poisson equation determines
the displacement of the string given appropriate
boundary conditions. Don't
overdo this experiment, in particular if
you have old strings.	
			\end{block}	
		\end{minipage}
	\end{frame}

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

\begin{frame}{Weak form}
Transform {\it strong} form into {\it weak} form by multiplying the equation with an arbitrary space-dependent test function that we denote as $v \rightarrow v(x)$. Then we integrate the equation on both sides over the entire physical domain $D$ with  $x \in D$
\begin{equation*}
		- \int_D \mu \partial_x^2 \; u \; v \; dx = \int_D \; f \; v \; dx  \ .
\end{equation*}
Performing an integration by parts of the left side: 
\begin{equation*}
- \int_D \mu \partial_x^2 \; u \; v \; dx = \ \int_D \mu \; \partial_x \; u \;  \partial_x \; v \; dx - \left[ \mu \partial_x u \; v  \right]_{x_{min}}^{x_{max}} \ . 
\end{equation*}
where the last term is an anti-derivative. 
\end{frame}

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

\begin{frame}{Free surface}
Free-surface condition $\Longrightarrow$ No stress at the boundaries. \\
As the anti-derivative is evaluated at the domain boundaries this implies that this term vanishes.
\begin{equation*}
	\mu \int_D \partial_x \; u \;  \partial_x \; v \; dx = \int f \; v \; dx \ 
\end{equation*}
which is still a description in the continuous world. To enter the discrete world we replace our exact solution $u(x)$ by a $\overline{u}$, a sum over some basis functions $\varphi_i$ that we do not yet specify
\begin{equation*}
	u \approx \overline{u} (x) = \sum_{i = 1}^N u_i \varphi_i \ .
\end{equation*}
Replacing $u$ by $\overline{u}$, we obtain
\begin{equation*}
	\mu \int_D  \partial_x \; \overline{u} \;  \partial_x \; v \; dx = \int f \; v \; dx 
\end{equation*}
\end{frame}

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

\begin{frame}{Basis functions}
As a choice for our test function $v(x)$ we use the same set of basis functions. Thus $v(x) \rightarrow \varphi_j (x)$.  \\
\textbf{What is the simplemost choice for our basis functions $\varphi_i$?} Denoting $x_i, i=1,2,...,N$ as the boundaries of our elements we define our basis functions such that $\varphi_i = 1, x = x_i$ and zero elsewhere. Inside the elements our solution field is described by a linear function:
\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}
\end{equation*}
\end{frame}

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

\begin{frame}{Basis functions}
	\begin{figure}
		\includegraphics[width=0.5\textwidth]{\figpath/fig_fe_basis.png}
	\end{figure}
	\footnotesize{Linear basis functions for the finite-element method. A 1D domain is discretized with $n-1$ elements having $n=10$ element boundaries (open circles). The basis functions $\varphi_i=1$ at $x=x_{i}$. With this basis an arbitrary function can be exactly interpolated at the element boundary points $x_i$.}
\end{frame}

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

\begin{frame}{Galerkin Principle}
	We are ready to assemble our discretized version of the weak form by replacing the continuous displacement field by its approximation and applying the Galerkin principle. We obtain
\begin{equation*}
	\begin{split}
		\mu \int_D \partial_x \left( \sum_{i=1}^N u_i \varphi_i  \right) \partial_x \; \varphi_j \; dx & = \int f \; \varphi_j \; dx \\
		\sum_{i=1}^N u_i \int_D \mu \; \partial_x \; \varphi_i \; \partial_x \; \varphi_j \; dx & = \int f \; \varphi_j \; dx \\
	\end{split}
\end{equation*}
which is a system of N equations as we project the solution on the basis functions $\varphi_j$ with $j=1,...,N$. In the second equation we switched the sequence of integration and sum. 
\end{frame}

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

\begin{frame}{Matrix-Vector Notation}
	The discrete system thus obtained can be written using matrix-vector notation. Defining the solution vector $\textbf{u}$ as 
\begin{equation*}
\textbf{u} = \left( \begin{array}{c} u_1 \\ u_2 \\ \vdots \\ u_N  \end{array} \right)  \;
\end{equation*}
the source vector $\textbf{f}$ as
\begin{equation*}
	 \textbf{f} = \left( \begin{array}{c} \int_D f \; \varphi_1 \; dx \\  \int_D f \; \varphi_2 \; dx \\ \vdots \\ \int_D f \; \varphi_N \; dx \\ \end{array}	 \right)
\end{equation*}
and the matrix containing the integral over the basis function derivatives as $\textbf{K}$
\begin{equation*}
 \textbf{K} \rightarrow K_{i j} = \mu \int_D \partial_x \; \varphi_i \;  \partial_x \; \varphi_j
\end{equation*}
\end{frame}

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

\begin{frame}{Solution}
System of equations can be written in component form as 
\begin{equation*}
	u_i \; K_{i j} = f_j
\end{equation*}
where we use the Einstein summation convention and in matrix-vector notation
\begin{equation*}
 \; \; \; \textbf{K}^T \; \textbf{u} = \textbf{f}
\end{equation*}
Note: Matrix $\textbf{K}$ is called the {\ it stiffness} matrix. \\
This system of equations has as many unknowns as equations. Provided that the matrix is positive definite we can determine its inverse:
\begin{equation*}
	\textbf{u} = (\textbf{K}^T)^{-1} \textbf{f}	
\end{equation*}
\end{frame}

% ---------------------------------------------------------
\subsection{Boundary Conditions}

\begin{frame}{Boundary Conditions}
In case you would like to invoke specific values at the boundaries the approximate solution becomes
\begin{equation*}
	\overline{u}\ = \ u_1 \varphi_1 + \sum_{i = 2}^{N-1} u_i \varphi_i   \ +   \ u_N \varphi_N \ 
\end{equation*}
where $u_1$ and $u_N$ are the boundary values. The weak form becomes
\begin{equation*}
	\begin{split}
	\sum_{i=2}^{N-1} \; \mu \; \int_D \partial_x \; \varphi_i \; \partial_x \; \varphi_j \; dx & = \int_D f \; \varphi_j \; dx \\
	& + u(x_{min}) \int_D \partial_x \; \varphi_1 \; \partial_x \; \varphi_j \; dx \\
	& + u(x_{max}) \int_D \partial_x \; \varphi_N \; \partial_x \; \varphi_j \; dx
	\end{split}
\end{equation*}
\end{frame}

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

\begin{frame}{Boundary Conditions}
		\begin{minipage}{0.45\textwidth}
			\begin{figure}
					\includegraphics[width=0.9\textwidth]{\figpath/fig_fe_bc.png}
			\end{figure}
		\end{minipage}
    \mbox{ }
    \hfill		
		\begin{minipage}{0.5\textwidth}
Graphical representation of the matrix-vector system with boundary conditions. The global system matrix has $N-2 \times N-2$ elements. The system {\it feels} the boundary conditions through a modified source terms. The red spot indicates the source location inside the medium.
		\end{minipage}
\end{frame}

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

\begin{frame}{Reference Element, Mapping}
Mapping the physical domain to a reference element. In our case we center the local coordinate system denoted as $\xi$ at point $x_i$ and obtain 
\begin{equation*}
	\begin{split}	
	\xi & = x - x_i \\
	h_i & = x_i - x_{i-1} \\
	\end{split}
\end{equation*}
where $h_i$ denotes the size of element $i$ defined in the interval $x\in [x_i, x_{i+1}]$. The local basis functions becomes 
\begin{equation*}
	\varphi_i (\xi) = 
	\begin{cases} \frac{\xi}{h} +1 \; \; \text{for} \; - h < \xi \leq 0 \\
	1 - \frac{\xi}{h} \; \;  \text{for} \;  0 < \xi < h \\
	\phantom{-} 0 \; \; \; \; \; \; \; \; \; \; \text{elsewhere}
	\end{cases}
\end{equation*}
and their derivatives
\begin{equation*}
	\partial_\xi \; \varphi_i (\xi) = 
	\begin{cases} \phantom{-} \frac{1}{h} \; \; \text{for} \; - h < \xi \leq 0 \\
	- \frac{1}{h} \; \;  \text{for} \;  0 < \xi < h \\
	\phantom{-} 0 \; \; \; \; \; \; \;  \text{elsewhere}
	\end{cases}
\end{equation*}
\end{frame}

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

\begin{frame}{Reference Element, Mapping}
		\begin{minipage}{0.45\textwidth}
			\begin{figure}
				\includegraphics[width=0.6\textwidth]{\figpath/fig_fe_basis_der.png}
			\end{figure}
		\end{minipage}
    \mbox{ }
    \hfill		
		\begin{minipage}{0.5\textwidth}
Basis functions and their derivatives. {\bf Top:} The basis function $\varphi_i$ (thick solid line) is shown along with the neighboring functions $\varphi_{i\pm1}$ (thin dotted lines). {\bf Bottom:} The same for their derivatives with respect to the space coordinate $ \xi$.
		\end{minipage}
\end{frame}

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

\begin{frame}{Stiffness Matrix}
	We can now proceed with calculating the elements of the stiffness matrix $K$ defined as
\begin{equation*}
	K_{i j} = \mu \int_D \partial_x  \varphi_i \; \partial_x \; \varphi_j \; dx \ 
\end{equation*}
with the corresponding expression in local coordinates $\xi$  
\begin{equation*}
	K_{i j} = \mu \int_{D_\xi} \partial_\xi  \varphi_i \; \partial_\xi \varphi_j \; d\xi \ . 
\end{equation*}
\end{frame}

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

\begin{frame}{Stiffness Matrix}
	 Let us calculate some of the elements of matrix $K_{ij}$ starting with the diagonal elements. For example, for $K_{11}$ we obtain
\begin{equation*}
	\begin{split}
		K_{11} & = \mu \int_D \partial_x \; \varphi_1 \; \partial_x \; \varphi_1 \; dx \\
		&= \mu \int_{0}^{h} \frac{-1}{h} \frac{-1}{h} d\xi \ =  \frac{\mu}{h^2} \int_0^{h} d\xi = \frac{\mu}{h} 
		 \end{split}
\end{equation*}
For diagonal element $A_{22}$ the derivatives overlap in element $1$ and $2$ implying the integration has to be performed for the interval $\xi \in [-h, h]$.  
\begin{equation*}
	\begin{split}	
		K_{22} & = \mu \int_D \partial_x \; \varphi_2 \; \partial_x \; \varphi_2 \; dx  \\
		& = \ \mu \int_{-h}^{0} \partial_\xi \; \varphi_2 \; \partial_\xi \; \varphi_2 \; d\xi + \mu \int_{0}^{h} \partial_\xi \; \varphi_2 \; \partial_\xi \; \varphi_2 \; d\xi \\
		& = \frac{\mu}{h^2}  \int_{-h}^0 d\xi + \frac{\mu}{h^2} \int_0^{h} = \frac{2\mu}{h} 
		\end{split}
\end{equation*}
\end{frame}

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

\begin{frame}{Stiffness Matrix}
	Equivalently, the off-diagonal terms overlap only in one element, for example
\begin{equation*}
\begin{split}
		K_{1 2} & = \mu \int_D \partial_x \; \varphi_1 \; \partial_x \; \varphi_2 \; dx \\
		&= \mu \int_{0}^{h} \partial_\xi \; \varphi_1 \; \partial_\xi \; \varphi_2 \; d\xi =  \mu \int_{0}^{h} \frac{-1}{h} \frac{1}{h} d\xi \\
		& = \frac{-\mu}{h^2} \int_0^{h} d\xi = \frac{-\mu}{h} \\
		K_{2 1} & = K_{1 2}\\
	\end{split}
\end{equation*}
\end{frame}


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

\begin{frame}{Stiffness Matrix}
Finally, the stiffness matrix for an elastic physical system with constant shear modulus $\mu$ and element size $h$ reads
\begin{equation*}
K_{i j} = \frac{\mu}{h}
	\begin{pmatrix}
1 & -1 &  &  & \\ 
-1 & 2 & -1 &  & \\ 
 &  & \ddots  &  & \\ 
 &  & -1 & 2 & -1\\ 
 &  &  & -1 & 1
	\end{pmatrix}
\end{equation*}
\textbf{Note:} Space-dependent terms in our linear system are proportional to the 3-point operator matrix for a 2nd finite-difference derivative	
\end{frame}

% ---------------------------------------------------------
\subsection{Simulation Example}

\begin{frame}{Simulation Example}
	\begin{minipage}{0.45\textwidth}
		\begin{table}
			\begin{tabular}[t]{ll}
Parameter & Value \\
\hline
$x_{max}$ & 1  \\
nx & 20 \\
$\mu$ & 1 \\
h & 0.0526 \\
u(0) & 0.15 \\
u(1) & 0.05
\end{tabular}
\end{table}				
	\end{minipage}
    \mbox{ }
    \hfill		
	\begin{minipage}{0.5\textwidth}
		The physical domain is defined in the interval $x\in  [0,1]$ and we apply a unit forcing at $x=0.75$ at one of the boundary element points.  	
	\end{minipage}
\end{frame}

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

\begin{frame}[fragile]{Simulation Example}
	\begin{minipage}{0.45\textwidth}
			\begin{lstlisting}
(...)
% Basic parameters
nx=20;   % Number of boundary points
u=zeros([1 nx]);  % solution vector
f=zeros([1 nx]);  % source vector
mu=1;  % constant shear modulus
% Element boundary points
x=linspace(0,1,nx); % x in [0, 1]
h=x(2)-x(1);   % constant element size

% Assemble stiffness matrix Kij
K=zeros(nx);

for i=2:nx-1,
   for j=2:nx-1,
      if i==j,
         K(i,j)=2*mu/h;
      elseif i==j+1
         K(i,j)=-mu/h;
      elseif i+1==j
         K(i,j)=-mu/h;
      else
         K(i,j)=0;
      end
   end
end
		\end{lstlisting}
	\end{minipage}
    \mbox{ }
    \hfill		
	\begin{minipage}{0.5\textwidth}
		\begin{lstlisting}
% source term is a spike at i=3*nx/4; 
f(3*nx/4)=1.;
% Boundary condition at x=0;
u(1)=0.15;   f(2) =u(1)/h;
% Boundary condition at x=1;
u(nx)=0.05; f(nx-1)=u(nx)/h;

% finite element solution
u(2:nx-1)=inv(K(2:nx-1,2:nx-1))*f(2:nx-1)';

(...)
		\end{lstlisting}
	\end{minipage}
\end{frame}

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

\begin{frame}{Simulation Example}
Starting with the Poisson equation $-\mu\partial^2_x u=f$, omitting space dependencies, we replace the l.h.s. with a finite-difference approximation and obtain
\begin{equation*}
-\mu \frac{u(x-h)-2u(x)+u(x+h)}{h^2} = f
\end{equation*}
and after rearranging
\begin{equation*}
u(x)=\frac{u(x-h)+u(x+h)}{2}-\frac{h^2}{2\mu} f \ . 
\end{equation*}
This equation can be used as an iterative procedure with an initial guess for the the unknown field $u$. 
\end{frame}

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

\begin{frame}[fragile]{Simulation Example}
With discretization $u_i=u(x_i)$ and iteration step $k$ this can be written as
\begin{equation*}
u_i^{k+1}=\frac{u^k_{i+1}-u^k_{i-1}}{2}-\frac{h^2}{2\mu} f_i
\end{equation*}
with initial guess $u_i^{k=1}=0$ - also called a {\it relaxation method}.
\begin{lstlisting}
(...)
% forcing
f(3*nx/4)=1/h;

for it = 1:nit, % iterate
% calculate the average of u (omit boundaries)
for i=2:nx-1, du(i)=u(i+1)+u(i-1); end
u=.5*( f*h^2 + du );
u(1)=0.15;  % Boundary at x=0;
u(nx)=.05;  % Boundary at x=1;
end
(...)
\end{lstlisting}
\end{frame}

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

\begin{frame}{Simulation Example}
	\begin{minipage}{0.45\textwidth}
		\begin{figure}
			\includegraphics[width=0.9\textwidth]{\figpath/fig_fe1d_static.png}
		\end{figure}			
	\end{minipage}
    \mbox{ }
    \hfill		
	\begin{minipage}{0.5\textwidth}
Simulation example comparing the finite-element solution (thick solid line) with a finite-difference based relaxation method (thin lines) that iteratively converges to the correct solution (see text for details). 500 iterations were employed for the relaxation method and the solution is shown after any 25 iterations.
	\end{minipage}
\end{frame}

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

%\begin{frame}
	%\begin{center}
		%\begin{block}{}
			%\centering \Huge The Finite-Element Method: 1D Elastic Wave Equation
		%\end{block}
	%\end{center}
%\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}



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

\end{document}