External libraries in JIT compiled expression

Question

I'd like to use the C (!) package cubature in a JIT compiled expression. I can't figure out how to do it. The package documentation reads:

The current version of the code can be downloaded from:

cubature-1.0.2.tgz

a gzipped tar file. This unpacks to a directory containing a README file with instructions and a stand-alone hcubature.c or pcubature.c file (along with a couple of private header files) that you can compile and link into your program for h-adaptive and p-adaptive integration, respectively, and a header file cubature.h that you #include.

In my view, the problem is that this package cannot be "installed".

The package compiles and runs when I use either gcc or cc. I ran a few tests to see if I can call it from another file and all works well. Here is a minimal example of my problem in FEniCS.

The following python file pyfile.py runs well except for a warning (see comment at bottom):

from dolfin import *

loc_file = open( "xpr.cpp" , 'r' )  
code = loc_file.read()
xpr = Expression( code ) 
xpr( (2.,2.) )

and the c++ file is xpr.cpp:

namespace dolfin {

    class Xpr : public Expression
    {
    public:
        Xpr() : Expression() { }

        void eval(Array<double>& values, const Array<double>& x) const
        {
            values[0]  = 1;
        }
    };
}

adding the following line to the beginning of xpr.cpp makes the compilation fail with a long error message.

#include "cubature.h"

I know this is the line I need because this is how I included the package in my tests that didn't involve python and dolfin.

Any help would be greatly appreciated!

UPDATE: I believe what I need to do is give the Expression the option include_dirs = [ "." ] and place the compiled cubature files in the same library as my code. I can't figure this out (yet).

PS 1: The warning mentioned above is Automatic determination of degree for Expressions has been deprecated in FEniCS version 2016.1. I'd love to get rid of it if anyone knows how.

PS 2: Why do I want to use this package? I want to calculate an integral of a function over a computatinal domain. Normally, this can be done via FEniCS. However, the function has a singularity on the boundary, so plan on (adaptively) integrating over every cell, and using that to calculate the domain integral. I think this might help other people who may use FEniCS in thee future. Also, I did try to increase the degree of the expression, but to no avail.

Yair Daon · Answer 1 · Jul 19, 2016

The most useful reference for me was chapter 14 in the FEniCS book (covering JIT compilation using Instant). The book is downloadable from launchpad via this link (hopefully will not be broken anytime soon). Here's the python file I use:

from dolfin import *
import instant 

import numpy as np


def det_A ( v ):
    v[6] =  abs( (v[2] - v[0]) * (v[5] - v[3]) - (v[4] - v[2]) * (v[3] - v[1]) )
tol = 1e-7

xpr_file = open( "xpr.cpp" , 'r' )  
xpr_code = xpr_file.read()
xpr_file.close()

xpr_code = xpr_code.replace( "FUNCTION_NAME", "exam" )

#instant.output.set_log_level("DEBUG")

xpr = instant.build_module(
code=xpr_code,
sources = ["helper.cpp" , "hcubature.c" ],
system_headers=["numpy/arrayobject.h"],
include_dirs=[np.get_include()],
init_code="import_array();",
local_headers = [ "header.h" , "cubature.h" ],
arrays=[[ "m", "vertices" ]]
)

arr = np.array( [
    -3, 0, # a
     6, 0, # b
     0, 3,  # c
     -1 # placeholder for det
     ] , dtype = np.float64 )
det_A( arr )

val = xpr.integrate( arr, tol, 1 )
print "\n\nExample from some exam\n"
print "Computed integral = "  + str (val/2        )  
print "Analytic solution = "  + str (13.5         )  
print "Difference        = "  + str (val/2 - 13.5 )  
print "Tolerance         = "  + str (tol          )

Of course, this needs also some C\C++ files, which I chose not to post because the python file really is the issue.
For some reason, the compiler doesn't appreciate the fact that it needs to do some includes (specifically, vwrapper.h and converged.h) so I copy-pasted them in place of the #include statement.

OxbowQuest · Answer 2 · Jul 15, 2016

Hello, this post by Hans Petter Langtangen should be helpful. Further, you could consider using scipy.integrate.dblquad and tplquad instead of cubature. They might be less performant but interfacing is trivial. As for PS1, you need to give the degree of the expression (e.g. for estimating the quadrature degree). So

xpr = Expression(code, degree=2)

gets rid of the error.

External libraries in JIT compiled expression

3 Answers