function EltArray = WG_RectQ0Q0RT0_MixProdVecGradValBasFxn(...
  VecFxn,RectMesh,GAUSSQUAD) 
%% WG(Q0,Q0;RT[0]) on a rectangular mesh: Computing a mix inner product 
% James Liu, ColoState; 2012/07--2017/02 

%% Mesh info 
NumEms = RectMesh.NumEms;
x1 = RectMesh.node(RectMesh.elem(:,1),1);
y1 = RectMesh.node(RectMesh.elem(:,1),2);
x2 = RectMesh.node(RectMesh.elem(:,3),1);
y2 = RectMesh.node(RectMesh.elem(:,3),2);
area = (x2-x1).*(y2-y1);
xc = 0.5*(x1+x2);
yc = 0.5*(y1+y2);

%% 
% NOTE: CDWGB has size NumEms*5*4 
CDWGB = WG_RectQ0Q0RT0_CofRT0NmlzBas_DiscWkGradBasFxn(RectMesh);

%% Computing the elementwise 5x1 array (for the entire mesh) 
V1 = zeros(NumEms,1);  V2 = zeros(NumEms,1);  
VX = zeros(NumEms,1);  VY = zeros(NumEms,1);
NumQuadPts = size(GAUSSQUAD.RECT,1);
GAUSSQUADRECT1 = ones(NumQuadPts,2) - GAUSSQUAD.RECT(:,1:2);
qp = zeros(NumEms,2);
for k=1:NumQuadPts
  qp(:,1) = GAUSSQUAD.RECT(k,1)*x1 + GAUSSQUADRECT1(k,1)*x2;
  qp(:,2) = GAUSSQUAD.RECT(k,2)*y1 + GAUSSQUADRECT1(k,2)*y2;
  vec = VecFxn(qp);
  X = qp(:,1) - xc;
  Y = qp(:,2) - yc;
  V1 = V1 + GAUSSQUAD.RECT(k,3)*vec(:,1);
  V2 = V2 + GAUSSQUAD.RECT(k,3)*vec(:,2);
  VX = VX + GAUSSQUAD.RECT(k,3)*(vec(:,1).*X);
  VY = VY + GAUSSQUAD.RECT(k,3)*(vec(:,2).*Y);
end
V1 = V1.*area;  V2 = V2.*area;  VX = VX.*area;  VY = VY.*area;
% 
EltArray = zeros(NumEms,5);
for i=1:5
  EltArray(:,i) = sum(squeeze(CDWGB(:,i,:)).*[V1,V2,VX,VY],2);
end

return;