How to build VectorFunctionSpace on part of a 2D mesh boundary?

Hi everyone,

I would like to build a VectorFunctionSpace on a part (only) of the boundary of a 2D mesh. I tried the following code:

mesh = UnitSquareMesh(20, 20)
bmesh = BoundaryMesh(mesh, “exterior”)

class Top(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], 20) and on_boundary

bmesh_top = SubMesh(bmesh, Top())
VectorSpace_Boundaries_Top = VectorFunctionSpace(bmesh_top, “Lagrange”, 1)

But got this RuntimeError at last line:
“*** Error: Unable to create mesh entity.
*** Reason: Mesh entity index 0 out of range [0, 0] for entity of dimension 1.”

Could you explain me why?
(Running VectorSpace_Boundaries = VectorFunctionSpace(bmesh, “Lagrange”, 1) works)

Many thanks,
Anne

The error indicates you are trying to create a function space on a mesh with zero elements.

UnitSquareMesh(n, m) creates an n \times m mesh of the unit square, i.e. the square extending from [0, 1] in the x direction and from [0, 1] in the y direction. Therefore, there are no facets satisfying your Top condition:

This is why you receive the error. You should change the above to

class Top(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 1.0) and on_boundary
1 Like

Thank you for your answer.

Actually, I also ran the code with a rectangle mesh:
mesh = RectangleMesh(Point(0.0, 0.0), Point(1.0, 0.5), 20, 20, “right/left”)

and with top boundary defined as follow:
class Top(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], 0.5) and on_boundary

and aslo got the same error. Could you tell me where it could come from?

Many thanks,
Anne

You need to remove the on_boundary from the test condition. For a given mesh, the boundary consists of all facets that are connected to a single cell. With respect to the original mesh, cells are triangles, and facets are lines. However, with respect to the bmesh, cells are lines and facets are points. Since the bmesh is composed of the contour bounding the rectangle (a closed contour), all facets (points) are connected to two cells (lines), and thus there are no facets (points) that are considered to be on the boundary. Hence your Top subdomain is empty. The following works for me:

from dolfin import *

mesh = RectangleMesh(Point(0.0, 0.0), Point(1.0, 0.5), 20, 20, "right/left")
bmesh = BoundaryMesh(mesh, "exterior")

class Top(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 0.5)
bmesh_top = SubMesh(bmesh, Top())

VectorSpace_Boundaries_Top = VectorFunctionSpace(bmesh_top, "Lagrange", 1)
1 Like

Many thanks conpierce8, that works for me.

1 Like