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

% Set number of Gauss points (only used for forcing term in this example)
NGf = 1;
if (NGf == 3),
  xiGf = [-sqrt(3/5);          0; sqrt(3/5)];
  aGf  = [       5/9;        8/9;       5/9];
elseif (NGf == 2),
  xiGf = [-1/sqrt(3); +1/sqrt(3)];
  aGf  = [       1.0;        1.0];
else, 
  NGf = 1;
  xiGf = [0.0];
  aGf  = [2.0];
end

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

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

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

  kn1 = 1 + 2*(ii-1);
  kn2 = 2 + 2*(ii-1);
  kn3 = 3 + 2*(ii-1);
  
  x1 = x(ii);
  x3 = x(ii+1);
  
  dx    = x3 - x1;
  dxidx = 2/dx;
  dxdxi = 1/dxidx;
  
  % Add contribution to kn1 weighted residual 
  K(kn1, kn1) = K(kn1, kn1) - dxidx*( 7/6);
  K(kn1, kn2) = K(kn1, kn2) - dxidx*(-4/3);
  K(kn1, kn3) = K(kn1, kn3) - dxidx*( 1/6);

  % Add contribution to kn2 weighted residual 
  K(kn2, kn1) = K(kn2, kn1) - dxidx*(-4/3);
  K(kn2, kn2) = K(kn2, kn2) - dxidx*( 8/3);
  K(kn2, kn3) = K(kn2, kn3) - dxidx*(-4/3);

  % Add contribution to kn3 weighted residual 
  K(kn3, kn1) = K(kn3, kn1) - dxidx*( 1/6);
  K(kn3, kn2) = K(kn3, kn2) - dxidx*(-4/3);
  K(kn3, kn3) = K(kn3, kn3) - dxidx*( 7/6);
  
  % Use Gaussian quadrature to evaluate forcing term integral  
  for nn = 1:NGf,

    % Get xi for Gauss point    
    xiG = xiGf(nn); 
    
    % Find N1, N2 and N3 (i.e. weighting/intepolants) at xiG
    N1 = -0.5*xiG*(1-xiG);
    N2 =  (1-xiG)*(1+xiG);
    N3 =  0.5*xiG*(1+xiG);
    
    % Find x for Gauss point
    xG = 0.5*(1-xiG)*x1 + 0.5*(1+xiG)*x3;
    
    % Find f for Gauss point
    fG = -50*exp(xG);
    
    % Evaluate integrand at Gauss point for weight functions at nodes
    gG1 = N1*fG*dxdxi;
    gG2 = N2*fG*dxdxi;
    gG3 = N3*fG*dxdxi;
    
    % Send to correct right-hand side term
    b(kn1) = b(kn1) + aGf(nn)*gG1;
    b(kn2) = b(kn2) + aGf(nn)*gG2;
    b(kn3) = b(kn3) + aGf(nn)*gG3;
  
  end
  
end


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


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

  
% Solve for solution
v = K\b;


% Plot it and compare.  Note: since even the finite element
% solution varies more than linearly across an element, we need to
% subdivide each element, evaluate the basis functions, and plot
% the FEM solution to see the higher order variations.

Nplot = 20; % Number of points per element to plot
nnn = 0;
for ii = 1:Ne,

  kn1 = 1 + 2*(ii-1);
  kn2 = 2 + 2*(ii-1);
  kn3 = 3 + 2*(ii-1);
  
  x1 = x(ii);
  x3 = x(ii+1);

  v1 = v(kn1);
  v2 = v(kn2);
  v3 = v(kn3);

  for nn = 1:Nplot,
    
    % Get xi for plot point    
    xiG = -1 + 2*(nn-1)/(Nplot-1);
    
    % Find N1, N2 and N3 (i.e. weighting/intepolants) at xiG
    N1 = -0.5*xiG*(1-xiG);
    N2 =  (1-xiG)*(1+xiG);
    N3 =  0.5*xiG*(1+xiG);
    
    % Find x and v for plot point
    xG = 0.5*(1-xiG)*x1 + 0.5*(1+xiG)*x3;
    vG = v1*N1 + v2*N2 + v3*N3;
    
    nnn = nnn + 1;
    xp(nnn) = xG;
    vp(nnn) = vG;
    up(nnn) = -50*exp(xG) + 50*xG*sinh(1) + 100 + 50*cosh(1);
    
  end
end

figure(1);
plot(xp,vp,'r');hold on;
plot(xp,up); hold off;
xlabel('x');
ylabel('u');

figure(2);
plot(xp,vp-up);
xlabel('x');
ylabel('Error');