Solve Poisson equation with non-linear boundary conditions on a submesh

variational formulation
2

The problem is that you are not setting mesh markers properly on your submesh. This can be seen by visualizing them after you have marked them.
For this particular case, the issue is that you are using an unstructured mesh, which means that the submesh boundary is not a straight line.

I’ve added a code using built in meshes, that gives the correct and expected results.

import matplotlib.pyplot as plt
from fenics import *
import numpy as np
from ufl import nabla_div
from ufl import sinh


f = Constant(0.0)

i0 = 12.6
r = Constant(2*i0)

iApp = -100
g = iApp

F = 96485.33
R = 8.314
T = 300
C = Constant(F/(2*R*T))

kappa = 1

mesh = RectangleMesh(Point(0,0), Point(100e-6, 400e-6), 100,100)
# from mshr import *
# domain = Rectangle(Point(0,0), Point(100e-6, 400e-6))
# mesh = generate_mesh(domain,100)

regions = MeshFunction('size_t', mesh, mesh.topology().dim())
boundaries = MeshFunction('size_t', mesh, mesh.topology().dim() - 1)

class A(SubDomain):
    def inside(self, x, on_boundary):
        return x[1] <= 200e-6
a = A()

class B(SubDomain):
    def inside(self, x, on_boundary):
        return x[1] >= 200e-6
b = B()

class BottomBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 0.0)
bottom = BottomBoundary()

class TopBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 200e-6)
top = TopBoundary()

regions.set_all(0)
b.mark(regions, 1)

boundaries.set_all(0)
bottom.mark(boundaries, 1)
top.mark(boundaries, 2)
# plot(regions, title="Regions")
# plot(boundaries, title="Boundaries")

# define a submesh, composed of region 0 and 1
new_region = MeshFunction('size_t', mesh, mesh.topology().dim())
new_region.set_all(0)
new_region.array()[regions.array() == 0] = 1
submesh = SubMesh(mesh, new_region, 1)

# # define new meshfunctions on this submesh with the same values as their original mesh
submesh_regions = MeshFunction('size_t', submesh, submesh.topology().dim())
submesh_regions.set_all(0)
submesh_boundaries = MeshFunction('size_t', submesh, submesh.topology().dim() - 1)
submesh_boundaries.set_all(0)

# plot(submesh_regions, title="Submesh regions")
# plot(submesh_boundaries, title="Submesh boundaries")

bottom.mark(submesh_boundaries, 1)
top.mark(submesh_boundaries, 2)
File("submesh_boundaries.pvd") << submesh_boundaries
ds = Measure("ds")[submesh_boundaries]

Vsub = FunctionSpace(submesh, 'P', 1)
u = TrialFunction(Vsub)
v = TestFunction(Vsub)
u = Function(Vsub)
dx = Measure("dx")[submesh_regions]

bcs = []


F = kappa*inner(nabla_grad(u), nabla_grad(v))*dx + g*v*ds(2) + r*sinh(C*u)*v*ds(1) - f*v*dx

J = derivative(F, u)
problem = NonlinearVariationalProblem(F, u, bcs, J)
solver  = NonlinearVariationalSolver(problem)
solver.solve()

File("u.pvd").write(u)
2 Likes