So you want to compute the following integral:
\begin{align*}
\frac{1}{(2\pi|\Sigma|)^{n/2}}\int_{V_k}\exp\left(-\frac{1}{2}(x-a)^T\Sigma^{-1}(x-a)\right)dx
\end{align*}
where $V_k\subset \mathbb{R}^k$ is a $k$-polytope, $a\in \mathbb{R}^k$ is the mean and $\Sigma$ is the $k\times k$ covariance matrix of $k$-variate normal vector.
My approach would be doing a change of variables $y=\Sigma^{1/2}(x-a)$ and then integrate by one variable at a time, since the integrand function will be a product of one-variable functions. If $\Sigma^{1/2}V_k$ is a $k$-cube then the integration is straightforward. For the general case it will not be easy. Take $k=2$, normal with zero mean and unit covariance matrix, and $V_2$ a triangle with vertices in $(0,0)$, $(1,0)$, $(0,1)$. Then your integral will be
\begin{align}
\frac{1}{2\pi}\int_{V_2}e^{-\frac{x_1^2}{2}}e^{-\frac{x_2^2}{2}}dx_1dx_2&=\frac{1}{2\pi}\int_{0}^{1}e^{-\frac{x_1^2}{2}}dx_1\int_{0}^{1-x_1}e^{-\frac{x_2^2}{2}}dx_2\\
&=\frac{1}{\sqrt{2\pi}}\int_0^1e^{-\frac{x_1^2}{2}}(\Phi(1-x_1)-\Phi(0))dx_1
\end{align}
which I do not think has closed formula. (I might be wrong).
As for approximation, I do not see why standard numerical integral solving methods would not work.