% FEm.PSE template script, axisymmetric 2-D steady heat transfer
% ref:  FE123, page 154.

% physical constants for this problem
k_pipe = 10;
k_insulation = 0.1;
h_flow = 20;
T_flow = 1500;
Tb = 306.8528;
h_air = 5;
T_air = 70;

% Generate node locations in polar coordinates.  The full mesh will be,
% in a sense, the outer product of a list of radial locations (r_seed)
% and a list of angular locations, in radians (theta_seed).
r_seed = linspace(1,2,5);
theta_seed = linspace(pi/2,0,5);

% generate the coordinate lists from the seeds
[r_grid,theta_grid] = meshgrid(r_seed,theta_seed);

% number of nodes in this problem
nnodes = prod(size(theta_grid));

% generate elements from the coordinate lists
elcon = delaunay(theta_grid,r_grid);

% transform the grid coordinates from polar to cartesian
[Xgrid,Ygrid] = pol2cart(theta_grid,r_grid);

% quick look at the mesh
figure(1)
trimesh(elcon,Xgrid,Ygrid,zeros(size(Xgrid)))
view(2)			        	% see a 2D view
axis square				% equal scaling of x and y axes

% calculate element zeta's and determinants.
[ZETAe,DETe] = do_ZETAe_2Dt(Xgrid,Ygrid,elcon);

% generate element diffusion matrices, [K]e
Ke = do_DIFFe_2Dt(ZETAe);

% assemble global diffusion matrix, [K].
K = asjac(k_pipe,[],[],-1,Ke,[],elcon,DETe,nnodes);

% implement the interior surface convective BC.  The nodes in the first
% column of Xgrid and Ygrid lie on the interior surface of the pipe.  Use
% that fact to build an ordered list of nodes for the interior surface.
ntheta = prod(size(theta_seed));	% number of theta values
flow_nodes = 1:ntheta;

% make a connectivity table for the line elements comprising the interior 
% surface
flow_linecon = do_linecon(flow_nodes);

% compute the lengths of the elements
le_flow = do_le2D(flow_linecon,Xgrid,Ygrid);

% assemble the interior surface convective BC contributions to [H] and {b}
load femlib;				% load 1D library matrices
H_flow = asjac(h_flow,[],[],1,A200L,[],flow_linecon,le_flow,nnodes);
TF = T_flow * ones(size(Xgrid));
b_flow = asres(h_flow,[],[],1,A200L,TF,flow_linecon,le_flow,nnodes);

% add convective BC contribution to [K]
KplusH = K + H_flow;

% Make a table to implement the Dirichlet BC.  The table has one row for
% each node.  The first column is a flag indicating whether the node is
% diri, the second column is the diri value.  The nodes in the last column
% of the grid arrays on the outer surface of the pipe.
diri_flag = zeros(size(Xgrid));
diri_flag(:,end) = ones(ntheta,1);
diri_value = diri_flag * Tb;
% stretch diri_flag and diri_value into columns, concatenate side by side
diri_table = [diri_flag(:), diri_value(:)];

% implement the diri BC
[KplusH_diri,b] = do_diri(KplusH,b_flow,diri_table);

% solve the linear system
Q = KplusH_diri \ b;

% calculate and display the energy semi-norm
format long
enorm = 1/2 * (Q' * KplusH * Q)
format short

% The linear solve returns the solution Q as a column matrix.  Reshape Q
% to a rectangular array with the same dimensions as the grid arrays.
% This form is more convienient for plotting the solution.
Q = reshape(Q,size(Xgrid));

% This plots Q as a function of radius for each value of theta in
% theta_seed.  These curves should all be identical.
figure(2)
plot(r_seed,Q');
xlabel('radius');
ylabel('temperature');
