%% JL20171206: TO BE FINISHED/REVISED BY Graham 
function [BFGrad,BFDiv] = CG_TriBR1_BasFxnGradDiv(TriMesh,alpha,beta,gamma)
%% FE: Triangular mesh: Bernardi-Raugel element (1st order)
% Basis function gradient & divergence 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);
area = 0.5*abs((x2-x1).*(y3-y1)-(x3-x1).*(y2-y1));
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);

%% NOTE: Lagrangian P1 bas.fxn.val. = Barycentric coordinates 
CP1BasVal = [alpha,beta,gamma];
CP1BasGrad = zeros(NumEms,3,2);
CP1BasGrad(:,1,:) = [(y2-y3)./(2*area), (x3-x2)./(2*area)];
CP1BasGrad(:,2,:) = [(y3-y1)./(2*area), (x1-x3)./(2*area)];
CP1BasGrad(:,3,:) = [(y1-y2)./(2*area), (x2-x1)./(2*area)];

%%
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 gradient of each basis function is a 2x2 matrix 
BFGrad = zeros(NumEms,9,2,2);
% For 3 vertices 
for j=1:3
  BFGrad(:,2*j-1,1,1) = CP1BasGrad(:,j,1);
  BFGrad(:,2*j-1,1,2) = CP1BasGrad(:,j,2);
  BFGrad(:,2*j,  2,1) = CP1BasGrad(:,j,1);
  BFGrad(:,2*j,  2,2) = CP1BasGrad(:,j,2);
end
% JL20171206: TO BE FINISHED BY Graham 
% For 3 edges
idx1_list = [2,3,1];
idx2_list = [3,1,2];
for j=1:3
  % GH20171219: this is the most efficient method I can think of to keep track of indices
  idx1 = idx1_list(j);
  idx2 = idx2_list(j);
  % gradient goes across rows, normal goes down columns
  BFGrad(:,j+6,1,1) = (CP1BasVal(idx1)*CP1BasGrad(:,idx2,1)+CP1BasVal(idx2)*CP1BasGrad(:,idx1,1)).*EmEgNml(:,j,1);
  BFGrad(:,j+6,1,2) = (CP1BasVal(idx1)*CP1BasGrad(:,idx2,2)+CP1BasVal(idx2)*CP1BasGrad(:,idx1,2)).*EmEgNml(:,j,1);
  BFGrad(:,j+6,2,1) = (CP1BasVal(idx1)*CP1BasGrad(:,idx2,1)+CP1BasVal(idx2)*CP1BasGrad(:,idx1,1)).*EmEgNml(:,j,2);
  BFGrad(:,j+6,2,2) = (CP1BasVal(idx1)*CP1BasGrad(:,idx2,2)+CP1BasVal(idx2)*CP1BasGrad(:,idx1,2)).*EmEgNml(:,j,2);
end

%% The divergence of each basis function is a scalar 
BFDiv = zeros(NumEms,9);
% For 3 vertices 
for j=1:3
  BFDiv(:,2*j-1) = CP1BasGrad(:,j,1);
  BFDiv(:,2*j  ) = CP1BasGrad(:,j,2);
end
% JL20171206: TO BE FINISHED BY Graham 
% For 3 edges 
% GH20171217: No need to loop since the formulas are simple
BFDiv(:,7) = gamma* sum(CP1BasGrad(:,2,:).*EmEgNml(:,1,:),3)  +  beta*  sum(CP1BasGrad(:,3,:).*EmEgNml(:,1,:),3);
BFDiv(:,8) = gamma* sum(CP1BasGrad(:,1,:).*EmEgNml(:,2,:),3)  +  alpha* sum(CP1BasGrad(:,3,:).*EmEgNml(:,2,:),3);
BFDiv(:,9) = beta*  sum(CP1BasGrad(:,1,:).*EmEgNml(:,3,:),3)  +  alpha* sum(CP1BasGrad(:,2,:).*EmEgNml(:,3,:),3);

return;