function kg = geometricStiffMtxUL(Anl,Elem,U)
    % Element properties
    A  = Elem.A;
    I  = Elem.I;
    L  = elemLength(Elem,U);

    % Internal forces
    P  = Elem.fi(4);
    M1 = Elem.fi(3);
    M2 = Elem.fi(6);

    % Simplifications
    L2  = L*L;
    L3  = L2*L;
    PI  = P*I;
    AL  = A*L;
    AL2 = A*L2;
    AL3 = A*L3;

    % First order terms
    kg1 = P/L *[ 1        0        0       -1        0        0;
                 0        6/5      L/10     0       -6/5      L/10;
                 0        L/10     2*L2/15  0       -L/10    -L2/30;
                -1        0        0        1        0        0;
                 0       -6/5     -L/10     0        6/5     -L/10;
                 0        L/10    -L2/30    0       -L/10    2*L2/15 ];
    % Second order terms
%     if (Anl.geom_mtx == 1)
        kg2 = zeros(6,6);

%     elseif (Anl.geom_mtx == 2)
%         kg2 = [ 0      0          -M1/L        0      0           -M2/L;
%                 0      12*PI/AL3   6*PI/AL2    0     -12*PI/AL3    6*PI/AL2;
%                -M1/L   6*PI/AL2    4*PI/AL     M1/L  -6*PI/AL2     2*PI/AL;
%                 0      0           M1/L        0      0            M2/L;
%                 0     -12*PI/AL3  -6*PI/AL2    0      12*PI/AL3   -6*PI/AL2;
%                -M2/L   6*PI/AL2    2*PI/AL     M2/L  -6*PI/AL2     4*PI/AL ];
%     end

    % Geometric matrix
    kg = kg1 + kg2;
end