% Lab 7 file: Transient Heat Conduction
clear all; close all;

% Number of element nodes
nen = 3;

% Mesh density
numel    = 200;
numnodes = (nen-1)*numel + 1;

% Properties and loading
kc  = 54.0; % Conductivity        [W/mK]
L   = 1;    % Half-span of domain [m]
% [COMMENT: This is a good place to set the other material properties]

% Set up isoparametric shape functions and their derivatives
n1 = @(xi) xi.*(xi-1)/2; dn1 = @(xi) xi-1/2;
n2 = @(xi) (1-xi.^2);    dn2 = @(xi) -2*xi;
n3 = @(xi) xi.*(xi+1)/2; dn3 = @(xi) xi+1/2;

% Store these functions in two cell arrays for convenience
N = {n1, n2, n3};  dN = {dn1, dn2, dn3};

% Array of physical nodal coordinates
X = linspace(-L,L,numnodes);

% LM array
lm = [1:(nen-1):numnodes-nen+1;
      2:(nen-1):numnodes-nen+2;
      3:(nen-1):numnodes];

% Allocate stiffness and force vector
K = spalloc(numnodes,numnodes,5*numnodes);
F = zeros(numnodes,1);
% [COMMENT: This is a good place to allocate memory for the thermal mass]

% Quadrature point locations and corresponding weights
xg    = [-sqrt(3/5),   0,   sqrt(3/5)];
wt    = [       5/9, 8/9,         5/9];
nint  = length(wt); % Number of quadrature points

% Loop the elements
% [COMMENT: You will need to modify this loop to compute the thermal mass terms]
for e = 1:numel

    % Get element's nodal coordinates from global array
    ll = lm(:,e);
    xl = X( ll );
    
    % Set up isoparametric map
    xe = @(xi) xl(1)*n1(xi) + xl(2)*n2(xi) + xl(3)*n3(xi);

    % Set up Jacobian of the isoparametric map
    je = @(xi) xl(1)*dn1(xi) + xl(2)*dn2(xi) + xl(3)*dn3(xi);

    % Initialize the element force and stiffness matrix
    f = zeros(nen,1);
    k = zeros(nen,nen);

    % Loop over the quadrature points to accumulate the values for the
    % needed integrals
    for l = 1:nint 

        for i = 1:nen
                
                % Element force vector for this problem is zero

            for j = 1:nen

                % Compute the stiffness matrix
                k(i,j) = k(i,j) + wt(l)*dN{i}(xg(l))*kc*dN{j}(xg(l))/je(xg(l));

            end

        end

    end
    
    F(ll)    = F(ll)    + f;
    K(ll,ll) = K(ll,ll) + k;
   
end


% [COMMENT: Add code for the transient time stepping]



