function EGM = CG_RectQ1_GramMat(RectMesh)
%% CG: Rect.mesh. Q1 elements: Elementwise Gram matrix (for all elements)
% James Liu, ColoState; 2012/07--2017/02 

%% Mesh info 
NumEms = RectMesh.NumEms;
% Element vertices are couterclockwise, starting from the lower-left corner
k1 = RectMesh.elem(:,1);  x1 = RectMesh.node(k1,1);  y1 = RectMesh.node(k1,2);
k2 = RectMesh.elem(:,3);  x2 = RectMesh.node(k2,1);  y2 = RectMesh.node(k2,2);

%% Auxiliary quantities 
Deltax = x2 - x1;  
Deltay = y2 - y1;
hxhy = (1/36) * Deltax .* Deltay;  

%% Setting up all entries 
EGM = zeros(NumEms,4,4);
% 
EGM(:,1,1) = 4*hxhy ;
EGM(:,1,4) = 2*hxhy;
EGM(:,1,2) = 2*hxhy;
EGM(:,1,3) =   hxhy;
%
EGM(:,4,1) = 2*hxhy;
EGM(:,4,4) = 4*hxhy;
EGM(:,4,2) =   hxhy;
EGM(:,4,3) = 2*hxhy; 
%
EGM(:,2,1) = 2*hxhy; 
EGM(:,2,4) =   hxhy; 
EGM(:,2,2) = 4*hxhy;
EGM(:,2,3) = 2*hxhy;
%
EGM(:,3,1) =   hxhy; 
EGM(:,3,4) = 2*hxhy; 
EGM(:,3,2) = 2*hxhy; 
EGM(:,3,3) = 4*hxhy;

% Another way 
% 
% EGM(:,1,1) = 4*hxhy ;
% EGM(:,1,2) = 2*hxhy;
% EGM(:,1,3) = 2*hxhy;
% EGM(:,1,4) =   hxhy;
% 
% EGM(:,2,1) = 2*hxhy;
% EGM(:,2,2) = 4*hxhy;
% EGM(:,2,3) =   hxhy;
% EGM(:,2,4) = 2*hxhy; 
% 
% EGM(:,3,1) = 2*hxhy; 
% EGM(:,3,2) =   hxhy; 
% EGM(:,3,3) = 4*hxhy;
% EGM(:,3,4) = 2*hxhy;
% 
% EGM(:,4,1) =   hxhy; 
% EGM(:,4,2) = 2*hxhy; 
% EGM(:,4,3) = 2*hxhy; 
% EGM(:,4,4) = 4*hxhy;

return;