%% JL20171206: TO BE FINISHED/REVISED BY Graham 
function [NumerDsplEmAvg,NumerStrsEmAvg,NumerDivEmAvg]...
  = LinElas_CG_TriBR1_Prsn(lambda,mu,TriMesh,sln,GAUSSQUAD)
%% Computing emws. avg. numer. displacement, stress, dilation for presentation
% Bernardi-Raugel 1st order elements on triangles 
% James Liu, Graham Harper, ColoState; 2012/07--2017/12 

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

%% Auxiliary quantities 
% Elementwise cofficients in 9 BR1 basis functions for numerical solution 
cof = zeros(NumEms,9);
for j=1:3  % Treating vertices and edges in one loop 
  k = TriMesh.elem(:,j);
  cof(:,2*j-1) = sln(2*k-1);
  cof(:,2*j  ) = sln(2*k  );
  k = TriMesh.elem2edge(:,j);
  cof(:,6+j) = sln(2*NumNds+k) .* TriMesh.EmEgSign(:,j);
end

%% Computing numerical displacement elementwise averages 
NumerDsplEmAvg = zeros(NumEms,2);
NumQuadPts = size(GAUSSQUAD.TRIG,1);
for k=1:NumQuadPts
  BFVal = CG_TriBR1_BasFxnVal(TriMesh,...
    GAUSSQUAD.TRIG(k,1), GAUSSQUAD.TRIG(k,2), GAUSSQUAD.TRIG(k,3));
  NumerDsplEmAvg(:,1) = NumerDsplEmAvg(:,1) + GAUSSQUAD.TRIG(k,4)*sum(cof.*BFVal(:,:,1),2);
  NumerDsplEmAvg(:,2) = NumerDsplEmAvg(:,2) + GAUSSQUAD.TRIG(k,4)*sum(cof.*BFVal(:,:,2),2);
end

%% Computing numerical dilation (div of displacement) elementwise averages 
% NumerDivEmAvg = zeros(NumEms,1);
BasFxnDivAvg = zeros(NumEms,9);
CP1BasGrad = CG_TriP1_BasFxnGrad(TriMesh);
for j=1:3
  BasFxnDivAvg(:,2*j-1) = CP1BasGrad(:,j,1);
  BasFxnDivAvg(:,2*j  ) = CP1BasGrad(:,j,2);
end
BasFxnDivAvg(:,7) = TriMesh.EgLen(TriMesh.elem2edge(:,1))./(6*area);
BasFxnDivAvg(:,8) = TriMesh.EgLen(TriMesh.elem2edge(:,2))./(6*area);
BasFxnDivAvg(:,9) = TriMesh.EgLen(TriMesh.elem2edge(:,3))./(6*area);
NumerDivEmAvg = sum(cof.*BasFxnDivAvg,2);

%% Computing numerical stress elementwise averages 
NumerStrsEmAvg = zeros(NumEms,2,2);
NumQuadPts = size(GAUSSQUAD.TRIG,1);
for k=1:NumQuadPts
  [BFGrad,BFDiv] = CG_TriBR1_BasFxnGradDiv(TriMesh,...
    GAUSSQUAD.TRIG(k,1), GAUSSQUAD.TRIG(k,2), GAUSSQUAD.TRIG(k,3));
  StrnPhi = zeros(NumEms,9,2,2);
  StrnPhi(:,:,1,1) = BFGrad(:,:,1,1);
  StrnPhi(:,:,2,2) = BFGrad(:,:,2,2);
  StrnPhi(:,:,1,2) = 0.5*(BFGrad(:,:,1,2)+BFGrad(:,:,2,1));
  StrnPhi(:,:,2,1) = StrnPhi(:,:,1,2);
  NumerStrsEmAvg(:,1,1) = NumerStrsEmAvg(:,1,1) ...
    + ((2*mu)*sum(cof.*squeeze(StrnPhi(:,:,1,1)),2)+lambda*sum(cof.*BFDiv,2))...
    * GAUSSQUAD.TRIG(k,3);
  NumerStrsEmAvg(:,2,2) = NumerStrsEmAvg(:,2,2) ...
    + ((2*mu)*sum(cof.*squeeze(StrnPhi(:,:,2,2)),2)+lambda*sum(cof.*BFDiv,2))...
    * GAUSSQUAD.TRIG(k,3);
  NumerStrsEmAvg(:,1,2) = NumerStrsEmAvg(:,1,2) ...
    + (2*mu)* sum(cof.*squeeze(StrnPhi(:,:,1,2)),2) ...
    * GAUSSQUAD.TRIG(k,3);
  NumerStrsEmAvg(:,2,1) = NumerStrsEmAvg(:,1,2);
end


return;