% FEm.PSE template - Lab01
% steady conduction with variable conductivity and source

clear all;
format short;

global X1;        % array for node coordinates

% load the FE matrix library
load femlib;

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

% Boundary Conditions
q  = 5;			% BC: Heat flux at x=0
Tb = 10;		% BC: Prescribed temperature at x=L=1

% average element data - conductivity
K = 0.1+0.001*X1';
Kbar = (K(1:nnodes-1)+K(2:nnodes))/2;

% nodal element data - source
S = 25*(1-X1.^2);

%-----------------------------------------------------------------------------

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

    % since [DIFF] is only term in LHS
    LHS = DIFF;
   
% Assemble all vectors on the RHS
	% assemble Vector for GWSe term SRC
	SRC = asres1D([],[],[],1,A200L,S);

	% form Vector for GWSe term BFLX
	BFLX = zeros(nnodes,1);  % an nnodes by 1 matrix of zeros
    % modify BFLX for flux BC at node 1
	BFLX(1,1) = q;

	% Build RHS
	RHS = SRC + BFLX;

% modify RHS for the Dirichlet boundary condition
RHSdiri = RHS;
RHSdiri(nnodes,1) = Tb;

% modify LHS for Dirichlet BC direct solve
LHSdiri = LHS; % copy DIFF
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;

Tleft = Q(1,1)

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

% create nice graphics
figure(1)
plot(X1,Q,'ko-')
xlabel('x')
ylabel('T')
title('Temperature Distribution - Lab 01')
text(0.1,30,'Zachariah Chambers')
text(0.1,20,'March 18, 20002')