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, iām not able to get it to work. Iām trying to use a submesh for the subdomain. 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 feature.

I have attached the working code in a single domain first, and then the code to solve in in only one subdomain (not working).

Please help me with this issue. Iām not able to resolve it.

Thanks in advance. @bmf @cdaversin

## Working (Single domain)

from **future** import print_function

from dolfin import *

import matplotlib.pyplot as plt

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*vds(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Ā§

plt.title(r"\Phi",fontsize=26)

plt.show()

Solving on subdomain (not working)

from **future** import print_function

from dolfin import *

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, 400e-6))

mesh = generate_mesh(domain,100)

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 OmegaMiddle(SubDomain):

def inside(self, x, on_boundary):

return x[1] <= 200e-6

omegaM = OmegaMiddle()

# Initialize mesh function for interior domains

subdomains = MeshFunction(āsize_tā, mesh, 2)

subdomains.set_all(0)

omegaM.mark(subdomains, 1)

dx = Measure(ādxā)[subdomains]

submesh_omegaL = SubMesh(mesh, subdomains, 1)

# define boundaries for the subdomain

class BottomBoundary(SubDomain):

def inside(self, x, on_boundary):

tol = 1E-14

return on_boundary and abs(x[1]) < tol

Bottom_boundary = BottomBoundary()

class TopBoundary(SubDomain):

def inside(self, x, on_boundary):

tol = 1E-14

return on_boundary and abs(x[1]-200e-6) < tol

Top_boundary = TopBoundary()

omegaL_boundaries = MeshFunction(āsize_tā, submesh_omegaL,1)

omegaL_boundaries.set_all(0)

Bottom_boundary.mark(omegaL_boundaries, 1)

Top_boundary.mark(omegaL_boundaries, 2)

ds_omegaL = Measure(ādsā)[omegaL_boundaries]

mydomains = MeshFunction(āsize_tā, submesh_omegaL,2)

mydomains.set_all(0)

dx_subdomain = Measure(ādxā)[mydomains]

V = FunctionSpace(submesh_omegaL, āPā, 1)

u = TrialFunction(V)

v = TestFunction(V)

u = Function(V)

bcs = []

F = kappa*inner(nabla_grad(u), nabla_grad(v))**dx_subdomain(0) + g*vds_omegaL(2) + r*sinh(C*u)*v*ds_omegaL(1) - f*v*dx_subdomain(0)

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Ā§

plt.show()