# Mixed-dimensional branch - Problem with time integration

Hi fenics community,

I use the mixed dimensional branch of @cdaversin and try to solve a time-dependent problem with the following reduced code:

``````from dolfin import *
from ufl.algorithms import expand_derivatives

# Create mesh
L = 0.20
H = 0.01

mesh = RectangleMesh(Point(0., 0.), Point(L, H), 160, 8, "right")

# Define boundaries and boundary conditions
marker = MeshFunction("size_t", mesh, mesh.topology().dim() - 1, 0)
marker.set_all(99)

class BoundaryLeft(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x, 0 + DOLFIN_EPS)

class BoundaryRight(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x, L - DOLFIN_EPS)  # H-h

boundary_left = BoundaryLeft()
boundary_left.mark(marker, 1)
boundary_mesh_left = MeshView.create(marker, 1)

boundary_right = BoundaryRight()
boundary_right.mark(marker, 2)

dV = Measure("dx", domain=mesh)
dS_l = Measure("dx", domain=boundary_mesh_left)
dS_r = Measure("ds", domain=mesh, subdomain_data=marker)

tau = Expression(('0.', 't <= tc ? 2*t/tc : 0'), t=0, tc=0.01, degree=1)
nue = Expression(('0', '0'), degree=1)

# Define function space
V_v = VectorElement('CG', mesh.ufl_cell(), 1, dim=2)
V_S = TensorElement('DG', mesh.ufl_cell(), 0, shape=(2, 2), symmetry=True)
MX_elem = MixedElement([V_v, V_S])
MX = FunctionSpace(mesh, MX_elem)

LM_left = VectorFunctionSpace(boundary_mesh_left, 'CG', 1, dim=2)

W = MixedFunctionSpace(MX, LM_left)

# Define trial and test functions
w_ = Function(W)
mx, l_left = w_.split()
v, S = split(mx)

dmx, dl_left = TestFunctions(W)
dv, dS = split(dmx)

# Variables of current time step
w_n = Function(W)
mx_n, l_left_n = w_n.split()
v_n, S_n = split(mx_n)

# Time-stepping parameters
T = 0.1
n = 1000
dt = float(T / n)

# Time-stepping
for i in range(n):

print("Time: ", dt * (i + 1))

# Update Neumann boundary conditions
tau.t = float(dt * (i + 1 / 2))

# Residual
v_mid = (v + v_n) / 2
S_mid = (S + S_n) / 2
l_left_mid = (l_left + l_left_n) / 2

R0 = dot(dv, (v - v_n) / dt) * dV + \
inner(grad(dv).T, S_mid) * dV - \
dot(dv, (l_left - l_left_n) / dt) * dS_l - dot(dv, tau) * dS_r(2)

R1 = dot(dl_left, (v - v_n) / dt) * dS_l - dot(dl_left, nue) * dS_l

R2 = inner(dS, (S - S_n) / dt) * dV - inner(dS, grad(v_mid)) * dV

R = R0 + R1 + R2

# Solve for next step
F_blocks = extract_blocks(R)
Jacobi = []
for Fi in F_blocks:
for ui in w_.split():
d = derivative(Fi, ui)
d = expand_derivatives(d)
Jacobi.append(d)

problem = MixedNonlinearVariationalProblem(F_blocks, w_.split(), [], J=Jacobi)
solver = MixedNonlinearVariationalSolver(problem)
solver.solve()

# Update field
mx_n.assign(w_.sub(0))
l_left_n.assign(w_.sub(1))
``````

Unfortunately, the program crashes after a few time iterations and generates a RuntimeError in solver.solve() with the error message:

*** -------------------------------------------------------------------------
*** Error: Unable to evaluate function at point.
*** Reason: The point is not inside the domain. Consider calling “Function::set_allow_extrapolation(true)” on this Function to allow extrapolation.
*** Where: This error was encountered inside Function.cpp.
*** Process: 0

*** DOLFIN version: 2019.2.0.dev0
*** Git changeset: 68d9ed51bee6723163e7c4fbfffe31e6e3079bca
*** -------------------------------------------------------------------------

It is very strange, because the error does not appear every time and sometimes the code works fine. I wonder why the error occurs randomly.

I appreciate any kind of help to solve my problem This was most likely fixed in: fenics-project / DOLFIN / commit / 444128a456b9 — Bitbucket

Hi dokken,

thank you very much for your respond. Is there a docker container including this last commit?
I pulled the docker container mixed dimensional branch of @cdaversin 2 weeks ago.

See for instance: Package fenics · GitHub

1 Like

Hi @dokken,

thank you once again, it seems to work but my code is running much slower now.

Are you running on a Mac? Then you should use the arm image. If not use the non-arm image

I am running on windows via pycharm and I use the arm image. I switched to the non-arm image. It is way better now, thank you .

Nonetheless, my original code does not converge now. Did I mess up something with the variable assignment (like in the last lines of my code)? Do I have to consider something special compared to the branch of @cdaversin?

I dont think so. Could you post an MWE?

Hi @dokken,

for example the following code does not converge at the second time step, when I use the new image.

``````from dolfin import *
import numpy as np
from ufl.algorithms import expand_derivatives

# Material properties
model = 'plane_stress'
E = 4.5e6
nu = 0.4
rho = 28
Lambda = E * nu / (1 + nu) / (1 - 2 * nu)
G = E / (2.0 * (1.0 + nu))
mu = G
if model == 'plane_stress':
Lambda = 2 * mu * Lambda / (Lambda + 2 * mu)

# Create mesh
L = 0.20
H = 0.01

mesh = RectangleMesh(Point(0., 0.), Point(L, H), 160, 8, "right")

# Define boundaries and boundary conditions
marker = MeshFunction("size_t", mesh, mesh.topology().dim() - 1, 0)
marker.set_all(99)

class BoundaryLeft(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x, 0 + DOLFIN_EPS)

class BoundaryRight(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x, L - DOLFIN_EPS)

boundary_left = BoundaryLeft()
boundary_left.mark(marker, 1)
boundary_mesh_left = MeshView.create(marker, 1)

boundary_right = BoundaryRight()
boundary_right.mark(marker, 2)

dV = Measure("dx", domain=mesh)
dS_l = Measure("dx", domain=boundary_mesh_left)
dS_r = Measure("ds", domain=mesh, subdomain_data=marker)

tau = Expression(('0.', 't <= tc ? 2e3*t/tc : 0'), t=0, tc=0.01, degree=1)
nue = Expression(('0', '0'), degree=1)

# Define function space
V_v = VectorElement('CG', mesh.ufl_cell(), 1, dim=2)
V_S = TensorElement('DG', mesh.ufl_cell(), 0, shape=(2, 2), symmetry=True)
V_F = TensorElement('DG', mesh.ufl_cell(), 0, shape=(2, 2))
V_u = VectorElement('CG', mesh.ufl_cell(), 1, dim=2)
MX_elem = MixedElement([V_v, V_S, V_F, V_u])
MX = FunctionSpace(mesh, MX_elem)

LM_left = VectorFunctionSpace(boundary_mesh_left, 'CG', 1, dim=2)

W = MixedFunctionSpace(MX, LM_left)

# Define trial and test functions
w_ = Function(W)
mx, l_left = w_.split()
v, S, F, u = split(mx)

dmx, dl_left = TestFunctions(W)
dv, dS, dF, du = split(dmx)

# Variables of current time step
w_n = Function(W)
mx_n, l_left_n = w_n.split()
v_n, S_n, F_n, u_n = split(mx_n)

MX_F = FunctionSpace(mesh, V_F)
assign((w_.sub(0)).sub(2), project(Identity(MX.mesh().geometry().dim()), MX_F))
assign((w_n.sub(0)).sub(2), project(Identity(MX.mesh().geometry().dim()), MX_F))

# Time-stepping parameters
T = 0.1
n = 1000
dt = float(T / n)

# Time-stepping
for i in range(n):

print("Time: ", dt * (i + 1))

# Update Neumann boundary conditions
tau.t = float(dt * (i + 1 / 2))

# Residual
v_mid = (v + v_n) / 2
S_mid = (S + S_n) / 2
F_mid = (F + F_n) / 2
l_left_mid = (l_left + l_left_n) / 2

R0 = dot(dv, rho * (v - v_n) / dt) * dV + \
inner(dot(grad(dv).T, F_mid), S_mid) * dV - \
dot(dv, (l_left - l_left_n) / dt) * dS_l - dot(dv, tau) * dS_r(2)

R1 = dot(dl_left, (v - v_n) / dt) * dS_l - dot(dl_left, nue) * dS_l

R2 = inner(dS, -(Lambda / (2 * mu * (3 * Lambda + 2 * mu))) * tr((S - S_n) / dt) * Identity(2) + 1 / (2 * mu) * (S - S_n) / dt) * dV - inner(dS, dot(F_mid.T, grad(v_mid))) * dV

R3 = inner(dF, (F - F_n) / dt) * dV - \

R4 = dot(du, (u - u_n) / dt) * dV - \
dot(du, v_mid) * dV

R = R0 + R1 + R2 + R3 + R4

# Solve for next step
F_blocks = extract_blocks(R)
Jacobi = []
for Fi in F_blocks:
for ui in w_.split():
d = derivative(Fi, ui)
d = expand_derivatives(d)
Jacobi.append(d)

problem = MixedNonlinearVariationalProblem(F_blocks, w_.split(), [], J=Jacobi)
solver = MixedNonlinearVariationalSolver(problem)
solver.solve()

# Update field
mx_n.assign(w_.sub(0))
l_left_n.assign(w_.sub(1))
``````

It is strange and I have probally mess up something with the variable assignment or the split function?

This is far from and MWE (minimal working example), and I do not have time to parse all this code.

One thing that stands out is:
Why do you redefine the whole variational form inside the loop?

I believe that should be split(w). This holds in all the different places you use split

Hi @dokken,

I am sorry, I tried to create a MWE.

Using the following code

``````from dolfin import *

# Create mesh
mesh = RectangleMesh(Point(0., 0.), Point(0.2, 0.01), 160, 8, "right")

# Define boundaries
marker = MeshFunction("size_t", mesh, mesh.topology().dim() - 1, 0)
marker.set_all(99)

class BoundaryLeft(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x, 0 + DOLFIN_EPS)

boundary_left = BoundaryLeft()
boundary_left.mark(marker, 1)
boundary_mesh_left = MeshView.create(marker, 1)

dV = Measure("dx", domain=mesh)
dS_l = Measure("dx", domain=boundary_mesh_left)

# Define function space
V_v = VectorElement('CG', mesh.ufl_cell(), 1, dim=2)
V_S = TensorElement('DG', mesh.ufl_cell(), 0, shape=(2, 2), symmetry=True)
MX_elem = MixedElement([V_v, V_S])
MX = FunctionSpace(mesh, MX_elem)

LM_left = VectorFunctionSpace(boundary_mesh_left, 'CG', 1, dim=2)

W = MixedFunctionSpace(MX, LM_left)

# Define trial and test functions
w_ = Function(W)
mx, l_left = w_.split()
v, S = split(mx)

dmx, dl_left = TestFunctions(W)
dv, dS = split(dmx)

# Weak form
tau = Expression(('0.', '20'), degree=1)

R0 = dot(dv, v) * dV + inner(grad(dv).T , S) * dV - \
dot(dv, l_left) * dS_l - dot(dv, tau) * dV

R1 = dot(dl_left, v) * dS_l

R2 = inner(dS, tr(S) * Identity(2) + S) * dV - \

R = R0 + R1 + R2

# Solution
solve(R==0, w_)
``````

works well using the branch of @cdaversin. But the solver does not converge in case of the new image.

Hi,

I tested the MWE posted above and I also cannot get it to work on ghcr.io/scientificcomputing/fenics:2023-02-17 or on ceciledc/fenics_mixed_dimensional:16-01-23. As FEM_TT says the Newton solver fails to converge. I will post a bug report on bitucket.

It does work and give the correct result using e.g. ceciledc/fenics_mixed_dimensional:21-06-22 - but then your original problem remains wrt the program crashing You could assemble the code using fenics_ii instead:

``````from xii import *
from xii.nonlin.jacobian import block_jacobian
from fenics import *

# Create mesh
mesh = RectangleMesh(Point(0., 0.), Point(0.2, 0.01), 160, 8, "right")

# Define boundaries
marker = MeshFunction("size_t", mesh, mesh.topology().dim() - 1, 0)
marker.set_all(99)

class BoundaryLeft(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x, 0 + DOLFIN_EPS)

boundary_left = BoundaryLeft()
boundary_left.mark(marker, 1)
boundary_mesh_left = EmbeddedMesh(marker, 1)

dV = Measure("dx", domain=mesh)
dS_l = Measure("dx", domain=boundary_mesh_left)

# Define function space
V_v = VectorFunctionSpace(mesh, 'CG', 1, 2)
V_S = TensorElement('DG', mesh.ufl_cell(), 0, shape=(2, 2), symmetry=True)
V_S = FunctionSpace(mesh, V_S)

LM_left = VectorFunctionSpace(boundary_mesh_left, 'CG', 1, dim=2)

W = [V_v, V_S, LM_left] # list of function spaces

# Define trial and test functions
w = ii_Function(W)
v, S, l_left = w # w is a list of functions

dv, dS, dl_left = map(TestFunction, W)

Tdv = Trace(dv, boundary_mesh_left)
Tv = Trace(v, boundary_mesh_left)

# Weak form
tau = Expression(('0.', '20'), degree=1)
tau = interpolate(tau, V_v) #fenics_ii later needs to know the "type" of tau, that is why we interpolate it

R0 = dot(dv, v) * dV + inner(grad(dv).T , S) * dV + \
- dot(Tdv, l_left) * dS_l \
- dot(dv, tau) * dV

R1 = dot(dl_left, Tv) * dS_l

R2 = inner(dS, tr(S) * Identity(2) + S) * dV - \

F = [R0, R2, R1]

# Solution
dF = block_jacobian(F, w)

# Newton
eps = 1.0
tol = 1.0E-10
niter = 0
maxiter = 125

dw = ii_Function(W)
while eps > tol and niter < maxiter:
niter += 1

A, b = list(map(ii_assemble, (dF, F)))
A, b = list(map(ii_convert, (A, b)))

solve(A, dw.vector(), b)

eps = sqrt(sum(x.norm('l2')**2 for x in dw.vectors()))

print('\t%d |du| = %.6E |A|= %.6E |b| = %.6E' % (niter, eps, A.norm('linf'), b.norm('l2')))

for i in range(len(W)):
w[i].vector().axpy(-1, dw[i].vector())
``````

You can install fenics_ii using pip (as described in the readme, feel free to message me if you need assistance) or run on the docker image mirok/emi-book-fem 2 Likes

thank you very much, that is very kind.
Unfortunately, I can not find the modul “xii” using the docker image mirok/emi-book-fem:latest. Do you know, what might be the problem?

Aha, the `mirok/emi-book-fem` image is a bit old and uses python2.

Can you try this image: `ghcr.io/ingeborggjerde/graphnics:v0.8.0` ?

Sorry for the late reply, I am travelling at the moment. I’ll update my post in a few days with pip installation instructions 2 Likes

thank you again, your image seems to work perfectly. Now, I will try to solve my original problem and I will let you know, if it works . I wish you a nice trip.