Numerical integration on a constrained surface

I am interested in the numerical calculation of a surface integral in the space, with a quadratic constraint determined by a non-diagonal matrix g, which effectively bounds the integral on the surface of an ellipsoid. So far, I tried to include the constraint under the sign of integral via the Gaussian approximation of the Dirac's delta function, but this slowed significantly the convergence.

Is there an optimal library to implement efficiently such integration? Possibly in Python?

A good approach is to first generate a triangulation of the domain, e.g., with pygalmesh

import pygalmesh
import numpy as np


class Surface(pygalmesh.DomainBase):
    def __init__(self):
        self.G = np.array([
            [1.0, 3.0, -1.0],
            [2.0, 1.0, 2.0],
            [0.0, 0.0, 1.0],
        ])
        super().__init__()

    def eval(self, x):
        return np.dot(x, self.G @ x) - 1.0

    def get_bounding_sphere_squared_radius(self):
        return 10.0


d = Surface()
mesh = pygalmesh.generate_surface_mesh(d, max_radius_surface_delaunay_ball=0.1)
mesh.write("out.vtk")

Then you can use quadpy to integrate any function over the triangulation.