examples/navier-stokes.md

@@ -142,3 +142,55 @@ where we have used the identity
 ```math
 \nabla\mathbf{u}:\nabla\mathbf{v} = \sum_i \nabla u_i\cdot\nabla v_i
 ```
+
+### External volume force
+The force may be implemented as a vector-valued function of the global coordinates or
+any other expression that leads to a local force vector at the quadrature points, e.g.
+```c++
+using WorldVector = typename Grid::template Codim<0>::Geometry::GlobalCoordinate;
+auto opForce = makeOperator(tag::testvec{},
+  [](WorldVector const& x) { return WorldVector{1.0, 0.0}; });
+prob.addVectorOperator(opForce, _v);
+```
+
+Numerical example
+-----------------
+We consider a square domain $`\Omega`$ with $`L_x=L_y=1`$ and Dirichlet boundary conditions
+```math
+\mathbf{g} = \begin{pmatrix}0 \\ x_1(1 - x_1)(1 - x_0)\end{pmatrix}
+```
+Thus, the boundary value is zero on top, right, and bottom boundary and parabolic
+on the left boundary.
+
+This is the lid-driven cavity problem, where we do not include any other external force
+term.
+
+The boundary condition is implemented, by first setting boundary ids for parts of
+the grid boundary and second, setting the Dirichlet value:
+```c++
+double vel = 1.0;
+Parameters::get("stokes->boundary velocity", vel);
+
+// define boundary values
+auto g = [vel](FieldVector<double,2> const& x)
+{
+  return FieldVector<double,2>{0.0, vel*x[1]*(1.0 - x[1])*(1.0 - x[0])};
+};
+
+// set boundary conditions for velocity
+prob.boundaryManager()->setBoxBoundary({1,1,1,1});
+prob.addDirichletBC(1, _v, _v, g);
+```
+
+### Simulation
+The time-stepping process can stationary problem in each iteration can be started
+using an `AdaptInstationary` wrapper class:
+```c++
+// set initial conditions
+prob.solution(_v).interpolate(g);
+
+// start simulation
+AdaptInfo adaptInfo("adapt");
+AdaptInstationary adapt("adapt", prob, adaptInfo, probInstat, adaptInfo);
+adapt.adapt();
+```
\ No newline at end of file