function EGGM = CG_RectQ1_GradGradMat(RectMesh)
%% CG: Rect.mesh. Q1 elements: grad-grad 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/6) * Deltax ./ Deltay;  
hyhx = (1/6) * Deltay ./ Deltax;

%% Setting up all entries 
EGGM = zeros(NumEms,4,4);
% 
EGGM(:,1,1) =  2*hxhy + 2*hyhx;
EGGM(:,1,4) = -2*hxhy +   hyhx;
EGGM(:,1,2) =    hxhy - 2*hyhx;
EGGM(:,1,3) =   -hxhy -   hyhx;
% 
EGGM(:,4,1) = -2*hxhy +   hyhx;
EGGM(:,4,4) =  2*hxhy + 2*hyhx;
EGGM(:,4,2) =   -hxhy -   hyhx;
EGGM(:,4,3) =    hxhy - 2*hyhx; 
% 
EGGM(:,2,1) =    hxhy - 2*hyhx; 
EGGM(:,2,4) =   -hxhy -   hyhx; 
EGGM(:,2,2) =  2*hxhy + 2*hyhx;
EGGM(:,2,3) = -2*hxhy +   hyhx;
% 
EGGM(:,3,1) =   -hxhy -   hyhx; 
EGGM(:,3,4) =    hxhy - 2*hyhx; 
EGGM(:,3,2) = -2*hxhy +   hyhx; 
EGGM(:,3,3) =  2*hxhy + 2*hyhx;

return;