--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test/test_all.m	Mon Jul 22 14:49:10 2013 +0200
@@ -0,0 +1,59 @@
+# We follow the dolfin example for the Poisson problem
+#    -div ( grad (u) ) = f      on omega
+#                    u = h      on gamma_d;
+#                du/dn = g      on gamma_n;
+# See (http://fenicsproject.org/documentation/dolfin/1.2.0/cpp/demo/pde/poisson/cpp/documentation.html#index-0)
+# we check if: 
+#   1) the classes created within fem-fenics
+#      like "mesh" and "functionspace" hold correctly the dolfin data
+#   2) the class "expression", which we derived from dolfin::Expression
+#      correctly sets up the value for the bc using a function_handle
+#   3) the class "boundarycondition", which handle a vecotr of pointer
+#      to dolfin::DirichletBC correctly stores the value for the bc
+
+
+pkg load msh
+addpath ("../src/")
+fem_init_env ();
+
+# create a unit square mesh using msh: labels for the boundary sides are 1,2,3,4
+# we can use only 2D mesh for the moment
+# if you want to try with a 3D mesh, you need to use tetrahedron instead of
+# triangle inside Laplace.ufl and recompile fem_fs.cpp
+msho = msh2m_structured_mesh (0:0.05:1, 0:0.05:1, 1, 1:4);
+mshd = fem_init_mesh (msho);
+
+V  = fem_fs (mshd);
+
+# fem_bc takes as input the functionspace V, a function handler f,
+# and the sides where we want to apply the condition
+# The value on each point of the boundary is computed using
+# the eval method available inside expression.h
+# if a side is not specified, Neumann conditions are applied
+# with g specified below
+f = @(x,y) 0;
+bc = fem_bc (V, f, [2, 4]);
+
+# fem_coeff takes as input a string and a function handler
+# and is used below to set the value of the coefficient of the rhs
+
+ff = @(x,y) 10*exp(-((x - 0.5).^2 + (y - 0.5).^2) / 0.02);
+f = fem_coeff ('f', ff);
+
+gg = @(x,y) sin (5.0 * x);
+g = fem_coeff ('g', gg);
+
+# fem_rhs and fem_lhs takes as input the functionspace V, and the
+# boundarycondition bc and solve the Poisson problem with
+# the velues specified inside f and g;
+A = fem_rhs (V, bc);
+
+b = fem_lhs (V, f, g, bc);
+
+u = A \ b;
+
+func = fem_func (V, u)
+
+fem_plot (func);
+
+fem_save (func, 'solution');