function [p,enew,tnew]=nmesh(p,e,t,n)
% Heikki Laitinen March -97
% Example: p = [ 0,0; 1,0; .5,1/sqrt(2)]';t = [1;2;3;1];
%          e = [1,2,3; 2,3,1; 0,0,0; 1,1,1; 1,1,1; 1,1,1; 1,2,3];
%          pdeplot(p,e,t); [p,e,t]=nmesh(p,e,t,8); figure;pdeplot(p,e,t)
N=size(p,2);
T=size(t,2);
E=size(e,2);
tcount=1;
for k=1:T
  p1=p(:,t(1,k));
  p2=p(:,t(2,k));
  p3=p(:,t(3,k));
  dp12=(p2-p1)/n;
  dp13=(p3-p1)/n;
  delta=min([norm(dp12) norm(dp13)])/(10*n);
  for i=1:n+1
    for j=1:n+2-i
      newnode=p1+(j-1)*dp12+(i-1)*dp13;
      index=find(abs(p(1,:)-newnode(1))<delta & abs(p(2,:)-newnode(2))<delta);
      if isempty(index)
        N=N+1;
        indices(i,j)=N;
        p(:,N)=newnode;
      else
        indices(i,j)=index;
      end;
    end;
  end;

  for j=1:n
    tnew(:,tcount)=[indices(1,j); indices(1,j+1); indices(2,j); t(4,k)];
    tcount=tcount+1;
  end;
  for i=2:n
    for j=1:n+1-i
      tnew(:,tcount)=[indices(i,j); indices(i-1,j+1); indices(i,j+1); t(4,k)];
      tcount=tcount+1;
      tnew(:,tcount)=[indices(i,j); indices(i,j+1); indices(i+1,j); t(4,k)];
      tcount=tcount+1;
    end;
  end;
end;

% Reunamatriisi e pitäisi vielä päivittää...
