Finite Element assembly

Question

I'm having serious problems in a simple example of fem assembly. I just want to assemble the Mass matrix without any coefficient. The geometry is simple:

conn=[1, 2, 3];
p = [0 0; 1 0; 0 1];

I made it like this so that the physical element will be equal to the reference one.

my basis functions:

phi_1 = @(eta) 1 - eta(1) - eta(2);
phi_2 = @(eta) eta(1);
phi_3 = @(eta) eta(2);
phi = {phi_1, phi_2, phi_3};

Jacobian matrix:

J = @(x,y) [x(2) - x(1), x(3) - x(1);
            y(2) - y(1), y(3) - y(1)];

The rest of the code:

M = zeros(np,np);
for K = 1:size(conn,1)
    l2g = conn(K,:); %local to global mapping 
    x = p(l2g,1); %node x-coordinate
    y = p(l2g,2); %node y-coordinate
    jac = J(x,y);
    loc_M = localM(jac, phi);
    M(l2g, l2g) = M(l2g, l2g) + loc_M; %add element masses to M
end

localM:

function loc_M = localM(J,phi)
    d_J = det(J);
    loc_M = zeros(3,3);
    for i = 1:3
        for j = 1:3

            loc_M(i,j) = d_J * quadrature(phi{i}, phi{j});

        end
    end

end

quadrature:

function value = quadrature(phi_i, phi_j)

    p = [1/3, 1/3;
         0.6, 0.2;
         0.2, 0.6;                                    
         0.2, 0.2];

    w = [-27/96, 25/96, 25/96, 25/96];

    res = 0;
    for i = 1:size(p,1)
        res = res + phi_i(p(i,:)) * phi_j(p(i,:))  * w(i);
    end

    value = res;  

end

For the simple entry (1,1) I obtain 0.833, while computing the integral by hand or on wolfram alpha I get 0.166 (2 times the result of the quadrature). I tried with different points and weights for quadrature, but really I do not know what I am doing wrong.