%-----------------------------------------------------------------
% Script 'Discretization Methods', Example 4.3
% (c) 12/2004 by A. Schmidt, IAM, Uni Stuttgart
%-----------------------------------------------------------------

%-----------------------------------------------------------------
% remark: the same problem given in Example 2.2 / 3.2 
%         this time solved using the FEM
%         difference: in Example 3.2, the heat flux q (heat flux
%         per unit cross section area in J/(s*m^2)) is given,
%         here the heat flux Q is an absolute value given in J/s.
%-----------------------------------------------------------------
clear all

% Variables
l = 1.0;                            % length of the rod
A = 1.0;                            % cross section area of the rod
T = 1.0;                            % temperature of the basin / right side
k = 1.0;                            % thermal conductivity of the rod
Q = 1.0;                            % heat flux into the rod
n = 6;                              % number of nodes (without ghost points)

%-----------------------------------------------------------------
% Calculation of the System (conductivity) matrix K_tt [Eq. (4.97)]
% first row (i=1)
K_tt(1,1) = 1.0;
K_tt(1,2) = -1.0;

% middle block
for i=2:1:n-2
  K_tt(i,i-1) = -1.0;
  K_tt(i,i) = 2.0;
  K_tt(i,i+1) = -1.0;
end

% last row (i=n-1)
K_tt(n-1,n-2) = -1.0;
K_tt(n-1,n-1) = 2.0;

K_tt = K_tt*(n-1)*k*A/l;

% Calculation of the right-hand side f [Eq. (4.97)]
f(1) = Q;
f(n-1) = T*(n-1)*k*A/l;
f = f';                             % make f a column vector

%-----------------------------------------------------------------
% Calculation of the vector of unknowns t (temperature): K_tt*t=f
t = K_tt\f;

%-----------------------------------------------------------------
% Calculation of the unknown heat flux Q_n at node n [Eq. (4.98)]
K_ft(n-1) = -1.0;
K_ft(n) = 1.0;
K_ft = K_ft*(n-1)*k*A/l;
t(n) = T;                           % given temperature at right side 
Q_n = K_ft*t;

%-----------------------------------------------------------------
% Output
% sample points
for i=1:1:n
  x(i) = (i-1)*l/(n-1);
end

% plot
A = ['heat flux at node n: Q_n = ' num2str(Q_n)];
figure(1)
plot(x,t,'k')
xlabel('Location / m')
ylabel('Temperature / °C')
title(A)