function PermK = Darcy_SmplnPerm_RectMesh(fxnK,RectMesh,GAUSSQUAD)
%% Darcy a rectangular mesh: Sampling permeability using a Gaussian quadrature
% PermK will be an elementwise constant 2x2 SPD matrix (as averages on elements) 
% James Liu, ColoState; 2012/07--2017/02 

% Mesh info 
NumEms = RectMesh.NumEms;
x1 = RectMesh.node(RectMesh.elem(:,1),1);
y1 = RectMesh.node(RectMesh.elem(:,1),2);
x2 = RectMesh.node(RectMesh.elem(:,3),1);
y2 = RectMesh.node(RectMesh.elem(:,3),2);

%% 
PermK = zeros(NumEms,2,2);
NumQuadPts = size(GAUSSQUAD.RECT,1);
GAUSSQUADRECT1 = ones(NumQuadPts,2) - GAUSSQUAD.RECT(:,1:2);
qp = zeros(NumEms,2);
for k=1:NumQuadPts
  qp(:,1) = GAUSSQUAD.RECT(k,1)*x1 + GAUSSQUADRECT1(k,1)*x2;
  qp(:,2) = GAUSSQUAD.RECT(k,2)*y1 + GAUSSQUADRECT1(k,2)*y2;
  PermK = PermK + GAUSSQUAD.RECT(k,3) * fxnK(qp);
end

return;