function [NumerPresEmAvg, NumerVelEmCntr] = ...
  Darcy_WG_PlygnP0P0CW0_Presn(PlygnMesh,NumerPresCofP0,NumerVelCofCW0)
%% Darcy: WG(P0,P0;CW0)Plygn: Quantities for graphical presentation  
% James Liu, ColoState; 2012/07--2018/12 

%% Mesh info 
NumEms = PlygnMesh.NumEms;

%% Numerical pressure elementwise 
NumerPresEmAvg = NumerPresCofP0;

%% Auxiliary  
% ConvMat: Cell NumEms*1; Each entry (n+1)*n matrix; n varies 
ConvMat = Hdiv_PlygnCW0_Wachspress_ConvMatBasFrm(PlygnMesh);

%% Numerical velocity at element center 
NumerVelEmCntr = zeros(NumEms,2);
for ie=1:NumEms
  n = PlygnMesh.ElemType(ie);
  vrtx = PlygnMesh.node(PlygnMesh.elem{ie}(:),:);
  cntr = PlygnMesh.EmCntr(ie,:);
  [lambda,lambdagrad,lambdacurl] = Plygn_Wachspress_ValGradCurl(vrtx,cntr);
  frm = zeros(n+1,2);
  for i=1:n
    frm(1+i,1) = lambdacurl(i,1);
    frm(1+i,2) = lambdacurl(i,2);
  end
  bas = zeros(n,2);
  for i=1:n
    bas(i,1) = ConvMat{ie}(i,:) * frm(:,1);
    bas(i,2) = ConvMat{ie}(i,:) * frm(:,2);
  end
  NumerVelEmCntr(ie,1) = NumerVelCofCW0{ie} * bas(:,1);
  NumerVelEmCntr(ie,2) = NumerVelCofCW0{ie} * bas(:,2);
end

return;