function [CDWGBQ,CDWGBT] = WG_HybridPQ0_CofRT0NmlzBas_DiscWkGradBasFxn(HybridMesh) 
%% Weak Galerkin on a hybrid mesh: Discrte weak gradients for P0/Q0 bas.fxns.
% CDWGBQ: Cofs. in RT[0] nmlz.bas. for disc.wk.grad. of 5 WG0 bas.fxns. on quadri.
% CDWGBT: Cofs. in RT0 nmlz.bas. for disc.wk.grad. of 4 WG0 bas.fxns. on trig.
% James Liu, Zhuoran Wang, ColoState; 2012/07--2017/02 

%% For the quadrilateral elements 
NumEmsQ = HybridMesh.NumEmsQ;
k1 = HybridMesh.elemQ(:,1);  k2 = HybridMesh.elemQ(:,2);  
k3 = HybridMesh.elemQ(:,3);  k4 = HybridMesh.elemQ(:,4);
x = zeros(NumEmsQ,5);  y = zeros(NumEmsQ,5);
x(:,1) = HybridMesh.node(k1,1);  y(:,1) = HybridMesh.node(k1,2);
x(:,2) = HybridMesh.node(k2,1);  y(:,2) = HybridMesh.node(k2,2);
x(:,3) = HybridMesh.node(k3,1);  y(:,3) = HybridMesh.node(k3,2);
x(:,4) = HybridMesh.node(k4,1);  y(:,4) = HybridMesh.node(k4,2);
x(:,5) = x(:,1);                 y(:,5) = y(:,1);
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);
CofA = zeros(NumEmsQ,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(NumEmsQ,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 a 5-pt Gaussian quadrature on [0,1]^2
SX = zeros(NumEmsQ,1);  SY = zeros(NumEmsQ,1);
SX2 = zeros(NumEmsQ,1);  SY2 = zeros(NumEmsQ,1);  % SXY = zeros(NumEmsQ,1);
GAUSSQUADRECT = SetGaussQuadRect(5,5);
NumQuadPts = size(GAUSSQUADRECT,1);
GAUSSQUADRECT1 = ones(NumQuadPts,2) - GAUSSQUADRECT(:,1:2);
qp = zeros(NumEmsQ,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 + GAUSSQUADRECT(k,3)*jac.*X;
  SY = SY + GAUSSQUADRECT(k,3)*jac.*Y;
  SX2 = SX2 + GAUSSQUADRECT(k,3)*jac.*(X.*X);
  SY2 = SY2 + GAUSSQUADRECT(k,3)*jac.*(Y.*Y);
  % SXY = SXY + GAUSSQUADRECT(k,3)*jac.*(X.*Y);
end
% Gram matrices for all elements 
GM = zeros(NumEmsQ,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(NumEmsQ,5,4);
RHS(:,1,3) = -area;
RHS(:,1,4) = -area;
Xm = zeros(NumEmsQ,4);  Ym = zeros(NumEmsQ,4);  NmlEg = zeros(NumEmsQ,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. 
CDWGBQ = zeros(NumEmsQ,5,4);
% 1 WG bas.fxn. for element interior 
% 4 WG bas.fxns. for 4 edges 
for ie=1:NumEmsQ
  EltGM = squeeze(GM(ie,:,:));
  for j=1:5
    CDWGBQ(ie,j,:) = EltGM \ squeeze(RHS(ie,j,:));
  end
end

%% For the triangular elements 
NumEmsT = HybridMesh.NumEmsT;
k1 = HybridMesh.elemT(:,1);  
k2 = HybridMesh.elemT(:,2);  
k3 = HybridMesh.elemT(:,3);
x1 = HybridMesh.node(k1,1);  y1 = HybridMesh.node(k1,2); 
x2 = HybridMesh.node(k2,1);  y2 = HybridMesh.node(k2,2);
x3 = HybridMesh.node(k3,1);  y3 = HybridMesh.node(k3,2);
area = 0.5*abs((x2-x1).*(y3-y1)-(x3-x1).*(y2-y1));
% Auxiliary integrals 
SX2SY2 = (1/36)*area.*( (x1-x2).^2 + (x2-x3).^2 + (x3-x1).^2 ...
                      + (y1-y2).^2 + (y2-y3).^2 + (y3-y1).^2 );
% Coeff. in RT0 nmlz.bas. for discrete weak gradients of WG0 bas.fxns.
% See notes for mathematical simplification 
% NOTE: CDGWBT has size NumEms*4*3, All elts., 4 WG0 bas.fxns., 3 RT0 bas.fxns.
% For the element interior basis functions 
CDWGBT(:,1,:) = [zeros(NumEmsT,1), zeros(NumEmsT,1), -2*area./SX2SY2];
% For all three edge basis functions 
CDWGBT(:,2,:) = [(y3-y2)./area, (x2-x3)./area, (2/3)*area./SX2SY2];
CDWGBT(:,3,:) = [(y1-y3)./area, (x3-x1)./area, (2/3)*area./SX2SY2];
CDWGBT(:,4,:) = [(y2-y1)./area, (x1-x2)./area, (2/3)*area./SX2SY2];

return;