Proper way to turn MeshFunction over cells into a Function?

Question

If I have a MeshFunction over cells, representing (for example) thermal conductivity in the heat equation, how can I convert that into a Function object, and use it in a form?

mesh = Mesh("mymesh.xml.gz")
V = FunctionSpace(mesh, "Lagrange", 1)
u = Function(V)
tcond = MeshFunction("size_t", mesh, "materials.xml.gz")

I was going to try something like this, but the number of nodes != number of cells …

tcond_vals = (10., 100.)
helper = numpy.asarray(tcond.array(), dtype=numpy.int32)
tcond.vector()[:] = numpy.choose(helper, tcond_vals)

Answer 1

Hi, a proper function space for your conductivity could be DG0 where you get one degree
of freedom per cell. Try

V = FunctionSpace(mesh, "DG", 0)
u = Function(V)

dm = V.dofmap()
for cell in cells(mesh):
   helper[dm.cell_dofs(cell.index())] = tcond[cell]

u.vector()[:] = numpy.choose(helper, tcond_vals)

Comments

That worked, thanks!
(going to go back and read some more, so I can ask less silly questions)

---

Have in mind that mesh entities like the cell number does not necessary correspond to the dof number in the vector. So even if the above code does not complain it most probably does not generate the correct result.

You need to map the cell indices through the dofmap from the FunctionSpace. Something like:

dm = V.dofmap()
for cell in cells(mesh):
    helper[dm.cell_dofs(cell.index())] = tcond[cell]
tcond.vector()[:] = helper

---

Hi, you are right of course. I have edited the answer following your comment.