solving PDE on submesh

Question

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.

Comments

This error message

*** 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.

actually means PETSc error: object in wrong state. But due to its common cause it is translated into some diagonal elements not preallocated for DOLFIN user.

---

Both code gives assertion fail

*** Error:   Unable to complete call to function operator[]().
*** Reason:  Assertion &entity.mesh() == _mesh.get() failed.
*** Where:   This error was encountered inside /usr/users/blechta/fenics/fenics/src/dolfin/dolfin/mesh/MeshFunction.h (line 686).

on line solve(a == L, ...). What about dx measure defined on another mesh?

---

Funny; I don't get this error message for either the codes, neither in 1.2.0 nor in the latest dev.
I'm not sure what should be wrong with dx(1). -- It's created as Measure('dx')[subdomains], and the subdomain with index 1 is omegaL.

---

That's probably because my DOLFIN is compiled with CMAKE_BUILD_TYPE:STRING=Developer and yours is with RelWithDebInfo or Release. Hence you have turned off dolfin_assert checks.

Rest in my answer.

Answer 1

Problem is that your measure dx holds different mesh than your form function space and cell indices does not match. With assertion checks turned off in DOLFIN this is not checked. Just try defining your measure dx using CellFunction on correct mesh submesh_omegaL.

This is quite dangerous and should be reported as a bug as wrong user code passed in one your case without error but computing non-sense.

Comments

Bug filed at https://bitbucket.org/fenics-project/dolfin/issue/88/solving-on-submesh-may-lead-to-wrong.

---

Using dx_subdomain(0) from

mydomains = CellFunction('size_t', submesh_omegaL)
mydomains.set_all(0)
dx_subdomain = Measure('dx')[mydomains]

doesn't help. -- The wrong behavior persists.

---

Are you sure? For me both codes with dx_subdomain(0) in place of dx(1) yields constant solution sol==0.421 which is correct.

Both solutions also looks reasonable for f = Constant(1.0) and fulfills all the BCs. In particular u-u0==0.5 on Robin boundary as required by necessary condition for Neumann problem.

---

It works now, thanks.