I need your help. My sample is a 3D cuboid. I want to take a 2D integration over a cross section parallel to the x-y plane and cuts the z-axis at z = z0. I only know Integrate(function, mesh) for an integration over the whole sample.
How can I do the cross section integration with NGSolve?

If the plane z=z0 is part of the mesh as a boundary you can integrate on the boundary using a proper region flag: defined=mesh.Boundaries("myplane").
Otherwise there is no core functionality in NGSolve. However, you can use ngs-xfem, which is a small add-on module for things like that ( see ngsxfemhttps://github.com/ngsxfem/ngsxfem ). There you have to provide a level set function where the zero level corresponds to your plane, i.e. this would look like