function CDWGB = WG_QuadriP02P02AC02_CofNmlzPiolaBas_DiscWkGradBasFxn( ... 
  QuadriMesh, GAUSSQUAD)
%% Cofs. in AC0^2 nmlz.Piola.bas.: CDWGB size NumEms*10*8 
% Discrete weak gradient for 10 WG(P0^2,P0^2) bas.fxns. on a quadri.mesh 
% James Liu, ColoState; 2017/02--2021/02 

%% Mesh info 
NumEms = QuadriMesh.NumEms; 
area = QuadriMesh.area;
xc = QuadriMesh.EmCntr(:,1);
yc = QuadriMesh.EmCntr(:,2);
% CofA = QuadriMesh.CofA;
% CofB = QuadriMesh.CofB;

%% Mesh info: More: Edge "normal" vector carrying edge length info 
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;
  NmlEg(:,j,1) = y(:,j1) - y(:,j);
  NmlEg(:,j,2) = x(:,j) - x(:,j1);
end
Xm = zeros(NumEms,4);  
Ym = zeros(NumEms,4);
for j=1:4
  j1 = mod(j,4) + 1;
  Xm(:,j) = 0.5*(x(:,j)+x(:,j1)) - xc;
  Ym(:,j) = 0.5*(y(:,j)+y(:,j1)) - yc;
end

%% JL20210207: Important! This approach implicitly use AC0 properties 
%% RHS used in the definition of discrete weak gradient 
RHS = zeros(NumEms,10,8);
RHS(:,1,3) = -2*area;  RHS(:,2,7) = -2*area;
for j=1:4
  RHS(:,2*j+1,1) = NmlEg(:,j,1);  
  RHS(:,2*j+1,2) = NmlEg(:,j,2);
  RHS(:,2*j+1,3) = Xm(:,j).*NmlEg(:,j,1) + Ym(:,j).*NmlEg(:,j,2);
  % 
  RHS(:,2*j+2,5) = NmlEg(:,j,1);
  RHS(:,2*j+2,6) = NmlEg(:,j,2);
  RHS(:,2*j+2,7) = Xm(:,j).*NmlEg(:,j,1) + Ym(:,j).*NmlEg(:,j,2);
end
RHS(:,5,4) = 1;  RHS(:,7,4) = -1;
RHS(:,6,8) = 1;  RHS(:,8,8) = -1;

% %%
% RHS = zeros(NumEms,10,8);
% RHS(:,1,3) = -2*area;  
% RHS(:,2,7) = -2*area;
% NumQuadPts = size(GAUSSQUAD.LINE,1);
% rxhat = zeros(4,NumQuadPts);
% ryhat = zeros(4,NumQuadPts);
% for k=1:NumQuadPts
%   rxhat(1,k) = GAUSSQUAD.LINE(k,2);  ryhat(1,k) = 0;  % Edge #1
%   rxhat(3,k) = GAUSSQUAD.LINE(k,2);  ryhat(3,k) = 1;  % Edge #3
%   rxhat(2,k) = 1;  ryhat(2,k) = GAUSSQUAD.LINE(k,2);  % Edge #2
%   rxhat(4,k) = 0;  ryhat(4,k) = GAUSSQUAD.LINE(k,2);  % Edge #4
% end
% for j=1:4
%   for k=1:NumQuadPts
%     FV = Hdiv_QuadriAC0_NmlzPiolaBas_FxnVal(QuadriMesh,rxhat(j,k),ryhat(j,k));
%     for i=1:4
%       RHS(:,2+2*j-1,i  ) = RHS(:,2+2*j-1,i  ) + GAUSSQUAD.LINE(k,3) * squeeze(sum(FV(:,i,:).*NmlEg(:,j,:),3));
%       RHS(:,2+2*j,  i+4) = RHS(:,2+2*j,  i+4) + GAUSSQUAD.LINE(k,3) * squeeze(sum(FV(:,i,:).*NmlEg(:,j,:),3));
%     end
%   end
% end

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

%% Solving now an 8x8 SPD system on each element 
CDWGB = zeros(NumEms,10,8);
for ie=1:NumEms
  EltGM = squeeze(GM(ie,:,:));
  for j=1:10
    EltRHS = squeeze(RHS(ie,j,:));
    CDWGB(ie,j,:) = EltGM\EltRHS;
  end 
end

return;