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_inner = ( inner(grad(ua), grad(va)) * dxa + inner(grad(ub), grad(vb)) * dxb
- (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"}})