%% %%%%%%%%%% Indirect Finite Element Method ( IFEM ) %%%%%%%%%%%%
% ***** Fractional Diffusion Equation: Piecewise constant coefficient *****
% Counter Example 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 < alpha <= 1
%  $$    u(0) = u_left, u(1) = u_right.
%
%  Use finite element method with piecewise linear basis functions to 
%  approximately solve 1-D steady-state diffusion equation
%           -D (K(x) Dwf) = f,  0 < x < 1,
%           wf(0) = u_left,  wf(1) = u_right.
%  The symmetric bilinear form a and linear functional "L" associated 
%  with this BVP are
%     a(u,v) = \int_0^1 kappa(x) du/dx dv/dx 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 )
%             wb(x) = 2*x/(1+kl)
%             DxBeta_wb_true = 2*x.^(1-beta)/(gamma(1-beta)*(1+Kl)*(1-beta))
%             f(x) = -2*Kl*x.^(beta-1)/GammaBeta
%            If x in (1/2,1]
%             u_true(x) = 2*(1-x)
%             K(x) = 1
%             wb(x) = (2*kl*x+1-kl)/(1+kl)
%             DxBeta_wb_true = (2*x.^(1-beta)+(2*Kl-2)*(x-0.5).^(-beta))
%                                       .../(gamma(1-beta)*(1+Kl)*(1-beta))
%             f(x) = -2*(x.^(beta-1)-2*(x-1/2).^(beta-1))/GammaBeta
%
%--------------------------------------------------------------------------
%                     Shengfeng Zhu, 2013, ECNU
%--------------------------------------------------------------------------
 clear;
 beta = input(' Parameter Beta ( in (0,1) ) = ');
 fprintf('Beta: %g\n', beta);
 GammaBeta = gamma(beta);
 
 %Kl = input(' Parameter Kl ( in (0,1] ) = ');
 Kl = 2^beta - 1.5;
 fprintf('Kl: %g\n', Kl);
 
 nlevels = input(' No. of grid levels = ');
 nvec = zeros(nlevels,1);
 Eh_vec = zeros(nlevels,1);
 Eh_wfvec = zeros(nlevels,1);
 hvec = zeros(nlevels,1);
 
 GammaBeta = gamma(beta);
 temp = 1/((beta+1)*beta*GammaBeta*(1+Kl)); %gamma(beta+2)*(1+Kl)
 
 a = @(Kl,x) Kl*(x<0.5) + 1*(x>0.5);
 f = @(Kl,beta,x) -2*Kl*x.^(beta-1)/GammaBeta.*(x<0.5)...
     + (-2*(x.^(beta-1)-2*(x-1/2).^(beta-1))/GammaBeta).*(x>0.5);

 Exact_u = @(x) 2*x.*(x <= 0.5) + (2*(1-x)).*(x>0.5);
 
 ExactDxBeta_wb = @(Kl,beta,x) (2*x.^(1-beta)/(gamma(1-beta)*(1+Kl)*(1-beta))).*(x<=0.5) +...
    ((2*x.^(1-beta)+(2*Kl - 2)*(x-0.5).^(1-beta))/(gamma(1-beta)*(1+Kl)*(1-beta))).*(x>0.5);

 Exact_wf = @(Kl,beta,x) ( temp*( 2*x.^(beta+1)*(1+Kl) +...
    (2^(1-beta)*((beta+3)-Kl*(beta+1))-4)*x ) ).* (x<=0.5) ...
     + ( temp*( (1+Kl)*(2*x.^(beta+1)-2^(1-beta)*(2*x-1).^(beta+1)+2^(1-beta)-2)...
     +((beta+1)-Kl*(beta+3)+2^(1+beta)*Kl)*2^(1-beta)*(1-x) ) ).*(x>0.5);
    
 Exact_wfx = @(Kl,beta,x) ( temp*( 2*x.^beta*(1+Kl)*(beta+1) + ...
     (2^(1-beta)*((beta+3)-Kl*(beta+1))-4) ) ).* (x<=0.5) ...
     + ( temp*( (1+Kl)*(2*(beta+1)*x.^beta-2^(2-beta)*(beta+1)*(2*x-1).^beta)...
     - 2^(1-beta)*((beta+1)-Kl*(beta+3)+2^(1+beta)*Kl) ) ).*(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.
     nvec(k) = n;
     h = 1/nelement;                    %  Mesh spacing.
     x = (0:h:1)';                      %  Mesh points.
     hvec(k) = h;
     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_true = Exact_u(x);
     DxBeta_wb_true = ExactDxBeta_wb(Kl,beta,x);
     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;   %  No. of grid points.
     %  Quadrature definitions.
     %
     quad_num = 2;
     abscissa(1) = -0.577350269189625764509148780502;
     abscissa(2) = +0.577350269189625764509148780502;
     weight(1) = 1.0;
     weight(2) = 1.0;
     %
     %  Make room for the matrix A and right hand side b.
     %
     A = zeros ( n, n );
     b = zeros ( n, 1 );
     b1 = 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, where:
             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 = a ( Kl,xq );
             fxq = f ( Kl,beta,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, where:
             %
             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 = a ( Kl,xq );
             fxq = f ( Kl,beta,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
     %%%%
     parameter = -2/(h*gamma(beta+2));
     for i = 2: (n-1)/2
         b1(i) = Kl*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
     b1((n+1)/2) = Kl*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
         b1(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
     %
     %  Equation N is the right boundary condition, U(1.0) = 0.0;
     %
     A(n,n) = 1.0;
     b(n) = 0.0;
     
     % Solve the linear system.
     wf_h = A \ b1;
     
     %  Compute max absolute error at mesh points, and record step size.
     wf_true = Exact_wf(Kl,beta,x);
     wf_error = wf_h-wf_true;
     indication = norm(wf_error,2)*h^.5;
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     %
     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);
         wf_l = wf_h(i);
         wf_r = wf_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      ) * wf_l ...
                 + (      xq - xl ) * wf_r ) ...
                 / ( xr      - xl );
             uxq = ( wf_r - wf_l ) / ( xr - xl );
             eq = Exact_wf ( Kl,beta,xq );
             exq = Exact_wfx ( Kl,beta,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 of wf: %d\n', l2_error);
     Eh_wfvec(k) = l2_error;
     %   fprintf('seminorm_error of wf: %d\n', seminorm_error);
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     %  Evaluate DxBeta_wf_h.
     DxBeta_wf_h = zeros(nelement+1,1);
     DxBeta_wf_h(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_wf_h(i+1) = 1/(gamma(2-beta)*h^beta)*(temp + (wf_h(i+1)-wf_h(i)));
     end
     %  Evaluate DxBeta_wf Numerically.
     DxBeta_wf = zeros(nelement+1,1);
     Dwfm = zeros(nelement,1); % Function values on the middle points of each element
     c1 = 2*(1+Kl)*(beta+1)/(gamma(beta+2)*(1+Kl));
     c2 = (2^(1-beta)*((beta+3)-Kl*(beta+1))-4)/(gamma(beta+2)*(1+Kl));
     d1 = c1;
     d2 = (-2^(2-beta)*(1+Kl)*(beta+1))/(gamma(beta+2)*(1+Kl));
     d3 = -2^(1-beta)*((beta+1)-Kl*(beta+3)+2^(1+beta)*Kl)/(gamma(beta+2)*(1+Kl));
     Dwfm(1:nelement/2) = c1*(x(2:nelement/2+1)-h/2).^(beta) + c2;
     Dwfm(nelement/2+1:end) = d1*(x(nelement/2+2:end)-h/2).^(beta) + ...
         d2*(2*(x(nelement/2+2:end)-h/2) - 1).^beta + d3;
     DxBeta_wf(2) = Dwfm(1)/gamma(2-beta)*h^(1-beta);
     for i = 2:nelement
         temp1 = (i+1-(2:i)'+1).^(1-beta) - (i+1-(2:i)').^(1-beta);
         temp = temp1.*Dwfm(2:i);
         DxBeta_wf(i+1) = h^(1-beta)/gamma(2-beta)*(sum(temp) + Dwfm(i));
     end
     ErrorDxBeta_wf =  DxBeta_wf_h - DxBeta_wf;
     uError = zeros(nelement+1,1);
     u_h(2:end-1) = DxBeta_wf_h(2:end-1) ...
         - DxBeta_wf_h(end)/DxBeta_wb_true(end)*DxBeta_wb_true(2:end-1);
     error = u_h - u_true;
     coeff_end = - DxBeta_wf_h(end)/DxBeta_wb_true(end);
     uError(2:end-1) = ErrorDxBeta_wf(2:end-1) + ...
      (- DxBeta_wf_h(end)+DxBeta_wf(end))/DxBeta_wb_true(end)*DxBeta_wb_true(2:end-1);
     uErrNorm = norm(uError,2)*h^(1/2);
     fprintf('Error Norm of u : %E\n',uErrNorm);
     Eh_vec(k) = uErrNorm;
 end
  %  Compute rate of convergence.
  conv_rate = diff(log(Eh_wfvec)) ./ diff(log(hvec));
  fprintf(' wf Convergence Order: %f\n',conv_rate);
  %
  conv_rate = diff(log(Eh_vec)) ./ diff(log(hvec));
  fprintf(' u Convergence Order: %f\n',conv_rate);
