function EGGMK = CG_RectQ1_GradGradMatK(RectMesh,MatK)
%% CG: Q1 rect.elements: grad-grad matrix with mat. MatK (for all elements)
% Assuming MatK is an elementwise constant 2x2 matrix 
% 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;
A = (1/6)*MatK(:,1,1).*(Deltay./Deltax);
D = (1/6)*MatK(:,2,2).*(Deltax./Deltay);
B = 0.25*MatK(:,1,2);
C = 0.25*MatK(:,2,1);

%% Setting up all entries 
EGGMK = zeros(NumEms,4,4);
% 
EGGMK(:,1,1) =  2*A + ( B+C) + 2*D;
EGGMK(:,1,2) = -2*A + (-B+C) +   D;
EGGMK(:,1,3) =   -A + (-B-C) -   D;
EGGMK(:,1,4) =    A + ( B-C) - 2*D;
% 
EGGMK(:,2,1) = -2*A + ( B-C) +   D;
EGGMK(:,2,2) =  2*A + (-B-C) + 2*D;
EGGMK(:,2,3) =    A + (-B+C) - 2*D;
EGGMK(:,2,4) =   -A + ( B+C) -   D;
% 
EGGMK(:,3,1) =   -A + (-B-C) -   D;
EGGMK(:,3,2) =    A + ( B-C) - 2*D;
EGGMK(:,3,3) =  2*A + ( B+C) + 2*D;
EGGMK(:,3,4) = -2*A + (-B+C) +   D;
% 
EGGMK(:,4,1) =    A + (-B+C) - 2*D;
EGGMK(:,4,2) =   -A + ( B+C) -   D;
EGGMK(:,4,3) = -2*A + ( B-C) +   D;
EGGMK(:,4,4) =  2*A + (-B-C) + 2*D;

return;