--- a/example/Elasticity/HyperElasticity.m	Mon Sep 09 22:45:03 2013 +0200
+++ b/example/Elasticity/HyperElasticity.m	Mon Sep 09 23:03:26 2013 +0200
@@ -1,35 +1,57 @@
-clear all
-pkg load msh
-pkg load fem-fenics
+pkg load fem-fenics msh
+problem = 'HyperElasticity';
+import_ufl_Problem (problem);
 
-fem_init_env ();
-problem = 'HyperElasticity';
-fem_create_all (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)];
 
-x = y = z = linspace (0, 1, 16);
-msho = msh3m_structured_mesh (x, y, z, 1, 1:6);
-mshd = Mesh (msho);
-
-V  = HyperElasticity_FunctionSpace (mshd);
+#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);
 
-bc1 = DirichletBC (V, @(x,y,z) [0; 0; 0], 1);
+# 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]);
 
-bc2 = DirichletBC (V, @(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)], 2);
-
-B = Expression ('B', @(x,y,z) [0.0; -0.5; 0.0]);
-T = Expression ('T', @(x,y,z) [0.1; 0.0; 0.0]);
+# Set material parameters
 E = 10.0;
 nu = 0.3;
-mu = Expression ('mu', @(x,y,z) E./(2*(1+nu)));
-lmbda = Expression ('lmbda', @(x,y,z) E*nu./((1+nu)*(1-2*nu)));
+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]);
 
-L = HyperElasticity_LinearForm (V, mu, lmbda, B, T, u);
-u0 = assemble (L, bc1, bc2);
-fs = @(xx) f (xx, V, bc1, bc2, B, T, mu, lmbda)
+# 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');
 
-func = Function ('u', V, x);
-fem_plot (func);
+# Plot solution
+plot (func);

--- a/example/Elasticity/f.m	Mon Sep 09 22:45:03 2013 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,14 +0,0 @@
-function [y, jac] = f (xx, V, bc1, bc2, B, T, mu, lmbda)
-
-  u = Function ('u', V, xx);
-
-  a = HyperElasticity_BilinearForm (V, mu, lmbda, u);
-  L = HyperElasticity_LinearForm (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

--- a/example/Mixed-Poisson/MixedPoisson.m	Mon Sep 09 22:45:03 2013 +0200
+++ b/example/Mixed-Poisson/MixedPoisson.m	Mon Sep 09 23:03:26 2013 +0200
@@ -1,31 +1,29 @@
-pkg load msh
-pkg load fem-fenics
-
-fem_init_env ();
+pkg load fem-fenics msh
+import_ufl_Problem ('MixedPoisson')
 
-problem = "MixedPoisson";
-fem_create_all (problem);
-x = y = linspace (0, 1, 32);
-msho = msh2m_structured_mesh (x, y, 1, 1:4);
-mshd = Mesh (msho);
+# Create mesh
+x = y = linspace (0, 1, 33);
+mesh = Mesh(msh2m_structured_mesh (x, y, 1, 1:4));
 
-V = MixedPoisson_FunctionSpace (mshd);
-V0 = SubSpace (V, 0);
-bc1 = DirichletBC (V0, @(x,y) [0; -sin(5.0*x)], 1);
-bc2 = DirichletBC (V0, @(x,y) [0;  sin(5.0*x)], 3);
+V = FunctionSpace('MixedPoisson', mesh);
 
 f = Expression ('f', @(x,y) 10*exp(-((x - 0.5)^2 + (y - 0.5)^2) / 0.02));
 
-a = MixedPoisson_BilinearForm (V);
-L = MixedPoisson_LinearForm (V, f);
+a = BilinearForm ('MixedPoisson', V);
+L = LinearForm ('MixedPoisson', V, f);
 
-[A, b] = assemble_system (a, L, bc1, bc2);
-
-u = A \ b;
+# 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);
 
-func = Function ('solution', V, u);
-sigma = Function ('sigma', func, 0);
-uu = Function ('u', func, 1);
+# Compute solution
+[A, b] = assemble_system (a, L, bc1, bc2);
+sol = A \ b;
+func = Function ('func', V, sol);
 
-fem_plot (sigma);
-fem_plot (uu);
+sigma = Function ('sigma', func, 1);
+u = Function ('u', func, 2);
+
+# Plot solution
+plot (sigma);
+plot (u);

--- a/example/NavierStokes/NS.m	Mon Sep 09 22:45:03 2013 +0200
+++ b/example/NavierStokes/NS.m	Mon Sep 09 23:03:26 2013 +0200
@@ -1,69 +1,102 @@
-pkg load msh
-pkg load fem-fenics
+pkg load fem-fenics msh
+import_ufl_Problem ("TentativeVelocity");
+import_ufl_Problem ("VelocityUpdate");
+import_ufl_Problem ("PressureUpdate");
 
-fem_init_env ();
-
-fem_create_all ("TentativeVelocity");
-fem_create_all ("VelocityUpdate");
-fem_create_all ("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");
 
-# We use the fem pkg and the script above to generate the L-shape domain
-run lshape-domain.m;
-mshd = Mesh (msho);
-
-V = VelocityUpdate_FunctionSpace (mshd);
-Q = PressureUpdate_FunctionSpace (mshd);
+fputs (fid,"Line (1)  = {5, 6};\n");
+fputs (fid,"Line (2) = {2, 3};\n");
 
-# no slip BC inside
-bcnoslip = DirichletBC (V, @(x,y) [0;0], [3,4]);
-# constant pressure at the outflow
-bcout = DirichletBC (Q, @(x,y) 0, 2);
+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));
 
-t = 0;
+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.;
 
-k = Expression ('k', @(x,y) dt);
-f = Expression ('f', @(x,y) [0; 0]);
+# 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]);
-p0 = Expression ('p1', @(x,y) 0);
 
-a1 = TentativeVelocity_BilinearForm (V, k);
-A1 = assemble (a1, bcnoslip);
-a3 = VelocityUpdate_BilinearForm (V);
-A3 = assemble (a3, bcnoslip);
+# 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);
 
-# solve the problem for each time step
-# We use the Chorin-Temam algorithm
-i = 0;
+# 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
-  t += dt;
-  # time dependent pressure at the inflow
-  bcin = DirichletBC (Q, @(x,y) sin(3.0*t), 1);
+
+  # Update pressure boundary condition
+  inflow = DirichletBC (Q, @(x,y) sin(3.0*t), 1);
 
-  L1 = TentativeVelocity_LinearForm (V, k, u0, f);
-  b1 = assemble (L1, bcnoslip);
+  # Compute tentative velocity step
+  L1 = LinearForm ('TentativeVelocity', V, k, u0, f);
+  b1 = assemble (L1, noslip);
   utmp = A1 \ b1;
- 
   u1 = Function ('u1', V, utmp);
 
-
-  # we need to assemble the system at each time-step because of time dependent BC
-  a2 = PressureUpdate_BilinearForm (Q);
-  L2 = PressureUpdate_LinearForm (Q,u1, k);
-  [A2, b2] = assemble_system (a2, L2, bcin, bcout);
-
+  # Pressure correction
+  L2 = LinearForm ('PressureUpdate', Q, u1, k);
+  [A2, b2] = assemble_system (a2, L2, inflow, outflow);
   ptmp = A2 \ b2;
   p1 = Function ('p1', Q, ptmp);
 
-  L3 = VelocityUpdate_LinearForm (V, k, u1, p1);
-  b3 = assemble (L3, bcnoslip);
+  # Velocity correction
+  L3 = LinearForm ('VelocityUpdate', V, k, u1, p1);
+  b3 = assemble (L3, noslip);
+  ut = A3 \ b3;
+  u1 = Function ('u0', V, ut);
 
-  ut = A3 \ b3;
-  u0 = Function ('u0', V, ut);
+  # Plot solution
+  plot (p1);
+  plot (u1);
 
-  # We use the fem_save () function to store the solution in vtu format
-  name = sprintf ("u_%3.3d", ++i);
-  fem_save (u0, name);
-  delete ([name ".pvd"]);
+  # 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/Poisson.m	Mon Sep 09 22:45:03 2013 +0200
+++ b/example/Poisson.m	Mon Sep 09 23:03:26 2013 +0200
@@ -1,28 +1,28 @@
-pkg load msh
-pkg load fem-fenics
-
-fem_init_env ();
-problem = 'Poisson';
-fem_create_all (problem);
+pkg load fem-fenics msh
+import_ufl_Problem ('Poisson')
 
-x = y = linspace (0, 1, 32);
-msho = msh2m_structured_mesh (x, y, 1, 1:4);
-mshd = Mesh (msho);
+# Create mesh and define function space
+x = y = linspace (0, 1, 33);
+mesh = Mesh(msh2m_structured_mesh (x, y, 1, 1:4));
 
-V  = Poisson_FunctionSpace (mshd);
+V = FunctionSpace('Poisson', mesh);
 
-bc = DirichletBC (V, @(x,y) 0, [2, 4]);
+# 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 = Poisson_BilinearForm (V);
-L = Poisson_LinearForm (V, f, g);
-
-[A, b] = assemble_system (a, L, bc);
+a = BilinearForm ('Poisson', V);
+L = LinearForm ('Poisson', V, f, g);
 
-u = A \ b;
+# Compute solution
+[A, b] = assemble_system (a, L, bc);
+sol = A \ b;
+u = Function ('u', V, sol);
 
-func = Function ('u', V, u);
+# Save solution in VTK format
+save (u, 'poisson')
 
-fem_plot (func);
+# Plot solution
+plot (u);