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