%                 Aysajan Abidin
%
% 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
%         H,S:    The same as in Solver_spde1.m
%
% Output
%     L2_error:   A 1-x-L vector, where 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.

function L2_error = L2_ErrorEst(u_true,u_est,H,S)
L = H./S;
l = size(u_est,2);   % Number of estimated solutions 
L2_error = zeros(1,l); 
p = 1/L(1):1/L(1):(1-1/L(1));
n = length(p); 
t = [1:(n-1); 2:n; ones(1,(n-1))];
e = [1, n; 1,  2];         
[stiff mass] = stiff_mass(p,e,t);  % Assemble stiffness and mass matrices 
for i = 1:l   % Loop over different estimated solutions
    clear u_e
    n = H/2^i - 1; 
    u_e = u_est(:,i); 
    u_e(2^i:2^i:end) = u_e(1:n);%Example: this makes u_e=[0, 1, 0, 2, 0, 3, 0]', 
                                %here only in this example we wrote other entries
                                %of u_e as 0, but it doen't affect the computation. 
                                % originally u_e=[1,2,3,0,0,0,0]'
    %%%find the values for the entries before the first nonzero entry.
    for j = 1 : 2^i-1
        u_e(j) = j/2^i * u_e(2^i);  
    end
    %%%finde the values for the entries with 0 between the first and last nonzero entry. 
    for j = 2 : L(i+1)-1
        h = 1;
        for jj = j*2^i+1 : (j+1)*2^i-1            
            u_e(jj) = h/2^i * abs(u_e((j-1)*2^i) - u_e(j*2^i));
            h = h+1;
        end
    end
    %%%finde the values for the entries after the last nonzero entry. 
    h = 1;
    for j = L(i+1)*2^i+1 : L(1)
        u_e(j) = h/2 * u_e(L(i)*2^i);
        h = h+1;
    end
    err = u_e-u_true;    
    L2_error(1,i) = L2_error(1,i) + err'*mass*err;     
end
L2_error = sqrt(L2_error); % Computes the L_2 error
