% Finite elements program
% adapted from "Applied Numerical Anaylsis" by L. Faussett, Prentice Hall.
% Solves u_xx + u_yy = 0 
% with fixed boundary conditions 
% A list of vertices V and a list of triangles T must be in memory.
% x and y must be row vectors containing coordinates of the vertices.
% The following variable must also be in memory:
% n is the number of vertices
% m is the number of internal vertices
% p is the number of triangles
% z input values at the vertices, including boundary values.

% Plot the boundary conditions and initial values of nodes
figure
trimesh(T,x,y,z)

%Define basis function for each node
for j = 1:n                       % for each node (finite element)  
    for k = 1:p                   % for each triangle
        for i = 1:3               % for each vertex of the triangle
            aa(i,1) = 1;
            aa(i,2) = V(T(k,i),1);
            aa(i,3) = V(T(k,i),2);
            if T(k,i) == j        % if the vertex goes with the  element
                bb(i,1) = 1;      % set the value to one.
            else
                bb(i,1) = 0;      % otherwise set it to zero.
            end
        end
        abc = aa\bb;      % solve for coefficients of function j on T_k.
        A(j,k) = abc(1);
        B(j,k) = abc(2);
        C(j,k) = abc(3);
    end
end

% Calculate the areas of the triangles
H = ones(1,p);  
for i = 1:p
   H(i) = .5*abs(det([x(T(i,1)), x(T(i,2)), x(T(i,3)); y(T(i,1)), y(T(i,2)), y(T(i,3)); 1, 1, 1]));
end

% This section solves for the values at the interior vertices
for i = 1:m
    for j = 1:m
        s(1:p) = B(i,1:p).*B(j,1:p) + C(i,1:p).*C(j,1:p);
        AAA(i,j) = s*H';
    end
end
for i = 1:m
    for j = m+1:n
        k = j - m;
        s(1:p) = B(i,1:p).*B(j,1:p) + C(i,1:p).*C(j,1:p);
        G(i,k) = s*H';
    end
end
d(1:m) = -G*z(m+1:n)';
z(1:m) = AAA\d';


% Plot the solution
figure
trimesh(T,x,y,z)