Hi Dokken,

I want to couple the biharmonic equation with the stokes equation. The displacement and the following derivative in time of the plate which is modelled by the biharmonic equation is supposed to be the boundary condition of the fluid system which is modelled by the stokes equation.

Following code produces quite good results for a evolution in time:

```
mesh = UnitSquareMesh(16,16)
X = FiniteElement('CG', mesh.ufl_cell(), 2)
Y = FiniteElement('CG', mesh.ufl_cell(), 2)
Z = X * Y
S = FunctionSpace(mesh, Z)
u = TrialFunctions(S)
phi = TestFunctions(S)
a = inner(grad(u[1]), grad(phi[0]))*dx + inner(grad(u[0]), grad(phi[1]))*dx
def boundary(x, on_boundary):
return on_boundary
bc = DirichletBC(S, Constant((0,0)), boundary)
L_dummy = Constant(0)*phi[0]*dx
m = u[0]*phi[0]*dx + u[1]*phi[1]*dx
A, _ = assemble_system(a,L_dummy,bc)
B = assemble(m)
eigensolver = SLEPcEigenSolver(as_backend_type(A), as_backend_type(B))
eigensolver.parameters['spectrum'] = 'target real'
eigensolver.parameters['tolerance'] = 1e-6
eigensolver.solve()
r, c, rx, cx = eigensolver.get_eigenpair(0)
print("Eigenvalue: ", r)
ef = Function(S)
ef.vector()[:] = rx
#plot(ef[0])
#
# %% Computation of time dependent problem
#
T = 1.5 # final time
num_steps = 5000 # number of time steps
dt = T / num_steps # time step size
t = 0
u_m = Function(S)
u_n = Function(S)
u = TrialFunctions(S)
phi = TestFunctions(S)
# Create VTK file for saving solution
vtkfile = File('solution/ResonatingPlate.pvd')
#define F and assemble unknown and known variables
a_evo = -dt*dt*(inner(grad(u[0]),grad(phi[0])) + inner(grad(u[1]),grad(phi[1])) + u[1]*phi[0])*dx + (u[0]*phi[1])*dx
f = Constant(0.0)
L = (2*u_n[0]*phi[1] - u_m[0]*phi[1])*dx + dt*dt*(inner(f*ef[0],phi[0]))*dx
l=Function(S)
for n in range(num_steps):
# Compute solution
solve(a_evo == L, l, bc)
# Save to file
vtkfile << (l, t)
#vtkfile2 << (w_temp,t)
u_m.assign(u_n)
# Update previous solution
u_n.assign(l)
# Update current time
t += dt
f = Constant(2*r*sin(r*t))
L = (2*u_n[0]*phi[1] - u_m[0]*phi[1])*dx + dt*dt*(inner(f*ef[0],phi[0]))*dx
# Update time of the expression
f.t = t
```

Now i need to find a way of using the results in the following codesection for the stokes system:

```
# Test for PETSc or Tpetra
if not has_linear_algebra_backend("PETSc") and not has_linear_algebra_backend("Tpetra"):
info("DOLFIN has not been configured with Trilinos or PETSc. Exiting.")
exit()
if not has_krylov_solver_preconditioner("amg"):
info("Sorry, this demo is only available when DOLFIN is compiled with AMG "
"preconditioner, Hypre or ML.")
exit()
if has_krylov_solver_method("minres"):
krylov_method = "minres"
elif has_krylov_solver_method("tfqmr"):
krylov_method = "tfqmr"
else:
info("Default linear algebra backend was not compiled with MINRES or TFQMR "
"Krylov subspace method. Terminating.")
exit()
height = 1.0
width = 1.0
length = 1.0
mesh = BoxMesh(Point(0.0, 0.0, 0.0), Point(length,width,height), 16, 16, 16)
# Build function space
P2 = VectorElement("Lagrange", mesh.ufl_cell(), 2)
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
TH = P2 * P1
W = FunctionSpace(mesh, TH)
def Top(x, on_boundary): return x[2] > height - DOLFIN_EPS
def Container(x, on_boundary):
return x[0] > length - DOLFIN_EPS or x[0] < DOLFIN_EPS or x[1] > width - DOLFIN_EPS or x[1] < DOLFIN_EPS #or x[2] < DOLFIN_EPS
# No-slip boundary condition for velocity on each facet
noslip = Constant((0.0, 0.0, 0.0))
bcContainer = DirichletBC(W.sub(0), noslip, Container)
```

Here it gets tricky, i basically want to use use the results as a cover wihich stimulates the fluid, the resulting pressure that works back at the plate is an additional task. Moreover i am right now not sure how to handle the preconditioner matrices.

```
# Coupling with plate
#inflow = Expression(('0.0, 0.0, displacement*x[2]'), degree = 2)
bcTop = DirichletBC(W.sub(0), (0.0, 0.0, displacement), Top)
# Collect boundary conditions
bcs = [bcTop, bcContainer]
# Define variational problem
(u, p) = TrialFunctions(W)
(v, q) = TestFunctions(W)
f = Constant((0.0, 0.0, 0.0))
a = inner(grad(u), grad(v))*dx + div(v)*p*dx + q*div(u)*dx
L = inner(f, v)*dx
# Form for use in constructing preconditioner matrix
b = inner(grad(u), grad(v))*dx + p*q*dx
# Assemble system
A, bb = assemble_system(a, L, bcs)
# Assemble preconditioner system
P, btmp = assemble_system(b, L, bcs)
# Create Krylov solver and AMG preconditioner
solver = KrylovSolver(krylov_method, "amg")
# Associate operator (A) and preconditioner matrix (P)
solver.set_operators(A, P)
# Solve
U = Function(W)
solver.solve(U.vector(), bb)
# Get sub-functions
u, p = U.split()
# Save solution in VTK format
ufile_pvd = File("CoupledSystem/Velocitiy/velocity.pvd")
ufile_pvd << u
pfile_pvd = File("CoupledSystem/Pressure/pressure.pvd")
pfile_pvd << p
```