Complex userExpression using C++ class

Errors
1

Hi all,
I have a question how to write C++ class for a userExpression that takes the solution u as an argument. Thanks to @dokken (How to difine UserExpression that takes the solution u), I have a Python code (Version ‘2019.1.0’) similar to this.

N = 100
X = 1.0
mesh = dl.IntervalMesh(N, 0, X)
degree = 1 # linear basis functions
Vh = dl.FunctionSpace(mesh, "Lagrange", degree)
u = dl.interpolate(dl.Expression("x[0]", degree = 1), Vh)
alpha = 2.0

class user_expr(dl.UserExpression):
    def set_alpha(self, u, alpha): # in the user_expr class
        self.u = u
        self.alpha = alpha

    def eval(self, value, x): # in the user_expr class
        if self.u(x) > 0.5:
            value[0] = 0.5
        else:  
            value[0] = self.alpha * self.u(x)* self.u(x)
        
    def value_shape(self):
        return ()

user_expr2 = user_expr(degree = 0)
user_expr2.set_alpha(u, alpha)

q = dl.interpolate(user_expr2, dl.FunctionSpace(mesh, "Lagrange", 1))
d_to_v = dl.dof_to_vertex_map(q.function_space())
print(q.vector()[:][d_to_v])

, and it works fine. I want to speed up this code by rewriting it into C++ codes. Regardless of my search, I could not make it work due to C++ compilation errors etc. the current C++ is like this.

cpp_code = """

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

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

class user_expr : public dolfin::Expression
{
public:

  std::shared_ptr<dolfin::Function> u;
  double alpha;
  
  user_expr(std::shared_ptr<dolfin::Function> u_) : dolfin::Expression(){
    u = u_;
  }

  void eval(Eigen::Ref<Eigen::VectorXd> value, Eigen::Ref<const Eigen::VectorXd> x) const override
  { 
    u->eval(value, x);
    if u > 0.5
        value[0] = 0.5
    else
        value[0] = alpha * u * u
  }

};

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

"""

with

profile = dl.compile_cpp_code(cpp_code)
expr = dl.CompiledExpression(profile.user_expr(u.cpp_object()), degree=1)
expr.alpha = alpha

Could anyone help me? I appreciate any help.

2

It would help if you posted your compiler error.

However the following looks incorrect

u->eval(value, x);
if u > 0.5

As u is a dolfin::Function, you cannot compare it with a floating point litteral. Also, the C++ syntax is incorrect.
What you probably meant was

u->eval(value, x);
if (value[0] > 0.5)
    value[0] = 0.5;
else
    value[0] = alpha * value[0] * value[0];

(Not entierly sure the above actually works, I have not checked the syntax or tried to run it)

3

You also do not initialize alpha, so you will probably need to change the constructor to

user_expr(std::shared_ptr<dolfin::Function> u_, double alpha_) : dolfin::Expression(){
    u = u_;
    alpha = alpha_;
  }

and pass in alpha from python to the profile.user_expr function call.

4

Thank you so much for your help! After some try and error, I could compile the cpp code and it works well (although I could not put alpha in the constructor…). The code below works fine. I appreciate your help.

 cpp_code = """

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

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

class user_expr : public dolfin::Expression
{
public:

  std::shared_ptr<dolfin::Function> u;
  double alpha;
  
  user_expr(std::shared_ptr<dolfin::Function> u_) : dolfin::Expression(){
    u = u_;
  }

  void eval(Eigen::Ref<Eigen::VectorXd> value, Eigen::Ref<const Eigen::VectorXd> x) const override
  { 
    u->eval(value, x);
    if (value[0] > 0.5){
        value[0] = 0.5;
    }
    else {
        double val = value[0];
        value[0] = alpha * val * val;
    }
  }
};

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

"""

with

profile = dl.compile_cpp_code(cpp_code)
expr = dl.CompiledExpression(profile.user_expr(u.cpp_object()), degree=1)
expr.alpha = 2.0
r = dl.interpolate(expr, dl.FunctionSpace(mesh, "Lagrange", 1))
dl.plot(r)