Creating an Expression Using a MeshFunction

1

I simulating fluid structure interaction problems, and as part of this I need an expression which is 1 in the solid and zero in the fluid. As part of my mesh creation I already make a mesh function that satisfies this, and I find it useful to be able to simply create an expression from this mesh function, rather than having to try again. I have found a way of doing this by adapting the tensor-weighted-poisson demo, which I show below in a minimal working example,

from dolfin import *
from mshr import *
import matplotlib.pyplot as plt

# create mesh (marking the solid part in the process)
L, H = 2.5, 0.41
l, h = 0.35, 0.02
C, r, SEGMENTS = (0.2,0.2), 0.05, 10
RESOLUTION = 30
domain = Rectangle(Point(0.0,0.0),Point(L,H)) - \
        Circle(Point(C[0],C[1]),r,SEGMENTS)
domain.set_subdomain(1,Rectangle(Point(C[0],C[1]-h/2),Point(C[0]+r+l,C[1]+h/2)))
mesh = generate_mesh(domain, RESOLUTION)

# create mesh function using domains already defined
mf = MeshFunction("size_t", mesh, mesh.geometric_dimension(), mesh.domains())

# c++ code to create an expression using a mesh function
cpp_code = """
#include <pybind11/pybind11.h>
#include <pybind11/eigen.h>
namespace py = pybind11;

#include <dolfin/function/Expression.h>
#include <dolfin/mesh/MeshFunction.h>

class MeshFunExpression : public dolfin::Expression
{
    public:

    // Create scalar expression
    MeshFunExpression() : dolfin::Expression() {}

    // Function for evaluating expression on each cell
    void eval(Eigen::Ref<Eigen::VectorXd> values,
        Eigen::Ref<const Eigen::VectorXd> x, const ufc::cell& cell)
        const override
    {
        const uint cell_index = cell.index;
        values[0] = (*mf)[cell_index];
    }

    // The data stored in mesh functions
    std::shared_ptr<dolfin::MeshFunction<size_t>> mf;

};

PYBIND11_MODULE(SIGNATURE, m)
{
    py::class_<MeshFunExpression, std::shared_ptr<MeshFunExpression>,
        dolfin::Expression> (m, "MeshFunExpression")
    .def(py::init<>())
    .def_readwrite("mf", &MeshFunExpression::mf);
}

"""

# compile the c++ code to make an expression
phi = CompiledExpression( \
        compile_cpp_code(cpp_code).MeshFunExpression(), \
        mf=mf, degree=0)

# plotting
Z = FunctionSpace(mesh,"DG",0)
phi = project(phi,Z)
plot(phi)
plt.savefig("tmp.png")

I would be curious to know if anyone has done this in a neater or better way.