%                 Aysajan Abidin

%                 January 9th 2006 - Asta Herdis Hall
%
% Function:    L2_error = L2_ErrorEst(u_true,u_est,p,e,t)
%
% 
% Input  
%      u_true:    The true solution of the equation, which is the FEM
%                 solution on the fine mesh.
%       u_est:    The estimated solution of the equation.
%                 An Np-x-L matrix, where L is the number of estimated 
%                 solutions that we are examining
%      TAG:       Represents the spatial dimension. 
%        
%        If TAG = 1, then
%        p:       1 x n vector of node coordinates.
%              
%        t:       3 x m matrix of element node numbers t(1,el)=nn1
%                 and t(2,el)=nn2, 1 <= el <= m, pointing to 
%                 corresponding node coordinates p(nn1) and p(nn2),
%                 and subdomain reference tags t(3,el), useful
%                 if data are given by different expresssions in
%                 different parts of the domain.
%
%        e:       2 x 2 matrix of boundary node numbers and 
%                 reference tags, with node number e(1,1)=nn pointing
%                 to the left boundary node coordinate p(nn), and e(2,2)
%                 is a reference tag (see example below).       
%
%        If TAG = 2, then
%        p,e,t:   Information about the space mesh (p=points, e=edges,
%                 t=triangles). For further information try: help initmesh

% Output
%     L2_error:   A Np-x-L matrix, where Np is the number of points in the
%                 space mesh and L is the number of estimated solutions that
%                 we are examining (i.e. length(u_est)).  It contains the
%                 L2 norm of the difference between the true solution
%                 and the estimated solution, taken over the domain
%                 specified by the data p,e,t.

function L2_error = L2_ErrorEst(u_true,u_est,p,e,t,TAG)

    Np = length(p);      % Number of points in space mesh
    L = size(u_est,2);   % Number of estimated solutions
    if TAG==1
    [stiff,mass] = Assemble_stiff_mass(p,e,t,TAG);
    L2_error = zeros(1,L);   
        for pp = 1:L   % Loop over different estimated solutions
            e = u_est(:,pp)-u_true;
            L2_error(1,pp) = L2_error(1,pp) + e'*mass*e;
        end
    L2_error = sqrt(L2_error); % computes the L_2 error.

    else
    Nt = size(t,2);      % Number of triangles in space mesh
    
    L2_error = zeros(1,L);
    
    for i = 1:Nt              % Looping over triangles in mesh
        triInd = t(1:3,i);     % Indices of triangle corner nodes
        triCoord = p(:,triInd);         % xy-coordinates of triangle corner nodes
        y = triCoord(2,:);x = triCoord(1,:);
        dx = [y(2)-y(3);y(3)-y(1);y(1)-y(2)];
        triSize = 0.5 * abs(x(1)*dx(1) + x(2)*dx(2) + x(3)*dx(3));
        
        for pp = 1:L   % Loop over different estimated solutions
            L2_error(1,pp) = L2_error(1,pp) + triSize/3* sum((u_est(triInd,pp)-u_true(triInd)).^2);
        end
    end
    L2_error = sqrt(L2_error);
end


    


