function [IntpltVelCofRT0] = Darcy_WG_QuadriQ0Q0RT0_IntpltVel(...
  EqnBC,QuadriMesh,GAUSSQUAD)
%% Darcy: WG(Q0,Q0,RT[0])Quadri: 
% Interpolating flux edgewise and then computing velocity elementwise 
% 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);
% Auxiliary quantities 
x = zeros(NumEms,5);
y = zeros(NumEms,5);
x(:,1) = QuadriMesh.node(k1,1);  y(:,1) = QuadriMesh.node(k1,2);
x(:,2) = QuadriMesh.node(k2,1);  y(:,2) = QuadriMesh.node(k2,2);
x(:,3) = QuadriMesh.node(k3,1);  y(:,3) = QuadriMesh.node(k3,2);
x(:,4) = QuadriMesh.node(k4,1);  y(:,4) = QuadriMesh.node(k4,2);
x(:,5) = x(:,1);  y(:,5) = y(:,1);
xc = 0.25*sum(x(:,1:4),2);
yc = 0.25*sum(y(:,1:4),2);
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

%% Computing the interpolated flux
% Interpolating flux edgewise 
IntpltFlux = zeros(NumEgs,1);
x1 = QuadriMesh.node(QuadriMesh.edge(:,1),1);
y1 = QuadriMesh.node(QuadriMesh.edge(:,1),2);
x2 = QuadriMesh.node(QuadriMesh.edge(:,2),1);
y2 = QuadriMesh.node(QuadriMesh.edge(:,2),2);
% Carrying edge-length info 
NmlEg = QuadriMesh.NmlEg .* [QuadriMesh.LenEg, QuadriMesh.LenEg];
NumQuadPts = size(GAUSSQUAD.LINE,1);
qp = zeros(NumEgs,2);
for k=1:NumQuadPts
  qp(:,1) = GAUSSQUAD.LINE(k,1)*x1 + GAUSSQUAD.LINE(k,2)*x2;
  qp(:,2) = GAUSSQUAD.LINE(k,1)*y1 + GAUSSQUAD.LINE(k,2)*y2;
  IntpltFlux = IntpltFlux + GAUSSQUAD.LINE(k,3)*dot(EqnBC.fxnu(qp),NmlEg,2);
end

%% Computing elementwise interpolated velocity based on the interpolated flux
% This is the mathematical interpolation operator \Pi_h 
IntpltMat = zeros(NumEms,4,4);  % For every element: four edges, four terms 
for j=1:4  % The j-th edge 
  IntpltMat(:,j,1) = EmEgNml(:,j,1);
  IntpltMat(:,j,2) = EmEgNml(:,j,2);
  IntpltMat(:,j,3) = EmEgNml(:,j,1).*Xm(:,j);
  IntpltMat(:,j,4) = EmEgNml(:,j,2).*Ym(:,j);
end
IntpltRHS = zeros(NumEms,4);
for j=1:4
  k = QuadriMesh.elem2edge(:,j);
  IntpltRHS(:,j) = IntpltFlux(k) .* QuadriMesh.SignEmEg(:,j);
end
% JL20161224: TO BE REVISED FOR EFFECIENCY 
% Solving for (interpolated velocity) coefficients in RT[0] normalized basis 
IntpltVelCofRT0 = zeros(NumEms,4);
for ie=1:NumEms
  EltMat = squeeze(IntpltMat(ie,:,:));
  EltRHS = squeeze(IntpltRHS(ie,:));
  IntpltVelCofRT0(ie,:) = EltMat\EltRHS';
end

return;