Combine a set of PDE with a fenics mechanics package

2

One approach that is straightforward, albeit perhaps difficult to optimize for large problems, is to just add the ODE system unknowns as “Real”-type components of a mixed space, i.e., global constants. Consider a contrived example (which I just made up for a code demonstration, rather than to be a reasonable model any physical system). The PDE problem is: Find u:\Omega\to\mathbb{R} such that

\begin{array}{lr} \partial_t u - c_1\Delta u = 0 & \text{ on }\Omega\text{ ,}\\ c_1\partial u/\partial\mathbf{n} = q & \text{ on }\partial\Omega\text{ .} \end{array}

Its Neumann BC data is the solution q to the ODE

\frac{d}{dt}q = -\frac{c_2}{\vert\partial\Omega\vert}\int_{\partial\Omega} u\text{ ,}

which is in turn driven by an integral of the PDE solution. This problem can be approximated in FEniCS as follows:

from dolfin import *

# Setup:
N = 64
mesh = UnitSquareMesh(N,N)
cell = mesh.ufl_cell()
uE = FiniteElement("CG",cell,1)
qE = FiniteElement("R",cell,0)
V = FunctionSpace(mesh,MixedElement([uE,qE]))
uq = TrialFunction(V)
uq_old = Function(V)
vr = TestFunction(V)
u,q = split(uq)
u_old,q_old = split(uq_old)
v,r = split(vr)
T = 5.0
Dt = Constant(T/N)
n = FacetNormal(mesh)

# PDE for u on interior, with flux BC in terms of ODE unknown, q:
u_t = (u-u_old)/Dt
c1 = Constant(1e-1)
F_v = (u_t*v + c1*dot(grad(u),grad(v)))*dx - q*v*ds

# ODE for q on boundary, driven by average value of PDE uknown, u:
q_t = (q-q_old)/Dt
c2 = Constant(1e0)
F_r = (q_t + c2*u)*r*ds

# Combined residual:
F = F_v + F_r

# Time stepping:
uqh = Function(V)
x = SpatialCoordinate(mesh)
uq_old.assign(project(as_vector([sin(pi*x[0])*sin(pi*x[1]),Constant(1.0)]),V))
f_out = File("u.pvd")
for i in range(0,N):
    solve(lhs(F)==rhs(F),uqh)
    uq_old.assign(uqh)
    u_out = uqh.split()[0]
    u_out.rename("u","u")
    f_out << u_out
1 Like