use python side quadrature and make grid function evaluation vectorize
There is for example a package quadpy https://github.com/nschloe/quadpy which we could support. For example
for e in grid.elements():
e.quadpy( lamnda point: ... , **kwargs)
Since the lambda could contain evaluations of grid functions it might be nice to make the evaluation vectorize. According to http://pybind11.readthedocs.io/en/stable/advanced/pycpp/numpy.html#vectorizing-functions the arguments have to be PODs so CoordinateWrapper would not work, which means we do need also the (element,hatx) call, which then could be vectorize for quadrature.