Very costly interpolation of UserExpression

It seems to me that the interpolation of UserExpression is very costly as soon as if-statements come into play. For instance, I am using the following implementation for a characteristic function of some set of class dolfin.SubDomain:

class G(dolfin.UserExpression):    
    def __init__(self, SubDomain, degree = 2, **kwargs):
        try:    super(G, self).__init__(**kwargs)
        except: pass
    
        self.degree = degree
        self.SubDomain = SubDomain
    
    def eval(self, value, x): 
        if self.SubDomain.inside(x, True):
            value[0] = 1.0
        else:
            value[0] = 0.0
    
    def value_shape(self):
        return ()

Is there a better way to implement a characteristic function such that the interpolation into a dolfin.FunctionSpace is not so costly? It is important for me that I obtain a dolfin.Function. Please find attached a small MWE which compares the interpolation with and without the if-statement. Any help is very much appreciated.

import dolfin
import time

class G(dolfin.UserExpression):    
    def __init__(self, SubDomain, degree = 2, **kwargs):
        try:    super(G, self).__init__(**kwargs)
        except: pass
    
        self.degree = degree
        self.SubDomain = SubDomain
    
    def eval(self, value, x): 
        if self.SubDomain.inside(x, True):
            value[0] = 1.0
        else:
            value[0] = 0.0
    
    def value_shape(self):
        return ()

class G_fast(dolfin.UserExpression):    
    def __init__(self, SubDomain, degree = 2, **kwargs):
        try:    super(G_fast, self).__init__(**kwargs)
        except: pass
    
        self.degree = degree
        self.SubDomain = SubDomain
    
    def eval(self, value, x): 
        value[0] = 1.0
    
    def value_shape(self):
        return ()

class Omega(dolfin.SubDomain):
    def inside(self, x, on_boundary):
        radius = dolfin.sqrt(x[0]**2 + x[1]**2)
        return radius <= 0.2


if __name__ == "__main__":

    nx = ny = 100

    mesh = dolfin.UnitSquareMesh(nx,ny)

    V = dolfin.FunctionSpace(mesh, "DG", 0)

    omega   = Omega()
    g       = G(omega)
    g_fast  = G_fast(omega)

    start_1 = time.process_time()
    g_inter = dolfin.interpolate(g,V)
    end_1   = time.process_time()

    start_2 = time.process_time()
    g_fast_inter = dolfin.interpolate(g_fast, V)
    end_2   = time.process_time()

    print('slow: {:5.10f}s'.format(end_1 - start_1))
    print('fast: {:5.10f}s'.format(end_2 - start_2))

Maybe take a look at these, and try to implement your code in C++, this will make it much faster than using a UserExpression, that uses python to do the computations

1 Like

Thank you for the links. It really helped me a lot! However, the code I put together still gives me an Error that I cannot resolve. I guess, the problem is the pybind for cpp_code_charFunc. Since I get the following error:

ImportError: generic_type: type "CharFunc" is already registered!

Please find below the code. Any further help is very much appreciated.

import dolfin
import numpy as np

cpp_code_charFunc = """

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

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

class CharFunc : public dolfin::Expression
{
public:
    
    std::shared_ptr<dolfin::Constant> jc;
    std::shared_ptr<dolfin::SubDomain> omega;

    CharFunc(std::shared_ptr<dolfin::Constant> jc_, std::shared_ptr<dolfin::SubDomain> omega_) : dolfin::Expression(){
        jc = jc_;
        omega = omega_;
    }

    void eval(Eigen::Ref<Eigen::VectorXd> values, Eigen::Ref<const Eigen::VectorXd> x) const override
    {
          if (omega->inside(x, 1>0)){
              values[0] = jc->values()[0];
          }
          else{
              values[0] = 0.0;
          }
    }

};

PYBIND11_MODULE(SIGNATURE, m)
{
  py::class_<CharFunc, std::shared_ptr<CharFunc>, dolfin::Expression>(m, "CharFunc")
    .def(py::init<std::shared_ptr<dolfin::Constant>, std::shared_ptr<dolfin::SubDomain>  >())
    .def_readwrite("jc", &CharFunc::jc)
    .def_readwrite("omega", &CharFunc::omega);
}

"""

cpp_code_Omega = """

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

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

class Omega : public dolfin::SubDomain
{
public:

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

  // Function for evaluating expression on each cell
  bool inside(Eigen::Ref<const Eigen::VectorXd> x, bool on_boundary) const override
  {
    return sqrt( pow(x[0], 2) + pow(x[1], 2)) < 0.2;
  }

};

PYBIND11_MODULE(SIGNATURE, m)
{
  py::class_<Omega, std::shared_ptr<Omega>, dolfin::SubDomain>
    (m, "Omega")
    .def(py::init<>());
}

"""

mesh = dolfin.UnitSquareMesh(10,10)
V = dolfin.FunctionSpace(mesh,"DG",0)

omega = dolfin.compile_cpp_code(cpp_code_Omega).Omega()

#omega = Omega()
jc = dolfin.Constant(2.0)
g = dolfin.compile_cpp_code(cpp_code_charFunc)
expr = dolfin.CompiledExpression(g.CharFunc(jc, omega), degree=1)

expr_inter = dolfin.interpolate(expr, V)

dolfin.File("expr.pvd") << expr_inter

I resolved the problem myself. Replace the line

with

expr = dolfin.CompiledExpression(g.CharFunc(jc.cpp_object(), omega), degree=1)