CompiledSubdomain using C++ class

My goal is to speed up the evaluation of the SubDomain Gamma

    from dolfin import *
    mesh =     UnitSquareMesh(50,50, 'left/right')
    PHI = FunctionSpace(mesh, 'CG', 1)
    phi_expr = Expression('2.0*(1.0 - x[0])*x[0] + 0.2 - x[1]', degree=2)
    phi = interpolate(phi_expr , PHI)
    File("phi_initial.pvd") << phi

    class Gamma(SubDomain):
        def inside(self, x, on_boundary):
            return phi(x) > 0.5

with a C++ expression as follows

omega_out_code = """

#include <pybind11/pybind11.h>
#include <pybind11/eigen.h>
#include <math.h>
namespace py = pybind11;

#include <dolfin/function/Function.h>
#include <dolfin/mesh/SubDomain.h>

class OmegaOut : public dolfin::SubDomain

  // Create expression with 3 components
  OmegaOut() : dolfin::SubDomain() {}

  // Function for evaluating expression on each cell
  bool inside(Eigen::Ref<const Eigen::VectorXd> x, bool on_boundary) const override
    Eigen::Array<double, 1, 1> values;
    phi.eval(values, x);
    return values[0] > 0.5;
  dolfin::Function phi;


  py::class_<OmegaOut, std::shared_ptr<OmegaOut>, dolfin::SubDomain>
    (m, "OmegaOut")
    .def_readwrite("phi", &OmegaOut::phi);


c = CompiledSubDomain(compile_cpp_code(omega_out_code).OmegaOut(), phi=phi, degree=0)

but I get the error:

Traceback (most recent call last):
  File "", line 140, in <module>
  File "", line 57, in main
    c = CompiledSubDomain(compile_cpp_code(omega_out_code).OmegaOut(), phi=phi, degree=0)
  File "/usr/local/lib/python3.6/dist-packages/dolfin/mesh/", line 125, in __new__
    return compile_subdomain(inside_code, properties)
  File "/usr/local/lib/python3.6/dist-packages/dolfin/mesh/", line 118, in compile_subdomain
    subdomain = compile_class(cpp_data, mpi_comm=mpi_comm)
  File "/usr/local/lib/python3.6/dist-packages/dolfin/jit/", line 145, in compile_class
    raise RuntimeError("Expression must be a string, or a list or tuple of strings")
RuntimeError: Expression must be a string, or a list or tuple of strings

Is this not supported for SubDomain yet?

The following way of defining c seems to work:

c = compile_cpp_code(omega_out_code).OmegaOut()
c.phi = phi.cpp_object()

# (Change phi_expr definition to just "x[0]" for simpler testing.)
1 Like