function GMK1 = Hdiv_QuadriRT0_PiolaBas_GramMatK1(QuadriMesh,MatK1)
%% Hdiv elements on quadrilaterals: 
% RT[0] basis based on Piola transformation; Its Gram matrix with inv.mat. MatK1
% Using a Gaussian quadrature implicitly 
% James Liu, ColoState; 2012/07--2016/12 

%% Mesh info 
NumEms = QuadriMesh.NumEms; 
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);
       
%% Coefficients 
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);

%% Gram matrix 
GMK1 = zeros(NumEms,4,4);
GAUSSQUADRECT = SetGaussQuadRect(5,5);
NumQuadPts = size(GAUSSQUADRECT,1);
GAUSSQUADRECT1 = ones(NumQuadPts,2) - GAUSSQUADRECT(:,1:2);
MatJ = zeros(NumEms,2,2);
MatK1J = zeros(NumEms,2,2);
MatJTK1J = zeros(NumEms,2,2);
MatA = zeros(NumEms,2,2);
Phihat = zeros(NumEms,4,2);
WO = ones(NumEms,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;
  MatK1J(:,1,1) = MatK1(:,1,1).*MatJ(:,1,1) + MatK1(:,1,2).*MatJ(:,2,1);
  MatK1J(:,1,2) = MatK1(:,1,1).*MatJ(:,1,2) + MatK1(:,1,2).*MatJ(:,2,2);
  MatK1J(:,2,1) = MatK1(:,2,1).*MatJ(:,1,1) + MatK1(:,2,2).*MatJ(:,2,1);
  MatK1J(:,2,2) = MatK1(:,2,1).*MatJ(:,1,2) + MatK1(:,2,2).*MatJ(:,2,2);
  MatJTK1J(:,1,1) = MatJ(:,1,1).*MatK1J(:,1,1) + MatJ(:,2,1).*MatK1J(:,2,1);
  MatJTK1J(:,1,2) = MatJ(:,1,1).*MatK1J(:,1,2) + MatJ(:,2,1).*MatK1J(:,2,2);
  MatJTK1J(:,2,1) = MatJ(:,1,2).*MatK1J(:,1,1) + MatJ(:,2,2).*MatK1J(:,2,1);
  MatJTK1J(:,2,2) = MatJ(:,1,2).*MatK1J(:,1,2) + MatJ(:,2,2).*MatK1J(:,2,2);
  for i=1:2
    for j=1:2
      MatA(:,i,j) = MatJTK1J(:,i,j)./jac;
    end
  end
  Phihat(:,4,1) = (1-xhat)*WO;  Phihat(:,2,1) = xhat*WO;
  Phihat(:,1,2) = (1-yhat)*WO;  Phihat(:,3,2) = yhat*WO;
  for i=1:4
    for j=1:4
      GMK1(:,i,j) = GMK1(:,i,j) + GAUSSQUADRECT(k,3)*...
        ( MatA(:,1,1).*Phihat(:,i,1).*Phihat(:,j,1)... 
        + MatA(:,1,2).*Phihat(:,i,1).*Phihat(:,j,2)... 
        + MatA(:,2,1).*Phihat(:,i,2).*Phihat(:,j,1)... 
        + MatA(:,2,2).*Phihat(:,i,2).*Phihat(:,j,2));
    end
  end
end

return;