%% JL20171206: TO BE FINISHED/REVISED BY Graham 
function [BFVal] = CG_TriBR1_BasFxnVal(TriMesh,alpha,beta,gamma)
%% FE: Triangular mesh: Bernardi-Raugel element (1st order)
% Basis function values at (alpha,beta,gamma) for all elements 
% James Liu, Graham Harper, ColoState; 2012/07--2017/12 

%% Mesh info 
NumEms = TriMesh.NumEms; 
k1 = TriMesh.elem(:,1);  k2 = TriMesh.elem(:,2);  k3 = TriMesh.elem(:,3);
x1 = TriMesh.node(k1,1);  y1 = TriMesh.node(k1,2); 
x2 = TriMesh.node(k2,1);  y2 = TriMesh.node(k2,2);
x3 = TriMesh.node(k3,1);  y3 = TriMesh.node(k3,2);
EgLen = zeros(NumEms,3);
EgLen(:,1) = sqrt((x2-x3).^2+(y2-y3).^2);
EgLen(:,2) = sqrt((x3-x1).^2+(y3-y1).^2);
EgLen(:,3) = sqrt((x1-x2).^2+(y1-y2).^2);

%% 
EmEgNml = zeros(NumEms,3,2);
EmEgNml(:,1,:) = [y3-y2, x2-x3]./[EgLen(:,1),EgLen(:,1)];
EmEgNml(:,2,:) = [y1-y3, x3-x1]./[EgLen(:,2),EgLen(:,2)];
EmEgNml(:,3,:) = [y2-y1, x1-x2]./[EgLen(:,3),EgLen(:,3)];

%% The value of each basis function is a 2-vector 
BFVal = zeros(NumEms,9,2);
% WO = ones(NumEms,1);
% For 3 vertices 
BFVal(:,1,1) = alpha; 
BFVal(:,2,2) = alpha; 
BFVal(:,3,1) =  beta;
BFVal(:,4,2) =  beta; 
BFVal(:,5,1) = gamma; 
BFVal(:,6,2) = gamma; 
% JL20171206: TO BE CHECKED BY Graham 
% For 3 edges
BFVal(:,7,:) =  beta * gamma * EmEgNml(:,1,:);
BFVal(:,8,:) = gamma * alpha * EmEgNml(:,2,:);
BFVal(:,9,:) = alpha * beta  * EmEgNml(:,3,:);

return;