% JL20210116: BEING REVISED By James 
function PermK = Darcy_SmplnNlinPerm_QuadriMesh( ... 
  EqnBIC,QuadriMesh,GAUSSQUAD,NumerDspl)
%% Darcy on a quadri.mesh: Sampling nonlinear permeability using a Gaussian quadrature
% PermK: An elementwise constant 2x2 SPD matrix (as averages on elements) 
% Zhuoran Wang, SYSU; James Liu, ColoState; 2017/07--2021/01 

%% Equation info 
lambda = EqnBIC.lambda;
    mu = EqnBIC.mu;
    c0 = EqnBIC.c0;  
  fxnK = EqnBIC.fxnK;
 alpha = EqnBIC.alpha;
parama = EqnBIC.parama;  % Used in "kappa = 1 + a*dilation" 

%% Mesh info 
NumEms = QuadriMesh.NumEms;
k1 = QuadriMesh.elem(:,1);  k2 = QuadriMesh.elem(:,2);  
k3 = QuadriMesh.elem(:,3);  k4 = QuadriMesh.elem(:,4);
x1 = QuadriMesh.node(k1,1);  y1 = QuadriMesh.node(k1,2);
x2 = QuadriMesh.node(k2,1);  y2 = QuadriMesh.node(k2,2);
x3 = QuadriMesh.node(k3,1);  y3 = QuadriMesh.node(k3,2);
x4 = QuadriMesh.node(k4,1);  y4 = QuadriMesh.node(k4,2);
area = 0.5*( (x2-x1).*(y3-y1) - (x3-x1).*(y2-y1)...
           + (x3-x1).*(y4-y1) - (x4-x1).*(y3-y1));

%% Coefficients of bilinear mapping 
CofA = zeros(NumEms,4);
CofB = zeros(NumEms,4);
CofA(:,1) = x1;
CofA(:,2) = x2-x1;
CofA(:,3) = x4-x1;
CofA(:,4) = (x3+x1)-(x2+x4);
CofB(:,1) = y1;
CofB(:,2) = y2-y1;
CofB(:,3) = y4-y1;
CofB(:,4) = (y3+y1)-(y2+y4);

%% NOTE: Permability is a 2x2 SPD matrix on each element 
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
  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);
%   MatK = fxnKt(qp,t);
  [NumerDivEmAvg, NumerDivEmQuad] = ...
   NonlinearPoroElas_GS_QuadriBR1_dilation(lambda,mu,QuadriMesh,NumerDspl,GAUSSQUAD);
  % for the case with kappa = 1+parameter*dilation
  % this will be changed for other cases.
  MatK_component = 1 + parama * NumerDivEmAvg;
  MatK(:,1,1) = MatK_component;
  MatK(:,2,2) = MatK_component;
  for i=1:2
    for j=1:2
      PermK(:,i,j) = PermK(:,i,j) + GAUSSQUAD.RECT(k,3)*(jac.*MatK(:,i,j));
    end
  end
end
for i=1:2
  for j=1:2
    PermK(:,i,j) = PermK(:,i,j)./area;
  end
end

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

return;