%[delta_U ,new_coordinate_sum ,K] = FEM_for_linear_tasks(D ,N ,M ,t ,Cn ,Cs ,overloaded_elements ,overloaded_contact_elements ,point_global_number_sum  ,coordinate_sum ,point_stress ,local_forse)

function [delta_U ,new_coordinate_sum ,K] = FEM_for_linear_tasks(D ,N ,M ,t ,Cn ,Cs ,overloaded_elements ,overloaded_contact_elements ,point_global_number_sum  ,coordinate_sum ,point_stress ,local_forse)

        % Create stiffness global matrix

            K = zeros(2*2*N*M);

        % Create stiffness local matrix
            for i = 1:(N-1)*(2*M-1)%length(overloaded_elements)

                if i > (N-1)*(M-1) && i <= M*(N-1)

                    % Overload check at contact element
                    j = i - (N-1)*(M-1);
                    if abs(overloaded_contact_elements(j ,2)) > 0
                        C_s = 0;
                        C_n = 0;
                    elseif abs(overloaded_contact_elements(j ,1)) > 0
                        C_s = 0;
                        C_n = Cn;
                    else
                        C_s = Cs;
                        C_n = Cn;
                    end

                    % Stiffness matrix (B'*C*B)
                        K_contact_local = 1/3*[ 2*C_s,      0,    C_s,      0,   -C_s,      0, -2*C_s,      0;
                                                    0,  2*C_n,      0,    C_n,      0,   -C_n,      0, -2*C_n;
                                                  C_s,      0,  2*C_s,      0, -2*C_s,      0,   -C_s,      0;
                                                    0,    C_n,      0,  2*C_n,      0, -2*C_n,      0,   -C_n;
                                                 -C_s,      0, -2*C_s,      0,  2*C_s,      0,    C_s,      0;
                                                    0,   -C_n,      0, -2*C_n,      0,  2*C_n,      0,    C_n;
                                               -2*C_s,      0,   -C_s,      0,    C_s,      0,  2*C_s,      0;
                                                    0, -2*C_n,      0,   -C_n,      0,    C_n,      0,  2*C_n];
                    K_local = K_contact_local;
                else

                    % Overload check at element
                    if i <= (N-1)*(M-1) && overloaded_elements(i) == 0
                        D_local = D;
                    elseif i <= (N-1)*(M-1) && overloaded_elements(i) > 0
                        D_local = D/2;
                    elseif i > M*(N-1) && overloaded_elements(i - (N-1)) == 0
                        D_local = D;
                    else
                        D_local = D;
                    end

                    K_local = zeros(8);

                    a = ( abs(coordinate_sum(point_global_number_sum(i,1) ,1) - coordinate_sum(point_global_number_sum(i,2) ,1)) + abs(coordinate_sum(point_global_number_sum(i,3) ,1) - coordinate_sum(point_global_number_sum(i,4) ,1)) )/4;
                    b = ( abs(coordinate_sum(point_global_number_sum(i,1) ,2) - coordinate_sum(point_global_number_sum(i,4) ,2)) + abs(coordinate_sum(point_global_number_sum(i,2) ,2) - coordinate_sum(point_global_number_sum(i,3) ,2)) )/4;

                    for k = 1:1:2
                        xi = 0 + ((-1)^k)*1/sqrt(3);
                        for l = 1:1:2
                            eta = 0 + ((-1)^l)*1/sqrt(3);
                            B = [-b*(1-eta),          0, b*(1-eta),         0, b*(1+eta),         0, -b*(1+eta),          0;
                                          0,  -a*(1-xi),         0, -a*(1+xi),         0,  a*(1+xi),          0,   a*(1-xi);
                                  -a*(1-xi), -b*(1-eta), -a*(1+xi), b*(1-eta),  a*(1+xi), b*(1+eta),   a*(1-xi), -b*(1+eta)]/(4*a*b);
                            K_local = K_local + t*a*b*(4*a*b)*B'*D_local*B;
                        end
                    end
                end

            % load K_local in K_global
                for k = 1:1:8
                    node1 = point_global_number_sum(i,fix((k+1)/2));
                    mod1 = mod(k,2);
                    for l = 1:1:8
                        node2 = point_global_number_sum(i,fix((l+1)/2));
                        mod2 = mod(l,2);
                        K(2*node1-mod1,2*node2-mod2) = K(2*node1-mod1,2*node2-mod2) + K_local(k,l);
                    end
                end
            end

        % Create fixed
%             bound = find(point_global_number_sum(:,2) == 0)
%             for i = 1:1:size(bound,1)
%                 for j = -1:1:0
%                     K(2*bound(i)+j, :) = 0;
%                     K(2*bound(i)+j,2*bound(i)+j) = 1;
%                 end
%             end
            for ni = 1:2*N
                for nj = 1:2*N
                    if ni == nj
                        K(ni ,nj) =1;
                    else
                        K(ni ,nj) =0;
                    end
                end
            end

%             for ni = 1:2*N
%                 K(ni ,:) = [];
%                 K(: ,ni) = [];
%             end

        % Create force

            F = zeros(2*2*N*M,1);
            for i = 1:length(point_stress)
                F(point_stress(i)*2 - 1) = local_forse(i ,1)*t;
                F(point_stress(i)*2) = local_forse(i ,2)*t;
            end

%             for ni = 1:2*N
%                 F(ni) = [];
%             end

        % Solve system

%             short_U = K\F;
%
%             delta_U = zeros(2*M*N,1);
%
%             for ni = 1:length(short_U)
%                 delta_U(ni + 2*N) = short_U(ni);
%             end
%

            delta_U = K\F;


        % New coordinate grid

            new_coordinate_sum = coordinate_sum;
            for i = 1:max(size(coordinate_sum))
                new_coordinate_sum(i ,1) = new_coordinate_sum(i ,1) + delta_U(2*i - 1);
                new_coordinate_sum(i ,2) = new_coordinate_sum(i ,2) + delta_U(2*i);
            end
end