# Initial conditions for Navier--Stokes equations

Hi there, I am starting with the test case on 2D incompressible flow around the cylinder. I went through the tutorial many times but did not figure out where the initial condtions are set.

To be more specific, I would like to do a real simulation of an uniform flow entering the rectangular domain by the inlet boundary. I have already modified the inlet boundary conditions. As for the initial conditions, (v,p) should be zero in the whole domain, and (u) is preferred to be initialized as u_in on the inlet boundary and zero elsewhere. Is that implemented in the tutorial? Otherwise is it possible to set up such an initial condition for (u) ?

As I am not sure what version of FEniCS you are using (legacy FEniCS or DOLFINx), and how you have implemented the flow simulation (direct coupled or via a splitting scheme) I will not provide specific code.

You can use `bc.apply(u.vector())` in the case of legacy fenics, and `dolfinx.fem.petsc.set_bc(u.vector, [bc])` in the case of DOLFINx to apply a boundary condition to a function.

I am using DOLFINx. petsc.set_bc works well. The simulation seem to diverge rapidly, but it’s an irrelevant topic to this one…Thank you anyway.

When I tried to transform the velocity from function to numpy arrays, I get strange array shapes. Indeed, `u_.x.array.shape` returns (25478,), but `p_.x.array.shape` and `mesh.geometry.x.shape` return (3252, ) and (3252, 3). I recognize 3252 as the number of nodes in my mesh. The returned value about velocity seem to be caused by the `CG2` function space for velocity. Do you know how can I get the numpy array of velocity fields ?

What do you mean by the numpy array of the velocity field?
When you use a CG-2 function space, you have 6 degrees of freedom (see: DefElement)

Thus it will not match the number of vertices of your mesh.
Also note that each component (x,y,z) of the velocity is a separate degree of freedom.
The array `u_.x.array` matches the coordinates of (`x=V_space.tabulate_dof_coordinates()`), where the i-th row of `x` corresponds to the `u_.x.array[2*i:2*(i+1)]` entries of `u_` (in the 2D setting).

I try to obtain the velocity values on mesh nodes, like what I get by calling `p_.x.array` for the pressure. `V_space.tabulate_dof_coordinates()` looks like a refined version of the original `mesh.geometry.x`, but I do not understand the transformation in between. The correspondence describe in DefElement doesn’t match the reality, at least the first 10 rows of `V_space.tabulate_dof_coordinates()` and `mesh.geometry.x`.
Doesn’t there exist a clean way to obtain velocity values on the mesh nodes?

The point is that there is no one-to-one correspondence between the dof coordinates and the mesh vertices, as illustrated by the plots on DefElement. A second order space has dofs on the facets and vertices. If you want to have only values at the nodes, you could interpolate your CG-2 solution into CG-1.
See for instance: Mapping of results to mesh nodes, debugging of dolfinx implentation of shell formulation - #6 by dokken and the previous posts in that thread

Ok I see. I will go for the interpolation method.