function N=SmoothMesh(N,T,bn);
% SMOOTHMESH smooth mesh closer to improve mesh quality
%   Nr=SmootheMesh(N,T,bn); modify node positions iteratively to get
%   closer to a Delaunay triangulation, by iteratively assigning to
%   nodes the averages of neighboring nodes. Node numbers in bn are
%   not allowed to have their position be modified

nn=size(N,1);
Nn=zeros(nn,1);
nt=size(T,1);                    % number of triangles
if nargin<3
  bn=[];
end
for i=1:nt                       % first find boundary nodes
  for j=1:3
    if T(i,3+j)==1 | T(i,mod(3+j+1,3)+4)==1 
      if isempty(find(bn==T(i,j)))
        bn=[bn T(i,j)];
      end
    end    
  end    
end    
for i=1:nt                       % construct neigboring node list
  for j=1:3
    if isempty(find(bn==T(i,j))) % not a boundary node
      for k=1:2                  % add neigboring nodes from triangle
        id=find(Nn(T(i,j),2:Nn(T(i,j),1)+1)==T(i,mod(j+k-1,3)+1));
        if isempty(id)
          Nn(T(i,j),1)=Nn(T(i,j),1)+1;
          Nn(T(i,j),Nn(T(i,j),1)+1)=T(i,mod(j+k-1,3)+1);
        end
      end
    end
  end
end  
for i=1:10
  for j=1:nn
    if Nn(j,1)>0
      N(j,:)=mean(N(Nn(j,2:Nn(j,1)+1),:));
    end
  end
end