What is the best way to set a subdomain to a given value?

Question

Having a domain that consists of two subdomains (eg. a smaller square inside a larger one), I would like the solution to be zero in the whole 'border' region, rather than just on a boundary.

DirichletBC(V, Constant(0.0), MeshFunction, nr)

does not seem to work when MeshFunction is defined over cells instead of facets. Is there an method to achieve this though? (eg getting all facets within a given region/subdomain or so?)

Thanks

Answer 1

Maybe you can try something like this (this is only one way to do what you are looking)

# Begin demo
from dolfin import *

# Create mesh and define function space
mesh = UnitSquareMesh(32, 32)
V = FunctionSpace(mesh, "Lagrange", 1)

# Exterior domain (PUT ATTENTION HERE!)
class Exterior(SubDomain):
    def inside(self, x, on_boundary):
        return x[0] > 0.75 - DOLFIN_EPS or x[0] < 0.25 + DOLFIN_EPS or x[1] > 0.75 - DOLFIN_EPS or x[1] < 0.25 + DOLFIN_EPS

# Define Dirichlet boundary (x = 0 or x = 1)
def boundary(x):
    return x[0] < DOLFIN_EPS or x[0] > 1.0 - DOLFIN_EPS

# Define boundary condition (AND HERE!)
exterior = Exterior()
u0 = Constant(0.0)
bc1 = DirichletBC(V, u0, boundary)
bc2 = DirichletBC(V, u0, exterior)
bcs = [bc1, bc2]

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)")
g = Expression("sin(5*x[0])")
a = inner(grad(u), grad(v))*dx
L = f*v*dx + g*v*ds

# Compute solution
u = Function(V)
solve(a == L, u, bcs)

# Plot solution
plot(u, interactive=True)

The above code is a modification of the poisson demo.

Comments

Thank you for your answer. I am a bit surprised that this works, but no longer does when I replace the 3 lines related to the boundary condition by

f = MeshFunction('size_t', mesh, mesh.topology().dim())
f.set_all(0)
exterior.mark(f, 1)
bcs = DirichletBC(V, u0, f, 1)

I suppose this is because DirichletBC expects a FacetFunction, not a CellFunction. Since I actually have a complex geometry, I cannot simply define a subdomain as in the example. I came up with some code that 'translates' the cellfunction into a facetfunction. Will post it below as reference. As usual, if such method already existed and I overlooked it, I'd be happy to hear:

# the original example, using a cellfunction to define regions
from dolfin import *
mesh = UnitSquareMesh(32, 32)
V = FunctionSpace(mesh, "Lagrange", 1)
class Exterior(SubDomain):
    def inside(self, x, on_boundary):
        return x[0] > 0.75 - DOLFIN_EPS or x[0] < 0.25 + DOLFIN_EPS or x[1] > 0.75 - DOLFIN_EPS or x[1] < 0.25 + DOLFIN_EPS
exterior = Exterior()
def boundary(x):
    return x[0] < DOLFIN_EPS or x[0] > 1.0 - DOLFIN_EPS
u0 = Constant(0.0)
cellfct = MeshFunction('size_t', mesh, mesh.topology().dim())
cellfct.set_all(0)
exterior.mark(cellfct, 1)

# so suppose we can only use cellfct to define a condition on a region?
dim = mesh.topology().dim()
facetfct = MeshFunction('size_t', mesh, dim - 1)
facetfct.set_all(0)
mesh.init(dim - 1, dim) # relates facets to cells
for f in facets(mesh):
    if 1 in cellfct.array()[f.entities(dim)]: # one of the neighbouring cells of the facet is part of region 1
        facetfct[f.index()] = 1
bcs = [DirichletBC(V, u0, facetfct, 1)]

# the rest of the example
u = TrialFunction(V)
v = TestFunction(V)
f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)")
g = Expression("sin(5*x[0])")
a = inner(grad(u), grad(v))*dx
L = f*v*dx + g*v*ds
u = Function(V)
solve(a == L, u, bcs)
plot(u, interactive=True)

---

your guess in the second paragraph is correct, but I do not know if exists an overloaded method to do what you're looking for

---

How can I achieve the same when the mesh function is defined over vertices instead of cells? http://fenicsproject.org/qa/8968/trouble-applying-function-boundary-definition-dirichletbc