function CDWGB = WG_QuadriQ0Q0RT0_CofRT0NmlzBas_DiscWkGradBasFxn(QuadriMesh,GAUSSQUAD) 
%% Weak Galerkin on quadri.mesh: Discrte weak gradients for (Q0,Q0) bas.fxns.
% Cofs. in RT[0] nmlz.bas. for disc.wk.grad. of 5 WG(Q0,Q0) bas.fxns. on quadri.
% James Liu, ColoState; 2012/07--2020/08 

%% 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 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);

%% Coefficients 
CofA = zeros(NumEms,4);
CofA(:,1) = x(:,1);
CofA(:,2) = x(:,2) - x(:,1);
CofA(:,3) = x(:,4) - x(:,1);
CofA(:,4) = (x(:,1)+x(:,3)) - (x(:,2)+x(:,4));
CofB = zeros(NumEms,4);
CofB(:,1) = y(:,1);
CofB(:,2) = y(:,2) - y(:,1);
CofB(:,3) = y(:,4) - y(:,1);
CofB(:,4) = (y(:,1)+y(:,3)) - (y(:,2)+y(:,4));

%% Computing auxiliary quantities using Gaussian quadrature on [0,1]^2
SX = zeros(NumEms,1);  SY = zeros(NumEms,1);
SX2 = zeros(NumEms,1);  SY2 = zeros(NumEms,1);  % SXY = zeros(NumEms,1);
NumQuadPts = size(GAUSSQUAD.RECT,1);
GAUSSQUADRECT1 = ones(NumQuadPts,2) - GAUSSQUAD.RECT(:,1:2);
qp = zeros(NumEms,2);
for k=1:NumQuadPts
  xhat = GAUSSQUADRECT1(k,1);
  yhat = GAUSSQUADRECT1(k,2);
  qp(:,1) = CofA(:,1) + CofA(:,2)*xhat + CofA(:,3)*yhat + CofA(:,4)*xhat*yhat;
  qp(:,2) = CofB(:,1) + CofB(:,2)*xhat + CofB(:,3)*yhat + CofB(:,4)*xhat*yhat;
  jac = (CofB(:,3)+CofB(:,4)*xhat) .* (CofA(:,2)+CofA(:,4)*yhat)...
      - (CofA(:,3)+CofA(:,4)*xhat) .* (CofB(:,2)+CofB(:,4)*yhat);
  X = qp(:,1) - xc;
  Y = qp(:,2) - yc;
  SX = SX + GAUSSQUAD.RECT(k,3)*jac.*X;
  SY = SY + GAUSSQUAD.RECT(k,3)*jac.*Y;
  SX2 = SX2 + GAUSSQUAD.RECT(k,3)*jac.*(X.*X);
  SY2 = SY2 + GAUSSQUAD.RECT(k,3)*jac.*(Y.*Y);
  % SXY = SXY + GAUSSQUAD.RECT(k,3)*jac.*(X.*Y);
end

%% Gram matrices for all elements 
GM = zeros(NumEms,4,4);
GM(:,1,1) = area;  GM(:,2,2) = area;
GM(:,1,3) = SX;  GM(:,2,4) = SY;
GM(:,3,1) = SX;  GM(:,4,2) = SY;
GM(:,3,3) = SX2;  GM(:,4,4) = SY2;

%% More auxiliary quantities 
RHS = zeros(NumEms,5,4);
RHS(:,1,3) = -area;
RHS(:,1,4) = -area;
Xm = zeros(NumEms,4);  Ym = zeros(NumEms,4);  NmlEg = zeros(NumEms,4,2);
for i=1:4
  Xm(:,i) = 0.5*(x(:,i)+x(:,i+1)) - xc;
  Ym(:,i) = 0.5*(y(:,i)+y(:,i+1)) - yc;
  % Normal vector with the length of the edge 
  NmlEg(:,i,1) = y(:,i+1) - y(:,i);
  NmlEg(:,i,2) = x(:,i) - x(:,i+1);
  RHS(:,i+1,1) = NmlEg(:,i,1);
  RHS(:,i+1,2) = NmlEg(:,i,2);
  RHS(:,i+1,3) = Xm(:,i).*NmlEg(:,i,1);
  RHS(:,i+1,4) = Ym(:,i).*NmlEg(:,i,2);
end

%% Computing 4 coeffs. in RT[0] 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;