function GM = Hdiv_TriRT0_EgBas_GramMat(TriMesh)
%% "Hdiv" local RT0 on tri.mesh: Gram matrices for edge-based basis 
% JL20131102: TO BE FURTHER REVISED FOR EFFICIENCY 
% James Liu, ColoState; 2012/07--2017/02 

%% Mesh info 
NumEms = TriMesh.NumEms;
% NumEgs = TriMesh.NumEgs;
% Assuming TriMesh.flag>=2 
e2g = TriMesh.elem2edge;
area = TriMesh.area;
LenEg = TriMesh.LenEg;

%% Gram matrices: locally for individual elements 
% Assuming vrtx is a 2-by-3 matrix 
GM = zeros(NumEms,3,3);
A = kron(ones(3,3)+eye(3),eye(2));
for ie=1:NumEms
  vrtx = (TriMesh.node(TriMesh.elem(ie,1:3),:))';
  B = [vrtx(:),vrtx(:),vrtx(:)] - [vrtx;vrtx;vrtx];
  C = diag(LenEg(e2g(ie,1:3)));
  GM(ie,1:3,1:3) = (1/48)/area(ie)*(C'*B'*A*B*C);
end

return;