function [meshnodnew,S]=sprngnod(meshele,meshnod,bndele,nrelaxes)
% [meshnodnew,S]=sprngnod(meshele,meshnod,bndele,[nrelaxes]);
% meshtools/sprngnod.m   Like nodes with springs between this will')
% relax all the nodes to centers of their connections')
% nrelaxes is optional, one default,  number of times through relaxation loop
%  Version 3 makes use of sparse matrix multiplications

if nargin==3
nrelaxes = 1;
end

disp(' meshtools/sprngnod.m   Like nodes with springs between this will')
disp('relax all the nodes to centers of their connections')




polys = meshele;    
polys(:,1)=meshele(:,4);   %   triangles close on themselves

[m,n]=size(polys);
c = ones(4*m,1);
c(:) = polys;
l = min(c):max(c);    % this will be harder    
mc = min(c)-1

xy=meshnod(l,2:3);
i1 = polys(:,1)-mc;
j1 = polys(:,2)-mc;
i2 = polys(:,2)-mc;
j2 = polys(:,3)-mc;
i3 = polys(:,3)-mc;
j3 = polys(:,4)-mc;
i4 = polys(:,4)-mc;
j4 = polys(:,1)-mc;
S = sparse([i1;i2;i3],[j1;j2;j3],1);
%spy(S,'g')
Srow= sum(S);
for k = 1:length(Srow)
S(k,:) = S(k,:)./Srow(k);
end



bnds = [ bndele(1,2) ; bndele(:,3) ];
meshnodnew=meshnod;

tic
for k = 1:nrelaxes
% disp(['sprngnod Relaxation number  ',int2str(k),'   ',num2str(toc)])
meshnodnew(:,2) = S*meshnodnew(:,2);
meshnodnew(:,3) = S*meshnodnew(:,3);
meshnodnew(bnds,2) = meshnod(bnds,2);
meshnodnew(bnds,3) = meshnod(bnds,3);
end
 disp(['sprngnod Relaxation number  ',int2str(k),'   ',num2str(toc)])

x3 =  meshele;
x3(:,1) = meshele(:,4);
x3(:) = meshnodnew(x3,2);
y3 =  meshele;
y3(:,1) = meshele(:,4);
y3(:) = meshnodnew(y3,3);
x3 = [ x3 NaN*meshele(:,1)]';
y3 = [ y3 NaN*meshele(:,1)]';
plot(x3,y3,'w')
axis('equal')
title(['Spring nodes Relaxation=',int2str(nrelaxes)])
