function[K]=computeStiffness(Mesh, Data)

% computeStiffness: Function which assembles the global stiffness matrix
% Assumes all elements are of the same type
% INPUT:
% Mesh:     Structure containing all the mesh parameters
% Data:     Structure containing all the data parameters
% OUTPUT:
% K:        Global stiffness matrix

% Mesh parameters
% Local number of degrees of freedom (dimension in solid mechanics)
d=Mesh.nb_dof_local;
% Total number of degrees of freedom
nb_dof=Mesh.nb_dof;
% Coordinate matrix
coord=Mesh.coord_mat;
% Connectivity matrix
connect=Mesh.connect_mat;
% Finite element order
order=Mesh.order;
% Number of nodes per element
nne=Mesh.nb_nodes_elem;
% Total number of elements
ne=Mesh.nb_elem;
% Matrix linking an element to the degrees of freedom associated to it.
Eq=Mesh.Eq;

% Data parameters
% Mu
f1=Data.Coefficient.mu.function;
% Lambda
f2=Data.Coefficient.lambda.function;


% Retrieve functions
[Functions]=getFunctions(order);
% Jacobian matrix
J_phi=Functions.J_phi;
% Define permutation matrix
Sdd=zeros(d*d);
Id=eye(d);
for k=1:d
    Sdd=Sdd+kron(Id(:,k)', kron(Id, Id(:,k)));
end

% Retrieve quadrature nodes and weights
switch order
    case 1
        [points, weights] = getQuadrature(1, 'bulk');
    case 2
        [points, weights] = getQuadrature(2, 'bulk');
    case 3
        [points, weights] = getQuadrature(3, 'bulk');
    otherwise
        error('Not yet implemented');
end

% Initialization
s=nne*d;
nq=length(weights);
Q=zeros(s^2, nq*d^4);

% Loop over the Gauss points
for i = 1:nq
    L=kron(J_phi(points(i,:)), Id);
    Q(:, (i-1)*d^4+1:i*d^4)=kron(L, L);
end

T=Eq';
G=T(:);

I=repmat(Eq, [1 nne*d])';
J=repmat(G', [nne*d 1]);

T=connect';
G=T(:);

% Coordinates of the nodes of the elements
coord_nodes=coord(G,:);
% 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,:);

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

V=[(vecnorm(c-a,2,2).^2) -sum((c-a).*(b-a),2) -sum((c-a).*(b-a),2) (vecnorm(b-a,2,2).^2)];
W=1./(detJK').^2.*(V');

V=reshape(W,d,d*ne);

JK_inv=zeros(d, d*ne);
JK_inv(:,1:d:d*ne)=1./detJK'.*[c(:,2)-a(:,2) -(b(:,2)-a(:,2))]';
JK_inv(:,2:d:d*ne)=1./detJK'.*[-(c(:,1)-a(:,1)) b(:,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, d^2, ne);

[S1]=product(V , repmat(Id, 1, ne), d, ne);
[S2]=product(JK_invT , JK_inv, d, ne);
S2=Sdd*S2;

A1=reshape(S1, d^4, ne);
A2=reshape(S2, d^4, ne);
A3=kr(vec_JK_invT, vec_JK_invT);


Lambda1=(abs(detJK).*f1(X).*weights)';
Lambda2=(abs(detJK).*f2(X).*weights)';

X=kr(Lambda1, A1+A2)+kr(Lambda2, A3);
K=Q*X;

% Assemblage
K=sparse(I(:), J(:), K(:), nb_dof, nb_dof);

    function[M]=product(A,B,d,ne)
        
        % Blockwise Kronecker product of two matrices
        % !! 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