# Coupling discontinuous functions from multiple submeshes

Hello,
I am trying to couple two discontinuous functions from different submeshes. I tried to couple them by using Neumann BCS which doesn’t work at all. Now I am interested to couple them by making use of the jump operator, which returns various errors. I’ll highly appreciate your valuable suggestions, remarks/ comments.

A minimal example is:

from dolfin import *
import matplotlib.pyplot as plt

class DirichletBoundary1(SubDomain):
def inside(self, x, on_boundary):
return x[0] <= 2.0 + DOLFIN_EPS or on_boundary

class DirichletBoundary2(SubDomain):
def inside(self, x, on_boundary):
return x[0] >= 2.0 - DOLFIN_EPS or on_boundary

class Left(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0.0) and on_boundary

class Right(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 4.0) and on_boundary

class Middle(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 2.0)

mesh = IntervalMesh(100, 0, 4)
marker = MeshFunction(“size_t”, mesh, mesh.topology().dim(), 0)

for c in cells(mesh):
marker[c] = c.midpoint().x() < 2.0

submesh1 = MeshView.create(marker, 1)
submesh2 = MeshView.create(marker, 0)
W1 = FunctionSpace(submesh1, “Lagrange”, 1)
W2 = FunctionSpace(submesh2, “Lagrange”, 1)
W = MixedFunctionSpace( W1, W2 )
f = Expression(“x[0]”, degree=1)
(v1,v2) = TestFunctions(W)
u = Function(W)
u1 = u.sub(0)
u2 = u.sub(1)

jump = u1(‘+’) - u2(‘-’)
dx1 = Measure(“dx”, domain=W.sub_space(0).mesh())
dx2 = Measure(“dx”, domain=W.sub_space(1).mesh())
ds = Measure(“ds”, domain=mesh, subdomain_data=marker)
F1 = ( inner ( grad ( v1 ), grad ( u1 ) ) + v1 * u1 ) * dx1
F2 = ( inner ( grad ( v2 ), grad ( u2 ) ) + v2 * u2 ) * dx2
F = F1 + F2 - jump(v1)jumpds(‘+’) - jump(v2)jumpds(‘-’)
L1 = fv1dx1
L2 = fv2dx2
L = L1 + L2
bc1 = DirichletBC(W1, Constant(-2.0), Middle())
bc2 = DirichletBC(W1, Constant(0.0), Left())
bc3 = DirichletBC(W2, Constant(0.0), Right())
bc4 = DirichletBC(W2, Constant(1.0), Middle())
bcs = [bc1,bc4, bc3,bc2]
solve(F - L == 0, u, bcs)
c1 = plot(u1)
plt.show()
c2 = plot(u2)

## and this returns the following error:

NotImplementedError Traceback (most recent call last)
/tmp/ipykernel_81779/1044633422.py in
61 F2 = ( inner ( grad ( v2 ), grad ( u2 ) ) + v2 * u2 ) * dx2
62
—> 63 F = F1 + F2 - jump(v1)jumpds(‘+’) - jump(v2)jumpds(‘-’)
64
65 L1 = fv1dx1

/usr/lib/python3/dist-packages/ufl/exproperators.py in _call(self, arg, mapping, component)
328 return _restrict(self, arg)
329 else:
→ 330 return _eval(self, arg, mapping, component)
331
332

/usr/lib/python3/dist-packages/ufl/exproperators.py in _eval(self, coord, mapping, component)
318 mapping = {}
319 index_values = StackDict()
→ 320 return f.evaluate(coord, mapping, component, index_values)
321
322

/usr/lib/python3/dist-packages/ufl/algebra.py in evaluate(self, x, mapping, component, index_values)
85
86 def evaluate(self, x, mapping, component, index_values):
—> 87 return sum(o.evaluate(x, mapping, component,
88 index_values) for o in self.ufl_operands)
89

/usr/lib/python3/dist-packages/ufl/algebra.py in (.0)
85
86 def evaluate(self, x, mapping, component, index_values):
—> 87 return sum(o.evaluate(x, mapping, component,
88 index_values) for o in self.ufl_operands)
89

/usr/lib/python3/dist-packages/ufl/algebra.py in evaluate(self, x, mapping, component, index_values)
191 tmp = 1
192 for o in ops:
→ 193 tmp *= o.evaluate(x, mapping, (), index_values)
194 return tmp
195

/usr/lib/python3/dist-packages/ufl/restriction.py in evaluate(self, x, mapping, component, index_values)
32
33 def evaluate(self, x, mapping, component, index_values):
—> 34 return self.ufl_operands[0].evaluate(x, mapping, component,
35 index_values)
36

/usr/lib/python3/dist-packages/ufl/core/terminal.py in evaluate(self, x, mapping, component, index_values, derivatives)
58 # If it has an ufl_evaluate function, call it
59 if hasattr(self, ‘ufl_evaluate’):
—> 60 return self.ufl_evaluate(x, component, derivatives)
61 # Take component if any
62 warning(“Couldn’t map ‘%s’ to a float, returning ufl object without evaluation.” % str(self))

/usr/lib/petsc/lib/python3/dist-packages/dolfin/function/function.py in ufl_evaluate(self, x, component, derivatives)
286 else:
287 # Scalar evaluation
→ 288 return self(*x)
289
290 def call(self, *args, **kwargs):

/usr/lib/python3/dist-packages/ufl/core/expr.py in len(self)
388 if len(s) == 1:
389 return s[0]
→ 390 raise NotImplementedError(“Cannot take length of non-vector expression.”)
391
392 def iter(self):

NotImplementedError: Cannot take length of non-vector expression.

Hi Hina,

This can be implemented using `multiphenics`, see this thread for the exact code.

Hope that helps, let me know if you have any questions.

Ingeborg

1 Like