function MyFEMheat()
global xmean
c = imread('cat.png');
cc = sum(c,3);
h = contour(cc,[1 1]);
% extract contours of face and eyes
c1 = h(:,2:6:1682)';  % face
c2 = h(:,2884:6:3080)'; % eye
c3 = h(:,3082:6:3278)'; % another eye
xmax = max(c1(:,1));
xmin = min(c1(:,1));
ymax = max(c1(:,2));
ymin = min(c1(:,2));
rx = xmax - xmin;
cx = 0.5*(xmax + xmin);
ry = ymax - ymin;
cy = 0.5*(ymax + ymin);
c1(:,1) = (c1(:,1) - cx)/rx + 0.5;  
c1(:,2) = (c1(:,2) - cy)/ry + 0.5;  
c2(:,1) = (c2(:,1) - cx)/rx + 0.5;  
c2(:,2) = (c2(:,2) - cy)/ry + 0.5;  
c3(:,1) = (c3(:,1) - cx)/rx + 0.5;  
c3(:,2) = (c3(:,2) - cy)/ry + 0.5;  
xmean = mean(c1,1);
lc1 = length(c1); 
lc2 = length(c2);
lc3 = length(c3);
% connect boundary points
a1 = [1 : lc1]';
e1 = [a1, circshift(a1,[-1 0])];
a2 = [1 : lc2]';
e2 = [lc1 + a2, lc1 + circshift(a2,[-1 0])];
a3 = [1 : lc3]';
e3 = [lc1 + lc2 + a3, lc1 + lc2 + circshift(a3,[-1 0])];
bdry = [c1; c2; c3];
bdry_connect = [e1; e2; e3];

% make mesh
[coordinates, elements3] = mesh2d(bdry,bdry_connect);

% Find boundary points with homogeneous Neumann BCs
% The number of rows in neumann is the number of boundary intervals with
% Neumann BCs
% Each row of neumann contains indices of endpoints of the corresponding
% boundary interval
ind = zeros(lc1,1);
for i = 1 : lc1
    ind(i) = find(coordinates(:,1) == c1(i,1) & coordinates(:,2) == c1(i,2));
end
neumann = [ind,circshift(ind,[-1, 0])];

% Find boundary points with Dirichlet BCs
% dirichlet is a column vector of indices of points with Dirichlet BCs
dirichlet = zeros(lc2 + lc3,1);
ind = zeros(lc2,1);
for i = 1 : lc2
    dirichlet(i) = find(coordinates(:,1) == c2(i,1) & coordinates(:,2) == c2(i,2));
end
ind = zeros(lc3,1);
for i = 1 : lc3
    dirichlet(lc2 + i) = find(coordinates(:,1) == c3(i,1) & coordinates(:,2) == c3(i,2));
end

% call FEM
fem2d_heat(coordinates,elements3,neumann,dirichlet);
end


%%
function fem2d_heat(coordinates,elements3,neumann,dirichlet)
%FEM2D_HEAT finite element method for two-dimensional heat equation. 2 %Initialisation
FreeNodes=setdiff(1:size(coordinates,1),unique(dirichlet));
A = sparse(size(coordinates,1),size(coordinates,1));
B = sparse(size(coordinates,1),size(coordinates,1)); 
T = 1; dt = 0.01; N = T/dt;
U = zeros(size(coordinates,1),N+1);
% Assembly
for j = 1:size(elements3,1)
	A(elements3(j,:),elements3(j,:)) = A(elements3(j,:), ...
    elements3(j,:)) + stima3_2(coordinates(elements3(j,:),:));
end
for j = 1:size(elements3,1)
	B(elements3(j,:),elements3(j,:)) = B(elements3(j,:), ...
    elements3(j,:)) + det([1,1,1;coordinates(elements3(j,:),:)'])... 
    *[2,1,1;1,2,1;1,1,2]/24;
end
% Initial Condition
U(:,1) = IC(coordinates); 
% time steps
for n = 2:N+1
    b = sparse(size(coordinates,1),1);
    % Volume Forces
    for j = 1:size(elements3,1)
        b(elements3(j,:)) = b(elements3(j,:)) + ... 
            det([1,1,1; coordinates(elements3(j,:),:)']) * ... 
            dt*myf(sum(coordinates(elements3(j,:),:))/3,n*dt)/6;
    end
    % Neumann conditions
    for j = 1 : size(neumann,1)
       b(neumann(j,:)) = b(neumann(j,:)) + ...
         norm(coordinates(neumann(j,1),:)-coordinates(neumann(j,2),:))*...
         dt*myg(sum(coordinates(neumann(j,:),:))/2,n*dt)/2;
    end
    % previous timestep
    b=b+B*U(:,n-1);
    % Dirichlet conditions
    u = sparse(size(coordinates,1),1);
    u(unique(dirichlet)) = myu_d(coordinates(unique(dirichlet),:),n*dt);
    b=b-(dt*A+B)*u;
    % Computation of the solution
    u(FreeNodes) = (dt*A(FreeNodes,FreeNodes)+ ...
            B(FreeNodes,FreeNodes))\b(FreeNodes);
    U(:,n) = u;
end
for k = 1 : 6
    figure;
    t = 0.2*(k - 1);
    p = ceil(t/dt) + 1;
    u = U(:,p);
    trisurf(elements3,coordinates(:,1),coordinates(:,2),full(u)','facecolor','interp')
    title(sprintf('Time = %d\n',t),'Fontsize',24);
    axis ij
    colorbar
    view(2)
    set(gca,'Fontsize',24);
end
end
%%
function u0 = IC(x)
global xmean
u0 = 10*exp(-4*((x(:,1) - 1.5*xmean(1)).^2 + (x(:,2) - 0.5*xmean(2)).^2));
end

%%
function DirichletBoundaryValue = myu_d(x,t)
xmin = min(x(:,1));
xmax = max(x(:,1));
midx = 0.5*(xmin + xmax);
%DirichletBoundaryValue =  0.5 * (sign(x(:,1) - midx) + 1);
DirichletBoundaryValue =  zeros(size(x,1),1);
end

%%
function Stress = myg(x,t)
Stress = zeros(size(x,1),1);
end

%%
function M = stima3_2(vertices)
d = size(vertices,2);
G = [ones(1,d+1);vertices'] \ [zeros(1,d);eye(d)];
M = det([ones(1,d+1);vertices']) * G * G' / prod(1:d);
end

%%
function heatsource = myf(x,t)
heatsource = zeros(size(x,1),1);
end