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

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

%% Gram matrix for AC1 
GM = zeros(NumEms,10,10);
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;
  Fval = Hdiv_QuadriAC1_NmlzPiolaBas_FxnVal(QuadriMesh,xhat,yhat);
  for i=1:10
    for j=1:10
      GM(:,i,j) = GM(:,i,j) + GAUSSQUAD.RECT(k,3) * jac ...
        .* ( Fval(:,i,1).*Fval(:,j,1) + Fval(:,i,2).*Fval(:,j,2) );
    end
  end
end

return;