--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test/test_coeff.cpp	Thu Jul 18 14:52:47 2013 +0200
@@ -0,0 +1,80 @@
+#include "Laplace.h"
+using namespace Laplace;
+
+#include <dolfin.h>
+#include "functionspace.h"
+#include "boundarycondition.h"
+#include "coefficient.h"
+
+DEFUN_DLD (test_coeff, 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 () &&
+          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 ();
+              LinearForm L (V);
+              std::size_t ncoef = L.num_coefficients ();
+              if (nargin < 2 + ncoef || nargin > 2 + ncoef)
+                error ("Wrong number of coefficient");
+              else
+                {
+                  for (octave_idx_type i = 0; i < ncoef; ++i)
+                    {
+                      const coefficient & cf =static_cast<const coefficient&> (args(i+2).get_rep ());
+                      std::size_t n = L.coefficient_number (cf.get_str ());
+                      const boost::shared_ptr<const expression> & pexp = cf.get_expr ();
+                      L.set_coefficient (n, pexp);
+                    }
+
+                  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])));
+
+                  //Now use dolfin, fem-fenics not yet available
+                  dolfin::Vector b;
+                  dolfin::assemble (b, L);
+
+                  for (std::size_t i = 0; i < bcu.size(); i++)
+                    bcu[i]->apply(b);
+
+                  BilinearForm a (V, V);
+                  dolfin::Matrix A;
+                  dolfin::assemble (A, a);
+
+                  for (std::size_t i = 0; i < bcu.size(); i++)
+                    bcu[i]->apply(A);
+
+                  dolfin::Function u(V);
+                  dolfin::solve(A, *u.vector(), b, "gmres", "default");
+
+                  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_coeff.m	Thu Jul 18 14:52:47 2013 +0200
@@ -0,0 +1,49 @@
+# 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
+# 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);
+
+f = @(x,y) 0;
+
+# 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
+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);
+
+# test_coeff takes as input the functionspace V, and the
+# boundarycondition bc and solve the Poisson problem with
+# the velues specified inside f and g;
+test_coeff (V, bc, f, g);