Hello,
I’m trying to solve this PDE:
\nabla^2 u = u - b, \,\,\, u\in \bf{R} \, ;
in a 2D domain \Omega, with Robin boundary conditions on \partial \Omega:
\left( \nabla u \cdot \hat{n} + c u \right)|_{\partial \Omega} = 0 .
The variational form is thus:
\int_\Omega \left( \nabla u \cdot \nabla v + u v \right) \, \text{d}^2 {\bf{x}} + c \int_{\partial \Omega} u v \, \text{d} {\bf{s}} = b \int_\Omega v \, \text{d}^2 {\bf{x}} .
I know how to solve this for a constant c over the boundary, i.e., using the pre-defined measures:
a = (dot(grad(u), grad(v)) + u*v)*dx + Constant(c)*u*v*ds
L = Constant(b)*v*dx
However, I’m getting nonsense when I try to set different c values for different boundaries. Here’s a MWE:
import dolfin as fem
# simple mesh
mesh = fem.UnitSquareMesh(20, 20)
# horizontal boundary
class HBoundary(fem.SubDomain):
def __init__(self, Y0):
self.Y0 = Y0
super().__init__()
def inside(self, x, on_boundary):
return on_boundary and fem.near(x[1], self.Y0, eps= 1.e-12)
# vertical boundary
class VBoundary(fem.SubDomain):
def __init__(self, X0):
self.X0 = X0
super().__init__()
def inside(self, x, on_boundary):
return on_boundary and fem.near(x[0], self.X0, eps= 1.e-12)
# a MeshFunction over facets
boundfn = fem.MeshFunction('size_t', mesh, mesh.topology().dim()-1)
boundfn.set_all(999) # dummy value for inner facets
# mark boundaries
left, right = VBoundary(0.), VBoundary(1.)
bottom, top = HBoundary(0.), HBoundary(1.)
bottom.mark(boundfn, 0)
right.mark(boundfn, 1)
top.mark(boundfn, 2)
left.mark(boundfn, 3)
# space, trial and test functions
V = fem.FunctionSpace(mesh, 'CG', 1)
u = fem.TrialFunction(V)
v = fem.TestFunction(V)
# redefine the ds measure
ds = fem.Measure('ds', domain= mesh, subdomain_data= boundfn)
# the variational form, volume terms
a = (fem.dot(fem.grad(u), fem.grad(v)) + u*v)*fem.dx
b = fem.Constant(1.); L = b*v*fem.dx
# the boundary terms
c = fem.Constant(1.) # for simplicity; what matters is the following
for i in sorted(list(set(boundfn.array())))[:-1]: # not to include inner facets
a += c*u*v*ds(i)
# if instead of using the previous loop, we use this:
#a += c*u*v*ds
# it works; however this will not allow different c values
# further, setting the boundary terms for all
# but one of the boundaries produces the expected solution
# for c = 0 along the unspecified boundary term (3 here, left boundary)
#a += c*u*v*ds(0)
#a += c*u*v*ds(1)
#a += c*u*v*ds(2)
# compute solution
m = fem.Function(V)
fem.solve(a == L, m)
I am lost here, I will appreciate your input
Cheers,
Miguel