09/06/2024 - This is an old post I wrote when I started the FEM project. This was an important literature for understanding the implementation of PML, but took me quite some time to reproduce in C++ since some of the equations in the paper are wrong. I’m posting it here in case someone like me needs to read it in detail.

Recently I’ve been trying to replicate Perfect Matched Layer described in the paper “Matzen, René. “An efficient finite element time‐domain formulation for the elastic second‐order wave equation: A non‐split complex frequency shifted convolutional PML.” International Journal for Numerical Methods in Engineering 88.10 (2011): 951-973.” As an important paper in the PML field, it turns out that there are a lot of notation inconsistencies and definition mistakes in this paper which misled me and wasted me quite some time. As I’ve sorted out the details, I hope my clarification here can save you some time.

2.1 Frequency domain equations

• Equation (16) \(D_0 = \begin{bmatrix}0 & 1 & 0 & 0 \\1 & 0 & 0& 0 \\ 0 & 0 & \color{red}{\frac{1}{2}} & 0 \\ 0 & 0& 0& \color{red}{-\frac{1}{2}} \end{bmatrix}\)

The new definition can be verified from equation (12) and equation (14).

2.2 Time domain equations

• Equation (17a) \(\nabla \cdot \sigma ^{\prime} + {\color{red}{\rho{\prime}}} \pmb{p} = \mathcal{L_0}(t) {\color{red}{\rho}} \pmb{u}\)

Since the source is in the original domain, in this case \(\rho = \rho{\prime}\)

• Equation (17c) \(\mathbf{C}\prime = \mathbf{C_0} + \mathcal{L_1}(t){\color{red}{\mathbf{C_1}}} + \mathcal{L_2}(t)\mathbf{C_2}\)

2.3 Finite element implementation

• Equation (24) The displacement within an element is expanded using the basis function \(u \approx \sum_{k=1}^{m_e}\mathit{N_k}\mathbf{I}\mathbf{d_k^e}=[\mathit{N_1}\mathbf{I}, \mathit{N_2}\mathbf{I}, ..., \mathit{N_{m_{e}}}\mathbf{I}] [\mathbf{d_1} , \mathbf{d_2}, ..., \mathbf{d_{m_{e}}} ]^{T}\)

where \(\mathbf{d_k} = (d_1, d_2)^{T}\).

It needs to be clarified that in equations (24), the \(\mathbf{N_i}\) and \(\mathbf{N_j}\) denote one \(\color{red}{column}\) of the \(\mathbf{N}\) matrix, rather than a \(2\times 2\) sub matrix.

Similarly, since \(\mathbf{B_i} = \partial{\mathbf{N_i}}\), then \(\mathbf{B_i}\) is only one \(\color{red}{column}\) of the classical strain-displacement matrix, rather than the whole matrix itself.

• Equation (25) \(\mathbf{f_i^e}=\int_{\Omega^{e}}{\color{red}{\rho}}\mathbf{N_i^T}\pmb{p}d\Omega\)

2.4 Time integration

It should be noted that the time integration scheme given in equation (30) is \(\color{red}{not}\) the classical Newmark-beta as it stated. When I coded the time marching as given in the paper, no wave is generated. The reader should refer to Wikipedia or other related papers for the exact update algorithm.