function md=Loadresults(md,fem,u_g);
%LOADRESULTS - load results onto model
%
%   The solution (vel,pressure,...) are on the g-set,
%   use index to get the value on each node
%
%   Usage:
%      md=Loadresults(md,fem,u_g)

if strcmpi(md.type,'2d')
	if fem.ismacayealpattyn,
		m_dh=fem.m_dh;
		indx=1:6:m_dh.gridset.gsize;
		indy=2:6:m_dh.gridset.gsize;
		indp=4:6:m_dh.gridset.gsize;
		indx=indx(m_dh.tpart);
		indy=indy(m_dh.tpart);
		indp=indp(m_dh.tpart);
	else
		m_dhu=fem.m_dhu;
		indx=1:6:m_dhu.gridset.gsize;
		indy=2:6:m_dhu.gridset.gsize;
		indp=4:6:m_dhu.gridset.gsize;
		indx=indx(m_dhu.tpart);
		indy=indy(m_dhu.tpart);
		indp=indp(m_dhu.tpart);
	end

	md.vx=u_g(indx)*md.yts;
	md.vy=u_g(indy)*md.yts;
	md.vel=sqrt(md.vx.^2+md.vy.^2);
	md.pressure=u_g(indp);

else
	m_dv=fem.m_dv;
	indx=1:6:m_dv.gridset.gsize;
	indy=2:6:m_dv.gridset.gsize;
	indz=3:6:m_dv.gridset.gsize;
	indp=4:6:m_dv.gridset.gsize;
	indx=indx(m_dv.tpart);
	indy=indy(m_dv.tpart);
	indz=indz(m_dv.tpart);
	indp=indp(m_dv.tpart);

	md.vx=u_g(indx)*md.yts;
	md.vy=u_g(indy)*md.yts;
	md.vz=u_g(indz)*md.yts;
	md.vel=sqrt(md.vx.^2+md.vy.^2+md.vz.^2);
	md.pressure=u_g(indp);
end
