--- a/doc/doc.bib	Thu Dec 19 09:42:38 2013 +0100
+++ b/doc/doc.bib	Fri Dec 20 09:09:56 2013 +0100
@@ -57,6 +57,19 @@
  title = {\url{http://gedeone-gsoc.blogspot.co.uk/2013/06/update-4.html}},
 }
 
+@Unpublished{mixedpois,
+ title = {\url{http://fenicsproject.org/documentation/dolfin/1.2.0/python/demo/pde/mixed-poisson/python/documentation.html}},
+}
+
+@Unpublished{navierstokes,
+ title = {\url{http://fenicsproject.org/documentation/dolfin/1.2.0/python/demo/pde/navier-stokes/python/documentation.html}},
+}
+
+@Unpublished{hyperelasticity,
+ title = {\url{http://fenicsproject.org/documentation/dolfin/1.2.0/python/demo/pde/hyperelasticity/python/documentation.html}},
+}
+
+
 @Article{Formaggia_smart,
  author = {Formaggia, Luca},
  title = {Advanced Programming for Scientific Computing. Lecture title: Smart Pointers},

Binary file doc/doc.pdf has changed

--- a/doc/doc.tex	Thu Dec 19 09:42:38 2013 +0100
+++ b/doc/doc.tex	Fri Dec 20 09:09:56 2013 +0100
@@ -15,6 +15,8 @@
 \usepackage[english]{babel}
 \usepackage[T1]{fontenc}
 \usepackage[framemethod=TikZ]{mdframed}
+\usepackage{parcolumns}
+\usepackage[strict]{changepage}
 
 \definecolor{dkgreen}{rgb}{0,0.6,0}
 \definecolor{gray}{rgb}{0.5,0.5,0.5}
@@ -39,15 +41,48 @@
   tabsize=2
 }
 
-\author{Marco Vassallo}
-\title{ \textbf{Fem-fenics} \\ \bigskip \textit{General purpose Finite Element library \\ for GNU-Octave}\\ 
-\bigskip
-       \textsc{ Work in progress \\(help and remarks are welcome)}}
-
-
+\def\changemargin#1#2{\list{}{\rightmargin#2\leftmargin#1}\item[]}
+\let\endchangemargin=\endlist
 \begin{document}
 
-\maketitle
+\begin{titlepage}
+\begin{center}
+
+% Upper part of the page. The '~' is needed because \\
+% only works if a paragraph has started.
+
+\textsc{\LARGE Google Summer of Code 2013}\\[0.5cm]
+
+\textsc{\Large GNU-Octave}\\[2.5cm]
+
+%\includegraphics[width=1.\textwidth]{./octave-header.png}~\\[1.5cm]
+
+{ \huge \bfseries Fem-fenics \\[0.5cm] }
+{ \Large \bfseries Genaral Purpose Finite Element Library for GNU-Octave \\[2.4cm] }
+% Author and supervisor
+%\begin{minipage}{0.4\textwidth}
+%\begin{flushleft} \large
+%\emph{Author:}\\
+Marco \textsc{Vassallo}\\[0.5cm]
+%matr. 780787
+%\end{flushleft}
+%\end{minipage}
+%\begin{minipage}{0.4\textwidth}
+%\begin{flushright} \large
+%\emph{Supervisor:} \\
+%Dr.~Carlo \textsc{de Falco}
+%\end{flushright}
+%\end{minipage}
+
+\vfill
+
+% Bottom of the page
+{\large Version 0.0 \\[0.5cm]}
+{\large \today}
+
+\end{center}
+\end{titlepage}
+
 \tableofcontents
 
 \chapter{Introduction}
@@ -1036,6 +1071,11 @@
 
 
 \chapter{More Advanced Examples}\label{exem}
+In this chapter more examples are provided.
+At the beginning of each section, the problem is briefly presented and then
+the Octave script for the resolution of the problem using Fem-fenics is presented alongside the code 
+written in C++ and/or the Python.
+For each problem, we refer the reader to the complete desciption on the FEniCS website.
 \iffalse
 In the following examples we can see directly in action the classes and the functions presented in the 
 chapters before. A comparison with DOLFIN is given only for the first example, while more extensive case can
@@ -1052,150 +1092,170 @@
     Advection-Diffusion: we solve a time dependent problem using the lsode () command and save 
     the solution using the pkg flp.
 
-For each problem, we refer the reader to the complete desciption on FEniCS or bim web-page, 
+or bim web-page, 
 while here we highlight only the  implementation detail relevant for our pkg.
 \fi
-\section{Navier-Stokes equation with Chorin-Temam projection algorithm}
-Navier-Stokes: we learn how to deal with a vector-field problem and how we can save the solution using the 
-\verb$fem_save$ () function. We also use the msh pkg to generate a mesh using gmesh.
-\paragraph{TentativeVelocity.ufl}
-\begin{lstlisting}[language=Octave]
-# Copyright (C) 2010 Anders Logg
-# Define function spaces (P2-P1)
-V = VectorElement("CG", triangle, 2)
-Q = FiniteElement("CG", triangle, 1)
-
-# Define trial and test functions
-u = TrialFunction(V)
-v = TestFunction(V)
-
-# Define coefficients
-k  = Constant(triangle)
-u0 = Coefficient(V)
-f  = Coefficient(V)
-nu = 0.01
-
-# Define bilinear and linear forms
-eq = (1/k)*inner(u - u0, v)*dx + inner(grad(u0)*u0, v)*dx + \
-    nu*inner(grad(u), grad(v))*dx - inner(f, v)*dx
-a  = lhs(eq)
-L  = rhs(eq) 
- 
-\end{lstlisting}
-
-\paragraph{PressureUpdate.ufl}
-\begin{lstlisting}[language=Octave]
- # Copyright (C) 2010 Anders Logg
- # Define function spaces (P2-P1)
-V = VectorElement("CG", triangle, 2)
-Q = FiniteElement("CG", triangle, 1)
-
-# Define trial and test functions
-p = TrialFunction(Q)
-q = TestFunction(Q)
-
-# Define coefficients
-k  = Constant(triangle)
-u1 = Coefficient(V)
-
-# Define bilinear and linear forms
-a = inner(grad(p), grad(q))*dx
-L = -(1/k)*div(u1)*q*dx 
- 
-\end{lstlisting}
-
-\paragraph{VelocityUpdate.ufl}
-\begin{lstlisting}[language=Octave]
- # Copyright (C) 2010 Anders Logg
-# Define function spaces (P2-P1)
-V = VectorElement("CG", triangle, 2)
-Q = FiniteElement("CG", triangle, 1)
+\section{Mixed Formulation for the Poisson Equation}
+In this example the Poisson equation is solved with a 
+''mixed approach'': it is usedthe stable FE space obtained using Brezzi-Douglas-Marini 
+polynomial of order 1 and Dicontinuos element of order 0.
+\begin{align*}
+-\mathrm{div}\ ( \mathbf{\sigma} (x, y) ) ) &= f (x, y) & \quad \mbox{ in } \Omega \\
+\sigma (x, y) &= \nabla u (x, y) & \quad \mbox{ in } \Omega \\
+u(x, y) &= 0 & \quad \mbox{ on } \Gamma_D \\
+(\sigma (x, y) )  \cdot \mathbf{n} &= \sin (5x) & \quad \mbox{ on } \Gamma_N
+\end{align*}
+
+A complete description of the problem is avilable on the Fenics website \cite{mixedpois}.
+\begin{changemargin}{-1.5cm}{-1.5cm}
+$\phantom {u}$
+\begin{parcolumns}[colwidths={1=0.65\textwidth,2=0.65\textwidth}]{2}
+\colchunk{\begin{lstlisting}[caption=Fem-fenics, language=Octave, numbers=none]{Name}
+pkg load fem-fenics msh
+import_ufl_Problem ('MixedPoisson')
+
+# Create mesh
+x = y = linspace (0, 1, 33);
+mesh = Mesh(msh2m_structured_mesh (x, y, 1, 1:4));
+
+# File MixedPoisson.ufl
+#  BDM = FiniteElement("BDM", triangle, 1)
+#  DG  = FiniteElement("DG", triangle, 0)
+#  W = BDM * DG
+V = FunctionSpace('MixedPoisson', mesh);
+
+# Define trial and test function
+# File MixedPoisson.ufl
+#  (sigma, u) = TrialFunctions(W)
+#  (tau, v)   = TestFunctions(W)
+#  CG = FiniteElement("CG", triangle, 1)
+#  f = Coefficient(CG)
+f = Expression ('f', 
+      @(x,y) 10*exp(-((x - 0.5)^2 + (y - 0.5)^2) / 0.02));
+
+# Define variational form
+# File MixedPoisson.ufl
+#  a = (dot(sigma, tau) + div(tau)*u + div(sigma)*v)*dx
+#  L = - f*v*dx
+a = BilinearForm ('MixedPoisson', V, V);
+L = LinearForm ('MixedPoisson', V, f);
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+# Define essential boundary
+
+bc1 = DirichletBC (SubSpace (V, 1), @(x,y) [0; -sin(5.0*x)], 1);
+bc2 = DirichletBC (SubSpace (V, 1), @(x,y) [0;  sin(5.0*x)], 3);
+
+# Compute solution
+[A, b] = assemble_system (a, L, bc1, bc2);
+sol = A \ b;
+func = Function ('func', V, sol);
+
+sigma = Function ('sigma', func, 1);
+u = Function ('u', func, 2);
+
+# Plot solution
+plot (sigma);
+plot (u);
+
+
+#
+\end{lstlisting}}
+
+\colchunk{\begin{lstlisting}[caption=Python, language=Python, numbers=none]{Name}
+from dolfin import *
+
+
+# Create mesh
+mesh = UnitSquareMesh(32, 32)
+
+
+# Define function spaces and mixed (product) space
+BDM = FunctionSpace(mesh, "BDM", 1)
+DG = FunctionSpace(mesh, "DG", 0)
+W = BDM * DG
+
+
 
 # Define trial and test functions
-u = TrialFunction(V)
-v = TestFunction(V)
-
-# Define coefficients
-k  = Constant(triangle)
-u1 = Coefficient(V)
-p1 = Coefficient(Q)
-
-# Define bilinear and linear forms
-a = inner(u, v)*dx
-L = inner(u1, v)*dx - k*inner(grad(p1), v)*dx
-\end{lstlisting}
-
-\paragraph{NS.m}
-\begin{lstlisting}[language=Octave]
-pkg load fem-fenics msh
-import_ufl_Problem ("TentativeVelocity");
-import_ufl_Problem ("VelocityUpdate");
-import_ufl_Problem ("PressureUpdate");
- 
-# We can either load the mesh from the file as in Dolfin but 
-# we can also use the msh pkg to generate the L-shape domain
-L-shape-domain;
-mesh = Mesh (msho);
- 
-# Define function spaces (P2-P1).
-V = FunctionSpace ('VelocityUpdate', mesh);
-Q = FunctionSpace ('PressureUpdate', mesh);
- 
-# Set parameter values and define coefficients
-dt = 0.01;
-T = 3.;
-k = Constant ('k', dt);
-f = Constant ('f', [0; 0]);
-u0 = Expression ('u0', @(x,y) [0; 0]);
-
-# Define boundary conditions
-noslip = DirichletBC (V, @(x,y) [0; 0], [3, 4]);
-outflow = DirichletBC (Q, @(x,y) 0, 2);
-
-# Assemble matrices
-a1 = BilinearForm ('TentativeVelocity', V, V, k);
-a2 = BilinearForm ('PressureUpdate', Q, Q);
-a3 = BilinearForm ('VelocityUpdate', V, V);
-A1 = assemble (a1, noslip);
-A3 = assemble (a3, noslip);
-
-# Time-stepping
-t = dt; i = 0;
-while t < T
- 
-  # Update pressure boundary condition
-  inflow = DirichletBC (Q, @(x,y) sin(3.0*t), 1);
- 
-  # Compute tentative velocity step
-  L1 = LinearForm ('TentativeVelocity', V, k, u0, f);
-  b1 = assemble (L1, noslip);
-  utmp = A1 \ b1;
-  u1 = Function ('u1', V, utmp);
- 
-  # Pressure correction
-  L2 = LinearForm ('PressureUpdate', Q, u1, k);
-  [A2, b2] = assemble_system (a2, L2, inflow, outflow);
-  ptmp = A2 \ b2;
-  p1 = Function ('p1', Q, ptmp);
- 
-  # Velocity correction
-  L3 = LinearForm ('VelocityUpdate', V, k, u1, p1);
-  b3 = assemble (L3, noslip);
-  ut = A3 \ b3;
-  u1 = Function ('u0', V, ut);
-  
-  # Save to file
-  save (p1, sprintf ("p_%3.3d", ++i));
-  save (u1, sprintf ("u_%3.3d", i));
- 
-  # Move to next time step
-  u0 = u1;
-  t += dt
- 
-end
-\end{lstlisting}
-
-\paragraph{L-shape-domain.m}
+(sigma, u) = TrialFunctions(W)
+(tau, v) = TestFunctions(W)
+
+
+f = Expression
+    ("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)")
+
+# Define variational form
+
+a = (dot(sigma, tau) + div(tau)*u + div(sigma)*v)*dx
+L = - f*v*dx
+
+
+
+# Define function G such that G \cdot n = g
+class BoundarySource(Expression):
+    def __init__(self, mesh):
+        self.mesh = mesh
+    def eval_cell(self, values, x, ufc_cell):
+        cell = Cell(self.mesh, ufc_cell.index)
+        n = cell.normal(ufc_cell.local_facet)
+        g = sin(5*x[0])
+        values[0] = g*n[0]
+        values[1] = g*n[1]
+    def value_shape(self):
+        return (2,)
+
+G = BoundarySource(mesh)
+
+# Define essential boundary
+def boundary(x):
+    return x[1] < DOLFIN_EPS or x[1] > 1.0 - DOLFIN_EPS
+
+bc = DirichletBC(W.sub(0), G, boundary)
+
+# Compute solution
+w = Function(W)
+solve(a == L, w, bc)
+
+
+(sigma, u) = w.split()
+
+
+# Plot sigma and u
+plot(sigma)
+plot(u)
+interactive()
+
+# Copyright 2011, The FEniCS Project
+\end{lstlisting}}
+
+\colplacechunks
+\end{parcolumns}
+\end{changemargin} 
+
+\section{Incompressible Navier-Stokes equation}
+In this example the incompressible Navier-Stokes equation
+\begin{align*}
+\dfrac{\partial u}{\partial t} + (\mathbf u \cdot \mathrm{\nabla})  \mathbf u - \nu \Delta \mathbf u 
+  + \nabla p &= f & \quad \mbox{ in } \Omega \\
+\mathrm{\nabla} \cdot \mathbf u &= 0 & \quad \mbox{ in } \Omega \\
+\end{align*}
+are solved using the Chorin-Temam algorithm. The L-shaped domain $\Omega$ can be obtained using
+the msh pkg.
 \begin{lstlisting}[language=Octave]
 name = [tmpnam ".geo"];
 fid = fopen (name, "w");
@@ -1219,8 +1279,653 @@
 unlink (canonicalize_file_name (name));
 \end{lstlisting}
 
-\section{A penalization method to take into account obstacles in incompressible viscous flows}
-
+The flow is driven by an oscillating pressure $p_{in}(t) = \sin 3t$ at the inflow 
+while the pressure is kept constant $p_{out} = 0$ at the outflow.
+A complete description of the problem is avilable on the Fenics website \cite{navierstokes}.
+
+\begin{changemargin}{-1.5cm}{-1.5cm}
+$\phantom {u}$
+\begin{parcolumns}[colwidths={1=0.65\textwidth,2=0.65\textwidth}]{2}
+\colchunk{\begin{lstlisting}[caption=Fem-fenics, language=Octave, numbers=none]{Name}
+pkg load fem-fenics msh
+import_ufl_Problem ("TentativeVelocity");
+import_ufl_Problem ("VelocityUpdate");
+import_ufl_Problem ("PressureUpdate");
+
+# We can either load the mesh from the file as in Dolfin but 
+# we can also use the msh pkg to generate the L-shape domain
+# as showed above
+
+mesh = Mesh ('lshape.xml');
+
+# Define function spaces (P2-P1). From ufl file
+#  V = VectorElement("CG", triangle, 2)
+#  Q = FiniteElement("CG", triangle, 1)
+V = FunctionSpace ('VelocityUpdate', mesh);
+Q = FunctionSpace ('PressureUpdate', mesh);
+
+# Define trial and test functions. From ufl file
+#  u = TrialFunction(V)
+#  p = TrialFunction(Q)
+#  v = TestFunction(V)
+#  q = TestFunction(Q)
+
+# Set parameter values. From ufl file
+#  nu = 0.01
+dt = 0.01;
+T = 3.;
+
+
+
+
+# Define boundary conditions
+noslip = DirichletBC (V, @(x,y) [0; 0], [3, 4]);
+
+
+
+
+
+
+outflow = DirichletBC (Q, @(x,y) 0, 2);
+
+
+
+
+# Create functions
+u0 = Expression ('u0', @(x,y) [0; 0]);
+
+
+
+# Define coefficients
+k = Constant ('k', dt);
+f = Constant ('f', [0; 0]);
+
+# Tentative velocity step. From ufl file
+#  eq = (1/k)*inner(u - u0, v)*dx + inner(grad(u0)*u0, v)*dx \
+#       + nu*inner(grad(u), grad(v))*dx - inner(f, v)*dx
+a1 = BilinearForm ('TentativeVelocity', V, V, k);
+
+# Pressure update. From ufl file
+#  a = inner(grad(p), grad(q))*dx
+#  L = -(1/k)*div(u1)*q*dx
+a2 = BilinearForm ('PressureUpdate', Q, Q);
+
+# Velocity update
+#  a = inner(u, v)*dx
+#  L = inner(u1, v)*dx - k*inner(grad(p1), v)*dx
+a3 = BilinearForm ('VelocityUpdate', V, V);
+
+# Assemble matrices
+A1 = assemble (a1, noslip);
+
+A3 = assemble (a3, noslip);
+
+
+
+
+
+
+
+
+
+
+# Time-stepping
+t = dt; i = 0;
+while t < T
+
+  # Update pressure boundary condition
+  inflow = DirichletBC (Q, @(x,y) sin(3.0*t), 1);
+
+  # Compute tentative velocity step
+  "Computing tentative velocity"
+  L1 = LinearForm ('TentativeVelocity', V, k, u0, f);
+  b1 = assemble (L1, noslip);
+  utmp = A1 \ b1;
+  u1 = Function ('u1', V, utmp);
+
+  # Pressure correction
+  "Computing pressure correction"
+  L2 = LinearForm ('PressureUpdate', Q, u1, k);
+  [A2, b2] = assemble_system (a2, L2, inflow, outflow);
+  ptmp = A2 \ b2;
+  p1 = Function ('p1', Q, ptmp);
+
+  # Velocity correction
+  "Computing velocity correction"
+  L3 = LinearForm ('VelocityUpdate', V, k, u1, p1);
+  b3 = assemble (L3, noslip);
+  ut = A3 \ b3;
+  u1 = Function ('u0', V, ut);
+
+  # Plot solution
+  plot (p1);
+  plot (u1);
+
+  # Save to file
+  save (p1, sprintf ("p_%3.3d", ++i));
+  save (u1, sprintf ("u_%3.3d", i));
+
+  # Move to next time step
+  u0 = u1;
+  t += dt
+
+end
+
+
+
+#
+\end{lstlisting}}
+
+\colchunk{\begin{lstlisting}[caption=Python, language=Python, numbers=none]{Name}
+from dolfin import *
+
+
+
+
+# Load mesh from file
+
+
+
+
+
+mesh = Mesh("lshape.xml")
+
+# Define function spaces (P2-P1)
+
+
+V = VectorFunctionSpace(mesh, "CG", 2)
+Q = FunctionSpace(mesh, "CG", 1)
+
+# Define trial and test functions
+
+u = TrialFunction(V)
+p = TrialFunction(Q)
+v = TestFunction(V)
+q = TestFunction(Q)
+
+# Set parameter values
+dt = 0.01
+T = 3
+nu = 0.01
+
+# Define time-dependent pressure BC
+p_in = Expression("sin(3.0*t)", t=0.0)
+
+# Define boundary conditions
+noslip  = DirichletBC(V, (0, 0),
+           "on_boundary && \
+           (x[0] < DOLFIN_EPS | x[1] < DOLFIN_EPS | \
+           (x[0] > 0.5 - DOLFIN_EPS && x[1] > 0.5 - DOLFIN_EPS))")
+inflow  = DirichletBC(Q, p_in, "x[1] > 1.0 - DOLFIN_EPS")
+outflow = DirichletBC(Q, 0, "x[0] > 1.0 - DOLFIN_EPS")
+bcu = [noslip]
+bcp = [inflow, outflow]
+
+# Create functions
+u0 = Function(V)
+u1 = Function(V)
+p1 = Function(Q)
+
+# Define coefficients
+k = Constant(dt)
+f = Constant((0, 0))
+
+# Tentative velocity step
+F1 = (1/k)*inner(u - u0, v)*dx + inner(grad(u0)*u0, v)*dx \
+     + nu*inner(grad(u), grad(v))*dx - inner(f, v)*dx
+a1 = lhs(F1)
+L1 = rhs(F1)
+
+# Pressure update
+a2 = inner(grad(p), grad(q))*dx
+L2 = -(1/k)*div(u1)*q*dx
+
+
+# Velocity update
+a3 = inner(u, v)*dx
+L3 = inner(u1, v)*dx - k*inner(grad(p1), v)*dx
+
+
+
+# Assemble matrices
+A1 = assemble(a1)
+A2 = assemble(a2)
+A3 = assemble(a3)
+
+# Use amg preconditioner if available
+prec = "amg" if has_krylov_solver_preconditioner("amg") 
+             else "default"
+
+# Create files for storing solution
+ufile = File("results/velocity.pvd")
+pfile = File("results/pressure.pvd")
+
+# Time-stepping
+t = dt
+while t < T + DOLFIN_EPS:
+
+    # Update pressure boundary condition
+    p_in.t = t
+
+    
+    # Compute tentative velocity step
+    begin("Computing tentative velocity")
+    b1 = assemble(L1)
+    [bc.apply(A1, b1) for bc in bcu]
+    solve(A1, u1.vector(), b1, "gmres", "default")
+    end()
+
+    # Pressure correction
+    begin("Computing pressure correction")
+    b2 = assemble(L2)
+    [bc.apply(A2, b2) for bc in bcp]
+    solve(A2, p1.vector(), b2, "gmres", prec)
+    end()
+
+    
+    # Velocity correction
+    begin("Computing velocity correction")
+    b3 = assemble(L3)
+    [bc.apply(A3, b3) for bc in bcu]
+    solve(A3, u1.vector(), b3, "gmres", "default")
+    end()
+
+    # Plot solution
+    plot(p1, title="Pressure", rescale=True)
+    plot(u1, title="Velocity", rescale=True)
+
+    # Save to file
+    ufile << u1
+    pfile << p1
+
+    # Move to next time step
+    u0.assign(u1)
+    t += dt
+    print "t =", t
+
+# Hold plot
+interactive()
+
+# Copyright 2011, The FEniCS Project
+\end{lstlisting}}
+
+\colplacechunks
+\end{parcolumns}
+\end{changemargin}
+
+\section{HyperElasticity}
+This time we compare the code with the c++ version of DOLFIN. 
+The problem for an elastic material can be expressed as a minimization problem
+\begin{align*}
+\min_{u \in V} \Pi\\
+\Pi &= \int_{\Omega} \psi(u) \, {\rm d} x - \int_{\Omega} B \cdot u \, {\rm d} x - \int_{\partial\Omega} T \cdot u \, 
+{\rm d} s\\
+\end{align*}
+where $\Pi$  is the total potential energy, $\psi$ is the elastic stored energy, $B$  is a body force and $T$
+is a traction force.
+
+A complete description of the problem is avilable on the Fenics website \cite{hyperelasticity}.
+The final solution will look like in figure \ref{Hyp}.
+\begin{figure}
+ \begin{center}
+  \includegraphics[height=6 cm,keepaspectratio=true]{./HyperElasticity.png}
+   \caption{Solution of the HyperElasticity problem}
+   \label{Hyp}
+  \end{center}
+\end{figure}
+
+\begin{lstlisting}[caption=UFL code, language=Python, numbers=none]
+# Function spaces
+element = VectorElement("Lagrange", tetrahedron, 1)
+
+# Trial and test functions
+du = TrialFunction(element)     # Incremental displacement
+v  = TestFunction(element)      # Test function
+
+# Functions
+u = Coefficient(element)        # Displacement from previous iteration
+B = Coefficient(element)        # Body force per unit volume
+T = Coefficient(element)        # Traction force on the boundary
+
+# Kinematics
+I = Identity(element.cell().d)  # Identity tensor
+F = I + grad(u)                 # Deformation gradient
+C = F.T*F                       # Right Cauchy-Green tensor
+
+# Invariants of deformation tensors
+Ic = tr(C)
+J  = det(F)
+
+# Elasticity parameters
+mu    = Constant(tetrahedron)
+lmbda = Constant(tetrahedron)
+
+# Stored strain energy density (compressible neo-Hookean model)
+psi = (mu/2)*(Ic - 3) - mu*ln(J) + (lmbda/2)*(ln(J))**2
+
+# Total potential energy
+Pi = psi*dx - inner(B, u)*dx - inner(T, u)*ds
+
+# First variation of Pi (directional derivative about u in the direction of v)
+F = derivative(Pi, u, v)
+
+# Compute Jacobian of F
+J = derivative(F, u, du)
+
+# Copyright 2011, The FEniCS Project
+\end{lstlisting}
+
+\begin{changemargin}{-1.5cm}{-1.5cm}
+$\phantom {u}$
+\begin{parcolumns}[colwidths={1=0.65\textwidth,2=0.65\textwidth}]{2}
+\colchunk{\begin{lstlisting}[caption=Fem-fenics, language=Octave, numbers=none]{Name}
+pkg load fem-fenics msh
+problem = 'HyperElasticity';
+import_ufl_Problem (problem);
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+Rotation = @(x,y,z) ...
+ [0; ...
+ 0.5*(0.5 + (y - 0.5)*cos(pi/3) - (z-0.5)*sin(pi/3) - y);...
+ 0.5*(0.5 + (y - 0.5)*sin(pi/3) + (z-0.5)*cos(pi/3) - z)];
+
+
+
+
+
+
+# Create mesh and define function space
+x = y = z = linspace (0, 1, 17);
+mshd = Mesh (msh3m_structured_mesh (x, y, z, 1, 1:6));
+V  = FunctionSpace (problem, mshd);
+
+
+
+
+
+
+
+# Create Dirichlet boundary conditions
+bcl = DirichletBC (V, @(x,y,z) [0; 0; 0], 1);
+bcr = DirichletBC (V, Rotation, 2);
+bcs = {bcl, bcr};
+
+
+# Define source and boundary traction functions
+
+B = Constant ('B', [0.0; -0.5; 0.0]);
+T = Constant ('T', [0.1; 0.0; 0.0]);
+
+
+
+
+# Set material parameters
+E = 10.0;
+nu = 0.3;
+mu = Constant ('mu', E./(2*(1+nu)));
+lmbda = Constant ('lmbda', E*nu./((1+nu)*(1-2*nu)));
+u = Expression ('u', @(x,y,z) [0; 0; 0]);
+
+# Create (linear) form defining (nonlinear) variational problem
+L = ResidualForm (problem, V, mu, lmbda, B, T, u);
+
+
+
+
+
+
+
+# Solve nonlinear variational problem F(u; v) = 0
+u0 = assemble (L, bcs{:});
+# Create function for the resolution of the NL problem
+function [y, jac] = f (problem, xx, V, bc1, bc2, B, T, mu, lmbda)
+  u = Function ('u', V, xx);
+  a = JacobianForm (problem, V, mu, lmbda, u);
+  L = ResidualForm (problem, V, mu, lmbda, B, T, u);
+  if (nargout == 1)
+    [y, xx] = assemble (L, xx, bc1, bc2);
+  elseif (nargout == 2)
+    [jac, y, xx] = assemble_system (a, L, xx, bc1, bc2);
+  endif
+endfunction
+
+fs = @(xx) f (problem, xx, V, bcl, bcr, B, T, mu, lmbda);
+[x, fval, info] = fsolve (fs, u0, optimset ("jacobian", "on"));
+func = Function ('u', V, x);
+
+# Save solution in VTK format
+save (func, 'displacement');
+
+
+# Plot solution
+plot (func);
+
+
+
+
+#
+\end{lstlisting}}
+
+\colchunk{\begin{lstlisting}[caption=C++, language=C++, numbers=none]{Name}
+#include <dolfin.h>
+#include "HyperElasticity.h"
+
+using namespace dolfin;
+
+// Sub domain for clamp at left end
+class Left : public SubDomain
+{
+  bool inside(const Array<double>& x, bool on_boundary) const
+  {
+    return (std::abs(x[0]) < DOLFIN_EPS) && on_boundary;
+  }
+};
+
+// Sub domain for rotation at right end
+class Right : public SubDomain
+{
+  bool inside(const Array<double>& x, bool on_boundary) const
+  {
+    return (std::abs(x[0] - 1.0) < DOLFIN_EPS) && on_boundary;
+  }
+};
+
+// Dirichlet boundary condition for clamp at left end
+class Clamp : public Expression
+{
+public:
+
+  Clamp() : Expression(3) {}
+
+  void eval(Array<double>& values, const Array<double>& x) const
+  {
+    values[0] = 0.0;
+    values[1] = 0.0;
+    values[2] = 0.0;
+  }
+
+};
+
+// Dirichlet boundary condition for rotation at right end
+class Rotation : public Expression
+{
+public:
+
+  Rotation() : Expression(3) {}
+
+  void eval(Array<double>& values, const Array<double>& x) const
+  {
+    const double scale = 0.5;
+
+    // Center of rotation
+    const double y0 = 0.5;
+    const double z0 = 0.5;
+
+    // Large angle of rotation (60 degrees)
+    double theta = 1.04719755;
+
+    // New coordinates
+    double y = y0 + (x[1]-y0)*cos(theta) - (x[2]-z0)*sin(theta);
+    double z = z0 + (x[1]-y0)*sin(theta) + (x[2]-z0)*cos(theta);
+
+    // Rotate at right end
+    values[0] = 0.0;
+    values[1] = scale*(y - x[1]);
+    values[2] = scale*(z - x[2]);
+  }
+
+};
+
+int main()
+{
+  // Create mesh and define function space
+  UnitCubeMesh mesh (16, 16, 16);
+  HyperElasticity::FunctionSpace V(mesh);
+
+  // Define Dirichlet boundaries
+  Left left;
+  Right right;
+
+  // Define Dirichlet boundary functions
+  Clamp c;
+  Rotation r;
+
+  // Create Dirichlet boundary conditions
+  DirichletBC bcl(V, c, left);
+  DirichletBC bcr(V, r, right);
+  std::vector<const BoundaryCondition*> bcs;
+  bcs.push_back(&bcl); bcs.push_back(&bcr);
+
+  // Define source and boundary traction functions
+  Constant B(0.0, -0.5, 0.0);
+  Constant T(0.1,  0.0, 0.0);
+
+  // Define solution function
+  Function u(V);
+
+  // Set material parameters
+  const double E  = 10.0;
+  const double nu = 0.3;
+  Constant mu(E/(2*(1 + nu)));
+  Constant lambda(E*nu/((1 + nu)*(1 - 2*nu)));
+
+
+  
+  // Create (linear) form defining (nonlinear) variational problem
+  HyperElasticity::ResidualForm F(V);
+  F.mu = mu; F.lmbda = lambda; F.B = B; F.T = T; F.u = u;
+
+  // Create jacobian dF = F' (for use in nonlinear solver).
+  HyperElasticity::JacobianForm J(V, V);
+  J.mu = mu; J.lmbda = lambda; J.u = u;
+
+  // Solve nonlinear variational problem F(u; v) = 0
+  solve(F == 0, u, bcs, J);
+
+
+
+
+
+
+
+
+
+
+
+
+  
+  
+  
+  
+  
+  
+  
+
+
+
+  // Save solution in VTK format
+  File file("displacement.pvd");
+  file << u;
+
+  // Plot solution
+  plot(u);
+  interactive();
+
+  return 0;
+}
+# Copyright 2011, The FEniCS Project
+\end{lstlisting}}
+
+\colplacechunks
+\end{parcolumns}
+\end{changemargin}
+
+\section{Fictitious Domain}
+A penalization method to take into account obstacles in incompressible viscous flows
 
 \newpage 
 \appendix