--- a/test/Makefile	Fri Jul 12 13:58:58 2013 +0200
+++ b/test/Makefile	Mon Jul 15 14:42:18 2013 +0200
@@ -1,29 +1,25 @@
 MKOCTFILE ?= mkoctfile
 
-OCTFILES=test_expr.oct
+OCTFILES=test_expr.oct test_bc.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)
+	$(MKOCTFILE) -s ../src/fem_init_env.o test_expr.o -o $@ $(LDFLAGS) $(LIBS)
+
+test_expr.o:  test_expr.cpp Laplace.h
+	$(MKOCTFILE) -c test_expr.cpp $(LDFLAGS) -o $@ -I../src/
 
-test_expr.o:  test_expr.cpp laplace
-	$(MKOCTFILE) $(CPPFLAGS) -c test_expr.cpp $(LDFLAGS) -o $@ -I../src/
+test_bc.oct: test_bc.o
+	$(MKOCTFILE) -s ../src/fem_init_env.o test_bc.o -o $@ $(LDFLAGS) $(LIBS)
 
-laplace:  Laplace.ufl
+test_bc.o:  test_bc.cpp Laplace.h
+	$(MKOCTFILE) -c test_bc.cpp $(LDFLAGS) -o $@ -I../src/
+
+Laplace.h: 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
-
-
-
+	-rm -f *.o core octave-core *.oct *~ *.xml *.pvd *.vtu Laplace.h

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test/test_bc.cpp	Mon Jul 15 14:42:18 2013 +0200
@@ -0,0 +1,77 @@
+#include <dolfin.h>
+#include "Laplace.h"
+#include "functionspace.h"
+#include "boundarycondition.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_bc, args, , "test_bc: functionspace V, boundarycondition bc")
+{
+
+  int nargin = args.length ();
+  octave_value retval=0;
+
+  if (nargin < 2 || nargin > 2)
+    print_usage ();
+  else
+    {
+      if (args(0).type_id () == functionspace::static_type_id () &&
+          args(1).type_id () == boundarycondition::static_type_id ())
+        {
+          const functionspace & fspo =
+            static_cast<const functionspace&> (args(0).get_rep ());
+          const boundarycondition & bc =
+            static_cast<const boundarycondition&> (args(1).get_rep ());
+
+          if (!error_state)
+            {
+              const dolfin::FunctionSpace V = fspo.get_fsp ();
+              dolfin::Mesh mesh = *(V.mesh());
+
+              const std::vector<boost::shared_ptr <const dolfin::DirichletBC> > &
+              pbc = bc.get_bc ();
+              std::vector<const dolfin::BoundaryCondition*> bcu;
+
+              for (octave_idx_type i = 0; i < pbc.size (); ++i)
+                bcu.push_back (& (* (pbc[i])));
+
+              //We use dolfin now, fem-fenics not yet available
+              Source f;
+              dUdN g;
+              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-bc.pvd");
+              file << u;
+
+              dolfin::plot (u);
+              dolfin::interactive ();
+            }
+        }
+    }
+  return retval;
+
+}

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test/test_bc.m	Mon Jul 15 14:42:18 2013 +0200
@@ -0,0 +1,37 @@
+# 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/")
+
+# 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;
+
+# fem_bc 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
+# if a side is not specified, Neumann conditions are applied
+# with g = sin(5*x)
+bc = fem_bc (V, f, [2,4]);
+
+# test_bc take as input the functionspace V, and the
+# boundarycondition bc and solve the Poisson problem with
+#     f = 10*exp(-(dx*dx + dy*dy) / 0.02);
+test_bc (V, bc);