function CDWGB = WG_QuadriP1P0P02_CofNatBas_DiscWkGradBasFxn(QuadriMesh) 
%% Cofs. in P0^2 natural basis for disc.wk.grad. of WG(P1,P0)Quadri. 7 bas.fxns.
% James Liu, ColoState; 2012/07--2017/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,5);  y = zeros(NumEms,5);
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 outward normal vectors carrying edge-length 
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)));
NmlEg = zeros(NumEms,4,2);
for j=1:4
  NmlEg(:,j,:) = [y(:,j+1)-y(:,j), x(:,j)-x(:,j+1)];  % Carrying edge length 
end

%% Discrete weak gradients of WG(P1,P0)Quadri. 7 basis functions 
CDWGB = zeros(NumEms,7,2);
for j=1:4
  CDWGB(:,3+j,:) = squeeze(NmlEg(:,j,:))./[area,area];
end

return;