Locally-Conformal Perfectly Matched Layer implementation in FEniCSx

where t_1 and t_2 are the principal curvatures of \Gamma_i.
The curvilinear transformation that brings from such system back to the cartesian coordinate system is: 

\mathbf{x}(\xi_{1},\xi_{2},\xi_{3}) = \mathbf{p}(\xi_{1},\xi_{2}) + \xi_{1} \mathbf{n} (\xi_{1},\xi_{2})

Since we need it later, we calculate the Jacobian matrix of this transformation:

\mathbf{J}_{x/\xi} = \frac{\partial ( x, y,z)}{\partial \left( \xi_{1},\xi_{2},\xi_{3} \right)} = \left[ \frac{\partial \mathbf{x}}{\partial \xi_{1}} \quad \frac{\partial \mathbf{x}}{\partial \xi_{2}} \quad \frac{\partial \mathbf{x}}{\partial \xi_{3}} \right] = \left[ \mathbf{n} \quad (1+\kappa_{2} \xi_{1} ) \mathbf{t_{2}} \quad (1+\kappa_{3}\xi_{1})\mathbf{t_{3}} \right] = \left[ \mathbf{n} \quad h_{2}\mathbf{t_{2}} \quad h_{3}\mathbf{t_{3}} \right] \quad \text{(1)}

The core of Locally-Conformal PML approach is that the complex coordinate stretch (that is, the conventional transformation that makes the PML absorb all the incoming waves) is applied only on the normal direction \xi_1:

\tilde{\xi}_{1} (\xi_{1}) = \xi_{1} + \frac{1}{jk}f(\xi_{1}) \quad \text{(2)}

For which the Jacobian matrix is:

\mathbf{J}_{\tilde{\xi}/\xi} = \frac{\partial ( \tilde{\xi_{1}},\tilde{\xi_{2}},\tilde{\xi_{3}})}{\partial \left( \xi_{1},\xi_{2},\xi_{3} \right)} = \begin{bmatrix} {\frac{d \tilde{\xi}_{1}}{d\xi_{1}}}& {0 }& {0 }\\ {0 }& {1 }& {0 }\\ {0 }& {0 }& {1 }\\ \end{bmatrix} = \begin{bmatrix} {1+ \frac{1}{jk}\sigma }& {0 }& {0 }\\ {0 }& {1 }& {0 }\\ {0 }& {0 }& {1 }\\ \end{bmatrix}= \begin{bmatrix} {s_{1} }& {0 }& {0 }\\ {0 }& {1 }& {0 }\\ {0 }& {0 }& {1 }\\ \end{bmatrix}

With this in mind, we can compute the total jacobian J_{PML}, that brings us from the complexified cartesian coordinates to the real cartesian coordinates.

Thanks to the chain rule, we can split it as follows:

\mathbf{J}_{PML} = \frac{\partial \left( \tilde{x},\tilde{y},\tilde{z} \right)}{\partial \left( x,y,z \right)}=\frac{\partial \left( \tilde{x},\tilde{y},\tilde{z} \right)}{\partial (\tilde{\xi}_{1},\tilde{\xi}_{2},\tilde{\xi}_{3} )} \frac{\partial ( \tilde{\xi}_{1},\tilde{\xi}_{2},\tilde{\xi}_{3} )}{\partial \left( \xi_{1},\xi_{2},\xi_{3} \right)} \frac{\partial \left( \xi_{1},\xi_{2},\xi_{3} \right)}{\partial \left( x,y,z \right)}

Where the steps represent: 

1. the transformation from Cartesian coordinates to curvilinear coordinates in the complex space
2. the use of the complex stretch to come back in the real space
3. the transformation from curvilinear coordinates to Cartesian coordinates in the real space

The first jacobian matrix can be computed from the previous formula, with the difference of complex \xi_1:

\mathbf{J}_{\tilde{x}/\tilde{\xi} } = \frac{\partial ( \tilde{x}, \tilde{y},\tilde{z} )}{\partial \left( \tilde{\xi}_{1},\tilde{\xi}_{2},\tilde{\xi}_{3} \right)} = \left[ \mathbf{n} \quad (1+\kappa_{2} \tilde{\xi}_{1} ) \mathbf{t_{2}} \quad (1+\kappa_{3} \tilde{\xi}_{1})\mathbf{t_{3}} \right] = \left[ \mathbf{n} \quad (1+\kappa_{2} \tilde{\xi}_{1} ) \mathbf{t_{2}} \quad (1+\kappa_{3} \tilde{\xi}_{1})\mathbf{t_{3}} \right] = \left[ \mathbf{n} \quad s_{2} h_{2}\mathbf{t_{2}} \quad s_{3} h_{3}\mathbf{t_{3}} \right]

\text{where} \quad s_{i} = 1-\frac{\kappa_{i}}{jkh_{i}}f(\xi_{1})

Since the second jacobian has been computed already, we are only missing the last one. Notice that is just the inverse of (1):

\mathbf{J}_{\xi/x} = \frac{\partial ( \xi_{1},\xi_{2},\xi_{3})}{\partial \left( x, y,z \right)} = \left[ \mathbf{n} \quad h_{2}^{-1}\mathbf{t_{2}} \quad h_{3}^{-1}\mathbf{t_{3}} \right]

We can now compute the PML Jacobian matrix:

\mathbf{J}_{PML} = \frac{\partial \left( \tilde{x},\tilde{y},\tilde{z} \right)}{\partial \left( x,y,z \right)}=\mathbf{J}_{\tilde{x}/\tilde{\xi} } \mathbf{J}_{\tilde{\xi}/ \xi}\mathbf{J}_{\xi/x}=

= \left[ \mathbf{n} \quad s_{2} h_{2}\mathbf{t_{2}} \quad s_{3} h_{3}\mathbf{t_{3}} \right] \begin{bmatrix} {s_{1} }& {0 }& {0 }\\ {0 }& {1 }& {0 }\\ {0 }& {0 }& {1 }\\ \end{bmatrix} \left[ \mathbf{n} \quad h_{2}^{-1}\mathbf{t_{2}} \quad h_{3}^{-1}\mathbf{t_{3}} \right]= s_{1} \mathbf{n} \mathbf{n^{T}} +s_{2} \mathbf{t_{2}} \mathbf{t_{2}^{T}} + s_{3}\mathbf{t_{3}} \mathbf{t_{3}^{T}}

and its inverse, that is:

\mathbf{J}_{PML}^{-1} = s_{1}^{-1} \mathbf{n} \mathbf{n^{T}} +s_{2}^{-1} \mathbf{t_{2}} \mathbf{t_{2}^{T}} + s_{3}^{-1}\mathbf{t_{3}} \mathbf{t_{3}^{T}}

It is possible to demonstrate that the physics inside the PML is described by an anisotropic version of the Helmholtz equation: 

\nabla \cdot \left( \mathbf{\Lambda}_{PML} \right)+ \alpha_{PML}k^{2}\tilde{p} = -j\omega\rho_{0} \tilde{q}_{m}

where

\mathbf{\Lambda}_{PML} = \mathbf{J}_{PML}^{-1}\mathbf{J}_{PML}^{-T} det(\mathbf{J}_{PML})=\frac{s_{2} s_{3}}{s_{1}}\mathbf{n} \mathbf{n^{T}} +\frac{s_{1} s_{3}}{s_{2}}\mathbf{t_{2}} \mathbf{t_{2}^{T}} + \frac{s_{1} s_{2}}{s_{3}}\mathbf{t_{3}} \mathbf{t_{3}^{T}}

and

\alpha_{PML} = det(\mathbf{J}_{PML})\,=\,s_{1}s_{2}s_{3}