function BaryCrd = TriMesh_BaryCrd(TriMesh,pt)
%% Triangular mesh: Barycentric coordinates of points
% Assuming one point inside each element 
% James Liu, ColoState; 2012/07--2017/02 

%% 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*((x2-x1).*(y3-y1)-(x3-x1).*(y2-y1));

%% Barycentric coordinates
BaryCrd = zeros(NumEms,3);
x = pt(:,1);  y = pt(:,2);
area1 = 0.5*((x2-x).*(y3-y)-(x3-x).*(y2-y));
area2 = 0.5*((x-x1).*(y3-y1)-(x3-x1).*(y-y1));
area3 = 0.5*((x2-x1).*(y-y1)-(x-x1).*(y2-y1));
BaryCrd(:,1) = area1./area;
BaryCrd(:,2) = area2./area;
BaryCrd(:,3) = area3./area;

return;