function [NumerVelEmAvg,NumerPresEmAvg] = ... 
  StokesDarcy_QuadriCGBR1WGAC0_1meshBT_Presn( ...
    EqnBC, QuadriMesh, sln, CofQhg, GAUSSQUAD) 
% Assuming 
% -- QuadriMesh is 1 mesh for the combined domains; 
% -- sln contains all lin.comb.coeffs. for bas.fxns. Stokes & Darcy 
% Graham Harper, James Liu, ColoState; 2014/09--2020/09 

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

%% Auxiliary mesh info 
LblEmDarcy  = QuadriMesh.DarcyEms;
NumNdsStokes = nnz(QuadriMesh.posStokesNd);
NumEgsStokes = nnz(QuadriMesh.posStokesEg);
NumEmsStokes = length(QuadriMesh.StokesEms);
%NumStokesDoFs = 2*NumNdsStokes + NumStokesEgs + NumStokesEms;
%NumDarcyNds = nnz(QuadriMesh.posDarcyNd);
NumEgsDarcy = nnz(QuadriMesh.posDarcyEg);
NumEmsDarcy = length(QuadriMesh.DarcyEms);
NumDOFsDarcy = NumEmsDarcy + NumEgsDarcy;

%% Elementwise cofficients in 12 BR1 basis functions for numerical solution 
cof = zeros(length(QuadriMesh.StokesEms),12);
for j=1:4  % Treating vertices and edges in one loop 
  k = QuadriMesh.posStokesNd(QuadriMesh.elem(QuadriMesh.StokesEms,j));
  cof(:,2*j-1) = sln(NumEmsDarcy+NumEgsDarcy+2*k-1);
  cof(:,2*j  ) = sln(NumEmsDarcy+NumEgsDarcy+2*k  );
  k = QuadriMesh.posStokesEg(QuadriMesh.elem2edge(QuadriMesh.StokesEms,j));
  cof(:,8+j) = sln(NumEmsDarcy+NumEgsDarcy+2*NumNdsStokes+k) .* QuadriMesh.EmEgSign(QuadriMesh.StokesEms,j);
end

%% Pressure first
DarcyPresInt = sln(1:length(QuadriMesh.DarcyEms));
StokesPresInt = sln((end-length(QuadriMesh.StokesEms)+1):end);

% % normalize for debugging
% DarcyPresInt = DarcyPresInt - min(DarcyPresInt);
% StokesPresInt = StokesPresInt - min(StokesPresInt);
% DarcyPresInt = DarcyPresInt / abs(max(DarcyPresInt));
% StokesPresInt = StokesPresInt / abs(max(StokesPresInt)) * sin(1);

NumerPresEmAvg = zeros(QuadriMesh.NumEms,1);
NumerPresEmAvg(QuadriMesh.DarcyEms) = DarcyPresInt;
NumerPresEmAvg(QuadriMesh.StokesEms) = StokesPresInt;

%% Numerical velocity elementwise averages 
NumerVelEmAvg = zeros(NumEms,2);

%% Stokes: Numerical velocity elementwise averages 
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(QuadriMesh.StokesEms,3)+CofB(QuadriMesh.StokesEms,4)*xhat) .* (CofA(QuadriMesh.StokesEms,2)+CofA(QuadriMesh.StokesEms,4)*yhat)...
      - (CofA(QuadriMesh.StokesEms,3)+CofA(QuadriMesh.StokesEms,4)*xhat) .* (CofB(QuadriMesh.StokesEms,2)+CofB(QuadriMesh.StokesEms,4)*yhat);
  BFVal = CG_QuadriBR1_BasFxnVal(QuadriMesh,xhat,yhat);
  NumerVelEmAvg(QuadriMesh.StokesEms,1) = NumerVelEmAvg(QuadriMesh.StokesEms,1) + GAUSSQUAD.RECT(k,3) ...
    * jac ./area(QuadriMesh.StokesEms) .* sum(cof.*squeeze(BFVal(QuadriMesh.StokesEms,:,1)),2);
  NumerVelEmAvg(QuadriMesh.StokesEms,2) = NumerVelEmAvg(QuadriMesh.StokesEms,2) + GAUSSQUAD.RECT(k,3) ...
    * jac ./area(QuadriMesh.StokesEms) .* sum(cof.*squeeze(BFVal(QuadriMesh.StokesEms,:,2)),2);
end

%% Darcy: Numerical velocity elementwise averages 
% Elementwise 5 coeffs. in WG bas.fxns. 
cof = zeros(NumEms,5);
cof(QuadriMesh.DarcyEms,1) = sln(1:NumEmsDarcy);
for j=1:4
  k = QuadriMesh.posDarcyEg(QuadriMesh.elem2edge(QuadriMesh.DarcyEms,j));
  cof(QuadriMesh.DarcyEms,1+j) = sln(NumEmsDarcy+k);
end
% Elementwise cofs. in AC0 nmlz.Piola.bas.: Disc.wk.grad. of WG0 bas.fxns.
% CDWGB = zeros(NumEms,5,4) 
CDWGB = WG_QuadriP0P0AC0_CofAC0NmlzPiolaBas_DiscWkGradBasFxn(... 
  QuadriMesh, GAUSSQUAD);
CDWGS = zeros(NumEms,4);
for j=1:4
  CDWGS(:,j) = sum(cof.*squeeze(CDWGB(:,:,j)),2);
end
PermK = Darcy_SmplnPerm_QuadriMesh(EqnBC.fxnK_Darcy, QuadriMesh, GAUSSQUAD);
% GM = zeros(NumEms,4,4);
% GMK = zeros(NumEms,4,4);
GM = Hdiv_QuadriAC0_NmlzPiolaBas_GramMat(QuadriMesh, GAUSSQUAD);
GMK = Hdiv_QuadriAC0_NmlzPiolaBas_GramMatK(QuadriMesh, PermK, GAUSSQUAD);
% Projection coeffs. in AC0 nmlz.Piola.bas. for permeability * bas.fxn.
ProjCof = zeros(NumEms,4,4);  % j-th fxn. projected onto k-th fxn. 
for ie=1:NumEms  % TO BE REVISED FOR EFFECIENCY! 
  EltGM = squeeze(GM(ie,:,:));
  for j=1:4
    EltRHS = squeeze(GMK(ie,:,j));
    ProjCof(ie,j,:) = EltGM\EltRHS';
  end
end
% Numerical velocity: Coeffs. in AC0 elementwise 
NumerVelCofAC0 = zeros(NumEms,4);
for k=1:4
  NumerVelCofAC0(:,k) = -sum(CDWGS.*squeeze(ProjCof(:,:,k)),2);  % - CofQhg;
end

% Darcy: Numerical velocity average "at element centers" 
xhat = 0.5;  yhat = 0.5;
FV = Hdiv_QuadriAC0_NmlzPiolaBas_FxnVal(QuadriMesh,xhat,yhat);
NumerVelEmAvg(QuadriMesh.DarcyEms,1) ... 
  = NumerVelCofAC0(QuadriMesh.DarcyEms,1) ... 
  + NumerVelCofAC0(QuadriMesh.DarcyEms,3) .* FV(QuadriMesh.DarcyEms,3,1) ...
  + NumerVelCofAC0(QuadriMesh.DarcyEms,4) .* FV(QuadriMesh.DarcyEms,4,1);
NumerVelEmAvg(QuadriMesh.DarcyEms,2) ... 
  = NumerVelCofAC0(QuadriMesh.DarcyEms,2) ... 
  + NumerVelCofAC0(QuadriMesh.DarcyEms,3) .* FV(QuadriMesh.DarcyEms,3,2) ... 
  + NumerVelCofAC0(QuadriMesh.DarcyEms,4) .* FV(QuadriMesh.DarcyEms,4,2);
% NumerVelEmAvg(QuadriMesh.DarcyEms,:) = ... 
%   NumerVelCofAC0(QuadriMesh.DarcyEms,1:2) - CofQhg(:,1:2);

return;