function grid = refine_grid(coarse_grid)
% Refines a given grid uniformly
%
% Usage:
%    grid = refine_grid(coarse_grid)
%
% Input:
%    coarse_grid   - the coarse grid
%
% Output:
%    grid          - the refined grid


n_nodes = size(coarse_grid.nodes, 2);
n_edges = size(coarse_grid.edges, 2);
n_elements = size(coarse_grid.elements, 2);

nodes = [coarse_grid.nodes zeros(size(coarse_grid.nodes,1), n_edges)];
elements = zeros(3, n_elements*4);

for k=1:n_edges
	e = coarse_grid.edges(:,k);
	nodes(1:2,n_nodes+k) = (nodes(1:2,e(1)) + nodes(1:2,e(2)))/2;
end


for k = 1:n_elements
	t = sort_indices(coarse_grid.elements(:, k), coarse_grid.edges);
	tn = t(1:3);
	te = t(4:6);

	new_n = t(4:6) + n_nodes;
	new_e = n_edges*2 + n_elements*(0:3) + k;

	% new node new_n(i) is on old edge(i)
	% hence the new opposite edge must be the
	% half of an old edge(j) with i!=j
	new_el = (k-1)*4 + (1:4);
	elements(1:3,new_el(1)) =[t(1); new_n(2); new_n(3)];

	elements(1:3,new_el(2)) =[t(2); new_n(1); new_n(3)];

	elements(1:3,new_el(3)) =[t(3); new_n(1); new_n(2)];

	elements(1:3,new_el(4)) =[new_n(1); new_n(2); new_n(3)];
end

grid = generate_grid(nodes, elements);

end




function t = sort_indices(t,edges)
	t = [sort(t(1:3)); t(4:6)];
	
	for i=1:2
		for j=(i+1):3
			if ((edges(1,t(3+i))==t(i)) | (edges(2,t(3+i))==t(i)))
				dummy = t(3+i);
				t(3+i) = t(3+j);
				t(3+j) = dummy;
			end
		end
	end
end

