%  ************** Direct Finite Element Method ( DFEM ) **************
%  for Fractional Diffusion Equation With Smoothly Varying 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(x) = x.^2.*(1-x).^2
%  K(x) = 1./(1+x)
%
%--------------------------------------------------------------------------
%                     Shengfeng Zhu, 2013, ECNU
%--------------------------------------------------------------------------
  clear;
%  
  disp(' ************* Direct Finite Element Method ( DFEM ) ************* ') 
  disp(' for Fractional Diffusion Equation With Smoothly Varying Coefficient')
%
  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)
%
  diffusion_coeff = @(x) 1./(1+x); 
%  
  right_hand_side = @(x) -1./(gamma(beta+4)*(1+x).^2).*( (24*beta+48)*x.^(3+beta)...
    +(3+beta)*(12-12*beta)*x.^(2+beta)+(3+beta)*(2+beta)*(2*beta-12)*x.^(1+beta)...
    + (6+2*beta)*(2+beta)*(1+beta)*x.^(beta) );
%  
  true_u = @(x) x.^2.*(1-x).^2; 
%
  true_ux = @(x) 2*x.*(1-x).^2 + 2*x.^2.*(x-1);
%
  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;                    %  No. of elements.
      h = 1/nelement;                    %  Mesh spacing.
      hvec(k) = h;                       %  Mesh size saving
      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
      %
      CoeffMid = diffusion_coeff(x_mid);
      parameter = 1/(h^(1-beta)*GammaBeta2);
      % Diagonal: j == i
      vector = CoeffMid(1:end-1) + (3-2^beta1)*CoeffMid(2:end);
      A(2:end-1,2:end-1) = diag(parameter*vector);
      % j <= i-2
      for i = 4: n-1
          for j = 2: i-2
              A(i,j) = parameter*((CoeffMid(i-1)+3*CoeffMid(i))*(i-j+1)^beta1...
                  -3*(CoeffMid(i-1)+CoeffMid(i))*(i-j)^beta1+(3*CoeffMid(i-1)+CoeffMid(i))*(i-j-1)^beta1...
                  -CoeffMid(i-1)*(i-j-2)^beta1-CoeffMid(i)*(i-j+2)^beta1);
          end
      end
      % First Column of A
      A(2,1) = parameter*(CoeffMid(1)+CoeffMid(2)*(2^beta1-2));
      index = 3:n-1;
      A(index,1) = parameter*(CoeffMid(i-1)*(2*(index-1-1).^beta1-(index-1).^beta1-(index-1-2).^beta1)+...
          CoeffMid(i)*((index-1-1).^beta1-2*(index-1).^beta1+(index-1+1).^beta1));
      % j == i-1
      for i = 3:n-1
          A(i,i-1) = parameter*(CoeffMid(i-1)*(2^beta1-3)+CoeffMid(i)*(3*2^beta1-3^beta1-3));
      end
      % j == i+1
      for i = 2:n-1
          A(i,i+1) = -parameter*CoeffMid(i);
      end
      %
      % RHS: Force Vector Assembling
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      %  Quadrature definitions.
      quad_num = 10;
      %
      %  Quadrature points
      abscissa(1:10) = [ ...
          -0.973906528517172; ...
          -0.865063366688985; ...
          -0.679409568299024; ...
          -0.433395394129247; ...
          -0.148874338981631; ...
          0.148874338981631; ...
          0.433395394129247; ...
          0.679409568299024; ...
          0.865063366688985; ...
          0.973906528517172 ];
      %
      % Weights
      weight(1:10) = [...
          0.066671344308688; ...
          0.149451349150581; ...
          0.219086362515982; ...
          0.269266719309996; ...
          0.295524224714753; ...
          0.295524224714753; ...
          0.269266719309996; ...
          0.219086362515982; ...
          0.149451349150581; ...
          0.066671344308688 ];
      %
      w1_left = 0;
      w1_right = 0;
      %
      %  Make room for the matrix A and right hand side b.
      b = zeros ( n, 1 );
      %
      %  Equation 1 is the left boundary condition, U(0.0) = 0.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 variable for B(I).
          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;
              vm = ( xq - xl ) / ( xm - xl );
              fxq = right_hand_side ( xq );
              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;
              vm = ( xr - xq ) / ( xr - xm );
              fxq = right_hand_side ( xq );
              bm = bm + wq * ( fxq * vm );
          end
          b(i)     = bm;
      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.
      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 = 10;
      %
      %  Quadrature points
      %
      abscissa(1:10) = [ ...
          -0.973906528517172; ...
          -0.865063366688985; ...
          -0.679409568299024; ...
          -0.433395394129247; ...
          -0.148874338981631; ...
          0.148874338981631; ...
          0.433395394129247; ...
          0.679409568299024; ...
          0.865063366688985; ...
          0.973906528517172 ];
      %
      % Weights
      %
      weight(1:10) = [...
          0.066671344308688; ...
          0.149451349150581; ...
          0.219086362515982; ...
          0.269266719309996; ...
          0.295524224714753; ...
          0.295524224714753; ...
          0.269266719309996; ...
          0.219086362515982; ...
          0.149451349150581; ...
          0.066671344308688 ];
      %
      %  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 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: %d\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);
   
 