%% JL20171223: TO BE FINISHED/REVISED 
function [BFVal] = CG_QuadriBR1_BasFxnVal(QuadriMesh,xhat,yhat)
%% CG: Quadriangular mesh: Bernardi-Raugel element (1st order)
% Basis function values at (xhat,yhat) for all elements 
% James Liu, ColoState; 2012/07--2017/12 

%% Mesh info 
NumEms = QuadriMesh.NumEms; 
NumEgs = QuadriMesh.NumEgs;
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)];

%% The value of each basis function is a 2-vector 
BFVal = zeros(NumEms,12,2);
% For 4 vertices 
BFVal(:,1,1) = (1-xhat)*(1-yhat);  BFVal(:,2,2) = (1-xhat)*(1-yhat);
BFVal(:,3,1) =    xhat *(1-yhat);  BFVal(:,4,2) =    xhat *(1-yhat);
BFVal(:,5,1) =    xhat *   yhat ;  BFVal(:,6,2) =    xhat *   yhat ;
BFVal(:,7,1) = (1-xhat)*   yhat ;  BFVal(:,8,2) = (1-xhat)*   yhat ;
% For 4 edges
BFVal(:, 9,:) = xhat * (1-xhat) *        (1-yhat) * EmEgNml(:,1,:);
BFVal(:,10,:) = xhat *            yhat * (1-yhat) * EmEgNml(:,2,:);
BFVal(:,11,:) = xhat * (1-xhat) * yhat            * EmEgNml(:,3,:);
BFVal(:,12,:) =        (1-xhat) * yhat * (1-yhat) * EmEgNml(:,4,:);

return;