% sample.tex
\documentclass[pdf,slideColor,nocolorBG]{prosper}
\hypersetup{colorlinks=true,anchorcolor=red,filecolor=red}

\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{amsbsy}
\usepackage{relsize}
\usepackage{psfrag}

\newcommand{\fenics}{\scalebox{1.3}{\textbf{\textsf{FE}}}\scalebox{1.6}{\textbf{\textsf{ni}}}\scalebox{1.3}{\textbf{\textsf{CS}}}}
\newcommand{\fenicssmall}{\scalebox{1.0}{\textbf{\textsf{FE}}}\scalebox{1.2}{\textbf{\textsf{ni}}}\scalebox{1.0}{\textbf{\textsf{CS}}}}
\newcommand{\FENICS}{\textbf{\textsf{\Large{FE}\Huge{ni}\Large{CS}}}}
\newcommand{\dolfin}{\scalebox{1.3}{\textbf{\textsf{DOLFIN}}}}
\newcommand{\dolfinsmall}{\scalebox{1.0}{\textbf{\textsf{DOLFIN}}}}
\newcommand{\fiat}{\scalebox{1.3}{\textbf{\textsf{FIAT}}}}
\newcommand{\ko}{\scalebox{1.3}{\textbf{\textsf{Ko}}}}
\newcommand{\puffin}{\scalebox{1.3}{\textbf{\textsf{Puffin}}}}
\newcommand{\puffinsmall}{\scalebox{1.0}{\textbf{\textsf{Puffin}}}}
\newcommand{\ferari}{\scalebox{1.3}{\textbf{\textsf{FErari}}}}
\newcommand{\ffc}{\scalebox{1.3}{\textbf{\textsf{FFC}}}}
\newcommand{\ffcsmall}{\scalebox{1.0}{\textbf{\textsf{FFC}}}}
\newcommand{\cgq}{$\mathrm{cG}(q)$}
\newcommand{\dgq}{$\mathrm{dG}(q)$}
\newcommand{\mcgq}{$\mathrm{mcG}(q)$}
\newcommand{\mdgq}{$\mathrm{mdG}(q)$}
\newcommand{\back}{\scalebox{1.3}{$\hookleftarrow$}}
\newcommand{\comp}[1]{{\scriptstyle [#1]}}
\newcommand{\ind}[1]{{\scriptstyle (\cdot)}}

\ptsize{8}

\def\movieroot{file:///root/data/txt/tex/CAM2004/movies/}

\hyperbaseurl{}


\title{\black Solving an elasto-plastic model using DOLFIN}
\subtitle{TTI 2005}
\author{Johan Jansson}
\institution{Chalmers University of Technology}
\email{johanjan@math.chalmers.se}



\begin{document}

\maketitle

\begin{slide}{Overview}

\bigskip

\begin{itemize}

\item
        Motivation

\item
        Previous work

\item
        Simple derivation of model

\item
        Plasticity

\item
        Discretization

\item
        \fenicssmall{}        

\item
        Implementation using FFC/DOLFIN

\item
        Performance results

\item
        Future work

\end{itemize}

\end{slide}


\begin{slide}{Motivation}

Computational elasticity useful in many fields:

\begin{itemize}
\item
Computer animation

\item
Computer games

\item
CAD

\item
...
\end{itemize}

BUT, established models are crude, ad hoc (outside of CAD).

\end{slide}

\begin{slide}{Motivation - examples}

State of the art computer games use rigid body motion with joints
(Half Life 2).

\bigskip

Motion pictures primarily use animation by hand, some cases of
mass-spring simulation (hair).

\bigskip

Why don't these applications use more advanced/general models?

\bigskip

Traditional elasticity models are difficult to understand $\Rightarrow$
difficult to apply, use effectively.

\bigskip

Attempt to find a simple model, attempt to automatize discretization
of model.

\end{slide}

\begin{slide}{Previous work - mass-spring model}

Point-masses connected by ``springs''.

\bigskip

The springs represent forces between the masses - Hooke elasticity ($F
= -k (e(x))$), damping ($F = -b \dot{e}(x)$) (e - elongation of spring).

\bigskip

Newton's second law: ($F = m \ddot{x}$).

\bigskip

Collision/contact handled by each mass having a radius - add a spring when
two masses intersect.

\begin{center}
\href{run:movies/ko-movie-collision.mpg}{\includegraphics[height=30mm]{movies/ko-movie-collision.eps}}
\end{center}

%\begin{center}
%\href{\movieroot/ko-movie-collision.mpg}{\includegraphics[height=30mm]{movies/ko-movie-collision.eps}}
%\end{center}

\end{slide}

\begin{slide}{Previous work - simulations}

\begin{center}
\href{run:movies/susp-physical.mpg}{\includegraphics[height=30mm]{movies/susp-physical.eps}}
\end{center}

\begin{center}
\href{run:movies/susp-geometric.mpg}{\includegraphics[height=30mm]{movies/susp-geometric.eps}}
\end{center}

\bigskip

Can we find an analogous PDE-model?

\end{slide}

\begin{slide}{Simple derivation of model}

Classical linear elasticity:

\begin{equation*}
\label{classicalelast}
\begin{split}
u &= x - X,\\
\dot{u}-v&=0\quad\mbox{in } \Omega^0,\\
\dot{v} - \nabla \cdot \sigma& =f \quad\mbox{in } \Omega^0,\\
\sigma &= E\epsilon (u) = E(\nabla u^\top +\nabla u)\\
E\epsilon &= \lambda\sum_k\epsilon_{kk} \delta_{ij} + 2\mu \epsilon,\\
v(0,\cdot ) &=v^0,\quad u(0,\cdot )=u^0\quad\mbox{in } \Omega^0.
\end{split}
\end{equation*}

\large

Only works for small displacements. Computations carried out on fixed
geometry $\Omega^0$. Why not use the deformed geometry $\Omega(t)$?

\end{slide}

\begin{slide}{Simple derivation of model}

$u$ depends on $\Omega^0$ (X is defined on $\Omega^0$). Let's compute in
$v$ instead.

\bigskip

Differentiate $\sigma$ wrt. $t$:

\begin{equation*}
\label{elasteq5}
\begin{split}
\sigma &= E\epsilon (u) = E(\nabla u^\top +\nabla u) \Rightarrow\\
\dot{\sigma} &= E\epsilon (v) = E(\nabla v^\top +\nabla v)
\end{split}
\end{equation*}

Rest of the model remains the same.

\end{slide}

\begin{slide}{The elastic model}

Now formulate the model in the deformed geometry $\Omega(t)$:

\begin{equation*}
\label{classicalelast}
\begin{split}
\dot{u} - v &= 0 \quad\mbox{in } \Omega(t),\\
\dot{v} - \nabla \cdot \sigma &= f \quad\mbox{in } \Omega(t),\\
\dot{\sigma} &= E\epsilon (v) = E(\nabla v^\top +\nabla v)\\
v(0,\cdot ) &=v^0,\quad u(0,\cdot )=u^0\quad\mbox{in } \Omega^0.
\end{split}
\end{equation*}

The model is a piecewise linear elastic model. Given some geometry
$\Omega_i$ we compute using the linear model (small displacements) for
a small time step and produce the geometry $\Omega_{i+1}$. The process
is then repeated.

\end{slide}

\begin{slide}{Examples}

\begin{center}
\href{run:movies/updelast-bar03-model1.mpg}{Elastic bar (Updated Lagrange)}
\end{center}

\begin{center}
\href{run:movies/massspring-bar03.mpg}{Elastic bar (Mass-spring)}
\end{center}

\begin{center}
\href{run:movies/ko-movie-elasticity-classic.mpg}{Elastic cube (Classical elasticity)}
\end{center}

\end{slide}

\begin{slide}{Viscosity}

\begin{equation*}
\begin{split}
\label{viscosity}
\dot{v} - \nabla \cdot \sigma - \nu \epsilon(v) &= f \quad\mbox{in } \Omega(t)
\end{split}
\end{equation*}

We add a simple viscous term to model viscosity in materials.

\end{slide}

\begin{slide}{Plasticity}

\begin{equation*}
\label{elastplast}
\begin{split}
\dot{\sigma} &= E(\epsilon (v) - \frac{1}{\nu_p}(\sigma -\pi\sigma) )\quad\mbox{in } \Omega (t),\\
\pi\sigma &= \frac{\sigma}{\Vert \sigma \Vert}, \Vert \sigma \Vert > Y_s\\
\pi\sigma &= \sigma, \Vert \sigma \Vert \leq Y_s
\end{split}
\end{equation*}

Visco-plastic model. $\pi \sigma$ is the projection on to the set of
admissible stresses. $Y_s$ is the yield stress of the material.

\end{slide}

\begin{slide}{Examples (Plasticity)}

\begin{center}
\href{run:movies/updelast-elast-diamond01.mpg}{Elastic diamond}
\end{center}

\begin{center}
\href{run:movies/updelast-plast-diamond01.mpg}{Plastic diamond}
\end{center}

\begin{center}
\href{run:movies/updelast-plast-bar01.mpg}{Plastic bar (no gravity)}
\end{center}

\end{slide}

\begin{slide}{Discretization of model}

Finite element discretization (elastic model):

\begin{equation*}
\label{classicalelast}
\begin{split}
\dot{v} - \nabla \cdot \sigma &= f \Rightarrow\\
\int_{\Omega} \dot{v} q\ dx - \int_{\Omega} \nabla \cdot \sigma q\ dx &= \int_{\Omega} f q\ dx
\end{split}
\end{equation*}

Ansatz:

\begin{equation*}
\label{classicalelast}
\begin{split}
v_h &= \sum_{j} \xi^v_j \phi_j\\
\sigma_h &= \sum_{j} \xi^{\sigma}_j \psi_j
\end{split}
\end{equation*}

$\phi$ piecewise linear. $\psi$ piecewise constant.

\end{slide}

\begin{slide}{Discretization of model}

Insert ansatz and choose test space as piecewise linear:

\begin{equation*}
\label{classicalelast}
\begin{split}
\int_{\Omega} \sum_{j} \dot{\xi}^v_j \phi_j \phi_i\ dx - \int_{\Omega} \sum_{j} \xi^{\sigma}_j \psi_j \nabla \phi_i\ dx &= \int_{\Omega} f \phi_i\ dx, \forall i \Rightarrow\\
M \xi^v - A \xi^{\sigma} &= b
\end{split}
\end{equation*}
\end{slide}

\begin{slide}{Discretization of model}

\begin{equation*}
\label{classicalelast}
\begin{split}
\dot{\sigma}^{K_i}_h &= E(\nabla v^\top +\nabla v)^{K_i}, \forall K_i \in T_\Omega(t)\\
\end{split}
\end{equation*}

\end{slide}


\begin{slide}{General examples}

\begin{center}
\href{run:movies/cow01.mpg}{Elastic cow}
\end{center}

\begin{center}
\href{run:movies/cow05.mpg}{Visco-elastic cow}
\end{center}

\begin{center}
\href{run:movies/cow06.mpg}{Plastic cow}
\end{center}


\end{slide}


\begin{slide}{The \fenicssmall{} project}

A free software project for the Automation of Computational
Mathematical Modeling (ACMM), developed at
\begin{itemize}
\item
  The Toyota Technological Institute at Chicago
\item
  The University of Chicago
\item
  Chalmers University of Technology (Göteborg, Sweden)
\item
  Argonne National Laboratory (Chicago)
\end{itemize}

More information at \texttt{http://www.fenics.org/}

\begin{center}
  \includegraphics[width=3cm]{eps/mesh-2d.eps}
  \includegraphics[width=3cm]{eps/mesh-3d.eps}
  \includegraphics[width=3cm]{eps/convdiff.eps}
\end{center}

\end{slide}
%------------------------------------------------------------------------------
\begin{slide}{Projects}

\begin{itemize}
\item
  \dolfin{}, the C++ interface of FEniCS
  \begin{itemize}
  \item
    Hoffman, Jansson, Logg, et al.
  \end{itemize}
\item
  \ferari{}, optimized form evaluation
  \begin{itemize}
  \item
    Kirby, Knepley, Scott
  \end{itemize}
\item
  \ffc{}, the FEniCS Form Compiler
  \begin{itemize}
  \item
    Logg
  \end{itemize}
\item
  \fiat{}, automatic generation of finite elements
  \begin{itemize}
  \item
    Kirby, Knepley
  \end{itemize}
\item
  \ko{}, simulation of mechanical systems
  \begin{itemize}
  \item
    Jansson
  \end{itemize}
\item
  \puffin{}, light-weight version for Octave/MATLAB
  \begin{itemize}
  \item
    Hoffman, Logg
  \end{itemize}
\end{itemize}

\end{slide}
%------------------------------------------------------------------------------

%------------------------------------------------------------------------------
\begin{slide}{}

\null\vspace{1.0cm}

\begin{center}

\LARGE

\textit{Solving a PDE with DOLFIN}

\normalsize

\end{center}

\end{slide}
%------------------------------------------------------------------------------
\begin{slide}{Implementing a solver}

\small\begin{verbatim}
void PoissonSolver::solve()
{
  Galerkin     fem;
  Matrix       A;
  Vector       x, b;
  Function     u(mesh, x);
  Function     f(mesh, "source");
  Poisson      poisson(f);
  KrylovSolver solver;
  File         file("poisson.m");
 
  fem.assemble(poisson, mesh, A, b);
  solver.solve(A, x, b);
 
  u.rename("u", "temperature");
  file << u;
}
\end{verbatim}

\end{slide}
%------------------------------------------------------------------------------
\begin{slide}{Automatic assembling}

\small\begin{verbatim}
class Poisson : public PDE {
  ...
  real lhs(const ShapeFunction& u, const ShapeFunction& v)
  {
    return (grad(u),grad(v)) * dK;
  }
   
  real rhs(const ShapeFunction& v)
  {
    return f*v * dK;
  }
  ...
};
\end{verbatim}

\end{slide}
%------------------------------------------------------------------------------
\begin{slide}{Automatic assembling}

\small\begin{verbatim}
class ConvDiff : public PDE {
  ...     
  real lhs(const ShapeFunction& u, const ShapeFunction& v)
  {
    return (u*v + k*((b,grad(u))*v + a*(grad(u),grad(v))))*dK;
  }
     
  real rhs(const ShapeFunction& v)
  {
    return (up*v + k*f*v) * dK;
  }
  ...
};
\end{verbatim}

\end{slide}
%------------------------------------------------------------------------------
\begin{slide}{Abstract assembly algorithm}

\begin{eqnarray*}
	(A_h)_{ij} &=& a(\varphi_j, \hat{\varphi}_i) = \int_{\Omega} A(\varphi_j) \hat{\varphi}_i \ dx \\
			 	  &=& \sum_{K\in\mathcal{T}} \int_K A(\varphi_j) \hat{\varphi}_i \ dx
              = \sum_{K\in\mathcal{T}} a(\varphi_j, \hat{\varphi}_i)_K.
\end{eqnarray*}
Iterate over all elements $K$ and for each element $K$
compute the contributions to all $(A_h)_{ij}$, for which $\varphi_j$ and $\hat{\varphi}_i$
are supported within $K$.

\bigskip

Assemble on reference element and map values to real element.

\end{slide}
%-----------------------------------------------------------------------------
\begin{slide}{The elasto-plastic model in DOLFIN}

\small\begin{verbatim}
class ElasticityUpdated : public PDE {

  ...

  real rhs(ShapeFunction::Vector& v)
  {
    ... Compute sigma ...

     return (k * (f(0) * v(0) + f(1) * v(1) + f(2) * v(2)) -
              k * (sigma(0, 0) * v(0).ddx() +
                   sigma(0, 1) * v(0).ddy() +
                   sigma(0, 2) * v(0).ddz() +
                   sigma(1, 0) * v(1).ddx() +
                   sigma(1, 1) * v(1).ddy() +
                   sigma(1, 2) * v(1).ddz() +
                   sigma(2, 0) * v(2).ddx() +
                   sigma(2, 1) * v(2).ddy() +
                   sigma(2, 2) * v(2).ddz())) * dx;
  }
};
\end{verbatim}

\end{slide}
%------------------------------------------------------------------------------
\begin{slide}{}

\null\vspace{1.0cm}

\begin{center}

\LARGE

\textit{Solving a PDE with FFC/DOLFIN}

\normalsize

\end{center}

\end{slide}
%------------------------------------------------------------------------------
\begin{slide}{FFC: the FEniCS Form Compiler}

\begin{itemize}
\item
  Automates a key step in the implementation of finite element methods
  for partial differential equations
\end{itemize}

\begin{itemize}
\item
  Input: a variational form and a finite element
\item
  Output: optimal C/C++
\end{itemize}

\vspace{1cm}

\begin{center}
  \psfrag{FFC}{\ \ffc{}}
  \psfrag{form}{$a(v, u) = \int_{\Omega} \nabla v(x) \cdot \nabla u(x) \, dx$}
  \includegraphics[width=11cm]{eps/ffc.eps}
\end{center}

\vspace{0.5cm}

\begin{verbatim}
  >> ffc [-l language] poisson.form
\end{verbatim}

\end{slide}
%------------------------------------------------------------------------------
\begin{slide}{Example: Classical Elasticity}
The form:

\begin{equation*}
\label{classicalelastform}
\begin{split}
a(u, v) &= \int_{\Omega} \sigma(u) \epsilon(v)\\
l(v) &= \int_{\Omega} f v
\end{split}
\end{equation*}

where (as before):

\begin{equation*}
\label{classicalelastform}
\begin{split}
\epsilon(u) &= \nabla u^\top +\nabla u\\
\sigma(u) &= E\epsilon (u)\\
E\epsilon &= \lambda\sum_k\epsilon_{kk} \delta_{ij} + 2\mu \epsilon,\\
\end{split}
\end{equation*}


\end{slide}
%------------------------------------------------------------------------------
\begin{slide}{Example: Classical Elasticity}

FFC representation ({\tt Elasticity.form}):

\small\begin{verbatim}
# The bilinear form for classical linear elasticity
# Compile this form with FFC: ffc Elasticity.form.

name = "Elasticity"
element = FiniteElement("Lagrange", "tetrahedron", 1, 3)

c1 = Constant()	# Lame coefficient
c2 = Constant()	# Lame coefficient
f = Function(element) # Source

v = BasisFunction(element)
u = BasisFunction(element)

a = (2.0 * c1 * u[i].dx(i) * v[j].dx(j) +
	c2 * (u[i].dx(j) + u[j].dx(i)) * (v[i].dx(j) + v[j].dx(i))) * dx
L = f[i] * v[i] * dx

\end{verbatim}

\end{slide}
%------------------------------------------------------------------------------
\begin{slide}{Example: Classical Elasticity}

FFC output ({\tt Elasticity.h}):

\small\begin{verbatim}
  BilinearForm(const real& c0, const real& c1) : ...

  bool interior(real* block) const
  {
    // Compute geometry tensors
    real G0_0_0_0_0 = det*c0*g00*g00;
    real G0_0_0_0_1 = det*c0*g00*g10;
    ...

    // Compute element tensor
    block[0] =
    3.333333333333329e-01*G0_0_0_0_0 + 3.333333333333329e-01*G0_0_0_0_1 +
    3.333333333333329e-01*G0_0_0_0_2 + 3.333333333333329e-01*G0_0_0_1_0 +
    3.333333333333329e-01*G0_0_0_1_1 + 3.333333333333329e-01*G0_0_0_1_2 +
    3.333333333333329e-01*G0_0_0_2_0 + 3.333333333333329e-01*G0_0_0_2_1 +
    3.333333333333329e-01*G0_0_0_2_2 + 1.666666666666664e-01*G1_0_0 +
    1.666666666666664e-01*G1_0_1 + 1.666666666666664e-01*G1_0_2 +
    1.666666666666664e-01*G1_1_0 + 1.666666666666664e-01*G1_1_1 +
    ...
\end{verbatim}

\end{slide}
%------------------------------------------------------------------------------
\begin{slide}{Benchmark}

\href{run:movies/ffc-linelast-beam01.mpg}{Beam example}

\bigskip

Task: assemble global matrix and load vector

\end{slide}
%------------------------------------------------------------------------------
\begin{slide}{Benchmark results}

\begin{itemize}
\item{Mesh 1}

24576 (25k) cells

\bigskip

old:\\
matrix - 19.2s\\
vector - 5.21s

\bigskip

ffc:\\
matrix - 0.770s\\
vector - 0.0900s

\bigskip

speedup:\\
matrix - 25\\
vector - 58
\end{itemize}

\end{slide}
%------------------------------------------------------------------------------
\begin{slide}{Benchmark results}

\begin{itemize}
\item{Mesh 2}

162000 (162k) cells

\bigskip

old:\\
matrix - 119s\\
vector - 33.2s\\

\bigskip

ffc:\\
matrix - 9.12s\\
vector - 0.560s

\bigskip

speedup:\\
matrix - 13\\
vector - 60


\end{itemize}

\end{slide}
%------------------------------------------------------------------------------



\begin{slide}{Future work}

Solve elasto-plastic model using FFC

\bigskip

Multi-adaptive time stepping

\bigskip

Space adaptivity

\bigskip

Interface to fluid mechanics (Navier-Stokes).

\end{slide}



\end{document}


