% h-adaptivity
% for a given mesh and solution we generate a new adapted mesh by splitting
% the elements with large errors into two new elements

function mesh = adapt_h(nodes,u,lambda)
    n_nodes = length(nodes);
    flux_jumps = jumps_flux(nodes,u,lambda);
    alpha = 0.5;    % parameter for choosing elements to refine
    crit_error = alpha*max(abs(flux_jumps));
    mesh = nodes;   % will be the new nodes/mesh
    % finding elements to refine
    for i=2:n_nodes-1
        if abs(flux_jumps(i)) > crit_error
            mesh = [mesh, nodes(i-1)+(nodes(i)-nodes(i-1))/2, nodes(i)+(nodes(i+1)-nodes(i))/2];
        end
    end
    mesh = unique(mesh);
end