Do I need to create an `on_boundary` function when using a mesh

General
1

Hello all, another noobie question I am sorry. I had a question about setting boundary conditions when I import a mesh file. I am using Fenics 2019.1.

In the SOLVING PDES IN PYTHON book p. 23, the authors indicate that a user should specify a boundary function, to indicate whether some coordinate is on or off the boundary. The boundary function then goes in the arguments for the DirichletBC function. That boundary function
can be very simple like:

def boundary(x, on_boundary):
    return on_boundary

Now, I was watching a video tutorial for using Meshes in Fenics, and in that presenter’s code he imports the mesh but does not use a boundary function when setting the DirichletBC. In this case, he is showing the transient heat equation on a mesh. The code for that tutorial is below–taken from the linked source.

My question is, should there be a boundary function in the code below? I am not clear on what the rule is. Like if there is a mesh, then should I never include a boundary function? Or if there is a mesh, then the boundary function is optional, and it should look like …" Is there some specific form for the boundary function when using a mesh? That is my basic question.

from fenics import *

t_end = 10.0
dt = 0.1
k = 300
u_in = 20
u_out = -20

xml_file = "helix.xml"
mesh = Mesh(xml_file)
fd = MeshFunction('size_t', mesh, "helix_facet_region.xml");

V = FunctionSpace(mesh, 'P', 1)

bc1 = DirichletBC(V, Constant(u_in), fd, 3)
bc2 = DirichletBC(V, Constant(u_out), fd, 2)
bc = [bc1, bc2]

u = TrialFunction(V)
v = TestFunction(V)
u_n = Function(V)

F = u*v*dx + dt*k*dot(grad(u), grad(v))*dx - u_n*v*dx
a, L = lhs(F), rhs(F)

u = Function(V)
t = 0
vtkfile = File('output/output.pvd')

num_steps = int(t_end/dt)
for n in range(num_steps):
    t += dt
    solve(a == L, u, bc)
    u_n.assign(u)
    vtkfile << (u, t)

Again the question is, is it good practice to include a boundary function when specifying a DirichletBC with a mesh? If yes, then what would that boundary function look like?
Thanks.

2

The boundary conditions here are set on the DoFs which are contained by the topology marked by the meshfuction.

The method you describe using an “on_boundary” function searches for DoFs geometrically as you state.

Use the right tool for the job. E.g. a simple 2d mesh will easily be searchable using geometry. A complex 3d mesh with many intricate parts will be near impossible to search geometrically and topological definitions should be employed as defined by the mesh generator.

3

@nate thanks for the response. Okay, so sounds like you are saying that the meshfunction already knows how to handle the topology and identify the boundaries. Seems like this is the preferred way, and that I should not use the on_boundary function unless I am using some toy example?

Also, does any of this change now that everyone seems to be moving towards the XDMF file types? The file structure is different–with the h5 backend, but I am not sure if the topology handling is also improved? But I will follow your recommendation.