clc
clear

%Solving not linear task to destroy contact at using finite element method.
%In program uses the rectangular multiplex element.

%
%
%Initial data (input data)

    %Parameters material contactive parts
        %сталь 20л/ steel 430A
        mu = 0.25;        %(0,24..0,28) Poisson's coeff.
        sigma_critical = 216e4; %sigma_0.2
        E = 200e9;        %Young's modulus (for this task elastic modulus)

    %Parameters contact
        Cs = 0.01*E;
        Cn = 0.01*E;

    %Body geometry
        L = 1;    %length
        H = 1/2;    %height parts body
        t = 1/10;

    %Grid options
        N = 11;   %points at length
        M = 11;   %points at height part body




%
%
%Program body

    %Create the initial state and the grid

            [point_global_number ,point_coordinate_in_global_numbers] = create_grid(N ,M ,L ,H);

        %Create the vector of tension
            tension = zeros(2*(N-1)*(M-1) ,3);

        %Create the indication overloaded elements
            overloaded_elements = zeros(2*(N-1)*(M-1) ,1);   %for control elements

        %Create the indication overloaded point on contact
            overloaded_contact_elements = zeros((N - 1) ,2);   %for overload check at basic elements

    %Parametrs finite element

        %Basis elements

        %Matrix elastic constants
            D = [        1, mu/(1-mu),                  0;
                 mu/(1-mu),         1,                  0;
                         0,         0, (1-2*mu)/(2*(1-mu)) ] *E*(1-mu)/((1+mu)*(1-2*mu));

    %Create the stress

        %local force

            point_stress = [203 220 201  232 221 202 240 241 233 234 237];      %the global number of the point
            local_forse = [-5e6 ,0; -5e7 ,0; -5e7 ,0; 5e7 ,0; 5e7 ,0; 5e6 ,0; 0 ,-5e7; 0 ,-5e7; 0 ,-5e7; 0 ,-5e7;  0 ,5e5;];     %[f_x ,f_y]


    %The task is solving at the initial data

%     [delta_U ,new_point_coordinate_in_global_numbers ,K] = FEM_for_linear_tasks(D ,N ,M ,t ,Cn ,Cs ,overloaded_elements ,overloaded_contact_elements ,point_global_number  ,point_coordinate_in_global_numbers ,point_stress ,local_forse);

    %Create parametrs for plots
    FaceColor = zeros(size(point_global_number,1),3);
    FaceColor(:,1) = 0.98824; FaceColor(:,2) = 0.74902; FaceColor(:,3) = 0.10588;
    EdgeColor = zeros(size(point_coordinate_in_global_numbers,1),3);
    EdgeColor(:,1) = 0.0; EdgeColor(:,2) = 0.0; EdgeColor(:,3) = 0.0;

    count = 1;

    if ~exist('plots/', 'dir')
        mkdir('plots/');
    end
    if ~exist('plots/grid/', 'dir')
        mkdir('plots/grid/');
    end
    if ~exist('plots/strains/', 'dir')
        mkdir('plots/strains/');
    end

    old_U = 0;
    while (true)
    if (count == 1)
        fig = domainPlot( point_global_number, point_coordinate_in_global_numbers, FaceColor, EdgeColor,'Iteration 0', false );
        saveas(fig, 'plots/grid/iteration0.png');

    end
    [delta_U ,new_point_coordinate_in_global_numbers ,K] = FEM_for_linear_tasks(D ,N ,M ,t ,Cn ,Cs ,overloaded_elements ,overloaded_contact_elements ,point_global_number  ,point_coordinate_in_global_numbers ,point_stress ,local_forse);
    [new_tension ,new_overloaded_elements ,new_overloaded_contact_elements ,virtual_tension,a ,b] = tension_in_element (M ,N ,Cs ,Cn ,sigma_critical ,D ,tension ,overloaded_elements ,point_global_number ,point_coordinate_in_global_numbers ,overloaded_contact_elements ,delta_U);

    if (new_overloaded_contact_elements == overloaded_contact_elements) | (new_overloaded_contact_elements == overloaded_contact_elements)
        break;
    else
        fig = domainPlot( point_global_number, new_point_coordinate_in_global_numbers, FaceColor, EdgeColor, strcat('Iteration',{' '},int2str(count)), false );
        saveas(fig, sprintf('plots/grid/iteration%d.png',count));

        for i = 1:1:3
            switch i
                case 1
                    title = sprintf('\\sigma_x %d',count);
                case 2
                    title = sprintf('\\sigma_y %d',count);
                otherwise
                    title = sprintf('\\tau_{xy} %d',count);
            end
            fig = domainPlot( point_global_number, new_point_coordinate_in_global_numbers, abs(new_tension(:,i)), EdgeColor, title, true );
            saveas(fig, sprintf('plots/strains/%s.png',regexprep(title,'[\\{}_ ]','')));

        end
        overloaded_elements = new_overloaded_elements;
        overloaded_contact_elements = new_overloaded_contact_elements;
        old_U = delta_U;
    end
    count = count + 1;
end