function CDWGB = WG_QuadriP0P0AC0_CofAC0NmlzPiolaBas_DiscWkGradBasFxn(...
  QuadriMesh, GAUSSQUAD) 
%% Weak Galerkin on quadri.mesh: Discrte weak gradients: (P0,P0) bas.fxns.
% Cofs. in AC0 nmlz.Piola.bas. for disc.wk.grad. of 5 WG(P0,P0) bas.fxns.
% James Liu, ColoState; 2017/02--2020/02 

%% Mesh info 
NumEms = QuadriMesh.NumEms;
k1 = QuadriMesh.elem(:,1);  k2 = QuadriMesh.elem(:,2);  
k3 = QuadriMesh.elem(:,3);  k4 = QuadriMesh.elem(:,4);
x = zeros(NumEms,4);  y = zeros(NumEms,4);
x(:,1) = QuadriMesh.node(k1,1);  y(:,1) = QuadriMesh.node(k1,2);
x(:,2) = QuadriMesh.node(k2,1);  y(:,2) = QuadriMesh.node(k2,2);
x(:,3) = QuadriMesh.node(k3,1);  y(:,3) = QuadriMesh.node(k3,2);
x(:,4) = QuadriMesh.node(k4,1);  y(:,4) = QuadriMesh.node(k4,2);
% x(:,5) = x(:,1);  y(:,5) = y(:,1);

%% Elementwise area and center 
% area = 0.5*( (x(:,2)-x(:,1)).*(y(:,3)-y(:,1)) - (x(:,3)-x(:,1)).*(y(:,2)-y(:,1))...
%            + (x(:,3)-x(:,1)).*(y(:,4)-y(:,1)) - (x(:,4)-x(:,1)).*(y(:,3)-y(:,1)));
% xc = 0.25*sum(x(:,1:4),2);
% yc = 0.25*sum(y(:,1:4),2);
area = QuadriMesh.area;
xc = QuadriMesh.EmCntr(:,1);
yc = QuadriMesh.EmCntr(:,2);

%% Normal vectors carrying edge lengths
LenEg = zeros(NumEms,4);
NmlEg = zeros(NumEms,4,2);
for j=1:4
  j1 = mod(j,4) + 1;
  tmpx = x(:,j1) - x(:,j);
  tmpy = y(:,j1) - y(:,j);
  LenEg(:,j) = sqrt(tmpx.^2+tmpy.^2);
  NmlEg(:,j,1) =  tmpy;
  NmlEg(:,j,2) = -tmpx;
end

%% RHS 
RHS = zeros(NumEms,5,4);
RHS(:,1,3) = -2*area;
for j=1:4
  j1 = mod(j,4) + 1;
  Xm = 0.5*(x(:,j)+x(:,j1)) - xc;
  Ym = 0.5*(y(:,j)+y(:,j1)) - yc;
  RHS(:,j+1,1) = NmlEg(:,j,1);
  RHS(:,j+1,2) = NmlEg(:,j,2);
  RHS(:,j+1,3) = Xm.*NmlEg(:,j,1) + Ym.*NmlEg(:,j,2);
end
RHS(:,3,4) =  1;  % Edge #2 
RHS(:,4,4) = -1;  % Edge #3 

%% Gram matrices for all elements 
GM = Hdiv_QuadriAC0_NmlzPiolaBas_GramMat(QuadriMesh,GAUSSQUAD);

%% Computing 4 coeffs. in AC0 basis for disc.wk.grad. of 5 WG bas.fxns. 
CDWGB = zeros(NumEms,5,4);
% 1 WG bas.fxn. for element interior 
% 4 WG bas.fxns. for 4 edges 
for ie=1:NumEms
  EltGM = squeeze(GM(ie,:,:));
  for j=1:5
    CDWGB(ie,j,:) = EltGM \ squeeze(RHS(ie,j,:));
  end
end

return;