function average=averaging(md,data,iterations)
%AVERAGING - smooths the input over the mesh
%
%   This routine takes a list over the elements or the grids in input
%   and return a list over the grids.
%   For each iterations it computes the average over each element (average 
%   of the vertices values) and then computes the average over each grid
%   by taking the average of the element around a grid weighted by the
%   elements volume
%
%   Usage:
%      md=averaging(md,data,iterations)
%
%   Examples:
%      velsmoothed=averaging(md,md.vel,4);
%      pressure=averaging(md,md.pressure,0);

if length(data)~=md.numberofelements & length(data)~=md.numberofgrids
	error('averaging error message: data not supported yet');
end

%initialization
weights=zeros(md.numberofgrids,1);
areas=zeros(md.numberofelements,1);
data=data(:);

%load some variables (it is much faster if the variab;es are loaded from md once for all)
index=md.elements;
numberofgrids=md.numberofgrids;
numberofelements=md.numberofelements;

%build some variables
line=index(:);
if strcmpi(md.type,'3d')
	rep=6;
else
	rep=3;
end
summation=1/rep*ones(rep,1);
linesize=rep*numberofelements;

%compute the volume of each element
areas=area(md);

%update weights that holds the volume of all the element holding the grid i
weights=sparse(line,ones(linesize,1),repmat(areas,rep,1),numberofgrids,1);

%initialization
if length(data)==numberofelements
	average_grid=sparse(line,ones(linesize,1),repmat(areas.*data,rep,1),numberofgrids,1);
	average_grid=average_grid./weights;
else
	average_grid=data;
end

%loop over iteration
for i=1:iterations
	average_el=average_grid(index)*summation;
	average_grid=sparse(line,ones(linesize,1),repmat(areas.*average_el,rep,1),numberofgrids,1);
	average_grid=average_grid./weights;
end

%return output as a full matrix (C code do not like sparse matrices)
average=full(average_grid);
