# 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);