A comparison of $C^0$ and $G^1$ continuous curved tetrahedral meshes for high-order finite element simulations

Daniel W. Zaide, Qiukai Lu, and Mark Shephard

Scientific Computation Research Center, Rensselaer Polytechnic Institute


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.

Surface Continuity

$C^0$ continuity
normals in different directions
different magnitude

$G^1$ continuity
normals in same direction,
different magnitude

$C^1$ continuity,
normals in same direction,
same magnitude

C0 Continuity
G1 Continuity
C1 Continuity

$4^{\mathrm{th}}$ Order Bézier Triangle

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

G1 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$ Continuous Triangle

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

G1 Continuity

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}) $$

Mesh Curving

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

Interpolation PointsControl Points Control Points

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.

Blended Tetrahedron

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*} $$

Hagen-Poiseuille Flow

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}$

Poisson's Equation

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.


Poisson's Equation

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}$

Hagen-Poiseulle Flow


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.

Hagen-Poiseulle Flow


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.

Other Lessons

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 on Surface Continuity

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.

Current Research

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.

Current Research

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.


A Meshing Cartoon

Happy Canadian Thanksgiving!

These slides are available as a PDF here: SLIDES