For a given tetrahedron E, I need to compute numerically the integral of a polynomial over one of its faces e:
$$\int_e f(x,y,z) ds$$ I know that in the finite elements context, the standard practice to compute volume integrals is to perform a change of variable and do the actual computations on the unit reference element using a Gaussian quadrature, but I can't seem to figure out exactly how that would work for a surface integral.
First, you create a diffeomorphism $\Phi$ that maps the triangle $T = \{(0,0,0),(1,0,0),(0,1,0)\}$ the face $e$.
Then, you employ the change of variables for integrals
$$\int_e f(x,y,z)\, ds = \int_T f(\Phi(\hat{x},\hat{y},\hat{z})) \vert \det D\Phi \vert \, d\hat{s}\,.$$
Finally, you approximate the right-hand side with a 2D Gauss quadrature extended by zero to the $\hat{z}$ component.