I am very new to FEM and FeniCS. I need to learn how to use FeniCS to solve NS equations for my work. I am following the tutorial on Naveir Stokes from here https://fenicsproject.org/pub/tutorial/html/._ftut1009.html and have a few questions about the weak form and the code. For the first step of IPCS (eq below), the velocity used for the tensor \sigma is the midpoint velocity \frac{1}{2} (u^n + u^{n+1}):
.
but in the code provided, I see this as the equation for solving the tentative u^*
F1 = rho*dot((u  u_n) / k, v)*dx + \
rho*dot(dot(u_n, nabla_grad(u_n)), v)*dx \
+ inner(sigma(U, p_n), epsilon(v))*dx \
+ dot(p_n*n, v)*ds  dot(mu*nabla_grad(U)*n, v)*ds \
 dot(f, v)*dx
a1 = lhs(F1)
L1 = rhs(F1)
where U is defined as
U = 0.5*(u_n + u)
 Why do I think this means that the U to evaluate \sigma is \frac{1}{2} (u^n + u^{*}) instead of the midpoint velocity?
 Is this step 1 a nonlinear problem as the term \sigma depends on the implicit u? why this first step can be solved as by a linear solve calling
solve()
:
b1 = assemble(L1)
[bc.apply(b1) for bc in bcu]
solve(A1, u_.vector(), b1)

Where can I find an example of using a fully implicit formulation to solve NS for a channel flow?

I want to apply a neumann BC \nabla\vec{u} \cdot \vec{n}p \vec{n} = \vec{g} as an outlet BC. how do I modify the surface integral/boundary terms in the weak form given this midpoint velocity and pressure at time n?
Thanks a lot. Truly appreciate your help.