function PermK = Darcy_SmplnPerm_PlygnMesh(fxnK, PlygnMesh, GAUSSQUAD)
%% Darcy on a plygn.mesh: Sampling permeability using a Gaussian quadrature
% PermK is an elementwise constant 2x2 SPD matrix (as averages) 
% James Liu, Zhuoran Wang, ColoState; 2012/07--2017/08 

%% Mesh info 
NumEms = PlygnMesh.NumEms;

%% Integration 
PermK = zeros(NumEms,2,2);
NumQuadPts = size(GAUSSQUAD.TRIG,1); 
for ie=1:NumEms
  n = PlygnMesh.ElemType(ie);
  for m=1:(n-2)
    x1 = PlygnMesh.node(PlygnMesh.elem{ie}(1),1);
    y1 = PlygnMesh.node(PlygnMesh.elem{ie}(1),2);
    x2 = PlygnMesh.node(PlygnMesh.elem{ie}(m+1),1);
    y2 = PlygnMesh.node(PlygnMesh.elem{ie}(m+1),2);
    x3 = PlygnMesh.node(PlygnMesh.elem{ie}(m+2),1);
    y3 = PlygnMesh.node(PlygnMesh.elem{ie}(m+2),2);
    ar = abs(0.5*((x2-x1).*(y3-y1)-(x3-x1).*(y2-y1)));
    for k=1:NumQuadPts
      qp = GAUSSQUAD.TRIG(k,1) * [x1,y1] ...
         + GAUSSQUAD.TRIG(k,2) * [x2,y2] ...
         + GAUSSQUAD.TRIG(k,3) * [x3,y3] ;
      PermK(ie,:,:) = PermK(ie,:,:) + GAUSSQUAD.TRIG(k,4) * fxnK(qp) * ar;
    end
  end
end

%% Computing averages 
for i=1:2
  for j=1:2
    PermK(:,i,j) = PermK(:,i,j) ./ PlygnMesh.area;
  end
end

return;