function CDWGB = ... 
  WG_QuadriP1P1AC1_CofAC1NmlzPiolaBas_DiscWkGradBasFxn( ...
  QuadriMesh,GAUSSQUAD) 
%% Weak Galerkin on quadri.mesh: Discrte weak gradients: (P1,P1) bas.fxns.
% Cofs. in AC1 nmlz.Piola.bas. for disc.wk.grad. of 11 WG(P1,P1) bas.fxns.
% James Liu, ColoState; 2017/02--2017/10 

%% Mesh info 
NumEms = QuadriMesh.NumEms;
k1 = QuadriMesh.elem(:,1);  k2 = QuadriMesh.elem(:,2);  
k3 = QuadriMesh.elem(:,3);  k4 = QuadriMesh.elem(:,4);
x = zeros(NumEms,1);  y = zeros(NumEms,1);
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);
xc = QuadriMesh.EmCntr(:,1);
yc = QuadriMesh.EmCntr(:,2);
CofA = QuadriMesh.CofA;
CofB = QuadriMesh.CofB;

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

%% Auxiliary quantities 
MatJ = zeros(NumEms,2,2);
S1 = zeros(NumEms,1);  SXY = zeros(NumEms,1);
SX = zeros(NumEms,1);  SX2 = zeros(NumEms,1);
SY = zeros(NumEms,1);  SY2 = 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);
  MatJ(:,1,1) = CofA(:,2) + CofA(:,4)*yhat;  
  MatJ(:,1,2) = CofA(:,3) + CofA(:,4)*xhat;
  MatJ(:,2,1) = CofB(:,2) + CofB(:,4)*yhat;
  MatJ(:,2,2) = CofB(:,3) + CofB(:,4)*xhat;
  jac = MatJ(:,1,1).*MatJ(:,2,2) - MatJ(:,1,2).*MatJ(:,2,1);
  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;
  X = qp(:,1) - xc;
  Y = qp(:,2) - yc;
  S1 = S1 + GAUSSQUAD.RECT(k,3) * jac;
  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

%% RHS 
RHS = zeros(NumEms,11,10);
% For interior WG basis functions 
RHS(:,1,2) = -S1;    RHS(:,2,2) = -SX;     RHS(:,3,2) = -SY;
RHS(:,1,6) = -S1;    RHS(:,2,6) = -SX;     RHS(:,3,6) = -SY;
RHS(:,1,7) = -3*SX;  RHS(:,2,7) = -3*SX2;  RHS(:,3,7) = -3*SXY;
RHS(:,1,8) = -3*SY;  RHS(:,2,8) = -3*SXY;  RHS(:,3,8) = -3*SY2;
% For edge WG basis functions 
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
    r = GAUSSQUAD.LINE(k,2) - 0.5;
    FV = Hdiv_QuadriAC1_NmlzPiolaBas_FxnVal(QuadriMesh,rxhat(j,k),ryhat(j,k));
    for i=1:10 
      RHS(:,3+2*j-1,i) = RHS(:,3+2*j-1,i) + GAUSSQUAD.LINE(k,3) * sum(FV(:,i,:).*NmlEg(:,j,:),3)    ; 
      RHS(:,3+2*j,  i) = RHS(:,3+2*j,  i) + GAUSSQUAD.LINE(k,3) * sum(FV(:,i,:).*NmlEg(:,j,:),3) * r; 
    end
  end
end

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

%% Computing 10 coeffs. in AC1 basis for disc.wk.grad. of 11 WG bas.fxns. 
CDWGB = zeros(NumEms,11,10);
% 3 WG bas.fxns. for element interior 
% 2 WG bas.fxns. for each of the 4 edges 
for ie=1:NumEms
  EltGM = squeeze(GM(ie,:,:));
  for j=1:11
    CDWGB(ie,j,:) = EltGM \ squeeze(RHS(ie,j,:));
  end
end

return;