--- a/example/Biharmonic/Biharmonic.m	Thu Sep 12 15:19:31 2013 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,24 +0,0 @@
-pkg load fem-fenics msh
-
-problem = 'Biharmonic';
-import_ufl_Problem (problem);
-
-# Create mesh and define function space
-x = y = linspace (0, 1, 33);
-mesh = Mesh(msh2m_structured_mesh (x, y, 1, 1:4));
-
-V = FunctionSpace(problem, mesh);
-
-bc = DirichletBC (V, @(x,y) 0, 1:4);
-
-f = Expression ('f', @(x,y) 4.0*pi^4.*sin(pi.*x).*sin(pi.*y));
-g = Expression ('alpha', @(x,y) 8);
-
-a = BilinearForm (problem, V, g);
-L = LinearForm (problem, V, f);
-
-[A, b] = assemble_system (a, L, bc);
-u = A \ b;
-
-func = Function ('u', V, u);
-plot (func);

--- a/example/Biharmonic/Biharmonic.ufl	Thu Sep 12 15:19:31 2013 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,26 +0,0 @@
-# Copyright (C) 2005-2009 Anders Logg
-
-# 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

--- a/example/Elasticity/HyperElasticity.m	Thu Sep 12 15:19:31 2013 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,57 +0,0 @@
-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
-
-# 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
-
-u0 = assemble (L, bcs{:});
-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);

--- a/example/Elasticity/HyperElasticity.ufl	Thu Sep 12 15:19:31 2013 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,63 +0,0 @@
-# Copyright (C) 2009-2010 Harish Narayanan and Garth N. Wells
-#
-# 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-09-29
-# Last changed: 2011-06-28
-#
-# The bilinear form a(du, v) and linear form L(v) for
-# a hyperelastic model.
-#
-# Compile this form with FFC: ffc -l dolfin -feliminate_zeros -fprecompute_basis_const -fprecompute_ip_const HyperElasticity.ufl
-
-# 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)

--- a/example/Evolution/Evolution.m	Thu Sep 12 15:19:31 2013 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,41 +0,0 @@
-pkg load fem-fenics msh
-
-problem = 'Evolution';
-import_ufl_Problem (problem);
-
-mesh = Mesh (msh2m_structured_mesh (0:0.05:1, 0:0.05:1, 1, 1:4));
-
-V = FunctionSpace(problem, mesh);
-
-bc = DirichletBC (V, @(x,y) 1, 1:4);
-
-t = 0;
-dt = 0.01;
-T = 0.3;
-
-k = Expression ('dt', @(x,y) dt);
-
-u0 = Expression ('u0', @(x,y) 10*exp(-((x - 0.5)^2 + (y - 0.5)^2) / 0.02));
-
-a = BilinearForm (problem, V, k);
-L = LinearForm (problem, V, k, u0);
-
-A = assemble (a, bc);
-
-# solve the problem for each time step
-# We need to update only the lhs
-while t < T
-  t += dt;
-
-  # we can pass u0 to the lhs indifferently as a fem_coeff or
-  # as a fem_func
-  L = LinearForm (problem, V, k, u0);
-  b = assemble (L, bc);
-
-  u = A \ b;
- 
-  u0 = Function ('u0', V, u);
-
-  #press Q to make the plot continue
-  plot (u0);
-end

--- a/example/Evolution/Evolution.ufl	Thu Sep 12 15:19:31 2013 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,12 +0,0 @@
-element = FiniteElement("Lagrange", triangle, 1)
-
-u = TrialFunction(element)
-v = TestFunction(element)
-
-u0 = Coefficient(element)
-
-dt = Constant(triangle)
-
-eq = (1/dt)*(u-u0)*v*dx + inner(grad(u), grad(v))*dx
-a = rhs (eq)
-L = lhs (eq)

--- a/example/Mixed-Poisson/MixedPoisson.m	Thu Sep 12 15:19:31 2013 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,29 +0,0 @@
-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));
-
-V = FunctionSpace('MixedPoisson', mesh);
-
-f = Expression ('f', @(x,y) 10*exp(-((x - 0.5)^2 + (y - 0.5)^2) / 0.02));
-
-a = BilinearForm ('MixedPoisson', 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);

--- a/example/Mixed-Poisson/MixedPoisson.ufl	Thu Sep 12 15:19:31 2013 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,12 +0,0 @@
-BDM = FiniteElement("BDM", triangle, 1)
-DG  = FiniteElement("DG", triangle, 0)
-W = BDM * DG
-
-(sigma, u) = TrialFunctions(W)
-(tau, v)   = TestFunctions(W)
-
-CG = FiniteElement("CG", triangle, 1)
-f = Coefficient(CG)
-
-a = (dot(sigma, tau) + div(tau)*u + div(sigma)*v)*dx
-L = - f*v*dx

--- a/example/NavierStokes/NS.m	Thu Sep 12 15:19:31 2013 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,102 +0,0 @@
-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
-name = [tmpnam ".geo"];
-fid = fopen (name, "w");
-fputs (fid,"Point (1)  = {0, 0, 0, 0.1};\n");
-fputs (fid,"Point (2)  = {1, 0, 0, 0.1};\n");
-fputs (fid,"Point (3)  = {1, 0.5, 0, 0.1};\n");
-fputs (fid,"Point (4)  = {0.5, 0.5, 0, 0.1};\n");
-fputs (fid,"Point (5) = {0.5, 1, 0, 0.1};\n");
-fputs (fid,"Point (6) = {0, 1, 0,0.1};\n");
-
-fputs (fid,"Line (1)  = {5, 6};\n");
-fputs (fid,"Line (2) = {2, 3};\n");
-
-fputs (fid,"Line(3) = {6,1,2};\n");
-fputs (fid,"Line(4) = {5,4,3};\n");
-fputs (fid,"Line Loop(7) = {3,2,-4,1};\n");
-fputs (fid,"Plane Surface(8) = {7};\n");
-fclose (fid);
-msho = msh2m_gmsh (canonicalize_file_name (name)(1:end-4),"scale", 1,"clscale", .2);
-unlink (canonicalize_file_name (name));
-
-mesh = Mesh (msho);
-
-# Define function spaces (P2-P1). From ufl file
-V = FunctionSpace ('VelocityUpdate', mesh);
-Q = FunctionSpace ('PressureUpdate', mesh);
-
-# Set parameter values
-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.
-a1 = BilinearForm ('TentativeVelocity', V, k);
-
-# Pressure update.
-a2 = BilinearForm ('PressureUpdate', Q);
-
-# Velocity update
-a3 = BilinearForm ('VelocityUpdate', 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
-  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);
-
-  # Plot solution
-  plot (p1);
-  plot (u1);
-
-  # Save to file
-  save (p1, sprintf ("p_%3.3d", ++i));
-  delete (sprintf ("p_%3.3d.pvd", i));
-  save (u1, sprintf ("u_%3.3d", i));
-  delete (sprintf ("u_%3.3d.pvd", i));
-
-  # Move to next time step
-  u0 = u1;
-  t += dt
-
-end
-
-

--- a/example/NavierStokes/PressureUpdate.ufl	Thu Sep 12 15:19:31 2013 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,41 +0,0 @@
-# Copyright (C) 2010 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-30
-# Last changed: 2010-09-01
-#
-# The bilinear form a(u, v) and linear form L(v) for the pressure
-# update step in Chorin's method for the incompressible Navier-Stokes
-# equations.
-#
-# Compile this form with FFC: ffc -l dolfin PressureUpdate.ufl.
-
-# 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

--- a/example/NavierStokes/TentativeVelocity.ufl	Thu Sep 12 15:19:31 2013 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,45 +0,0 @@
-# Copyright (C) 2010 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-30
-# Last changed: 2011-06-30
-#
-# The bilinear form a(u, v) and linear form L(v) for the tentative
-# velocity step in Chorin's method for the incompressible
-# Navier-Stokes equations.
-#
-# Compile this form with FFC: ffc -l dolfin TentativeVelocity.ufl.
-
-# 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)

--- a/example/NavierStokes/VelocityUpdate.ufl	Thu Sep 12 15:19:31 2013 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,42 +0,0 @@
-# Copyright (C) 2010 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-30
-# Last changed: 2010-09-01
-#
-# The bilinear form a(u, v) and linear form L(v) for the velocity
-# update step in Chorin's method for the incompressible Navier-Stokes
-# equations.
-#
-# Compile this form with FFC: ffc -l dolfin TentativeVelocity.ufl.
-
-# 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)
-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

--- a/example/NavierStokes/lshape-domain.m	Thu Sep 12 15:19:31 2013 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,20 +0,0 @@
-name = [tmpnam ".geo"];
-fid = fopen (name, "w");
-fputs (fid,"Point (1)  = {0, 0, 0, 0.1};\n");
-fputs (fid,"Point (2)  = {1, 0, 0, 0.1};\n");
-fputs (fid,"Point (3)  = {1, 0.5, 0, 0.1};\n");
-fputs (fid,"Point (4)  = {0.5, 0.5, 0, 0.1};\n");
-fputs (fid,"Point (5) = {0.5, 1, 0, 0.1};\n");
-fputs (fid,"Point (6) = {0, 1, 0,0.1};\n");
-
-fputs (fid,"Line (1)  = {5, 6};\n");
-fputs (fid,"Line (2) = {2, 3};\n");
-
-fputs (fid,"Line(3) = {6,1,2};\n");
-fputs (fid,"Line(4) = {5,4,3};\n");
-fputs (fid,"Line Loop(7) = {3,2,-4,1};\n");
-fputs (fid,"Plane Surface(8) = {7};\n");
- fclose (fid);
- msho = msh2m_gmsh (canonicalize_file_name (name)(1:end-4),"scale", 1,"clscale", .2);
- trimesh (msho.t(1:3,:)', msho.p(1,:)', msho.p(2,:)');
- unlink (canonicalize_file_name (name));

--- a/example/Poisson/Poisson.m	Thu Sep 12 15:19:31 2013 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,28 +0,0 @@
-pkg load fem-fenics msh
-import_ufl_Problem ('Poisson')
-
-# Create mesh and define function space
-x = y = linspace (0, 1, 33);
-mesh = Mesh(msh2m_structured_mesh (x, y, 1, 1:4));
-
-V = FunctionSpace('Poisson', mesh);
-
-# Define boundary condition
-bc = DirichletBC(V, @(x, y) 0.0, [2;4]);
-
-f = Expression ('f', @(x,y) 10*exp(-((x - 0.5)^2 + (y - 0.5)^2) / 0.02));
-g = Expression ('g', @(x,y) sin (5.0 * x));
-
-a = BilinearForm ('Poisson', V);
-L = LinearForm ('Poisson', V, f, g);
-
-# Compute solution
-[A, b] = assemble_system (a, L, bc);
-sol = A \ b;
-u = Function ('u', V, sol);
-
-# Save solution in VTK format
-save (u, 'poisson')
-
-# Plot solution
-plot (u);

--- a/example/Poisson/Poisson.ufl	Thu Sep 12 15:19:31 2013 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,10 +0,0 @@
-# Copyright (C) 2005-2009 Anders Logg
-element = FiniteElement("Lagrange", triangle, 1)
-
-u = TrialFunction(element)
-v = TestFunction(element)
-f = Coefficient(element)
-g = Coefficient(element)
-
-a = inner(grad(u), grad(v))*dx
-L = f*v*dx + g*v*ds

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/inst/example/Biharmonic/Biharmonic.m	Thu Sep 12 15:22:19 2013 +0200
@@ -0,0 +1,24 @@
+pkg load fem-fenics msh
+
+problem = 'Biharmonic';
+import_ufl_Problem (problem);
+
+# Create mesh and define function space
+x = y = linspace (0, 1, 33);
+mesh = Mesh(msh2m_structured_mesh (x, y, 1, 1:4));
+
+V = FunctionSpace(problem, mesh);
+
+bc = DirichletBC (V, @(x,y) 0, 1:4);
+
+f = Expression ('f', @(x,y) 4.0*pi^4.*sin(pi.*x).*sin(pi.*y));
+g = Expression ('alpha', @(x,y) 8);
+
+a = BilinearForm (problem, V, g);
+L = LinearForm (problem, V, f);
+
+[A, b] = assemble_system (a, L, bc);
+u = A \ b;
+
+func = Function ('u', V, u);
+plot (func);

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/inst/example/Biharmonic/Biharmonic.ufl	Thu Sep 12 15:22:19 2013 +0200
@@ -0,0 +1,26 @@
+# Copyright (C) 2005-2009 Anders Logg
+
+# 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

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/inst/example/Elasticity/HyperElasticity.m	Thu Sep 12 15:22:19 2013 +0200
@@ -0,0 +1,57 @@
+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
+
+# 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
+
+u0 = assemble (L, bcs{:});
+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);

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/inst/example/Elasticity/HyperElasticity.ufl	Thu Sep 12 15:22:19 2013 +0200
@@ -0,0 +1,63 @@
+# Copyright (C) 2009-2010 Harish Narayanan and Garth N. Wells
+#
+# 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-09-29
+# Last changed: 2011-06-28
+#
+# The bilinear form a(du, v) and linear form L(v) for
+# a hyperelastic model.
+#
+# Compile this form with FFC: ffc -l dolfin -feliminate_zeros -fprecompute_basis_const -fprecompute_ip_const HyperElasticity.ufl
+
+# 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)

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/inst/example/Evolution/Evolution.m	Thu Sep 12 15:22:19 2013 +0200
@@ -0,0 +1,41 @@
+pkg load fem-fenics msh
+
+problem = 'Evolution';
+import_ufl_Problem (problem);
+
+mesh = Mesh (msh2m_structured_mesh (0:0.05:1, 0:0.05:1, 1, 1:4));
+
+V = FunctionSpace(problem, mesh);
+
+bc = DirichletBC (V, @(x,y) 1, 1:4);
+
+t = 0;
+dt = 0.01;
+T = 0.3;
+
+k = Expression ('dt', @(x,y) dt);
+
+u0 = Expression ('u0', @(x,y) 10*exp(-((x - 0.5)^2 + (y - 0.5)^2) / 0.02));
+
+a = BilinearForm (problem, V, k);
+L = LinearForm (problem, V, k, u0);
+
+A = assemble (a, bc);
+
+# solve the problem for each time step
+# We need to update only the lhs
+while t < T
+  t += dt;
+
+  # we can pass u0 to the lhs indifferently as a fem_coeff or
+  # as a fem_func
+  L = LinearForm (problem, V, k, u0);
+  b = assemble (L, bc);
+
+  u = A \ b;
+ 
+  u0 = Function ('u0', V, u);
+
+  #press Q to make the plot continue
+  plot (u0);
+end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/inst/example/Evolution/Evolution.ufl	Thu Sep 12 15:22:19 2013 +0200
@@ -0,0 +1,12 @@
+element = FiniteElement("Lagrange", triangle, 1)
+
+u = TrialFunction(element)
+v = TestFunction(element)
+
+u0 = Coefficient(element)
+
+dt = Constant(triangle)
+
+eq = (1/dt)*(u-u0)*v*dx + inner(grad(u), grad(v))*dx
+a = rhs (eq)
+L = lhs (eq)

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/inst/example/Mixed-Poisson/MixedPoisson.m	Thu Sep 12 15:22:19 2013 +0200
@@ -0,0 +1,29 @@
+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));
+
+V = FunctionSpace('MixedPoisson', mesh);
+
+f = Expression ('f', @(x,y) 10*exp(-((x - 0.5)^2 + (y - 0.5)^2) / 0.02));
+
+a = BilinearForm ('MixedPoisson', 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);

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/inst/example/Mixed-Poisson/MixedPoisson.ufl	Thu Sep 12 15:22:19 2013 +0200
@@ -0,0 +1,12 @@
+BDM = FiniteElement("BDM", triangle, 1)
+DG  = FiniteElement("DG", triangle, 0)
+W = BDM * DG
+
+(sigma, u) = TrialFunctions(W)
+(tau, v)   = TestFunctions(W)
+
+CG = FiniteElement("CG", triangle, 1)
+f = Coefficient(CG)
+
+a = (dot(sigma, tau) + div(tau)*u + div(sigma)*v)*dx
+L = - f*v*dx

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/inst/example/NavierStokes/NS.m	Thu Sep 12 15:22:19 2013 +0200
@@ -0,0 +1,102 @@
+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
+name = [tmpnam ".geo"];
+fid = fopen (name, "w");
+fputs (fid,"Point (1)  = {0, 0, 0, 0.1};\n");
+fputs (fid,"Point (2)  = {1, 0, 0, 0.1};\n");
+fputs (fid,"Point (3)  = {1, 0.5, 0, 0.1};\n");
+fputs (fid,"Point (4)  = {0.5, 0.5, 0, 0.1};\n");
+fputs (fid,"Point (5) = {0.5, 1, 0, 0.1};\n");
+fputs (fid,"Point (6) = {0, 1, 0,0.1};\n");
+
+fputs (fid,"Line (1)  = {5, 6};\n");
+fputs (fid,"Line (2) = {2, 3};\n");
+
+fputs (fid,"Line(3) = {6,1,2};\n");
+fputs (fid,"Line(4) = {5,4,3};\n");
+fputs (fid,"Line Loop(7) = {3,2,-4,1};\n");
+fputs (fid,"Plane Surface(8) = {7};\n");
+fclose (fid);
+msho = msh2m_gmsh (canonicalize_file_name (name)(1:end-4),"scale", 1,"clscale", .2);
+unlink (canonicalize_file_name (name));
+
+mesh = Mesh (msho);
+
+# Define function spaces (P2-P1). From ufl file
+V = FunctionSpace ('VelocityUpdate', mesh);
+Q = FunctionSpace ('PressureUpdate', mesh);
+
+# Set parameter values
+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.
+a1 = BilinearForm ('TentativeVelocity', V, k);
+
+# Pressure update.
+a2 = BilinearForm ('PressureUpdate', Q);
+
+# Velocity update
+a3 = BilinearForm ('VelocityUpdate', 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
+  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);
+
+  # Plot solution
+  plot (p1);
+  plot (u1);
+
+  # Save to file
+  save (p1, sprintf ("p_%3.3d", ++i));
+  delete (sprintf ("p_%3.3d.pvd", i));
+  save (u1, sprintf ("u_%3.3d", i));
+  delete (sprintf ("u_%3.3d.pvd", i));
+
+  # Move to next time step
+  u0 = u1;
+  t += dt
+
+end
+
+

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/inst/example/NavierStokes/PressureUpdate.ufl	Thu Sep 12 15:22:19 2013 +0200
@@ -0,0 +1,41 @@
+# Copyright (C) 2010 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-30
+# Last changed: 2010-09-01
+#
+# The bilinear form a(u, v) and linear form L(v) for the pressure
+# update step in Chorin's method for the incompressible Navier-Stokes
+# equations.
+#
+# Compile this form with FFC: ffc -l dolfin PressureUpdate.ufl.
+
+# 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

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/inst/example/NavierStokes/TentativeVelocity.ufl	Thu Sep 12 15:22:19 2013 +0200
@@ -0,0 +1,45 @@
+# Copyright (C) 2010 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-30
+# Last changed: 2011-06-30
+#
+# The bilinear form a(u, v) and linear form L(v) for the tentative
+# velocity step in Chorin's method for the incompressible
+# Navier-Stokes equations.
+#
+# Compile this form with FFC: ffc -l dolfin TentativeVelocity.ufl.
+
+# 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)

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/inst/example/NavierStokes/VelocityUpdate.ufl	Thu Sep 12 15:22:19 2013 +0200
@@ -0,0 +1,42 @@
+# Copyright (C) 2010 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-30
+# Last changed: 2010-09-01
+#
+# The bilinear form a(u, v) and linear form L(v) for the velocity
+# update step in Chorin's method for the incompressible Navier-Stokes
+# equations.
+#
+# Compile this form with FFC: ffc -l dolfin TentativeVelocity.ufl.
+
+# 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)
+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

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/inst/example/NavierStokes/lshape-domain.m	Thu Sep 12 15:22:19 2013 +0200
@@ -0,0 +1,20 @@
+name = [tmpnam ".geo"];
+fid = fopen (name, "w");
+fputs (fid,"Point (1)  = {0, 0, 0, 0.1};\n");
+fputs (fid,"Point (2)  = {1, 0, 0, 0.1};\n");
+fputs (fid,"Point (3)  = {1, 0.5, 0, 0.1};\n");
+fputs (fid,"Point (4)  = {0.5, 0.5, 0, 0.1};\n");
+fputs (fid,"Point (5) = {0.5, 1, 0, 0.1};\n");
+fputs (fid,"Point (6) = {0, 1, 0,0.1};\n");
+
+fputs (fid,"Line (1)  = {5, 6};\n");
+fputs (fid,"Line (2) = {2, 3};\n");
+
+fputs (fid,"Line(3) = {6,1,2};\n");
+fputs (fid,"Line(4) = {5,4,3};\n");
+fputs (fid,"Line Loop(7) = {3,2,-4,1};\n");
+fputs (fid,"Plane Surface(8) = {7};\n");
+ fclose (fid);
+ msho = msh2m_gmsh (canonicalize_file_name (name)(1:end-4),"scale", 1,"clscale", .2);
+ trimesh (msho.t(1:3,:)', msho.p(1,:)', msho.p(2,:)');
+ unlink (canonicalize_file_name (name));

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/inst/example/Poisson/Poisson.m	Thu Sep 12 15:22:19 2013 +0200
@@ -0,0 +1,28 @@
+pkg load fem-fenics msh
+import_ufl_Problem ('Poisson')
+
+# Create mesh and define function space
+x = y = linspace (0, 1, 33);
+mesh = Mesh(msh2m_structured_mesh (x, y, 1, 1:4));
+
+V = FunctionSpace('Poisson', mesh);
+
+# Define boundary condition
+bc = DirichletBC(V, @(x, y) 0.0, [2;4]);
+
+f = Expression ('f', @(x,y) 10*exp(-((x - 0.5)^2 + (y - 0.5)^2) / 0.02));
+g = Expression ('g', @(x,y) sin (5.0 * x));
+
+a = BilinearForm ('Poisson', V);
+L = LinearForm ('Poisson', V, f, g);
+
+# Compute solution
+[A, b] = assemble_system (a, L, bc);
+sol = A \ b;
+u = Function ('u', V, sol);
+
+# Save solution in VTK format
+save (u, 'poisson')
+
+# Plot solution
+plot (u);

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/inst/example/Poisson/Poisson.ufl	Thu Sep 12 15:22:19 2013 +0200
@@ -0,0 +1,10 @@
+# Copyright (C) 2005-2009 Anders Logg
+element = FiniteElement("Lagrange", triangle, 1)
+
+u = TrialFunction(element)
+v = TestFunction(element)
+f = Coefficient(element)
+g = Coefficient(element)
+
+a = inner(grad(u), grad(v))*dx
+L = f*v*dx + g*v*ds