function EltArrayUpwnd = WG_TriP0P0RT0_Upwind(VecFxn,TriMesh,GAUSSQUAD) 
% WG(P0,P0;RT0) on a triangular mesh: Elementwise array for upwinding  
% NOTE: EltAryUpwnd has size zeros NumEms*3 
% James Liu, ColoState; 2012/07--2017/02 

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

%% Computing edgewise normal fluxes using a Gaussian quadrature 
EgNmlFlux = zeros(NumEgs,1);
x1 = TriMesh.node(TriMesh.edge(:,1),1);
y1 = TriMesh.node(TriMesh.edge(:,1),2);
x2 = TriMesh.node(TriMesh.edge(:,2),1);
y2 = TriMesh.node(TriMesh.edge(:,2),2);
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;
  EgNmlFlux = EgNmlFlux + GAUSSQUAD.LINE(k,3)...
    * dot(VecFxn(qp),TriMesh.NmlEg,2);
end
EgNmlFlux = EgNmlFlux .* TriMesh.LenEg;

%% Converting to elementwise normal fluxes 
NmlFlux = zeros(NumEms,3);
for j=1:3
  NmlFlux(:,j) = TriMesh.SignEmEg(:,j)...
    .* EgNmlFlux(TriMesh.elem2edge(:,j));
end

%% Upwinding 
EltArrayUpwnd = zeros(NumEms,3);
for j=1:3
  idx = find(NmlFlux(:,j)>0);
  EltArrayUpwnd(idx,j) = -NmlFlux(idx,j);
end

return;