% INTERNODES method for contact mechanics
% Case study: Hertzian contact problem with Dirichlet boundary conditions
% applied to both bodies. Elastic contact between a sphere and a half
% space.

% Limitations: 2D
clc
clear variables
close all

% Read GMSH mesh file (.msh)
[Mesh] = readGMSH('../../GMSH/Meshes/contact10_p2_0_05.msh');

%% Parameter prescription
% Mesh parameters
[Mesh]=getParameters(Mesh);
% Initialize model
[Data]=setModel('model', 'static_solid', 'submodel', 'contact', 'type', 'linear');

%% Set material parameters
% Elastic modulus [N\m^2]
E=30e9;
% Poisson ratio [-]
nu=0.2;
% Specific mass [kg/m^3]
rho=2500;
% Thickness [m]
t=1;
% Amplification factor for the stiffness of the upper body
f=10;
% Radius of the upper sphere
R=0.5;
% Prescribed displacement
d=0.15;

% Plane stress conditions
lambda=E*nu/(1-nu^2);
mu=E/(2*(1+nu));

[Data]=setCoefficient(Data, 'rho', t*rho);
[Data]=setCoefficient(Data, 'lambda', @(x) t*(lambda*(x(2)<=2)+f*lambda*(x(2)>2)));
[Data]=setCoefficient(Data, 'mu', @(x) t*(mu*(x(2)<=2)+f*mu*(x(2)>2)));
[Data]=setCoefficient(Data, 'f', @(x) [0; 0]);

%% Set boundary conditions

% Used for testing with invertible stiffness matrix
BC={{'Type', 'Dirichlet', 'Edges', 9, 'Function', @(x) [0; 0]},...
    {'Type', 'Dirichlet', 'Edges', 14, 'Function', @(x) [0; -d]}};

p2={'BC', BC};    

[Data]=setBoundaryConditions(Mesh, Data, p2);
%% Set numerical parameters
% For more information, type 'help setNumericalParameters'

[Data]=setNumericalParameters(Mesh, Data, 'Contact', 'rescaling', E, 'radius', 0.2, 'maxiter', 30);
[Data]=setNumericalParameters(Mesh, Data, 'PostProcessing', 'ampl_fact', 1, 'quantity', {'stress', 'strain'});

%% Assemblage of the right-hand side
[Data]=computeRHS(Mesh, Data);

%% Solve the contact problem
[Mesh, Data, Solution, A, b]=solve(Mesh, Data, {9:12, 13:14});

%% Computation of the stresses and strains
[Mesh, Solution]=postProcess(Mesh, Data, Solution);

%% Visualization of the solution
% For more information, type 'help plotSolution'
plotSolution(Mesh, Data, Solution, 'disp', 'standard', 'stress', 'sigma_yy');

%% Comparison between numerical and analytical solutions

E_star=(1/E*(1-nu^2)+1/(10*E)*(1-nu^2))^(-1);
F=4/3*E_star*sqrt(R)*d^(3/2);
% Maximum pressure
p0=1/pi*(6*F*E_star^2/R^2)^(1/3);

lambda=Solution.lambda;
lambda_max=max(abs(lambda));