%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% MATLAB routine of the counterexample presented in: 
%
% "A. Bemporad, M. Paggi: Optimization algorithms for the solution of the frictionless normal contact between rough surfaces"
% submitted to Int. J. Sol. Structures
%
% Release date: 01/11/2014
%
% Authors: Alberto Bemporad and Marco Paggi
%
% Affiliation: IMT Institute for Advanced Studies
%              Piazza San Francesco 19
%              55100 Lucca - Italy
%              emails: alberto.bemporad@imtlucca.it; marco.paggi@imtlucca.it
%              http://musam.imtlucca.it 
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Purpose: show counterexamples where the Greedy algorithm fails in solving the linear complementarity problem (LCP).
%          By repeating the subroutine, the vector \bar_{u} is randomly generated.
%          In some cases (about 40% of the cases for the values of the parameters reported below and almost regardless of n=#I_C)
%          w^* is negative valued in some elements, violating the LCP condition.
%
% Input:   Young modulus, E
%          Grid spacing, delta
%          Number of boundary elements per side, N (square grid)
%          Tolerance of the algorithm, feastol
%          \bar_{u} randomly generated by repeating the execution of the routine
%
% Output:  \bar_{u} vs. boundary element number, Fig.1, displaying the input data of the problem. 
%          w^* vs. boundary element number, Fig.2. If nonnegative, the LCP is solved exactly.
%          p^*/E vs. boundary element number, Fig.3, displaying the contact tractions.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all
clc
clf

N=10;   
E=0.01;
delta=1;
feastol=1e-5;

n=N*N; 

% generation of the matrices containing the barycentric coordinates of the boundary elements

for i=1:N
    for j=1:N
       xm(i,j)=i*delta;
       ym(i,j)=j*delta;
    end
end

xv=reshape(xm,n,1);
yv=reshape(ym,n,1);

radius=delta/2;
C=1/(E*pi*radius);
Q=diag(C*ones(n,1));

for i=1:n
    for j=i:n
        if i~=j,
           r(i,j)=sqrt((xv(j)-xv(i))^2+(yv(j)-yv(i))^2);
           Q(i,j)=C*asin(radius/r(i,j));
           Q(j,i)=Q(i,j); 
        end
    end
end

% random assignment of \bar{u}

ub=rand(n,1);  

% application of the Greedy algorithm in conjunction with the built-in mldivide function of MATLAB

k=0;
kept=true(n,1);
keepgoing=true;

while keepgoing,
    p(kept)=Q(kept,kept)\ub(kept);
    k=k+1;
    f=find(p<-feastol);
    keepgoing=any(f);
    p(f)=0;
    kept(f)=false;
    w=Q*p'-ub;
    minw(k)=min(w); % store of the minimum of w during the iterations to check the correctness of the solution
end

numiter=1:k;

% minimum value of w vs. number of Greedy iterations, k 
% The last value corresponds to w^*, which should be always nonnegative

[numiter',minw'] 

% Figure 1: \bar{u} vs. boundary element number

figure(1)
plot(ub,'k')
xlabel('Boundary element','fontsize',20,'Interpreter','latex')
ylabel('$$\bar{u}$$','fontsize',20,'Interpreter','latex')

% Figure 2: w^* vs. boundary element number

figure(2)
plot(w,'r')
xlabel('Boundary element','fontsize',20,'Interpreter','latex')
ylabel('$$w^*$$','fontsize',20,'Interpreter','latex')

% Figure 3: p/E vs. boundary element number

figure(3)
plot(p/E,'k')
xlabel('Boundary element','fontsize',20,'Interpreter','latex')
ylabel('$$p/E$$','fontsize',20,'Interpreter','latex')

