# Very costly interpolation of UserExpression

**URL:** <https://fenicsproject.discourse.group/t/very-costly-interpolation-of-userexpression/1987>\
**Category:** Uncategorized\
**Created:** [November 26, 2019, 11:09am UTC](https://fenicsproject.discourse.group/t/very-costly-interpolation-of-userexpression/1987 "2019-11-26T11:09:17Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![malte\_winckler](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/malte_winckler/32/643_2.png) [@malte\_winckler](https://fenicsproject.discourse.group/u/malte_winckler)\
**Post date:** [November 26, 2019, 11:09am UTC](https://fenicsproject.discourse.group/t/very-costly-interpolation-of-userexpression/1987/1 "2019-11-26T11:09:17Z")

</div>

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))
```

---

<div class="post-metadata">

**Author:** ![plugged](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/plugged/32/5115_2.png) [@plugged](https://fenicsproject.discourse.group/u/plugged)\
**Post date:** [November 26, 2019, 12:49pm UTC](https://fenicsproject.discourse.group/t/very-costly-interpolation-of-userexpression/1987/2 "2019-11-26T12:49:26Z")

</div>

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

> [@How to pass numpy array to compiled expression member variable?](https://fenicsproject.discourse.group/t/how-to-pass-numpy-array-to-compiled-expression-member-variable/963/9):
>
> Thanks Kamensky! Just to add: Using the Eigen::Ref\<\> template instead of py::array\_t\<\> , pybind11 casts numpy arrays directly and also gives access to the Eigen library methods. i.e Change npArray definition to typedef Eigen::Ref\<Eigen::VectorXd\> npArray; arr declaration to npArray arr; and the constructor method to test\_exp(npArray a) : dolfin::Expression(), arr(a) {} This can help in using Eigen Methods e.g arr.sum() in the eval method.

> [@Pass Function to C++ expression](https://fenicsproject.discourse.group/t/pass-function-to-c-expression/1081/3):
>
> Thanks a lot, it works

> [@CompiledSubdomain using C++ class](https://fenicsproject.discourse.group/t/compiledsubdomain-using-c-class/918):
>
> 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 \<pybin…

---

<div class="post-metadata">

**Author:** ![malte\_winckler](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/malte_winckler/32/643_2.png) [@malte\_winckler](https://fenicsproject.discourse.group/u/malte_winckler)\
**Post date:** [November 27, 2019, 6:03pm UTC](https://fenicsproject.discourse.group/t/very-costly-interpolation-of-userexpression/1987/3 "2019-11-27T18:03:32Z")

</div>

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
```

---

<div class="post-metadata">

**Author:** ![malte\_winckler](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/malte_winckler/32/643_2.png) [@malte\_winckler](https://fenicsproject.discourse.group/u/malte_winckler)\
**Post date:** [November 28, 2019, 9:05am UTC](https://fenicsproject.discourse.group/t/very-costly-interpolation-of-userexpression/1987/4 "2019-11-28T09:05:52Z")

</div>

I resolved the problem myself. Replace the line

> [@malte\_winckler](#):
>
> expr = dolfin.CompiledExpression(g.CharFunc(jc, omega), degree=1)

with

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

```