%***************** Indirect Finite Element Method ( IFEM ) ****************
%****  Fractional Diffusion Equation With Smoothly Varying Coefficient  ***
%
% 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.
% Test Case: 
%  $$ K(x) = 1./(1+x)
%  $$ u(x) = x.^2.*(1-x).^2
%--------------------------------------------------------------------------
%                     Shengfeng Zhu, 2013, ECNU
%--------------------------------------------------------------------------
clear;
beta = input(' Parameter Beta ( in (0,1) ) = ');
diffusion_coeff = @(x) 1./(1+x);
right_hand_side = @(x) -1./(gamma(beta+4)*(1+x).^2).*( (24*beta+48)*x.^(3+beta)...
    +(3+beta)*(12-12*beta)*x.^(2+beta)+(3+beta)*(2+beta)*(2*beta-12)*x.^(1+beta)...
    + (6+2*beta)*(2+beta)*(1+beta)*x.^(beta) );
true_u = @(x) x.^2.*(1-x).^2;

fprintf(' Parameter Beta: %g\n',beta);
disp(' Two-point Gauss formula to assemble RHS force vector.');
nlevels = input(' No. of grid levels = ');
nvec = 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;
    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.
    
    n = n + 2;  % Grid number
    %  Quadrature definitions.
    %
    quad_num = 2;
    abscissa(1) = -0.577350269189625764509148780502;
    abscissa(2) = +0.577350269189625764509148780502;
    weight(1) = 1.0;
    weight(2) = 1.0;
    w1_left = 0;
    w1_right = 0;
    %
    %  Make room for the matrix A and right hand side b.
    %
    A = zeros ( n, n );
    b = zeros ( n, 1 );
    %
    %  Equation 1 is the left boundary condition, U(0.0) = 0.0;
    %
    A(1,1) = 1.0;
    b(1) = 0.0;
    
    for i = 2 : n - 1
        %
        %  Get the left, right and middle coordinates.
        %
        xl = x(i-1);
        xm = x(i);
        xr = x(i+1);
        %
        %  Make temporary variables for A(I,I-1), A(I,I), A(I,I+1) and B(I).
        %
        al = 0.0;
        am = 0.0;
        ar = 0.0;
        bm = 0.0;
        %
        %  We approximate the integrals by using a weighted sum of
        %  the integrand values at quadrature points.
        %
        for q = 1 : quad_num
            %
            %  Integrate over the LEFT interval, between XL and XM
            xq = ( ( 1.0 - abscissa(q) ) * xl   ...
                + ( 1.0 + abscissa(q) ) * xm ) ...
                / 2.0;
            
            wq = weight(q) * ( xm - xl ) / 2.0;
            
            %      vl = ( xm - xq ) / ( xm - xl );
            vlp =     - 1.0  / ( xm - xl );
            
            vm = ( xq - xl ) / ( xm - xl );
            vmp = +1.0       / ( xm - xl );
            
            %     vr =  0.0;
            vrp = 0.0;
            
            axq = diffusion_coeff ( xq );
            fxq = right_hand_side ( xq );
            
            al = al + wq * ( axq * vlp * vmp );
            am = am + wq * ( axq * vmp * vmp );
            ar = ar + wq * ( axq * vrp * vmp );
            bm = bm + wq * ( fxq * vm );
            %
            %  Integrate over the RIGHT interval, between XM and XR, 
            xq = ( ( 1.0 - abscissa(q) ) * xm   ...
                + ( 1.0 + abscissa(q) ) * xr ) ...
                / 2.0;
            
            wq = weight(q) * ( xr - xm ) / 2.0;
            
            vl = 0.0;
            vlp = 0.0;
            
            vm = ( xr - xq ) / ( xr - xm );
            vmp =     - 1.0  / ( xr - xm );
            
            vr = ( xq - xm ) / ( xr - xm );
            vrp =  1.0       / ( xr - xm );
            
            axq = diffusion_coeff ( xq );
            fxq = right_hand_side ( xq );
            
            al = al + wq * ( axq * vlp * vmp );
            am = am + wq * ( axq * vmp * vmp );
            ar = ar + wq * ( axq * vrp * vmp );
            bm = bm + wq * ( fxq * vm );
            
        end
        A(i,i-1) = al;
        A(i,i)   = am;
        A(i,i+1) = ar;
        b(i)     = bm;
    end
    %
    %  Equation N is the right boundary condition, U(1.0) = 0.0;
    A(n,n) = 1.0;
    b(n) = 0.0;
    
    %  Then solve for the coefficients of u_h(x) and append BC's.
    w1_h = A \ b;
        
    %  Compute error at mesh points, and record step size.
    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 = 1/3*(x.^2+2*x);
    DerBeta_w2 = 2/(3*gamma(3-beta))*( (2-beta)*x.^(1-beta) + x.^(2-beta) );
    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,n,x,w1_h,2);
    fprintf(' Error Norm: %E\n',uErrNorm);
    Eh_vec(k) = uErrNorm;
    
    %  Plot true and approximate solutions in figure 1.
    figure(1)
    plot(x,u_h,'o-', xfine,u_trueFine,'--')
    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')
    disp(' Hit any key to continue.');
    pause
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)|')
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);
%************************************************************************************
