Genotype to phenotype map and the G-matrix
From the pheno-geno map and the genotypes frequencies, you have the whole distribution of phenotypes in your population. The mean of the phenoype $n$, $P_n$ is the 

$$\bar P_n = \sum f_{G_i} P_{G_i}$$

, where $f_{G_i}$ is the proportions of individuals having genotype $G_i$ and $P_{G_i}$ is the phenotype of the individuals with genotype $G_i$. Therefore, the variance of $P_n$ is 

$$var(P_n) = \sum f_{G_i} (P_{G_i} - \bar P_n)^2$$.

The covariance between the phenotypic traits $n$ and $m$ is 

$$cov(P_n,P_m) = \sum f_{G_i} (P_{G_i} - \bar P_n)(P_{G_i} - \bar P_m)$$

The only assumptions you need is that there is no environmental variance. I don't think the above calculations required any other assumption. $G$ can be defined at the haploid or diploid stage, it doesn't change anything.

The above calculations however give you the genetic variance and covariance but not the additive genetic variance covariance. I think (but I might be mistaken) that the additive genetic variance would be equal to 

$$\sigma_D^2 = \sum_i f_i^2(1-f_i)^2(2 \cdot P_{n,i,12} - P_{n,i,11} - P_{n,i,22})^2$$

$$\sigma_A^2 = \sum_i 2f_i(1-f_i)(f_iP_{n,i,11}+(1-2f_i)P_{n,i,12} - (1-f_i)P_{n,i,22})^2$$

, where $f_i$ is the frequency of the allele $A$ at the locus $i$ (you should get those data by reformatting your geno-pheno map matrix) and $P_{n,i,11}$, $P_{n,i,12}$ and 

$$P_{n,i,22} =\frac{1}{l-1} \sum_j P_{n,i,j,22}$$

is the $n^{th}$ phenotypic trait for the genotype $22$ and locus $i$ averaged over all other genotypes. There are $l$ genotypes in total. Of course, this works only for biallelic loci. I don't know what it would look like for more alleles.