How to express time-dependent compiled (cpp) expression with FEniCS 2018.1.0.dev0/pybind11?


356
views
1
5 months ago by
I'm trying to create a minimal example of how to define a compiled expression with time represented as a Constant inspired by https://bitbucket.org/fenics-project/dolfin/src/master/python/demo/documented/tensor-weighted-poisson/demo_tensor-weighted-poisson.py?fileviewer=file-view-default and https://www.allanswered.com/post/jrmxq/expression-with-cpp-code-behaviour-changed-with-2018-1-0-dev0pybind11/. The compiled expression code currently looks like this:

f_code = """

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

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

class F : public dolfin::Expression
{
public:

  // Create expression
  F() : dolfin::Expression() {}

  // Function for evaluating expression on each cell
  void eval(Eigen::Ref<Eigen::VectorXd> values, Eigen::Ref<const Eigen::VectorXd> x) const override
  {
    values[0] = 1.0;//time->double();
  }

  // The data stored in Constant
  std::shared_ptr<dolfin::Constant> time;

};

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

"""​


and is invoked with

 time = Constant(0.0)
 f = CompiledExpression(compile_cpp_code(f_code).F(),
                           time=time, degree=0)


which yields this error:

Traceback (most recent call last):

File "/usr/local/lib/python3.5/dist-packages/dolfin/function/expression.py", line 319, in __init__
setattr(self._cpp_object, k, val)
TypeError: (): incompatible function arguments. The following argument types are supported:
1. (self: dolfin_cpp_module_b5bfb8417c3f881dcc37206004268f99.F, arg0: dolfin.cpp.function.Constant) -> None

Invoked with: <dolfin_cpp_module_b5bfb8417c3f881dcc37206004268f99.F object at 0x7fbca7029ed8>, Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 8)

Any tips on how to get this working?

Community: FEniCS Project

3 Answers


2
5 months ago by
Thanks for the answers! Posting working example comparing the two/three different approaches here for future reference. Sample output:

Time elapsed (Python) = 1.1004652976989746
Time elapsed (cpp_code) = 0.0013682842254638672
Time elapsed (CompiledExpression) = 0.0011909008026123047

import time

from dolfin import *

N = 8
mesh = UnitCubeMesh(N, N, N)
delta = 2.0
t = Constant(0.1)
V = FunctionSpace(mesh, 'DG', 1)

# -----------------------------------------------------------------------------
# Version 1: Just Python
# -----------------------------------------------------------------------------
class MyFunc(UserExpression):
    def __init__(self, delta, t, *arg, **kwargs):
        super().__init__(*arg, **kwargs)
        self.delta = delta
        self.t = t
        
    def eval(self, values, x):
        if (abs(x[0]) <= 0.5) and (abs(x[1]) <= 0.5):
            values[0] = self.delta*(sin(2.0*pi*self.t))
        else:
            values[0] = 0.0

    def value_shape(self):
        return ()

# Test and time evaluation of Python Expression
f = MyFunc(delta, t, degree=0)
t_start = time.time()
p = interpolate(f, V)
t_stop = time.time()
print("Time elapsed (Python) = ", (t_stop - t_start))

# -----------------------------------------------------------------------------
# Version 2: with cpp lambda function:
# -----------------------------------------------------------------------------
cpp_code = """
[&]() {
    if (std::abs(x[0]) <= 0.5 && std::abs(x[1]) <= 0.5) {
        return delta*(sin(2.0*pi*t));
    } else {
        return 0.0;
    }
}()
"""

# Test and time evaluation of cpp_code Expression
e = Expression(cpp_code, delta=delta, t=t, degree=0)
t_start = time.time()
u = interpolate(e, V)
t_stop = time.time()
print("Time elapsed (cpp_code) = ", (t_stop - t_start))
assert(u.vector().sum() == p.vector().sum())

# -----------------------------------------------------------------------------
# Version 3: with Compiled Expression function:
# -----------------------------------------------------------------------------

f_code = """

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

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

class F : public dolfin::Expression
{
public:

  // Create expression
  F() : dolfin::Expression() {}

  // Function for evaluating expression on each cell
  void eval(Eigen::Ref<Eigen::VectorXd> values, Eigen::Ref<const Eigen::VectorXd> x) const override
  {
    if (std::abs(x[0]) <= 0.5 && std::abs(x[1]) <= 0.5) {
        values[0] = delta*(sin(2.0*M_PI*t->values()[0]));
    } else {
        values[0] =  0.0;
    }

  }

  // Constant representing time
  std::shared_ptr<dolfin::Constant> t;

  // Other variables
  double delta;

};

PYBIND11_MODULE(SIGNATURE, m)
{
  py::class_<F, std::shared_ptr<F>, dolfin::Expression>
    (m, "F")
    .def(py::init<>())
    .def_readwrite("t", &F::t)
    .def_readwrite("delta", &F::delta);
}
"""
f = CompiledExpression(compile_cpp_code(f_code).F(), t=t.cpp_object(),
                       delta=delta, degree=0)
t_start = time.time()
g = interpolate(f, V)
t_stop = time.time()
print("Time elapsed (CompiledExpression) = ", (t_stop - t_start))
assert(g.vector().sum() == p.vector().sum())
​
1
5 months ago by
As a quick fix invoke the compiled code with
time = Constant(0.0)
f = CompiledExpression(compile_cpp_code(f_code).F(),
                       time=time.cpp_object(), degree=0)​


Explanation:
You created bindings for getters and setters for member variable std::shared_pointer<dolfin::Constant> time. Actually, pythonic Constant() is not pybind11's "version" of dolfin::Constant , there is an additional layer in between, see https://bitbucket.org/fenics-project/dolfin/src/9eba551bd3ed5ab7bb61cd63436d0f23e224a685/python/dolfin/function/constant.py?at=master&fileviewer=file-view-default

Therefore you must call .cpp_object() to get an actual object that pybind11 recognizes as valid input for &F::time.

Thanks Michal, will test!
written 5 months ago by Marie E. Rognes  
But you must use line values[0] = double(*time); The double is an operator in Constant.cpp.

Or better do values[0] = time->values()[0]; and stay dolfinX compatible ;)
written 5 months ago by Michal Habera  
1
5 months ago by
Maybe something like this:
from dolfin import *

f_code = """

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

#include <dolfin/function/Expression.h>

class F : public dolfin::Expression
{
public:

  // Create expression
  F() : dolfin::Expression() {}

  // Function for evaluating expression on each cell
  void eval(Eigen::Ref<Eigen::VectorXd> values, Eigen::Ref<const Eigen::VectorXd> x) const override
  {
    values[0] = t;
  }

  double t;

};

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

"""

F = compile_cpp_code(f_code).F()

f = CompiledExpression(F, degree=0)

F.t = 0.123

print(f(0.0,0.0))​

Not using "Constant". As Michal says, you need to use cpp_object() on some classes.

Thanks Chris, however using the Constant is the key point . In particular, so that the Expression is "automatically" updated when the constant is updated.
written 5 months ago by Marie E. Rognes  
I think, if you can get away with a one-liner "Expression" instead of "CompiledExpression", then some of the wrapping will be done for you, see "python/test/unit/function/test_expression.py". Maybe Tormod's trick with a lambda function would be enough?
written 5 months ago by Chris Richardson  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »