function [NumerDispEmAvg,NumerStrsEmAvg,NumerDivEmAvg] = ...
  LinElas_CG_QuadriEQ1_Presn(lambda,mu,QuadriMesh,sln,GAUSSQUAD)
%% LinElas.: CG: Enriched Lagrangian Q1 on quadri.mesh (w/ red.integration)
% Computing em.avg. numer. displacement, stress, dilation for presentation
% James Liu, ColoState; 2012/07--2020/08 

%% Mesh info 
NumNds = QuadriMesh.NumNds;
NumEms = QuadriMesh.NumEms;
% NumEgs = QuadriMesh.NumEgs;
area = QuadriMesh.area;
CofA = QuadriMesh.CofA;
CofB = QuadriMesh.CofB;

%% Retrieving coefficients 
% Elementwise cofficients in 12 BR1 basis functions for numerical solution 
cof = zeros(NumEms,12);
for j=1:4  % Treating vertices and edges in one loop 
  k = QuadriMesh.elem(:,j);
  cof(:,2*j-1) = sln(2*k-1);
  cof(:,2*j  ) = sln(2*k  );
  k = QuadriMesh.elem2edge(:,j);
  cof(:,8+j) = sln(2*NumNds+k) .* QuadriMesh.EmEgSign(:,j);
end

%% Computing numerical displacement elementwise averages
NumerDispEmAvg = zeros(NumEms,2);
NumQuadPts = size(GAUSSQUAD.RECT,1);
GAUSSQUADRECT1 = ones(NumQuadPts,2) - GAUSSQUAD.RECT(:,1:2);
for k=1:NumQuadPts
  xhat = GAUSSQUADRECT1(k,1);
  yhat = GAUSSQUADRECT1(k,2);
  jac = (CofB(:,3)+CofB(:,4)*xhat) .* (CofA(:,2)+CofA(:,4)*yhat)...
      - (CofA(:,3)+CofA(:,4)*xhat) .* (CofB(:,2)+CofB(:,4)*yhat);
  BFVal = CG_QuadriBR1_BasFxnVal(QuadriMesh,xhat,yhat);
  NumerDispEmAvg(:,1) = NumerDispEmAvg(:,1) + GAUSSQUAD.RECT(k,3) ...
    * jac ./area .* sum(cof.*squeeze(BFVal(:,:,1)),2);
  NumerDispEmAvg(:,2) = NumerDispEmAvg(:,2) + GAUSSQUAD.RECT(k,3) ...
    * jac ./area .* sum(cof.*squeeze(BFVal(:,:,2)),2);
end

%% Computing numerical dilation (div of displacement) elementwise averages 
NumerDivEmAvg = zeros(NumEms,1);
NumQuadPts = size(GAUSSQUAD.RECT,1);
GAUSSQUADRECT1 = ones(NumQuadPts,2) - GAUSSQUAD.RECT(:,1:2);
for k=1:NumQuadPts
  xhat = GAUSSQUADRECT1(k,1);
  yhat = GAUSSQUADRECT1(k,2);
  jac = (CofB(:,3)+CofB(:,4)*xhat) .* (CofA(:,2)+CofA(:,4)*yhat)...
      - (CofA(:,3)+CofA(:,4)*xhat) .* (CofB(:,2)+CofB(:,4)*yhat);
  [BFGrad,BFDiv] = CG_QuadriBR1_BasFxnGradDiv(QuadriMesh,xhat,yhat);
  NumerDivEmAvg(:,1) = NumerDivEmAvg(:,1) + GAUSSQUAD.RECT(k,3) ...
    * jac ./area .* sum(cof.*BFDiv,2);
end

%% Computing numerical stress elementwise averages 
NumerStrsEmAvg = zeros(NumEms,2,2);
NumQuadPts = size(GAUSSQUAD.RECT,1);
GAUSSQUADRECT1 = ones(NumQuadPts,2) - GAUSSQUAD.RECT(:,1:2);
for k=1:NumQuadPts
  xhat = GAUSSQUADRECT1(k,1);
  yhat = GAUSSQUADRECT1(k,2);
  jac = (CofB(:,3)+CofB(:,4)*xhat) .* (CofA(:,2)+CofA(:,4)*yhat)...
      - (CofA(:,3)+CofA(:,4)*xhat) .* (CofB(:,2)+CofB(:,4)*yhat);
  [BFGrad,BFDiv] = CG_QuadriBR1_BasFxnGradDiv(QuadriMesh,xhat,yhat);
  StrnPhi = zeros(NumEms,12,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))...
    .* jac ./ area * GAUSSQUAD.RECT(k,3);
  NumerStrsEmAvg(:,2,2) = NumerStrsEmAvg(:,2,2) ...
    + ((2*mu)*sum(cof.*squeeze(StrnPhi(:,:,2,2)),2)+lambda*sum(cof.*BFDiv,2))...
    .* jac ./ area * GAUSSQUAD.RECT(k,3);
  NumerStrsEmAvg(:,1,2) = NumerStrsEmAvg(:,1,2) ...
    + (2*mu)* sum(cof.*squeeze(StrnPhi(:,:,1,2)),2) ...
    .* jac ./ area * GAUSSQUAD.RECT(k,3);
  NumerStrsEmAvg(:,2,1) = NumerStrsEmAvg(:,1,2);
end

return;