Way too late for this answer to be useful. In fact, way too late even for it to have been useful for me, when I had to solve the same problem, back in maybe 1978 or so?
But what I would do now, is tessellate the region of interest, dissecting into tetrahedrae. In the case where the mode of the Gaussian is inside the region, I would insure that point was one of the nodes of said dissection. If the Gaussian mode is not inside the region, then just dissect it arbitrarily. You can dissect a cube into either 5 or 6 tetrahedrae - your choice. The dissection into 6 pieces is nice in a way, in that they are all similar tetrahedrae.
Then inside every simplex, use a Gaussian style quadrature, designed to work on a simplex. I recall finding some that apply to 3-d. Perhaps better yet would be to use a simple rule, and then split each simplex into a set of smaller similar simplexes. Use a Richardson estrapolation to measure convergence and speed it up, subdividing only those pieces as needed.
I would add the above scheme works reasonably well for 2-d problems.