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')
        ke_s = element_data{4}(6);
        K = [ke_s, -ke_s;
            -ke_s, ke_s];
        
    elseif strcmp(element_data{1}, '2dModElasticBeam')
        E=element_data{5};
        A=element_data{6};
        Ie=element_data{11};
        L = norm(element_data{4} - element_data{3});
        S22mod=element_data{7};
        S23mod=element_data{8};
        S32mod=element_data{9};
        S33mod=element_data{10};
        K=computeInitialLocalMaterialElemStiffnessMatrix_2dModElasticBeam(E,A,S22mod,S23mod,S32mod,S33mod,Ie,L);
        
        
    elseif strcmp(element_data{1}, '2dDispBeamColumn')
        E=element_data{7}(1);
        L = norm(element_data{4} - element_data{3});
        coordinates_fibers=element_data{8};
        area_fibers=element_data{9};
        A=0;
        I=0;
        for iFiber=1:size(coordinates_fibers,1)
            A=A+area_fibers(iFiber);
            I=I+area_fibers(iFiber)*(coordinates_fibers(iFiber,2))^2;
        end
        K=computeInitialLocalElemStiffnessMatrix_2dBeam(E,A,I,L);
        
    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


function Kloc=computeInitialLocalMaterialElemStiffnessMatrix_2dModElasticBeam(E,A,S22mod,S23mod,S32mod,S33mod,Ie,l)
% Computes the local element stiffness matrix for a 2D elastic beam
% element with stiffness modifiers

% Tangent stiffness matrix in the basic reference frame
KBar = [
    E * A / l, 0, 0;
    0, S22mod * E * Ie / l, S23mod * E * Ie / l;
    0, S32mod * E * Ie / l, S33mod * E * Ie / l
    ];

% Linear transformation matrix
LTransfo = [
    -1, 0, 0, 1, 0, 0;
    0, 1/l, 1, 0, -1/l, 0;
    0, 1/l, 0, 0, -1/l, 1
    ];

% Tangent stiffness matrix in the local reference frame
Kloc = LTransfo' * KBar * LTransfo;

end