function GM = Hdiv_TriBDM1_NmlzBas_GramMat(TriMesh)
%% "Hdiv" local BDM1 on a triangular mesh: Gram matrces for the normalized basis
% [1;0], [0;1], [X;0], [0;X], [Y;0], [0;Y] 
% where X=x-xc, Y=y-yc, (xc,yc) is the element center 
% 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*abs((x2-x1).*(y3-y1)-(x3-x1).*(y2-y1));

%% 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).*(2*(x2-x1).*(y2-y1) - (x2-x1).*(y3-y1) ...
                  - (x3-x1).*(y2-y1) + 2*(x3-x1).*(y3-y1) );

%% Gram matrices for all elements in the triangular mesh 
% For each element, Gram matrix is a 6x6 mat. 
GM = zeros(NumEms,6,6);
GM(:,1,1) = area;  GM(:,2,2) = area;
GM(:,3,3) = SX2;   GM(:,4,4) = SX2;
GM(:,5,5) = SY2;   GM(:,6,6) = SY2;
GM(:,3,5) = SXY;   GM(:,5,3) = SXY;
GM(:,4,6) = SXY;   GM(:,6,4) = SXY;

return;