--- a/example/Biharmonic.m	Mon Aug 05 11:45:54 2013 +0200
+++ b/example/Biharmonic.m	Mon Aug 05 13:20:12 2013 +0200
@@ -7,7 +7,7 @@
 
 x = y = linspace (0, 1, 32);
 msho = msh2m_structured_mesh (x, y, 1, 1:4);
-mshd = fem_init_mesh (msho);
+mshd = Mesh (msho);
 
 V  = Biharmonic_FunctionSpace (mshd);
 

--- a/example/Evolution.m	Mon Aug 05 11:45:54 2013 +0200
+++ b/example/Evolution.m	Mon Aug 05 13:20:12 2013 +0200
@@ -7,21 +7,21 @@
 fem_create_all (problem);
 
 msho = msh2m_structured_mesh (0:0.05:1, 0:0.05:1, 1, 1:4);
-mshd = fem_init_mesh (msho);
+mshd = Mesh (msho);
 
-V  = fem_fs_Evolution (mshd);
+V  = Evolution_FunctionSpace (mshd);
 
-bc = fem_bc (V, @(x,y) 1, 1:4);
+bc = DirichletBC (V, @(x,y) 1, 1:4);
 
 t = 0;
 dt = 0.01;
 T = 0.3;
 
-k = fem_coeff ('dt', @(x,y) dt);
+k = Expression ('dt', @(x,y) dt);
 
-u0 = fem_coeff ('u0', @(x,y) 10*exp(-((x - 0.5)^2 + (y - 0.5)^2) / 0.02));
+u0 = Expression ('u0', @(x,y) 10*exp(-((x - 0.5)^2 + (y - 0.5)^2) / 0.02));
 
-A = fem_rhs_Evolution (V, bc, k);
+A = Evolution_BilinearForm (V, bc, k);
 
 # solve the problem for each time step
 # We need to update only the lhs
@@ -30,11 +30,11 @@
 
   # we can pass u0 to the lhs indifferently as a fem_coeff or
   # as a fem_func
-  b = fem_lhs_Evolution (V, k, u0, bc);
+  b = Evolution_LinearForm (V, k, u0, bc);
 
   u = A \ b;
  
-  u0 = fem_func ('u0', V, u);
+  u0 = Function ('u0', V, u);
 
   #press Q to make the plot continue
   fem_plot (u0);

--- a/example/MixedPoisson.m	Mon Aug 05 11:45:54 2013 +0200
+++ b/example/MixedPoisson.m	Mon Aug 05 13:20:12 2013 +0200
@@ -6,20 +6,23 @@
 fem_create_all (problem);
 x = y = linspace (0, 1, 32);
 msho = msh2m_structured_mesh (x, y, 1, 1:4);
-mshd = fem_init_mesh (msho);
+mshd = Mesh (msho);
 
-V = fem_fs_MixedPoisson (mshd);
+V = MixedPoisson_FunctionSpace (mshd);
 
-bc1 = fem_bc (V, @(x,y) [0; -sin(5.0*x); 0], 1);
-bc2 = fem_bc (V, @(x,y) [0;  sin(5.0*x); 0], 3);
+bc1 = DirichletBC (V, @(x,y) [0; -sin(5.0*x); 0], 1);
+bc2 = DirichletBC (V, @(x,y) [0;  sin(5.0*x); 0], 3);
 
-f = fem_coeff ('f', @(x,y) 10*exp(-((x - 0.5)^2 + (y - 0.5)^2) / 0.02));
+f = Expression ('f', @(x,y) 10*exp(-((x - 0.5)^2 + (y - 0.5)^2) / 0.02));
 
-A = fem_rhs_MixedPoisson (V, bc1, bc2);
-b = fem_lhs_MixedPoisson (V, f, bc1, bc2);
+A = MixedPoisson_BilinearForm (V, bc1, bc2);
+b = MixedPoisson_LinearForm (V, f, bc1, bc2);
 
 u = A \ b;
 
 func = fem_func ('u', V, u);
+uu = fem_func ('u', func, 0);
+sigma = fem_func ('sigma', func, 1);
 
-fem_plot (func);
+fem_plot (sigma);
+fem_plot (uu);

--- a/example/NavierStokes/NS.m	Mon Aug 05 11:45:54 2013 +0200
+++ b/example/NavierStokes/NS.m	Mon Aug 05 13:20:12 2013 +0200
@@ -8,46 +8,46 @@
 fem_create_all ("PressureUpdate");
 
 run lshape-domain.m;
-mshd = fem_init_mesh (msho);
+mshd = Mesh (msho);
 
-V = fem_fs_VelocityUpdate (mshd);
-Q = fem_fs_PressureUpdate (mshd);
+V = VelocityUpdate_FunctionSpace (mshd);
+Q = PressureUpdate_FunctionSpace (mshd);
 
-bcnoslip = fem_bc (V, @(x,y) [0;0], [3,4]);
-bcout = fem_bc (Q, @(x,y) 0, 2);
+bcnoslip = DirichletBC (V, @(x,y) [0;0], [3,4]);
+bcout = DirichletBC (Q, @(x,y) 0, 2);
 
 t = 0;
 dt = 0.01;
 T = 3.;
 
-k = fem_coeff ('k', @(x,y) dt);
-f = fem_coeff ('f', @(x,y) [0; 0]);
-u0 = fem_coeff ('u0', @(x,y) [0; 0]);
-p0 = fem_coeff ('p1', @(x,y) 0);
-A1 = fem_rhs_TentativeVelocity (V, k, bcnoslip);
-A3 = fem_rhs_VelocityUpdate (V, bcnoslip);
+k = Expression ('k', @(x,y) dt);
+f = Expression ('f', @(x,y) [0; 0]);
+u0 = Expression ('u0', @(x,y) [0; 0]);
+p0 = Expression ('p1', @(x,y) 0);
+A1 = TentativeVelocity_BilinearForm (V, k, bcnoslip);
+A3 = VelocityUpdate_BilinearForm (V, bcnoslip);
 # solve the problem for each time step
 # We need to update only the lhs
 while t < T
   t += dt;
-  bcin = fem_bc (Q, @(x,y) sin(3.0*t), 1);
+  bcin = DirichletBC (Q, @(x,y) sin(3.0*t), 1);
 
-  b1 = fem_lhs_TentativeVelocity (V, k, u0, f, bcnoslip);
+  b1 = TentativeVelocity_LinearForm (V, k, u0, f, bcnoslip);
 
   utmp = A1 \ b1;
  
-  u1 = fem_func ('u1', V, utmp);
+  u1 = Function ('u1', V, utmp);
 
-  A2 = fem_rhs_PressureUpdate (Q, bcin, bcout);
-  b2 = fem_lhs_PressureUpdate (Q, bcin, bcout, u1, k);
+  A2 = PressureUpdate_BilinearForm (Q, bcin, bcout);
+  b2 = PressureUpdate_LinearForm (Q, bcin, bcout, u1, k);
 
   ptmp = A2 \ b2;
-  p1 = fem_func ('p1', Q, ptmp);
+  p1 = Function ('p1', Q, ptmp);
 
-  b3 = fem_lhs_VelocityUpdate (V, k, u1, p1, bcnoslip);
+  b3 = VelocityUpdate_LinearForm (V, k, u1, p1, bcnoslip);
 
   ut = A3 \ b3;
-  u0 = fem_func ('u0', V, ut);
+  u0 = Function ('u0', V, ut);
 
   #press Q to make the plot continue
   fem_plot (u0);

--- a/example/Poisson.m	Mon Aug 05 11:45:54 2013 +0200
+++ b/example/Poisson.m	Mon Aug 05 13:20:12 2013 +0200
@@ -5,21 +5,19 @@
 problem = 'Poisson';
 fem_create_all (problem);
 
- 
-
 x = y = linspace (0, 1, 32);
 msho = msh2m_structured_mesh (x, y, 1, 1:4);
-mshd = fem_init_mesh (msho);
+mshd = Mesh (msho);
 
-V  = fem_fs_Poisson (mshd);
+V  = Poisson_FunctionSpace (mshd);
 
-bc = fem_bc (V, @(x,y) 0, [2, 4]);
+bc = DirichletBC (V, @(x,y) 0, [2, 4]);
 
-f = fem_coeff ('f', @(x,y) 10*exp(-((x - 0.5)^2 + (y - 0.5)^2) / 0.02));
-g = fem_coeff ('g', @(x,y) sin (5.0 * x));
+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 = fem_rhs_Poisson (V, bc);
-b = fem_lhs_Poisson (V, f, g, bc);
+A = Poisson_BilinearForm (V, bc);
+b = Poisson_LinearForm (V, f, g, bc);
 
 u = A \ b;
 
@@ -27,5 +25,3 @@
 func = fem_func ('u', V, u);
 
 fem_plot (func);
-
-fem_save (func, problem);