# We follow the dolfin example for the Biharmonic problem
# See (http://fenicsproject.org/documentation/dolfin/1.2.0/cpp/demo/pde/biharmonic/cpp/documentation.html#index-0)

# initialize the environment and create function related to the problem
# defined inside Biharmonic.ufl
pkg load msh
pkg load fem-fenics
fem_init_env ();
problem = 'Biharmonic';
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
x = y = linspace (0, 1, 32);
msho = msh2m_structured_mesh (x, y, 1, 1:4);
mshd = fem_init_mesh (msho);

V  = fem_fs_Biharmonic (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
f = @(x,y) 0;
bc = fem_bc (V, f, 1: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) 4.0*pi^4.*sin(pi.*x).*sin(pi.*y);
f = fem_coeff ('f', ff);

gg = @(x,y) 8;
alpha = fem_coeff ('alpha', 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 alpha;
A = fem_rhs_Biharmonic (V, alpha, bc);

b = fem_lhs_Biharmonic (V, f, bc);

u = A \ b;

func = fem_func ('u', V, u);

fem_plot (func);

fem_save (func, problem);