Commits

Anders Logg  committed aa291d8 Merge with conflicts

Merge branch 'dolfin-demo-documentation'

Conflicts:
cmake/scripts/download-demo-data
demo/documented/mesh-generation/python/demo_mesh_generation.py
demo/undocumented/mesh-generation/cpp/main.cpp
demo/undocumented/mesh-generation/python/demo_mesh-generation.py
demo/undocumented/mesh-generation/python/demo_mesh_generation.py
test/regression/test.py

  • Participants
  • Parent commits 251e291, 61504fb

Comments (0)

Files changed (303)

File cmake/scripts/download-demo-data

     fi
 }
 
-download_mesh demo/la/eigenvalue                             box_with_dent.xml.gz
+download_mesh demo/documented/eigenvalue                     box_with_dent.xml.gz
+download_mesh demo/undocumented/meshfunction                 unitsquare_2_2.xml.gz
+download_mesh demo/undocumented/meshfunction                 unitsquare_2_2_subdomains.xml.gz
 download_mesh demo/undocumented/mixed-poisson-sphere         sphere_16.xml.gz
 download_mesh demo/undocumented/advection-diffusion          dolfin_fine.xml.gz
 download_mesh demo/undocumented/advection-diffusion          dolfin_fine_velocity.xml.gz
 download_mesh demo/undocumented/elasticity                   gear.xml.gz
 download_mesh demo/undocumented/lift-drag                    dolfin_fine.xml.gz
 download_mesh demo/undocumented/lift-drag                    dolfin_fine_pressure.xml.gz
-download_mesh demo/undocumented/tensor-weighted-poisson      unitsquare_32_32.xml.gz
-download_mesh demo/undocumented/tensor-weighted-poisson      unitsquare_32_32_c01.xml.gz
-download_mesh demo/undocumented/tensor-weighted-poisson      unitsquare_32_32_c00.xml.gz
-download_mesh demo/undocumented/tensor-weighted-poisson      unitsquare_32_32_c11.xml.gz
-download_mesh demo/undocumented/stokes-taylor-hood/          dolfin_fine.xml.gz
-download_mesh demo/undocumented/stokes-taylor-hood           dolfin_fine_subdomains.xml.gz
+download_mesh demo/documented/tensor-weighted-poisson        unitsquare_32_32.xml.gz
+download_mesh demo/documented/tensor-weighted-poisson        unitsquare_32_32_c01.xml.gz
+download_mesh demo/documented/tensor-weighted-poisson        unitsquare_32_32_c00.xml.gz
+download_mesh demo/documented/tensor-weighted-poisson        unitsquare_32_32_c11.xml.gz
+download_mesh demo/documented/stokes-taylor-hood/            dolfin_fine.xml.gz
+download_mesh demo/documented/stokes-taylor-hood             dolfin_fine_subdomains.xml.gz
 download_mesh demo/undocumented/stokes-taylor-hood           dolfin_fine.xml.gz
 download_mesh demo/undocumented/elastodynamics               dolfin_fine.xml.gz
 download_mesh demo/undocumented/auto-adaptive-navier-stokes  channel_with_flap.xml.gz
-download_mesh demo/undocumented/stokes-mini                  dolfin_fine.xml.gz
-download_mesh demo/undocumented/stokes-mini                  dolfin_fine_subdomains.xml.gz
-download_mesh demo/undocumented/subdomains                   dolfin_fine.xml.gz
-download_mesh demo/undocumented/stokes-stabilized            dolfin_fine.xml.gz
-download_mesh demo/undocumented/stokes-stabilized            dolfin_fine_subdomains.xml.gz
+download_mesh demo/documented/stokes-mini                    dolfin_fine.xml.gz
+download_mesh demo/documented/stokes-mini                    dolfin_fine_subdomains.xml.gz
+download_mesh demo/documented/subdomains                     dolfin_fine.xml.gz
+download_mesh demo/documented/stokes-stabilized              dolfin_fine.xml.gz
+download_mesh demo/documented/stokes-stabilized              dolfin_fine_subdomains.xml.gz
 download_mesh demo/undocumented/plot                         dolfin_fine.xml.gz
 download_mesh demo/undocumented/dg-advection-diffusion       unitsquare_64_64.xml.gz
 download_mesh demo/undocumented/dg-advection-diffusion       unitsquare_64_64_velocity.xml.gz
-download_mesh demo/undocumented/bcs                          aneurysm.xml.gz
-download_mesh demo/pde/navier-stokes                         lshape.xml.gz
+download_mesh demo/documented/bcs                            aneurysm.xml.gz
+download_mesh demo/documented/navier-stokes                  lshape.xml.gz
 download_mesh test/unit/fem/python                           aneurysm.xml.gz
 download_mesh test/unit/io                                   snake.xml.gz
 download_mesh test/unit/nls/python                           doughnut.xml.gz

File demo/documented/auto-adaptive-poisson/common.txt

+In this demo we will use goal oriented adaptivity and error control
+which applies a duality technique to derive error estimates taken
+directly from the computed solution which then are used to weight
+local residuals. To this end, we derive an :math:`\textit{a
+posteriori}` error estimate and error indicators. We define a goal
+functional :math:`\mathcal{M} : V \rightarrow \mathbb{R}`, which
+expresses the localized physical properties of the solution of a
+simulation. The objective of goal oriented adaptive error control is
+to minimize computational work to obtain a given level of accuracy in
+:math:`\mathcal{M}`.
+
+We will thus illustrate how to:
+
+* Solve a linear partial differential equation with automatic adaptive mesh refinement
+* Define a goal functional
+* Use :py:class:`AdaptiveLinearVariationalSolver <dolfin.cpp.fem.AdaptiveLinearVariationalSolver>`
+
+The two solutions for u in this demo will look as follows, where the
+first is the unrefined while the second is the refined solution:
+
+.. image:: ../u_unrefined.png
+    :scale: 75 %
+
+.. image:: ../u_refined.png
+    :scale: 75 %
+
+
+Equation and problem definition
+-------------------------------
+
+The Poisson equation is the canonical elliptic partial differential
+equation. For a domain :math:`\Omega \subset \mathbb{R}^n` with
+boundary :math:`\partial \Omega = \Gamma_{D} \cup \Gamma_{N}`, the
+Poisson equation with particular boundary conditions reads:
+
+
+.. math::
+
+    - \nabla^{2} u &= f \quad {\rm in} \ \Omega, \\
+    u &= 0 \quad {\rm on} \ \Gamma_{D}, \\
+    \nabla u \cdot n &= g \quad {\rm on} \ \Gamma_{N}. \\
+
+Here, :math:`f` and :math:`g` are input data and n denotes the outward
+directed boundary normal. The variational form of Poisson equation
+reads: find :math:`u \in V` such that
+
+.. math::
+
+	a(u, v) = L(v) \quad \forall \ v \in \hat{V},
+
+which we will call the continous primal problem, where :math:`V`,
+:math:`\hat{V}` are the trial- and test spaces and
+
+
+.. math::
+
+    a(u, v) &= \int_{\Omega} \nabla u \cdot \nabla v \, {\rm d} x, \\
+    L(v)    &= \int_{\Omega} f v \, {\rm d} x + \int_{\Gamma_{N}} g v \, {\rm d} s.
+
+The expression :math:`a(u, v)` is the bilinear form and :math:`L(v)`
+is the linear form. It is assumed that all functions in :math:`V`
+satisfy the Dirichlet boundary conditions (:math:`u = 0 \ {\rm on} \
+\Gamma_{D}`).
+
+The above definitions is that of the continuous problem. In the actual
+computer implementation we use a descrete formulation which reads:
+find :math:`u \in V_h` such that
+
+
+.. math::
+
+    a(u_h, v) = L(v) \quad \forall \ v \in \hat{V}_h.
+
+We will refer to the above equation as the discrete primal
+problem. Here, :math:`V_h` and :math:`\hat{V_h}` are finite
+dimensional subspaces.
+
+The weak residual is defined as
+
+.. math::
+
+    r(v) = L(v) - a(u_h, v).
+
+By the Galerkin orthogonality, we have
+
+
+.. math::
+
+    r(v) = L(v) - a(u_h, v) = a(u_h, v) -  a(u_h, v) = 0\,\, \forall v \in \hat{V}_{h},
+
+which means that the residual vanishes for all functions in
+:math:`\hat{V}_{h}`. This property is used further in the derivation
+of the error estimates. wish to compute a solution :math:`u_h` to the
+discrete primal problem defined above such that for a given tolerance
+:math:`\mathrm{TOL}` we have
+
+.. math::
+
+    \eta = \left| \mathcal{M}(u_h) - \mathcal{u_h} \right| \leq \mathrm{TOL}.
+
+Next we derive an :math:`\textit{a posteriori}` error estimate by
+defining the discrete dual variational problem: find :math:`z_h \in
+V_h^*` such that
+
+.. :math::
+
+    a^*(z_h,v) = \mathcal{v}, \quad \forall v \in \hat{V}_h^*.
+
+Here :math:`V^*, \hat{V}_h^*` are the dual trial- and test spaces and
+:math:`a^* : V^* \times \hat{V}^* \rightarrow \mathbb{R}` is the
+adjoint of :math:`a` such that :math:`a^*(v,w) = a(w,v)`. We find that
+
+.. math::
+
+    \mathcal{M}(u - \mathcal{u_h}) = a^*(z, u - u_h) = a(u - u_h, z) = a(u,z) - a(u_h,z) = L(z) - a(u_h,z) = r(z)
+
+and by Galerkin orthogonality we have :math:`r(z) = r(z - v_h)\,\,
+\forall v_h \in \hat{V}_h`. Note that the residual vanishes if
+:math:`z \in \hat{V}_h^*` and has to either be approximated in a
+higher order element space or one may use an extrapolation. The choice
+of goal functional depends on what quantity you are interested in.
+Here, we take the goal functional to be defined as
+
+.. math::
+
+    \mathcal{M}(u) = \int_{\Omega}  u dx.
+
+
+We use :math:`D\ddot{o}rfler` marking as the mesh marking procedure.
+
+In this demo, we shall consider the following definitions of the input functions, the domain, and the boundaries:
+
+* :math:`\Omega = [0,1] \times [0,1]\,` (a unit square)
+* :math:`\Gamma_{D} = \{(0, y) \cup (1, y) \subset \partial \Omega\}\,` (Dirichlet boundary)
+* :math:`\Gamma_{N} = \{(x, 0) \cup (x, 1) \subset \partial \Omega\}\,` (Neumann boundary)
+* :math:`g = \sin(5x)\,` (normal derivative)
+* :math:`f = 10\exp(-((x - 0.5)^2 + (y - 0.5)^2) / 0.02)\,` (source term)
+

File demo/documented/auto-adaptive-poisson/cpp/AdaptivePoisson.ufl

+# Copyright (C) 2010 Anders Logg and Marie E. Rognes
+#
+# This file is part of DOLFIN.
+#
+# DOLFIN is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# DOLFIN is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with DOLFIN. If not, see <http://www.gnu.org/licenses/>.
+#
+# First added:  2010-08-19
+# Last changed: 2012-12-05
+#
+# Compile this with: ffc -l dolfin -e AdaptivePoisson.ufl
+
+element = FiniteElement("CG", triangle, 1)
+u = TrialFunction(element)
+v = TestFunction(element)
+
+f = Coefficient(element)
+g = Coefficient(element)
+
+a = dot(grad(u), grad(v))*dx()
+L = f*v*dx() + g*v*ds()
+M = u*dx()

File demo/documented/auto-adaptive-poisson/cpp/documentation.rst

+.. Documentation for the auto adaptive poisson demo from DOLFIN.
+
+.. _demo_pde_auto-adaptive-poisson_cpp_documentation:

File demo/documented/auto-adaptive-poisson/cpp/main.cpp

+// Copyright (C) 2010-2012 Anders Logg and Marie E. Rognes
+//
+// This file is part of DOLFIN.
+//
+// DOLFIN is free software: you can redistribute it and/or modify
+// it under the terms of the GNU Lesser General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// DOLFIN is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with DOLFIN. If not, see <http://www.gnu.org/licenses/>.
+//
+// First added:  2010-08-19
+// Last changed: 2012-11-14
+
+#include <dolfin.h>
+#include "AdaptivePoisson.h"
+
+using namespace dolfin;
+
+// Source term (right-hand side)
+class Source : public Expression
+{
+  void eval(Array<double>& values, const Array<double>& x) const
+  {
+    double dx = x[0] - 0.5;
+    double dy = x[1] - 0.5;
+    values[0] = 10*exp(-(dx*dx + dy*dy) / 0.02);
+  }
+};
+
+// Normal derivative (Neumann boundary condition)
+class dUdN : public Expression
+{
+  void eval(Array<double>& values, const Array<double>& x) const
+  {
+    values[0] = sin(5*x[0]);
+  }
+};
+
+// Sub domain for Dirichlet boundary condition
+class DirichletBoundary : public SubDomain
+{
+  bool inside(const Array<double>& x, bool on_boundary) const
+  {
+    return x[0] < DOLFIN_EPS or x[0] > 1.0 - DOLFIN_EPS;
+  }
+};
+
+int main()
+{
+  // Create mesh and define function space
+  UnitSquareMesh mesh(8, 8);
+  AdaptivePoisson::BilinearForm::TrialSpace V(mesh);
+
+  // Define boundary condition
+  Constant u0(0.0);
+  DirichletBoundary boundary;
+  DirichletBC bc(V, u0, boundary);
+
+  // Define variational forms
+  AdaptivePoisson::BilinearForm a(V, V);
+  AdaptivePoisson::LinearForm L(V);
+  Source f;
+  dUdN g;
+  L.f = f;
+  L.g = g;
+
+  // Define Function for solution
+  Function u(V);
+
+  // Define goal functional (quantity of interest)
+  AdaptivePoisson::GoalFunctional M(mesh);
+
+  // Define error tolerance
+  double tol = 1.e-5;
+
+  // Solve equation a = L with respect to u and the given boundary
+  // conditions, such that the estimated error (measured in M) is less
+  // than tol
+  LinearVariationalProblem problem(a, L, u, bc);
+  AdaptiveLinearVariationalSolver solver(problem, M);
+  solver.parameters("error_control")("dual_variational_solver")["linear_solver"] = "cg";
+  solver.solve(tol);
+
+  solver.summary();
+
+  // Plot final solution
+  plot(u.root_node(), "Solution on initial mesh");
+  plot(u.leaf_node(), "Solution on final mesh");
+  interactive();
+
+  return 0;
+}

File demo/documented/auto-adaptive-poisson/python/demo_auto-adaptive_poisson.py

+# Copyright (C) 2011-2012 Marie E. Rognes and Anders Logg
+#
+# This file is part of DOLFIN.
+#
+# DOLFIN is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# DOLFIN is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with DOLFIN. If not, see <http://www.gnu.org/licenses/>.
+#
+# First added:  2010-08-19
+# Last changed: 2012-11-14
+
+# Begin demo
+
+from dolfin import *
+
+# Create mesh and define function space
+mesh = UnitSquareMesh(8, 8)
+V = FunctionSpace(mesh, "Lagrange", 1)
+
+# Define boundary condition
+u0 = Function(V)
+bc = DirichletBC(V, u0, "x[0] < DOLFIN_EPS || x[0] > 1.0 - DOLFIN_EPS")
+
+# Define variational problem
+u = TrialFunction(V)
+v = TestFunction(V)
+f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)",
+               degree=1)
+g = Expression("sin(5*x[0])", degree=1)
+a = inner(grad(u), grad(v))*dx()
+L = f*v*dx() + g*v*ds()
+
+# Define function for the solution
+u = Function(V)
+
+# Define goal functional (quantity of interest)
+M = u*dx()
+
+# Define error tolerance
+tol = 1.e-5
+
+# Solve equation a = L with respect to u and the given boundary
+# conditions, such that the estimated error (measured in M) is less
+# than tol
+problem = LinearVariationalProblem(a, L, u, bc)
+solver = AdaptiveLinearVariationalSolver(problem, M)
+solver.parameters["error_control"]["dual_variational_solver"]["linear_solver"] = "cg"
+solver.solve(tol)
+
+solver.summary()
+
+# Plot solution(s)
+plot(u.root_node(), title="Solution on initial mesh")
+plot(u.leaf_node(), title="Solution on final mesh")
+interactive()

File demo/documented/auto-adaptive-poisson/python/documentation.rst

+.. Documentation for the auto adaptive poisson demo from DOLFIN.
+
+.. _demo_pde_auto-adaptive-poisson_python_documentation:
+
+
+Auto adaptive Poisson equation
+==============================
+
+This demo is implemented in a single Python file,
+:download:`demo_auto-adaptive_poisson.py`, which contains both the
+variational forms and the solver.
+
+.. include:: ../common.txt
+
+Implementation
+--------------
+
+This description goes through the implementation (in
+:download:`demo-auto-adaptive_poisson.py`) of a solver for the above
+described Poisson equation step-by-step.
+
+First, the dolfin module is imported:
+
+.. code-block:: python
+
+	from dolfin import *
+
+We begin by defining a mesh of the domain and a finite element
+function space V relative to this mesh. We used the built-in mesh
+provided by the class :py:class:`UnitSquareMesh
+<dolfin.cpp.mesh.UnitSquareMesh>`. In order to create a mesh
+consisting of 8 x 8 squares with each square divided into two
+triangles, we do as follows:
+
+
+.. code-block:: python
+
+    # Create mesh and define function space
+    mesh = UnitSquareMesh(8, 8)
+    V = FunctionSpace(mesh, "Lagrange", 1)
+
+The second argument to :py:class:`FunctionSpace
+<dolfin.cpp.function.FunctionSpace>`, "Lagrange", is the finite
+element family, while the third argument specifies the polynomial
+degree. Thus, in this case, our space V consists of first-order,
+continuous Lagrange finite element functions (or in order words,
+continuous piecewise linear polynomials).
+
+Next, we want to consider the Dirichlet boundary condition. In our
+case, we want to say that the points (x, y) such that x = 0 or x = 1
+are inside on the inside of :math:`\Gamma_D`. (Note that because of
+rounding-off errors, it is often wise to instead specify :math:`x <
+\epsilon` or :math:`x > 1 - \epsilon` where :math:`\epsilon` is a
+small number (such as machine precision).)
+
+
+.. code-block:: python
+
+    # Define boundary condition
+    u0 = Function(V)
+    bc = DirichletBC(V, u0, "x[0] < DOLFIN_EPS || x[0] > 1.0 - DOLFIN_EPS")
+
+Next, we want to express the variational problem. First, we need to
+specify the trial function u and the test function v, both living in
+the function space V. We do this by defining a
+:py:class:`TrialFunction <dolfin.functions.function.TrialFunction>`
+and a :py:class:`TestFunction
+<dolfin.functions.function.TestFunction>` on the previously defined
+:py:class:`FunctionSpace <dolfin.cpp.function.FunctionSpace>` V.
+
+Further, the source f and the boundary normal derivative g are
+involved in the variational forms, and hence we must specify
+these. Both f and g are given by simple mathematical formulas, and can
+be easily declared using the :py:class:`Expression
+<dolfin.cpp.function.Expression>` class. Note that the strings
+defining f and g use C++ syntax since, for efficiency, DOLFIN will
+generate and compile C++ code for these expressions at run-time.
+
+With these ingredients, we can write down the bilinear form a and the
+linear form L (using UFL operators). In summary, this reads
+
+.. code-block:: python
+
+    # Define variational problem
+    u = TrialFunction(V)
+    v = TestFunction(V)
+    f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)",
+                   degree=1)
+    g = Expression("sin(5*x[0])", degree=1)
+    a = inner(grad(u), grad(v))*dx()
+    L = f*v*dx() + g*v*ds()
+
+Now, we have specified the variational forms and can consider the
+solution of the variational problem. First, we need to define a
+:py:class:`Function <dolfin.cpp.function.Function>` u to represent the
+solution. (Upon initialization, it is simply set to the zero
+function.) A Function represents a function living in a finite element
+function space.
+
+.. code-block:: python
+
+    # Define function for the solution
+    u = Function(V)
+
+Then define the goal functional:
+
+.. code-block:: python
+
+    # Define goal functional (quantity of interest)
+    M = u*dx()
+
+Next we specify the error tolerance for when the refinement shall stop:
+
+.. code-block:: python
+
+    # Define error tolerance
+    tol = 1.e-5
+
+Now, we have specified the variational forms and can consider the
+solution of the variational problem. First, we define the
+:py:class:`LinearVariationalProblem
+<dolfin.cpp.fem.LinearVariationalProblem>` function with the arguments
+a, L, u and bc. Next we send this problem to the
+:py:class:`AdaptiveLinearVariationalSolver
+<dolfin.cpp.fem.AdaptiveLinearVariationalSolver>` together with the
+goal functional. Note that one may also choose several adaptations in
+the error control. At last we solve the problem with the defined
+tolerance:
+
+.. code-block:: python
+
+    # Solve equation a = L with respect to u and the given boundary
+    # conditions, such that the estimated error (measured in M) is less
+    # than tol
+    problem = LinearVariationalProblem(a, L, u, bc)
+    solver = AdaptiveLinearVariationalSolver(problem, M)
+    solver.parameters["error_control"]["dual_variational_solver"]["linear_solver"] = "cg"
+    solver.solve(tol)
+
+    solver.summary()
+
+    # Plot solution(s)
+    plot(u.root_node(), title="Solution on initial mesh")
+    plot(u.leaf_node(), title="Solution on final mesh")
+    interactive()
+
+Complete code
+-------------
+
+.. literalinclude:: demo_auto-adaptive_poisson.py
+	:start-after: # Begin demo

File demo/documented/auto-adaptive-poisson/u_refined.png

Added
New image

File demo/documented/auto-adaptive-poisson/u_unrefined.png

Added
New image

File demo/documented/bcs/common.txt

+This demo illustrates how to:
+
+* Use meshes with boundary indicators
+* Solve the Poisson equation with Dirichlet boundary values
+
+
+The solution for :math:`u` in this demo will look as follows:
+
+.. image:: ../demo_bcs.png
+    :scale: 75 %
+
+Equation and the problem definition
+-----------------------------------
+We will use the Poisson equation as a model problem when we demonstrate
+how to set boundary conditions for an imported mesh that includes boundary indicators.
+The Poisson equation is the canonical elliptic partial differential equation.
+For a domain :math:`\Omega \subset \mathbb{R}^3` with boundary :math:`\partial \Omega = \cup_{i = 0}^3 \Gamma_{D, i}`,
+the Poisson equation with particular boundary conditions reads:
+
+.. math::
+	- \nabla^{2} u &= f \quad {\rm in} \ \Omega, \\
+             u &= u_0 \quad {\rm on} \ \Gamma_{D, 0}, \\
+             u &= u_1 \quad {\rm on} \ \Gamma_{D, 1}, \\
+             u &= u_2 \quad {\rm on} \ \Gamma_{D, 2}, \\
+             u &= u_3 \quad {\rm on} \ \Gamma_{D, 3}. \\
+
+Here, :math:`f` is some given input data and :math:`u_0`, :math:`u_1`, :math:`u_2` and :math:`u_3`
+are the prescribed values of :math:`u` at the boundaries.
+The variational form of the Poisson equation reads: find :math:`u \in V` such that
+
+.. math::
+	a(u, v) = L(v) \quad \forall \ v \in V,
+
+where :math:`V` is a suitable function space and
+
+.. math::
+	a(u, v) &= \int_{\Omega} \nabla u \cdot \nabla v \, {\rm d} x, \\
+	L(v)    &= \int_{\Omega} f v \, {\rm d} x.
+
+The expression :math:`a(u, v)` is the bilinear form and :math:`L(v)` is the linear form.
+It is assumed that all functions in :math:`V` satisfy the Dirichlet boundary conditions.
+
+In this demo we shall consider the domain :math:`\Omega` to be a model
+of blood vessels with an aneurysm.  It has one inlet vessel and two
+outlet vessels.  We define noslip boundary conditions on the walls of
+the vessels and aneurysm, that is :math:`u = u_0 = 0.0`.  We let
+:math:`u = u_1 = 1.0` be the Dirichlet condition on the inlet, the
+outlets will have the prescribed values :math:`u = u_2 = 2.0` and
+:math:`u = u_3 = 3.0`.  In summary, we have:
+
+* :math:`u = u_0 = 0.0 \text{ on } \Gamma_{D, 0}`  (noslip boundary)
+* :math:`u = u_1 = 1.0 \text{ on } \Gamma_{D, 1}`  (inlet)
+* :math:`u = u_2 = 2.0 \text{ on } \Gamma_{D, 2}`  (outlet 1)
+* :math:`u = u_3 = 3.0 \text{ on } \Gamma_{D, 3}`  (outlet 2)
+* :math:`f = 0.0`  (source term)

File demo/documented/bcs/cpp/Poisson.ufl

+# Copyright (C) 2008 Anders Logg
+#
+# This file is part of DOLFIN.
+#
+# DOLFIN is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# DOLFIN is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with DOLFIN. If not, see <http://www.gnu.org/licenses/>.
+#
+# First added:  2008-05-23
+# Last changed: 2011-03-09
+#
+# The bilinear form a(u, v) and linear form L(v) for
+# Poisson's equation.
+#
+# Compile these forms with FFC: ffc -l dolfin Poisson.ufl
+
+element = FiniteElement("Lagrange", tetrahedron, 1)
+
+v = TestFunction(element)
+u = TrialFunction(element)
+f = Coefficient(element)
+
+a = dot(grad(u), grad(v))*dx
+L = f*v*dx

File demo/documented/bcs/cpp/documentation.rst

+.. Documentation for the bcs demo from DOLFIN.
+
+.. _demo_pde_bcs_cpp_documentation:
+

File demo/documented/bcs/cpp/main.cpp

+// Copyright (C) 2008 Anders Logg
+//
+// This file is part of DOLFIN.
+//
+// DOLFIN is free software: you can redistribute it and/or modify
+// it under the terms of the GNU Lesser General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// DOLFIN is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with DOLFIN. If not, see <http://www.gnu.org/licenses/>.
+//
+// First added:  2008-05-23
+// Last changed: 2012-07-05
+//
+// This demo illustrates how to set boundary conditions for meshes
+// that include boundary indicators. The mesh used in this demo was
+// generated with VMTK (http://www.vmtk.org/).
+
+#include <dolfin.h>
+#include "Poisson.h"
+#include <boost/assign/list_of.hpp>
+
+using namespace dolfin;
+
+int main()
+{
+  // Create mesh and finite element
+  Mesh mesh("../aneurysm.xml.gz");
+
+  // Define variational problem
+  Constant f(0.0);
+  Poisson::FunctionSpace V(mesh);
+  Poisson::BilinearForm a(V, V);
+  Poisson::LinearForm L(V);
+  L.f = f;
+
+  // Define boundary condition values
+  Constant u0(0.0);
+  Constant u1(1.0);
+  Constant u2(2.0);
+  Constant u3(3.0);
+
+  // Define boundary conditions
+  DirichletBC bc0(V, u0, 0);
+  DirichletBC bc1(V, u1, 1);
+  DirichletBC bc2(V, u2, 2);
+  DirichletBC bc3(V, u3, 3);
+  std::vector<const DirichletBC*> bcs = boost::assign::list_of(&bc0)(&bc1)(&bc2)(&bc3);
+
+  // Compute solution
+  Function u(V);
+  solve(a == L, u, bcs);
+
+  // Write solution to file
+  File file("u.pvd");
+  file << u;
+
+  // Plot solution
+  plot(u);
+  interactive();
+
+  return 0;
+}

File demo/documented/bcs/demo_bcs.png

Added
New image

File demo/documented/bcs/python/demo_bcs.py

+"""This demo illustrates how to set boundary conditions for meshes
+that include boundary indicators. The mesh used in this demo was
+generated with VMTK (http://www.vmtk.org/)."""
+
+# Copyright (C) 2008-2012 Anders Logg
+#
+# This file is part of DOLFIN.
+#
+# DOLFIN is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# DOLFIN is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with DOLFIN. If not, see <http://www.gnu.org/licenses/>.
+#
+# Modified by Martin Alnaes 2012
+#
+# First added:  2008-05-23
+# Last changed: 2012-10-16
+# Begin demo
+
+from dolfin import *
+
+# Create mesh and define function space
+mesh = Mesh("../aneurysm.xml.gz")
+V = FunctionSpace(mesh, "CG", 1)
+
+# Define variational problem
+u = TrialFunction(V)
+v = TestFunction(V)
+f = Constant(0.0)
+a = dot(grad(u), grad(v))*dx
+L = f*v*dx
+
+# Define boundary condition values
+u0 = Constant(0.0)
+u1 = Constant(1.0)
+u2 = Constant(2.0)
+u3 = Constant(3.0)
+
+# Define boundary conditions
+bc0 = DirichletBC(V, u0, 0)
+bc1 = DirichletBC(V, u1, 1)
+bc2 = DirichletBC(V, u2, 2)
+bc3 = DirichletBC(V, u3, 3)
+
+# Compute solution
+u = Function(V)
+solve(a == L, u, [bc0, bc1, bc2, bc3])
+
+# Write solution to file
+File("u.pvd") << u
+
+# Plot solution
+plot(u, interactive=True)

File demo/documented/bcs/python/documentation.rst

+.. Documentation for the bcs demo from DOLFIN.
+
+.. _demo_pde_bcs_documentation:
+
+
+Set boundary conditions for meshes that include boundary indicators
+===================================================================
+
+This demo is implemented in a single python file, :download:`demo_bcs.py`,
+which contains both the variational form and the solver.
+
+
+.. include:: ../common.txt
+
+Implementation
+--------------
+
+This description goes through the implementation (in
+:download:`demo_bcs.py`) of a solver for the above described Poisson
+equation and how to set boundary conditions for a mesh that includes
+boundary indicators.
+
+First, the :py:mod:`dolfin` module is imported:
+
+.. code-block:: python
+
+	from dolfin import *
+
+Then, we import the mesh and create a finite element function space
+:math:`V` relative to this mesh. In this case we create a
+:py:class:`FunctionSpace<dolfin.functions.functionspace.FunctionSpace>`
+``V`` consisting of continuous piecewise linear polynomials.
+
+.. code-block:: python
+
+	# Create mesh and define function space
+	mesh = Mesh("../aneurysm.xml.gz")
+	V = FunctionSpace(mesh, "CG", 1)
+
+Now, we define the trial function u and the test function v, both living in the function space ``V``. We also define our variational problem, a and L. u and v are defined using the classes
+:py:class:`TrialFunction <dolfin.functions.function.TrialFunction>` and
+:py:class:`TestFunction<dolfin.functions.function.TrialFunction>`,
+respetively, on the :py:class:`FunctionSpace<dolfin.functions.functionspace.FunctionSpace>` ``V``.
+The source :math:`f` may be defined as a :py:class:`Constant <dolfin.functions.constant.Constant>`.
+The bilinear and linear forms, ``a`` and ``L`` respectively, are defined using UFL operators.
+Thus, the definition of the variational problem reads:
+
+.. code-block:: python
+
+	# Define variational problem
+	u = TrialFunction(V)
+	v = TestFunction(V)
+	f = Constant(0.0)
+	a = dot(grad(u), grad(v))*dx
+	L = f*v*dx
+
+Before we can solve the problem we must specify the boundary conditions.
+We begin with specifying the values of the boundary conditions as :py:class:`Constant <dolfin.functions.constant.Constant>`s.
+Then we use the class :py:class:`DirichletBC <dolfin.fem.bcs.DirichletBC>` to define the Dirichlet boundary conditions:
+
+.. code-block:: python
+
+	# Define boundary condition values
+	u0 = Constant(0.0)
+	u1 = Constant(1.0)
+	u2 = Constant(2.0)
+	u3 = Constant(3.0)
+
+	# Define boundary conditions
+	bc0 = DirichletBC(V, u0, 0)
+	bc1 = DirichletBC(V, u1, 1)
+	bc2 = DirichletBC(V, u2, 2)
+	bc3 = DirichletBC(V, u3, 3)
+
+:py:class:`DirichletBC <dolfin.fem.bcs.DirichletBC>` takes three arguments, the first one is our function space ``V``,
+the next is the boundary condition value and the third is the subdomain indicator which is information stored in the mesh.
+
+At this point we are ready to create a :py:class:`Function <dolfin.cpp.function.Function>` ``u`` to store the solution and call the solve function with the arguments
+``a == L``, ``u``, ``[bc0, bc1, bc2, bc3]``, as follows:
+
+.. code-block:: python
+
+	# Compute solution
+	u = Function(V)
+	solve(a == L, u, [bc0, bc1, bc2, bc3])
+
+When we have solved the problem, we save the solution to file and plot it.
+
+.. code-block:: python
+
+	# Write solution to file
+	File("u.pvd") << u
+
+	# Plot solution
+	plot(u, interactive=True)
+
+
+Complete code
+-------------
+
+.. literalinclude:: demo_bcs.py
+   :start-after: # Begin demo
+

File demo/documented/biharmonic/biharmonic_u.png

Added
New image

File demo/documented/biharmonic/common.txt

+
+This demo illustrates how to:
+
+* Solve a linear partial differential equation
+* Use a discontinuous Galerkin method
+* Solve a fourth-order differential equation
+
+The solution for :math:`u` in this demo will look as follows:
+
+.. image:: ../biharmonic_u.png
+    :scale: 75 %
+
+Equation and problem definition
+-------------------------------
+
+The biharmonic equation is a fourth-order elliptic equation. On the domain
+:math:`\Omega \subset \mathbb{R}^{d}`, :math:`1 \le d \le 3`, it reads
+
+.. math::
+   \nabla^{4} u = f \quad {\rm in} \ \Omega,
+
+where :math:`\nabla^{4} \equiv \nabla^{2} \nabla^{2}` is the biharmonic
+operator and :math:`f` is a prescribed source term. To formulate a complete
+boundary value problem, the biharmonic equation must be complemented by
+suitable boundary conditions.
+
+Multiplying the biharmonic equation by a test function and integrating
+by parts twice leads to a problem second-order derivatives, which
+would requires :math:`H^{2}` conforming (roughly :math:`C^{0}`
+continuous) basis functions.  To solve the biharmonic equation using
+Lagrange finite element basis functions, the biharmonic equation can
+be split into two second-order equations (see the Mixed Poisson demo
+for a mixed method for the Poisson equation), or a variational
+formulation can be constructed that imposes weak continuity of normal
+derivatives between finite element cells.  The demo uses a
+discontinuous Galerkin approach to impose continuity of the normal
+derivative weakly.
+
+Consider a triangulation :math:`\mathcal{T}` of the domain :math:`\Omega`,
+where the union of interior facets is denoted by :math:`\Gamma`.
+Functions evaluated on opposite sides of a facet are indicated by the
+subscripts ':math:`+`' and ':math:`-`'.
+Using the standard continuous Lagrange finite element space
+
+.. math::
+    V_ = \left\{v \in H^{1}_{0}(\Omega): \ v \in P_{k}(K) \ \forall \ K \in \mathcal{T} \right\}
+
+and considering the boundary conditions
+
+.. math::
+   u            &= 0 \quad {\rm on} \ \partial\Omega \\
+   \nabla^{2} u &= 0 \quad {\rm on} \ \partial\Omega
+
+a weak formulation of the biharmonic reads: find :math:`u \in V` such that
+
+.. math::
+   \sum_{K \in \mathcal{T}} \int_{K} \nabla^{2} u \nabla^{2} v \, dx \
+  - \int_{\Gamma} \left<\nabla^{2} u \right> \llbracket\nabla v \rrbracket \, ds
+  - \int_{\Gamma} \llbracket\nabla u \rrbracket \left<\nabla^{2} v \right>  \, ds
+  + \int_{\Gamma} \frac{\alpha}{h} \llbracket \nabla u \rrbracket \llbracket \nabla v \rrbracket \, ds
+  = \int_{\Omega} vf \, dx \quad \forall \ v \in V
+
+where :math:`\left< u \right> = (1/2) (u_{+} + u_{-})`, :math:`\llbracket w
+\rrbracket = w_{+} \cdot n_{+} + w_{-} \cdot n_{-}`, :math:`\alpha \ge 0`
+is a penalty term and :math:`h` is a measure of the cell size.  For the
+implementation, it is useful to identify the bilinear form
+
+.. math::
+   a(u, v) = \sum_{K \in \mathcal{T}} \int_{K} \nabla^{2} u \nabla^{2} v \, dx \
+  - \int_{\Gamma} \left<\nabla^{2} u \right> \llbracket\nabla v \rrbracket \, ds
+  - \int_{\Gamma} \llbracket\nabla u \rrbracket \left<\nabla^{2} v \right>  \, ds
+  + \int_{\Gamma} \frac{\alpha}{h} \llbracket \nabla u \rrbracket \llbracket \nabla v \rrbracket \, ds
+
+and the linear form
+
+.. math::
+  L(v) = \int_{\Omega} vf \, dx
+
+The input parameters for this demos are defined as follows:
+
+* :math:`\Omega = [0,1] \times [0,1]` (a unit square)
+* :math:`\alpha = 8.0` (penalty parameter)
+* :math:`f = 4.0 \pi^4\sin(\pi x)\sin(\pi y)` (source term)

File demo/documented/biharmonic/cpp/Biharmonic.ufl

+# Copyright (C) 2009 Kristian B. Oelgaard, Garth N. Wells and Anders Logg
+#
+# This file is part of DOLFIN.
+#
+# DOLFIN is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# DOLFIN is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with DOLFIN. If not, see <http://www.gnu.org/licenses/>.
+#
+# First added:  2009-06-26
+# Last changed: 2011-05-11
+#
+# The bilinear form a(u, v) and linear form L(v) for
+# Biharmonic equation in a discontinuous Galerkin (DG)
+# formulation.
+#
+# Compile this form with FFC: ffc -l dolfin Biharmonic.ufl
+
+# Elements
+element = FiniteElement("Lagrange", triangle, 2)
+
+# Trial and test functions
+u = TrialFunction(element)
+v = TestFunction(element)
+f = Coefficient(element)
+
+# Normal component, mesh size and right-hand side
+n  = element.cell().n
+h = 2.0*triangle.circumradius
+h_avg = (h('+') + h('-'))/2
+
+# Parameters
+alpha = Constant(triangle)
+
+# Bilinear form
+a = inner(div(grad(u)), div(grad(v)))*dx \
+  - inner(avg(div(grad(u))), jump(grad(v), n))*dS \
+  - inner(jump(grad(u), n), avg(div(grad(v))))*dS \
+  + alpha('+')/h_avg*inner(jump(grad(u), n), jump(grad(v),n))*dS
+
+# Linear form
+L = f*v*dx

File demo/documented/biharmonic/cpp/documentation.rst

+.. Documentation for the biharmonic demo from DOLFIN.
+
+.. _demo_pde_biharmonic_cpp_documentation:
+
+
+Biharmonic equation
+===================
+
+.. include:: ../common.txt
+
+
+Implementation
+--------------
+
+The implementation is split in two files, a form file containing the definition
+of the variational forms expressed in UFL and the solver which is implemented
+in a C++ file.
+
+Running this demo requires the files: :download:`main.cpp`,
+:download:`Biharmonic.ufl` and :download:`CMakeLists.txt`.
+
+UFL form file
+^^^^^^^^^^^^^
+
+First we define the variational problem in UFL in the file called
+:download:`Biharmonic.ufl`.
+
+In the UFL file, the finite element space is defined:
+
+.. code-block:: python
+
+    # Elements
+    element = FiniteElement("Lagrange", triangle, 2)
+
+On the space ``element``, trial and test functions, and the source term
+are defined:
+
+.. code-block:: python
+
+    # Trial and test functions
+    u = TrialFunction(element)
+    v = TestFunction(element)
+    f = Coefficient(element)
+
+Next, the outward unit normal to cell boundaries and a measure of the
+cell size are defined. The average size of cells sharing a facet will
+be used (``h_avg``).  The UFL syntax ``('+')`` and ``('-')`` restricts
+a function to the ``('+')`` and ``('-')`` sides of a facet,
+respectively.  The penalty parameter ``alpha`` is made a
+:cpp:class:`Constant` so that it can be changed in the program without
+regenerating the code.
+
+.. code-block:: python
+
+    # Normal component, mesh size and right-hand side
+    n  = element.cell().n
+    h = 2.0*triangle.circumradius
+    h_avg = (h('+') + h('-'))/2
+
+    # Parameters
+    alpha = Constant(triangle)
+
+Finally the bilinear and linear forms are defined. Integrals over
+internal facets are indicated by ``*dS``.
+
+.. code-block:: python
+
+    # Bilinear form
+    a = inner(div(grad(u)), div(grad(v)))*dx \
+      - inner(avg(div(grad(u))), jump(grad(v), n))*dS \
+      - inner(jump(grad(u), n), avg(div(grad(v))))*dS \
+      + alpha('+')/h_avg*inner(jump(grad(u), n), jump(grad(v),n))*dS
+
+    # Linear form
+    L = f*v*dx
+
+
+C++ program
+^^^^^^^^^^^
+
+The DOLFIN interface and the code generated from the UFL input is included,
+and the DOLFIN namespace is used:
+
+.. code-block:: c++
+
+  #include <dolfin.h>
+  #include "Biharmonic.h"
+
+  using namespace dolfin;
+
+A class ``Source`` is defined for the function :math:`f`, with the
+function ``Expression::eval`` overloaded:
+
+.. code-block:: c++
+
+  // Source term
+  class Source : public Expression
+  {
+  public:
+
+    void eval(Array<double>& values, const Array<double>& x) const
+    {
+    values[0] = 4.0*std::pow(DOLFIN_PI, 4)*
+      std::sin(DOLFIN_PI*x[0])*std::sin(DOLFIN_PI*x[1]);
+    }
+
+  };
+
+A boundary subdomain is defined, which in this case is the entire boundary:
+
+.. code-block:: c++
+
+  // Sub domain for Dirichlet boundary condition
+  class DirichletBoundary : public SubDomain
+  {
+    bool inside(const Array<double>& x, bool on_boundary) const
+    {
+      return on_boundary;
+    }
+  };
+
+The main part of the program is begun, and a mesh is created with 32 vertices
+in each direction:
+
+.. code-block:: c++
+
+  int main()
+  {
+    // Create mesh
+    UnitSquareMesh mesh(32, 32);
+
+
+The source function, a function for the cell size and the penalty term
+are declared:
+
+.. code-block:: c++
+
+    // Create functions
+    Source f;
+    Constant alpha(8.0);
+
+A function space object, which is defined in the generated code, is created:
+
+.. code-block:: c++
+
+    // Create function space
+    Biharmonic::FunctionSpace V(mesh);
+
+The Dirichlet boundary condition on :math:`u` is constructed by
+defining a :cpp:class:`Constant` which is equal to zero, defining the boundary
+(``DirichletBoundary``), and using these, together with ``V``, to create
+``bc``:
+
+.. code-block:: c++
+
+    // Define boundary condition
+    Constant u0(0.0);
+    DirichletBoundary boundary;
+    DirichletBC bc(V, u0, boundary);
+
+Using the function space ``V``, the bilinear and linear forms
+are created, and function are attached:
+
+.. code-block:: c++
+
+    // Define variational problem
+    Biharmonic::BilinearForm a(V, V);
+    Biharmonic::LinearForm L(V);
+    a.alpha = alpha; L.f = f;
+
+A :cpp:class:`Function` is created to hold the solution and the
+problem is solved:
+
+.. code-block:: c++
+
+    // Compute solution
+    Function u(V);
+    solve(a == L, u, bc);
+
+The solution is then written to a file in VTK format and plotted to
+the screen:
+
+.. code-block:: c++
+
+    // Save solution in VTK format
+    File file("biharmonic.pvd");
+    file << u;
+
+    // Plot solution
+    plot(u);
+    interactive();
+
+    return 0;
+  }
+
+Complete code
+-------------
+
+Complete UFL file
+^^^^^^^^^^^^^^^^^
+
+.. literalinclude:: Biharmonic.ufl
+   :start-after: # Compile
+   :language: python
+
+Complete main file
+^^^^^^^^^^^^^^^^^^
+
+.. literalinclude:: main.cpp
+   :start-after: // using
+   :language: c++

File demo/documented/biharmonic/cpp/main.cpp

+// Copyright (C) 2009 Kristian B. Oelgaard
+//
+// This file is part of DOLFIN.
+//
+// DOLFIN is free software: you can redistribute it and/or modify
+// it under the terms of the GNU Lesser General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// DOLFIN is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with DOLFIN. If not, see <http://www.gnu.org/licenses/>.
+//
+// Modified by Anders Logg, 2011
+//
+// First added:  2009-06-26
+// Last changed: 2012-11-12
+//
+// This demo program solves the Biharmonic equation,
+//
+//     - nabla^4 u(x, y) = f(x, y)
+//
+// on the unit square with source f given by
+//
+//     f(x, y) = 4 pi^4 sin(pi*x)*sin(pi*y)
+//
+// and boundary conditions given by
+//
+//     u(x, y)         = 0
+//     nabla^2 u(x, y) = 0
+//
+// using a discontinuous Galerkin formulation (interior penalty method).
+
+#include <dolfin.h>
+#include "Biharmonic.h"
+
+using namespace dolfin;
+
+// Source term
+class Source : public Expression
+{
+public:
+
+  void eval(Array<double>& values, const Array<double>& x) const
+  {
+    values[0] = 4.0*std::pow(DOLFIN_PI, 4)*
+      std::sin(DOLFIN_PI*x[0])*std::sin(DOLFIN_PI*x[1]);
+  }
+
+};
+
+// Sub domain for Dirichlet boundary condition
+class DirichletBoundary : public SubDomain
+{
+  bool inside(const Array<double>& x, bool on_boundary) const
+  {
+    return on_boundary;
+  }
+};
+
+int main()
+{
+  // Create mesh
+  UnitSquareMesh mesh(32, 32);
+
+  // Create functions
+  Source f;
+  Constant alpha(8.0);
+
+  // Create function space
+  Biharmonic::FunctionSpace V(mesh);
+
+  // Define boundary condition
+  Constant u0(0.0);
+  DirichletBoundary boundary;
+  DirichletBC bc(V, u0, boundary);
+
+  // Define variational problem
+  Biharmonic::BilinearForm a(V, V);
+  Biharmonic::LinearForm L(V);
+  a.alpha = alpha; L.f = f;
+
+  // Compute solution
+  Function u(V);
+  solve(a == L, u, bc);
+
+  // Save solution in VTK format
+  File file("biharmonic.pvd");
+  file << u;
+
+  // Plot solution
+  plot(u);
+  interactive();
+
+  return 0;
+}

File demo/documented/biharmonic/python/demo_biharmonic.py

+"""This demo program solves the Biharmonic equation,
+
+    nabla^4 u(x, y) = f(x, y)
+
+on the unit square with source f given by
+
+    f(x, y) = 4 pi^4 sin(pi*x)*sin(pi*y)
+
+and boundary conditions given by
+
+    u(x, y)         = 0
+    nabla^2 u(x, y) = 0
+
+using a discontinuous Galerkin formulation (interior penalty method).
+"""
+
+# Copyright (C) 2009 Kristian B. Oelgaard
+#
+# This file is part of DOLFIN.
+#
+# DOLFIN is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# DOLFIN is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with DOLFIN. If not, see <http://www.gnu.org/licenses/>.
+#
+# Modified by Anders Logg, 2011
+#
+# First added:  2009-06-26
+# Last changed: 2012-11-12
+
+# Begin demo
+
+from dolfin import *
+
+# Optimization options for the form compiler
+parameters["form_compiler"]["cpp_optimize"] = True
+parameters["form_compiler"]["optimize"] = True
+
+# Create mesh and define function space
+mesh = UnitSquareMesh(32, 32)
+V = FunctionSpace(mesh, "CG", 2)
+
+# Define Dirichlet boundary
+class DirichletBoundary(SubDomain):
+    def inside(self, x, on_boundary):
+        return on_boundary
+
+class Source(Expression):
+    def eval(self, values, x):
+        values[0] = 4.0*pi**4*sin(pi*x[0])*sin(pi*x[1])
+
+# Define boundary condition
+u0 = Constant(0.0)
+bc = DirichletBC(V, u0, DirichletBoundary())
+
+# Define trial and test functions
+u = TrialFunction(V)
+v = TestFunction(V)
+
+# Define normal component, mesh size and right-hand side
+h = CellSize(mesh)
+h_avg = (h('+') + h('-'))/2.0
+n = FacetNormal(mesh)
+f = Source()
+
+# Penalty parameter
+alpha = Constant(8.0)
+
+# Define bilinear form
+a = inner(div(grad(u)), div(grad(v)))*dx \
+  - inner(avg(div(grad(u))), jump(grad(v), n))*dS \
+  - inner(jump(grad(u), n), avg(div(grad(v))))*dS \
+  + alpha('+')/h_avg*inner(jump(grad(u),n), jump(grad(v),n))*dS
+
+# Define linear form
+L = f*v*dx
+
+# Solve variational problem
+u = Function(V)
+solve(a == L, u, bc)
+
+# Save solution to file
+file = File("biharmonic.pvd")
+file << u
+
+# Plot solution
+plot(u, interactive=True)

File demo/documented/biharmonic/python/documentation.rst

+.. Documentation for the biharmonic demo from DOLFIN.
+
+.. _demo_pde_biharmonic_python_documentation:
+
+Biharmonic equation
+===================
+
+This demo is implemented in a single Python file,
+:download:`demo_biharmonic.py`, which contains both the variational
+forms and the solver.
+
+.. include:: ../common.txt
+
+
+Implementation
+--------------
+
+This demo is implemented in the :download:`demo_biharmonic.py` file.
+
+First, the :py:mod:`dolfin` module is imported:
+
+.. code-block:: python
+
+    from dolfin import *
+
+Next, some parameters for the form compiler are set:
+
+.. code-block:: python
+
+    # Optimization options for the form compiler
+    parameters["form_compiler"]["cpp_optimize"] = True
+    parameters["form_compiler"]["optimize"] = True
+
+A mesh is created, and a quadratic finite element function space:
+
+.. code-block:: python
+
+    # Create mesh and define function space
+    mesh = UnitSquareMesh(32, 32)
+    V = FunctionSpace(mesh, "CG", 2)
+
+A subclass of :py:class:`SubDomain <dolfin.cpp.SubDomain>`,
+``DirichletBoundary`` is created for later defining the boundary of
+the domian:
+
+.. code-block:: python
+
+    # Define Dirichlet boundary
+    class DirichletBoundary(SubDomain):
+        def inside(self, x, on_boundary):
+            return on_boundary
+
+A subclass of :py:class:`Expression
+<dolfin.functions.expression.Expression>`, ``Source`` is created for
+the source term :math:`f`:
+
+.. code-block:: python
+
+    class Source(Expression):
+        def eval(self, values, x):
+            values[0] = 4.0*pi**4*sin(pi*x[0])*sin(pi*x[1])
+
+The Dirichlet boundary condition is created:
+
+.. code-block:: python
+