function [ag,bg,cg] = NC_TriP1_CofBasFxn(TriMesh)
%% Nonconfirming FEs: Triangular mesh: P1: Elementwise bas.fxn.coeffs.
% 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*abs((x2-x1).*(y3-y1)-(x3-x1).*(y2-y1));

%% Generating info for edge midpoints 
% Elementwise edge midpoints x-, y- coordinates 
xm1 = 0.5*(x2+x3);  xm2 = 0.5*(x3+x1);  xm3 = 0.5*(x1+x2);
ym1 = 0.5*(y2+y3);  ym2 = 0.5*(y3+y1);  ym3 = 0.5*(y1+y2);
% Determinants for Lagrangian basis functions 
% mdet = 0.25*area;
mdet = (xm2.*ym3-xm3.*ym2) + (xm3.*ym1-xm1.*ym3) + (xm1.*ym2-xm2.*ym1);

%% Computing coefficients of elementwise 3 basis functions: a+bx+cy 
% For each element: 3 edge-oriented bas.fxns., 3 coefficients a,b,c 
ag = zeros(NumEms,3);  % for const. term 
bg = zeros(NumEms,3);  % for x term 
cg = zeros(NumEms,3);  % for y term 
ag(:,1) = (xm2.*ym3-xm3.*ym2)./mdet;
bg(:,1) = (ym2-ym3)./mdet;
cg(:,1) = (xm3-xm2)./mdet;
ag(:,2) = (xm3.*ym1-xm1.*ym3)./mdet;
bg(:,2) = (ym3-ym1)./mdet;
cg(:,2) = (xm1-xm3)./mdet;
ag(:,3) = (xm1.*ym2-xm2.*ym1)./mdet;
bg(:,3) = (ym1-ym2)./mdet;
cg(:,3) = (xm2-xm1)./mdet;

return;