%%
%% FEM code for 1D, 2nd order steady state diffusion equation with source
%% Marcel Frehner, ETH Zurich, 2012
%%

% GENERAL STUFF
clear all
close all
clc
    
% GEOMETRICAL PARAMETERS
Lx = 10;
    
% PHYSICAL PARAMETERS
kappa	 = 1;
source = -1;
    
% NUMERICAL PARAMETERS
el_tot   = 10;
n_tot    = el_tot+1;
n_per_el = 2;

% calculated from above
dx      = Lx/el_tot;
GCOORD  = 0:dx:Lx;
    
% LOCAL-TO-GLOBAL MAPPING
EL_N = [ 1:n_tot-1 ; 2:n_tot ];
    
% BOUNDARY CONDITIONS
bc_dof = [ 1 n_tot ];
bc_val = [ 0 0 ];
    
% INITIALIZE GLOBAL STIFFNESS MATRIX AND GLOBAL FORCE VECTOR
KG = zeros(n_tot,n_tot);
FG = zeros(n_tot,1);
    
% ELEMENT LOOP
for iel = 1:el_tot 

	% GET NODE INDICES FOR ELEMENT iel
	el_nodes   =   EL_N(:,iel);

	% DEFINE ELEMENT STIFFNESS MATRIX AND VECTOR
	Kloc    =   kappa/dx*[1 -1;-1 1];
	Floc    =   source*dx/2;
        
	% INSERT/SUM LOCAL ELEMENT STIFFNESS MATRICES INTO GLOBAL STIFFNESS MATRIX
	KG(el_nodes,el_nodes)  = KG(el_nodes,el_nodes) + Kloc;
	FG(el_nodes)           = FG(el_nodes)          + Floc;

end

% APPLY BOUNDARY CONDITIONS
for ibc = 1:length(bc_dof)

	% GET NODE INDICES FOR BOUNDARY CONDITION AT NODE ibc 
	bc_idx              =   bc_dof(ibc);

	% ZERO ROW AND INSERT 1 ON THE DIAGONAL
	KG(bc_idx,:)        =   0;
	KG(bc_idx,bc_idx)   =   1;

	% FORCE DIRICHLET VALUE INTO RHS VECTOR
	FG(bc_idx)          =   bc_val(ibc);

end
    
% SOLVE
T = KG\FG;

% COMPUTE ANALYTICAL SOLUTION
x_ana = 0:Lx/1000:Lx;
T_ana = -1/2*source/kappa.*x_ana.^2 + 1/2*source*Lx/kappa*x_ana;

% PLOT SOLUTIONS
plot(x_ana,T_ana,'k-',GCOORD,T,'ro');
xlim([0 Lx])
xlabel('x-coordinate [m]')
ylabel('T [°]')
legend('analytical solution','fem solution')
drawnow
