function [ProjPresEm,ProjVelCofRT0] = Darcy_WG_QuadriQ0Q0RT0_ProjPresVel(...
  EqnBC,QuadriMesh,GAUSSQUAD)
%% Darcy: WG(Q0,Q0,RT[0])Quadri: 
% Projecting pressure elementwise and edgewise 
% Projecting velocity elementwise 
% % and then computing the fluxes on 4 edges 
% James Liu, Zhuoran Wang, ColoState; 2012/07--2017/02 

%% Mesh info 
NumEms = QuadriMesh.NumEms;
% NumEgs = QuadriMesh.NumEgs;
% Assuming TriMesh.flag>=2 
k1 = QuadriMesh.elem(:,1);  k2 = QuadriMesh.elem(:,2);  
k3 = QuadriMesh.elem(:,3);  k4 = QuadriMesh.elem(:,4);
x1 = QuadriMesh.node(k1,1);  y1 = QuadriMesh.node(k1,2);
x2 = QuadriMesh.node(k2,1);  y2 = QuadriMesh.node(k2,2);
x3 = QuadriMesh.node(k3,1);  y3 = QuadriMesh.node(k3,2);
x4 = QuadriMesh.node(k4,1);  y4 = QuadriMesh.node(k4,2);
% Coefficients for the bilinear mapping 
CofA = zeros(NumEms,4);       CofB = zeros(NumEms,4);
CofA(:,1) = x1;               CofB(:,1) = y1;
CofA(:,2) = x2-x1;            CofB(:,2) = y2-y1;
CofA(:,3) = x4-x1;            CofB(:,3) = y4-y1;
CofA(:,4) = (x1+x3)-(x2+x4);  CofB(:,4) = (y1+y3)-(y2+y4);
x = zeros(NumEms,5);
y = zeros(NumEms,5);
x(:,1) = x1;  x(:,2) = x2;  x(:,3) = x3;  x(:,4) = x4;  x(:,5) = x1;
y(:,1) = y1;  y(:,2) = y2;  y(:,3) = y3;  y(:,4) = y4;  y(:,5) = y1;
xc = 0.25*sum(x(:,1:4),2);
yc = 0.25*sum(y(:,1:4),2);
area = QuadriMesh.area;
Xm = zeros(NumEms,4);
Ym = zeros(NumEms,4);
for j=1:4
  Xm(:,j) = 0.5*(x(:,j)+x(:,j+1)) - xc;
  Ym(:,j) = 0.5*(y(:,j)+y(:,j+1)) - yc;
end
% Elementwise edge normal vector with length info 
EmEgNml = zeros(NumEms,4,2);
for j=1:4
  EmEgNml(:,j,1) = y(:,j+1) - y(:,j);
  EmEgNml(:,j,2) = x(:,j) - x(:,j+1);
end

%% Projecting pressure elementwise 
ProjPresEm = zeros(NumEms,1);
NumQuadPts = size(GAUSSQUAD.RECT,1);
GAUSSQUADRECT1 = ones(NumQuadPts,2) - GAUSSQUAD.RECT(:,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;
  jac = (CofB(:,3)+CofB(:,4)*xhat) .* (CofA(:,2)+CofA(:,4)*yhat)...
      - (CofA(:,3)+CofA(:,4)*xhat) .* (CofB(:,2)+CofB(:,4)*yhat);
  ProjPresEm = ProjPresEm + GAUSSQUAD.RECT(k,3) * jac .* EqnBC.fxnp(qp);
end
ProjPresEm = ProjPresEm./area;

%% Computing local L2 projection of velocity and then the related flux 
% Computing projection Gram matrix 
SX = zeros(NumEms,1);  SX2 = zeros(NumEms,1);
SY = zeros(NumEms,1);  SY2 = zeros(NumEms,1);
NumQuadPts = size(GAUSSQUAD.RECT,1);
GAUSSQUADRECT1 = ones(NumQuadPts,2) - GAUSSQUAD.RECT(:,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;
  jac = (CofB(:,3)+CofB(:,4)*xhat) .* (CofA(:,2)+CofA(:,4)*yhat)...
      - (CofA(:,3)+CofA(:,4)*xhat) .* (CofB(:,2)+CofB(:,4)*yhat);
  X = qp(:,1) - xc;
  Y = qp(:,2) - yc;
  SX = SX + GAUSSQUAD.RECT(k,3)*jac.*X;
  SY = SY + GAUSSQUAD.RECT(k,3)*jac.*Y;
  SX2 = SX2 + GAUSSQUAD.RECT(k,3)*jac.*(X.*X);
  SY2 = SY2 + GAUSSQUAD.RECT(k,3)*jac.*(Y.*Y);
end
ProjVelGM = zeros(NumEms,4,4);
ProjVelGM(:,1,1) = area;  ProjVelGM(:,2,2) = area;
ProjVelGM(:,1,3) = SX;    ProjVelGM(:,2,4) = SY;
ProjVelGM(:,3,1) = SX;    ProjVelGM(:,4,2) = SY;
ProjVelGM(:,3,3) = SX2;   ProjVelGM(:,4,4) = SY2;
% Computing projection RHS 
ProjVelRHS = zeros(NumEms,4);
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;
  jac = (CofB(:,3)+CofB(:,4)*xhat) .* (CofA(:,2)+CofA(:,4)*yhat)...
      - (CofA(:,3)+CofA(:,4)*xhat) .* (CofB(:,2)+CofB(:,4)*yhat);
  X = qp(:,1) - xc;
  Y = qp(:,2) - yc;
  vel = EqnBC.fxnu(qp);
  ProjVelRHS(:,1) = ProjVelRHS(:,1) + GAUSSQUAD.RECT(k,3)*jac.*vel(:,1);
  ProjVelRHS(:,2) = ProjVelRHS(:,2) + GAUSSQUAD.RECT(k,3)*jac.*vel(:,2);
  ProjVelRHS(:,3) = ProjVelRHS(:,3) + GAUSSQUAD.RECT(k,3)*jac.*vel(:,1).*X;  
  ProjVelRHS(:,4) = ProjVelRHS(:,4) + GAUSSQUAD.RECT(k,3)*jac.*vel(:,2).*Y;
end
% JL20161223: TO BE REVISED FOR EFFECIENCY
% Solving for (velocity) projection coefficients in RT[0] normalized basis 
ProjVelCofRT0 = zeros(NumEms,4);
for ie=1:NumEms
  EltGM = squeeze(ProjVelGM(ie,:,:));
  EltRHS = squeeze(ProjVelRHS(ie,:));
  ProjVelCofRT0(ie,:) = EltGM\EltRHS';
end

% %% Computing fluxes of the projected velocity: For each element: 4 edges
% ProjFlux = zeros(NumEms,4);
% for j=1:4
%   ProjFlux(:,j)...
%     = (ProjVelCofRT0NB(:,1)+ProjVelCofRT0NB(:,3).*Xm(:,j)).*EmEgNml(:,j,1)...
%     + (ProjVelCofRT0NB(:,2)+ProjVelCofRT0NB(:,4).*Ym(:,j)).*EmEgNml(:,j,2);  
% end

return;