# Combine a set of PDE with a fenics mechanics package

Morning, everybody.
I’m trying to solve a blood flow problem using the lumped parameters model. The flow equations are governed by ODEs and I used a 3D model to represent my left ventricle. Is it possible to solve ODEs coupled with the fenics mechanics package? If yes how do you have ? Has anyone ever solved a problem ??
I would also like to modify the constitutive laws to introduce a new factor in the equations (the phenomenon of growth of the cardiac muscle).
Please can someone help me with the steps or a bit of code?

Thank you in advance, I hope I have been clear in explaining my problem. On this I remain open for more informations

renethierrydjoumessi@yahoo.com

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)

# 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

thank you so much kamensky, let me try and give you a feedback