function LMCR = Parab_IE_CG_TriP1_LMCR(EqnBIC,TriMesh,tyme,sln,GAUSSQUAD)
%% Computing elementwise local-mass-conservation residual 
% for Implicit Euler + Continuous Galerkin on a tri.mesh with P1 elements 
% James Liu, ColoState; 2012/07--2017/02 

%% EqnBIC 
D = EqnBIC.D;

%% Mesh info 
NumEms = TriMesh.NumEms;
k1 = TriMesh.elem(:,1);  k2 = TriMesh.elem(:,2);  k3 = TriMesh.elem(:,3);
LenEg = TriMesh.LenEg;
NmlEg = TriMesh.NmlEg;
e2g = TriMesh.elem2edge;

%% Timing info 
NT = tyme.NT;  Deltat = tyme.Deltat;  tm = tyme.tm;

%% Computing elementwise concentration averages 
EltMass = zeros(NumEms,NT+1);
for j=1:(NT+1)
  c1 = sln(k1,j);  c2 = sln(k2,j);  c3 = sln(k3,j);
  EltMass(:,j) = (c1+c2+c3)/3 .* TriMesh.area;
end

%% Computing element-edgewise outward normal fluxes 
% CP1BasGrad = zeros(NumEms,3,2);
CP1BasGrad = CG_TriP1_BasFxnGrad(TriMesh);
EltEgNmlFlux = zeros(NumEms,3,NT+1);
for j=1:(NT+1)
  c1 = sln(k1,j);  c2 = sln(k2,j);  c3 = sln(k3,j);
  % NumerSlnGrad = zeros(NumEms,2);
  NumerSlnGrad = [c1,c1] .* squeeze(CP1BasGrad(:,1,:))...
               + [c2,c2] .* squeeze(CP1BasGrad(:,2,:))... 
               + [c3,c3] .* squeeze(CP1BasGrad(:,3,:));
  for k=1:3
    EltEgNmlFlux(:,k,j) = Deltat * dot(-D*NumerSlnGrad, ...
      NmlEg(e2g(:,k),:).*[LenEg(e2g(:,k)),LenEg(e2g(:,k))], 2);
  end
end
    
%% Computing elementwise source contribution at all time moments 
EltSrcCtrbn = zeros(NumEms,NT+1);
NumQuadPts = size(GAUSSQUAD.TRIG,1);
for j=1:(NT+1)
  EltSrc = zeros(NumEms,1);
  for k=1:NumQuadPts
    qp = GAUSSQUAD.TRIG(k,1) * TriMesh.node(TriMesh.elem(:,1),:)...
       + GAUSSQUAD.TRIG(k,2) * TriMesh.node(TriMesh.elem(:,2),:)...
       + GAUSSQUAD.TRIG(k,3) * TriMesh.node(TriMesh.elem(:,3),:);
    EltSrc = EltSrc + GAUSSQUAD.TRIG(k,4)*EqnBIC.fxnf(qp,tm(j));
  end
  EltSrc = EltSrc .* TriMesh.area;
  EltSrcCtrbn(:,j) = Deltat * EltSrc; 
end

%% Local-mass-conservation residual 
LMCR = zeros(NumEms,NT);
for j=1:NT 
  LMCR(:,j) = EltMass(:,j+1) - EltMass(:,j+1) ...
    + sum(EltEgNmlFlux(:,:,j+1),2) - EltSrcCtrbn(:,j+1);
end

return;