%% JL20210207: TO BE FINISHED 
function [NumerDispEmAvg,NumerStrsEmAvg,NumerDivEmAvg] ...
  = LinElas_WG_QuadriP02P02AC02P0_Presn(lambda,mu,QuadriMesh,sln,GAUSSQUAD) 
%% Comput.quant. for preseentation: Numer.disp., stress, div. 
% James Liu, ColoState; 2017/02--2021/02 

%% Mesh info 
NumEms = QuadriMesh.NumEms;

%% Elementwise coefficients in WG(P0^2,P0^2) basis functions 
cof = zeros(NumEms,10);
cof(:,1) = sln(2*(0:NumEms-1)+1);
cof(:,2) = sln(2*(0:NumEms-1)+2);
for j=1:4
  k = QuadriMesh.elem2edge(:,j);
  cof(:,2+2*j-1) = sln(2*NumEms+2*k-1);
  cof(:,2+2*j  ) = sln(2*NumEms+2*k  );
end

%% "Computing"/extracting numerical displacement element "averages" 
NumerDispEmAvg = zeros(NumEms,2);
NumerDispEmAvg(:,1) = cof(:,1);
NumerDispEmAvg(:,2) = cof(:,2);

%% JL20210207: TO BE FINISHED 
%% Computing numerical stress at element centers  
NumerStrsEmAvg = zeros(NumEms,2,2);

%% Computing numer.div.(em.avg.) in P0; Note CDWDB size NumEms*10 
CDWDB = WG_QuadriP02P02P0_CofNmlzBas_DiscWkDivBasFxn(QuadriMesh);
NumerDivEmAvg = sum(cof.*CDWDB,2);

return;