function [NumerPresEmAvg,NumerPresEmEg,NumerVelEmAvg] = ...
  Darcy_WG_TriP0P0RT0_Prsn_mandel(TriMesh,PermK,sln,n)
%% Darcy: WG(P0,P0;RT0)Tri: 
% Computing (numer. pressure and velocity elementwise averages for presentation 
% Assuming PermK is an elementwise constant 2x2 SPD matrix 
% James Liu, ColoState; 2017/02--2018/01 

%% Mesh info 
NumEms = TriMesh.NumEms;
NumEgs = TriMesh.NumEgs;
nx = n;
ny = n;
% Assuming TriMesh.flag>=2

%% "Computing" numerical pressure elementwise averages 
NumerPresEmAvg = sln(1:NumEms);

% For mandel's problem, get edge solns.
% For trig. meshes, this only stores vertical edges, not the diagonal
% edges. Need to check whether they are needed.
NumerPresEmEg = sln(NumEms+1:ny:NumEms+(nx+1)*ny);

%% Auxiliary quantities 
% Elementwise 5 coeffs. in WG bas.fxns. 
cof = zeros(NumEms,4);
cof(:,1) = sln(1:NumEms);
for j=1:3
  cof(:,1+j) = sln(NumEms+TriMesh.elem2edge(:,j));
end
% Elementwise cofs. in RT0 nmlz.bas.: Disc.wk.grad. of WG0 bas.fxns.
% CDWGB = zeros(NumEms,5,4) 
CDWGB = WG_TriP0P0RT0_CofRT0NmlzBas_DiscWkGradBasFxn(TriMesh);
CDWGS = zeros(NumEms,3);
for j=1:3  % The j-th RT0 bas.fxn.
  CDWGS(:,j) = CDWGB(:,1,j) .* sln(1:NumEms) ...
             + CDWGB(:,2,j) .* sln(NumEms+TriMesh.elem2edge(:,1)) ...
             + CDWGB(:,3,j) .* sln(NumEms+TriMesh.elem2edge(:,2)) ...
             + CDWGB(:,4,j) .* sln(NumEms+TriMesh.elem2edge(:,3)) ;
end
% GM = zeros(NumEms,4,4);
% GMK = zeros(NumEms,4,4);
% GM = Hdiv_TriRT0_NmlzBas_GramMat(TriMesh);
% GMK = Hdiv_TriRT0_NmlzBas_GramMatK(TriMesh,PermK);
% Projection coeffs. in RT0 nmlz.bas. for permeability * bas.fxn.
% ProjCof = zeros(NumEms,3,3);  % j-th fxn. projected onto k-th fxn. 
% for ie=1:NumEms  % TO BE REVISED FOR EFFECIENCY! 
%   EltGM = squeeze(GM(ie,:,:));
%   for j=1:3
%     EltRHS = squeeze(GMK(ie,:,j));
%     ProjCof(ie,j,:) = EltGM\EltRHS';
%   end
% end
ProjCof = Hdiv_TriRT0_NmlzBas_ProjMatK(TriMesh,PermK);

%% Numerical velocity: Coeffs. in RT0 local bases 
NumerVelCofRT0 = zeros(NumEms,3);
for k=1:3
  NumerVelCofRT0(:,k) = -sum(CDWGS.*squeeze(ProjCof(:,:,k)),2);
end

%% JL20180119: TO BE REVISED FOR GENERALITY 
NumerVelEmAvg = NumerVelCofRT0(:,1:2);

return;