%% Solution 1
load week5

% Solve the eigenvalue problem.
[evecs,evals] = eigs(stiffness,mass,1,'smallestabs');
evecs = evecs/sqrt(evecs'*mass*evecs); % normalize the eigenfunction
sol = zeros(ncoord,1); % initialize the solution vector
sol(interior) = -evecs; % plug in the eigenfunction values solved in the interior of the domain

% Plot the results over the actual FE mesh (note the use of "element" in trisurf!)
figure;
trisurf(element,nodes(:,1),nodes(:,2),sol)

%% Solution 2
% Solve the same problem using the FE mesh constructed in task 3!
clear all;
level = 5;
h = 1/2^level;
n1 = 2^(level+1)+1;

% FE nodes
[X,Y] = meshgrid(0:h:2);
nodes = [X(:),Y(:)];
ind = nodes(:,1)<=1 | nodes(:,2)<=1;
nodes = nodes(ind,:);

% Elements
element = [];
for ii = 1:n1-1
    for jj = 1:(n1-1)/2+1
        if ii <= (n1-1)/2 || jj <= (n1-1)/2
            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];
        end
        if ii <= (n1-1)/2 && jj <= (n1-3)/2
            element = [element;
                (jj-1)*(n1+1)/2+[n1*(n1+1)/2+ii,n1*(n1+1)/2+ii+(n1-1)/2+1,n1*(n1+1)/2+ii+1];
                (jj-1)*(n1+1)/2+[n1*(n1+1)/2+ii+1,(n1+1)*(n1+1)/2+ii,(n1+1)*(n1+1)/2+ii+1]];
        end
    end
end
ncoord = size(nodes,1);
nelem = size(element,1);
boundary1 = find(sum(nodes == 2,2)>0);
boundary2 = find(sum(nodes == 0,2)>0);
boundary3 = find((nodes(:,1)==1).*nodes(:,2)>=1);
boundary4 = find((nodes(:,1)>1).*nodes(:,2)==1);
boundary = [boundary1;boundary2;boundary3;boundary4];
interior = setdiff((1:ncoord)',boundary);
centers = [];
for ii = 1:nelem
    centers = [centers;mean(nodes(element(ii,:),:))];
end

% Note that the measurement geometry does not affect
% the assembly of the FE matrices!

% 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
mass = mass(interior,interior);
% 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
stiffness = reshape(grad*ones(nelem,1),ncoord,ncoord);
stiffness = stiffness(interior,interior);

[evecs,evals] = eigs(stiffness,mass,1,'smallestabs');
evecs = evecs/sqrt(evecs'*mass*evecs);
sol = zeros(ncoord,1);
sol(interior) = -evecs;
figure;
trisurf(element,nodes(:,1),nodes(:,2),sol)