%% JL20171222: TO BE FINISHED/REVISED 
function [NumerDivEmAvg, NumerDivEmQuad] = ...
  NonlinearPoroElas_GS_QuadriBR1_dilation(lambda,mu,QuadriMesh,sln,GAUSSQUAD)
%% Computing numer. dilation for sampling nonlinear permeability
% Bernardi-Raugel 1st order elements on quadrilaterals 
% James Liu, ColoState; 2012/07--2017/12 

%% Mesh info 
NumNds = QuadriMesh.NumNds;
NumEms = QuadriMesh.NumEms;
NumEgs = QuadriMesh.NumEgs;
area = QuadriMesh.area;
CofA = QuadriMesh.CofA;
CofB = QuadriMesh.CofB;


%% Auxiliary quantities 
% Elementwise cofficients in 12 BR1 basis functions for numerical solution 
cof = zeros(NumEms,12);
for j=1:4  % Treating vertices and edges in one loop 
  k = QuadriMesh.elem(:,j);
  cof(:,2*j-1) = sln(2*k-1);
  cof(:,2*j  ) = sln(2*k  );
  k = QuadriMesh.elem2edge(:,j);
  cof(:,8+j) = sln(2*NumNds+k) .* QuadriMesh.EmEgSign(:,j);
end


%% Computing numerical dilation (div of displacement) elementwise averages 
NumerDivEmAvg = zeros(NumEms,1);
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);
  jac = (CofB(:,3)+CofB(:,4)*xhat) .* (CofA(:,2)+CofA(:,4)*yhat)...
      - (CofA(:,3)+CofA(:,4)*xhat) .* (CofB(:,2)+CofB(:,4)*yhat);
  [BFGrad,BFDiv] = CG_QuadriBR1_BasFxnGradDiv(QuadriMesh,xhat,yhat);
  NumerDivEmAvg(:,1) = NumerDivEmAvg(:,1) + GAUSSQUAD.RECT(k,3)...
    * jac .* sum(cof.*BFDiv,2) ./area;
end

%% Computing numerical dilation (div of displacement) at quadrature points
NumQuadPts = size(GAUSSQUAD.RECT,1);
NumerDivEmQuad = zeros(NumEms,NumQuadPts);
GAUSSQUADRECT1 = ones(NumQuadPts,2) - GAUSSQUAD.RECT(:,1:2);
for k=1:NumQuadPts
  xhat = GAUSSQUADRECT1(k,1);
  yhat = GAUSSQUADRECT1(k,2);
%   jac = (CofB(:,3)+CofB(:,4)*xhat) .* (CofA(:,2)+CofA(:,4)*yhat)...
%       - (CofA(:,3)+CofA(:,4)*xhat) .* (CofB(:,2)+CofB(:,4)*yhat);
  [BFGrad,BFDiv] = CG_QuadriBR1_BasFxnGradDiv(QuadriMesh,xhat,yhat);
  NumerDivEmQuad(:,k) = sum(cof.*BFDiv,2) ;
end

return;