Legendre spectral finite elements

The GEBT equations are discretized with Legendre spectral finite elements (LSFEs), which are defined by \(P\) nodes and \((P-1)\)th order Lagrangian-interpolant basis functions. Nodes are located at the \(P\) Gauss-Legendre-Lobatto (GLL) points (in the element reference domain), \(\xi_i\) for \(i\in\{1,2, \ldots,P\}\), which are the solutions to the equation

\[\begin{aligned} \left( 1-\xi^2\right) \frac{\partial L_{P-1}}{\partial \xi} = 0 \end{aligned}\]

for \(\xi \in[-1,1]\), where \(L_{P-1}(\xi)\) is the Legendre polynomial of degree \((P-1)\) [@Deville-etal:2002]. The beam reference line is approximated as

\[\begin{aligned} \underline{x}^{0,h}(\xi) = \sum_{\ell=1}^{P} \phi_\ell(\xi) \underline{x}^0_\ell \end{aligned}\]

where \(\phi_\ell(\xi)\) is the Lagrangian-interpolant of the \(\ell^{th}\) node. For LSFEs, those can be written as

\[\begin{aligned} \phi_\ell(\xi) = \frac{-1}{P (P-1)} \frac{\left(1-\xi^2\right) \frac{\partial L_{P-1}}{\partial \xi}}{(\xi - \xi_\ell) L_{P-1}(\xi_\ell)} \end{aligned}\]

for \(\xi \in [-1,1]\). The weak form of the residual of the governing equations Eq. (9), for the \(i^\mathrm{th}\)-node, can be written (after integration by parts) as

(11)\[\underline{R}_i = \int_{-1}^{1} \phi_i \underline{\mathcal{R}}\, J d\xi = \int_{-1}^{1} \left[ \frac{\partial \phi_i}{\partial \xi} \left(\underline{\mathcal{F}}^\mathrm{E1} + \underline{\mathcal{F}}^\mathrm{D1} \right)J^{-1} + \phi_i \left( \underline{\mathcal{F}}^\mathrm{I} + \underline{\mathcal{F}}^\mathrm{E2} + \underline{\mathcal{F}}^\mathrm{D2} - \underline{\mathcal{F}}^\mathrm{ext} \right) \right] J d\xi\]

\(\underline{R}\in\mathbb{R}^6\), for all \(i\in\{1,\ldots,P\}\), where we have mapped the domain \(s\in[0,L]\) to \(\xi \in[-1,1]\), \(J(\xi) \in \mathbb{R}\) is the Jacobian of the mapping,

\[\begin{aligned} J(\xi) = \sqrt{\frac{\partial \underline{x}^{0,h}}{\partial \xi}^T \frac{\partial \underline{x}^{0,h}}{\partial \xi} } \end{aligned}\]

We remark that standard linear (i.e., 2-node) and quadratic (i.e., 3-node) elements are a subset of LSFEs.

As described above, the generalized displacement, \(\underline{q}\) is in \(\mathbb{R}^3 \times \mathrm{SO(3)}\). However, rotations are represented as quaternions for storage and interpolation, and are stored at nodes. In this form, nodal degrees of freedom are denoted, for the \(i^\mathrm{th}\) node,

\[\begin{split}\begin{aligned} \underline{q}_i = \begin{bmatrix} \underline{u}_i \\ \hat{q}_i \\ \end{bmatrix} \end{aligned}\end{split}\]

where \(\underline{u}_i \in \mathbb{R}^3\) and \(\hat{q}_i \in \mathbb{R}^4\). The generalized displacement in \(\mathbb{R}^3\times\mathrm{SO(3)}\) along the beam reference line is then given as

\[\begin{split}\begin{aligned} \underline{q}(s,t) = \begin{bmatrix} \underline{u}^h \\ \underline{\underline{R}} \left( \widehat{q}^h \right) \end{bmatrix} \end{aligned}\end{split}\]

where displacement is interpolated in the normal manner, i.e.,

\[\begin{aligned} \underline{u}^h = \sum_{j=1}^{p} \phi_j\underline{u}_j \end{aligned}\]

but quaternion interpolation requires normalization, i.e.,

\[\begin{aligned} \widehat{q}^h = \frac{ \sum_{j=1}^{p} \phi_j \hat{q}_j } {|| \sum_{j=1}^{p} \phi_j \hat{q}_j ||} \end{aligned}\]

For a given quaternion, the associated rotation matrix is calculated as

\[\begin{aligned} \underline{\underline{R}}\left(\hat{q}\right) = \underline{\underline{I}} + q \widetilde{q} + 2 \widetilde{q} \widetilde{q} \end{aligned}\]

Introducing a quadrature scheme with \(n^Q\) points with locations and weights, \(\xi_k^Q\), \(w_k^Q\), \(k\in \{1, \ldots, n^Q\}\), respectively, the approximate form of the residual, Eq. (11), can be written

\[\begin{split}\begin{aligned} \underline{R} = \begin{bmatrix} \underline{F}^\mathrm{I}_1 + \underline{F}^\mathrm{E}_1 + \underline{F}^\mathrm{D}_1 - \underline{F}^\mathrm{ext}_1 \\ \underline{F}^\mathrm{I}_2 + \underline{F}^\mathrm{E}_2 + \underline{F}^\mathrm{D}_2 - \underline{F}^\mathrm{ext}_2 \\ \vdots \\ \underline{F}^\mathrm{I}_P + \underline{F}^\mathrm{E}_P + \underline{F}^\mathrm{D}_P - \underline{F}^\mathrm{ext}_P \end{bmatrix} \end{aligned}\end{split}\]

where \(\underline{R} \in \mathbb{R}^{6 P}\) and

\[\begin{split}\begin{aligned} \underline{F}_i^\mathrm{I} &= \sum_{k=1}^{n^Q} J(\xi^Q_k) \phi_i(\xi^Q_k) \underline{\mathcal{F}}^I(\xi^Q_k) w^Q_k\, \\ \underline{F}_i^\mathrm{E} &= \sum_{k=1}^{n^Q} \left[\phi_i^\prime(\xi^Q_k) {\underline{\mathcal{F}}^{E1}}(\xi^Q_k)+ J(\xi^Q_k) \phi_i(\xi^Q_k) \underline{\mathcal{F}}^\mathrm{E2}(\xi^Q_k) \right] w^Q_k\, \\ \underline{F}_i^\mathrm{D} &= \sum_{k=1}^{n^Q} \left[\phi_i^\prime(\xi^Q_k) {\underline{\mathcal{F}}^{D1}}(\xi^Q_k)+ J(\xi^Q_k) \phi_i(\xi^Q_k) \underline{\mathcal{F}}^\mathrm{D2}(\xi^Q_k) \right] w^Q_k\, \\ \underline{F}_i^\mathrm{ext} &= \sum_{k=1}^{n^Q} \phi_i (\xi^Q_k) \underline{F}^{ext}(\xi^Q_k) J(\xi^Q_k) w^Q_k \end{aligned}\end{split}\]

and

\[\phi_i^\prime(\xi^Q_k) = \left .\frac{\partial \phi_i}{\partial \xi}\right|_{\xi=\xi^Q_k}\]

The matrices required for the time-integration iteration matrix in Eq. (4) are constructed from the following:

\[\begin{split}\begin{aligned} \underline{\underline{M}}_{ij} = \sum_{k=1}^{n^Q} & \phi_i(\xi^Q_k) \underline{\underline{M}}(\xi^Q_k) \phi_j(\xi^Q_k) J(\xi^Q_k) w^Q_k \\ % \underline{\underline{G}}_{ij} = \sum_{k=1}^{n^Q} & \Big\{ \phi'_i(\xi^Q_k) \underline{\underline{D}}(\xi^Q_k) \phi'_j(\xi^Q_k) \frac{1}{J(\xi^Q_k)} \\ &+ \phi_i^\prime (\xi^Q_k) \underline{\underline{\mathcal{G}}}^\mathrm{D1}(\xi^Q_k) \phi_j(\xi^Q_k) \\ & +\phi_i(\xi^Q_k) \underline{\underline{D}}^\mathrm{D2}(\xi^Q_k) \phi^\prime_j(\xi^Q_k) \\ & \phi_i(\xi^Q_k) \left[\underline{\underline{\mathcal{G}}}^\mathrm{I}(\xi^Q_k) + \underline{\underline{\mathcal{G}}}^\mathrm{D2}(\xi^Q_k) \right] \phi_j(\xi^Q_k) J(\xi^Q_k) \Big\}w^Q_k \\ % \underline{\underline{K}}_{ij} = \sum_{k=1}^{n^Q} \Big\{& \phi'_i(\xi^Q_k) \left[ \underline{\underline{C}}(\xi^Q_k) + \underline{\underline{D}}^\mathrm{D1}(\xi^Q_k) \right] \phi'_j(\xi^Q_k) \frac{1}{J(\xi^Q_k)} \\ & + \phi'_i(\xi^Q_k) \left[ \underline{\underline{K}}^\mathrm{E1}(\xi^Q_k) + \underline{\underline{K}}^\mathrm{D1}(\xi^Q_k) \right] \phi_j(\xi^Q_k) \\ & \phi_i(\xi^Q_k) \left[ \underline{\underline{\mathcal{P}}}^\mathrm{E2}(\xi^Q_k) + \underline{\underline{\mathcal{P}}}^\mathrm{D2}(\xi^Q_k) \right] \phi'_j(\xi^Q_k) \\ & + \phi_i(\xi^Q_k) \left[ \underline{\underline{\mathcal{K}}}^\mathrm{E2}(\xi^Q_k) +\underline{\underline{\mathcal{K}}}^\mathrm{D2}(\xi^Q_k) \right]\phi_j(\xi^Q_k) J(\xi^Q_k) \Big\} w^Q_k \end{aligned}\end{split}\]

for all \(i,j \in\{1,2, \ldots, P\}\) and \(\underline{\underline{M}}_{ij},\underline{\underline{G}}_{ij},\underline{\underline{K}}_{ij} \in \mathbb{R}^{6 \times 6}\). These matrices define the full matrices for a single beam LSFE. For example, the mass matrix is assembled as

\[\begin{split}\underline{\underline{M}} = \begin{bmatrix} \underline{\underline{M}}_{11}& \underline{\underline{M}}_{12}& \ldots & \underline{\underline{M}}_{1P}\\ \underline{\underline{M}}_{21} & \underline{\underline{M}}_{22}& \ldots & \underline{\underline{M}}_{2P}\\ \vdots & \vdots & \vdots & \vdots \\ \underline{\underline{M}}_{P1}& \underline{\underline{M}}_{P2}& \ldots & \underline{\underline{M}}_{PP}\\ \end{bmatrix}\end{split}\]

where \(\underline{\underline{M}} \in \mathbb{R}^{6 P \times 6 P}\); similarly for \(\underline{\underline{G}}, \underline{\underline{K}} \in \mathbb{R}^{6 P \times 6 P}\).

Deville, M. O., P. F. Fischer, and E. H. Mund. 2002. High-Order Methods for Incompressible Fluid Flow. Cambridge University Press.