function KMaterial_initialGlobal = assembleInitialGlobalMaterialStiffnessMatrix(nDofTot, numEquations, element_data)
    KMaterial_initialGlobal = 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_elemInitial_inLocalRef = computeInitialLocalMaterialElemStiffnessMatrix(elem);

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

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