function FxnAvg = Hdiv_QuadriAC02_NmlzPiolaBas_FxnValAvg(QuadriMesh,xhat,yhat)
%% Quadrilateral mesh: 
% James Liu, ColoState; 2017/02--2017/11 

%% Mesh info 
NumEms = QuadriMesh.NumEms; 

%% Vector version 
% Fval = zeros(NumEms,4,2);
% Fdiv = zeros(NumEms,4);
[Fval,Fdiv] = Hdiv_QuadriAC0_NmlzPiolaBas_FxnValDiv(QuadriMesh,xhat,yhat);

%% Matrix version: Basis function values 
FxnVal = zeros(NumEms,8,2,2);
for i=1:4
  FxnVal(:,i,  1,:) = Fval(:,i,:);
  FxnVal(:,i+4,2,:) = Fval(:,i,:);
end

%% Matrix version: Basis function value averages  
FxnAvg = zeros(NumEms,8,2,2);
for i=1:8
  FxnAvg(:,i,1,1) = FxnVal(:,i,1,1);
  FxnAvg(:,i,2,2) = FxnVal(:,i,2,2);
  FxnAvg(:,i,1,2) = 0.5 * ( FxnVal(:,i,1,2) + FxnVal(:,i,2,1) );
  FxnAvg(:,i,2,1) = FxnAvg(:,i,1,2);
end

return;