function BasFxnVal = CG_RectQ1_BasFxnVal(RectMesh,pt)
%% CG: Rect.mesh. Q1 elements: 4 basis function values (for all elements) 
% Assuming pt has size NumEms*2 (one point in each element) 
% James Liu, ColoState; 2012/07--2017/02 

%% Mesh info 
NumEms = RectMesh.NumEms;
k1 = RectMesh.elem(:,1);  k3 = RectMesh.elem(:,3);
x1 = RectMesh.node(k1,1);  y1 = RectMesh.node(k1,2);
x2 = RectMesh.node(k3,1);  y2 = RectMesh.node(k3,2);

%% Mapping back to the reference rectangle (unit square)
x = pt(:,1);  X = (x-x1)./(x2-x1);
y = pt(:,2);  Y = (y-y1)./(y2-y1);

%% 4 corners/vertices: Starting from lower-left, going counterclockwise 
BasFxnVal = zeros(NumEms,4);
BasFxnVal(:,1) = (1-X) .* (1-Y);
BasFxnVal(:,2) =     X .* (1-Y);
BasFxnVal(:,3) =     X .* Y;
BasFxnVal(:,4) = (1-X) .* Y;

return;