Mesh movement calculation for FSI - boundary conditions

In a fluid structure interaction solver (written in c++, latest stable version of FEniCS available from docker) I need to solve two problems, one for the continuum (i.e. the fluid and the solid), and one for the mesh, in which the amount by which the mesh needs to move in the fluid domain is calculated (in the solid domain, the mesh moves with the material). I am trying to solve these two problems in a single call to the solve function.

To do this I need to apply a boundary condition to the mesh displacement such that at the boundary between the solid and the fluid the mesh moves with the solid, and within the fluid domain I solve (for example) div(grad(um)) == 0, where um is the mesh displacement.

If u is the displacement of the solid, then I need to apply a boundary condition um == u on the boundary between the fluid and solid domains. My problem is that if I try to do this using DirichletBC(W.sub(3), u, fluid_solid_bdry) then it does not treat u in the boundary condition as a variable, but simply takes the value that is stored within the function u from the previous time-step (note that W.sub(3) is the part of the mixed function space corresponding to um).

Does anyone know a way of applying a boundary condition to um such that it is “slaved” to u at the boundary between the fluid and the solid, treating u as a variable to be calculated? I am not sure if this is possible as it is important that the mesh problem remains “secondary” to the continuum problem, and does not in any way influence the dynamics of the fluid structure interaction.

You can actually set up a variational problem in which there is a single continuous displacement field. There is a complete formulation of this type (albeit using a somewhat difficult notation) in Thomas Wick’s paper here (see (15)–(17)):

https://doi.org/10.1016/j.compstruc.2011.02.019

The main “trick” I would anticipate needing to implement this in FEniCS would be assembling the fluid–solid interface term on the appropriate interior facets, since it requires knowing which of the "+" or "-" sides is the fluid subdomain, to use the correct normal vector. This can be done using @MiroK’s answer here:

(I’ve been meaning to sit down and implement this for some time, as an instructive example, but still haven’t gotten around to it.)