%% %%%%%%%%%% Indirect Finite Element Method ( IFEM ) %%%%%%%%%%%%
% Fractional Diffusion Equation with Constant Coefficient & Smooth Solution 
%
% Associated with the paper:
% H. Wang, D. Yang, S. Zhu, SIAM J. Numer. Anal. 52 (2014) 1292-1310
%
%  $$    - D ( K(x) * D^(-beta)D u ) = f,  0 =< beta < 1;
%  $$    u(0) = 0, u(1) = 0.
% Constant Coefficient Case: K(x) = 1
% Smooth True Solution: u(x) = x^2*(1-x)^2
%****************** Accurately assemble RHS force vector ******************
%--------------------------------------------------------------------------
%                     Shengfeng Zhu, 2013, ECNU
%--------------------------------------------------------------------------
clear;
beta = input(' Parameter Beta in [0,1] = ');
%  "Anonymous functions" used by MATLAB to define a, f, and u_true.

diffusion_coeff = @(x) x./x; % constant coefficient: a(x) = 1
right_hand_side = @(x) -2*x.^beta./((2+beta)*(1+beta)*beta*...   % f
    gamma(beta)).*(12*x.^2 - 6*(2+beta)*x + (2+beta)*(1+beta));
% True Solution of u:
true_u = @(x) x.^2.*(1-x).^2;
% True Solution of DerBeta_w2: (w2 = x)
true_DerBeta_w2 = @(x) 1/gamma(2-beta)*x.^(1-beta);
% True Solution of w1;
true_w1 = @(x) 2/gamma(beta+5)*(12*x.^(4+beta)-6*(4+beta)*x.^(3+beta)...
    +(4+beta)*(3+beta)*x.^(2+beta)-(beta^2+beta)*x);
%
disp(' ************* Indirect Finite Element Method ( IFEM ) ************* ')
disp(' for Fractional Diffusion Equation with Cons. Coeff. & Smooth Sol.')
fprintf(' Parameter Beta: %g\n',beta);
disp(' Accurately assemble RHS force vector.');
nlevels = input(' No. of grid levels = ');
%
nvec = zeros(nlevels,1);
w1Err_vec = zeros(nlevels,1);
Eh_vec = zeros(nlevels,1);
hvec = zeros(nlevels,1);

for k = 1:nlevels
    
    %  Set up finite element mesh for level k.
    
    n = 4 * 2^k - 1;                   %  No. of interior grid points.
    nelement = n+1;
    nnodes = n+2;
    nvec(k) = n;
    h = 1/(n+1);                        %  Mesh spacing.
    fprintf(' Mesh size: %f\n',h);
    x = [0:h:1]';                      %  Mesh points.
    xfine = [0:h/8:1]';                %  Finer grid for true solution.
    x_int = x(2:n+1);                  %  Interior nodes
    x_mid = [h/2; x_int+h/2];          %  Midpoints of subintervals.
    u_h = zeros(nelement+1,1);
    
    %  Evaluate forcing function f(x) and diffusion coefficient kappa
    %  at subinterval midpoints. Set up left and right boundary conditions.
    
    kappa_mid = diffusion_coeff(x_mid);
    w1_left = 0;
    w1_right = 0;
    
    %  Compute FEM load vector b. 
    xm = x(2:end-1);
    xl = x(1:end-2);
    xr = x(3:end);
    
    b = 2/(h*gamma(beta+3))*( 12/((3+beta)*(4+beta))*(2*xm.^(4+beta)-xl.^(4+beta)...
        -xr.^(4+beta))-6/(3+beta)*(2*xm.^(3+beta)-xl.^(3+beta)-xr.^(3+beta))...
        +(2*xm.^(2+beta)-xl.^(2+beta)-xr.^(2+beta)) );
    
    %  Compute FEM stiffness matrix A. First set up main diagonal and
    %  subdiagonal of A. Since A is symmetric, the sub- and superdiagonals
    %  are the same.
    
    Adiag = (kappa_mid(1:n) + kappa_mid(2:n+1)) / h;
    Asub = -kappa_mid(2:n) / h;
    Asuper = Asub;
    
    %  Store nonzero components of A using MATLAB's sparse matrix
    %  storage mode.
    
    A = spdiags([[Asub;0],Adiag,[0;Asub]],[-1 0 1],n,n);
    
    %  Then solve for the coefficients of u_h(x) and append BC's.
    
    w1_interior = A \ b;
    w1_h = [w1_left; w1_interior; w1_right];
    
    %  Compute error at mesh points, and record step size.
    
    w1_true = true_w1(x);
    w1Error = w1_true - w1_h;
    
    w1Err_vec(k) = compute_l2_error ( n, x, w1_h, true_w1 );
    
    hvec(k) = h;
    u = true_u(x);
    u_trueFine = true_u(xfine);
    
    %  Evaluate DxBeta_wf_h.
    
    DerBeta_w1h = zeros(nelement+1,1);
    DerBeta_w1h(2) = w1_h(2)/(gamma(2-beta)*h^beta);
    for i = 2:nelement
        temp1 = ( (i+1-(2:i)'+1).^(1-beta) - (i+1-(2:i)').^(1-beta) ).*w1_h(2:i);
        temp2 = ( - (i+1-(2:i)'+1).^(1-beta) + (i+1-(2:i)').^(1-beta) ).*w1_h(1:i-1);
        temp = sum(temp1) + sum(temp2);
        DerBeta_w1h(i+1) = 1/(gamma(2-beta)*h^beta)*(temp + (w1_h(i+1)-w1_h(i)));
    end
    
    w2 = x;
    DerBeta_w2 = true_DerBeta_w2(x);
    coeff_end = - DerBeta_w1h(end)/DerBeta_w2(end);
    u_h(2:end-1) = DerBeta_w1h(2:end-1) + coeff_end*DerBeta_w2(2:end-1);
    uError = u - u_h;
    uErrNorm = compute_frac_l2_error(beta,coeff_end,nnodes,x,w1_h,1);
    fprintf(' Error Norm: %E\n',uErrNorm);
    Eh_vec(k) = uErrNorm;
    
    %  Plot true and approximate solutions in figure 1.
    
    figure(1)
    plot(x,u_h,'.-r', xfine,u_trueFine,'-b')
    xlabel('x axis')
    ylabel('u(x)')
    title_str = ['True solution u and FEM solution for  h = ', num2str(h)];
    title(title_str)
    legend('FEM Solution','True Solution')
    
end

%  Plot error vector E_h vs. h on a log-log scale.
figure(2)
loglog(hvec,Eh_vec,'o', hvec,Eh_vec)
xlabel('Mesh spacing h')
%ylabel('E_h = max_i |u_{true}(x_i) - u_h(x_i)|')
ylabel('E_h = ||u - u_h||_{L^2}')
title('Log-log plot of solution error E_h vs. mesh spacing h')

%  Compute rate of convergence.
conv_rate = diff(log(Eh_vec)) ./ diff(log(hvec));
fprintf(' Convergence Order: %f\n',conv_rate);
