Defining a boundary throughout a 3D subdomain rather than on a 2D facet.

Question

Hi,

I regularly use a facet_region.xml file when defining boundary conditions on a 3D mesh. E.g.

mesh = Mesh('mesh.xml')
subdomains = MeshFunction('size_t',mesh,'mesh_physical_region.xml')
boundaries = MeshFunction('size_t',mesh,'mesh_facet_region.xml')
V = FunctionSpace(mesh,'Lagrange',1)

bc = DirichletBC(V,Constant(0.0),boundaries,1)

Instead of using the facet_region file, is it possible to set a boundary throughout the entire volume of a subdomain? I naively tried

bc = DirichletBC(V,Constant(0.0),subdomains,1)

and

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

But I got the error message "User MeshFunction is not a facet MeshFunction (dimension is wrong)" and "Found no facets matching domain for boundary condition" respectively.

Thanks

Comments

Hi, you have to manually mark the facets of cells where you want to set the boundary condition, e.g.

from dolfin import *

mesh = UnitSquareMesh(20, 20)
h = mesh.hmin()

sub_cells = CellFunction('size_t', mesh, 0)
AutoSubDomain(lambda x: x[0] < 0.5 + h).mark(sub_cells, 1)

sub_facets = FacetFunction('size_t', mesh, 0)
mesh.init(2, 1)
for cell in SubsetIterator(sub_cells, 1):
    for facet in facets(cell): 
        sub_facets[facet] = 1

V = FunctionSpace(mesh, 'CG', 1)
f = interpolate(Expression('x[0]*x[0]+x[1]*x[1]+1'), V)
bc = DirichletBC(V, Constant(1), sub_facets, 1)
bc.apply(f.vector())

plot(f)
interactive()