function A = aloc_gradgrad_L1rect(k,x)

% A = aloc_gradgrad_L1rect(tau,x)
%             calculates a local contribution to the system matrix
%             that corresponds to the bilinear form:
%              a(v,u) = k * \int_{K^e} grad(v(x))*grad(u(x)) dx
%             L1rect stands for the linear Lagrange (i.e. nodal)
%             rectangular elements.
%
%             The nodal degrees of freedom are nodal values at the following
%             nodes:
%               x1 (reference is (0,0)), x2 (reference is (1,0)),
%               x3 (reference is (1,1)), x4 (reference is (0,1)).

x1 = x(:,1);
x2 = x(:,2);
x3 = x(:,3);
x4 = x(:,4);

R = [x2-x1 x4-x1];
invRT = inv(R)';

wGP = 0.25*ones(1,4);
c1 = (sqrt(3)-1)/(2*sqrt(3));
c2 = (sqrt(3)+1)/(2*sqrt(3));
xGP = [[c1;c1] [c1;c2] [c2;c1] [c2;c2]];

A = zeros(4,4);

for i=1:length(wGP)
  B = invRT * [-(1-xGP(2,i)),(1-xGP(2,i)),xGP(2,i),-xGP(2,i);-(1-xGP(1,i)),-xGP(1,i),xGP(1,i),(1-xGP(1,i))];
  A = A + wGP(i)*B'*B;
end

A = k*A*abs(det(R));
