%******************** Direct Finite Element Method (DFEM) ***************** 
%* for Fractional Diffusion Equation with Piecewise 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,
%  $$    K(x) = K1 in (0,1/2), K(x) = K2 in (1/2,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;   
%              If x in (0,1/2]
%                u_true(x) = 2*x
%                K(x) = Kl ( 0 < kl < 2^(beta+1)-3 )
%                f(x) = -2*Kl*x.^(beta-1)/GammaBeta
%              If x in (1/2,1]
%                u_true(x) = 2*(1-x)
%                K(x) = 1
%                f(x) = -2*(x.^(beta-1)-2*(x-1/2).^(beta-1))/GammaBeta
%--------------------------------------------------------------------------
%                     Shengfeng Zhu, 2013, ECNU
%--------------------------------------------------------------------------
clear;
Kl = 0.1;%2^beta - 1.5;%
fprintf('Kl: %g\n', Kl);
K1 = Kl;
K2 = 1;
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)
disp(' ************* Direct Finite Element Method ( DFEM ) ************* ')
disp(' Fractional Diffusion Equation with Piecewise Cons. Coeff. & Linear Sol.')
nlevels = input('No. of grid levels = ');
nvec = zeros(nlevels,1);
Eh_vec = zeros(nlevels,1);
hvec = zeros(nlevels,1);
Exact_u = @(x) 2*x.*(x <= 0.5) + (2*(1-x)).*(x>0.5);
Exact_ux = @(x) 2*ones(length(x),1).*(x <= 0.5) + (-2*ones(length(x),1).*(x > 0.5));
    
 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] = Exact_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 = K1/(h^(1-beta)*GammaBeta2);
     parameter2 = K2/(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
     parameter = -2/(h*GammaBeta2);
     for i = 2: (n-1)/2
         b(i) = K1*parameter*( x(i+1)^(beta+1) + x(i-1)^(beta+1) + 2*beta*x(i)^(beta+1) ...
             -(beta+1)*x(i-1)*x(i)^beta - (beta+1)*x(i)^beta*x(i+1) );
     end
     %
     b((n+1)/2) = K1*parameter*((0.5-h)^(beta+1)-0.5^(beta+1)+h*(beta+1)*0.5^beta)...
         + parameter*(-0.5^beta*(beta+1)*h+(0.5+h)^(beta+1)-0.5^(beta+1)-2*h^(beta+1));
     %
     for i = (n+1)/2 + 1 : n-1
         b(i) = parameter*( x(i+1)^(beta+1) + x(i-1)^(beta+1) - 2*x(i)^(beta+1) ...
             +4*(x(i)-0.5)^(beta+1) -2*(x(i-1)-0.5)^(beta+1) -2*(x(i+1)-0.5)^(beta+1) );
     end
     %
     % 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 = Exact_u ( xq );
             exq = Exact_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);
   
 