Finite element methods are sensitive to the accuracy and smoothness of the mesh that approximates the curved domain boundaries.

Increasing interelement continuity has demonstrated increased accuracy on quadrilateral and hexahedral meshes.

We perform a preliminary study of the effects of surface triangle continuity, comparing a $C^0$ continuous surface with a $G^1$ continuous surface representation.

normals in different directions

different magnitude

$G^1$ continuity

normals in same direction,

different magnitude

$C^1$ continuity,

normals in same direction,

same magnitude

$C^0$ continuity obtained from the Bézier triangle, $4^{\mathrm{th}}$ order being the minimum needed for $G^1$ continuity.

$$ \mathbf{B}(u,v,w) = \displaystyle\sum_{i+j+k=4} \frac{4!}{i!j!k!}u^iv^jw^k\mathbf{P}_{ijk} $$

$G^1$ continuity obtained from the $4^{\mathrm{th}}$ order triangle of Walton and Meek.

And can be computed using the $4^{\mathrm{th}}$ order Bézier equation, blending pairs of points to form the internal Bézier points, e.g. $$ \mathbf{P}_{112} = \frac{1}{v+w}(v\mathbf{G}_{12} + w\mathbf{G}_{21}) $$

$G^1$ continuity obtained from the $4^{\mathrm{th}}$ order triangle of Walton and Meek.

Interpolation points are from surface parametrization using the optimal locations of Chen and Babuška.

For $G^1$ control points, the algorithm of Walton and Meek is used except on CAD edges, where edge tangents are used instead of surface normals.

The tetrahedral mapping we use the blending approach of Dey et al., to account for differences in exterior and interior shape functions, specifically from $G^1$ surface triangles. $$ \begin{eqnarray*} \mathbf{x}_i(\boldsymbol\xi) &=& (1-\xi_1)^{2} F_1(\xi_{F_1}') + (1-\xi_2)^{2} F_2(\xi') \\ &+& (1-\xi_3)^{2} F_3(\xi') (1-\xi_4)^{2} F_4(\xi') \\ &-& (1 - \xi_1 - \xi_2)^{2} E_1 (\xi') - (1 - \xi_1 - \xi_3)^{2} E_2 (\xi') \\ &-& (1 - \xi_1 - \xi_4)^{2} E_3 (\xi') - (1 - \xi_2 - \xi_3)^{2} E_4 (\xi') \\ & -& (1 - \xi_2 - \xi_4)^{2} E_5 (\xi') - (1 - \xi_3 - \xi_4)^{2} E_6 (\xi') \\ & +& \xi_1^{2} V_1 + \xi_2^{2} V_2 + \xi_3^{2} V_3 + \xi_4^{2} V_4 \end{eqnarray*} $$ where $\boldsymbol\xi = (\xi_1,\xi_2,\xi_3,\xi_4)$ are barycentric coordinates and $F_i$, $E_i$, and $V_i$ correspond to faces, edges, and vertices.

Curved elements implemented in Nektar++, a high-order finite element code.

Meshes are identical except for boundary surfaces and tets.

Accuracy is measured in the $L_2$ and $L_\infty$ norm of solution and $L_2$ norm of solution gradient. $$ \begin{eqnarray*} L_2(u^{h,p}) &=& \left(\int_{\Omega_M} (u^{h,p}-u_{exact})^2 \;\mathrm{d}\Omega_M \right)^{1/2} \\ L_\infty(u^{h,p}) &=& \max_{\Omega_M} \left|u^{h,p}-u_{exact}\right| \\ L_2(\nabla u^{h,p})& =& \left(\int_{\Omega_M} \nabla(u^{h,p}-u_{exact})\cdot \nabla(u^{h,p}-u_{exact}) \; \mathrm{d}\Omega_M \right)^{1/2} \end{eqnarray*} $$

Fully developed, laminar viscous flow through a cylindrical pipe of unit radius and length of 10 in the $z$- direction $$ u_z(\mathbf{x}) = -\frac{1}{4\nu}\left(\frac{\mathrm{d}p}{\mathrm{d}z}\right)(R^2-r^2), \quad p(\mathbf{x}) = 1-\frac{z}{10}, \quad r^2 = x^2 + y^2 $$ with $u_x = u_y = 0$, $\nu = 0.025$, such that $u_z = 1-r^2$.

interp. error. | ||
---|---|---|

$C^0$ | $1.19\times 10^{-5}$ | |

$G^1$ | $6.32\times 10^{-3}$ |

We first consider Poisson's equation on a sphere of unit radius. $$ -\Delta u(\mathbf{x}) = f(\mathbf{x})$$
Two curved tetrahedral meshes with 94 and 2942 tetrahedra are used.

We choose $u(\mathbf{x}) = (1-r^6)$ such that $u(\mathbf{x}) = 0$ on the boundary.

5$^{\text{th}}$ order solutions

interp. error. | $L_2(p^{h,p})$ | $L_\infty(p^{h,p}) $ | $L_2(\nabla p^{h,p})$ | ||
---|---|---|---|---|---|

$C^0$ | $6.71\times 10^{-4}$ | $8.88\times 10^{-4}$ | $2.10\times 10^{-3}$ | $1.29\times 10^{-2}$ | |

$G^1$ | $4.88\times 10^{-3}$ | $1.18\times 10^{-2}$ | $1.52\times 10^{-2}$ | $3.21\times 10^{-2}$ | |

$C^0$ | $2.01\times 10^{-7}$ | $1.96\times 10^{-6}$ | $1.90\times 10^{-5}$ | $1.86\times 10^{-4}$ | |

$G^1$ | $1.74\times 10^{-5}$ | $4.36\times 10^{-5}$ | $1.35\times 10^{-4}$ | $2.56\times 10^{-4}$ |

$x$-velocity

p | $L_2(p^{h,p})$ | $L_\infty(p^{h,p}) $ | $L_2(\nabla p^{h,p})$ | ||
---|---|---|---|---|---|

$C^0$ | 5 | $8.57\times 10^{-4}$ | $3.68\times 10^{-3}$ | $6.34\times 10^{-3}$ | |

$G^1$ | 5 | $1.67\times 10^{-3}$ | $1.99\times 10^{-2}$ | $1.36\times 10^{-2}$ | |

$C^0$ | 6 | $3.10\times 10^{-4}$ | $2.43\times 10^{-3}$ | $3.08\times 10^{-4}$ | |

$G^1$ | 6 | $1.65\times 10^{-3}$ | $1.56\times 10^{-3}$ | $1.30\times 10^{-2}$ |

Accuracy of $G^1$ dominated by geometric error.

$y-$velocity results are similar.

Pressure

p | $L_2(p^{h,p})$ | $L_\infty(p^{h,p}) $ | $L_2(\nabla p^{h,p})$ | ||
---|---|---|---|---|---|

$C^0$ | 5 | 1.71$\times 10^{-3}$ | $2.40\times 10^{-2}$ | $5.50\times 10^{-2}$ | |

$G^1$ | 5 | 4.51$\times 10^{-3}$ | $\mathbf{5.25\times 10^{-3}}$ | $\mathbf{2.75\times 10^{-2}}$ | |

$C^0$ | 6 | 4.40$\times 10^{-4}$ | $6.85\times 10^{-3}$ | $2.12\times 10^{-2}$ | |

$G^1$ | 6 | 4.47$\times 10^{-3}$ | $\mathbf{4.61\times 10^{-3}}$ | $\mathbf{1.77\times 10^{-2}}$ |

$G^1$ has lower error in several norms, but far from conclusive.

For a simple flow problem, $G^1$ continuous surface patches do not significantly improve solution accuracy compared to similar $C^0$ continuous patches.

The increased interpolation error of $G^1$ triangles appears to outweigh the potential benefits.

Geometric continuity is far from the only factor; there are relevant metrics to investigate, such as the scaled Jacobian.

The effect of curvature may also play a role, to achieve $G^1$ continuity, the derivatives may have to change rapidly to match direction on triangle edges.

Future work includes applying $G^1$ continuous boundary representations to other problems, examining its effect on forces, stresses, and other engineering quantities.

If there is shown to be a significant advantage, fifth order and higher $G^1$ triangles will be investigated.

Determining curved element validity and analyzing quality using Bézier properties.

Previous work (Remacle and others) introduced an adaptive check using Bézier subdivision.

We are investigating an elevation based adaptive check, studying and comparing the two algorithms.

With a robust and efficient validity check in place, we can extend local mesh modification operators for mesh adaptation and quality improvement:

- Refinement (splits), with subdivision
- Coarsening (collapses), with degree reduction
- Reshaping (control point movement), with convex hull property

Investigate variable-order meshes - i.e. highest order on the boundary, lowest order on in the interior.

Happy Canadian Thanksgiving!

These slides are available as a PDF here: SLIDES