function[Mesh, Solution]=postProcessDisp(Mesh, Data, Solution)

% postProcessDisp: Function used to compute strains and stresses once the 
% displacements are known.
% INPUT:
% Mesh:     Structure containing all mesh parameters
% Data:     Structure containing all data parameters
% Solution: Structure containing the computed solutions
% OUTPUT:
% Mesh:     Updated mesh structure with evaluation points coordinates
% Solution: Updated solution structure including strains and stresses


% Finite element order
order=Mesh.order;
% Number of local degrees of freedom
d=Mesh.nb_dof_local;
% Number of elements
ne=Mesh.nb_elem;
% Number of nodes per element
nne=Mesh.nb_nodes_elem;
% Coordinate matrix
coord=Mesh.coord_mat;
% Connectivity matrix
connect=Mesh.connect_mat;
% Retrieve basis functions
[Functions]=getFunctions(order);
% Retrieve Jacobian matrix
J_phi=Functions.J_phi;
% Number of evaluation points per element
nq=Data.PostProcessing.eval_points;
% Retrieve evaluation points
[points, ~] = getQuadrature(nq, 'interpolation');
% Number of independent strain and stress components
Mesh.nb_comp_ind=1/2*d*(d+1);
% Mu
f1=Data.Coefficient.mu.function;
% Lambda
f2=Data.Coefficient.lambda.function;
% Retrieve solution
U_sort=Solution.U_sort;
% Compute stresses and strains using current position
coord=coord+U_sort;

% Define permutation matrix
Sdd=zeros(d*d);
Id1=eye(d);
Id2=eye(d^2);

for k=1:d
    Sdd=Sdd+kron(Id1(:,k)', kron(Id1, Id1(:,k)));
end

% Reduce size of matrices to avoid unnecessary computations
% sigma_xx, sigma_xy, sigma_yy
h=[1 2 4];

I1=(Id2+Sdd);
I2=Id1(:);

% Truncate the matrices
I1=I1(h,:);
I2=I2(h);

% Initialization
Q=zeros(nq*d^2, nne*d);

% Loop over the Gauss points
for i = 1:nq
    L=kron(J_phi(points(i,:))', Id1);
    Q((i-1)*d^2+1:i*d^2, :)=L;
end
Q=sparse(Q);
T=connect';
G=T(:);

% Coordinates of the nodes of the elements
coord_nodes=coord(G,:);
% Displacements of the nodes of the elements
W=reshape(U_sort(G,:)', nne*d, ne);
% Vertices of the elements
a=coord_nodes(1:nne:end,:);
b=coord_nodes(2:nne:end,:);
c=coord_nodes(3:nne:end,:);
% Interpolated coordinates
P=a(:)+[b(:)-a(:) c(:)-a(:)]*points';
X(:,:,1)=P(1:ne,:);
X(:,:,2)=P(ne+1:end,:);
coord_inter=reshape(P', nq*ne, d);
% Update mesh
Mesh.coord_inter=coord_inter;

% Computation of geometric quantities
% Determinants of Jacobian matrices
detJK=(b(:,1)-a(:,1)).*(c(:,2)-a(:,2))-(b(:,2)-a(:,2)).*(c(:,1)-a(:,1));

JK_invT=zeros(d, d*ne);
JK_invT(:,1:d:d*ne)=1./detJK'.*[c(:,2)-a(:,2) -(c(:,1)-a(:,1))]';
JK_invT(:,2:d:d*ne)=1./detJK'.*[-(b(:,2)-a(:,2)) b(:,1)-a(:,1)]';

vec_JK_invT=reshape(JK_invT, 1, d^2*ne);

% Coefficient evaluation
F1=f1(X);
F2=f2(X);
if size(F1,1)>1
F1=reshape(F1, 1, nq, ne);
end
if size(F2,1)>1
F2=reshape(F2, 1, nq, ne);
end

% Matrix product computations
W=reshape(Q*W, d^2, nq, ne);

B1=reshape(product(JK_invT,repmat(Id1, 1, ne),d,ne), d^2, d^2, ne);
B2=reshape(vec_JK_invT, 1, d^2, ne);

% Computation of strains
if any(contains(Data.PostProcessing.quantity, 'strain'))
    E=pagemtimes(B1, W);
    E=0.5*I1*reshape(E, d^2, nq*ne);
    Solution.epsilon=E;
end

% Computation of stresses
if any(contains(Data.PostProcessing.quantity, 'stress'))
    S1=pagemtimes(B1, F1.*W);
    S2=pagemtimes(B2, F2.*W);
    S1=I1*reshape(S1, d^2, nq*ne);
    S2=I2*reshape(S2, 1, nq*ne);
    S=S1+S2;
    Solution.sigma=S;
end


    function[M]=product(A,B,d,ne)
        
        % !! Only for a two-dimensional case !!
        M=zeros(d^2, d^2*ne);
        
        M(:,1:d^2:end)=kr(A(:,1:d:end), B(:,1:d:end));
        M(:,2:d^2:end)=kr(A(:,1:d:end), B(:,2:d:end));
        M(:,3:d^2:end)=kr(A(:,2:d:end), B(:,1:d:end));
        M(:,4:d^2:end)=kr(A(:,2:d:end), B(:,2:d:end));
    end
end