Cpp-based Expression

1

Hi,

I am having some trouble with defining an expression via cpp code. The following minimal working example does not run

from fenics import *

mesh = UnitSquareMesh(8, 8)
V0 = FunctionSpace(mesh, 'DG', 0)


funcstring = '''
#include <dolfin/function/Expression.h>

class Function : public dolfin::Expression
{
	void eval(Array<double>& values, constArray <double>& x) const
	{
		if (x[0] >= 0.25 && x[0] < 0.5)
		{
			if (x[1] >= 0.5 && x[1] < 0.625)
			{
				values[0] = 1.0
			}
			else if (x[1] >= 0.625 && x[1] < 0.75)
			{
				values[0] = -1.0
			}
			else
			{
				values[0] = 0.0
			}
		}
		else 
		{
			values[0] = 0.0
		}
	}'''

function = project(Expression(funcstring, degree=0), V0)

This returns

------------------- Start compiler output ------------------------
/tmp/tmpsg02dhul/dolfin_expression_48aa633c4be2d56887d6fc2cebc19a1f.cpp: In member function ‘virtual void dolfin::dolfin_expression_48aa633c4be2d56887d6fc2cebc19a1f::eval(Eigen::Ref<Eigen::Matrix<double, -1, 1> >, Eigen::Ref<const Eigen::Matrix<double, -1, 1> >) const’:
/tmp/tmpsg02dhul/dolfin_expression_48aa633c4be2d56887d6fc2cebc19a1f.cpp:64:1: error: expected primary-expression before ‘class’
class Function : public dolfin::Expression
^~~~~
/tmp/tmpsg02dhul/dolfin_expression_48aa633c4be2d56887d6fc2cebc19a1f.cpp:92:8: error: a function-definition is not allowed here before ‘{’ token
{
^
/tmp/tmpsg02dhul/dolfin_expression_48aa633c4be2d56887d6fc2cebc19a1f.cpp:98:8: error: a function-definition is not allowed here before ‘{’ token
{
^
/tmp/tmpsg02dhul/dolfin_expression_48aa633c4be2d56887d6fc2cebc19a1f.cpp:105:8: error: a function-definition is not allowed here before ‘{’ token
{
^
/tmp/tmpsg02dhul/dolfin_expression_48aa633c4be2d56887d6fc2cebc19a1f.cpp:111:8: error: a function-definition is not allowed here before ‘{’ token
{
^
/tmp/tmpsg02dhul/dolfin_expression_48aa633c4be2d56887d6fc2cebc19a1f.cpp: At global scope:
/tmp/tmpsg02dhul/dolfin_expression_48aa633c4be2d56887d6fc2cebc19a1f.cpp:119:8: error: expected unqualified-id before string constant
extern “C” DLL_EXPORT dolfin::Expression * create_dolfin_expression_48aa633c4be2d56887d6fc2cebc19a1f()
^~~
/tmp/tmpsg02dhul/dolfin_expression_48aa633c4be2d56887d6fc2cebc19a1f.cpp:122:1: error: expected ‘}’ at end of input
}
^

------------------- End compiler output ------------------------
Compilation failed! Sources, command, and errors have been written to: /home/jitfailure-dolfin_expression_48aa633c4be2d56887d6fc2cebc19a1f
Traceback (most recent call last):
File “/usr/lib/python3/dist-packages/dolfin/jit/jit.py”, line 167, in compile_class
mpi_comm=mpi_comm)
File “/usr/lib/python3/dist-packages/dolfin/jit/jit.py”, line 47, in mpi_jit
return local_jit(*args, **kwargs)
File “/usr/lib/python3/dist-packages/dolfin/jit/jit.py”, line 103, in dijitso_jit
return dijitso.jit(*args, **kwargs)
File “/home/philipp/.local/lib/python3.6/site-packages/dijitso/jit.py”, line 217, in jit
% err_info[‘fail_dir’], err_info)
dijitso.jit.DijitsoError: Dijitso JIT compilation failed, see ‘/home/jitfailure-dolfin_expression_48aa633c4be2d56887d6fc2cebc19a1f’ for details

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
File “cppfunction.py”, line 35, in
function = project(Expression(funcstring, degree=0), V0)
File “/usr/lib/python3/dist-packages/dolfin/function/expression.py”, line 400, in init
self._cpp_object = jit.compile_expression(cpp_code, params)
File “/usr/lib/python3/dist-packages/dolfin/function/jit.py”, line 158, in compile_expression
expression = compile_class(cpp_data, mpi_comm=mpi_comm)
File “/usr/lib/python3/dist-packages/dolfin/jit/jit.py”, line 170, in compile_class
raise RuntimeError(“Unable to compile C++ code with dijitso”)
RuntimeError: Unable to compile C++ code with dijitso

2

The difficulty is that the string expected by the Python Expression constructor is something that you should be able to paste in as an R-value to complete values[0] = ... , not a full definition of a C++ Expression subclass. For the conditional logic you have in your class’s eval method, you can compress it all into one line of C++ with the ternary operator, as follows:

from fenics import *

mesh = UnitSquareMesh(8, 8)
V0 = FunctionSpace(mesh, 'DG', 0)

funcstring = """((x[0] >= 0.25 && x[0] < 0.5))?
                 ((x[1] >= 0.5 && x[1] < 0.625)?
                   1.0 : ((x[1] >= 0.625 && x[1] < 0.75)? 
                          -1.0 : 0.0)) : 0.0;"""
function = project(Expression(funcstring, degree=0), V0)

from matplotlib import pyplot as plt
plot(function)
plt.show()

If you really do need to define a C++ Expression subclass, I think the way to go would be something more like what @miguel did with SubDomain here, using pybind11 and compile_cpp_code: