% Written by Vesa Kaarnioja, 2023
% Demonstration on how to estimate the FE error for two FE solutions computed using two different FE meshes

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; % store FE nodes of dense mesh for later use
mass_dense = mass; % store mass matrix of dense mesh for later use

% Find FE solution on a coarse mesh
[mass,grad,nodes,element,interior,centers,ncoord,nelem] = FEMdata(4);
nodes_coarse = nodes;
A = UpdateStiffness(grad,ones(length(centers),1));
rhs = mass(interior,:) * f(nodes);
sol_coarse = zeros(ncoord,1);
tmp = A(interior,interior)\rhs;
sol_coarse(interior) = tmp;

% Interpolate the coarse FE results onto the dense FE mesh
sol_interp = griddata(nodes(:,1),nodes(:,2),sol_coarse,nodes_dense(:,1),nodes_dense(:,2),'linear');

% Compute the L2 error between the two FE solutions
sqrt((sol_dense-sol_interp)'*mass_dense*(sol_dense-sol_interp))
