% FEM code for 1D 2nd order steady state equation
% Marcel Frehner, ETH Zurich, 2017

% GENERAL STUFF
    clear                                           % clear the current Workspace
    close all                                       % close all figure windows
    clc                                             % clear the Command Window
    
% GEOMETRICAL PARAMETERS
    Lx	=	10;                                     % length of model [m]
    
% PHYSICAL PARAMETERS
    kappa	=	1;                                  % thermal diffusivity [m2/s]
    source	=	1;                                  % heat source [K/s]
    
% NUMERICAL PARAMETERS
    el_tot      =	9;                              % #elements total
    n_tot       =	i_dont_know;                   	% #nodes total
    n_per_el	=	2;                              % #nodes per element
    
% CREATE NUMERICAL GRID
    dx      =   i_dont_know;                     	% distance between two nodes
    GCOORD  =	i_dont_know;                    	% array of coordinates
	
% LOCAL-TO-GLOBAL MAPPING
    EL_N    =   i_dont_know;                        % relates local to global node numbers per element
	
% BOUNDARY CONDITIONS
    bc_dof  =   [   1  n_tot];                      % dof's to which Dirichlet bc's are assigned
    bc_val  =   [   0  0    ];                      % value for these dof's
    
% INITIALIZATION OF ALL KINDS OF STUFF
    KG	=	zeros(n_tot,n_tot);                     % global stiffness matrix
    FG	=	zeros(n_tot,1);                         % global force vector
    
for iel = 1:el_tot % ELEMENT LOOP
    n_now   =   i_dont_know;                        % which nodes are in the current element?
    
    Kloc    =   i_dont_know;                        % local stiffness matrix
    Floc    =   i_dont_know;                        % local force vector
    
    % add local matrices to global ones
        KG(n_now,n_now)     =   i_dont_know;
        FG(n_now)           =   i_dont_know;
end % ELEMENT LOOP

% APPLY BOUNDARY CONDITIONS
    KG(bc_dof,:)        =   0;                      % set bc-rows to zero
    KG(bc_dof,bc_dof)   =   i_dont_know;            % put 1 onto the main diagonal
    FG(bc_dof)          =   i_dont_know;          	% set bc-value to the rhs
    
% SOLVER
    T	= KG\FG;
    
% ANALYTICAL SOLUTION
    x_ana 	=   0:Lx/1000:Lx;
    T_ana  	=   -1/2*source/kappa.*x_ana.^2 + ...
                (1/2*source*Lx/kappa+(bc_val(2)-bc_val(1))/Lx)*x_ana + bc_val(1);
    
% PLOTTING
    plot(x_ana,T_ana,'k-','LineWidth',2)            % plot analytical solution as thick black line
    hold on                                         % set the figure status to hold, to plot more stuff on top
    plot(GCOORD,T,'ro')                             % plot numerical solution as red circles
    
    set(gca,'FontSize',14)                          % set the font size of the figure to 14
    xlim([0 Lx])                                    % left and right limits of the x-axis
    ylim([0 14])                                    % bottom and top limits of the y-axis
    xlabel('x-coordinate [m]','FontSize',14)        % label of the x-axis
    ylabel('T [°]','FontSize',14)                   % label of the y-axis
    legend('Analytical','FE-solution')              % add a legend to the figure
    drawnow