Python: Surface Integral over Mesh

Question

How can I compute the surface integral of a specified function of the outer surface of a mesh in python?

Code

I tried to create a minimal workable example. In this example, I try to compute the surface area of a box, which I define using trimesh.

I only compute the surface area because it can be easily checked using the trimesh package (mesh.area), but in the future I will have much more complex functions, therefore, I cannot use the build-in trimesh.integrate.

This is why I have turned to quadpy for the surface integration.

Here is the code:

import numpy as np
import quadpy
import trimesh


# CREATE THE BOX
# object (box)
box = trimesh.creation.box(extents=[0.3, 0.3, 0.3])

# show
box.show(viewer='gl')


# SURFACE INTEGRATION 

# 1 GET TRIANGLES of trimesh mesh object
triangles = box.triangles
# 2 SHAPE the triangles array into the shape accepted by quadpy
# trimesh -> (num_traingles, corners, xyz coordinates)
# quadpy ->  (corners, num_triangles, xyz_coords)
triangles_re = np.reshape(triangles,(
                       np.shape(triangles)[1],
                       np.shape(triangles)[0],
                       np.shape(triangles)[2])
                 )
# 3 DEFINE INTEGRATION FUNCTION 
def f(x):
    # -> get area of surface Integral[1, Area] = Area
    return np.ones(np.shape(x))

sol, error_estimate = quadpy.t2.integrate_adaptive(f, triangles_re,  1.0e-10)

Issue

I am getting the following error:

assert all(areas > minimum_triangle_area)

AssertionError

from Python: Surface Integral over Mesh