function CofK = Poisson_SmplnCofK_QuadriMesh(fxnK,QuadriMesh,GAUSSQUAD)
%% Darcy on quadri.mesh: Sampling permeability using a Gaussian quadrature
% CofK: elementwise constant 2x2 SPD matrix (as averages on elements) 
% James Liu, Zhuoran Wang, ColoState; 2012/07--2021/03 

%% Mesh info 
NumEms = QuadriMesh.NumEms;
% Assuming QuadriMesh.flag>=2 
area = QuadriMesh.area;
CofA = QuadriMesh.CofA;
CofB = QuadriMesh.CofB;

%% NOTE: Coefficient K is a 2x2 SPD matrix (pointwise) 
CofK = zeros(NumEms,2,2);
NumQuadPts = size(GAUSSQUAD.RECT,1);
GAUSSQUADRECT1 = ones(NumQuadPts,2) - GAUSSQUAD.RECT(:,1:2);
for k=1:NumQuadPts
  xhat = GAUSSQUADRECT1(k,1);
  yhat = GAUSSQUADRECT1(k,2);
  qp = [CofA(:,1) + CofA(:,2)*xhat + CofA(:,3)*yhat + CofA(:,4)*xhat*yhat,...
        CofB(:,1) + CofB(:,2)*xhat + CofB(:,3)*yhat + CofB(:,4)*xhat*yhat];
  jac = (CofB(:,3)+CofB(:,4)*xhat) .* (CofA(:,2)+CofA(:,4)*yhat)...
      - (CofA(:,3)+CofA(:,4)*xhat) .* (CofB(:,2)+CofB(:,4)*yhat);
  MatK = fxnK(qp);
  for i=1:2
    for j=1:2
      CofK(:,i,j) = CofK(:,i,j) + GAUSSQUAD.RECT(k,3) * jac.*MatK(:,i,j);
    end
  end
end
for i=1:2
  for j=1:2
    CofK(:,i,j) = CofK(:,i,j)./area;
  end
end

return;