function [Fval,Fdiv] = Hdiv_QuadriAC0_NmlzPiolaBas_FxnValDiv(QuadriMesh,xhat,yhat)
%% Quadrilateral mesh: 
% James Liu, ColoState; 2017/02--2017/11 

%% 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 
Fval = zeros(NumEms,4,2);
Fval(:,1,1) = 1;
Fval(:,2,2) = 1;
Fval(:,3,1) = X;  Fval(:,3,2) = Y;
Fval(:,4,1) = MatP(:,1,1) * xhat - MatP(:,1,2) * yhat;
Fval(:,4,2) = MatP(:,2,1) * xhat - MatP(:,2,2) * yhat;

%% Basis function divergence  
Fdiv = zeros(NumEms,4);
Fdiv(:,3) = 2;

return;