function GMD = CG_TriP1_GlbMatDiff(TriMesh)
%% TriMesh: Continuous Galerkin, Lagrange P1, Global mat. for diffusion problem
% For all elements 
% James Liu, ColoState; 2012/07--2017/02 

%% Mesh info 
NumEms = TriMesh.NumEms;
NumNds = TriMesh.NumNds;
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));

%% Edge quantities 
edge = zeros(NumEms,3,2);
edge(:,1,:) = [x3-x2, y3-y2];  
edge(:,2,:) = [x1-x3, y1-y3];  
edge(:,3,:) = [x2-x1, y2-y1];

%% 
EltMat = zeros(NumEms,3,3);
for i=1:3
  for j=1:3
    EltMat(:,i,j) = squeeze(sum(edge(:,i,:).*edge(:,j,:),3)) ./ (4*area);
  end
end

%% Global matrix for diffusion 
GMD = sparse(NumNds,NumNds);
for i=1:3
  II = TriMesh.elem(:,i);
  for j=1:3
    JJ = TriMesh.elem(:,j);
    GMD = GMD + sparse(II,JJ,EltMat(:,i,j),NumNds,NumNds);
  end
end

return;