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


pkg load msh
pkg load fem-fenics

# initialize the environment and create function related to the problem
# defined inside Laplace.ufl
fem_init_env ();
problem = 'Laplace';
fem_create_fs (problem);
fem_create_rhs (problem);
fem_create_lhs (problem);

# 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_Laplace (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 Lapalce problem with
# the velues specified inside f and g;
A = fem_rhs_Laplace (V, bc);

b = fem_lhs_Laplace (V, f, g, bc);

u = A \ b;

# fem_func create a function object which can be plotted, saved or
#  used to set the coefficient of the lhs/rhs
func = fem_func ('u', V, u)

fem_plot (func);

fem_save (func, problem);