function[Data]=preCompute(Mesh, Data)

% preCompute: Function which pre-computes quantities needed in a single 
% iteration of the contact algorithm
% INPUT:
% Mesh:     Structure containing all the mesh parameters
% Data:     Structure containing all the data parameters
% OUTPUT:
% Data:     Updated data structure

d=Mesh.nb_dof_local;

% Retrieve bodies
body=Data.body;
% Retrieve local degrees of freedom
indx1=body{1}.indx_local;
indx2=body{2}.indx_local;
l1=length(indx1);
l2=length(indx2);

Data.l1=l1;
Data.l2=l2;

b2=Data.b2;

L22=Data.L22;

i2_old=Data.i2;
i2=[indx1; indx2];
% Beware of the ordering in the linear operators
[~, index, reor]=intersect(i2_old, i2, 'stable');
Data.index=index;

%% Modified LU

Phi11=Data.Phi11;
Phi21=Data.Phi21;
Phi12=Data.Phi12;
Phi22=Data.Phi22;

r1=size(Phi11,1);
r2=size(Phi22,1);

q1=ones(r1,1);
q2=ones(r2,1);

s1=Phi12*(Phi22\q2);
s2=Phi21*(Phi11\q1);

% Define initial sparsity
S1=num2cell(1:r1);
S2=num2cell(1:r2);
% Maximum fill in per column of M per iteration
Param.s=5;
Param.maxit=20;
Param.tol=0.2;

[Phi11ia]=SPAI(Phi11, S1, Param);
[Phi22ia]=SPAI(Phi22, S2, Param);

R12a=Phi12*Phi22ia;
R21a=Phi21*Phi11ia;

% Rescaling
R12a=R12a./s1;
R21a=R21a./s2;

M1=kron(body{1}.interface_mass, speye(d));
M2=kron(body{2}.interface_mass, speye(d));

R12a=kron(R12a, speye(d));
R21a=kron(R21a, speye(d));

Q1a=-diag(M2).*R21a./(diag(M1)');
Q2a=-R12a';

Y1a=[eye(l1); Q1a];
Y2a=[eye(l1); Q2a];

Z1=sparse(b2, l1);
Z2=sparse(b2, l1);

% Z1(index,:)=-alpha*Y1a(reor,:);
% Z2(index,:)=Y2a(reor,:);

Z1(index,:)=Y1a(reor,:);
Z2(index,:)=Y2a(reor,:);

Lv=L22;
Uv=L22;

[i,j]=find(Z1);
[I]=findNNZc(i,j);

Lv(:,I)=Z1;
Uv(:,I)=Z2;

Data.Lv=Lv;
Data.Uv=Uv';
Data.orl=1:b2;
Data.orc=1:b2;

end