function K_global = assembleGlobalStiffnessMatrix(nDofTot, numEquations, element_data)
    K_global = zeros(nDofTot, nDofTot);  % Initialize the global stiffness matrix

    for elemIdx = 1:length(element_data)
        elem = element_data{elemIdx};

        % Compute element stiffness matrices in the local reference system
        K_elem_inLocalRef = computeLocalElemStiffnessMatrix(elem);

        % Compute element stiffness matrices in the global reference system
        R = computeRotationMatrix(elem);
        K_elem_inGlobalRef = R' * K_elem_inLocalRef * R;  

        % Assemble the global stiffness matrix
        dof_elem = numEquations{elemIdx};  % Get the degrees of freedom for the element
        K_global(dof_elem, dof_elem) = K_global(dof_elem, dof_elem) + K_elem_inGlobalRef;  % Assemble into the global matrix
    end
end