% Written by Vesa Kaarnioja, 2022
function [mass,grad,nodes,element,interior,centers,ncoord,nelem] = FEMdata(level)
n1 = 2^level+1;
h=1/2^level;

% Topology: FE nodes, mesh elements, interior
[X,Y]=meshgrid(0:h:1);
nodes = [X(:),Y(:)];
element = []; interior = [];
for ii = 1:n1-1
    for jj = 1:n1-1
        element = [element;
            (jj-1)*n1+ii,jj*n1+ii,(jj-1)*n1+ii+1;
            jj*n1+ii,jj*n1+ii+1,(jj-1)*n1+ii+1];
        if ii<n1-1 && jj<n1-1
            interior = [interior;jj*n1+1+ii];
        end
    end
end
ncoord = size(nodes,1);
nelem = size(element,1);

% For convenience, let's also find the element centers
centers = [];
for ii = 1:nelem
    centers = [centers;mean(nodes(element(ii,:),:))];
end

% Sanity check
% If PDE toolbox is installed, you can run the following:
% figure;pdeplot(nodes',element')
% hold on;
% plot(centers(:,1),centers(:,2),'r+')
% hold off;
%
% If PDE toolbox is _not_ installed, you can use the following.
% figure;hold on;
% xlim([0,1]), ylim([0,1])
% for ii = 1:nelem
%     plot(polyshape(nodes(element(ii,:),:)));
% end
% plot(centers(:,1),centers(:,2),'r+')
% hold off;

% Mass matrix assembly
mass = sparse(ncoord,ncoord);
for k = 1:nelem
  ind = element(k,:);
  g = nodes(ind,:);
  detB = abs(det([g(2,:)-g(1,:); g(3,:)-g(1,:)]));
  local = [1/12,1/24,1/24;
           1/24,1/12,1/24;
           1/24,1/24,1/12]*detB;
  mass(ind,ind) = mass(ind,ind)+local;
end

% Stiffness "tensor" assembly
grad = sparse(ncoord*ncoord,nelem);
for k = 1:nelem
  ind = element(k,:);
  g = nodes(ind,:);
  Dt = [g(3,:)-g(2,:); g(1,:)-g(3,:); g(2,:)-g(1,:)];
  % Use shoelace formula to find area of the triangle
  triarea = abs(det([g(1,:);g(2,:)])+det([g(2,:);g(3,:)])+det([g(3,:);g(1,:)]))/2;
  local = Dt*Dt'/(4*triarea);
  dummy = sparse(ncoord,ncoord); dummy(ind,ind) = local;
  grad(:,k) = dummy(:);
end