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()