Pass Function to C++ expression

variational formulation
1

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.

2

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)
6 Likes
3

Thanks a lot, it works

4

Hey, David!

It seems like your code is not reproducible in the latest dolfin version. Some import errors appear.

I am currently doing things closely related to this topic: passing python list of Functions to c++ function. However, I get TypeError: incompatible function arguments on a simple example with just one Function:

# P1.ufl
# lies in ./cutils/tests/
element = FiniteElement("Lagrange", triangle, 1)

u = Coefficient(element)

M = u*dx
forms = [M]

# MWE
from dolfin import *

cpp_code = """
#include <dolfin.h>
#include <pybind11/pybind11.h>
#include "P1.h"

using namespace dolfin;
namespace py = pybind11;

double simple(std::shared_ptr<Function> u) {
  auto V = u->function_space();
  P1::Form_M mass_form(V->mesh(), u);
  return assemble(mass_form);
}
  

PYBIND11_MODULE(SIGNATURE, m) {
    m.def("simple", &simple, py::arg("u"));
}
""" 

module = compile_cpp_code(cpp_code, include_dirs=['cutils/tests'])
# ./cutils/tests - is the location of P1.h compiled with `ffc -l dolfin P1.ufl`

mesh = UnitSquareMesh(32, 32)
V = FunctionSpace(mesh, 'P', 1)
u = Function(V)

module.simple(u.cpp_object())

Execution of the last piece of code raises:

TypeError: simple(): incompatible function arguments. The following argument types are supported:
    1. (u: dolfin::Function) -> float

Invoked with: <dolfin.cpp.function.Function object at 0x7f79885703b0>

I tried to follow David’s reply, but it didn’t help me.
Does anyone know how to handle this?

2 Likes
5

Yep, the problem is in the dolfin’s version. In 2018.1.0 it works, in the latest docker image of 2019.1.0, it doesn’t

2 Likes
6

Were you able to figure this out? I am having the same problem.

8

For anyone reading this later: You have to includesome header files :slight_smile:
This worked for me. This code simply reads the values from the function/Expression/MeshFunction and adds the values. Hope this helps anyone :slight_smile: :upside_down_face:

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

#include <dolfin/function/Expression.h> // for expressions
#include <dolfin/mesh/MeshFunction.h> //  for MeshFunctions
#include <dolfin/function/Function.h> // for functions



class your_Function_Name : public dolfin::Expression
{
   public:
       your_Function_Name() : dolfin::Expression(3){} // 3 Dimensional vector as output, use Expression() for scalar
       void eval(Eigen::Ref<Eigen::VectorXd> values, Eigen::Ref<const Eigen::VectorXd> x, const ufc::cell& cell) const override
       {

           func1->eval(values, x);
           const double val1 = values[0];
           func1->eval(values, x);
           const double val2 = values[0];
           values[0] = val1 + val2; // do your stuff here

           exp1->eval(values, x);
           const double val3 = values[0];
           exp2->eval(values, x);
           const double val4 = values[0];
           values[1] = val3 + val4; // do your stuff here

           const uint cell_index = cell.index;
           values[2] = (*MshFunc1)[cell_index] + (*MshFunc2)[cell_index]; // do your stuff here

       }

       std::shared_ptr<dolfin::Function> func1;
       std::shared_ptr<dolfin::Function> func2;
       std::shared_ptr<dolfin::Expression> exp1;
       std::shared_ptr<dolfin::Expression> exp2;
       std::shared_ptr<dolfin::MeshFunction<size_t>> MshFunc1;
       std::shared_ptr<dolfin::MeshFunction<size_t>> MshFunc2;
};
PYBIND11_MODULE(SIGNATURE, m)
{
   py::class_<your_Function_Name, std::shared_ptr<your_Function_Name>, dolfin::Expression>(m, "your_Function_Name")
   .def(py::init<>())
   .def_readwrite("func1", &your_Function_Name::func1)
   .def_readwrite("func2", &your_Function_Name::func2)
   .def_readwrite("exp1", &your_Function_Name::exp1)
   .def_readwrite("exp2", &your_Function_Name::exp2)
   .def_readwrite("MshFunc1", &your_Function_Name::MshFunc1)
   .def_readwrite("MshFunc2", &your_Function_Name::MshFunc2);
}
3 Likes