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.