# Internal boundary condition with mixed_dimensional branch

Hi everyone,

I would like to solve a problem consisting of two different equations on different subdomains coupled by nonlinear boundary condition of Neumann type. Using the mixed_dimensional branch (based on the answer to this question Coupled Equations Over Multiple Domains) I tried to add the internal boundary condition to the code, but this gives me the error : “Error: Unable to evaluate function at point.”, even though I only evaluate at common points of both submeshes. Maybe someone could give my a hint on how to correctly couple equation on different domains?
Thank you!

Code:

``````from dolfin import *
import matplotlib.pyplot as plt

# Width of total domain
width = 1.0
# Height of total domain
height = 1.0
# Total resolution
resolution = 16
# Forcing function
f = Constant(10)
# Boundary values for the two equations
a_0 = 0
b_0 = 0

# Define SubDomain
class LeftDomain(SubDomain):
def inside(self, x, on_boundary):
return between(x[0], (0, width/2))

# Define Boundary
class LeftBoundary(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], 0.0)

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

class Interface(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], width/2)

# Mesh for the entire domain
mesh_full = RectangleMesh(Point(0, 0),
Point(width, height),
resolution,
resolution)

# Define domains
domains = MeshFunction("size_t", mesh_full, 2)
domains.set_all(0)
left_domain = LeftDomain()
left_domain.mark(domains, 1)

# Create mesh for the left domain
mesh_left = MeshView.create(domains, 1)
mesh_right = MeshView.create(domains, 0)

# Define boundaries

left_boundary = LeftBoundary()
right_boundary = RightBoundary()

boundary_right = MeshFunction("size_t", mesh_right, 1)
right_boundary.mark(boundary_right, 1)
boundary_left = MeshFunction("size_t", mesh_left, 1)
left_boundary.mark(boundary_left, 1)

interf = Interface()
boundary_int = MeshFunction("size_t", mesh_full, 1)
boundary_int.set_all(0)
interf.mark(boundary_int, 1)

mesh_int = MeshView.create(boundary_int, 1)

# Create the function spaces
W_a = FunctionSpace(mesh_right, "Lagrange", 1)
W_b = FunctionSpace(mesh_left, "Lagrange", 1)
W = MixedFunctionSpace(W_a, W_b)

#dS = Measure("dS", domain=mesh_full, subdomain_data=boundary_int)

ds_right = Measure("ds", domain=mesh_right,  subdomain_data=boundary_int)
ds_int = ds_right(1)
ds_int = Measure("dx", domain=mesh_int)

# Define boundary conditions
bc_a = DirichletBC(W.sub_space(0), a_0, boundary_right, 1)
bc_b = DirichletBC(W.sub_space(1), b_0, boundary_left, 1)
bcs = [bc_a, bc_b]

# Create all of the weak form variables
u = Function(W)
ua, ub = u.sub(0), u.sub(1)
va, vb = TestFunctions(W)

# Need to manually specify these.
dxa = Measure("dx", domain=W.sub_space(0).mesh())
dxb = Measure("dx", domain=W.sub_space(1).mesh())

g = Expression("1.0", degree=2)
c = Constant("10.0")

# The weak formulation
- (f * va * dxa + f * vb * dxb )
)

F_coupling = (va-vb)*(c + (ua - ub)**2)*ds_int

F = F_inner + F_coupling

solve(F==0, u, bcs, solver_parameters= {"nonlinear_solver": "newton",
"newton_solver": {"linear_solver": "gmres"}})``````
2 Likes

Hi,

Sorry for the late reply. My last commit in the mixed-dimensional branch should fix the issue. (If you are using the Docker container, the image `ceciledc/fenics_mixed_dimensional:latest` is up-to-date)

1 Like

Hi @cdaversin,
thank you very for your reply and you effort to fix the problem. Unfortunately I am still not able to impose the internal boundary condition. How do I correctly chose the measure of the internal interface? If I use the corresponding part of the right mesh:

``````ds_right = Measure("ds", domain=mesh_right,  subdomain_data=boundary_int)
ds_int = ds_right(1)
``````

and a simplified boundary term:

``````F_coupling =  va*c*ds_int + vb*c*ds_int
``````

only the solution of the right mesh is affected. Additionally, if I combine trial and test functions from both meshes, e. g.

``````F_coupling =  va*ub*ds_int
``````

my kernel dies without any error message. Any idea how to solve this?

The measure on the internal interface should be with type “dx”, as you did in the code you sent previously :

``````ds_int = Measure("dx", domain=mesh_int)
``````

Indeed if you use `ds_int` defined from `ds_right`, it will affect only the right mesh since `ds_right` is defined from `mesh_right` (and if you combine trial and test functions with both meshes, one of them is not defined on `mesh_right`).
If `ds_int` is defined from `mesh_int` instead, the appropriate mappings between the interface and the left and right meshes are built to take into account the contributions from both sides.

Hi @cdaversin,
thank you for the explanation. However, using the suggested measure does not solve the problem, since the kernel reliably dies when functions of both meshes are combined (e.g. “F_coupling = va * ub * ds_int”). Any idea why this happens?

I was able to reproduce the issue, and it turned out that in this particular case one of the additional mappings needed for the coupling term was not properly built.
It should work fine now with the last versions of the mixed-dimensional branch / latest tagged Docker container. Thanks for your feedback, that’s helpful to make the framework more robust.

Hi @cdaversin and @causem,

I tried the code posted above by @causem using the latest mixed dimensional branch and I received the following error:

``````Traceback (most recent call last):
File "Main.py", line 111, in <module>
"newton_solver": {"linear_solver": "gmres"}})
File "/usr/local/lib/python3.6/dist-packages/fenics_dolfin-2019.2.0.dev0-py3.6-linux-x86_64.egg/dolfin/fem/solving.py", line 233, in solve
_solve_varproblem(*args, **kwargs)
File "/usr/local/lib/python3.6/dist-packages/fenics_dolfin-2019.2.0.dev0-py3.6-linux-x86_64.egg/dolfin/fem/solving.py", line 298, in _solve_varproblem
solver.solve()
RuntimeError:

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** 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
*** -------------------------------------------------------------------------
``````

When I get rid of the F_coupling it works. Any idea why I am getting this error?