function CDWGB = WG_TriP0P0RT0_CofRT0NmlzBas_DiscWkGradBasFxn(TriMesh)
%% Weak Galerkin on tri.mesh: Discrte weak gradients for (P0,P0) bas.fxns.
% Cofs. RT0 nmlz.bas. for disc.wk.grad. of 4 WG(P0,P0) bas.fxns. on triangles 
% 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 
SX2SY2 = (1/36)*area.*( (x1-x2).^2 + (x2-x3).^2 + (x3-x1).^2 ...
                      + (y1-y2).^2 + (y2-y3).^2 + (y3-y1).^2 );

%% Coeff. in RT0 nmlz.bas. for discrete weak gradients of WG(P0,P0) bas.fxns.
% See notes for mathematical simplification 
CDWGB = zeros(NumEms,4,3);  % All elements, 4 WG0 bas.fxns., 3 RT0 bas.fxns.
% For the element interior basis function 
CDWGB(:,1,:) = [zeros(NumEms,1), zeros(NumEms,1), -2*area./SX2SY2];
% For the three edge basis functions 
CDWGB(:,2,:) = [(y3-y2)./area, (x2-x3)./area, (2/3)*area./SX2SY2];
CDWGB(:,3,:) = [(y1-y3)./area, (x3-x1)./area, (2/3)*area./SX2SY2];
CDWGB(:,4,:) = [(y2-y1)./area, (x1-x2)./area, (2/3)*area./SX2SY2];

return;