Here’s the code I’ve come up with so far.

```
from fenics import *
from mshr import Rectangle, Circle, generate_mesh
# Create refined mesh for fluid domain
lx = ly = 10
cx = lx / 2.0
cy = 0
radius = 2.5
rectangle_middle = Rectangle(Point(0, -ly/2), Point(lx, ly/2))
hole = Circle(Point(cx, cy), radius)
res = 40
domain = rectangle_middle - hole
mesh_fluid = generate_mesh(domain, res)
# "solve" NS equations
u_fluid = Expression(("cos(x[0])", "sin(x[1])"), degree=2)
V_fluid = VectorFunctionSpace(mesh_fluid, "CG", 1, 2)
u_fluid = project(u_fluid, V_fluid) # this is our solution
# mesh the whole domain (fluid + solid)
thickness = 2
rectangle_bottom = Rectangle(Point(0, -ly/2 - thickness), Point(lx, -ly/2))
rectangle_top = Rectangle(Point(0, ly/2 + thickness), Point(lx, ly/2))
rectangle = rectangle_middle + rectangle_bottom + rectangle_top
domain = rectangle - hole
domain.set_subdomain(1, rectangle_top + rectangle_bottom)
domain.set_subdomain(2, rectangle_middle)
res = 30
mesh_total = generate_mesh(
domain, res)
vm = MeshFunction("size_t", mesh_total, mesh_total.topology().dim(), mesh_total.domains())
dx = Measure("dx", subdomain_data=vm)
# create formulation mimicking HT simulation...
V_total = FunctionSpace(mesh_total, "DG", 1)
u_total = Function(V_total)
v_total = TestFunction(V_total)
F = u_total*v_total*dx - 3*v_total*dx #- u_fluid*v_total*dx(2)
solve(F == 0, u_total, bcs=[])
XDMFFile("out.xdmf").write(u_total)
```