%******************** Direct Finite Element Method (DFEM) ***************** 
%******* for Fractional Diffusion Equation With Constant Coefficient ******
%
% Associated with the paper:
% H. Wang, D. Yang, S. Zhu, SIAM J. Numer. Anal. 52 (2014) 1292-1310
%
%  $$    - D ( K(x) * Dx^(-beta) Du ) = f,  0 < x < 1,
%  $$    u(0) = 0, u(1) = 0.
%  $$    0 <= beta < 1,
%  The symmetric bilinear form a and linear functional "L" associated 
%  with this BVP are
%     a(u,v) = \int_0^1 K(x) Dx^(1-beta)u Dv dx  
%     L(v)   = \int_0^1 f(x) v(x) dx.
%
%  Test case:  u(0) = u(1) = 0;   
%  u_true(x) = x^2*(1-x)^2
%  K(x)= 1
%  f(x)=-2*x^beta*(12/gamma(3+beta)*x.^2-6/gamma(2+beta)*x+1/gamma(1+beta))
%--------------------------------------------------------------------------
%                     Shengfeng Zhu, 2013, ECNU
%--------------------------------------------------------------------------
clear;
disp('***** DFEM Case: Constant Coefficient and Smooth Solution ****** ');
beta = input(' Parameter Beta in (0,1) = ');
beta1 = beta + 1;                   % beta+1
GammaBeta = gamma(beta);            % Gamma(beta)
GammaBeta1 = GammaBeta*beta;        % Gamma(beta+1)
GammaBeta2 = GammaBeta1*(beta+1);   % Gamma(beta+2)

nlevels = input(' No. of grid levels = ');
nvec = zeros(nlevels,1);
Eh_vec = zeros(nlevels,1);
hvec = zeros(nlevels,1);

% True Solution of u:
true_u = @(x) x.^2.*(1-x).^2;
% True Solution of u:
true_ux = @(x) 2*x.*(x-1).^2 + 2*x.^2.*(x-1);
% Right Force function f
right_hand_side = @(x) -2*x.^beta./((2+beta)*(1+beta)*beta*...
    GammaBeta).*(12*x.^2 - 6*(2+beta)*x + (2+beta)*(1+beta));

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;                    %  No. of elements.
    h = 1/nelement;                    %  Mesh spacing.
    hvec(k) = 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.
    n = n + 2;                         %  No. of grid points.
    [u_true] = true_u(x);
    %
    %  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;
    %
    %  Equation N is the right boundary condition, U(1.0) = 0.0;
    %
    A(n,n) = 1.0;
    b(n) = 0.0;
    %
    % Coefficient Matrix Assembling
    %
    parameter1 = 1/(h^(1-beta)*GammaBeta2);
    parameter2 = 1/(h^(1-beta)*GammaBeta2);
    % Diagonal: j == i
    vector = (4-2^beta1)*ones((n-3)/2,1);
    center = parameter1+parameter2*(3-2^beta1);
    A(2:end-1,2:end-1) = diag([parameter1*vector;center;parameter2*vector]);
    % j <= i-2
    for i = 4: (n-1)/2
        for j = 2: i-2
            A(i,j) = parameter1*(4*(i-j+1)^beta1-6*(i-j)^beta1+4*(i-j-1)^beta1...
                -(i-j-2)^beta1-(i-j+2)^beta1);
        end
    end
    index = 2:((n+1)/2-2);
    A((n+1)/2,index) = parameter1*(((n+1)/2-index+1).^beta1-3*((n+1)/2-index).^beta1+...
        3*((n+1)/2-index-1).^beta1-((n+1)/2-index-2).^beta1)+...
        parameter2*(3*((n+1)/2-index+1).^beta1-3*((n+1)/2-index).^beta1+...
        ((n+1)/2-index-1).^beta1-((n+1)/2-index+2).^beta1);
    for i = (n+1)/2 + 1 : n-1
        for j = 2: i-2
            A(i,j) = parameter2*(4*(i-j+1)^beta1-6*(i-j)^beta1+4*(i-j-1)^beta1...
                -(i-j-2)^beta1-(i-j+2)^beta1);
        end
    end
    % First Column of A
    A((n+1)/2,1) = parameter1*(2*((n+1)/2-1-1)^beta1-((n+1)/2-1-2)^beta1-((n+1)/2-1)^beta1)+...
        parameter2*(-2*((n+1)/2-1)^beta1+((n+1)/2-1-1)^beta1+((n+1)/2-1+1)^beta1);
    A(2,1) = parameter1*(2^beta1-3);
    index = 3:(n-1)/2;
    A(index,1) = parameter1*((index-1+1).^beta1-3*(index-1).^beta1+...
        3*(index-1+1).^beta1-(index-1-2).^beta1);
    index = (n+1)/2+1:n-1;
    A(index,1) = parameter2*((index-1+1).^beta1-3*(index-1).^beta1+...
        3*(index-1+1).^beta1-(index-1-2).^beta1);
    % j == i-1
    for i = 3:(n-1)/2
        A(i,i-1) = parameter1*(4*2^beta1-3^beta1-6);
    end
    A((n+1)/2,(n-1)/2) = parameter1*(2^beta1-3) + parameter2*(3*2^beta1-3-3^beta1);
    for i = (n+1)/2+1:n-1
        A(i,i-1) = parameter2*(4*2^beta1-3^beta1-6);
    end
    % j == i+1
    for i = 2:(n-1)/2
        A(i,i+1) = -parameter1;
    end
    for i = (n+1)/2:n-1
        A(i,i+1) = -parameter2;
    end
    %
    % RHS: Force Vector Assembling ( Accurate Integration )
    %
    xm = x(2:end-1);
    xl = x(1:end-2);
    xr = x(3:end);
    %
    b = zeros(n,1);
    b(2:end-1) = 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)) );
    %
    % Solve the linear system.
    u_h = A \ b;
    %  Compute max absolute error at mesh points, and record step size.
    error = u_h - u_true;
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %
    l2_error = 0.0;
    seminorm_error = 0.0;
    %
    %  Quadrature definitions.
    %
    quad_num = 2;
    abscissa(1) = -0.577350269189625764509148780502;
    abscissa(2) = +0.577350269189625764509148780502;
    weight(1) = 1.0;
    weight(2) = 1.0;
    %
    %  Integrate over each interval.
    %
    for i = 1 : n - 1
        xl = x(i);
        xr = x(i+1);
        u_l = u_h(i);
        u_r = u_h(i+1);
        for q = 1 : quad_num
            xq = ( ( 1.0 - abscissa(q) ) * xl   ...
                + ( 1.0 + abscissa(q) ) * xr ) ...
                /   2.0;
            wq = weight(q) * ( xr - xl ) / 2.0;
            %
            %  Use the fact that U is a linear combination of piecewise linears.
            %
            uq = ( ( xr - xq      ) * u_l ...
                + (      xq - xl ) * u_r ) ...
                / ( xr      - xl );
            uxq = ( u_r - u_l ) / ( xr - xl );
            eq = true_u ( xq );
            exq = true_ux ( xq );
            l2_error = l2_error + wq * ( uq - eq )^2;
            seminorm_error = seminorm_error + wq * ( uxq - exq )^2;
        end
    end
    l2_error = sqrt ( l2_error );
    seminorm_error = sqrt ( seminorm_error );
    fprintf('Mesh Number: %d\n', nelement);
    fprintf('L2_error: %E\n', l2_error);
    %fprintf('seminorm_error: %d\n', seminorm_error);
    Eh_vec(k) = l2_error;
end

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