MeshFunction arguments

mesh
1

I don’t understand the meaning/order of the arguments in MeshFunction, for examples this code gives error when uncommenting one of the two lines:

from dolfin import *

Len = 10.0
W = 2.0
H = 2.0

mesh = BoxMesh(Point(0.0, 0.0, 0.0), Point(Len, W, H), 10, 2, 2)
V = VectorFunctionSpace(mesh, "Lagrange", 1)

class Omega0(SubDomain):
    def inside(self, x, on_boundary):
        return True if x[0] <= Len / 2 else False

class Omega1(SubDomain):
    def inside(self, x, on_boundary):
        return True if x[0] > Len / 2 else False

subdomains = MeshFunction("size_t", mesh, mesh.topology().dim(), 0) # GIVES ERROR
#subdomains = MeshFunction("size_t", mesh, 0) # THIS WORKS

subdomain0 = Omega0().mark(subdomains, 0)
subdomain1 = Omega1().mark(subdomains, 1)

V1 = FunctionSpace(mesh, "CG", 1)                
u = Function(V1)

v2d = vertex_to_dof_map(V1)
u.vector()[v2d] = subdomains.array()

#from vtkplotter.dolfin import plot
#plot(subdomains)

Traceback (most recent call last):
File “test.py”, line 28, in
u.vector()[v2d] = subdomains.array()
IndexError: Index mismatch

How do I set correctly the u.vector() in the first case when I get the error message? Thanks!

2

The third argument to the MeshFunction constructor is the dimension of the mesh entities for which values of type "size_t" are stored. The fourth, optional, argument is the value to which the MeshFunction is initialized. (By default, it is zero.) The reason the first version gives an error is because u has values at vertices, whereas the MeshFunction has values at entities of the topological dimension of the mesh, i.e., cells. There is no one-to-one correspondence between cells and vertices, which is why there is an index mismatch. The second version works because the MeshFunction associates values with entities of dimension zero, i.e., vertices.

2 Likes
3

Thanks for your answer, so what would be the way to map the values from the entities to the vertices?

4

You can read The FEniCS Book on chapter 16. Entities of dimension 0 have the same order with vertices. I think you got the correct result already if using the commented line.
I have two other ways to do the same job:

subdomains = MeshFunction("size_t", mesh, 0, 0)
# subdomain0 = Omega0().mark(subdomains, 0)
subdomain1 = Omega1().mark(subdomains, 1)

subdomains = MeshFunction("size_t", mesh, 0, 1)
subdomain0 = Omega0().mark(subdomains, 0)
# subdomain1 = Omega1().mark(subdomains, 1)

Is it interesting?

1 Like