clear all;clc
% ME452w lab2
% Perform a mesh refinement study for a 1D steady state conduction problem. 
% A constant flux is applied on the left boundary and a fixed temperature on the 
% right boundary.

global X1 elcon;
for n=1:2
    nnodes=2^n+1
%*********** Geometry ***********************************
XL = [1];
XR = [nnodes];
%*********** Physical Parameters ***********************
k = 1;
X1 = linspace(0,1,nnodes);
Tb = 100;
fn = 100;
s = ones(nnodes,1)*100;
%*********** Iteration Parameters ***********************
[Xnodes Ynodes] = size(X1);
q0 = zeros(Xnodes.*Ynodes,1);

% Applies the BC on the intial condition
for ii=1:1:length(XR)
q0(XR(ii)) = +Tb;
end
q = q0;

%********************* Equations *******************************
%EQUATION_SET  0:
b0 = + asres1D(k, [], [], -1, 'A211L', q)...
- asres1D(1, [], [], 1, 'A200L', s);

% A constant flux(Neumann BC) is applied on the left end
b0(1)= b0(1)-fn;
% Dirichlet BC is enforced on the Right end
for ii=1:1:length(XR)
b0(XR(ii)) = 0;
end

%***************** Jacobians ***************************
%Equation 1 by null
jac0 = + asjac1D(k, [], [], -1, 'A211L', []);
jac0(1) = jac0(1) ;

DIFF = jac0;

for ii=1:1:length(XR)
jac0(XR(ii),:) = 0;
jac0(:,XR(ii)) = 0;
jac0(XR(ii),XR(ii)) = 1;
end

%***************** Solver ********************************************
dq = jac0\(-b0);

%************* Solve Iterate *****************************************
q = q + dq
% Note :- The arrays of q and dq will be identical since the Newton process 
% converges in one iteration (linear problem), as q0 = 0, and q=q+dq=0+dq=dq.
format long
enorm = 0.5*q'*DIFF*q
figure(2)
hold on
plot(X1,q,'md-')
title('Steady 1-D heat conduction', 'FontSize', 16)
xlabel('X1 (m)','FontSize', 16), ylabel('T (K)','FontSize', 16)
end