%==========================================================================
function angle = elemAngle(Elem,U)
    % Element nodes
    N1 = Elem.n1;
    N2 = Elem.n2;

    % Nodal displacements
    dx1 = U(N1.dof(1));
    dy1 = U(N1.dof(2));
    dx2 = U(N2.dof(1));
    dy2 = U(N2.dof(2));

    % New nodal coordinates
    x1 = N1.x + dx1;
    y1 = N1.y + dy1;
    x2 = N2.x + dx2;
    y2 = N2.y + dy2;

    % Element angle with global axes
    angle = atan2(y2-y1,x2-x1);
end