forAll(rho, celli)
{
	vector radial_direction(0, 0, 1);
	double z_coordinate = mesh.C()[celli] & radial_direction;
	if(z_coordinate < 3.6*0.001)
	{
		U[celli].component(0) = 0.0;
		U[celli].component(1) = 49.6;
		U[celli].component(2) = 0.0;
	}
	else if(z_coordinate >= 3.6*0.001 && z_coordinate < 9.1*0.001)
	{
                U[celli].component(0) = 0.0;
                U[celli].component(1) = 11.4;
                U[celli].component(2) = 0.0;
	}
	else if(z_coordinate >= 9.1*0.001 && z_coordinate < 150*0.001)
	{
                U[celli].component(0) = 0.0;
                U[celli].component(1) = 0.9;
                U[celli].component(2) = 0.0;
	}
	else
	{
	}
}