# MeshFunction arguments

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!

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.

1 Like

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

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?