function [p,e,t]=qmesh(p,e,t);
%QMESH Modify a pde triangular mesh to be used with quadratic basis functions
%
%	[P,E,T]=qmesh(P,E,T) returns a triangular mesh with side midpoints
%	appended into the point matrix P. In the new Edge matrix E, the second
%	row contains midpoint indices and rows 3-8 contain rows 2-7 from the
%	old edge matrix. In the triangle matrix T, rows 4-6 contain midpoint
%	indices and row 7 contains the subdomain number (row 4 from the old
%	triangle matrix).
%	Input arguments P, E and T can be created with pdetool

% Heikki Laitinen, May -97
% This function was created by modifying JvP's tofour.m

t(7,:)=t(4,:);  % Make room for midpoint indices
e(3:8,:)=e(2:7,:);

n = size(p,2);
midp = sparse(n,n);  % Midpoint indices; the midpoint of the side p(:,i)-
		     % p(:,j) will be stored as p(:,midp(i,j))

for ti=1:size(t,2); % scan through all triangles
   
%  Bisect sides and calculate new points
  for ppair = [t([2:3,1],ti), t([3,1:2],ti)]'
% Bisect edge, if not already done
    if ~midp(ppair(1),ppair(2))
       n=n+1;
       p(:,n)=p(:,ppair)*[1;1]/2;
       midp(ppair(1),ppair(2)) = n;
       midp(ppair(2),ppair(1)) = n;
    end;
  end;  % for ppair

  t(4:6,ti)=[midp(t(2,ti),t(3,ti)), midp(t(1,ti),t(3,ti)), ...
	     midp(t(1,ti),t(2,ti))]';

end % for ti 

% Update edge matrix
for ei=1:size(e,2) % scan through edges
  e(2,ei)=midp(e(1,ei),e(3,ei));
end;
