function [p,enew,tnew]=halfmesh(p,e,t)
% 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]=halfmesh(p,e,t); figure;pdeplot(p,e,t)
%
n=size(p,2);
T=size(t,2);
E=size(e,2);
for k=1:T
  p1=p(:,t(1,k));
  p2=p(:,t(2,k));
  p3=p(:,t(3,k));
  p12=(p1+p2)/2;
  p23=(p2+p3)/2;
  p31=(p3+p1)/2;
  i12=find(p(1,:)==p12(1) & p(2,:)==p12(2));
  if isempty(i12)
    n=n+1;
    i12=n;
    p(:,i12)=p12;
  end;
  i23=find(p(1,:)==p23(1) & p(2,:)==p23(2));
  if isempty(i23)
    n=n+1;
    i23=n;
    p(:,i23)=p23;
  end;
  i31=find(p(1,:)==p31(1) & p(2,:)==p31(2));
  if isempty(i31)
    n=n+1;
    i31=n;
    p(:,i31)=p31;
  end;

% Each triangle is divided into four new triangles.
% (p1,p12,p31)
  tnew(1,4*k-3)=t(1,k);
  tnew(2,4*k-3)=i12;
  tnew(3,4*k-3)=i31;
  tnew(4,4*k-3)=t(4,k);
% (p12,p2,p23)
  tnew(1,4*k-2)=i12;
  tnew(2,4*k-2)=t(2,k);
  tnew(3,4*k-2)=i23;
  tnew(4,4*k-2)=t(4,k);
% (p31,p23,p3)
  tnew(1,4*k-1)=i31;
  tnew(2,4*k-1)=i23;
  tnew(3,4*k-1)=t(3,k);
  tnew(4,4*k-1)=t(4,k);
% (p12,p23,p31)
  tnew(1,4*k)=i12;
  tnew(2,4*k)=i23;
  tnew(3,4*k)=i31;
  tnew(4,4*k)=t(4,k);
end;

% Update the edge matrix
for k=1:E
  i1=e(1,k);
  i2=e(2,k);
  x12=(p(1,i1)+p(1,i2))/2;
  y12=(p(2,i1)+p(2,i2))/2;
  i12=1;
  while ~(p(1,i12)==x12 & p(2,i12)==y12)
    i12=i12+1;
  end;
  enew(1:2,2*k-1)=[i1; i12];
  enew(1:2,2*k)=[i12; i2];
  par12=(e(3,k)+e(4,k))/2;
  enew(3:4,2*k-1)=[e(3,k); par12];
  enew(3:4,2*k)=[par12; e(4,k)];
  enew(5:7,2*k-1)=e(5:7,k);
  enew(5:7,2*k)=e(5:7,k);
end;
