%[new_tension ,new_overloaded_elements ,new_overloaded_contact_elements] = tension_in_element (M ,N ,Cs ,Cn ,sigma_critical ,D ,tension ,overloaded_elements ,point_global_number_sum ,coordinate_sum ,overloaded_contact_elements ,delta_U)

function [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_sum ,coordinate_sum ,overloaded_contact_elements ,delta_U)

        % Create the vector of tension
            new_tension = zeros(size(tension ,1) ,size(tension ,2));

        % Create the local vector of tension

            virtual_tension = zeros(length(point_global_number_sum),3);
            for i = 1:length(point_global_number_sum)

                % Overload check at element
                    % 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
                % Tension determ
                    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;

                    U = zeros(8,1);
                    for k = 1:1:4
                        for l = 1:1:2
                            U(2*k-mod(l,2)) = delta_U(2*point_global_number_sum(i,k)-mod(l,2));
                        end
                    end
                    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);
                            virtual_tension_local = D_local*B*U;
                        end
                    end

                    virtual_tension(i ,:) = virtual_tension_local(:);
            end

            for i = 1:length(tension)
                if i <= length(tension)/2
                    new_tension(i ,:) =  virtual_tension(i ,:);
                else
                    k = i + ( length(point_global_number_sum) - length(tension) );
                    new_tension(i ,:) =  virtual_tension(k ,:);
                end
            end

        % Create new vector for overload check at element
            new_overloaded_elements = overloaded_elements;

            for i = 1:length(overloaded_elements)
                if abs(new_tension(i ,:)) >= abs(sigma_critical)
                    new_overloaded_elements(i) = 1;
                end
            end

        % Create new vector for overload check at element
            new_overloaded_contact_elements = overloaded_contact_elements;

            for i = 1:2*( N -1 )
                %index for contacts elemetns
                if i > ( N -1 )
                    j = i - ( N -1 );
                else
                    j = i;
                end

                if new_tension((N - 1)*(M - 1) - (N - 1) + i ,2) >= Cn || new_overloaded_contact_elements(j ,2) > 0

                    new_overloaded_contact_elements(j ,2) = 1;
                elseif abs(new_tension((N - 1)*(M - 1) - (N - 1) + i ,1)) >= Cs || new_overloaded_contact_elements(j ,1) > 0

                    new_overloaded_contact_elements(j ,1) = 1;
                end
            end
%         % Search the point of gap
%             gap_points = new_overloaded_contact_elements(: ,2);
%             for i = 1:lenght(new_overloaded_contact_elements)
%                 if i == 1 || i == lenght(new_overloaded_contact_elements)
%                     if gap_points(i) > 1
%                         new_overloaded_contact_elements(i ,2) = 1;
%                     end
%                 elseif gap_points(i) > 0 && gap_points(i - 1) 0
%             end

end