%% JL20170122: TO BE FINISHED 
function [NumerPresEm, NumerVelCofRT0, NumerFlux, LMCR, FluxDscp] = ...
  Darcy_WG_HybridPQ0RT0_PresVelFlux(HybridMesh,PermK,sln,GAUSSQUAD)
%% Darcy:  
% Computing (numerical) pressure, velocity, flux, etc. 
% Assuming PermK is an elementwise constant 2x2 SPD matrix 
% James Liu, Zhuoran Wang, ColoState; 2012/07--2017/02 

%% Mesh info 
NumEms = HybridMesh.NumEms;
NumEmsQ = HybridMesh.NumEmsQ;
NumEmsT = HybridMesh.NumEmsT;
NumEgs = HybridMesh.NumEgs;
LblQ = HybridMesh.labelemQ;
LblT = HybridMesh.labelemT;

%% Auxiliary quantities 

%% "Computing" the numerical pressure 
NumerPresEm = sln(1:NumEms);

%% Computing the numerical velocity (elementwise 4/3 coeffs. in RT[0]/RT0) 
[CDWGBQ,CDWGBT] = WG_HybridPQ0_CofRT0NmlzBas_DiscWkGradBasFxn(HybridMesh);
[ProjCofQ,ProjCofT] = Hdiv_HybridRT0_NmlzBas_ProjMatK(HybridMesh,PermK,GAUSSQUAD);
NumerVelCofRT0.quadri = zeros(NumEmsQ,4);
NumerVelCofRT0.tri    = zeros(NumEmsT,3);
% For quadrilaterals 
CofWGB = zeros(NumEmsQ,5);  % Coefficients in the WG basis functions 
CofWGB(:,1) = sln(LblQ);
for j=1:4
  CofWGB(:,1+j) = sln(NumEms+HybridMesh.elemQedge(:,j));
end
for k=1:4
  for j=1:4
    tmp = zeros(NumEmsQ,1);
    for i=1:5
      tmp = tmp + CofWGB(:,i).*squeeze(CDWGBQ(:,i,j));
    end
    NumerVelCofRT0.quadri(:,k) = NumerVelCofRT0.quadri(:,k)...
      - tmp.*squeeze(ProjCofQ(:,j,k));
  end
end
% For triangles 
CofWGB = zeros(NumEmsT,4);  % Coefficients in the WG basis functions 
CofWGB(:,1) = sln(LblT);
for j=1:3
  CofWGB(:,1+j) = sln(NumEms+HybridMesh.elemTedge(:,j));
end
for k=1:3
  for j=1:3
    tmp = zeros(NumEmsT,1);
    for i=1:4
      tmp = tmp + CofWGB(:,i).*squeeze(CDWGBT(:,i,j));
    end
    NumerVelCofRT0.tri(:,k) = NumerVelCofRT0.tri(:,k)...
      - tmp.*squeeze(ProjCofT(:,j,k));
  end
end

%% Computing the numerical fluxes (elementwise on 4 or 3 edges) 
NumerFlux.quadri = zeros(NumEmsQ,4);
NumerFlux.tri    = zeros(NumEmsT,3);

%% "Computing" local-mass-conservation residuals for all elements 
LMCR = zeros(NumEms,1);
% For quadrilaterals 
% For triangles 

%% "Computing" (normal) flux discrepancies across all edges 
FluxDscp = zeros(NumEgs,1);

return;