%% JL20170206: Old code for 2-phase problems, TO BE REVISED !!
function [NumerPres,NumerVel,NumerFlux] = ...
  Darcy_WG_RectQ0Q0RT0s_PresVelFlux(BndDescMat,RectMesh,PermKs,sln)
%% Darcy: WG(Q0,Q0;RT[0]): Computing numerical pressure, velocity, flux, etc. 
% Assuming PermKs is an elementwise constant scalar 
% James Liu, ColoState; 2012/07--2017/02

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

%% "Computing" the numerical pressure 
NumerPres.em = sln(1:NumEms);
NumerPres.eg = sln(NumEms+(1:NumEgs)');

%% Comput. elt.coeff. of disc.wk.grad. of numer.sln. in RT[0] nmlz.bas. 
CofDWGsln = WG_RectQ0Q0RT0_CofRT0NmlzBas_DiscWkGradSln(RectMesh,sln);

%% Computing the numerical velocity (a structure) 
% 
% .EmEgMP: At elementwise edge midpoint 
NumerVel.EmEgMP = zeros(NumEms,3,2);
for j=1:4
  X = RectMesh.EgMidPt(RectMesh.elem2edge(:,j),1) - RectMesh.EmCntr(:,1);
  Y = RectMesh.EgMidPt(RectMesh.elem2edge(:,j),2) - RectMesh.EmCntr(:,2);
  DWG_EgMP = [CofDWGsln(:,1), CofDWGsln(:,2)] ...
    + [CofDWGsln(:,3).*X, CofDWGsln(:,4).*Y];
  NumerVel.EmEgMP(:,j,:) = -[PermKs,PermKs].*DWG_EgMP;
end
% 
% .EmCntr: At element centers, size NumEms-by-2 
NumerVel.EmCntr = -[PermKs,PermKs].*[CofDWGsln(:,1),CofDWGsln(:,2)];

%% Computing NumerFlux: the numerical flux (a structure) 
%
% .EmEgNml: Elementwise edge normal flux (edge length already in) 
NumerFlux.EmEgNml = zeros(NumEms,4);
for j=1:4
  X = RectMesh.EgMidPt(RectMesh.elem2edge(:,j),1) - RectMesh.EmCntr(:,1);
  Y = RectMesh.EgMidPt(RectMesh.elem2edge(:,j),2) - RectMesh.EmCntr(:,2);
  DWG_EgMP = [CofDWGsln(:,1), CofDWGsln(:,2)]...
    + [CofDWGsln(:,3).*X, CofDWGsln(:,4).*Y];
  NumeVelEgMP = -[PermKs,PermKs].*DWG_EgMP;
  NumerFlux.EmEgNml(:,j) = RectMesh.SignEmEg(:,j)...
    .*dot(NumeVelEgMP,RectMesh.NmlEg(RectMesh.elem2edge(:,j),:),2)...
    .*RectMesh.LenEg(RectMesh.elem2edge(:,j));  % edge length already in 
end
% 
% .BndryFlux: Boundary outfluxes 
NumBndryPcs = size(BndDescMat,1);
NumerFlux.BndryFlux = zeros(NumBndryPcs,1);
% JL20140629: TO BE REVISED FOR EFFICIENCY
for ig=1:NumEgs
  k = RectMesh.BndryEdge(ig);
  if k>0 
    ie1 = RectMesh.edge2elem(ig,1);
    lbl = RectMesh.WhichEdge(ig,1);
    NumerFlux.BndryFlux(k) = NumerFlux.BndryFlux(k)...
      + NumerFlux.EmEgNml(ie1,lbl);
  end
end

return;