function GMK = Hdiv_QuadriAC0_NmlzPiolaBas_GramMatK( ... 
  QuadriMesh, MatK, GAUSSQUAD)
%% Quadrilateral mesh: 
% James Liu, ColoState; 2017/02--2020/09 

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

%% Computing Gram matrix with weight MatK for AC0 using Gaussian quadrature 
GMK = 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
      tmp1 = MatK(:,1,1).*FV(:,j,1) + MatK(:,1,2).*FV(:,j,2);
      tmp2 = MatK(:,2,1).*FV(:,j,1) + MatK(:,2,2).*FV(:,j,2);
      GMK(:,i,j) = GMK(:,i,j) + (FV(:,i,1).*tmp1 + FV(:,i,2).*tmp2) .* jac * GAUSSQUAD.RECT(k,3);
    end
  end
end

% 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 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);  SXY = zeros(NumEms,1);
% SU2 = zeros(NumEms,1);  SV2 = zeros(NumEms,1);  SUV = zeros(NumEms,1);
% SXU  = zeros(NumEms,1);  SXV  = zeros(NumEms,1);
% SYU  = zeros(NumEms,1);  SYV  = 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);
%   SXY = SXY + GAUSSQUAD.RECT(k,3) * jac .* (X.*Y);
%   SU2 = SU2 + GAUSSQUAD.RECT(k,3) * jac .* (U.^2);
%   SV2 = SV2 + GAUSSQUAD.RECT(k,3) * jac .* (V.^2);
%   SUV = SUV + GAUSSQUAD.RECT(k,3) * jac .* (U.*V);
%   SXU = SXU + GAUSSQUAD.RECT(k,3) * jac .* (X.*U);
%   SXV = SXV + GAUSSQUAD.RECT(k,3) * jac .* (X.*V);
%   SYU = SYU + GAUSSQUAD.RECT(k,3) * jac .* (Y.*U);
%   SYV = SYV + GAUSSQUAD.RECT(k,3) * jac .* (Y.*V);
% end

% %% Gram matrix for AC0 with coefficient matrix MatK 
% GMK = zeros(NumEms,4,4);
% GMK(:,1,1) = MatK(:,1,1).*S1;  GMK(:,1,2) = MatK(:,1,2).*S1;
% GMK(:,2,1) = MatK(:,2,1).*S1;  GMK(:,2,2) = MatK(:,2,2).*S1;
% GMK(:,1,3) = MatK(:,1,1).*SX + MatK(:,1,2).*SY;  GMK(:,3,1) = GMK(:,1,3);
% GMK(:,2,3) = MatK(:,2,1).*SX + MatK(:,2,2).*SY;  GMK(:,3,2) = GMK(:,2,3);
% GMK(:,3,3) = MatK(:,1,1).*SX2 + MatK(:,1,2).*SXY...
%            + MatK(:,2,1).*SXY + MatK(:,2,2).*SY2;
% GMK(:,1,4) = MatK(:,1,1).*SU + MatK(:,1,2).*SV;  GMK(:,4,1) = GMK(:,1,4);
% GMK(:,2,4) = MatK(:,2,1).*SU + MatK(:,2,2).*SV;  GMK(:,4,2) = GMK(:,2,4);
% GMK(:,3,4) = MatK(:,1,1).*SXU + MatK(:,1,2).*SXV...
%            + MatK(:,2,1).*SYU + MatK(:,2,2).*SYV;
% GMK(:,4,3) = GMK(:,3,4);
% GMK(:,4,4) = MatK(:,1,1).*SU2 + MatK(:,1,2).*SUV...
%            + MatK(:,2,1).*SUV + MatK(:,2,2).*SV2;

return;