function PermK = Darcy_SmplnPerm_HybridMesh(fxnK,HybridMesh,GAUSSQUAD)
%% Darcy a hybrid mesh (quadrilaterals & triangles): Sampling permeability 
% PermK will be an elementwise constant 2x2 SPD matrix (as averages on elements) 
% James Liu, Zhuoran Wang, ColoState; 2012/07--2017/02 

%% Mesh info 
NumEmsQ = HybridMesh.NumEmsQ;
NumEmsT = HybridMesh.NumEmsT;
LblQ = HybridMesh.labelemQ;
LblT = HybridMesh.labelemT;

%% Elementwise constant 2x2 matrix (averages) 
PermK = zeros(HybridMesh.NumEms,2,2); 

%% Auxiliary quantities for quadrilateral elements
k1 = HybridMesh.elemQ(:,1);  k2 = HybridMesh.elemQ(:,2);  
k3 = HybridMesh.elemQ(:,3);  k4 = HybridMesh.elemQ(:,4);
x1 = HybridMesh.node(k1,1);  y1 = HybridMesh.node(k1,2);
x2 = HybridMesh.node(k2,1);  y2 = HybridMesh.node(k2,2);
x3 = HybridMesh.node(k3,1);  y3 = HybridMesh.node(k3,2);
x4 = HybridMesh.node(k4,1);  y4 = HybridMesh.node(k4,2);
area = 0.5*( (x2-x1).*(y3-y1) - (x3-x1).*(y2-y1)...
           + (x3-x1).*(y4-y1) - (x4-x1).*(y3-y1));
% Coefficients for the bilinear mapping for quadrilaterals
CofA = zeros(NumEmsQ,4);      CofB = zeros(NumEmsQ,4);
CofA(:,1) = x1;               CofB(:,1) = y1;
CofA(:,2) = x2-x1;            CofB(:,2) = y2-y1;
CofA(:,3) = x4-x1;            CofB(:,3) = y4-y1;
CofA(:,4) = (x1+x3)-(x2+x4);  CofB(:,4) = (y1+y3)-(y2+y4);

%% JL20160913: TO BE REVISED FOR EFFECIENCY 
PermK(LblQ,:,:) = zeros(NumEmsQ,2,2);
NumQuadPts = size(GAUSSQUAD.RECT,1);
GAUSSQUADRECT1 = ones(NumQuadPts,2) - GAUSSQUAD.RECT(:,1:2);
qp = zeros(NumEmsQ,2);
JAC = zeros(NumEmsQ,2,2);
for k=1:NumQuadPts
  xhat = GAUSSQUADRECT1(k,1);
  yhat = GAUSSQUADRECT1(k,2);
  qp(:,1) = CofA(:,1) + CofA(:,2)*xhat + CofA(:,3)*yhat + CofA(:,4)*xhat*yhat;
  qp(:,2) = 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);
  for i=1:2
    for j=1:2
      JAC(:,i,j) = jac;
    end
  end
  PermK(LblQ,:,:) = PermK(LblQ,:,:) + GAUSSQUAD.RECT(k,3)*(JAC.*fxnK(qp));
end
for i=1:2
  for j=1:2
    PermK(LblQ,i,j) = PermK(LblQ,i,j)./area;
  end
end

%% JL20160517: ALT.: Splitting ea.quadri. into 2 triangles for integration 

%% For triangular elements, MUCH SIMPLER ! 
PermK(LblT,:,:) = zeros(NumEmsT,2,2);
NumQuadPts = size(GAUSSQUAD.TRIG,1);
for k=1:NumQuadPts
  qp = GAUSSQUAD.TRIG(k,1) * HybridMesh.node(HybridMesh.elemT(:,1),:) ...
     + GAUSSQUAD.TRIG(k,2) * HybridMesh.node(HybridMesh.elemT(:,2),:) ...
     + GAUSSQUAD.TRIG(k,3) * HybridMesh.node(HybridMesh.elemT(:,3),:) ;
  PermK(LblT,:,:) = PermK(LblT,:,:) + GAUSSQUAD.TRIG(k,4)*fxnK(qp);
end

return;