% %%%%%%%%%%%%%%% Indirect Finite Element Method ( IFEM ) %%%%%%%%%%%%%%%%
% Fractional Diffusion Equation with Constant Coeffi & NonSmooth Solution 
%
% Associated with the paper:
% H. Wang, D. Yang, S. Zhu, SIAM J. Numer. Anal. 52 (2014) 1292-1310
%
%  $$    - D ( a(x) * D^(-beta)D u ) = f,  0 =< beta < 1;
%  $$    u(0) = 0, u(1) = 0.
% Constant Coefficient Case: i.e. a(x) = 1;
% True Solutoin not very Smooth: u(x) = x^lamda - x with lamda < 1.5

% The example tested here comes from Ervin's paper 2005
%  "Anonymous functions" used by MATLAB to define a, f, and u_true.
%--------------------------------------------------------------------------
%                     Shengfeng Zhu, 2013, ECNU
%--------------------------------------------------------------------------
  clear;
  beta = input(' Parameter Beta in (0,1) = ');
  lamda = 1.1;
  diffusion_coeff = @(x) x./x; % constant coefficient: a(x) = 1
  right_hand_side = @(x) -gamma(lamda+1)*x.^(lamda+beta-2)./gamma(lamda+beta-1)...
                         + x.^(beta-1)./gamma(beta); % f(x)
  true_u = @(x) x.^lamda - x; % True Solution
  
  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.

    f_mid = right_hand_side(x_mid);
    kappa_mid = diffusion_coeff(x_mid);  
    w1_left = 0;
    w1_right = 0;

    %  Compute FEM load vector b. Apply midpoint quadrature to the 
    %  forcing function, and then append left and right BC's.
    
    b = h/2 * (f_mid(1:n) + f_mid(2:n+1));
    b(1) = b(1) + w1_left / h;          
    b(n) = b(n) + w1_right / h;
    
    %  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 max absolute error at mesh points, and record step size.
    
    %w1_true = true_w1(x);
    %tempw1 = w1_true;
    %w1_true = true_w1(xfine);  %  Evaluate true soln on fine grid.
    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 = 1/gamma(2-beta)*x.^(1-beta);
    u_h(2:end-1) = DerBeta_w1h(2:end-1) - DerBeta_w1h(end)/DerBeta_w2(end)*DerBeta_w2(2:end-1);
    uError = u - u_h;
    uErrNorm = norm(uError,inf);
    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);
