function GMK = Hdiv_QuadriAC1_NmlzPiolaBas_GramMatK(...
  QuadriMesh,MatK,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 
GMK = zeros(NumEms,10,10);
MatJ = zeros(NumEms,2,2);
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);
  MatJ(:,1,1) = CofA(:,2) + CofA(:,4)*yhat;  
  MatJ(:,1,2) = CofA(:,3) + CofA(:,4)*xhat;
  MatJ(:,2,1) = CofB(:,2) + CofB(:,4)*yhat;
  MatJ(:,2,2) = CofB(:,3) + CofB(:,4)*xhat;
  jac = MatJ(:,1,1).*MatJ(:,2,2) - MatJ(:,1,2).*MatJ(:,2,1);
  FV = Hdiv_QuadriAC1_NmlzPiolaBas_FxnVal(QuadriMesh,xhat,yhat);
  KV = zeros(NumEms,10,2);
  for j=1:10
    KV(:,j,1) = MatK(:,1,1).*FV(:,j,1) + MatK(:,1,2).*FV(:,j,2);
    KV(:,j,2) = MatK(:,2,1).*FV(:,j,1) + MatK(:,2,2).*FV(:,j,2);
  end
  for i=1:10
    for j=1:10
      GMK(:,i,j) = GMK(:,i,j) + GAUSSQUAD.RECT(k,3) * jac ...
        .* ( FV(:,i,1).*KV(:,j,1) + FV(:,i,2).*KV(:,j,2) );
    end
  end
end

return;