When solving a PDE on a submesh $\Omega_0\subset\Omega$ of a given mesh $\Omega$, how do I apply boundary conditions to the PDE?
I'm able to correctly extract the submesh and set up the FunctionSpaces. I suppose it's something about the Measures that I'm missing here.
The following code tries to setup the solution of Poisson's equation on $\Omega_1 = {(x,y)\in(0,1)^2: x<0.5}\subset \Omega = (0,1)^2$. On the left boundary of $\Omega_1$, I'd like to enforce $n\cdot\nabla u = \alpha (u_0 - u)$ with some given $\alpha, u_0\in\mathbb{R}$.
Dolfin solves something, but the suggested solution just looks funky.
from dolfin import *
mesh = UnitSquareMesh(20, 20, 'left/right')
# -------------------------------------------------
class OmegaLeft(SubDomain):
def inside(self, x, on_boundary):
return x[0] <= 0.5
omegaL = OmegaLeft()
# Initialize mesh function for interior domains
subdomains = CellFunction('size_t', mesh)
subdomains.set_all(0)
omegaL.mark(subdomains, 1)
dx = Measure('dx')[subdomains]
submesh_omegaL = SubMesh(mesh, subdomains, 1)
# -------------------------------------------------
# define boundaries for the subdomain
class LeftBoundary(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and abs(x[0]) < DOLFIN_EPS
left_boundary = LeftBoundary()
omegaL_boundaries = FacetFunction('size_t', submesh_omegaL)
omegaL_boundaries.set_all(0)
left_boundary.mark(omegaL_boundaries, 1)
ds_omegaL = Measure('ds')[omegaL_boundaries]
# -------------------------------------------------
# Solve Poisson's equation.
f = Constant(0.0)
V = FunctionSpace(submesh_omegaL, 'CG', 1)
u = TrialFunction(V)
v = TestFunction(V)
# Boundary conditions of the form
# n.grad(u) = alpha * (u0 - u)
alpha = 1.0
u0 = 0.421
F = dot(grad(u), grad(v)) * dx(1) \
- alpha * (u0-u) * v * ds_omegaL(1) \
- f * v * dx(1)
a = lhs(F)
L = rhs(F)
sol = Function(V)
solve(a == L, sol)
plot(sol)
interactive()
# -------------------------------------------------
Even simpler, when only applying Dirichlet boundary conditions, e.g.,
from dolfin import *
mesh = UnitSquareMesh(4, 4, 'left/right')
# -------------------------------------------------
class OmegaLeft(SubDomain):
def inside(self, x, on_boundary):
return x[0] <= 0.5
omegaL = OmegaLeft()
# Initialize mesh function for interior domains
subdomains = CellFunction('size_t', mesh)
subdomains.set_all(0)
omegaL.mark(subdomains, 1)
dx = Measure('dx')[subdomains]
submesh_omegaL = SubMesh(mesh, subdomains, 1)
# -------------------------------------------------
# define boundaries for the subdomain
class LeftBoundary(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and abs(x[0]) < DOLFIN_EPS
left_boundary = LeftBoundary()
# -------------------------------------------------
# Solve Poisson's equation.
f = Constant(0.0)
V = FunctionSpace(submesh_omegaL, 'CG', 1)
u = TrialFunction(V)
v = TestFunction(V)
F = dot(grad(u), grad(v)) * dx(1) \
- f * v * dx(1)
bcs = DirichletBC(V, 0.421, 'on_boundary')
a = lhs(F)
L = rhs(F)
sol = Function(V)
solve(a == L, sol, bcs=bcs)
plot(sol)
interactive()
# -------------------------------------------------
one gets
*** Error: Unable to set given rows to identity matrix.
*** Reason: some diagonal elements not preallocated (try assembler option keep_diagonal).
*** Where: This error was encountered inside PETScMatrix.cpp.
