function [NumerPresCofP0,NumerVelCofCW0,NumerFlux,FluxDscp,LMCR] = ...
  Darcy_WG_PlygnP0P0CW0_PresVelFlux(PlygnMesh,PermK,sln,GAUSSQUAD)
%% Darcy: WG(P0,P0;CW0) on a polygonal mesh: 
% Assuming PermK is an elementwise constant 2x2 SPD matrix 
% James Liu, Zhuoran Wang, ColoState; 2012/07--2018/12 

%% Mesh info 
NumEms = PlygnMesh.NumEms;
NumEgs = PlygnMesh.NumEgs;

%% Numerical pressure: Coeff. in P0 for the interior part  
NumerPresCofP0 = sln(1:NumEms);

%% Auxiliary 
cof = cell(NumEms,1);
for ie=1:NumEms
  n = PlygnMesh.ElemType(ie);
  cof{ie} = zeros(1,n+1);
  cof{ie}(1) = sln(ie);
  for j=1:n
    cof{ie}(1+j) = sln(NumEms+PlygnMesh.elem2edge{ie}(j));
  end
end
% CDWGB cell NumEms*1; Each entry (n+1)*n matrix; n=ElemType varies 
CDWGB = WG_PlygnP0P0CW0_Wachspress_DiscWkGradBasFxn(PlygnMesh,GAUSSQUAD);
% 
CDWGS = cell(NumEms,1);
for ie=1:NumEms
  n = PlygnMesh.ElemType(ie);
  CDWGS{ie} = zeros(1,n);
  for j=1:n
    CDWGS{ie}(j) = cof{ie} * CDWGB{ie}(:,j);
  end
end
% 
GM = Hdiv_PlygnCW0_Wachspress_GramMat(PlygnMesh,GAUSSQUAD);
GMK = Hdiv_PlygnCW0_Wachspress_GramMatK(PlygnMesh,PermK,GAUSSQUAD);
% 
% JL20170808: TO BE REVISED FOR EFFICIENCY 
ProjCof = cell(NumEms,1);
for ie=1:NumEms
  n = PlygnMesh.ElemType(ie);
  ProjCof{ie} = zeros(n,n);
  EltGM = GM{ie};
  for j=1:n
    EltRHS = GMK{ie}(:,j);
    ProjCof{ie}(j,:) = EltGM\EltRHS;
  end
end

%% Numerical velocity: Coeffs. in CW0 elementwise 
NumerVelCofCW0 = cell(NumEms,1);
for ie=1:NumEms
  n = PlygnMesh.ElemType(ie);
  NumerVelCofCW0{ie} = zeros(1,n);
  for j=1:n 
    NumerVelCofCW0{ie}(j) = - CDWGS{ie} * ProjCof{ie}(:,j);  
  end
end

%% JL20170808: TO BE FINISHED BY ... 
%% Numerical (normal) flux: Elementwise n edges 
NumerFlux = cell(NumEms,1);
for ie=1:NumEms
  n = PlygnMesh.ElemType(ie);
  NumerFlux{ie} = zeros(1,n);
end

%% JL20170808:  THEORETICALLY 0 
%% Flux (normal) discrepancy cross all edges 
FluxDscp = zeros(NumEgs,1);

%% JL20170808: THEORETICALLY 0 
%% Elementwise Local-Mass-Conservation-Residual (LMCR) 
LMCR = zeros(NumEms,1);

return;