pkg load msh
pkg load fem-fenics

fem_init_env ();

fem_create_all ('DS');
fem_create_functional ('Err_v');
fem_create_functional ('Err_eps');
fem_create_functional ('Err_p');

uex = @(x, y) [-2*pi*sin(pi*x)*sin(pi*x)*cos(pi*y)*sin(pi*y);...
                2*pi*sin(pi*x)*cos(pi*x)*sin(pi*y)*sin(pi*y)];

pex = @(x, y) sin(pi*x);

gradp = @(x, y) [pi*cos(pi*x); 0];

laplu = @(x, y) [-4*pi^3*(1 - 4*sin(pi*x)*sin(pi*x))*sin(pi*y)*cos(pi*y);...
                  4*pi^3*(1 - 4*sin(pi*y)*sin(pi*y))*sin(pi*x)*cos(pi*x)];

eps_vector = [1 2^-2 2^-4 2^-8 0]';
nel_vector   = [4 8 16 32 64]';
h_vector = 1./nel_vector;

for ii = 1:numel(nel_vector)

  nel = nel_vector(ii);
  x = y = linspace (0, 1, nel + 1);
  msho = msh2m_structured_mesh (x, y, 1, 1:4);
  mshd = Mesh (msho);

  v = Expression ('v', uex);
  p = Expression ('p', pex);

  V  = DS_FunctionSpace (mshd);
  V0 = SubSpace (V, 0);

  bc = DirichletBC (V0, @(x,y) [0, 0], 1:4);

  for jj = 1:numel(eps_vector)

    epsi = eps_vector(jj);
    ff = @(x, y) uex (x, y) - epsi * epsi * laplu (x, y) - gradp (x, y);
    f = Expression ('f', ff);
    g = Expression ('g', @(x,y) 0);
    ep = Expression ('ep', @(x, y) epsi);

    a = DS_BilinearForm (V, ep);
    L = DS_LinearForm (V, f, g);

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

    sol = A \ b;

    func = Function ('uh', V, sol);
    vh = Function ('vh', func, 0);
    ph = Function ('ph', func, 1);

    M = Err_v_Functional (mshd, v, vh);
    err_v_L2(ii,jj) = sqrt (assemble (M));

    const = fem_eval (ph, [0, 0]);
    pex = @(x, y) sin(pi*x) + const;
    p = Expression ('p', pex);
    M = Err_p_Functional (mshd, p, ph);
    err_p_L2(ii,jj) = sqrt (assemble (M));

    M = Err_eps_Functional (mshd, v, vh, ep);
    err_v_eps(ii,jj) = sqrt (assemble (M));

  end

end

err_v_L2
err_p_L2
err_v_eps
pause ()

rate_v_L2 = (log (err_v_L2 (1:end-1,:)) - log (err_v_L2 (2:end,:)))./((log (h_vector (1:end-1)) - log (h_vector (2:end))));
conv_v_L2 = mean (rate_v_L2)

rate_p_L2 = (log (err_p_L2 (1:end-1,:)) - log (err_p_L2 (2:end,:)))./((log (h_vector (1:end-1)) - log (h_vector (2:end))));
conv_p_L2 = mean (rate_p_L2)

rate_v_eps = (log (err_v_eps (1:end-1,:)) - log (err_v_eps (2:end,:)))./((log (h_vector (1:end-1)) - log (h_vector (2:end))));
conv_v_eps = mean (rate_v_eps)