function GMST = Hdiv_QuadriAC02_NmlzPiolaBas_GramMatST(QuadriMesh,GAUSSQUAD)
%% Hdiv: AC0^2 spaces on a quadrilateral mesh: Gram matrix 
% Colon products of one basis function and the transpose of another 
% Self and Transpose 
% James Liu, Graham Harper, ColoState; 2017/02--2021/02 

%% Mesh info 
NumEms = QuadriMesh.NumEms; 
xc = QuadriMesh.EmCntr(:,1);
yc = QuadriMesh.EmCntr(:,2);
CofA = QuadriMesh.CofA;
CofB = QuadriMesh.CofB;

%% 
GMST = zeros(NumEms,8,8);
MatJ = zeros(NumEms,2,2);
MatP = zeros(NumEms,2,2);
NumQuadPts = size(GAUSSQUAD.RECT,1); 
GAUSSQUADRECT1 = ones(NumQuadPts,2) - GAUSSQUAD.RECT(:,1:2);
bas = zeros(NumEms,4,2);
bas(:,1,1) = 1;  bas(:,2,2) = 1;
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 = [CofA(:,1) + CofA(:,2)*xhat + CofA(:,3)*yhat + CofA(:,4)*xhat*yhat,...
        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);
  bas(:,3,1) = X;  bas(:,3,2) = Y;
  bas(:,4,1) = U;  bas(:,4,2) = V;
  for i=1:4
    for j=1:4
      GMST(:,i,  j  ) = GMST(:,i,  j  ) + GAUSSQUAD.RECT(k,3) * bas(:,i,1).*bas(:,j,1) .* jac;
      GMST(:,i,  j+4) = GMST(:,i,  j+4) + GAUSSQUAD.RECT(k,3) * bas(:,i,2).*bas(:,j,1) .* jac;
      GMST(:,i+4,j  ) = GMST(:,i+4,j  ) + GAUSSQUAD.RECT(k,3) * bas(:,i,1).*bas(:,j,2) .* jac;
      GMST(:,i+4,j+4) = GMST(:,i+4,j+4) + GAUSSQUAD.RECT(k,3) * bas(:,i,2).*bas(:,j,2) .* jac;
    end
  end
end

return;