% Model solution to task 1 of week 11
clear all; close all;

% Find FE solution on a dense mesh
[mass,grad,nodes,element,interior,centers,ncoord,nelem] = FEMdata(5);
A = UpdateStiffness(grad,ones(length(centers),1));
f = @(x) x(:,1);
rhs = mass(interior,:)*f(nodes);

sol_dense = zeros(ncoord,1);
tmp = A(interior,interior)\rhs;
sol_dense(interior) = tmp;
nodes_dense = nodes;
mass_dense = mass;

errors = [];
for i = 1:4
    [mass,grad,nodes,element,interior,centers,ncoord,nelem] = FEMdata(i);
    % Find FE solution on a coarse mesh
    A = UpdateStiffness(grad,ones(length(centers),1));
    rhs = mass(interior,:)*f(nodes);
    sol = zeros(ncoord,1);
    tmp = A(interior,interior)\rhs;
    sol(interior) = tmp;
    % Intepolate onto dense FE mesh
    sol_interp = griddata(nodes(:,1),nodes(:,2),sol,nodes_dense(:,1),nodes_dense(:,2),'linear');
    errors = [errors;sqrt((sol_dense-sol_interp)'*mass_dense*(sol_dense-sol_interp))];
end

% Least squares fit
h = (1/2).^(1:4)';
A = ones(4,2);
A(:,2) = log(h);
lsq = A\log(errors);
lsq(1) = exp(lsq(1));

% Visualize the results
figure;
loglog(h,errors,'r.','MarkerSize',18), hold on
loglog(h,lsq(1)*h.^lsq(2),'b','LineWidth',2), hold off
set(gca,'FontSize',14,'FontWeight','bold')
legend('FE errors',['slope: ',num2str(lsq(2))]);
title('Finite element error')
xlabel('h');
ylabel('L^2 error')