Coupled System with a single weak formulation

Hey everyone. I have a coupled system in the following form:
Gleichungssystem
which is to be calculated on the following region:
Volumen
The biharmonic equation is supposed to be zero on the boundary and shall only solved on the top, of the cube, the volume under that is supposed to be described by the heat equation.

I tried to couple them by iterative solving: i first solved the biharmonic equation, then the heat equation and transfered the values from one equation with its own functionspace to the other and so on, however the results dont match as i get different results every time…if i even get results.
My question is now, if there is a way to formulate the problem in one equation, which i can solve over time.
If someone can help me, i would be very happy.
Thanks in regard!

If I understand you correctly, then you want to use a mixed finite element method to solve your problem.
I would recommend taking a look at these two demos

https://bitbucket.org/fenics-project/dolfin/src/master/python/demo/documented/mixed-poisson/demo_mixed-poisson.py.rst

https://bitbucket.org/fenics-project/dolfin/src/master/python/demo/documented/stokes-taylor-hood/documentation.rst

This should give you an idea, how to use mixed elements in FEniCS.
But be aware that the first boundary condition in your system cannot be implemented strongly via a DirichletBC in this case, but that you will have to implement it weakly, e.g., using Nitsche’s method

Hi,

yes a mixed finite element method came to my mind, however is still struggle with finding a single formulation to describe to whole system and applying the equations on the right places. Namely the biharmonic equation on the top, with its boundaries and the rest of the cube with the heat equation, with its boundary.
I know about submeshes, which i used for the iterative approach, but a mixed functionspace expects me to define it over one mesh alone.

I am personally not very familiar with solving PDEs defined on different domains. The functionality for this is implemented in the development version of fenics, with the MeshView feature.

Maybe take a look here: https://arxiv.org/pdf/1911.01166.pdf
and here MeshView and MixedFunctionSpace

Also, I would suggest that you supply a minimal working example (e.g. using your iterative solver)

Hey, i am making small steps in the mixed element approach, however i can provide the iterative version, it works in terms of not resulting in erroers, however the results are not stable. It is not possible to reproduce and the chances are high, that the array, which contains the values is filled with nan entries.

The complete Code:

### General Declarations ###
from __future__ import print_function
from fenics import *
import matplotlib.pyplot as plt
import numpy as np

parameters['allow_extrapolation'] = True

### declare mesh and submesh ###

height = 1.0
width = 1.0
length = 1.0
mesh = BoxMesh(Point(0.0, 0.0, 0.0), Point(length,width,height), 9, 9, 9)
## Create submesh ##

marker = MeshFunction("size_t", mesh, mesh.topology().dim() - 1, 0)
for f in facets(mesh):
    marker[f] = height - DOLFIN_EPS < f.midpoint().z() < height + DOLFIN_EPS
    
submesh = MeshView.create(marker, 1)    

### Computation of Eigenvalues ###

X = FiniteElement('CG', submesh.ufl_cell(), 2)
Y = FiniteElement('CG', submesh.ufl_cell(), 2)
Z = X * Y
S = FunctionSpace(submesh, Z)

#Computation Eigenvalues
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
#%% ### The Coupled Problem ###
num_steps = 20    # number of time steps
dt = 0.005         # time step size
T = dt * num_steps            # final time
t = 0

################ Functions/TrialFunctions/TestFunctions of the Plate ################
u_m_derivative_Plate = Function(S)
u_n_derivative_Plate = Function(S)
u_m_Plate = Function(S)
u_n_Plate = Function(S)
u_Plate = TrialFunctions(S)
phi_Plate = TestFunctions(S)
normalDerivative_Damping = Function(S)
################ FunctionSpace Function/TrialFunctions/TestFunctions damping ################
V_Damping = FunctionSpace(mesh, 'CG', 2)

v_Damping = TrialFunction(V_Damping)
psi_Damping = TestFunction(V_Damping)
boundaryFunction_Damping = Function(V_Damping)
Solution_Damping = Function(V_Damping)

#for plotting
valueArray_damping = np.array([])
t_plot = np.linspace(0, T, num = num_steps)

############## Boundaries Damping ################
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 

boundaryValue_Damping = Constant(0.0)
bcTop = DirichletBC(V_Damping, boundaryValue_Damping, Top)
bcContainer = DirichletBC(V_Damping, boundaryValue_Damping, Container)

bcs_Damping = [bcTop, bcContainer]

v_initial_Damping = interpolate(boundaryValue_Damping, V_Damping)

############## Boundary Plate ################
def boundary(x, on_boundary):
    return on_boundary

bc_Plate = DirichletBC(S, Constant((0,0)), boundary)

############## bilinear/linear Form Plate and load ################
a_Plate = -dt*dt*(inner(grad(u_Plate[0]),grad(phi_Plate[0])) + inner(grad(u_Plate[1]),grad(phi_Plate[1])) + u_Plate[1]*phi_Plate[0])*dx + (u_Plate[0]*phi_Plate[1])*dx
f_Plate = Constant(0.0)
L_Plate = (2*u_n_Plate[0]*phi_Plate[1] - u_m_Plate[0]*phi_Plate[1])*dx + dt*dt*(inner(f_Plate*ef[0],phi_Plate[1]))*dx 

solution_Plate=Function(S)
l_derivative_Plate = Function(S)
############## bilinear/linear Form Damping and load ################
f_Damping = Constant(0.0)
a_Damping = (inner(v_Damping,psi_Damping) + 2*dt * inner(grad(v_Damping), grad(psi_Damping)))*dx
L_Damping = (inner(v_initial_Damping,psi_Damping) + dt*inner(f_Damping,psi_Damping))*dx

############# Start Calculations ################
for n in range(num_steps): 
    #compute solution Damping    
    solve(a_Damping == L_Damping, Solution_Damping, bcs_Damping, solver_parameters={"linear_solver":"direct"})
    v_initial_Damping.assign(Solution_Damping)
    Solution_DampingZDerivative = Solution_Damping.dx(2)
    plt.figure()
    plot(Solution_DampingZDerivative)
    normalDerivative_Damping = project(Solution_DampingZDerivative, S.sub(0).collapse(), solver_type="mumps")
    # Compute solution Plate
    solve(a_Plate == L_Plate, solution_Plate, bc_Plate, solver_parameters={"linear_solver":"direct"})
    #save to plot
    pointValue = solution_Plate(0.5,0.5,height)   
    valueArray_damping = np.append(valueArray_damping, pointValue[0])    
    u_m_Plate.assign(u_n_Plate)
    u_m_derivative_Plate.assign(u_m_Plate)
    #Update previous solution
    u_n_Plate.assign(solution_Plate)
    u_n_derivative_Plate.assign(solution_Plate)
    
    l_derivative_Plate.assign((u_n_derivative_Plate - u_m_derivative_Plate)/2*dt)
    
    speed, speedAppendix = l_derivative_Plate.split()    
    boundaryFunction_Damping = interpolate(speed, V_Damping) 
    bcTop = DirichletBC(V_Damping, boundaryFunction_Damping, Top)
    bcContainer = DirichletBC(V_Damping, boundaryValue_Damping, Container)
    
    bcs_Damping = [bcTop, bcContainer]
    #Update current time
    t += dt
    
    f_Plate = Constant(2*r*sin(r*t))
    L_Plate = (2*u_n_Plate[0]*phi_Plate[1] - u_m_Plate[0]*phi_Plate[1])*dx + dt*dt*(inner(f_Plate*ef[0],phi_Plate[1]) - inner(normalDerivative_Damping, phi_Plate[0]))*dx 
    #Update time of the expression
    f_Plate.t = t
############# Plotting ################    
plt.figure()
plt.plot(t_plot, valueArray_damping, 'b-')
plt.grid(axis='y')
plt.xlabel('t')
plt.ylabel('u',rotation=0)
plt.legend(['Entwicklung Auslenkung'], loc='upper left')
plt.savefig('GekoppeltesSystem.pdf')