Area and volume of P2 elements

Context:

For a fluid solver, I need to compute areas and volumes of curved elements. My curved elements are quadratic triangles and quadratic tetrahedra, defined by their Lagrange nodes.

Question:

1. Consider a 3d quadratic $$\mathbb{P}_2$$ tetrahedron defined by the positions of its 10 Lagrange nodes. Its volume is the integral of the Jacobian. Since the Jacobian is a degree 3 polynomial, a degree 3 quadrature will give the exact answer. Is there a smart way of obtaining an expression that is as simple as possible (and as fast as possible to compute) ?

2. Consider now a 3d quadratic $$\mathbb{P}_2$$ triangle defined by the positions of its 6 Lagrange nodes. Now the Jacobian to be integrated has a square root (length of the cross product of the two gradient vectors). Is there another way of doing the computation that gives a polynomial ? In my specific case, the $$\mathbb{P}_2$$ triangle is on a sphere, can we make use of that to simplify the computation ?

As you may be aware, this is a fairly common issue in finite element models, since these are all about calculating integrals on shapes including n-simplices.

For your first question, the typical general approach for the codes I'm used to is simply to differentiate the reference to physical transformation, $$\mathbf{x} = \sum_i \mathbf{x}_i N_i(\xi,\eta, \zeta),$$ where the $$\mathbf{x}_i$$s are the positions of the Lagrange nodes and the $$N_i$$s are the equivalent shape functions, which for P2 look like $$\xi(2\xi-1)$$, $$4\xi\zeta$$, $$\ldots$$ etc. The full set can be fairly trivially generated by introducing a fourth coordinate, $$\sigma = 1-\xi-\eta-\zeta$$, then taking the four versions of the form $$a(2a-1)$$ and 6 combinations $$4ab$$.

The Jacobian matrix is then $$\mathcal{J}=\left(\begin{array}\mathcal{J}_{x\xi} & \mathcal{J}_{x\eta} & \mathcal{J}_{x\zeta}\\ \mathcal{J}_{y\xi} & \mathcal{J}_{y\eta} & \mathcal{J}_{y\zeta} \\ \mathcal{J}_{y\xi} & \mathcal{J}_{z\eta} & \mathcal{J}_{z\zeta}\end{array}\right)$$ consists of terms such as $$\mathcal{J}_{x\xi}=\sum_i x_i\frac{\partial N_i}{\partial \xi}$$.

If hand optimising, this can be calculated more quickly by noting that eg. only the shape functions involving $$\xi$$ or $$\sigma$$ vary with $$\xi$$, so only 6 terms appear in each sum, which are linear functions of the reference coordinates. Typically the most complicated part is making sure that the numbering is consistent everywhere. The final operation is calculating the determinant.

None of this is particularly cheap, so it makes sense (unless your mesh is changing every tilmestep, and frequently even then) to cache the results. If you do that, then the cost is amortised across your time steps, and there isn't actually that much to be gained by trying to be much cleverer.

For your second question, if you are genuinely working on spherical triangles and really only want the areas, then you can probably get a significant time saving using Girard's Theorem, $$\mbox{Area of spherical triangle} = \mbox{radius}^2\left(\sum \mbox{angles}-\pi\right).$$

Here the angle at vertex $$\mathbf{a}$$ can be found using $$\cos\alpha = \frac{|\mathbf{a}\times\mathbf{b}\cdot\mathbf{a}\times\mathbf{c}|}{\|\mathbf{a}\times\mathbf{b}\|\|\mathbf{a}\times\mathbf{c}\|} .$$

All told, you'll have three cross products, three normalisations, three dot products and three trig functions, which I suspect will still beat trying to use high order quadrature and affine transformations.

• Thank you for the detailed answer. My problem is that my triangles are not spherical triangles (they are on a sphere but edges are not great circle, but they are quadratic). – BrunoLevy Mar 6 at 16:46
• @BrunoLevy You mean the triangles aren't actually on the sphere, just the nodes are? In that case you may be out of luck on neat short cuts. – origimbo Mar 6 at 17:58
• The triangles are on the sphere. The edges of the triangles are quadratic curves included in the sphere (but they are not great circles). – BrunoLevy Mar 6 at 18:14
• @BrunoLevy. You may want to double check that. For the edges to be embedded in the sphere implies a set of constraints on the points that it's not immediately obvious have nice non-trivial solutions. – origimbo Mar 6 at 20:03
• take a plane that intersects the sphere but that does not contain the center of the sphere. The intersection with the sphere is a circle that is not a great circle. Pick three points on that circle. The edges of my triangles are like that. – BrunoLevy Mar 6 at 20:45