Hi everyone,
I try to implement the Bidomain equations coupled to a surrounding domain. They are to model electrical heart dynamics. I obtain a RuntimeError: Cannot find common parent mesh error which I do not know how to track down…:
File "simple_mesh_example.py", line 93, in <module>
experiment()
File "simple_mesh_example.py", line 89, in experiment
df.solve(a == L, sol, bc)
File "[...]/lib/python3.7/site-packages/dolfin/fem/solving.py", line 233, in solve
_solve_varproblem(*args, **kwargs)
File "[...]/.local/lib/python3.7/site-packages/dolfin/fem/solving.py", line 259, in _solve_varproblem
form_compiler_parameters=form_compiler_parameters)
File "[...]/.local/lib/python3.7/site-packages/dolfin/fem/problem.py", line 126, in __init__
cpp.fem.MixedLinearVariationalProblem.__init__(self, a_list, L_list, u_comps, bcs)
RuntimeError: Cannot find common parent mesh
I just give the code in case you are familiar with the problem, below that there will be some background on the equations and the variational form I used.
import dolfin as df
# mesh creation taken from https://github.com/cdaversin/mixed-dimensional-examples/blob/master/KNPEMI/KNPEMI.py
mesh = df.UnitSquareMesh(1, 1)
interior = df.CompiledSubDomain("std::max(fabs(x[0] - 0.5), fabs(x[1] - 0.5)) < 0.25")
subdomains = df.MeshFunction("size_t", mesh, mesh.topology().dim(), 0)
# mark domain
for cell in df.cells(mesh):
x = cell.midpoint().array()
subdomains[cell] = int(interior.inside(x, False))
assert sum(1 for _ in df.SubsetIterator(subdomains, 1)) > 0
# create interior mesh
mesh_heart = df.MeshView.create(subdomains, 1)
# id torso:0, id heart: 1,
dx_heart = df.Measure("dx", mesh, 1, subdomain_data=subdomains)
dx_torso = df.Measure("dx", mesh, 0, subdomain_data=subdomains)
# heart space
H = df.FunctionSpace(mesh_heart, "Lagrange", 1)
# heart and Torso function space
HuT = df.FunctionSpace(mesh, "Lagrange", 1) # u is german 'und', sorry...
W = df.MixedFunctionSpace(H, HuT)
psi_H, psi_HuT = df.TestFunctions(W)
v_mt, phi_et = df.TrialFunctions(W)
# time discretisation and ionic current
v_mtm1 = df.Function(H)
I_ion = df.Function(H)
# conductivities
M_i = df.Constant(1)
M_e = df.Constant(1)
M_t = df.Constant(1)
# eq1
a1 = (
M_i * df.inner(df.grad(v_mt), df.grad(psi_H)) * dx_heart
+ M_i * df.inner(df.grad(phi_et), df.grad(psi_H)) * dx_heart
)
L1 = -v_mt * psi_H * dx_heart + v_mtm1 * psi_H * dx_heart - I_ion * psi_H * dx_heart
# eq2
a2 = (
M_i * df.inner(df.grad(v_mt), df.grad(psi_HuT)) * dx_heart
+ (M_i + M_e) * df.inner(df.grad(phi_et), df.grad(psi_HuT)) * dx_heart
+ M_t * df.inner(df.grad(phi_et), df.grad(psi_HuT)) * dx_torso
)
L2 = df.Constant(0) * psi_HuT * df.Measure("dx", mesh)
a = a1 + a2
L = L1 + L2
# some Boundary condition to test integration (phi otherwise wit non-zero nullspace)
u_D = df.Expression("1 + x[0]*x[0] + 2*x[1]*x[1]", degree=2)
def boundary(x, on_boundary):
return on_boundary
bc = df.DirichletBC(HuT, u_D, boundary)
sol = df.Function(W)
df.solve(a == L, sol, bc)
Running this gives the error above. I would really appreciate any help!
My dolfin.__version__: 2019.2.0.dev0
Background
The domain looks like this:
With backwards Euler time discretisation, the corresponding equations are
In H:
(1) \frac{v^n - v^{n-1}}{\Delta t} + I_\text{Ion}^{n-1} =\nabla M_i \nabla v^n + \nabla M_i \nabla \phi^n ,
(2) \nabla M_i \nabla v^n + \nabla (M_i + M_e) \nabla \phi^n =0
In T:
(3) \nabla M_t \nabla \phi_t^n = 0
with continuity condition for \phi == \phi_t:
(B.1) \phi = \phi_t on \partial H
(B.2) \vec{n}_h (M_i \nabla v^n+ (M_i + M_e) \nabla \phi^n )= 0 on \partial H
(B.3) \vec{n}_h M_e \nabla \phi = \vec{n}_t M_t \nabla \phi_t on \partial H
The variational forms are (boundary conditions already applied…):
(1) \int_H M_i \nabla v^{n} \nabla \psi \:\text{dx} + \int_H M_i \nabla \phi^{n} \nabla \psi \: \text{dx} = - \int_H \frac{v_n- v_{n-1}}{\Delta t} \psi \: \text{dx}- \int_H I_\text{Ion}^{n-1}\psi\: \text{dx},\: \:\:\:\forall \psi \in V(H)
(2) \int_H M_i \nabla v^{n} \nabla \psi \:\text{dx} + \int_H (M_i + M_e) \nabla \phi^{n} \nabla \psi \:\text{dx} + \int_{\partial H} \vec{n}_h M_e \nabla \phi \psi \text{ds}=0 ,\: \:\:\:\forall \psi \in V(H).
(3) \int_T M_t \nabla \phi_t^{n} \nabla \xi \:\text{dx} + \int_{\partial H} \vec{n}_t M_t \nabla \phi \xi \text{ds} =0,\: \:\:\:\forall \xi \in V(T)
because \phi and \phi_t are continous and \vec{n}_t = - \vec{n}_h, the variational form can be formulated with two functionspaces, one for the Heart domain, one for the whole domain (H+T):
(1) \int_H M_i \nabla v^{n} \nabla \psi \:\text{dx} + \int_H M_i \nabla \phi^{n} \nabla \psi \: \text{dx} = - \int_H \frac{v_n- v_{n-1}}{\Delta t} \psi \: \text{dx}- \int_H I_\text{Ion}^{n-1}\psi\: \text{dx},\: \:\:\:\forall \psi \in V(H)
(2) \int_H M_i \nabla v^{n} \nabla \eta \:\text{dx} + \int_H (M_i + M_e) \nabla \phi^{n} \nabla \eta \:\text{dx} + \int_T M_t \nabla \phi_t^{n} \nabla \eta \:\text{dx} =0 ,\: \:\:\:\forall \eta \in V(H \cup T).
This follows a deduction from Sundnes here with a different time approximation scheme and this is what I used for the implementation above.