function GM = Hdiv_QuadriRT0_NmlzBas_GramMat(QuadriMesh,GAUSSQUAD)
%% "Hdiv" local RT[0] on quadrilateral mesh: Gram matrix for normalized basis 
% Using implicitly a 5-point Gaussian quadrature for quadrilaterals 
% 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);
x1 = QuadriMesh.node(k1,1);  y1 = QuadriMesh.node(k1,2);
x2 = QuadriMesh.node(k2,1);  y2 = QuadriMesh.node(k2,2);
x3 = QuadriMesh.node(k3,1);  y3 = QuadriMesh.node(k3,2);
x4 = QuadriMesh.node(k4,1);  y4 = QuadriMesh.node(k4,2);

%% Elementwise area and center 
area = 0.5*( (x2-x1).*(y3-y1) - (x3-x1).*(y2-y1)...
           + (x3-x1).*(y4-y1) - (x4-x1).*(y3-y1));
xc = 0.25*(x1+x2+x3+x4);
yc = 0.25*(y1+y2+y3+y4);

%% Coefficients for the bilinear mapping 
CofA = zeros(NumEms,4);       CofB = zeros(NumEms,4);
CofA(:,1) = x1;               CofB(:,1) = y1;
CofA(:,2) = x2-x1;            CofB(:,2) = y2-y1;
CofA(:,3) = x4-x1;            CofB(:,3) = y4-y1;
CofA(:,4) = (x1+x3)-(x2+x4);  CofB(:,4) = (y1+y3)-(y2+y4);

%% Computing integrals using a 5-point Gaussian quadrature on [0,1]^2
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);
  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);
end

%% Gram matrix 
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;

return;