function f = elastic_force_nonregular(F, gradN, nodeNum, edofMat, lx, u, l)

    f = sparse(nodeNum*2, 1);
    jac = lx^2/4;

    for ele = 1:size(edofMat, 1)
        fe = zeros(2, 4);

        % element stiffness
        for gp = 1:4 
            Fe = F{gp, ele};
            Pe = PK1(u, l, Fe);
            Be = gradN{gp};
            
            fe = fe + Pe * Be' * jac;
        end

        edof = edofMat(ele, :);
        f(edof) = f(edof) - fe(:); % NOTE: MINUS!
    end
end

function PP = PK1(u,l,F)
    % Neo hookean
    JJ = log(det(F));
    Finv = inv(F)';
    PP = u*(F-Finv) + l*JJ*Finv;

    % Linear
    % I = eye(3,3);
    % PP = u*(F + F'-2*I) + l * (trace(F)-3) * I;
end