## 2 Basic Theory of HSDTPEs
### Fictitious Thin Layer between Adjacent Triangular Elements  
According to the definition of the displacement approximation in Equation (5), the deformation along the share boundary of the adjacent elements \(e_{i}\) and \(e_{j}\) is discontinuous, which means that \(\mathbf{u}^{e_{i}}\big{(}x_{p},y_{p},z_{p}\big{)}\
eq\mathbf{u}^{e_{j}}\big{(} x_{p},y_{p},z_{p}\big{)}\) for an arbitrary point \(\mathrm{p}\) along the share boundary shown in Figure 2, where point \(\mathrm{p}\) has local coordinates \(\big{(}s_{p},n_{p},z_{p}\big{)}\) and global coordinates \(\big{(}x_{p},y_{p},z_{p}\big{)}\). Here, we introduce a fictitious thin layer \(e_{l}\) shown in Figure 3 to enforce the continuous condition over the share boundary of the elements. The geometry dimensions of \(e_{l}\) are the length \(l\), width \(d\), and height \(h\), where \(d\ll l\), \(d\ll h\), and \(h\) is the thickness of the plate in the transverse direction. Because \(d\ll l\) and \(d\ll h\), the strain-displacement relations \(\varepsilon^{l}=[\gamma_{ns},\varepsilon_{n},\gamma_{nz}]^{\mathrm{T}}\) in thin layer \(e_{l}\) can be simplified as:  
\[\gamma_{ns}=\frac{\partial\overline{u}}{\partial n}\approx\frac{\overline{u} \overline{v}^{p^{i}}-\overline{u}\overline{v}^{p^{i}}}{d},\varepsilon_{n}= \frac{\partial\overline{v}}{\partial n}\approx\frac{\overline{v}\overline{v} ^{p^{j}}-\overline{v}\overline{v}^{p^{j}}}{d},\gamma_{nz}=\frac{\partial \overline{u}}{\partial n}\approx\frac{\overline{v}\overline{v}^{p^{j}}- \overline{v}\overline{v}^{p^{j}}}{d} \tag{17}\]  
where \(\overline{u}\overline{v}^{p^{j}}=\left[\overline{u}\overline{v}^{p^{j}}, \overline{v}\overline{v}^{p^{j}}\right]^{\mathrm{T}}\) is the displacement of point \(\mathrm{p}\) in the local coordinate \((s,n,z)\) computed by the approximation of triangular element \(e_{i}\), and \(\overline{u}\overline{v}^{p^{j}}=\left[\overline{u}\overline{v}^{p^{j}}, \overline{v}\overline{v}^{p^{j}}\right]^{\mathrm{T}}\) is the displacement of point \(\mathrm{p}\) in local coordinate \((s,n,z)\) computed by the approximation of triangular element \(e_{j}\). \(\overline{u}\overline{v}^{p^{j}}\) and \(\overline{u}\overline{v}^{p^{j}}\) can be calculated using Equation (5). Substitution of Equation (5) into Equation (17) yields:  
\[\varepsilon^{l}=\frac{1}{d}\mathbf{N}^{l}\mathbf{a}^{l} \tag{18}\]  
where  
\[\mathbf{N}^{l}=\mathbf{\Lambda}^{l}\left[-\mathbf{N}^{e_{i}}\left(x_{p},y_{p},z_{p}\right)\quad\mathbf{N}^{e_{j}}\left(x_{p},y_{p},z_{p}\right)\right] \tag{19}\]  
where \(\mathbf{N}^{e_{i}}\left(x_{p},y_{p},z_{p}\right)\) is the shape function of point \(p\) in triangular element \(e_{i}\) and \(\mathbf{N}^{e_{j}}\left(x_{p},y_{p},z_{p}\right)\) is the shape function of point \(p\) in triangular element \(e_{j}\). \(\mathbf{\Lambda}^{l}\) is the transformation matrix of point \(p\) from the global coordinate \(\left(x_{p},y_{p},z_{p}\right)\) to the local coordinate \(\left(s_{p},n_{p},z_{p}\right)\), where:  
\[\mathbf{\Lambda}^{l}=\left[ \begin{array}{ccc}\cos\theta&\sin\theta&0\\ -\sin\theta&\cos\theta&0\\ 0&0&1\end{array} \right] \tag{20}\]  
and  
\[\mathbf{a}^{l}=\left[ \begin{array}{c}\mathbf{a}^{e_{i}}\\ \mathbf{a}^{e_{j}}\end{array} \right] \tag{21}\]  
where \(\mathbf{a}^{e_{i}}\) is the DOFs of element \(e_{i}\) and \(\mathbf{a}^{e_{j}}\) is the DOFs of element \(e_{j}\). Figure 3: A fictitious thin layer between adjacent triangular elements. The stress-strain relations in fictitious thin layer \(e_{l}\) are then given by:  
\[\mathbf{\sigma}^{l}=\mathbf{D}^{l}\varepsilon^{l} \tag{22}\]  
where \(\mathbf{\sigma}^{l}=\left[\tau_{\text{ns}},\sigma_{n},\tau_{nz}\right]^{\text{T}}\) and  
\[\mathbf{D}^{l}=\left[ \begin{array}{ccc}G_{0}&0&0\\ 0&E_{0}&0\\ 0&0&G_{0}/k\end{array} \right] \tag{23}\]  
where \(G_{0}=\frac{E}{2(1+\bar{v})}\) and \(E_{0}=\frac{E}{1-\bar{v}^{2}}\).