%% JL20170121: TO BE FINISHED/REVISED 
function [NumerPresEm,NumerVelCofRT0,NumerFlux,LMCR,FluxDscp] = ...
  Darcy_WG_QuadriQ0Q0RT0_PresVelFlux(QuadriMesh,PermK,sln)
%% Darcy: WG(Q0,Q0,RT[0])Quadri: 
% 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 = QuadriMesh.NumEms;
NumEgs = QuadriMesh.NumEgs;
% Assuming QuadriMesh.flag>=2
CofA = QuadriMesh.CofA;
CofB = QuadriMesh.CofB;

%% Auxiliary quantities 
% Elementwise 5 coeffs. in WG bas.fxns. 
CofWGB = zeros(NumEms,5);
CofWGB(:,1) = sln(1:NumEms);
for j=1:4
  CofWGB(:,1+j) = sln(NumEms+QuadriMesh.elem2edge(:,j));
end
% Elementwise coeffis. in RT[0] nmlz.bas. for disc.wk.grad. of 5 WG0 bas.fxns.
% CDWGB = zeros(NumEms,5,4) 
CDWGB = WG_QuadriQ0Q0RT0_CofRT0NmlzBas_DiscWkGradBasFxn(QuadriMesh);
% Computing auxiliary integrals 
S1 = zeros(NumEms,1);  SX2 = zeros(NumEms,1);
SX = zeros(NumEms,1);  SY2 = zeros(NumEms,1);
SY = zeros(NumEms,1);  SXY = zeros(NumEms,1);
GAUSSQUADRECT = SetGaussQuadRect(5,5);
NumQuadPts = size(GAUSSQUADRECT,1);
GAUSSQUADRECT1 = ones(NumQuadPts,2) - GAUSSQUADRECT(:,1:2);
qp = zeros(NumEms,2);
for k=1:NumQuadPts
  xhat = GAUSSQUADRECT1(k,1);
  yhat = GAUSSQUADRECT1(k,2);
  qp(:,1) = CofA(:,1) + CofA(:,2)*xhat + CofA(:,3)*yhat + CofA(:,4)*xhat*yhat;
  qp(:,2) = CofB(:,1) + CofB(:,2)*xhat + CofB(:,3)*yhat + CofB(:,4)*xhat*yhat;
  X = qp(:,1) - QuadriMesh.EmCntr(:,1);
  Y = qp(:,2) - QuadriMesh.EmCntr(:,2);
  jac = (CofB(:,3)+CofB(:,4)*xhat) .* (CofA(:,2)+CofA(:,4)*yhat)...
      - (CofA(:,3)+CofA(:,4)*xhat) .* (CofB(:,2)+CofB(:,4)*yhat);
  S1 = S1 + GAUSSQUADRECT(k,3)*jac;  % S1=area
  SX = SX + GAUSSQUADRECT(k,3)*jac.*X;
  SY = SY + GAUSSQUADRECT(k,3)*jac.*Y;
  SX2 = SX2 + GAUSSQUADRECT(k,3)*jac.*(X.*X);
  SY2 = SY2 + GAUSSQUADRECT(k,3)*jac.*(Y.*Y);
  SXY = SXY + GAUSSQUADRECT(k,3)*jac.*(X.*Y);
end
% Gram matrix for RT[0] 
GM = zeros(NumEms,4,4);
GM(:,1,1) = S1;  GM(:,2,2) = S1;  % S1=area
GM(:,1,3) = SX;  GM(:,2,4) = SY;
GM(:,3,1) = SX;  GM(:,4,2) = SY;
GM(:,3,3) = SX2;  GM(:,4,4) = SY2;
% Coeffis. in quadri.RT[0]nmlz.bas. for projection of permeability * bas.fxn.
CofProjKBasRT0 = zeros(NumEms,4,4);  % j-th fxn. projected onto k-th fxn. 
CofProjKBasRT0(:,1,1) = PermK(:,1,1);
CofProjKBasRT0(:,1,2) = PermK(:,2,1);
CofProjKBasRT0(:,2,1) = PermK(:,1,2);
CofProjKBasRT0(:,2,2) = PermK(:,2,2);
for ie=1:NumEms  % TO BE REVISED FOR EFFECIENCY! 
  EltGM = squeeze(GM(ie,:,:));
  % For the 3rd bas.fxn.
  EltRHS = [PermK(ie,1,1)*SX(ie);...
            PermK(ie,2,1)*SX(ie);...
            PermK(ie,1,1)*SX2(ie);...
            PermK(ie,2,1)*SXY(ie)];
  CofProjKBasRT0(ie,3,:) = EltGM\EltRHS;
  % For the 4th bas.fxn.
  EltRHS = [PermK(ie,1,2)*SY(ie);...
            PermK(ie,2,2)*SY(ie);...
            PermK(ie,1,2)*SXY(ie);...
            PermK(ie,2,2)*SY2(ie)];
  CofProjKBasRT0(ie,4,:) = EltGM\EltRHS;
end

% Elementwise 4 edge midpoints
MidPtEg = zeros(NumEms,4,2);
for j=1:3
  MidPtEg(:,j,:) = 0.5*(QuadriMesh.node(QuadriMesh.elem(:,j+1),:)...
                      + QuadriMesh.node(QuadriMesh.elem(:,j),:));
end
for j=4:4
  MidPtEg(:,4,:) = 0.5*(QuadriMesh.node(QuadriMesh.elem(:,1),:)...
                      + QuadriMesh.node(QuadriMesh.elem(:,4),:));
end
% Elementwise 4 edge outward (non-unit, with edge length) normal vectors 
NmlEg = zeros(NumEms,4,2);
for j=1:3
  tmp = QuadriMesh.node(QuadriMesh.elem(:,j+1),:)...
      - QuadriMesh.node(QuadriMesh.elem(:,j),:);
  NmlEg(:,j,:) = [tmp(:,2), -tmp(:,1)];
end
for j=4:4
  tmp = QuadriMesh.node(QuadriMesh.elem(:,1),:)...
      - QuadriMesh.node(QuadriMesh.elem(:,4),:);
  NmlEg(:,j,:) = [tmp(:,2), -tmp(:,1)];
end
% Elementwise 4 RT[0](quadrilateral) bas.fxns. 4 bulk normal fluxes 
BasFxnNmlFlux = zeros(NumEms,4,4);
for j=1:4  % j-th edge 
  BasFxnNmlFlux(:,1,j) = NmlEg(:,j,1);
  BasFxnNmlFlux(:,2,j) = NmlEg(:,j,2);
  BasFxnNmlFlux(:,3,j) = (MidPtEg(:,j,1)-QuadriMesh.EmCntr(:,1))...
                       .* NmlEg(:,j,1);
  BasFxnNmlFlux(:,4,j) = (MidPtEg(:,j,2)-QuadriMesh.EmCntr(:,2))...
                       .* NmlEg(:,j,2);
end

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

%% Computing the numerical velocity: Coefficients in RT[0] nmlz.bas. 
NumerVelCofRT0 = zeros(NumEms,4);
for ie=1:NumEms  % TO BE REVISED FOR EFFECIENCY 
  NumerVelCofRT0(ie,1:4) = ...
    - squeeze(CofWGB(ie,1:5)) ...
    * squeeze(CDWGB(ie,1:5,1:4)) ...
    * squeeze(CofProjKBasRT0(ie,1:4,1:4));
end

%% Computing elementwise 4 edge (total) outward normal fluxes 
NumerFlux = zeros(NumEms,4);
for j=1:4
  NumerFlux(:,j) = sum(NumerVelCofRT0.*squeeze(BasFxnNmlFlux(:,:,j)),2);
end

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

%% Computing (normal) flux discrepancies across all edges 
FluxDscp = zeros(NumEgs,1);
for j=1:4
  k = QuadriMesh.elem2edge(:,j);
  tmp = FluxDscp(k);
  FluxDscp(k) = tmp - NumerFlux(:,j);
end
% Correction for boundary edges 
k = find(QuadriMesh.BndryEdge>0);
FluxDscp(k) = 0;

return;