Pass Function to C++ expression

Hi everyone,

I would like to speed up my code and I have the following problem: I want to create a more complicated Expression based on a finite element function. I can create the functionality I need with an UserExpression, which looks like this:

class profile(fenics.UserExpression):
def __init__(self, u, *arg, **kwargs):
    super().__init__(*arg, **kwargs)
    self.u = u

def eval(self, values, x):
    x_local = np.mod(x[0], 1.5e-3)
    values[0] = 0
    values[1] = self.u(x_local)


def value_shape(self):
    return (2,)

where u is just a solution of a (either 1 or 2 dimensional) PDE computed in fenics.

I am not an Expert for C++, but I thought that I can also generate this result using a C++ Expression, which should be way faster to interpolate. I know that I can use the std::fmod command in order to do the modulus operation I need, but I am not sure how to pass a fenics function to C++ code, and how to evaluate it there.

So far I tried the following:

cpp_code = """

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

#include <dolfin/function/Expression.h>
#include <dolfin/function/Function.h>

class Profile : public dolfin::Expression
{
public:

 Profile() : dolfin::Expression(2) {}

  void eval(Eigen::Ref<Eigen::VectorXd> values, Eigen::Ref<const Eigen::VectorXd> x) const override
  {
    u->eval(values, x);
    const double val = values[0];
    values[0] = 0.0;
    values[1] = val;
  }

  std::shared_ptr<dolfin::Function> u;
};

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

"""

which I can also compile using

profile = fenics.compile_cpp_code(cpp_code)

However, once I try passing the function to the expression, or compiling it, with

expr = fenics.CompiledExpression(profile.Profile(), degree=1, u=u)

I get the following error:

TypeError: (): incompatible function arguments. The following argument types are supported:
1. (self: dolfin_cpp_module_d2657709e7d4b986fae5165a8cc3b358.Profile, arg0: dolfin.cpp.function.Function) -> None

Invoked with: <dolfin_cpp_module_d2657709e7d4b986fae5165a8cc3b358.Profile object at 0x7fe56413bae8>, Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement(‘Lagrange’, interval, 1), dim=1), 140), FiniteElement(‘Lagrange’, interval, 2)), 147)

Any help for this would be really appreciated, thank you very much.

The following modified version appears to work:

from dolfin import *
import numpy as np

cpp_code = """

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

#include <dolfin/function/Expression.h>
#include <dolfin/function/Function.h>

class Profile : public dolfin::Expression
{
public:

  std::shared_ptr<dolfin::Function> u;

  Profile(std::shared_ptr<dolfin::Function> u_) : dolfin::Expression(2){
    u = u_;
  }

  void eval(Eigen::Ref<Eigen::VectorXd> values, Eigen::Ref<const Eigen::VectorXd> x) const override
  {
    u->eval(values, x);
    const double val = values[0];
    values[0] = 0.0;
    values[1] = val;
  }

};

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

"""

mesh = UnitSquareMesh(10,10)
V = FunctionSpace(mesh,"CG",1)
u = Function(V)
profile = compile_cpp_code(cpp_code)
expr = CompiledExpression(profile.Profile(u.cpp_object()), degree=1)

u.interpolate(Constant(7.0))
values = np.zeros(2)
x = np.array([0.5,0.5])
expr.eval(values,x)
print(values)
3 Likes

Thanks a lot, it works