%% %%%%%%%%%% Indirect Finite Element Method ( IFEM ) %%%%%%%%%%%%
% ********* Inhomogeneous Dirichlet Fractional Diffusion Equation *********
% ********************* Piecewise Constant 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 < x < 1, 0 < beta < 1;
%  $$    u(0) = ul, u(1) = ur.
% True Solution: u(x) = 1 + x.^alpha
%  ul = 1; ur = 2

% Step 0: Homogenization of boundary condition leads to
% $$     - D ( K(x) * D^(-beta)D uf ) = f(x) + p(x),  0 < x < 1,
% $$     uf(0) = 0,  uf(1) = 0.
% Perturbated RHS: 
% p(x) = (ur-ul)*gamma(m+1)/gamma(m+beta)*(K'*x^(m+beta-1)+K*(m+beta-1)*x^(m+beta-2))
% Solution Composition: 
%  u = uf + ( ul*gl + ur*gr ) with gl = 1 - x^m, gr = x^m
%
% Step 1:
%  Linear FEM to solve diffusion equation with perturbated rhs:
% $$     - D (K(x) Dwf) = f(x) + p(x),  0 < x < 1,
% $$     wf(0) = 0,  wf(1) = 0.
%
% Step 2:
% $$     - D (K(x) Dwb) = 0,  0 < x < 1,
% $$     wb(0) = 0, wb(1) = 1.
%  uf(x) = D^beta wf(x) + D^beta wf(1) / D^beta wb(1) * D^beta wb(x)
%
% Test Case: Piecewise Constant Difussion Coefficient 
%  K = K1 in (0,1/2); K = K2 in (1/2,1)
%  wb(x) = x;   D^beta wb = 1/gamma(2-beta)*x^(1-beta)
%  f(x) = - gamma(alpha+1)/gamma(alpha+beta-1)*x.^(alpha+beta-2)
%  p(x) = gamma(m+1)/gamma(m+beta-1)*x.^(m+beta-2)

%****************** Accurately assemble RHS force vector ******************
%--------------------------------------------------------------------------
%                     Shengfeng Zhu, 2013, ECNU
%--------------------------------------------------------------------------
  clear; 
  alpha = 2;
  beta = input(' Parameter Beta in (0,1] = ');
  m = alpha - beta/2;
  K1 = 1/2;
  K2 = 1;
  
  %  "Anonymous functions" used by MATLAB to define a, f, and u_true.

  diffusion_coeff = @(x)  K1*(x<=0.5) + K2*(x>0.5); % Piecewise constant coefficient
  right_hand_side = @(x) (-gamma(alpha+1)/gamma(alpha+beta-1)*x.^(alpha+beta-2)...
   + gamma(m+1)/gamma(m+beta-1)*x.^(m+beta-2)).*(K1*(x<=0.5) + K2*(x>0.5)); % f + p
% True Solution of u:
  true_u = @(x) x.^alpha + 1; 
% ub for Homogenization:
  true_ub = @(x) x.^m + 1;
% True Solution of DxBeta_wb: (wb = x) 
  true_DxBeta_wb = @(x) 2*K2/((K1+K2)*gamma(2-beta))*x.^(1-beta).*(x<=0.5)...
      + 2/((K1+K2)*gamma(2-beta))*(K2*x.^(1-beta)+(K1-K2)*(x-0.5).^(1-beta)).*(x>0.5); 
% True Solution of wf;
  true_wf = @(x) (gamma(alpha+1)/gamma(alpha+beta+1)*x.^(alpha+beta)...
      - gamma(m+1)/gamma(m+beta+1)*x.^(m+beta)+(-gamma(alpha+1)/gamma(alpha+beta+1)...
      +gamma(m+1)/gamma(m+beta+1)+(K2-K1)/(K2+K1)*(gamma(alpha+1)/gamma(alpha+beta)...
      *(0.5^(alpha+beta-1)-1/(alpha+beta))-gamma(m+1)/gamma(m+beta)*(0.5^(m+beta-1)-1/(m+beta))))*x).*(x<=0.5)...
      +(gamma(alpha+1)/gamma(alpha+beta+1)*x.^(alpha+beta)...
      - gamma(m+1)/gamma(m+beta+1)*x.^(m+beta)+(-gamma(alpha+1)/gamma(alpha+beta+1)...
      +gamma(m+1)/gamma(m+beta+1))*x+(K2-K1)/(K2+K1)*(gamma(alpha+1)/gamma(alpha+beta)...
      *(0.5^(alpha+beta-1)-1/(alpha+beta))-gamma(m+1)/gamma(m+beta)*(0.5^(m+beta-1)-1/(m+beta)))*(1-x)).*(x>0.5);
% True Solution of Dwf;
  true_Dwf = @(x) (gamma(alpha+1)/gamma(alpha+beta)*x.^(alpha+beta-1)...
      - gamma(m+1)/gamma(m+beta)*x.^(m+beta-1)+(-gamma(alpha+1)/gamma(alpha+beta+1)...
      +gamma(m+1)/gamma(m+beta+1)+(K2-K1)/(K2+K1)*(gamma(alpha+1)/gamma(alpha+beta)...
      *(0.5^(alpha+beta-1)-1/(alpha+beta))-gamma(m+1)/gamma(m+beta)*(0.5^(m+beta-1)-1/(m+beta))))).*(x<=0.5)...
      +(gamma(alpha+1)/gamma(alpha+beta+1)*x.^(alpha+beta)...
      - gamma(m+1)/gamma(m+beta+1)*x.^(m+beta)+(-gamma(alpha+1)/gamma(alpha+beta+1)...
      +gamma(m+1)/gamma(m+beta+1))*x+(K2-K1)/(K2+K1)*(gamma(alpha+1)/gamma(alpha+beta)...
      *(0.5^(alpha+beta-1)-1/(alpha+beta))-gamma(m+1)/gamma(m+beta)*(0.5^(m+beta-1)-1/(m+beta)))*(1-x)).*(x>0.5);
  disp(' Accurately assemble RHS force vector.'); 
  disp('**** PARAMETERS ****:'); 

  fprintf(' Beta = %g.\n',beta);
  fprintf(' K1 = %g;   ',K1);
  fprintf(' K2 = %g.\n',K2);

  nlevels = input(' No. of grid levels = ');
  nvec = zeros(nlevels,1);
  wfErr_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.
    uf_h = zeros(nnodes,1);
    
    %  Evaluate forcing function f(x) and diffusion coefficient kappa
    %  at subinterval midpoints. Set up left and right boundary conditions.

    k_mid = diffusion_coeff(x_mid);  
    wf_left = 0;
    wf_right = 0;
    
    %  Exactly compute FEM load vector b. 
    xm = x(2:end-1);
    xl = x(1:end-2);
    xr = x(3:end);
    
    b = 1/h*( gamma(alpha+1)/gamma(alpha+beta+1)*(2*xm.^(alpha+beta)...
        -xl.^(alpha+beta)-xr.^(alpha+beta)) - gamma(m+1)/gamma(m+beta+1)...
        *(2*xm.^(m+beta)-xl.^(m+beta)-xr.^(m+beta)) );
    b(1:(n+1)/2-1) = K1*b(1:(n+1)/2-1);
    b((n+1)/2+1:end) = K2*b((n+1)/2+1:end);
    b((n+1)/2) = 1/h*( gamma(alpha+1)/gamma(alpha+beta+1)*((K1+K2)*0.5^(alpha+beta)...
        -K1*(0.5-h)^(alpha+beta)-K2*(0.5+h)^(alpha+beta)) - gamma(m+1)/gamma(m+beta+1)...
        *((K1+K2)*0.5^(m+beta)-K1*(0.5-h)^(m+beta)-K2*(0.5+h)^(m+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 = (k_mid(1:n) + k_mid(2:n+1)) / h;
    Asub = -k_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. 
  
    wf_interior = A \ b;
    wf_h = [wf_left; wf_interior; wf_right];
    
    %  Compute error at mesh points, and record step size.
    
    wf_true = true_wf(x);
    wfError = wf_true - wf_h;    

    wfErr_vec(k) = compute_l2_error ( n, x, wf_h, true_wf );

    hvec(k) = h;
    u_true = true_u(x);  
    u_trueFine = true_u(xfine);

    %  Evaluate DxBeta_wfh.
    
    DxBeta_wfh = zeros(nelement+1,1);
    DxBeta_wfh(2) = wf_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) ).*wf_h(2:i);
        temp2 = ( - (i+1-(2:i)'+1).^(1-beta) + (i+1-(2:i)').^(1-beta) ).*wf_h(1:i-1);
        temp = sum(temp1) + sum(temp2);
        DxBeta_wfh(i+1) = 1/(gamma(2-beta)*h^beta)*(temp + (wf_h(i+1)-wf_h(i)));
    end
    
    %ErrorDxBeta_wf = DxBeta_wfh - DxBeta_wf;
    uError = zeros(nelement+1,1);
    DxBeta_wb = true_DxBeta_wb(x);
    coeff_end = - DxBeta_wfh(end)/DxBeta_wb(end);    
    uf_h(2:end-1) = DxBeta_wfh(2:end-1) + coeff_end*DxBeta_wb(2:end-1);
    ub = true_ub(x);
    u_h = uf_h + ub;
    error = u_h - u_true;

    uErrNorm = compute_frac_l2error_nonHom ( m, alpha, beta, K1, K2, coeff_end, nnodes, x, wf_h , 2);  
    fprintf('Error Norm of u : %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);
%**************************************************************************


