This model considers an SPDE over a domain \(\Omega\) which is partitioned into \(k\) subdomains \(\Omega_d\), \(d\in\{1,\ldots,k\}\), where \(\cup_{d=1}^k\Omega_d=\Omega\). A common marginal variance is assumed but the range can be particular to each \(\Omega_d\), \(r_d\).
From Bakka et al. (2019), the precision matrix is \[ \mathbf{Q} = \frac{1}{\sigma^2}\mathbf{R}\mathbf{\tilde{C}}^{-1}\mathbf{R} \textrm{ for } \mathbf{R}_r = \mathbf{C} + \frac{1}{8}\sum_{d=1}^kr_d^2\mathbf{G}_d , \;\;\; \mathbf{\tilde{C}}_r = \frac{\pi}{2}\sum_{d=1}^kr_d^2\mathbf{\tilde{C}}_d \] where \(\sigma^2\) is the marginal variance. The Finite Element Method - FEM matrices: \(\mathbf{C}\), defined as \[ \mathbf{C}_{i,j} = \langle \psi_i, \psi_j \rangle = \int_\Omega \psi_i(\mathbf{s}) \psi_j(\mathbf{s}) \partial \mathbf{s},\] computed over the whole domain, while \(\mathbf{G}_d\) and \(\mathbf{\tilde{C}}_d\) are defined as a pair of matrices for each subdomain \[ (\mathbf{G}_d)_{i,j} = \langle 1_{\Omega_d} \nabla \psi_i, \nabla \psi_j \rangle = \int_{\Omega_d} \nabla \psi_i(\mathbf{s}) \nabla \psi_j(\mathbf{s}) \partial \mathbf{s}\; \textrm{ and }\; (\mathbf{\tilde{C}}_d)_{i,i} = \langle 1_{\Omega_d} \psi_i, 1 \rangle = \int_{\Omega_d} \psi_i(\mathbf{s}) \partial \mathbf{s} . \]
In the case when \(r = r_1 = r_2 = \ldots = r_k\) we have \(\mathbf{R}_r = \mathbf{C}+\frac{r^2}{8}\mathbf{G}\) and \(\mathbf{\tilde{C}}_r = \frac{\pi r^2}{2}\mathbf{\tilde{C}}\) giving \[ \mathbf{Q} = \frac{2}{\pi\sigma^2}( \frac{1}{r^2}\mathbf{C}\mathbf{\tilde{C}}^{-1}\mathbf{C} + \frac{1}{8}\mathbf{C}\mathbf{\tilde{C}}^{-1}\mathbf{G} + \frac{1}{8}\mathbf{G}\mathbf{\tilde{C}}^{-1}\mathbf{C} + \frac{r^2}{64}\mathbf{G}\mathbf{\tilde{C}}^{-1}\mathbf{G} ) \] which coincides with the stationary case in Lindgren and Rue (2015), when using \(\tilde{\mathbf{C}}\) in place of \(\mathbf{C}\).
In practice we define \(r_d\) as \(r_d = p_d r\), for known \(p_1,\ldots,p_k\) constants. This gives \[ \mathbf{\tilde{C}}_r = \frac{\pi r^2}{2}\sum_{d=1}^kp_d^2\mathbf{\tilde{C}}_d = \frac{\pi r^2}{2} \mathbf{\tilde{C}}_{p_1,\ldots,p_k} \textrm{ and } \frac{1}{8}\sum_{d=1}^kr_d^2\mathbf{G}_d = \frac{r^2}{8}\sum_{d=1}^kp_d^2\mathbf{\tilde{G}}_d = \frac{r^2}{8}\mathbf{\tilde{G}}_{p_1,\ldots,p_k} \] where \(\mathbf{\tilde{C}}_{p_1,\ldots,p_k}\) and \(\mathbf{\tilde{G}}_{p_1,\ldots,p_k}\) are pre-computed.