function K = computeInitialLocalMaterialElemStiffnessMatrix(element_data)
    if strcmp(element_data{1}, '2dTruss')
        E = element_data{5};
        A = element_data{6};
        L = norm(element_data{4} - element_data{3});
        K = computeInitialLocalElemStiffnessMatrix_2dTruss(E, A, L);
        
    elseif strcmp(element_data{1}, '2dElasticBeam')
        E = element_data{5};
        A = element_data{6};
        I = element_data{7};
        L = norm(element_data{4} - element_data{3});
        K = computeInitialLocalElemStiffnessMatrix_2dBeam(E, A, I, L);
        
    elseif strcmp(element_data{1}, 'spring')
        K = element_data{4};  
    end
end


function K = computeInitialLocalElemStiffnessMatrix_2dTruss(E, A, L)
    % Computes the local element stiffness matrix for a 2D truss element
    
    K = (E * A / L) * [1, 0, -1, 0;
                       0, 0,  0, 0;
                      -1, 0,  1, 0;
                       0, 0,  0, 0];
end


function K = computeInitialLocalElemStiffnessMatrix_2dBeam(E, A, I, L)
    % Computes the local element stiffness matrix for a 2D beam element
    
    K = [ E * A / L,         0,          0,       -E * A / L,         0,           0;
             0,   12 * E * I / L^3,   6 * E * I / L^2,         0, -12 * E * I / L^3,   6 * E * I / L^2;
             0,    6 * E * I / L^2,    4 * E * I / L,         0,  -6 * E * I / L^2,    2 * E * I / L;
         -E * A / L,         0,          0,         E * A / L,         0,           0;
             0,  -12 * E * I / L^3,  -6 * E * I / L^2,         0,  12 * E * I / L^3,  -6 * E * I / L^2;
             0,    6 * E * I / L^2,    2 * E * I / L,         0,  -6 * E * I / L^2,    4 * E * I / L];
end