pkg load msh
pkg load fem-fenics
fem_init_env ();

problem = "MixedPoisson";
fem_create_all (problem);
x = y = linspace (0, 1, 32);
msho = msh2m_structured_mesh (x, y, 1, 1:4);
mshd = Mesh (msho);

V = MixedPoisson_FunctionSpace (mshd);

bc1 = DirichletBC (V, @(x,y) [0; -sin(5.0*x); 0], 1);
bc2 = DirichletBC (V, @(x,y) [0;  sin(5.0*x); 0], 3);

f = Expression ('f', @(x,y) 10*exp(-((x - 0.5)^2 + (y - 0.5)^2) / 0.02));

a = MixedPoisson_BilinearForm (V);
L = MixedPoisson_LinearForm (V, f);

[A, b] = assemble_system (a, L, bc1, bc2);

u = A \ b;

func = Function ('u', V, u);
uu = Function ('u', func, 0);
sigma = Function ('sigma', func, 1);

fem_plot (sigma);
fem_plot (uu);