function [Nr,Tr,Pr]=RefineMesh(N,T,P);
% REFINEMESH refines the mesh by a factor of 4
%   [Nr,Tr,Pr]=RefineMesh(N,T,P); refines the mesh given by the nodes N and
%   the triangles T and physical property P on each triangle by cutting each 
%   triangle into four smaller ones. The boundary is traced so that the new 
%   triangles still contain the information that they are touching the 
%   boundary.

Nr=N;                            % new node list starts with old one
nn=size(N,1);
Tr=zeros(size(T,1)*4,6);         % triangles start from scratch
nt=0;
if nargin==3
  Pr=zeros(1,length(P)*4);
else 
  Pr=[];
end;
NewNid=sparse(nn,nn);
Nid=[1 2;1 3;2 3];
for j=1:size(T,1),               % quadrisect all old triangles
  i=T(j,1:3); n=N(i,:);          % old nodes of current triangle
  for k=1:3
    i1=Nid(k,1);i2=Nid(k,2);
    if i(i1)>i(i2), tmp=i1; i1=i2; i2=tmp; end;   % to only store once  
    if NewNid(i(i1),i(i2))==0
      nn=nn+1; Nr(nn,:)=(n(i1,:)+n(i2,:))/2;
      i(k+3)=nn; NewNid(i(i1),i(i2))=nn;          % store node
    else
      i(k+3)=NewNid(i(i1),i(i2));
    end  
  end
  Tr(nt+1,:)=[i(1) i(4) i(5) T(j,4) 0 T(j,6)];   % insert 4 new triangles
  Tr(nt+2,:)=[i(5) i(4) i(6) 0 0 0];  
  Tr(nt+3,:)=[i(6) i(4) i(2) 0 T(j,4) T(j,5)];  
  Tr(nt+4,:)=[i(6) i(3) i(5) T(j,5) T(j,6) 0];  
  if nargin==3
    Pr(nt+1:nt+4)=P(j);
  end  
  nt=nt+4;
end;