%Non-linear Heat Conduction
%Lab04 ME 422 - Brett Werner 7/12/02

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=32
nnodes=33;
X1 = linspace(XL,XR,nnodes);

% Boundary Conditions
Tl = 400;
Tr = 1000;

% coeficcients of conductivity
a= 73.2;
b= -0.04;
c= -0.000003;

% nodal element data - source
S = 10*linspace(1,1,nnodes);


%-----------------------------------------------------------------------------
% SOLVE FOR TEMPERATURE
% outer iteration loop
Q = pi*ones(nnodes,1); % initialize Q
Q(1,1) = Tl; % impose Dirichlet data
Q(nnodes,1) = Tr;
dQ = pi*ones(nnodes,1); % initialize dQ
iter_num = 1; % iteration counter
Q_big(:,iter_num) = Q; % big matrix of nodal Q

newt_crit=1e-6;
twoc=2*c;

while max(dQ)>=newt_crit

Qsq=Q.^2;    
    
% Assemble all contributions to [JAC]
jac = asjac1d(a,[],[],-1,A211L, []);
%jac = jac + asjac1d(b,[],Q ,-1,A3011L,[]);
jac = jac + asjac1d(b,[],Q ,-1,A3110L,[]);
jac = jac + asjac1d(c,[],Qsq ,-1,A3011L,[]);
jac = jac + asjac1d(twoc,[],Q ,-1,A3110L,[Q]);

% Assemble all contributions to {FQ}
res = asres1d(a, [],[],-1,A211L, Q);
res = res + asres1d(b, [],Q ,-1,A3011L,Q);
res = res + asres1d(c, [],Qsq ,-1,A3011L,Q);
res = res + asres1d(-1,[],[],1, A200L, S);

% Modify [JAC] for Dirichlet data
jac(1,:) = zeros(1,nnodes);
jac(1,1) = 1;
jac(nnodes,:) = zeros(1,nnodes);
jac(nnodes,nnodes) = 1;

% Modify {RES} for Dirichlet data
res(1,1) = 0;
res(nnodes,1) = 0;

% Solve for dQ
dQ = jac \ (-res);

% Update Q
Q = Q + dQ;
iter_num = iter_num+1; % increment iteration counter
Q_big(:,iter_num) = Q; % append current value of Q
dQ_max(iter_num-1,1) = max(dQ); % store max(dQ)
max(dQ) % echo max(dQ)
end


Q

% calculate log10 of dQ for convergence verification
log_dQ_x = log10(dQ_max(1:end-1,1))
log_dQ_y = log10(dQ_max(2:end,1))
conv = polyfit(log_dQ_x,log_dQ_y,1)