function GMK = Hdiv_TriRT0_NmlzBas_GramMatK(TriMesh,MatK)
%% "Hdiv" local RT0 on a triangular mesh: Gram mat. for MatK for nmlz.bas. 
% Assuming MatK is an elementwise constant 2x2 matrix 
% [1;0], [0;1], [X;Y] where X=x-xc, Y=y-yc and (xc,yc) is element center 
% James Liu, ColoState; 2012/07--2016/06

%% 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*abs((x2-x1).*(y3-y1)-(x3-x1).*(y2-y1));
% xc = (1/3)*(x1+x2+x3);
% yc = (1/3)*(y1+y2+y3);

%% Auxiliary integrals 
SX2 = (area/36).*( (x1-x2).^2 + (x2-x3).^2 + (x3-x1).^2 );
SY2 = (area/36).*( (y1-y2).^2 + (y2-y3).^2 + (y3-y1).^2 );
SXY = (area/36).*((x1-x2).*(y1-y2) + (x2-x3).*(y2-y3) + (x3-x1).*(y3-y1));
% % Applying a special Gaussian quadrature on triangles 
% GAUSSQUADTRIG = zeros(3,4);
% GAUSSQUADTRIG(:,1:3) = [4,1,1; 1,4,1; 1,1,4]/6;
% GAUSSQUADTRIG(:,4) = [1; 1; 1]/3;
% SXY = zeros(NumEms,1);
% NumQuadPts = size(GAUSSQUADTRIG,1);
% for k=1:NumQuadPts
%   qp = GAUSSQUADTRIG(k,1)*[x1,y1]...
%      + GAUSSQUADTRIG(k,2)*[x2,y2]...
%      + GAUSSQUADTRIG(k,3)*[x3,y3];
%   X = qp(:,1) - xc;  
%   Y = qp(:,2) - yc;
%   SXY = SXY + GAUSSQUADTRIG(k,4)*(X.*Y).*area;
% end

%% Gram matrices for all elements in the triangular mesh 
% For each element, Gram matrix is a 3x3 mat 
GMK = zeros(NumEms,3,3);
GMK(:,1,1) = MatK(:,1,1) .* area;  
GMK(:,1,2) = MatK(:,2,1) .* area;
GMK(:,2,1) = MatK(:,1,2) .* area;
GMK(:,2,2) = MatK(:,2,2) .* area;
GMK(:,3,3) = MatK(:,1,1).*SX2 + MatK(:,1,2).*SXY...
           + MatK(:,2,1).*SXY + MatK(:,2,2).*SY2;

return;