One can frequently read that the product of two multivariate Gaussian pdfs, $f_1(x)$*$f_2(x)$, is itself a Gaussian function, with parameters as defined for example in:
http://www2.imm.dtu.dk/pubdb/views/edoc_download.php/3274/pdf/imm3274.pdf (section 8.1.8 page 42)
But what if the parameters (mean vector, covariance matrix) for $f_2$ are not given for all elements of $x$? I.e. if the mean vector of $f_2$ contains values only for some of the random elements of $x$, and the $f_2$ covariance matrix - which is here assumed to be diagonal - also contains only values for this subset of the elements of $x$? Is the product then still a Gaussian and how can the parameters be calculated?
This is an incorrect statement no matter how many times or in how many different places you have read it; indeed it does not hold even for univariate Gaussian random variables. If the link you have provided does indeed say exactly what you have written, then I suggest that you apply the instructions given by Dorothy Parker in a famous book review: "This is not a book to be put down lightly; it should be thrown with great force."
What is true is that if $X$ and $Y$ are independent (univariate) Gaussian random variables with densities $f_X(x)$ and $f_Y(y)$ respectively, then the product $f_X(x)f_Y(y)$ equals the joint density function $f_{X,Y}(x,y)$. Note that the joint density is a function of two real variables which have different names $x$ and $y$, and it is not claimed that $f_X(x)f_Y(x)$ is a Gaussian density (note the lack of $y$).
More generally,if $\mathbf X = (X_1,X_2,\ldots, X_m)$ and $\mathbf Y = (Y_1, Y_2, \ldots, Y_n)$ are random vectors that are independent of each other and $\mathbf X$ and $\mathbf Y$ have multivariate Gaussian densities $$\begin{align} f_{\mathbf X}(\mathbf x) &= \frac{1}{(2\pi)^{m/2}\sqrt{|\operatorname{det}(\Sigma_X)|}}\exp\left(-\frac 12(\mathbf x - \mathbf m_X)\Sigma_X^{-1} (\mathbf x - \mathbf m_X)^T\right)\\ f_{\mathbf Y}(\mathbf y) &= \frac{1}{(2\pi)^{n/2}\sqrt{|\operatorname{det}(\Sigma_Y)|}}\exp\left(-\frac 12(\mathbf y - \mathbf m_Y)\Sigma_Y^{-1} (\mathbf y - \mathbf m_Y)^T\right) \end{align}$$ where $\Sigma_X$ and $\Sigma_Y$ are $m\times m$ and $n \times n$ covariance matrices, and $\mathbf m_X$ and $\mathbf m_Y$ are the mean vectors, of $\mathbf X$ and $\mathbf Y$ respectively. Note that these densities are functions of $m$ and $n$ real variables $(x_1, x_2, \ldots, x_m)$ and $(y_1, y_2,\ldots, y_n)$ respectively. Then $$f_{\mathbf X, \mathbf Y}\big ((\mathbf x, \mathbf y)\big) = f_{\mathbf X}(\mathbf x)f_{\mathbf Y}(\mathbf y)$$ where $(\mathbf x, \mathbf y) = (x_1, x_2, \ldots, x_m, y_1, y_2,\ldots, y_n)$ is a vector of $m+n$ real numbers. Plugging and chugging, we have that $$\begin{align} f_{\mathbf X, \mathbf Y}\big ((\mathbf x, \mathbf y)\big) &= \frac{\exp\left(-\frac 12(\mathbf x - \mathbf m_X)\Sigma_X^{-1} (\mathbf x - \mathbf m_X)^T-\frac 12(\mathbf y - \mathbf m_Y)\Sigma_Y^{-1} (\mathbf y - \mathbf m_Y)^T\right)}{(2\pi)^{(m+n)/2} \sqrt{|\operatorname{det}(\Sigma_X)|\cdot|\operatorname{det}(\Sigma_Y)|}}\tag{1}\\ &= \frac{\exp\left(-\frac 12(\mathbf z - \mathbf m)\Sigma^{-1} (\mathbf z - \mathbf m)^T\right)}{(2\pi)^{(m+n)/2} \sqrt{|\operatorname{det}(\Sigma)|}}\tag{2} \end{align}$$ where $\mathbf z = (\mathbf x, \mathbf y) = (x_1, x_2, \ldots, x_m, y_1, \ldots, y_n)$, $\mathbf m = (\mathbf m_X,\mathbf m_Y)$ is the mean vector of the $(m+n)$-dimensional multivariate Gaussian vector $\mathbf Z = (\mathbf X, \mathbf Y)$, and $$\Sigma = \left[\begin{matrix}\Sigma_X&\mathbf 0\\\mathbf 0& \Sigma_Y\end{matrix}\right]\tag{3}$$ is the $(m+n)\times (m+n)$ covariance matrix. The derivations are greatly simplified by observing the block structure $(3)$ of the covariance matrix.