% FEm.PSE template, steady conduction with variable conductivity and no source

clear all;
global X1;        % array for node coordinates

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

% mesh parameters
XL = 0;
XR = 1;
% uniform discretization, M=1
nnodes = 9;
X1 = linspace(XL,XR,nnodes);

% nodal element data - cond
COND = X1' * 10 + 10;

%-----------------------------------------------------------------------------
% load the FE matrix library
load femlib;

% Assemble all matrices on the LHS:
	% assemble Matrix for GWSe term [DIFF]
	DIFF = asjac1D([],[],COND,-1,A3011L,[]);

	% complete Matrix for GWSe term [HBC]
	HBC = zeros(nnodes,nnodes);
    HBC(1,1) = h;
   
   % Build LHS
   LHS = DIFF + HBC;
   
% Assemble all vectors on the RHS
	% Assemble Vector for GWSe term HTR
   HTR = zeros(nnodes,1);
   HTR(1,1) = h*Tr;
   
   % Build RHS
   RHS = HTR;
   
% Modify RHS for dirichlet data
RHSdiri = RHS;
RHSdiri(nnodes,1) = Tb;

% Modify LHS for dirichlet data
LHSdiri = LHS;									% Copy LHS
LHSdiri(nnodes,:)=zeros(1,nnodes);		% Zero out the last row
LHSdiri(nnodes,nnodes)=1;					% i.e., 1*QM = Tb

% solve the linear system
Q =  LHSdiri \ RHSdiri

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