function [GMKQ,GMKT] = Hdiv_HybridRT0_NmlzBas_GramMatK(HybridMesh,MatK)
%% Gram matrix w/ K for RT0 in normalized basis for hybrid mesh of quadri. & tri.
% This is not an Hdiv-subspace (Hdiv named abused, actually broken "Hdiv")  
%   for quadrilaterals: [1;0], [0;1], [X;0], [0;Y]  
%   for triangles:      [1;0], [0;1], [X;Y] 
%   where X=x-xc, Y=y-yc and (xc,yc) is element center
% Using implicitly a Gaussian quadrature for quadrilaterals 
% James Liu, Zhuoran Wang, ColoState; 2012/07--2017/02 

%% Mesh info 
NumEmsQ = HybridMesh.NumEmsQ;
NumEmsT = HybridMesh.NumEmsT;

%% More  
MatKQ = MatK(HybridMesh.labelemQ,:,:);
MatKT = MatK(HybridMesh.labelemT,:,:);

%% Auxiliary quantities for quadrilateral elements
%
k1 = HybridMesh.elemQ(:,1);  k2 = HybridMesh.elemQ(:,2);  
k3 = HybridMesh.elemQ(:,3);  k4 = HybridMesh.elemQ(:,4);
%
x1 = HybridMesh.node(k1,1);  y1 = HybridMesh.node(k1,2);
x2 = HybridMesh.node(k2,1);  y2 = HybridMesh.node(k2,2);
x3 = HybridMesh.node(k3,1);  y3 = HybridMesh.node(k3,2);
x4 = HybridMesh.node(k4,1);  y4 = HybridMesh.node(k4,2);
%
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 for quadrilaterals
CofA = zeros(NumEmsQ,4);      CofB = zeros(NumEmsQ,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 integrals using Gaussian quadrature on [0,1]^2
SX = zeros(NumEmsQ,1);  SY = zeros(NumEmsQ,1);
SX2 = zeros(NumEmsQ,1);  SY2 = zeros(NumEmsQ,1);  SXY = zeros(NumEmsQ,1);
GAUSSQUADRECT = SetGaussQuadRect(5,5);
NumQuadPts = size(GAUSSQUADRECT,1);
GAUSSQUADRECT1 = ones(NumQuadPts,2) - GAUSSQUADRECT(:,1:2);
qp = zeros(NumEmsQ,2);
for k=1:NumQuadPts
  xhat = GAUSSQUADRECT1(k,1);
  yhat = GAUSSQUADRECT1(k,2);
  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;
  jac = (CofB(:,3)+CofB(:,4)*xhat) .* (CofA(:,2)+CofA(:,4)*yhat)...
      - (CofA(:,3)+CofA(:,4)*xhat) .* (CofB(:,2)+CofB(:,4)*yhat);
  X = qp(:,1) - xc;
  Y = qp(:,2) - yc;
  SX = SX + GAUSSQUADRECT(k,3)*jac.*X;
  SY = SY + GAUSSQUADRECT(k,3)*jac.*Y;
  SX2 = SX2 + GAUSSQUADRECT(k,3)*jac.*(X.*X);
  SY2 = SY2 + GAUSSQUADRECT(k,3)*jac.*(Y.*Y);
  SXY = SXY + GAUSSQUADRECT(k,3)*jac.*(X.*Y);
end

%% Gram matrix w/ K for all quadrilateral elements in the mixed mesh 
GMKQ = zeros(NumEmsQ,4,4);
for i=1:2
  for j=1:2
    GMKQ(:,i,j) = area.*MatKQ(:,i,j);
  end
end
GMKQ(:,1,3) = SX.*MatKQ(:,1,1);  GMKQ(:,1,4) = SY.*MatKQ(:,1,2);
GMKQ(:,2,3) = SX.*MatKQ(:,2,1);  GMKQ(:,2,4) = SY.*MatKQ(:,2,2);
% 
GMKQ(:,3,1) = SX.*MatKQ(:,1,1);  GMKQ(:,3,2) = SX.*MatKQ(:,1,2);
GMKQ(:,4,1) = SY.*MatKQ(:,2,1);  GMKQ(:,4,2) = SY.*MatKQ(:,2,2);
% 
GMKQ(:,3,3) = SX2.*MatKQ(:,1,1);  GMKQ(:,3,4) = SXY.*MatKQ(:,1,2);
GMKQ(:,4,3) = SXY.*MatKQ(:,2,1);  GMKQ(:,4,4) = SY2.*MatKQ(:,2,2);

%% Auxiliary quantities for triangular elements
% 
k1 = HybridMesh.elemT(:,1);  
k2 = HybridMesh.elemT(:,2);  
k3 = HybridMesh.elemT(:,3);
% 
x1 = HybridMesh.node(k1,1);  y1 = HybridMesh.node(k1,2); 
x2 = HybridMesh.node(k2,1);  y2 = HybridMesh.node(k2,2);
x3 = HybridMesh.node(k3,1);  y3 = HybridMesh.node(k3,2);
% 
area = 0.5*abs((x2-x1).*(y3-y1)-(x3-x1).*(y2-y1));
% 
SX2 = (area/36).*( (x1-x2).^2 + (x2-x3).^2 + (x3-x1).^2 );
SY2 = (area/36).*( (y1-y2).^2 + (y2-y3).^2 + (y3-y1).^2 );
SXY = (area/36).*((x1-x2).*(y1-y2) + (x2-x3).*(y2-y3) + (x3-x1).*(y3-y1));

%% Gram matrix w/ K for all triangular elements in the mixed mesh 
% For each element, Gram matrix is a 3x3 mat 
GMKT = zeros(NumEmsT,3,3);
GMKT(:,1,1) = MatKT(:,1,1) .* area;  
GMKT(:,1,2) = MatKT(:,2,1) .* area;
GMKT(:,2,1) = MatKT(:,1,2) .* area;
GMKT(:,2,2) = MatKT(:,2,2) .* area;
GMKT(:,3,3) = squeeze(MatKT(:,1,1)).*SX2 + squeeze(MatKT(:,1,2)).*SXY...
            + squeeze(MatKT(:,2,1)).*SXY + squeeze(MatKT(:,2,2)).*SY2;

return;