The SPDE model with transparent barriers
In INLAspacetime: Spatial and Spatio-Temporal Models using 'INLA'

knitr::opts_chunk$set(echo = TRUE)

The transparent barrier model

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 @bakka2019barrier, 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 @lindgren2015bayesian, when using $\tilde{\mathbf{C}}$ in place of $\mathbf{C}$.

Implementation

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.