function EltArray = WG_TriP0P0RT0_MixProdVecGradValBasFxn(...
  VecFxn,TriMesh,GAUSSQUAD) 
%% WG(Q0,Q0;RT0) on a triangular mesh: Computing a mix inner product 
% James Liu, ColoState; 2012/07--2017/02 

%% Mesh info 
NumEms = TriMesh.NumEms;
k1 = TriMesh.elem(:,1);  k2 = TriMesh.elem(:,2);  k3 = TriMesh.elem(:,3);
x1 = TriMesh.node(k1,1);  y1 = TriMesh.node(k1,2); 
x2 = TriMesh.node(k2,1);  y2 = TriMesh.node(k2,2);
x3 = TriMesh.node(k3,1);  y3 = TriMesh.node(k3,2);
area = 0.5*abs((x2-x1).*(y3-y1)-(x3-x1).*(y2-y1));
xc = (1/3)*(x1+x2+x3);
yc = (1/3)*(y1+y2+y3);

%% 
% NOTE: CDWGB has size NumEms*4*3 
CDWGB = WG_TriP0P0RT0_CofRT0NmlzBas_DiscWkGradBasFxn(TriMesh);

%% Computing the elementwise 4x1 array (for the entire mesh) 
V1 = zeros(NumEms,1);  V2 = zeros(NumEms,1);  VXY = zeros(NumEms,1);
NumQuadPts = size(GAUSSQUAD.TRIG,1);
% qp = zeros(NumEms,2);
for k=1:NumQuadPts
  qp = GAUSSQUAD.TRIG(k,1) * TriMesh.node(TriMesh.elem(:,1),:)...
     + GAUSSQUAD.TRIG(k,2) * TriMesh.node(TriMesh.elem(:,2),:)...
     + GAUSSQUAD.TRIG(k,3) * TriMesh.node(TriMesh.elem(:,3),:);
  vec = VecFxn(qp);
  X = qp(:,1) - xc;
  Y = qp(:,2) - yc;
  V1 = V1 + GAUSSQUAD.TRIG(k,4)*vec(:,1);
  V2 = V2 + GAUSSQUAD.TRIG(k,4)*vec(:,2);
  VXY = VXY + GAUSSQUAD.TRIG(k,4)*(vec(:,1).*X+vec(:,2).*Y);
end
V1 = V1.*area;  V2 = V2.*area;  VXY = VXY.*area;
% 
EltArray = zeros(NumEms,4);
for i=1:4
  EltArray(:,i) = sum(squeeze(CDWGB(:,i,:)).*[V1,V2,VXY],2);
end

return;