function CDWGB = WG_TriP0P1BDM1_CofBDM1EgBas_DiscWkGradBasFxn(TriMesh) 
%% WG(P0,P1;BDM1) on triangular mesh: 
% Computing coeffs. in BDM1 eg.bas. for disc.wk.grad. of WG(P0,P1) bas.fxns.
% 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 
eg1 = [x3,y3] - [x2,y2];
eg2 = [x1,y1] - [x3,y3];
eg3 = [x2,y2] - [x1,y1];
lg1 = sqrt(eg1(:,1).^2+eg1(:,2).^2);
lg2 = sqrt(eg2(:,1).^2+eg2(:,2).^2);
lg3 = sqrt(eg3(:,1).^2+eg3(:,2).^2);
tan = zeros(NumEms,3,2);
tan(:,1,:) = eg1./[lg1,lg1];
tan(:,2,:) = eg2./[lg2,lg2];
tan(:,3,:) = eg3./[lg3,lg3];
% cosine theta 
ct1 = (lg2.^2+lg3.^2-lg1.^2)./(2*lg2.*lg3);
ct2 = (lg3.^2+lg1.^2-lg2.^2)./(2*lg3.*lg1);
ct3 = (lg1.^2+lg2.^2-lg3.^2)./(2*lg1.*lg2);
% sine theta 
st1 = sqrt(1-ct1.^2);
st2 = sqrt(1-ct2.^2);
st3 = sqrt(1-ct3.^2);

%% Gradients of the Lagrangian P1 bas.fxns. (= Barycentric coordinates) 
% NOTE: CP1BasGrad has size NumEms*3*2 
CP1BasGrad = CG_TriP1_BasFxnGrad(TriMesh);

% Divergence of the BDM1 bas.fxns.
DivBDM1Bas = zeros(NumEms,6);  % all elts., 3 edges, 3 ends 
DivBDM1Bas(:,1) =  dot(CP1BasGrad(:,2,:),tan(:,3,:),3)./st2;
DivBDM1Bas(:,2) = -dot(CP1BasGrad(:,3,:),tan(:,2,:),3)./st3;
DivBDM1Bas(:,3) =  dot(CP1BasGrad(:,3,:),tan(:,1,:),3)./st3;
DivBDM1Bas(:,4) = -dot(CP1BasGrad(:,1,:),tan(:,3,:),3)./st1;
DivBDM1Bas(:,5) =  dot(CP1BasGrad(:,1,:),tan(:,2,:),3)./st1;
DivBDM1Bas(:,6) = -dot(CP1BasGrad(:,2,:),tan(:,1,:),3)./st2;

%% Gram matrices of BDM1 eg.bas.bas.fxns for a triangular mesh 
GMBDM1 = Hdiv_TriBDM1_EgBas_GramMat(TriMesh);

%% "RHS" 
RHS = zeros(NumEms,6,7); 
% For all six edge basis functions 
RHS(:,1:2,1) = [2*lg1, lg1]/6;
RHS(:,1:2,2) = [lg1, 2*lg1]/6;
RHS(:,3:4,3) = [2*lg2, lg2]/6;
RHS(:,3:4,4) = [lg2, 2*lg2]/6;
RHS(:,5:6,5) = [2*lg3, lg3]/6;
RHS(:,5:6,6) = [lg3, 2*lg3]/6;
% For one element interior basis functions 
for j=1:6 
  RHS(:,j,7) = -DivBDM1Bas(:,j).*area; 
end

%% All coeffs. for the discrete weak gradients in BDM1 eg.bas.bas.fxns. 
% See notes for mathematical simplification 
CDWGB = zeros(NumEms,6,7);  % all elts., 6 BDM1 bas.fxns., 7 WG bas.fxns. 
for ie=1:NumEms
  EltGM = squeeze(GMBDM1(ie,1:6,1:6));
  for k=1:7
    EltRHS = RHS(ie,1:6,k)';
    CDWGB(ie,1:6,k) = (EltGM\EltRHS)';
  end
end

return;