function [A,b] = MyFEMCat()
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
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]; % coordinates of boundary points
bdry_connect = [e1; e2; e3]; % indices of endpoints of boundary segments

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

% coordinates is a N-by-2 array with coordinated of the mesh points
% elements3 is a Ntriag-by-3 array of indices of the triangular elements

%%
% 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
[A,b] = fem2d(coordinates,elements3,neumann,dirichlet);
end

%% fem2d_2
function [A,b] = fem2d(coordinates,elements3,neumann,dirichlet)
FreeNodes = setdiff(1:size(coordinates,1),dirichlet); %mesh points with unknown values of u
A = sparse(size(coordinates,1),size(coordinates,1));
b = sparse(size(coordinates,1),1);

%% Assembly
%% The Stiffness matrix
for j = 1:size(elements3,1) % for all triangles    
  A(elements3(j,:),elements3(j,:)) = A(elements3(j,:),elements3(j,:)) ...
      + stima3(coordinates(elements3(j,:),:));
  % stima3 computes M = 0.5*|T_j|*G*G';
end

%% The Right-hand side, i.e., the load vector
% Volume Forces
for j = 1:size(elements3,1)
  b(elements3(j,:)) = 0;  % for the case where f = 0
end

% Neumann conditions
for j = 1 : size(neumann,1)
  b(neumann(j,:)) = b(neumann(j,:)) + norm(coordinates(neumann(j,1),:)- ...
      coordinates(neumann(j,2),:)) * myg(sum(coordinates(neumann(j,:),:))/2)/2;
end

% Dirichlet conditions 
u = sparse(size(coordinates,1),1);
u(dirichlet) = myu_d(coordinates(dirichlet,:));
b = b - A * u;

% Computation of the solution
u(FreeNodes) = A(FreeNodes,FreeNodes) \ b(FreeNodes);

% graphic representation
trisurf(elements3,coordinates(:,1),coordinates(:,2),full(u)','facecolor','interp')
hold on
axis ij
end

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

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

%%
function M = stima3(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