function CDWDB = WG_QuadriP02P02P0_CofNmlzBas_DiscWkDivBasFxn( ...
  QuadriMesh, GAUSSQUAD)
%% Weak Galerkin: Coefficients in P0 normalized basis function (1) 
% for discrete weak divergence of 10 WG(P0^2,P0^2) basis functions 
% on a quadrilateral mesh 
% James Liu, ColoState; 2017/07--2021/02 

%% Mesh info 
NumEms = QuadriMesh.NumEms;
  area = QuadriMesh.area;

%% Mesh info: More 
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);
NmlEg = zeros(NumEms,4,2);
for j=1:4
  j1 = mod(j,4) + 1;
  tmpx = x(:,j1) - x(:,j);
  tmpy = y(:,j1) - y(:,j);
  NmlEg(:,j,1) =  tmpy;  % Carrying length info 
  NmlEg(:,j,2) = -tmpx;
end

%% Computing RHS needed for 10 WG bas.fxns.
RHS = zeros(NumEms,10);
% Nonzero only for the WG basis functions on the 4 edges 
for j=1:4
  RHS(:,2+2*j-1) = NmlEg(:,j,1);
  RHS(:,2+2*j  ) = NmlEg(:,j,2);
end

%% Computing coeff. in P0 nmlz.bas. for disc.wk.div. of 10 WG bas.fxns.
CDWDB = zeros(NumEms,10);
for j=1:10
  CDWDB(:,j) = RHS(:,j)./area;
end

return;