Given a vector-valued PDE, I'd like to enforce the boundary conditions $$ \vec{n}\cdot u = g\\ \vec{n}\cdot \nabla (\vec{t}\cdot u) = 0 $$ on the solution $\vec{u}$. If the boundary happens to align with one of the coordinate axes, I could use

bcs = DirichletBC(V.sub(0), g, 'on_boundary')

as Dirichlet condition. What to do for other geometries though?


2 Answers 2


You can use a Nitsche-type method for this. See the following reference:

J. Freund, R. Stenberg. On weakly imposed boundary conditions for second order problems. Proceedings of the Ninth Int. Conf. Finite Elements in Fluids, Venice 1995. M. Morandi Cecchi et al., Eds. pp. 327-336.

I have implemented this a while ago in some simple FEniCS-Code to deal with free-slip boundaries in general 2d geometries. Find the demo code here and the mesh here (tested with FEniCS 1.3).

Cylinder bump: free-slip boundary conditions via Nitsche's method

Update: Meanwhile, others have made similar implementations available. See for instance this github repository.

  • $\begingroup$ The link to the demo code is out of date. $\endgroup$ Apr 7, 2018 at 20:37
  • $\begingroup$ @NicoSchlömer: thanks. The file was still up, but Dropbox seems to have forgotten that it was supposed to be shared. Hope pastebin lasts longer. $\endgroup$ Apr 11, 2018 at 3:08

This maybe raises more questions than it answers, but here is what I think one should do.

Lets assume that at the boundary point p, the velocity u is defined as

$$ u(p) = u_1 \phi_k(p) + u_2 \phi_{k+1} (p), $$

where $$\phi_k = \begin{bmatrix} \phi \\ 0 \end{bmatrix}, \quad \phi_{k+1} = \begin{bmatrix} 0 \\ \phi \end{bmatrix} $$ are corresponding nodal basis functions in the corresponding 2D VectorFunctionSpace.

Let $$ n = \begin{bmatrix} n_1 \\ n_2 \end{bmatrix} \quad \text{and} \quad n_\perp = \begin{bmatrix} n_1 \\ -n_2 \end{bmatrix} $$ be the normal vector and the tangent vector at the boundary at point $p$. Then the TrialFunction that corresponds to the normal v component would be $$ \phi_n = n_1\phi_1 + n_2\phi_2 $$ and $$ \phi_t = n_1\phi_1 - n_2\phi_2 $$ would correspond to the tangential component.

I'm not sure whether FEniCS allows for redefinitions of TrialFunctions or TestFunctions and a corresponding adaptation of how the DirichletBCs are treated.

If one exports the assembled operators to matrices, one can now apply, say, normal and tangential boundary conditions at $p$ as follows.

  1. Replace the k-th column, by $n_1$ times the k-th plus $n_2$-times the k+1-st column (this is as if we had $\phi_n$ in the trial space instead of $\phi_k$).

  2. Replace the k+1-st column, by $n_1$ times the k-th minus $n_2$-times the k+1-st column (this is as if we had $\phi_t$ in the trial space instead of $\phi_{k+1}$).

  3. Then the k-th degree of freedom v_k can be fixed to the normal component of the boundary condition or the k+1-st (v_{k+1}) to the tangential component.

  4. If symmetry is an issue, one can perform the same actions on the columns of the matrices.

  5. The solution in the actual coordinates is then established as $$u(p) = v_kn(p) + v_{k+1}n_\perp (p).$$


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.