function FV = Hdiv_QuadriAC1_NmlzPiolaBas_FxnVal(QuadriMesh,xhat,yhat)
%% Quadrilateral mesh: 
% James Liu, ColoState; 2017/02--2017/10 

%% Mesh info 
NumEms = QuadriMesh.NumEms; 
xc = QuadriMesh.EmCntr(:,1);
yc = QuadriMesh.EmCntr(:,2);
CofA = QuadriMesh.CofA;
CofB = QuadriMesh.CofB;

%% Jacobian and Piola 
MatJ = zeros(NumEms,2,2);
MatP = zeros(NumEms,2,2);
MatJ(:,1,1) = CofA(:,2) + CofA(:,4)*yhat;  
MatJ(:,1,2) = CofA(:,3) + CofA(:,4)*xhat;
MatJ(:,2,1) = CofB(:,2) + CofB(:,4)*yhat;
MatJ(:,2,2) = CofB(:,3) + CofB(:,4)*xhat;
jac = MatJ(:,1,1).*MatJ(:,2,2) - MatJ(:,1,2).*MatJ(:,2,1);
for i=1:2
  for j=1:2
    MatP(:,i,j) = MatJ(:,i,j)./jac;
  end
end

%% Normalized coordinates 
qp = zeros(NumEms,2);
qp(:,1) = CofA(:,1) + CofA(:,2)*xhat + CofA(:,3)*yhat + CofA(:,4)*xhat*yhat;
qp(:,2) = CofB(:,1) + CofB(:,2)*xhat + CofB(:,3)*yhat + CofB(:,4)*xhat*yhat;
X = qp(:,1) - xc;
Y = qp(:,2) - yc;

%% Basis function values 
FV = zeros(NumEms,10,2);
% Part I: Polynomial P1^2 
FV(:,1,1) = 1;  FV(:,2,1) = X;  FV(:,3,1) = Y;
FV(:,4,2) = 1;  FV(:,5,2) = X;  FV(:,6,2) = Y;
% Part II: Homogeneous polynomial [X,Y]*P1tilde 
FV(:,7,1) = X.^2;  FV(:,8,1) = X.*Y;  
FV(:,7,2) = X.*Y;  FV(:,8,2) = Y.^2;
% Part III: Supplementary S1 
FV(:,9,1) = MatP(:,1,1)*(1-xhat^2) + MatP(:,1,2)*(2*xhat*yhat);
FV(:,9,2) = MatP(:,2,1)*(1-xhat^2) + MatP(:,2,2)*(2*xhat*yhat);
FV(:,10,1) = MatP(:,1,1)*(2*xhat*yhat) + MatP(:,1,2)*(1-yhat^2);
FV(:,10,2) = MatP(:,2,1)*(2*xhat*yhat) + MatP(:,2,2)*(1-yhat^2);

return;