Von Neumann boundary condition on internal boundary not working

Hi,

I’m having issues with solving a Poisson equation with a von Neumann boundary at the middle of the mesh. My problem consists of a rectangle, with a plate at y=0. I’m generating the mesh with mshr and my implementation is as follows:

from scipy.special import erf
from collections import *
from dolfin import *
from mshr import *
import matplotlib.pyplot as plt
import matplotlib.tri as tri
import numpy as np
import math
parameters["allow_extrapolation"]= True
parameters["refinement_algorithm"] = "plaza_with_parent_facets"
#Parameters!

VSiO2 = -8.52
Vgraphite = 1.236
B = 1.45

size = 1
Lx = 1*size
H = 0.5*size
d_bg = 0.25*size
d_hbn_bt = 0.035*size
d_hbn_tp = 0.002*size
hbn_height = 0.1*size

domain = Rectangle(Point(-Lx/2, -d_bg), Point(Lx/2, H))
domain.set_subdomain(1, Rectangle(Point(-Lx/2, -d_bg), Point(Lx/2, -d_hbn_bt)))
domain.set_subdomain(2, Rectangle(Point(-Lx/2, -d_hbn_bt), Point(Lx/2, 0)))
domain.set_subdomain(3, Rectangle(Point(-Lx/2, 0.0), Point(Lx/2, d_hbn_tp)))
domain.set_subdomain(4, Rectangle(Point(-Lx/2, d_hbn_tp), Point(Lx/2, H)))

meshr = generate_mesh(domain, 50, “cgal”)
mesh = refine(refine(meshr))
V = FunctionSpace(mesh, “Lagrange”, 1)

def boundary_graph(x):
return near(x[1], -d_hbn_bt, tolc) and between(x[0],(0,Lx/2))

def boundary_bg(x):
return near(x[1], -d_bg, tolc)

bcs = [DirichletBC(V, Constant(VSiO2),boundary_bg),
DirichletBC(V, Constant(Vgraphite), boundary_graph)]

boundary_markers = MeshFunction(“size_t”, mesh, mesh.topology().dim() - 1)
boundary_markers.set_all(0)

class BoundaryBackGate(SubDomain):
tol = 1E-14
def inside(self, x, on_boundary):
return on_boundary and near(x[1], -d_bg, tol)

class BoundaryGraphite(SubDomain):
tol = 1E-14
def inside(self, x, on_boundary):
return on_boundary and near(x[1], -d_hbn_bt, tol) and between(x[0],(0,Lx/2))

class BoundaryX0(SubDomain):
tol = 1E-14
def inside(self, x, on_boundary):
return on_boundary and near(x[0], -Lx/2, tol)

class BoundaryX1(SubDomain):
tol = 1E-14
def inside(self, x, on_boundary):
return on_boundary and near(x[0],Lx/2, tol)

class BoundaryY1(SubDomain):
tol = 1E-14
def inside(self, x, on_boundary):
return on_boundary and near(x[1], H, tol)

class BoundaryPlate(SubDomain):
tol = 1E-14
def inside(self, x, on_boundary):
return near(x[1], 0.0, tol)

bbg = BoundaryBackGate()
bbg.mark(boundary_markers, 0)

bgraph = BoundaryGraphite()
bgraph.mark(boundary_markers, 1)

bx0 = BoundaryX0()
bx0.mark(boundary_markers, 1)

bx1 = BoundaryX1()
bx1.mark(boundary_markers, 2)

by1 = BoundaryY1()
by1.mark(boundary_markers, 3)

bgraphe = BoundaryPlate()
bgraphe.mark(boundary_markers, 4)

dx = Measure(‘dx’, domain = mesh)
ds = Measure(‘ds’, domain=mesh, subdomain_data = boundary_markers)

u = TrialFunction(V)
v = TestFunction(V)

F = dot(grad(u), grad(v))dx+4e*(B/phi0)uuvds(4)
u_ = Function(V) # the most recently computed solution
F = action(F, u_)
J = derivative(F, u_, u)

problem = NonlinearVariationalProblem(F, u_, bcs, J)
solver = NonlinearVariationalSolver(problem)

solver.solve()

If I use any other surface, the von Neumann boundary condition is applied, but not if I pick the y=0 plate.

Can you set up a minimum working code example? Your code does not run as there are typos even in the variational form.