%Brett Werner ME 422 Lab 05 Matlab Code 
%- Transient Heat Conduction

clear all;
clf reset;
format long e;		% double precision exponent output

global X1;        % array for node coordinates
load femlib;		% load FEA libarary

% constant data
alpha = 17.7E-06;	% alpha for carbon steel (m^2/s)

% spatial data
XL = 0;				% left boundary
XR = 1;				% right boundary
nnodes = 33;		% 32 elements
X1 = linspace(XL,XR,nnodes);

% temporal data
theta = 0;					% implicitness factor
max_its = 5;					% max number of Newton iterations
newt_crit = 0.000001;			% Newton convergence criteria
to = 0;							% initial time (s)
tf = 2000;						% final time (s)
DELTA_T = 100;					% time step size
nsteps = (tf-to)/DELTA_T;	% number of time steps

% boundary conditions
TL = 273;			% left dirichlet temperature (K)
TR = 273;			% right dirichlet temperature (K)

%initial condition - sin distribution
T_ini = .1*sin(X1/(XR-XL)*pi)'+273;
T_ini(1) = TL;		% make sure BC's are exactly enforced
T_ini(end) = TR;

Q_n = T_ini;						% Initialize Qn with IC data
Q_np1 = T_ini;						% Initialize Qnp1 with guess of IC data
DELTA_Q = Q_np1 - Q_n;			% Form temporal change in Q
delta_Q = pi*ones(nnodes,1);	% Initialize delta Q
num_its = 1;

Q_big(:,1) = T_ini;
flops(0);
%**************************************************************************
%OUTER TIME STEPPING LOOP
for t_loop=1:nsteps
      
   while ((max(abs(delta_Q))>newt_crit)&(num_its <= max_its))

      % Form temporal MASS contribution to total residual FQp using DELTA_Q
       RES_MASS = asres1D([],[],[],1,A200L,DELTA_Q);
      
      % Form spatial residual contributions to total residual FQp
   	% {RES}_np1 term
   	   RES_np1 = asres1d(alpha,[],[],-1,A211L,Q_np1);
   	% {RES)_n term
       RES_n = asres1d(alpha,[],[],-1,A211L,Q_n);
      
      % Form total residual term FQp
       FQp = RES_MASS + DELTA_T*(theta*RES_np1 + (1-theta)*RES_n);
      
      %*******************************************************************
      % Form temporal MASS contribution to jacobian JACp
        JAC_MASS = asjac1d([],[],[],1,A200L,[]);
      
      % Form spatial residual contributions to jacobian JACp
        JAC_np1 = asjac1d(alpha,[],[],-1,A211L,[]);
      % Form total jacobian term JACp
        JACp = JAC_MASS + DELTA_T*theta*JAC_np1;
      
      %*******************************************************************
      % Modify JACp for dirichlet data
      JACp(1,:) = zeros(1,nnodes);
      JACp(1,1) = 1;
      JACp(nnodes,:) = zeros(1,nnodes);
      JACp(nnodes,nnodes) = 1;
      
      % Modify FQp for Dirichlet data
      FQp(1,1) = 0;
      FQp(nnodes,1) = 0;
      
      % Solve for incremental change in Q_np1
       delta_Q = JACp \ -FQp;
      
      % echo stuff
      status = [t_loop num_its max(abs(delta_Q))]
      
      % Update Q_np1 with incremental change
        Q_np1 = Q_np1 + delta_Q;
      
      % Recover new value of DELTA_Q
        DELTA_Q_np1 = Q_np1 - Q_n;
      
      % increment numits counter
      num_its = num_its + 1;
   end
   
   % Moving to new time station, the current Q now becomes the previous Q
   Q_big(:,t_loop+1) = Q_np1;
   Q_n = Q_np1;
   DELTA_Q = Q_np1 - Q_n;			% Form temporal change in Q
   delta_Q = pi*ones(nnodes,1);	% Initialize delta Q
   num_its = 1;

end
flops
Q_big(:,nsteps)
plot(X1,Q_big(:,1),X1,Q_big(:,end))