Part of the answer is in the comments of the question. In order to compute the directivity of an array antenna first you need to compute the electrical field produced by one of the antenna elements. Once you have computed it, you will need to compute the array factor, which is given by
$$AF(\theta, \phi) = \sum_{m = 0}^{N - 1} I_m e^{j(k \hat{r} \cdot \mathbf{r_m})}$$
with $$\mathbf{r_m}$$ being the position vector the element m of the antenna, $$I_m$$ the complex amplitude of the excitation of the element m, and
$$\hat{r} = \sin{\theta} \; cos{\phi} \; \hat{x} + \sin{\theta} \; \sin{\phi} \hat{y} + \cos{\theta} \; \hat{z}$$
The total array pattern is given by the multiplication of the array factor and the element radiation pattern, the previously computed electrical field of the antenna element. Finally, the directivity is given by
$$D = \frac{U(\theta, \phi)}{\frac{1}{4\pi} \int_0^{2\pi} \int_0^{\pi} \frac{U(\theta, \phi)}{\max{\{U(\theta, \phi)\}}} \sin{\theta} d\theta d\phi}$$
where $$U(\theta, \phi)$$ is the radiation pattern. The radiation pattern is given by
$$U(\theta, \phi) = |B(\theta, \phi)|^2$$
with $$B(\theta, \phi)$$ being the normalized total array pattern.
Please note that the spherical coordinates are given by figure 1.
Figure 1