% FEM solver for d2T/dx2 + q = 0 where q = 50 exp(x)
%
% BC's: T(-1) = 100 and T(1) = 100.
%
% Gaussian quadrature is used in evaluating the forcing integral.
%
% Note: the finite element degrees of freedom are
%       stored in the vector, v.  

% Number of elements
Ne = 10;
x = linspace(-1,1,Ne+1);

% Set quadrature rule 
Nq = 2;
if (Nq == 1),
  alphaq(1) = 2.0; xiq(1) = 0.0;
elseif (Nq == 2),
  alphaq(1) = 1.0; xiq(1) = -1/sqrt(3);
  alphaq(2) = 1.0; xiq(2) =  1/sqrt(3);
else
  fprintf('Error: Unknown quadrature rule (Nq = %i)\n',Nq);
  return;
end

% Zero stiffness matrix
K = zeros(Ne+1, Ne+1);
b = zeros(Ne+1, 1);

% Loop over all elements and calculate stiffness and residuals
for ii = 1:Ne,

  kn1 = ii;
  kn2 = ii+1;
  
  x1 = x(kn1);
  x2 = x(kn2);
  
  dx    = x2 - x1;
  
  % Add contribution to kn1 weighted residual due to kn1 function
  K(kn1, kn1) = K(kn1, kn1) - (1/dx);

  % Add contribution to kn1 weighted residual due to kn2 function
  K(kn1, kn2) = K(kn1, kn2) + (1/dx);
  
  % Add contribution to kn2 weighted residual due to kn1 function
  K(kn2, kn1) = K(kn2, kn1) + (1/dx);

  % Add contribution to kn2 weighted residual due to kn2 function
  K(kn2, kn2) = K(kn2, kn2) - (1/dx);

  % Evaluate forcing term using quadrature
  for nn = 1:Nq,

    % Get xi location of quadrature point
    xi = xiq(nn);
    
    % Calculate x location of quadrature point
    xq = x1 + 0.5*(1+xi)*dx;
    
    % Calculate q
    qq = 50*exp(xq);
    
    % Calculate phi1 and phi2
    phi1 = 0.5*(1-xi);
    phi2 = 0.5*(1+xi);
    
    % Subtract forcing term to kn1 weighted residual 
    b(kn1) = b(kn1) - alphaq(nn)*0.5*phi1*qq*dx;
  
    % Add forcing term to kn2 weighted residual
    b(kn2) = b(kn2) - alphaq(nn)*0.5*phi2*qq*dx;
    
  end
  
end


% Set Dirichlet conditions at x=0
kn1 = 1;
K(kn1,:)    = zeros(size(1,Ne+1));
K(kn1, kn1) = 1.0;
b(kn1)      = 100.0;


% Set Dirichlet conditions at x=1
kn1 = Ne+1;
K(kn1,:)    = zeros(size(1,Ne+1));
K(kn1, kn1) = 1.0;
b(kn1)      = 100.0;

  
% Solve for solution
v = K\b;


% Plot solution
figure(1);
plot(x,v,'*-');
xlabel('x');
ylabel('T');

% For the exact solution, we need to use finer spacing to plot 
% it correctly.  If we only plot it at the nodes of the FEM mesh,
% the exact solution would also look linear between the nodes.  To
% make sure there is always enough resolution relative to the FEM
% nodes, the size of the vector for plotting the exact solution is
% set to be 20 times the number of FEM nodes.
Npt = 20*Ne+1;
xe = linspace(-1,1,Npt);
Te = -50*exp(xe) + 50*xe*sinh(1) + 100 + 50*cosh(1);
hold on; plot(xe,Te); hold off;

% Plot the error.  To do this, calculate the error on the same
% set of points in which the exact solution was plot.  This
% requires that the location of the point xx(i) be found in the
% FEM mesh to construct the true solution at this point by linearly
% interpolating between the two nodes of the FEM mesh.

verr(1) = v(1) - Te(1);
h = x(2)-x(1);
for i = 2:Npt-1,
  xxi = xe(i);
  Ti  = Te(i);
  j = floor((xxi-xe(1))/h) + 1;
  x0 = x(j);
  x1 = x(j+1);
  v0 = v(j);
  v1 = v(j+1);
  xi = 2*(xxi - x0)/(x1-x0)-1;  % This gives xi between +/-1
  vi = 0.5*(1-xi)*v0 + 0.5*(1+xi)*v1;
  verr(i) = vi - Ti;
end  
verr(Npt) = v(Ne+1) - Te(Npt);
  
figure(2);
plot(xe,verr);
xlabel('x');
ylabel('Error');