% Poisson FEM based on MIT lectures by Qiqi Wang
% Export the mesh data [p,e,t] in the wokspace

% Compute the gradient of the basis functions
grad = zeros(size(t,2),3,2);

% loop over number of triangles
for k = 1: size(t,2)
    % Store the coordinates of the triangles
    p1 = p(:,t(1,k));
    p2 = p(:,t(2,k));
    p3 = p(:,t(3,k));

    % Compute the gradient
    % the coefficient matrix: differences of coordinates
    C1 = [(p1-p2)';(p1-p3)'];
    % The RHS. Note that the basis function has the value 1 at p1 (say)
    r1 = [1;1];
    grad(k,1,:) = C1  r1;
    % Do it for the other nodes in a similar fashion
    C2 = [(p2-p1)';(p2-p3)'];
    grad(k,2,:) = C2  r1;
    C3 = [(p3-p1)';(p3-p2)'];
    grad(k,3,:) = C3  r1;
end
% Define the coefficient matrix
C = zeros(size(p,2));
% Note: the size of A must be the number of nodes. Ideally...
% number of interior nodes.

% Compute the area of each finite elements (triangles)
% loop over number of triangles
A = zeros (size(t,2),1);
for k = 1:size(t,2)
    % Store the coordinates of the triangles
    p1 = p(:,t(1,k));
    p2 = p(:,t(2,k));
    p3 = p(:,t(3,k));
    % Area(k) = det (coefficient matrix of basis) /2
    A(k) = abs(det([(p1-p2)';(p1-p3)']))/2;
end

% Compute the stiffness matrix
for k = 1:size(t,2)
    cr = t(1:3,k)';
    for i = 1:3
        for j = 1:3
            C(cr(i) , cr(j)) = C(cr(i) , cr(j)) +...
            A(k) * dot(grad(k,i,:),grad(k,j,:));
        end
    end
end

% Compute the global load vector (RHS)
b =  zeros(size(p,2),1);
% Define the RHS function
f=1;
for i=1:size(t,2)
    xcell=num2cell(t(1:3,i));
    [i1,i2,i3] = deal(xcell{:});
    % Compute the centroid of the triangle
    centroid = (p(:,i1)+p(:,i2)+p(:,i3))/3;
    % The integral is appoximated by f(centroid)*area/3
    b(i1) = b(i1)+A(i)/3;
    b(i2) = b(i2)+A(i)/3;
    b(i3) = b(i3)+A(i)/3;
end

M=size(p,2);
% Define the solution matrix
sol = diag(zeros(M));
% Get the number of interior nodes
dof=setdiff(1:M,e);
% Get the solution (Finally !)
sol(dof) = C(dof,dof)b(dof);
% Plot the solution
figure
a=trisurf(t',p(1,:)',p(2,:)',sol);