% INTERNODES method for contact mechanics
% Case study: Hertzian contact problem with Dirichlet boundary conditions
% applied to the lower body and Neumann boundary conditions applied to the
% upper body.

% Limitations: 2D
clc
clear variables
close all

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

%% Parameter prescription
% Mesh parameters
[Mesh]=getParameters(Mesh);
% Initialize model
[Data]=setModel('model', 'static_solid', '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;

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

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

%% Set boundary conditions

% Used for testing with singular stiffness matrix
BC={{'Type', 'Dirichlet', 'Edges', 3, 'Function', @(x) [0; 0]},...
    {'Type', 'Neumann', 'Edges', 10, 'Function', @(x) [0; -1e9]}};

p2={'BC', BC};    

[Data]=setBoundaryConditions(Mesh, Data, p2);

%% Set numerical parameters
% For more information, type 'help setNumericalParameters'

[Data]=setNumericalParameters(Mesh, Data, 'Contact', 'rescaling', E, 'maxiter', 30, 'radius', inf, 'tol', 0.9);
[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]=solve(Mesh, Data, {3:6, 7:10});

%% 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');