function EGGM = CG_TriP1_GradGradMat(TriMesh)
%% CG: Tri.mesh. P1 elements: Elementwise grad-grad matrix (for all elements)
% James Liu, ColoState; 2012/07--2017/02 

%% Mesh info 
NumEms = TriMesh.NumEms;
k1 = TriMesh.elem(:,1);  k2 = TriMesh.elem(:,2);  k3 = TriMesh.elem(:,3);
x1 = TriMesh.node(k1,1);  y1 = TriMesh.node(k1,2); 
x2 = TriMesh.node(k2,1);  y2 = TriMesh.node(k2,2);
x3 = TriMesh.node(k3,1);  y3 = TriMesh.node(k3,2);
area = 0.5*((x2-x1).*(y3-y1)-(x3-x1).*(y2-y1));

%% NOTE: Lagrangian P1 bas.fxn.val. = Bary crd.
P1BasGrad = zeros(NumEms,3,2);
P1BasGrad(:,1,:) = [(y2-y3)./(2*area), (x3-x2)./(2*area)];
P1BasGrad(:,2,:) = [(y3-y1)./(2*area), (x1-x3)./(2*area)];
P1BasGrad(:,3,:) = [(y1-y2)./(2*area), (x2-x1)./(2*area)];

%% Setting up all entries
EGGM = zeros(NumEms,3,3);
for i=1:3
  for j=1:3
    EGGM(:,i,j) = squeeze(dot(P1BasGrad(:,i,:),P1BasGrad(:,j,:),3)).* area;
  end
end

return;