% assembly of stiffness matrix
function [L,M] = StiffMassMat(nodes,triangles,lambda,mu);

    N = size(nodes,1);
    D = [-ones(2, 1), eye(2)];        % derivative operator of the base functions on the reference triangle

% index and value vectors for the stiffness and mass matrices
    indicesT = ones(3,1)*(1:3);
    indices = indicesT';
    I = zeros( 3 * 3 * size(triangles,1), 1 );
    J = I;
    VL11 = I;
    VL12 = I;
    VL21 = I;
    VL22 = I;
    VM = I;
    counter = 1;

% assemble system matrix
    for k = 1:size(triangles,1)         % loop through all elements {T}
        i = triangles(k, 1:3);          % get global indices
        x = nodes(i, 1:2);              % global nodes (as matrix rows)
        A = (D * x)';                   % Jacobian of coordinate transform from T_ref to T.
        AinvD = (D'/A);                 % gradient operator
        vol = abs(det(A)/2);            % triangle-volume

        % local anisotropy matrices
        B11 = [2*mu(k)+lambda(k),0;0,mu(k)];
        B12 = [0,lambda(k);mu(k),0];
        B21 = [0,mu(k);lambda(k),0];
        B22 = [mu(k),0;0,2*mu(k)+lambda(k)];

        % assembly of local stiffness matrices
        LT11 = vol*AinvD*B11*AinvD';
        LT12 = vol*AinvD*B12*AinvD';
        LT21 = vol*AinvD*B21*AinvD';
        LT22 = vol*AinvD*B22*AinvD';

        % assembly of local mass matrix
        MT = [2,1,1;1,2,1;1,1,2] * (vol/12);

        % update the global stiffness and mass matrices
        VL11(counter:counter+8) = LT11(:);
        VL12(counter:counter+8) = LT12(:);
        VL21(counter:counter+8) = LT21(:);
        VL22(counter:counter+8) = LT22(:);
        VM(counter:counter+8) = MT(:);
        I(counter:counter+8) = i(indices(:));
        J(counter:counter+8) = i(indicesT(:));
        counter = counter + 9;
    end

    L11 = sparse(I,J,VL11);
    L12 = sparse(I,J,VL12);
    L21 = sparse(I,J,VL21);
    L22 = sparse(I,J,VL22);

    L = [L11,L12;L21,L22];

    M11 = sparse(I,J,VM);

    M = [M11,sparse(N,N);sparse(N,N),M11];

    display('matrices assembled');
end
