% FEm.PSE template, steady conduction with variable conductivity


global X1;        % array for node coordinates

% constant data
h = 20;           % convection coefficient
Tr = 1500;        % reference temperature
Tb = 306.85282;   % fixed temperature

% uniform discretization, M=4
nnodes = 5;
X1 = linspace(0,1,nnodes);

% distributed data at the nodes
k = X1' * 10 + 10;

% load the FE matrix library
load femlib;

% assemble Matrix for WSe term [DIFF]
DIFF = asjac1D(1,[],k,-1,'A3011L',[]);

% complete Matrix for [HBC] on element 1
DIFFHBC = DIFF;   % copy DIFF
DIFFHBC(1,1) = DIFFHBC(1,1) + h; % modify the 1,1 entry

% set up data matrix (b)
b = zeros(nnodes,1);  % an nnodes by 1 matrix of zeros

% modify b for convective BC at node 1
b(1,1) = h * Tr;

% modify b for the Dirichlet boundry condition
b(nnodes,1) = Tb;

% modify DIFFHBC for Dirichlet BC direct solve
DHdiri = DIFFHBC; % copy DIFFHBC
DHdiri(nnodes,:) = zeros(1,nnodes); % zero out the last row
DHdiri(nnodes,nnodes) = 1;   % i.e., 1*QM = Tb

% solve the linear system
Q =  DHdiri \ b

% compute Energy Norm using [DIFF + HBC] (without Dirichlet BC)
format long; % display Enorm in long format
Enorm = 0.5*Q'*DIFFHBC*Q
format short;
