Hello,
I’m trying to solve a poisson equation with non-linear boundary conditions in one subdomain of a domain.
First, I tried solving the poisson equation in the full domain to check if it is working. The code works fine.
But when i try to solve it only in a subdomain using a submesh, i’m not able to get it to work. There is even no change in the geometry of the domains in the first and second case. I’m making some mistake while integrating the poisson solver with the submesh functionality.
I have attached the working code in a single domain first, and then the code to solve in in only one subdomain (not working).
Thanks a lot for taking time to look into this !
Solving on a whole domain (Working)
from __future__ import print_function
from dolfin import *
import matplotlib.pyplot as plt
# from __future__ import print_function
from fenics import *
from mshr import *
import numpy as np
from ufl import nabla_div
from ufl import sinh
import random
domain = Rectangle(Point(0,0), Point(100e-6, 200e-6))
mesh = generate_mesh(domain,100)
V = FunctionSpace(mesh, 'P', 1)
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
class BottomBoundary(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14
return on_boundary and abs(x[1]) < tol
class TopBoundary(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14
return on_boundary and abs(x[1] - 200e-6) < tol
bottom_boundary = BottomBoundary()
top_boundary = TopBoundary()
boundary_parts = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundary_parts.set_all(0)
bottom_boundary.mark(boundary_parts, 1)
top_boundary.mark(boundary_parts, 2)
ds = Measure("ds")[boundary_parts]
u = TrialFunction(V)
v = TestFunction(V)
u = Function(V)
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()
import matplotlib.pyplot as plt
p = plot(u, mode='color')
plt.colorbar(p)
plt.title(r"$\Phi$",fontsize=26)
plt.show()
Solving on a submesh (Not Working)
from __future__ import print_function
from dolfin import *
import matplotlib.pyplot as plt
# from __future__ import print_function
from fenics import *
from mshr import *
import numpy as np
from ufl import nabla_div
from ufl import sinh
import random
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
# define regions and boundaries
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)
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()
import matplotlib.pyplot as plt
p = plot(u, mode='color')
plt.colorbar(p)
plt.title(r"$\Phi$",fontsize=26)
plt.show()