function GM = Hdiv_QuadriAC0_NmlzPiolaBas_GramMat(QuadriMesh,GAUSSQUAD)
%% Quadrilateral mesh: 
% James Liu, ColoState; 2017/02--2020/09 

%% Mesh info 
NumEms = QuadriMesh.NumEms; 
CofA = QuadriMesh.CofA;
CofB = QuadriMesh.CofB;

%% Computing GramMat 
GM = zeros(NumEms,4,4);
NumQuadPts = size(GAUSSQUAD.RECT,1);
GAUSSQUADRECT1 = ones(NumQuadPts,2) - GAUSSQUAD.RECT(:,1:2);
for k=1:NumQuadPts
  xhat = GAUSSQUADRECT1(k,1);
  yhat = GAUSSQUADRECT1(k,2);
  J11 = CofA(:,2) + CofA(:,4)*yhat;  J12 = CofA(:,3) + CofA(:,4)*xhat;
  J21 = CofB(:,2) + CofB(:,4)*yhat;  J22 = CofB(:,3) + CofB(:,4)*xhat;
  jac = J11.*J22 - J12.*J21;
  FV = Hdiv_QuadriAC0_NmlzPiolaBas_FxnVal(QuadriMesh,xhat,yhat);
  for i=1:4
    for j=1:4
      GM(:,i,j) = GM(:,i,j) + (FV(:,i,1).*FV(:,j,1) + FV(:,i,2).*FV(:,j,2)) .* jac * GAUSSQUAD.RECT(k,3);
    end
  end
end

% %% Computing auxiliary quantities using a default Gaussian quadrature 
% MatJ = zeros(NumEms,2,2);
% MatP = zeros(NumEms,2,2);
% S1 = zeros(NumEms,1);
% SX  = zeros(NumEms,1);  SY  = zeros(NumEms,1);
% SU  = zeros(NumEms,1);  SV  = zeros(NumEms,1);
% SX2 = zeros(NumEms,1);  SY2 = zeros(NumEms,1);
% SXUYV = zeros(NumEms,1);  SU2V2 = 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);
%   J11 = CofA(:,2) + CofA(:,4)*yhat;  J12 = CofA(:,3) + CofA(:,4)*xhat;
%   J21 = CofB(:,2) + CofB(:,4)*yhat;  J22 = CofB(:,3) + CofB(:,4)*xhat;
%   jac = J11.*J22 - J12.*J21;
%   MatJ(:,1,1) = J11;  MatJ(:,1,2) = J12;
%   MatJ(:,2,1) = J21;  MatJ(:,2,2) = J22;
%   for i=1:2
%     for j=1:2
%       MatP(:,i,j) = MatJ(:,i,j)./jac;
%     end
%   end
%   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;
%   U = MatP(:,1,1) .* xhat - MatP(:,1,2) .* yhat;
%   V = MatP(:,2,1) .* xhat - MatP(:,2,2) .* yhat;
%   S1 = S1 + GAUSSQUAD.RECT(k,3) * jac;
%   SX = SX + GAUSSQUAD.RECT(k,3) * jac .* X;
%   SY = SY + GAUSSQUAD.RECT(k,3) * jac .* Y;
%   SU = SU + GAUSSQUAD.RECT(k,3) * jac .* U;
%   SV = SV + GAUSSQUAD.RECT(k,3) * jac .* V;
%   SX2 = SX2 + GAUSSQUAD.RECT(k,3) * jac .* (X.^2);
%   SY2 = SY2 + GAUSSQUAD.RECT(k,3) * jac .* (Y.^2);
%   SXUYV = SXUYV + GAUSSQUAD.RECT(k,3) * jac .* (X.*U+Y.*V);
%   SU2V2 = SU2V2 + GAUSSQUAD.RECT(k,3) * jac .* (U.^2+V.^2);
% end
% 
% %% Gram matrix for AC0 
% GM = zeros(NumEms,4,4);
% % The part for normalized basis  
% GM(:,1,1) = S1;  GM(:,2,2) = S1;  
% GM(:,1,3) = SX;  GM(:,3,1) = SX;
% GM(:,2,3) = SY;  GM(:,3,2) = SY;
% GM(:,3,3) = SX2 + SY2;
% % The part involving Piola basis 
% GM(:,1,4) = SU;  GM(:,4,1) = SU;
% GM(:,2,4) = SV;  GM(:,4,2) = SV;
% GM(:,3,4) = SXUYV;  GM(:,4,3) = SXUYV;
% GM(:,4,4) = SU2V2;

return;