function NmlFlux = Hdiv_TriRT0_NmlzBas_NmlFlux(TriMesh)
%% "Hdiv" local RT0 on tri.mesh: Normal fluxes for 3 nmlz.bas.fxns. & 3 edges 
% [1;0],[0;1],[X;Y] where X=x-xc, Y=y-yc, (xc,yc) is 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));
% area = TriMesh.area;

%% All elements 3 nmlz.bas.fxns. 3 edge bulk outward normal fluxes 
NmlFlux = zeros(NumEms,3,3);
% For the basis function [1;0]
NmlFlux(:,1,1) = y3 - y2;
NmlFlux(:,1,2) = y1 - y3;
NmlFlux(:,1,3) = y2 - y1;
% For the basis function [0;1]
NmlFlux(:,2,1) = x2 - x3;
NmlFlux(:,2,2) = x3 - x1;
NmlFlux(:,2,3) = x1 - x2;
% For the basis function [X;Y]
NmlFlux(:,3,1) = (2.0/3)*area;
NmlFlux(:,3,2) = (2.0/3)*area;
NmlFlux(:,3,3) = (2.0/3)*area;

return;