Slip (or similar) boundary conditions

[MER] We are interested in implementing functionality for supporting in
particular strong slip boundary conditions.

A way of supporting this would be to introduce something like

  bc = DirectionalBC(V, direction, value, domain)

or maybe

  bc = RotatedDirichletBC(V, direction, value, domain)

We would also provide an averaged normal

  n = AveragedFacetNormal(mesh)

formed by averaging facet normals appropriately.

The naming is very open for discussion.

A typical use case would be (for slip boundary conditions for Stokes)

  V = VectorFunctionSpace(mesh, "CG", 2)
  Q = FunctionSpace(mesh, "CG", 2)
  W = V*Q

  A = assemble(...)
  b = assemble(...)
  n = AveragedFacetNormal(mesh)
  slip_bc = DirictionalBC(W.sub(0), n, 0.0, "on_boundary")
  slip_bc.apply(A, b)

The implementation of the above would typically involve

  * Creating a rotation matrix R from the directional argument
  * Forming A' = R A R^T
  * Applying "standard" essential boundary conditions to A'
  * Rotate back: A = R^T A' R

Note that this would involve a global matrix-matrix-multiply operation
(possibly some backends support this better than others)

An alternative would be to do this locally during assembly; in the
same vein as assemble_system today. Essentially, you would rotate the
(dense) local element matrix, apply boundary conditions locally, and
rotate back prior to global assembly. So, instead of the above, one
would extend assemble_system and do

    assemble_system(..., ..., slip_bc)

Any opinions or other suggestions?

[Johan Hoffman] This is already implemented in Unicorn as a class SlipBC as a subclass of dolfin::BoundaryCondition (see pp.351 of the FEniCS book, and see also the thesis of Murtazo Nazarov from 2011), where we employ a local version. We have used it for several years with very good results. To include this functionality in dolfin I suggest to have a look at the Unicorn code, which to a large extent could probably be used as it is.

[MER] Thanks for the links Johan. I have taken a look at the Unicorn implementation. As far as I can see there the slip boundary condition is applied after global assembly, so that would correspond to the global approach I guess? Also, there are (at least) two classes there: SlipBC and NewSlipBC, which of these is the relevant one?

[JH] Yes you are right; it is applied after assembly so it should be global in your terminology (but with the rotation matrix formed locally for each node at a time). You should look at "SlipBC" ("NewSlipBC" is connected to a Newton formulation of the system as far as I understand).