function Ae = alok_gradgrad_2d(ke,Pe)

Re = [Pe(:,2)-Pe(:,1),Pe(:,3)-Pe(:,1)];
Be = inv(Re')*[-1 1 0;-1 0 1];
Ae = ke * Be'*Be * abs(det(Re))/2;