function GMK = Hdiv_RectRT0_NmlzBas_GramMatK(RectMesh,MatK)
%% Hdiv local RT[0] on a rectangular mesh: Gram matrix for nmlz.bas. with MatK 
% MatK is an elementwise constant 2x2 matrix 
% This is not an Hdiv-subspace (Hdiv named abused, actually broken "Hdiv") 
%   Normalized basis: [1;0], [0;1], [X;0], [0;Y]
%   where X=x-xc, Y=y-yc and (xc,yc) is element center
% James Liu, ColoState; 2012/07--2017/02 

%% Mesh info and auxiliary quantties 
NumEms = RectMesh.NumEms;
k1 = RectMesh.elem(:,1);  k2 = RectMesh.elem(:,3);
x1 = RectMesh.node(k1,1);  y1 = RectMesh.node(k1,2); 
x2 = RectMesh.node(k2,1);  y2 = RectMesh.node(k2,2);
Deltax = x2 - x1;  
Deltay = y2 - y1;
area = Deltax.*Deltay;

%% Gram matrices for all elements in the rectangular mesh 
% For each element, Gram matrix is a 4x4 matrix 
GMK = zeros(NumEms,4,4);
GMK(:,1,1) = MatK(:,1,1) .* area;
GMK(:,1,2) = MatK(:,1,2) .* area;
GMK(:,2,1) = MatK(:,2,1) .* area;
GMK(:,2,2) = MatK(:,2,2) .* area;
GMK(:,3,3) = MatK(:,1,1) .* area .* (Deltax.^2)/12;
GMK(:,4,4) = MatK(:,2,2) .* area .* (Deltay.^2)/12;

return;