%% JL2020117: Need to finish computing NumerStrsCofP022 
function [NumerDispCofP12, NumerDivCofP0, NumerStrsCofP022] ...
  = LinElas_CG_TriP12_DispDivStrs(TriMesh, sln)
%% Computing numerical displacement, divergence, stress in resp. bases 
% James Liu, ColoState; 2012/07--2020/11 

%% Mesh info 
% NumNds = TriMesh.NumNds;
NumEms = TriMesh.NumEms;

%% Obtaining elementwise coefficients for the numerical solution 
cof = zeros(NumEms,6);
for j=1:3
  k = TriMesh.elem(:,j);
  cof(:,2*j-1) = sln(2*k-1);
  cof(:,2*j  ) = sln(2*k  );
end

%% Obtaining auxiliary quantities 
% CP1BasGrad = zeros(NumEms,3,2);
CP1BasGrad = CG_TriP1_BasFxnGrad(TriMesh);

%% Computing numerical displacement: Emws.cofs. in P1^2 nmlz.bas.: 1,X,Y 
NumerDispCofP12 = zeros(NumEms,6);
NumerDispCofP12(:,1) = (1/3)*sum(sln(2*TriMesh.elem(:,[1,2,3])-1),2);
NumerDispCofP12(:,2) = (1/3)*sum(sln(2*TriMesh.elem(:,[1,2,3])),  2);

%% Computing numerical divergence/dilation: Emws.cofs. in P0 
% NumerDivCofP0 = zeros(NumEms,1);
NumerDivCofP0 = sum(cof(:,[1,3,5]).*squeeze(CP1BasGrad(:,:,1)),2)...
              + sum(cof(:,[2,4,6]).*squeeze(CP1BasGrad(:,:,2)),2);

%% Computing numerical stress: 
% Elementwise coefficients in P0^{2x2} 
NumerStrsCofP022 = zeros(NumEms,4);

return;