%% JL20171223: TO BE FINISHED/REVISED 
function [BFGrad,BFDiv] = CG_QuadriBR1_BasFxnGradDiv(QuadriMesh,xhat,yhat)
%% CG: Quadrilateral mesh: Bernardi-Raugel elements (1st order)
% Basis function gradient and divergence at (xhat,yhat) for all elements 
% James Liu, ColoState; 2012/07--2017/12 

%% Mesh info 
NumEms = QuadriMesh.NumEms; 
CofA = QuadriMesh.CofA;
CofB = QuadriMesh.CofB;

%% 
k1 = QuadriMesh.elem(:,1);  k2 = QuadriMesh.elem(:,2);  
k3 = QuadriMesh.elem(:,3);  k4 = QuadriMesh.elem(:,4);  
x1 = QuadriMesh.node(k1,1);  y1 = QuadriMesh.node(k1,2); 
x2 = QuadriMesh.node(k2,1);  y2 = QuadriMesh.node(k2,2);
x3 = QuadriMesh.node(k3,1);  y3 = QuadriMesh.node(k3,2);
x4 = QuadriMesh.node(k4,1);  y4 = QuadriMesh.node(k4,2);
EgLen = zeros(NumEms,4);
EgLen(:,1) = sqrt((x2-x1).^2+(y2-y1).^2);
EgLen(:,2) = sqrt((x3-x2).^2+(y3-y2).^2);
EgLen(:,3) = sqrt((x4-x3).^2+(y4-y3).^2);
EgLen(:,4) = sqrt((x1-x4).^2+(y1-y4).^2);
EmEgNml = zeros(NumEms,4,2);
EmEgNml(:,1,:) = [y2-y1, x1-x2]./[EgLen(:,1),EgLen(:,1)];
EmEgNml(:,2,:) = [y3-y2, x2-x3]./[EgLen(:,2),EgLen(:,2)];
EmEgNml(:,3,:) = [y4-y3, x3-x4]./[EgLen(:,3),EgLen(:,3)];
EmEgNml(:,4,:) = [y1-y4, x4-x1]./[EgLen(:,4),EgLen(:,4)];

%% Auxiliary quantities for Jacobian and its inverse transpose 
MatJ = 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);
JIT = zeros(NumEms,2,2);  % Jacobian Inverse Transpose 
JIT(:,1,1) =  MatJ(:,2,2)./jac;  JIT(:,1,2) = -MatJ(:,2,1)./jac;
JIT(:,2,1) = -MatJ(:,1,2)./jac;  JIT(:,2,2) =  MatJ(:,1,1)./jac;

%% Auxiliary quantities for standard basis functions: hatphi, hatpsi 
gradhatphi = zeros(NumEms,4,2);
gradhatphi(:,1,1) = -1 + yhat;  gradhatphi(:,1,2) = -1 + xhat;
gradhatphi(:,2,1) =  1 - yhat;  gradhatphi(:,2,2) =    - xhat;
gradhatphi(:,3,1) =      yhat;  gradhatphi(:,3,2) =      xhat;
gradhatphi(:,4,1) =    - yhat;  gradhatphi(:,4,2) =  1 - xhat;
gradhatpsi = zeros(NumEms,4,2);
gradhatpsi(:,1,1) = (1-2*xhat)*(1-yhat);  gradhatpsi(:,1,2) = -xhat*(1-xhat);
gradhatpsi(:,2,1) = yhat*(1-yhat);        gradhatpsi(:,2,2) = xhat*(1-2*yhat);
gradhatpsi(:,3,1) = (1-2*xhat)*yhat;      gradhatpsi(:,3,2) = xhat*(1-xhat);
gradhatpsi(:,4,1) = -yhat*(1-yhat);       gradhatpsi(:,4,2) = (1-xhat)*(1-2*yhat);

%% Auxiliary quantities for basis functions: phi, psi
gradphi = zeros(NumEms,4,2);
for j=1:4
  gradphi(:,j,1) = sum(JIT(:,1,:).*gradhatphi(:,j,:),3);
  gradphi(:,j,2) = sum(JIT(:,2,:).*gradhatphi(:,j,:),3);
end
gradpsi = zeros(NumEms,4,2);
for j=1:4
  gradpsi(:,j,1) = sum(JIT(:,1,:).*gradhatpsi(:,j,:),3);
  gradpsi(:,j,2) = sum(JIT(:,2,:).*gradhatpsi(:,j,:),3);
end

%% The gradient of each basis function is a 2x2 matrix 
BFGrad = zeros(NumEms,12,2,2);
% For 4 vertices 
for j=1:4
  BFGrad(:,2*j-1,1,1) = gradphi(:,j,1);  
  BFGrad(:,2*j-1,1,2) = gradphi(:,j,2);
  BFGrad(:,2*j,  2,1) = gradphi(:,j,1);  
  BFGrad(:,2*j,  2,2) = gradphi(:,j,2);  
end
% For 4 edges
for j=1:4
  % k = QuadriMesh.elem2edge(:,j);
  BFGrad(:,8+j,1,1) = squeeze(EmEgNml(:,j,1)).*squeeze(gradpsi(:,j,1));
  BFGrad(:,8+j,1,2) = squeeze(EmEgNml(:,j,1)).*squeeze(gradpsi(:,j,2));
  BFGrad(:,8+j,2,1) = squeeze(EmEgNml(:,j,2)).*squeeze(gradpsi(:,j,1));
  BFGrad(:,8+j,2,2) = squeeze(EmEgNml(:,j,2)).*squeeze(gradpsi(:,j,2));
end

%% The divergence of each basis function is just a scalar 
BFDiv = zeros(NumEms,12);
% For 4 vertices 
for j=1:4
  BFDiv(:,2*j-1) = gradphi(:,j,1);
  BFDiv(:,2*j  ) = gradphi(:,j,2);  
end
% For 4 edges
for j=1:4 
  % k = QuadriMesh.elem2edge(:,j);
  BFDiv(:,8+j) = sum(squeeze(EmEgNml(:,j,:)).*squeeze(gradpsi(:,j,:)),2);
end

return;