Integrate 3D solution along an axis with dolfinx

I’m looking for a fast method to integrate a solution u(x,y,z) of a Poisson equation in a 3D domain along one of its axis such that I can obtain a 2D map: g(x,y) = \int_{z_0}^{z_1} u(x,y,z) dz . Currently I am evaluating u(x,y,z) on a regular grid using the compute_colliding_cell method as in this membrane deflection example , then just numerically integrated over this regular grid. However, this currently takes a while to evaluate u over the large 3d grid so I wonder if there is a more efficient way to perform this integration.

In 3D problems based in the discussion here there doesn’t seem to be another way yet, unless I totally misunderstand. Still, if anyone has more input I am also in the answer to this.