function GMC = CG_TriP1_GlbMatConv(TriMesh,CofRT0EgBas)
%% TriMesh: Continuous Galerkin, Lagrange P1, Global mat. for convection prob.
% For all elements 
% James Liu, ColoState; 2012/07--2017/02 

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

%% 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];
LenEg = zeros(NumEms,3);
LenEg(:,1) = sqrt((x3-x2).^2+(y3-y2).^2);
LenEg(:,2) = sqrt((x1-x3).^2+(y1-y3).^2);
LenEg(:,3) = sqrt((x2-x1).^2+(y2-y1).^2);

%% 
ArrayA = zeros(NumEms,3,9);
ZO = zeros(NumEms,1);
ArrayA(:,1,:) = [ZO,2*area,2*area, -2*area,ZO,ZO, -2*area,ZO,ZO];
ArrayA(:,2,:) = [ZO,-2*area,ZO, 2*area,ZO,2*area, ZO,-2*area,ZO];
ArrayA(:,3,:) = [ZO,ZO,-2*area, ZO,ZO,-2*area, 2*area,2*area,ZO];
% 
ArrayB = zeros(NumEms,3,9);
ArrayB(:,1,:) = [2*area,area,area, 2*area area area, 2*area area area]/12;
ArrayB(:,2,:) = [area,2*area,area, area 2*area area, area 2*area area]/12;
ArrayB(:,3,:) = [area,area,2*area, area area 2*area, area area 2*area]/12;
% 
ArrayC = zeros(NumEms,3,9);
for k=1:3
  for l=1:3
    ArrayC(:,1,3*(k-1)+l) = -0.25./(area.^2) .* LenEg(:,k) .* CofRT0EgBas(:,k);
    ArrayC(:,2,3*(k-1)+l) = -0.25./(area.^2) .* LenEg(:,k) .* CofRT0EgBas(:,k);
    ArrayC(:,3,3*(k-1)+l) = -0.25./(area.^2) .* LenEg(:,k) .* CofRT0EgBas(:,k);
  end
end

%% JL20151119: REDUNDANCE IN .*area, ./area
EltMat = (ArrayA .* ArrayB) .* ArrayC;

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

return;