Dirichlet problem on connected subset of mesh

Question

Given a subset of cells (of some mesh) the closed union of which forms a simply connected domain, I'd like to solve a problem with Dirichlet bcs using a restricted function space. However, I don't know how to define the boundary indicator function (SubDomain.inside()) to identity the boundary of my subdomain.
Thanks for advise,
Martin

Comments

A minimal code example would help us help you.

---

Sure, here goes:

from dolfin import *

# create mesh and select patch of some cell -> C
mesh = UnitSquareMesh(4,4)
selected_index = 3
selected_cell = Cell(mesh, selected_index)
C = [c.index() for c in cells(selected_cell)]
C.append(selected_cell.index())

cf = CellFunction('size_t', mesh)
cf.set_all(1)
for i in C:
    cf[i] = 0

# create boundary element measure
dx = dx[cf]

# define restricted space
restriction = Restriction(cf, 0)
V = FunctionSpace(restriction, "CG", 1)

# define functions
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(1.0)

# define forms
a = dot(grad(u),grad(v))*dx()
L = f*v*dx()

# define boundary condition
# BUT: how???
#zero = Constant(0.0)
#boundary = Boundary()
#bc = DirichletBC(V, zero, boundary)

# solve system
u = Function(V)
solve(a == L, u)

# plot cell patch and solution
plot(cf)
plot(u, interactive=True)

Now, this restricted problem should have Dirichlet boundary conditions. How can one define the boundary?

Answer 1

Hi,

If I understand your question correctly, here is a possible solution:

# define boundary condition
# BUT: how???
ff = FacetFunction("size_t", mesh, 0)
mesh.init()
# Set value of ff to 1 on boundary of restriction
for cell in cells(mesh):
    if cf[cell.index()] == 0:
        for facet in facets(cell):
            if facet.num_entities(2) == 1: # exterior boundary
                ff[facet.index()] = 1
            else:
                c0, c1 = facet.entities(2)
                if cf[int(c0)] != cf[int(c1)]:   # interior boundary
                    ff[facet.index()] = 1

zero = Constant(0.0)
bc = DirichletBC(V, zero, ff, 1)

# solve system
u = Function(V)
solve(a == L, u, bcs=bc)

Comments

Exactly what I need. Thanks!