function [ error ] = Errornorm2D(u_e, u_h, p, t)
%UNTITLED Computes the L^2 error
%   Takes the following input:
%     u_e  --  exact solution
%     u_h  --  nodal values of numerical solution
%     p    --  nodal coordinates
%     t    --  the triangles in the mesh

    % Reference element quadrature points
    q = [0.816847572980459, 0.091576213509771; ...
         0.091576213509771, 0.816847572980459; ...
         0.091576213509771, 0.091576213509771; ...
         0.108103018168070, 0.445948490915965; ...
         0.445948490915965, 0.108103018168070; ...
         0.445948490915965, 0.445948490915965]';
     
    w = zeros(6,1);
    w(1:3) = 0.109951743655322;
    w(4:6) = 0.223381589678011;
    w = w/2;
    
    m = length(w);
    
    % Tabulate reference basis
    basis = {@(x) 1-x(1)-x(2); @(x) x(1); @(x) x(2)};
    B = zeros(m, 3);
    for i=1:m
        for j=1:3
            B(i, j) = basis{j}(q(:, i));
        end
    end

    % Loop over triangles:
    error2 = 0;
    for k = 1:length(t)
        % Determine trinagle coordinates
        vertices = t(1:3, k);
        coordinates = p(:,vertices);
        
        % Obtain the affine map from reference simplex
        A = coordinates(:,2:end) - repmat(coordinates(:,1),1,2);
        b = coordinates(:,1);
        J = abs(det(A)); % Jacobian
        
        % Map reference quadrature points
        Q = repmat(b,1,6) + A * q;

        % evaluate the numerical solution
        u_h_k = B * u_h(vertices);
        
        % evaluate exact solution
        u_e_k = zeros(m,1);
        for j=1:m
            u_e_k(j) = u_e(Q(:,j));
        end
        
        % Sum local L^2 error contribution
        error2 = error2 + J * sum(w.*(u_e_k - u_h_k).^2); 
    end
    error = sqrt(error2);
end

