--- a/test/Makefile	Sat Jul 20 23:49:23 2013 +0200
+++ b/test/Makefile	Sat Jul 20 23:53:31 2013 +0200
@@ -1,6 +1,6 @@
 MKOCTFILE ?= mkoctfile
 
-OCTFILES=test_expr.oct test_bc.oct test_coeff.oct
+OCTFILES=test_expr.oct test_bc.oct test_coeff.oct test_rlhs.oct
 
 LIBS += -ldolfin
 
@@ -24,6 +24,12 @@
 test_coeff.o:  test_coeff.cpp Laplace.h
 	$(MKOCTFILE) -c test_coeff.cpp $(LDFLAGS) -o $@ -I../src/
 
+test_rlhs.oct: test_rlhs.o
+	$(MKOCTFILE) -s ../src/fem_init_env.o test_rlhs.o -o $@ $(LDFLAGS) $(LIBS)
+
+test_rlhs.o:  test_rlhs.cc Laplace.h
+	$(MKOCTFILE) -c test_rlhs.cc $(LDFLAGS) -o $@ -I../src/
+
 
 Laplace.h: Laplace.ufl
 	ffc -l dolfin Laplace.ufl

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test/test_rlhs.cc	Sat Jul 20 23:53:31 2013 +0200
@@ -0,0 +1,46 @@
+#include "Laplace.h"
+using namespace Laplace;
+
+#include <dolfin.h>
+#include "functionspace.h"
+#include "boundarycondition.h"
+#include "coefficient.h"
+
+DEFUN_DLD (test_rlhs, args, , "test_bc: functionspace V, boundarycondition bc")
+{
+
+  int nargin = args.length ();
+  octave_value retval=0;
+  
+  if (nargin < 2)
+    print_usage ();
+  else
+    {
+      if (args(0).type_id () == functionspace::static_type_id ())
+        {
+          const functionspace & fspo =
+            static_cast<const functionspace&> (args(0).get_rep ());
+          Array <double> myu = args(1).array_value ();
+
+          if (!error_state)
+            {
+              const boost::shared_ptr<const dolfin::FunctionSpace> & V = fspo.get_pfsp ();
+
+              dolfin::Vector du(myu.length ());
+              for (std::size_t i = 0; i < myu.length (); ++i)
+                du.setitem (i, myu(i));
+
+              boost::shared_ptr<dolfin::Vector> uu (new dolfin::Vector(du));
+              dolfin::Function u(V, uu);
+
+              dolfin::File file ("fem-fenics-rlhs.pvd");
+              file << u;
+
+              dolfin::plot (u);
+              dolfin::interactive ();
+
+            }
+        }
+    }
+  return retval;
+}

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test/test_rlhs.m	Sat Jul 20 23:53:31 2013 +0200
@@ -0,0 +1,55 @@
+# 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;
+
+test_rlhs (V, u);