--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test/Makefile	Fri Jul 12 13:57:20 2013 +0200
@@ -0,0 +1,29 @@
+MKOCTFILE ?= mkoctfile
+
+OCTFILES=test_expr.oct
+
+
+CPPFLAGS += -DHAVE_DOLFIN_H
+LIBS += -ldolfin
+
+all: $(OCTFILES)
+
+test_expr.oct: test_expr.o
+	$(MKOCTFILE) $(CPPFLAGS) -s ../src/fem_init_env.o test_expr.o -o $@ $(LDFLAGS) $(LIBS)
+
+test_expr.o:  test_expr.cpp laplace
+	$(MKOCTFILE) $(CPPFLAGS) -c test_expr.cpp $(LDFLAGS) -o $@ -I../src/
+
+laplace:  Laplace.ufl
+	ffc -l dolfin Laplace.ufl
+
+clean:
+	-rm -f *.o core octave-core *.oct *~ *.xml *.pvd *.vtu
+
+cleanall:
+	-rm -f *.o core octave-core *.oct *~ *.xml *.status *.log octave-workspace configure
+	-rm -r autom4te.cache
+	-rm -f Makefile
+
+
+

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test/test_expr.cpp	Fri Jul 12 13:57:20 2013 +0200
@@ -0,0 +1,76 @@
+#include <dolfin.h>
+#include "Laplace.h"
+#include "expression.h"
+#include "functionspace.h"
+
+
+class Source : public dolfin::Expression
+{
+  void eval(dolfin::Array<double>& values, const dolfin::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);
+  }
+};
+
+class dUdN : public dolfin::Expression
+{
+  void eval(dolfin::Array<double>& values, const dolfin::Array<double>& x) const
+  {
+    values[0] = sin(5*x[0]);
+  }
+};
+
+
+DEFUN_DLD (test_expr, args, , "test_expr: functionspace V, fcn_handle f, mesh_label label")
+{
+
+  int nargin = args.length ();
+  octave_value retval=0;
+
+  if (nargin < 3 || nargin > 3)
+    print_usage ();
+  else
+    {
+      if (args(0).type_id () == functionspace::static_type_id ())
+        {
+          const functionspace & fspo = static_cast<const functionspace&> (args(0).get_rep ());
+          octave_fcn_handle * fh = args(1).fcn_handle_value ();
+          Array<int> side = args(2).array_value ();
+
+          if (!error_state)
+            {
+              const dolfin::FunctionSpace V = fspo.get_fsp ();
+              dolfin::Mesh mesh = *(V.mesh());
+              expression ff (*fh);
+
+              //We use dolfin now, fem-fenics not yet available
+
+              std::vector<const dolfin::BoundaryCondition*> bcu;
+              Source f;
+              dUdN g;
+
+              for (octave_idx_type i = 0; i < side.length (); ++i)
+                bcu.push_back (new dolfin::DirichletBC (V, ff, side(i)));
+
+
+              Laplace::BilinearForm a (V, V);
+              Laplace::LinearForm L (V);
+              L.f = f;
+              L.g = g;
+
+              dolfin::Function u(V);
+              dolfin::solve(a == L, u, bcu);
+
+              dolfin::File file("fem-fenics.pvd");
+              file << u;
+
+              dolfin::plot(u);
+              dolfin::interactive();
+            }
+        }
+    }
+  return retval;
+
+}

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test/test_expr.m	Fri Jul 12 13:57:20 2013 +0200
@@ -0,0 +1,27 @@
+
+# We follow the dolfin example for the Poisson problem 
+# (http://fenicsproject.org/documentation/dolfin/1.2.0/cpp/demo/pde/poisson/cpp/documentation.html#index-0)
+# to 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
+
+
+pkg load msh
+addpath("../src/")
+
+# create a unit square mesh using msh: labels for the boundary sides are 1,2,3,4
+msho = msh2m_structured_mesh (0:0.05:1, 0:0.05:1, 1, 1:4);
+
+mshd = fem_init_mesh (msho);
+V  = fem_fs (mshd);
+
+f = @(x,y) 0;
+
+# test_expr take 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
+
+test_expr (V, f, [2,4])