function ProjCof = Hdiv_TriRT0_NmlzBas_ProjMatK(TriMesh, MatK)
%% "Hdiv" local RT0 on tri.mesh: Projection of normalized bas.fxns. 
% Projection coefficients of (MatK * nmlz.bas.fxns.) back in nmlz.bas.fxns. 
% Assuming MatK is an elementwise constant 2x2 matrix 
% ProjCof has the size #Elements*3*3 (in the order [1;0],[0;1],[X;Y])
% 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 quantities 
% JL20160101: Auxiliary matrix [2 -1 -1; -1 2 -1; -1 -1 2] could be used 
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));
C3 = ( MatK(:,1,1).*SX2 + MatK(:,1,2).*SXY ... 
     + MatK(:,2,1).*SXY + MatK(:,2,2).*SY2)./(SX2+SY2);

%% Projection coefficients 
ProjCof = zeros(NumEms,3,3);
ZO = zeros(NumEms,1);
ProjCof(:,1,:) = [MatK(:,1,1), MatK(:,2,1), ZO];
ProjCof(:,2,:) = [MatK(:,1,2), MatK(:,2,2), ZO];
ProjCof(:,3,:) = [ZO, ZO, C3];

return;