Expression with cpp code behaviour changed with 2018.1.0.dev0/pybind11


362
views
0
5 months ago by
Sa Wu  
With my older FEniCS code i have segments similar to
from dolfin import *

mesh = UnitSquareMesh(10,10)
basedim = mesh.geometry().dim()

cpp_string = """
class MyCppExpression : public Expression {{
    public:
        double xx, yy, {} rr;

        MyCppExpression(): Expression(), xx(0), yy(0), {} rr(1) {{}}

        void eval(Array<double>& values, const Array<double>& arg) const {{
            double dx = (arg[0]-xx)/rr, dy = (arg[1]-yy)/rr{};
            values[0] = dx*dx+dy*dy;
        }}
}};
""".format(*(' zz, ', ' zz(0), ', ', dz = (arg[2]-zz)/rr') if basedim == 3 else ('', '', ''))
exp = Expression(cpp_string, degree = 1)​

With fenics from git, 2018.1.0.dev0 with the pybind11 interface this produces the following for jit, which then obviously fails in compilation:

namespace dolfin
{
  class dolfin_expression_7b17744efe1d7d5f390ad591e52004e2 : public Expression
  {
     public:
       

       dolfin_expression_7b17744efe1d7d5f390ad591e52004e2()
       {
            
       }

       void eval(Eigen::Ref<Eigen::VectorXd> values, Eigen::Ref<const Eigen::VectorXd> x) const override
       {
          values[0] = 
class MyCppExpression : public Expression {
    public:
        double xx, yy,  rr;

        MyCppExpression(): Expression(), xx(0), yy(0),  rr(1) {}

        void eval(Array<double>& values, const Array<double>& arg) const {
            double dx = (arg[0]-xx)/rr, dy = (arg[1]-yy)/rr;
            values[0] = dx*dx+dy*dy;
        }
};
;

...

Is constructor that I'd need, with the cpp code defining the class not just the scalar for values[0], not part of the pybind11 bindings?

Or did just the syntax, datatypes, etc change?

Community: FEniCS Project

2 Answers


4
5 months ago by
I saw a reference to this question on the FEniCS slack discussion board and I'm posting a simpler solution which works for 90% of the cases where you would want a multi-line multi-statement C++ expression.

It works by creating a C++ lambda function and then immediately evaluating it to get the expression value. It cannot store custom class attributes etc like the full blown solution of Michal, but it works in many cases:

from dolfin import *


mesh = UnitSquareMesh(10, 10)
V = FunctionSpace(mesh, 'CG', 1)


# Create a standard expression and function
e1 = Expression('x[0]', degree=1)
u = interpolate(e1, V)


# Create an multi-line expression that uses the value of the
# function u by use of a C++ lambda
cpp_code = """
[&]() {
    if (x[0] < 0.5) {
        return u;
    } else {
        return 0.5;
    }
}()
"""
e2 = Expression(cpp_code, u=u, degree=1)


# Plot the results
import matplotlib
matplotlib.use('Agg')
from matplotlib import pyplot

c = plot(e2, mesh=mesh)
pyplot.colorbar(c)
pyplot.savefig('test.png')​


The resulting plot looks like this:

Thanks. I only need multiline cpp Code for some polar and other coordinate transformations before plugging it into some sympy generated code. So, this covers all cases I need. Also, this approach works for lists of expressions.
written 11 weeks ago by Sa Wu  
2
5 months ago by

This is the new C++ Expression compilation in a nutshell.
Essentially, you just subclass C++ dolfin::Expression in a standard way and include pybind11 bindings.
In the example below I exposed just the constructor method. For more advanced stuff just combine proper Expression subclassing and proper pybind11 bindings.

import dolfin as d

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

#include <dolfin/function/Expression.h>
#include <dolfin/function/Constant.h>
class MyCppExpression : public dolfin::Expression
{
public:
  MyCppExpression() : dolfin::Expression() {}

  void eval(Eigen::Ref<Eigen::VectorXd> values, Eigen::Ref<const Eigen::VectorXd> x) const
  {
    values[0] = 1.0;
  }
  };

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

'''

f = d.CompiledExpression(d.compile_cpp_code(code).MyCppExpression(), degree=1)

Thanks a lot Michal.

Any hints/clues/ideas how to do this for a list of expressions? I used to generate lists of precompiled cpp Expressions with sympy. With the proposed way I don't know how to replace

f = d.CompiledExpression(d.compile_cpp_code(code).MyCppExpression(), degree=1)

such that it works with a list of "code". Kind of obviously just directly building a loop around this yields

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

when dijitso loads the second expression.

I'll try to make an abridged example later. Thanks a lot.

written 4 months ago by Sa Wu  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »